1 Introduction

Paralichthys olivaceus is a commercially important teleost fish that is widely farmed across the world, especially in China, Japan, and Korea (Xiu et al., 2019). Bacteria, parasites and viruses are well-known infectious pathogens (Kim et al., 2010), which cause mixed infectious illnesses, countless fatalities, and massive economic losses (Xu and Zhang, 2013). Edwardsiella tarda is a Gram-negative bacterium that can cause intestinal or extraintestinal infections in fish. In previous studies, Danio rerio (Liu et al., 2014), Sebastes schlegelii (Cao et al., 2020), tilapia (Wu et al., 2018), Scophthalmus maximus L. (Chen et al., 2020b) and P. olivaceus (Tang et al., 2019) with intestinal inflammation are more susceptible to E. tarda. Meanwhile, the antimicrobial response during E. tarda infection might inhibit the growth of intestinal microbiota. Particularly, P. olivaceus infected with E. tarda shows various intestinal and extraintestinal diseases, causing huge losses in the aquaculture industry (Xiu et al., 2019, 2021). Therefore, researches on the innate immune regulation mechanisms of flounder may contribute to the development of disease-control measures and enhance the development of flounder aquaculture.

MicroRNAs (miRNAs) are a class of non-coding short RNA molecules with 22 – 24 nucleotides which play the significant role in transcriptional regulation through blocking mRNA translating or encouraging mRNA breakdown (Hatfield et al., 2005; Wienholds and Plasterk, 2005). Previous studies have suggested that miRNAs were attractive tools and targets for novel therapeutic approaches in cancer. In human, miR-144 is involved in the pathological development and progress of many malignancy tumors (Shi et al., 2019; Kooshkaki et al., 2020). Moreover, miR-144 has the ability to promote the pathological progress of some cancers through suppressing or accelerating tumor cell proliferation, migration and invasion (Cui and Wang, 2015; Cheng et al., 2020). On the other hand, miR-144 also plays a vital role in regulating the maturity, differentiation and immune functions of immune cells (Li et al., 2018). Particularly, ectopic expression of miR-144 can effectively block cell barrier disruption and hyperpermeability to regulate activation of inflammatory pathway in response to proinflammatory agents (Siddiqui et al., 2019).

In fish, substantial progress was made in recognizing the role of miR-144 during the last decade. On the one hand, with the progress and development of high-throughput sequencing technologies, miR-144 was differentially expressed in response to external stimulus in several fish (Wu et al., 2021). Moreover, miR-144 regulated a variety of physiological activities in vivo by targeting different genes. For example, miR-144 regulats erythropoiesis via repressing Dicer (Kretov et al., 2020) as well as hematopoiesis and vascular development by targeting meis1 (Su et al., 2014) in zebra fish. Additionally, miR-144 also mediates high fatinduced changes of cholesterol metabolism via targeting C/EBPα in yellow catfish (Chen et al., 2020a). Moreover, in terms of antibacterial function, it has been observed that in response to the infection of different bacteria including S. iniae (Liu et al., 2020b), E. tarda (Li et al., 2019), and megalocytivirus (Zhang et al., 2014), pol-miR-144-5p expression varied dramatically, indicating miR-144-5p might be a regulator of antibacterial immune response. In innate immunity, miR-144 were proposed to influence innate immunity in fish by binding immunity-related genes. Xu' group have reported that miR-144 played regulatory roles in inflammatory responses by regulating the NF-κB signaling pathway via targeting NOD1 (Chu et al., 2021) and IL-1β (Yan et al., 2019) in miiuy croaker. Thus, miR-144 plays important roles with antimicrobial function as well as regulation of immune and inflammatory responses in innate immunity.

In this study, we aim to investigate the underlying regulation mechanism of miR-144-5p in intestinal mucosal immunity of flounder. Firstly, we effectively inhibited miR-144-5p expression by miR-144-5p interference in vivo. Then, we detected the differentially expressed genes (DEGs) by high-throughput sequencing technology. Next, in order to confirm the immunity-related signal pathways influenced by miR-144-5p, DEGs were analyzed by Gene Ontology (GO) functionality enrichment, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. At last, we predicted protein-protein interaction networks. These findings will contribute to a better knowledge of the immunological regulation mechanisms of miR-144-5p in intestine of flounder and other fishes.

2 Materials and Methods 2.1 Experimental Fish and miRNA Antagomir

All healthy flounders were offered by the Huanghai Aquaculture Company (Haiyang, Shandong, China). The fish were approximately 24 weeks old, with a body length of 12 cm ± 4 cm (mean ± SD) and a weight of 10 g ± 3 g (mean ± SD). Before the trials, they have been habituated in a recirculating water system at 17℃ ± 1℃ for 7 days. miR-144-5p antagomir, and antagomir-NC (negative control of miR-144-5p antagomir) utilized in this study were synthesized by Genepharma, China.

2.2 In vivo Injection and qRT-PCR

Healthy P. olivaceus were randomly separated into three groups (10 fish per group). Two experimental groups were interfered by 2 μg g⁻¹ miR-144-5p antagomir or antagomir-NC, while the control group was without any treatment. Twenty four hours after injection, the intestine tissue was collected under aseptic conditions. To examine the expression level of miR-144-5p, the collected intestines of flounders was used for RNA isolation and qRT-PCR.

The PrimeScript 1st Strand cDNA Synthesis Kit (Takara, Japan) has been used to create the cDNA. Then, based on the sequence of miR-144-5p, appropriate primers were constructed, with Po5s serving as an internal reference. The CFX96 real-time fluorescent quantitative PCR system (Bio-Rad, USA) with TB GreenTM Premix Ex TaqTM Ⅱ was used for the qRT-PCR (TaKaRa, Japan). The amplifying procedure included one cycle at 95℃ for 30 s, followed by 40 cycles of 5 s at 95℃ and 1 min at 60℃, then a melted curve from 60 to 95℃. The data were presented as means ± SE for three variables, and statistical analysis was conducted through SigmaPlot 11.0.

2.3 Transcriptome Sequencing Library Construction

TRIzol Reagent (Invitrogen, USA) has been utilized to isolate total RNA according with the manufacturer guidelines. And 1% agarose gel electrophoresiswas employed to determine the quality and quantity of extracted RNA. Furthermore, the NanoPhotometer ® spectrophotometer (IMP-LEN, CA, USA), the RNA Assay Kit in the Qubit ® 2.0 Flurometer (Life Technologies, CA, USA) and the Bioanalyzer 2100 system's RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA) has been utilized to measure RNA purity, concentration and integrity, respectively. Approximately 22.335 μg of total RNA was subjected to isolate Poly(A) mRNA with poly-T oligo attached magnetic beads. Cracking has been performed in First Strand Synthesis Reaction Buffer (5X) utilizing divalent cations at high temperatures. RNaseH has been used to breakdown the RNA after the first strand cDNA was generated using a randomized hexamer primer and M-MuLV Reverse Transcriptase. Following that, second strand cDNA synthesis has been carried out utilizing DNA Polymerase Ⅰ and dNTP. The library fragments were reverse-transcribed to identify cDNA segments with 370 – 420 bp length (Beckman Coulter, Beverly, USA) and then confirmed the quality of the database. Finally, the libraries were analyzed on an Illumina NovaSeq 6000 platform, yielding 150 bp paired-end reads following cluster generation.

2.4 Transcriptome Assembly and Novel Transcripts Prediction

Raw data were handled and then clean data were acquired via using in-house perl programs including remove adapter-containing readings, N base-containing reads, and low-quality reads in original data. Meanwhile, the Q20, Q30, and GC contents of clean data were determined. Afterwards, standard genome and genetic modeling annotation data were obtained directly from the genomic site (ftp: /ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/970/005/GCF 001970005.1 Flounder ref guided V1.0/). Hisat2 (v2.0.5) was used to create the reference genome index, and mapped clean reads to the reference genome. Furthermore, the mapped reads of each sample were assembled by String-Tie (v1.3.3b) in a reference-based approach.

2.5 Different Expression Genes Analysis and miRNA-mRNA Prediction

The featureCounts v1.5.0-p3 program was used to check the amount of reads that matched each transcript. The FPKM of every transcript was computed based on its length and the number of reads linked to it. FPKM stands for the assumed number of fragments per kilobase of transcript sequence per million base pairs sequenced, and it takes into account the impact of assembling depth and genomic length on read count simultaneously, which was presently used to estimate gene expression levels as the most commonly method. The DESeq2 R package (1.20.0) was used to analyze the different expressions of the genes in the three groups (3 biological replicas per condition). Padj ≤ 0.05 and |log₂ (Fold change)| ≥ 1 were chosen as the cutoffs for significantly different expression. The potential target genes of miR-144-5p were predicted using the TargetScan 8.0 database. Total of 3003 differentially expressed transcripts were submitted and further predicted.

2.6 GO and KEGG Enrichment Analysis of Differentially Expressed Genes

GO enrichment assessment and KEGG pathway studies have been conducted on all of the DEGs to better understand their functional relevance. GO enrichment analysis of differentially expressed genes was implemented by the ClusterProfiler R package (3.8.1). KEGG is a database for understanding key functions and utilities of organisms, including the cell, embryo, and ecosphere, through the molecular-level data, especially large-scale molecular datasets produced by genetic analysis and other high-throughput experimental techniques. The ClusterProfiler R package (3.8.1) was used to examine the statistical enrichment of pathways in KEGG networks.

2.7 Establishing Protein-Protein Interaction Networks

Protein-protein interaction networks were generated to evaluate the connection of genes in the immune regulatory pathway by using PPI (protein-protein interaction) and WGCNA (Weighted correlation network analysis). To further explore the relationship between the genes in the immunological regulatory system, protein-protein interaction networks were generated using Cytoscape 3.9.0 (https://github-releases.org/) with default settings.

3 Results 3.1 Effect of miR-144-5p Interference

To examine the effect of miR-144-5p interference, the expression of miR-144-5p was determined by qRT-PCR. miR-144-5p in flounder can be knocked down by miR-144-5p antagomir. The qRT-PCR results showed that in the presence of miR-144-5p antagomir, the expression of miR-144-5p was significantly reduced (Fig. 1A). In contrast, the presence of antagomir-NC had no effect on miR-144-5p expression (Fig. 1B).

3.2 Results of Sequences Assembly

The Illumina-based RNA-Seq was performed on three groups each with three replicates in intestine of flounder 24 hours after the injection of antagomir (named as Ant_ 24H_1, Ant_24H_2, Ant_24H_3), antagomir-NC (named as Ant_NC_24H_1, Ant_NC_24H_2, Ant_NC_24H_3) or without any treatment (named as Control_1, Control_2, and Control_3). We acquired 137343982, 141233938 and 142285934 raw reads from the Ant_24H, Ant_NC_24H and control group, respectively. Then, these raw data were further filtered to obtain 134891344, 137509356, 139274542 clean reads. All the databases have high quality, with clean base values ≥ 6.5 G, error rates ≤ 0.03, Q20 ≥ 97.62% and Q30 ≥ 93.66%. As the result, all of the libraries were determined to be suitable for future research. The RNA-Seq results are summarized in Table 1. Then, from each group, a total of 118712851, 122177868 and 124245564 clean reads which accounted for approximately 88%, 88.85% and 89.21% of the total clean reads were matched to the reference sequence (Table 2). Moreover, on average 93.46%, 93.49%, 92.72% of the clean reads were found in the exon regions of the genome from three groups, as shown in Fig. 2. Finally, in three groups, a total of 24387 genes were identified for additional quantitative analysis and functional annotation.

3.3 Differentially Expressed Genes Identification and Analysis

In contrast to the control group, there were 2914 DEGs as a result of miR-144-5p interference (P < 0.05), including 1828 significantly downregulated genes and 1086 significantly upregulated genes (Fig. 3A). Meanwhile, there were 2690 DEGs following antagomir-NC interference (P < 0.05), including 1397 significantly downregulated genes and 1293 significantly upregulated genes (Fig. 3B). According to the Venn diagram, 2047 transcripts were not affected in the Ant_NC_24H but changed in Ant_24H, accounting for about 39.3% of the identified genes (Fig. 3C). We utilized heatmap analyses of hierarchical clustering to identify the characteristics of the DEGs in this study. As shown in Fig. 3D, the heatmap analysis revealed that DEGs in the control and PomiR_Ant_NC_24H groups exhibited comparable clustering and expression patterns, whereas the Po-miR_Ant_24H group has distinct expression features. These findings suggest that miR-144-5p interference affected gene expression in flounder intestinal tissue, indicating that the impact of miR-144-5p on gene expression regulation should not be overlooked.

3.4 Analysis of miRNA-mRNA Prediction

In our study, a total of 2117 potential targeted genes were obtained. We selected 155 potential targeted genes which were stable in the Ant_NC_24H group but changed in Ant_ 24H group (Table 3), and further screened out 22 potential targeted genes in enriched immunity-related pathways. For example, apaf1 can participate in the regulation of various immune signaling pathways including p53 signaling pathway, Apoptosis and Herpes simplex virus 1 infection. LOC109641154 can be involved in the immune response of Focal adhesion, ErbB signaling pathway, C-type lectin receptor signaling pathway, MAPK signaling pathway and Salmonella infection, as shown in Table 3. These results suggested that miR-144-5p might regulate the immune response of flounder through targeting immunity-related genes.

3.5 Functional Annotation of DEGs

In order to investigate the potential role of DEGs, we performed GO enrichment and KEGG pathway analysis. DEGs, which were stable in the Ant_NC_24H group but changed in the Ant_24H group, were divided into three groups based on their GO functional enrichment analysis. A total of 2318 genes were matched to biological process (BP), 618 genes were matched to cellular component (CC), and 2084 genes were matched to molecular function (MF). The top 30 GO terms in these three categories are displayed in Fig. 4A. Most of the differential expression genes were enriched into cellular response to endogenous stimulus in biological process, endomembrane system in cellular component, and zinc ion binding in molecular function. Otherwise, several immunity-related GO terms, such as cell adhesion, regulation of cell death, regulation of apoptotic process, Wnt signaling pathway, cell-cell signaling by wnt, immune systemcomponents and immune response were significantly enriched in the biological process category (Table 4).

The KEGG enrichment analysis enabled us to get a greater comprehension of the biological functions of DEGs. In this study, 1355 DEGs among 2047 DEGs that were stable in the Ant_NC_24H group but changed in Ant_24H group were enriched in 147 signaling pathways. In order to depict the KEGG pathway enrichment data, the top 20 most significantly enriched pathways were shown in Fig. 4B. According to the KEGG pathway enrichment data, the DEGs have a substantial association with immunity-related signaling pathways, including NOD-like receptor, MAPK, Wnt, mTOR, Toll-like receptor signaling pathway. A summary of the 33 key immunity-related signaling pathways were showed in Table 5. In this study, function analysis of DEGs via GO and KEGG reveals that an abundance of immunity-related genes were regulated by miR-144-5p.

3.6 The Immunity-Related Protein-Protein Interaction Networks Construction

In vivo, proteins are the material underpinnings of physiological actions. The creation of immunity-related protein-protein interaction networks might aid in the discovery of key genes implicated in immunoreaction. Protein-protein function interactions were the core of regular physiological processes in P. olivaceus. The immunity-related protein-protein association networks have been constructed utilizing 67 DEGs from 33 immunity-related signaling pathways which were unaltered in the Ant_NC_24H but changed in Ant_24H (Table 5). As shown in Fig. 5, some nodes with several edges are thought to represent immune response hub genes and we opted to further investigate and validate these hub genes.

As a result, 5 core genes with multiple connections (degree ≥ 8) were discovered, including platelet-derived growth factor alpha polypeptide b (pdjfa), vinculin (vcl), vascular endothelial growth factor D (vegfd), C-X-C motif chemokine receptor 4 (cxcr4) and CF transmembrane conductance regulator (cftr). The top 20 DEGs most associated with immunity identified by WGCNA are presented in Table 6.

4 Discussion

In teleost fishes, miRNAs play a crucial role in immune cell differentiation, immune reaction, as well as the development of immunological disorders (Chu and Xu, 2020). The intestine of flounder represent an important organ for immunological response (Xiu et al., 2019). In order to research the immune regulation function of miR-144-5p in the intestine of flounder, we screened the differentially expressed genes by transcriptome sequencing after miR-144-5p interference. In this study, 2047 DEGs and 1823 DEGs were identified when comparing the fish in Ant_24H group with the control fish, and the fish in Ant_NC_24H group with the control fish, respectively. To further investigate the immunological response of miR-144-5p targeted DEGs, we selected 2047 genes that were unaffected in the Ant_ NC_24H group but changed in Ant_24H group. We discovered DEGs in the control and Ant_NC_24H groups had similar clustering and expression patterns, while the Ant_ 24H group has distinguishable expression features as a result of heatmap studies of hierarchical clustering. This result suggests that the immunological regulation of miR-144-5p should not be overlooked.

Based on the DEGs, GO and KEGG functional analyses display that there are many immunity-related genes in the intestinal transcriptome that respond to miR-144-5p interference, showing a complex immunological defensive response in P. olivaceus intestine. The thorough study of these GO enrichment and KEGG signaling pathways will help the research of the miR-144-5p immunoregulation in intestines of flounder, as well as give some significant hints for prevention and treatment of disease. We found 503 DEGs were highly enriched in 33 immunity-related pathways which were unaffected in the Ant_NC_24H group but changed in Ant_24H group.

Several signaling pathways were identified in our research, such as the NOD-like and Toll-like receptor signaling pathways, the phagosome, the intestinal immune network for IgT (IgA) production, as well as the p53 signaling pathway. In previous studies, miR-144 was discovered to be involved in the regulation of immune responses via regulating immunity-related signaling pathways. For example, in human, miR-144 promotes IL-1β-induced chondrocyte pyroptosis by regulating the NOD-like signaling pathway (Jiang et al., 2021). Additionally, miR-144-5p inhibits cancer cells proliferation by regulating the Toll-like signaling pathway via targeting Toll-like receptor 2 (TLR2) (Liu et al., 2020a). Moreover, miR-144 induced apoptosis and autophagy in lung cancer cells via regulating P53 signaling pathway by targeting TIGAR (Chen et al., 2015). In mouse, miR-144-5p agomirs mitigated M1 macrophage-associated inflammation by regulating the NF-κB signaling pathway via targeting TLR2 (Shi et al., 2020). In fish, miR-144 play regulatory roles in inflammatory responses after bacterial infection by regulating the NF-κB signaling via targeting NOD1 (Chu et al., 2021) and IL-1β (Yan et al., 2019) in miiuy croaker. Overall, we suggested that the miR-144 not only regulate cancer cells proliferation but also plays several roles including antimicrobial function and regulation of inflammatory responses in immune regulation.

In the miRNA-mediated immune response, functional connections between miRNA and genes are crucial. We discovered dynamic variations in the expression of a wide variety of immunity-related genes. On the basis of the 67 immunity-related DEGs, we constructed protein-protein interaction networks. As shown in Fig. 4, 67 immunity-related DEGs yielded 5 hypothesized hub genes, each of which was engaged in several protein-protein connections or KEGG pathways.

Depending on protein-protein association data, we determined that pdgfa, vegfd, and vcl are the three top hub genes (Fig. 4). Pdgfa was involved in cell survival, immunity and wound healing via regulating focal adhesion by interacting with vcl and vegfd (Wang et al., 2019). In our study, pdgfa, vegfd, and vcl were significantly down-regulated 24 h after miR-144-5p interference. Similar to previous observations, miR-93, miR-155 and miR-10b inhibitors suppress the expression of pdgfa, vegfd and vcl, respectively (Ibrahim et al., 2012; Chang et al., 2019; Feng et al., 2019). In a recent research, pdgfa connected with proinflammatory cytokines (TNF-α, IL-6, IL-1, IL-8) and chemokines (CXCL-1 and CCL2) was regulated by miR-29a (Wang et al., 2019). In zebrafish, vegfd is involved in angiogenesis and lymphangiogenesis through regulation by miR-194-5p (Song et al., 2007). Otherwise, deficiency of miR-663b-mediated vcl can promote epicardial hyperplasia and coronary vessel disruption (You et al., 2021). Our findings suggest that in fish, miR-144-5p can regulate the expressions of pdgfa, vegfd, and vcl which might play an important role to mitigate inflammatory responses against invading pathogens. Other hub genes, such as cftr and cxcr4, can regulate embryonic organizer formation in zebrafish (Liu et al., 2020c) and activate chemokine signalling pathway in immune responses to infection in channel catfish (Yeh and Klesius, 2010), respectively. Furthermore, miR-200b participates in the immunoreaction of hosts by regulating Wnt pathway via targeting cftr (Liu et al., 2020c) and TFA (Xu et al., 2018). Cxcr4 can regulate leukocyte chemotaxis, inflammatory response and numerous autoimmune diseases through the regulation by miR-23a (Pan et al., 2018) and participate in the intestinal immune network for IgT (IgA) production pathway (Mousavi, 2020). The findings in this research will be utilized to further research the miRNA-mediated immunological response to bacterial infections in flounder. The regulation function of miRNA need to be further studied in immune system of the flounder in the future.

5 Conclusions

In our study, RNA-Seq was utilized to detect transcriptional expression in the intestine of flounder in response to miR-144-5p interference. Primarily, we detected 2047 DEGs which were unaffected in the Ant_NC_24H group but changed in Ant_24H group. Then, 503 DEGs were identified to be enriched in 33 immunity-related pathways, demonstrating a broad immune response mediated by miR-144-5p interference. Following that analysis, a protein association interaction network was constructed, and 5 immunity-related hub genes were discovered to be engaged in immune regulation. Our findings emphasize the regulation roles of miR-144-5p in the flounder immune reaction on inflammatory response, antibacterial function and autoimmune diseases.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 32002421), the Advanced Talents Foundation of QAU (No. 6651118016), the Natural Science Foundation of Shandong Province (No. ZR20 19BC009), the Fish Innovation Team of Shandong Agriculture Research System (No. SDAIT-12-06), the Aquatic Animal Immunologic Agents Engineering Research Center of Shandong Province, Key Research and Development Program in Shandong (No. 2018YFJH0703), the Consulting Research Project of Shandong Research Institute, China Engineering Science and Technology Development Strategy (No. 202103SDYB08), the 'First Class Fishery Discipline' Programme in Shandong Province, and the special top talent plan 'One Thing One Decision (YishiYiyi)' Programme in Shandong Province, China.