Integrative Transcriptomic and Small RNA Sequencing Reveals Immune-Related miRNA–mRNA Regulation Network for Soybean Meal-Induced Enteritis in Hybrid Grouper, Epinephelus fuscoguttatus♀ × Epinephelus lanceolatus♂

A 10-week feeding experiment was conducted to reveal the immune mechanism for soybean meal-induced enteritis (SBMIE) in hybrid grouper, Epinephelus fuscoguttatus ♀ × Epinephelus lanceolatus ♂. Four isonitrogenous and isolipidic diets were formulated by replacing 0, 10, 30, and 50% fish meal protein with soybean meal (namely FM, SBM10, SBM30, and SBM50, respectively). The weight gain rate of the SBM50 group was significantly lower than those of the other groups. Plica height, muscular layer thickness, and goblet cells of the distal intestine in the SBM50 group were much lower than those in the FM group. The intestinal transcriptomic data, including the transcriptome and miRNAome, showed that a total of 6,390 differentially expressed genes (DEGs) and 92 DEmiRNAs were identified in the SBM50 and FM groups. DEmiRNAs (10 known and 1 novel miRNAs) and their DE target genes were involved in immune-related phagosome, natural killer cell-mediated cytotoxicity, Fc gamma R-mediated phagocytosis, and the intestinal immune network for IgA production pathways. Our study is the first to offer transcriptomic and small RNA profiling for SBMIE in hybrid grouper. Our findings offer important insights for the understanding of the RNA profile and further elucidation of the underlying molecular immune mechanism for SBMIE in carnivorous fish.


INTRODUCTION
In recent years, the global contribution of fish meal to aquafeeds has sharply reduced (1). As a result, soybean meal, which is considered one of the most hopeful candidates for fish meal replacement, can partially or fully replace the fish meal but introduces many anti-nutritional factors. Consequently, fish enteritis induced by plant proteins has become one of the main challenges for sustainable aquaculture (2), and this occurs in a dose-dependent manner (1,3).
Soybean meal-induced enteritis (SBMIE) has been found in many commercial fish species such as Atlantic salmon (Salmo salar) (4), grass carp (Ctenopharyngodon idella) (5), and turbot (Scophthalmus maximus L.) (6,7), and mainly occurs in carnivorous fish (4,6,7). The symptoms of SBMIE are most apparent in the posterior/distal intestine of turbot (6,7) and grass carp (5), the most important mucosal immune organ (8). Fish SBMIE has been found to be accompanied by a decrease in the height of villi and microvilli, downregulation of tight junction protein claudin-4, occluding, and ZO-1 mRNA levels, and upregulation of pro-inflammatory cytokine genes, including TNF-α, IL-1β, IL-8, and IL-16 (9)(10)(11)(12), influencing the innate immune response. Zebrafish SBMIE is T celldependent and has a T helper (Th) 17 cytokine profile (13). What are the underlying immune mechanisms in fish SBMIE? The new omics technologies, including genomics, proteomics, and transcriptomics, have great potential for investigating and explaining the complex relationship between fish nutrition and immunity, both in intestine health and disease (14). At the genomic level, most components associated with T lymphocyte function have been identified in fish, suggesting that gut-associated lymphoid tissue has similar functionalities between fish and mammalian T lymphocytes (15). At the transcriptomic level, immune-related pathways of fish SBMIE have been gradually reported, showing that cytokinecytokine receptor interaction, NOD-like receptor interaction, the intestinal immune network for IgA production, and the NF-kB, Jak-STAT, T-cell receptor, and TNF signaling pathways played key roles in response to SBM stress (3,16). Fish meal replacement in fish diet by alternative protein sources could change the fish intestine proteome, including innate immune proteins (17).
The especially interesting fact is that miRNAs are involved in regulating intestine function, including epithelial cell growth (18), mucosal barrier function (19), and the development of gastroenteric diseases (20)(21)(22). An important aspect is that miRNA can also regulate mRNA expression in fish at a transcriptional level. Recently, studies of miRNAome in turbot intestinal function have reported that miRNAs contributed to the intestinal immune responses, preventing host infection, which the potential target genes of differentially expressed miRNAs were involved in multiple functional categories, including the RIG-I signaling pathway, immune defense/evasion, the toll-like receptor signaling pathway, and inflammatory responses (23). Also, it was found via small RNA sequencing that fish diet could affect the expression of intestinal miRNAs and target genes and immune-related pathways, including cell adhesion molecules, ECM-receptor interaction, the apoptosis signaling pathway, cytokine-cytokine receptor, and the VEGF signaling pathway (3,24). However, there has been a lack of investigation of the underlying immune response by combining both transcriptomes and small RNA sequencing, and this requires further elucidation, especially in carnivorous fish.
Hybrid grouper (Epinephelus fuscoguttatus♀ × Epinephelus lanceolatus♂) is a carnivorous fish species that is the main farmed species in China due to its outstanding delicious taste and better growth rate and survival compared with the broodfish. Information on how nutrition influences the intestine health of hybrid grouper has been gradually accumulated in recent years (25)(26)(27)(28)(29)(30). Our previous research used metabolomics technology to identify 17 potential markers of SBMIE in hybrid grouper (31). However, this study aims to reveal the immunerelated miRNA-mRNA regulation network for SBMIE in hybrid grouper by integrative transcriptome and small RNA profiling from the perspective of molecular immunology and to provide another important insight to further solve the problem of fish intestinal health.

Experimental Diets
The use of hybrid grouper juveniles was approved by the Animal Research and Ethics Committees of Guangdong Ocean University, China. Four isonitrogenous and isolipidic diets were formulated to contain 0, 7.41, 22.24, and 37.07% of soybean meal (SBM) by replacing 0% (FM, control), 10% (SBM10), 30% (SBM30), and 50% (SBM50) of fish meal (FM) protein, respectively. The formulation of the basic experimental feeds is presented in Table 1. All ingredients were systematically mixed with lipid sources such as fish oil, soybean oil, and soybean lecithin and then purified water was added to produce a homogenous mixture. The dough was pelleted through a double helix extrusion mechanism (F-75, South China University of Technology, China). Feeds (2.5 mm diameter) were air-dried and then stored at −20 • C until feeding.

Feeding Trial
Hybrid grouper Juveniles were obtained from a native species farm (Zhanjiang, China). All fish were adapted under the feeding system for 2 weeks by feeding with a commercial diet. Uniformly sized fish (mean initial weight ± SE = 17.01 ± 0.04 g) were randomly divided into four groups in triplicate, with 30 individuals in each fiberglass tank (300 L). The fish were slowly fed twice a day at 08:00 and 17:00 for 10 weeks. During the experiment, the water temperature fluctuated from 28 to 30 • C, the dissolved oxygen concentration was kept at >7 mg/L, and ammonia and nitrate were kept at <0.03 mg/L.

Sample Collection
Before the termination of the 10-week feeding experiment, the fish per tank were fasted for 24 h before collecting samples and were then counted and weighed to determine growth indexes, including weight gain rate, feed conversion ratio, and survival rate. After weighing, distal intestines of two fish per tank were collected and instantly transferred to 4% paraformaldehyde solution for histological examination. At the same time, distal intestines of another three fish per tank were collected as a single sample and instantly frozen in liquid nitrogen, then stored at −80 • C for RNA extraction. Based on growth performance (see section Growth Performance) and histological examination, small RNA and transcriptome analyses were performed on distal intestine samples of the FM and SBM50 groups to ensure the maximum difference between samples. Thereby, the probability of detecting differential expression was increased between samples.

Intestinal Morphology
The fixed distal intestine samples from FM and SBM50 groups were dehydrated in a series of graded ethanol and embedded in paraffin. Distal intestine sections (7 µm thick) from each sample were cut and then stained with hematoxylin/eosin. The sections were observed under an inverted microscope (Nikon, Japan), and 10 plicas and muscle layer thicknesses (MLT) were randomly selected per slice. Plica height (PH), plica width (PW), MLT, and the number of goblet cells (GC) per slice were measured using the image acquisition software (NIS Elements, version 4.60, Nikon, Japan).

Transcriptome Sequencing and de novo Assembly
One microgram total RNA from the FM and SBM50 treatment groups was used for transcriptome library preparation. Total RNA was purified by beads containing oligo (dT). First-strand cDNA was then generated in a First-Strand Reaction System by PCR, and the second-strand cDNA was generated as well. The cDNA fragments with adapters were amplified by PCR, and the products were purified using AMPure XP Beads. The library was validated on the Agilent Technologies 2100 bioanalyzer for quality control. Transcriptome sequencing was carried out on a BGISEQ-500 platform (BGI-Shenzhen, China). Trinity (32) was used to achieve de novo assembly with clean reads, and Tgicl (33) was then used to cluster transcripts to Unigenes. The expected number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) was used to calculate the mRNA gene expression. Differentially expressed genes (DEGs) of two groups (SBM50 and FM) were quantified using the cutoff |log2FC| >1, P < 0.001 by DESeq R package (34,35). After assembly, All-Unigenes were searched and annotated against the publicly available protein databases, including Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), KOG (EuKaryotic Orthologous Groups), Swiss-Prot, and GO (Gene Ontology). The pathway assignments were performed by sequence searches against the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. KEGG terms with corrected P-values (Q-values) ≤ 0.05 were considered significant. Transcriptome (de novo assembly) sequencing data were deposited into the NCBI SRA database with the accession number SUB7020170.

Construction and Sequencing of Small RNA Libraries
Six small RNA libraries were constructed from the FM and SBM50 treatment groups. Total RNA extraction was performed from the distal intestine using Trizol Reagent (Invitrogen, USA). Subsequently, 1 µg of total RNA per sample was used for small RNA sequencing. The quality of RNA samples was evaluated using the Agilent 2100 Bioanalyzer. Small RNA fractions were ligated to 3 ′ and 5 ′ adapter. Quantitative reverse transcription PCR (RT-PCR) was carried out on the adaptor-ligated small RNAs. PCR products were purified by QIAquick Gel Extraction Kit (Qiagen, Germany) and used for sequencing on the BGISEQ-500 platform (BGI-Shenzhen, China). Clean reads were obtained by cleaning low-quality tags, removing adapter sequences, and filtering adaptor-ligated contaminants and sequences fewer than 18 nucleotides (nt). The final reads were mapped to the Hypoplectrus puella (GCA_900610375.1) reference genome by Bowtie2 (36). Clean reads were compared against small RNAs (rRNA, scRNA, snoRNA, snRNA, tRNA, and mRNA) using the Rfam database to annotate small RNA sequences. Finally, the miRBase20.0 was used to look for know miRNA. The hairpin structures were used to predict novel miRNAs using miRDeep2 software (37). miRNA expression levels were compared between the SBM50 and FM groups to identify differentially expressed miRNAs (DEmiRNAs). Firstly, data were normalized to obtain transcripts per million (TPM) values using the following formula: normalized expression = actual miRNA count/(total reads) ×1,000,000 (38). Fold-change values were then calculated based on log 2 (SBM50/FM) expression. The corrected P-value corresponds to the differential gene expression test using the Bonferroni method (39). Differential miRNA expression between two groups was analyzed with DESeq software based on the following thresholds: P-value ≤ 0.01 and |log 2 ratio| ≥ 1. RNAhybrid (40) and miRanda (41) software predicted the potential target genes of miRNA candidates, as described previously elsewhere (42,43). The DAVID gene annotation tool was used for the KEGG pathway annotation of predicted miRNA targets. Small RNA sequencing data were deposited into the NCBI SRA database with the accession number SUB7175134.

Network Analysis of DEmiRNA and DEG Interaction
Pearson's correlation coefficients between DEmiRNAs and their target genes were calculated using the correlation function in RStudio. To obtain the positive and negative correlations between two groups, the potential target genes of miRNAs were overlapped with the identified upregulated or downregulated DEGs, respectively. All of the relationship pairs between DEmiRNAs and their DE target genes were used to construct the interaction network using Cytoscape v.3.7.2 software.

Validation by Real-Time Quantitative PCR (RT-qPCR)
RT-qPCR validation was carried out on the same samples used for transcriptome sequencing (n = 3). Primers were designed from the candidate gene sequences by premier 5.0 software and the online Primer-BLAST program. Primers used in this study are provided in Table S1. One microgram total RNA for RNA sequencing was reverse transcribed into cDNA. Real-time PCR assays were conducted on a CFX96 real-time PCR Detection System (Bio-Rad, Hercules, CA) with 5 µL SYBR Green Master Mix (Takara, China). β-actin was selected as the reference gene according to a previous study (26). The small RNA of the same samples used for sequencing was extracted using an RNAiso for Small RNA Kit (Takara, China) according to the manufacturer's protocol. Subsequently, the first-strand cDNA was synthesized for mature miRNA expression analysis by a Mir-X TM miRNA First-Strand Synthesis Kit (Code No. 638315, Takara, China). The qPCR was carried out using a miRNA SYBR Green RT-qPCR Kit (Takara, China) with the provided miRNA reference gene (U6). Relative quantitative levels were calculated based on the 2 − CT method (44).

Statistical Analysis
The normal distribution and the homogeneity of the variance of growth indexes were tested, followed by one-way analysis of variance and Tukey s test. Morphological analysis between two groups was assessed by two-tailed unpaired Student's t-test (GraphPad). For statistically significant differences, P < 0.05 was required. All statistical analyses were carried out using SPSS 24.0 software. The barplot was generated by Graphpad Prism 8.0.1 software.

Growth Performance
The survival rate (SR) was not affected by dietary treatment levels (P > 0.05, Table 2). The weight gain rate (WGR) of the SBM50 group was significantly lower than that of other groups (P < 0.05). The feed conversion ratio (FCR) of the SBM50 group was significantly higher than that of the other groups (P < 0.05).

Histological Examinations and Intestinal Morphometry
More swelling of the lamina propria (LP) of the SBM50 group was observed compared with the FM group, and the intestinal villi of the SBM50 group showed signs of shedding (Figures 1A,B). The plica height (PH), muscular layer thickness (MLT), and number of goblet cells (GC) of the SBM50 group were much lower than those of the FM group ( Figure 1C, P < 0.01). There was no significant difference in plica width (PW) between the FM and SBM50 groups (P > 0.05).

Analysis of mRNA Sequencing
A total of six qualified libraries from the FM and SBM50 groups, with three biological replicates per treatment, were sequenced. An overview of the sequencing and assembly data is presented in Table 3. Approximately 33.49 and 33.6 Gb of clean reads were obtained in the FM and SBM50 groups. More than 87.1% of the reads had Q-scores at the Q30 level, and more than 75.8% of the clean reads were aligned. The length distribution of the Unigene in all six libraries is shown in Figure S1.

Analysis of Small RNA Sequencing
Six small RNA libraries were constructed, with three biological replicates per treatment. A total of 85,789,507 and 87,445,155 raw reads were found in the FM and SBM50 groups, respectively ( Table 4). Also, 77,115,163 and 81,360,125 clean reads were found in the FM and SBM50 groups, respectively. The total mapped tags to the reference genome in the FM and SBM50 groups were 71,413,180 (92.60%) and 76,824,617 (94.43%), respectively. Most small RNAs were 21-23 nt in length in all six libraries, with 22 nt being the most frequent length (Figure S2), and more than 61.7% were miRNA in the catalog of small RNA in all six libraries ( Figure S3). A total of 682 mature miRNAs (Table S3) and 29 novel miRNAs (Table S4) were identified in these six small RNA libraries.

Integration Analysis of the DEmiRNAs and DEGs
A total of 244 miRNA-mRNA interactions were identified in the FM and SBM50 groups, with the involvement of 92 DEmiRNAs and 211 DEGs (Figure 6, Table S7). A positively correlated expression pattern was seen for 180 mRNA-miRNA pairs. Most miRNAs had multiple possible target genes, while different miRNAs could regulate the same target. For instance, miR-196 was the regulator of CL4654.Contig2_All, CL7786.Contig3_All, and Unigene16944_All, whereas miR-20a-5p and miR-459-5p_1 could regulate the expression of CL8081.Contig7_All.

KEGG Enrichment Analysis of Differential Target Genes
Differential target genes were predicted and annotated to 241 pathways. The top 5 of the top 20 KEGG pathways were phagosome, tuberculosis, osteoclast differentiation, natural killer cell-mediated cytotoxicity, and mineral absorption (ko04978), respectively (Figure 7, Table S8).
To further explain the possible immune response for SBMIE, combined with the KEGG enrichment pathway of DEGs found that 15 miRNA-mRNA pairs involved in the phagosome, 8 miRNA-mRNA pairs involved in natural killer cell-mediated    Table 5.

Real-Time Quantitative PCR Validation
In total, 10 miRNAs and 10 mRNAs were selected to test their expression and the results suggested that DEmiRNAs except for miR-194-5p and DEGs showed a similar expression pattern to the high-throughput sequencing data (Figure 8).

DISCUSSION
The nutritional trial of this study revealed that a 22.4% SBM, i.e., 30% substitution of FM protein, did not significantly influence the growth of hybrid grouper compared to the FM diet. However, the effect of SBM on growth was substitution-related, and the growth performance became progressively compromised with increasing SBM substitution level, which is consistent with the previous studies (45,46). In this study, growth performance was significantly reduced when dietary SBM content reached 370 g/kg, i.e., 50% substitution of the FM protein. Hybrid grouper fed with SBM replacing 50% of the FM protein diet showed swelling of the lamina propria and reduction of plica height, muscular layer thickness, and number of goblet cells in the distal intestine of hybrid grouper, indicating that enteritis and intestinal injury had appeared. The severity of histopathological changes observed under SBM application relies on the level of soybean inclusion. SBMIE is now commonly used as a model for the study of intestinal inflammation in fish (6,7,47). The histopathological changes of SBMIE have been widely researched and are characterized by a swelling of the subepithelial mucosa and lamina propria, a reduced mucosal fold height, a profound infiltration of various inflammatory cells, and loss of normal enterocyte supranuclear absorptive vacuolization (13,31,48). Liu et al. (10) found that turbot fed with SBM replacing 40% of the fish meal protein diet showed obvious enteropathy, including a reduction of absorptive surface and obvious infiltration of mixed leukocytes in the lamina propria. To reveal the underlying immune response of fish SBMIE, 6,390 DEGs in the distal intestine of hybrid grouper were identified in this study. Also, the enhanced gene expressions were found to involved in immune-and inflammatory-related pathways, including phagosome, natural killer cell-mediated cytotoxicity, the intestinal immune network for IgA production, Fc gamma R-mediated phagocytosis, and the NF-kappa B signaling pathway. The above results suggested that these pathways played a vital role in fish SBMIE. Similar immunerelated pathways have been reported in SBMIE in other carnivorous fish, such as salmon (16,49) and turbot (6,7). It is worth noting that the downregulated gene expression was in KEGG pathways mainly involved in lipid metabolisms such as biosynthesis of unsaturated fatty acids, linoleic acid metabolism, cholesterol metabolism, fat digestion, absorption, and arachidonic acid metabolism. This may indicate that impaired lipid metabolism could be a consequence of "tissue malfunction" (47,49). The results enable a better understanding of why LC-PUFA biosynthesis, cholesterol biosynthesis, lipid digestion, and the PPAR signaling pathway of distal intestine were influenced when Atlantic salmon ingested feed in which the fish meal was partially replaced by soybean meal (1,16).
Despite the major contribution of mRNAs, miRNAs also play a key role during immune processes in fish SBMIE. In this study, predicted target genes for DEmiRNAs are annotated to 340 signaling pathways. Immune-related pathways, including ECM-receptor interaction, the NF-kappa B signaling pathway, and the IL-17 signaling pathway, were enriched. The downregulation of miR-192-3p and miR-212-5p expression involved in the regulation of the ECM-receptor interaction pathway. The down-regulation of genes involved in the ECMreceptor interaction pathway in response to SBM stress was also reported in Grass carp (5). To further introduce miRNAs in detail, enteritis-related miRNAs have been found both in human and mammals. The present miRNA results include inflammatory bowel diseases (IBD) related-miRNAs, such as miR-124, miR-24, miR-221, and miR-132. The results of this study showed that upregulation of miR-124 expression and downregulation of miR-24, miR-221, and miR-132 expression were observed in the SBM50 group (Table S5). Similar expression patterns were found in colon tissues of children with active ulcerative colitis (UC), where decreased levels of miR-124 appeared to enhance expression and activity of STAT3, which could induce inflammation and pathogenesis (50). There were elevated levels of miR-24, miR-221, and miR-132 in colonic biopsies from UC, which suggested that they are an important regulator of the intestine barrier that may be essential in the pathogenesis of IBD (51). The above results indicated that miR-124, miR-24, miR-221, and miR-132 may play important roles in SBMIE of hybrid grouper juveniles, which may suggest their therapeutic potential.
To further explain the possible immune response for SBMIE through integrative transcriptomic and small RNA sequencing, the KEGG enrichment pathways of differential target genes were analyzed. Combined with the KEGG enrichment pathways of DEGs, it was concluded that immune-related signaling pathways such as phagosome, natural killer cell-mediated cytotoxicity, Fc gamma R-mediated phagocytosis, and the intestinal immune network for IgA production were enriched. The most enhanced gene expression being in the phagosome pathway suggested the involvement of a macrophage as the main intestine phagocyte during enteritis (52). A total of 14 miRNA-mRNA pairs also suggested that the phagosome pathway could play a key role in intestinal inflammation. In addition, this study identified 92 DEGs related to the natural killer cell-cytotoxicity pathway. The role of this pathway in the immune response to pathogens has also been reported in different fish species such as large yellow croaker (Larimichthys crocea) (53) and half-smooth tongue sole (Cynoglossus semilaevis) (54). A previous study reported that mammalian natural killer (NK) cells mediated cytotoxic activity via two distinct pathways (55). NK cells can release cytotoxic granules, including perforin and granzymes, on the surface of diseased cells. Granzymes can then stimulate caspase activation, mitochondrial dysfunction, or apoptosis (55). In this study, the upregulation of granzyme (gene id: Unigene24888_All) and perforin (gene id: Unigene31830_All, Blast Nr annotated in Table 5) gene expressions indicated that granule-mediated cytotoxicity may be triggered by the targeted release of lytic granules toward a locally attached target cell. Wu et al. (5) reported that the intestinal immune network for the IgA production pathway was upregulated in the early stages of Grass carp in response to high-SBM-content stress. Similar results were also reported in the present study: most genes in pathways of Fc gamma receptormediated phagocytosis and the intestinal immune network for IgA production were significantly upregulated. The role of the two pathways in immune response to pathogens has also been reported in fish such as half-smooth tongue sole (53), large yellow croaker (54), and darkbarbel catfish (Pelteobagrus vachellii) (56). Previous studies reported that the target genes of the DEmiRNAs were involved in pathways of Fcγ R-mediated phagocytosis and the intestinal immune network for IgA production (57,58). In this study, related miRNA-mRNA pairs were also enriched in the two pathways mentioned above (Table 5). Thus, our results suggested that phagosome, natural killer cell-mediated cytotoxicity, Fc gamma R-mediated phagocytosis, and the intestinal immune network for IgA production pathways may play a vital role in SBMIE of carnivorous fish.
In conclusion, our study is the first to offer the transcriptomic and small RNA profiles for SBMIE in hybrid grouper. Overall, 6,390 mRNAs and 92 miRNAs were differentially expressed under dietary SBM stress. Our findings support the notion that DEmiRNAs and their target mRNAs play an important role in immune regulation. Also, investigation of KEGG enrichment pathways by integrative transcriptomic and small RNA profiling revealed that the immune mechanism for SBMIE in hybrid grouper may be associated with the phagosome, natural killer cell-mediated cytotoxicity, Fc gamma R-mediated phagocytosis, and the intestinal immune network for IgA production pathways. Our findings offer important insights for the understanding of the RNA profiles and further elucidation of the underlying molecular immune mechanism for SBMIE in carnivorous fish.