Comparative Transcriptomic and Metagenomic Analysis Reveal the Key Factors on the Growth Rate in Red Swamp Craysh (Procambarus Clarkii)

Background: Although red swamp craysh (Procambarus clarkii) is one of the most important species in aquaculture, the factors that inuence differences in growth and development between individual siblings are still not fully understood. To address this lack of knowledge, we designed experiments to elucidate factors contributing to individual differences by comparing the hepatopancreatic transcriptome and gut ora structure of individuals with differences in sibling craysh size under the same rearing conditions and attempted to nd links between gut ora and host transcriptome information. Result: In total, 300691028 high-quality reads were obtained that were used to assemble 60637958 unigenes. Comparison of the expression proles of the hepatopancreas of craysh of different sizes revealed 497 differentially expressed genes (P<0.05). A total of 32 KEGG signaling pathways were found to be enriched after the KEGG database and GO functional annotation analyses, with the highest number of unigene enrichments related to organismal metabolism. Additionally, we found that proteobacteria, tenericutes, actinobacteria, and bacteroidetes made up a larger portion of the microbiome of larger individuals, suggesting that craysh microbiota adapts to rapid growth, which may promote accelerated development by regulating the expression level of relevant genes. Conclusions: This work has accumulated data to support the potential impact of structural alterations in gut ora under equivalent feeding conditions in explaining the differential expression of transcriptomic genes among differentially developing individuals of sibling craysh. These results further elucidate how the intestinal environment affects the rapid development of invertebrate crustaceans and provide a reference for further understanding of the regulatory mechanisms of the host's in vivo environment.


Background
The red swamp cray sh (Procambarus clarkii) is native to the south-central United States (primarily in Louisiana) and northeastern Mexico, but has been transplanted around the world [1][2][3]. P. clarkii is a freshwater cray sh, which is economically signi cant due to its culinary uses [4,5]. Additionally, it is omnivorous and has a rapid growth rate and high environmental adaptability [6]. In recent years, P. clarkii has become an important economic farming species in China.
The intestinal tract is the most important digestive and absorption organ of crustaceans, including shrimps and crabs, and hosts a large number of complex microorganisms that are interdependent and mutually constrained by their hosts, forming a unique intestinal microecosystem [7]. The relationship between gut microbes and their hosts has been shown to play a crucial role in growth rate and disease susceptibility in several different species, due to its impact on nutrient processing, energy balance, as well as growth and development in the host [8]. Shrimp and cray sh lack an acquired immune system and rely on various types of non-speci c immune factors to recognize foreign objects and pathogens [9]. The structure and function of the normal ora in the intestine of shrimp and sh has been shown to have a signi cant impact on their immune systems due to its ability to form a powerful barrier against pathogenic bacteria [10][11][12].
The hepatopancreas is an important multifunctional organ of crustaceans, which plays a role in metabolism, nutrient absorption, and immune function [13]. Studies have also found that changes in the hepatopancreas of cray sh can be caused by the accumulation of heavy metals and affect the ability of species to cope with various environmental stresses [14]. Recent studies have sequenced and analyzed the transcriptomes of ocular stalk, hepatopancreas, muscle, gills, intestine, blood cells, and stomach of cray sh in order to better understand environmental adaptation and response to disease, but there is a lack of transcriptomic data involving the role of intestinal ora in growth and development among siblings [15,16]. To ll this current knowledge gap, this study employed transcriptomic data to analyze the role that the gut ora and the host transcriptome play in affecting the growth rate of P. clarkii. These studies further elucidate how the intestinal environment affects the rapid development of invertebrate crustaceans, providing insights to further understand the regulatory mechanisms of the host's internal environment and to understand the safety of the host's internal environmental homeostasis.

Results
Transcriptome sequencing and assembly Hepatopancreatic mRNA from differentially developing sibling individuals was detected by RNA-Seq, with three replicates per group. Members of the smaller group were termed SI1, SI2, and SI3, while members of the larger group were termed LI1, LI2, and LI3. After the removal of low-quality reads, clean reads were assembled into reference sequences for subsequent analyses using Trinity [18]. The LI group sequencing reads produced 158810616 unigenes, while the SI group yielded 141880412 unigenes. The error rate of the samples was restricted between 0.02 and 0.03, with Q20 accounting for more than 97% (Table   1). The reference sequence set included 135061303 transcripts and 60637958 unigenes, with unigene (gene) N50 values of 2468 and N90 values of 511 in the samples, respectively (Table 2).

Identi cation of DEGs
After quanti cation of gene expression levels, differential expression analysis was performed using screening criteria (| log2 (fold change) |> 1 and adjusted pvalue <0.05) to identify genes with signi cant differences between different samples. In brief, 497 DEGs were found in the LI vs SI group, of which 221 were up-regulated and 276 were down-regulated (Fig 2A). The analysis of the signi cance of the DEGs is shown in Figure 2B. Signi cantly DEGs were plotted into clustered heat maps based on FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced)) values, as shown in Figure 3.
The annotated genes were then categorized into ve branches according to KEGG metabolic pathways (level 2): Cellular Processes, Environmental Information Processing, Genetic Information Processing, Genomic Information Processing, Metabolism and Organic systems. In cellular processes, there were 398 unigenes in transport and catabolism and 229 unigenes in cell community signaling. In environmental information processing, the number of unigenes in signal transduction (601) was greater than the number in signaling molecules and interacting transport. The number of unigenes annotated in translation, folding, sorting, and degradation based on information processing was 320, and the number of unigenes in lipid metabolism and carbohydrate metabolism was approximately 250. Among the ve branches of KEGG, the metabolic pathways were the most numerous. Among the immune signaling pathways, the endocrine system (368) and the immune system (251) were found in the most and second most unigenes, respectively. The organic system ranked second among the ve branches. Both metabolic and organic systems were closely related to organismal growth (Fig. 4).

Validation of RNA-Seq with RT-qPCR
To further validate the accuracy of our RNA-Seq data, we selected 13 genes for RT-qPCR (Tables S1 and S2). These genes included members involved in apoptosis, immunity and organismal metabolism. Figure 5 shows that all RNA sequences are consistent with the RT-qPCR results. Therefore, the RNA-Seq data were considered to be accurate and reliable in this study. The RT-qPCR results also showed that there were signi cant differences in the expression of genes associated with metabolism and immunity across different individuals.
Taxonomic composition of intestine microbiota At the end of the feeding cycle, high-throughput sequence analysis of bacterial 16S ribosomal RNA v3-v4 region was performed on the intestinal contents of P. clarkii. As shown in Supplemental Table S3, the top four most abundant bacteria in the LI group were Proteobacteria (49.22%), Tenericutes (4.66%), Actinobacteria (3.11%), and Bacteroidetes (2.97%). In the SI group, Proteobacteria, Firmicutes, Actinobacteria, and Chloro exi were the four most abundant phyla, accounting for 39.57%, 5.83%, 5.35%, and 5.34%, respectively.
Based on the sequence information available for each sample, clustering was performed at both the species and sample levels to create heat maps ( Fig 6A).
To further investigate the phylogenetic relationships among bacterial genera, representative sequences of the top 100 genera were compared by multiple sequence alignment (Fig. 7). The species richness of Vogesella, Bacillus, Aeromonas, Bacteroides, De uviicoccus, Acinetobacter, Gemmobacter, and Rhodobacter were further compared at the genus level. The results showed signi cant differences (P<0.05) in LI vs SI in Gemmobacter, Acinetobacter, Candidatus-Bacilloplasma, Vogesella, Rhodobacter, Paracoccus, and Bosea.

Microbial Function Prediction and Functional Enrichment of the Hepatopancreas Transcriptome
Annotation results of gut microbial function showed that cray sh appeared to have a higher proportion of KEGG enrichment in the environmental information and metabolism categories (Fig. 8A). In the metabolism category, carbohydrate metabolism and amino acid metabolism were the most signi cantly enriched categories. Membrane transport was found to be most abundant in the environmental information processing category.
A scatterplot was used to generate a graphical representation of the KEGG enrichment of the transcriptome DEGs from previous analyses, with the 20 most enriched pathway entries displayed (Fig. 8B). This list included fatty acid biosynthesis, phototransduction, prion diseases, ribo avin metabolism, fructose, mannose metabolism, ether lipid metabolism, tyrosine metabolism, propanoate metabolism, PPAR signaling pathway, glycolysis/gluconeogenesis, cysteine and methionine metabolism, amyotrophic lateral sclerosis (ALS), drug metabolism, adipocytokine signaling pathway, melanogenesis, citrate cycle (TCA cycle), AMPK signaling pathway, vascular smooth muscle contraction, insulin signaling pathway, and lysosome. The majority of these terms were associated with nutritional metabolism or immune processes FigureS1.

Intestinal biopsy analysis and humoral factor testing
The results of the humoral immune factor assay showed that the activity of SOD was signi cantly higher in the LI group than in the SI group (P<0.05), and the activity of LI Ppo and Lyz were higher than in the SI group ( Fig 9A). Gls, AMS, and lipase in body uids of the LI group was signi cantly more active than those of the SI group (P<0.05) (Fig 9B).
The developmental differences in the intestinal tissue sections of the P. clarkii (Fig 10) revealed that the mucosal folds of the LI group were intact (Fig10A), with a higher height, tighter arrangement, larger folds, and epithelial cells in the mucosal folds arranged in an orderly manner with smooth striated edges on the surface (Fig10B); the intestinal villi of the SI group were sparsely arranged and short, and the depth of the intestinal crypt was deepened.

Discussion
In this study, we used RNA-Seq and 16S ribosomal RNA-Seq to analyze the hepatopancreatic transcriptome and the intestinal micro ora structure of P. clarkii, to identify factors that affect the growth and development of the cray sh ( Fig S2). Our ndings provided a basis for better understanding of the biology of this species at the molecular level and the intestinal micro ora level. In this study, RNA-Seq and bioinformatics analyses were used to study the dry pancreas transcriptome information of different individuals with the same genetic background and the same feeding conditions. There were signi cant differences in the expression of 497 genes between large and small groups (276 up-regulated and 221 down-regulated genes). The GO annotation and KEGG pathway analyses showed that the DEGs were enriched in catalytic activity and binding of metabolic processes in the biological processes category, and single-cell, cellular processes, and cellular components of the cellular and membrane in the molecular functions category. Analysis of the transcriptomics data also revealed many key genes and important pathways related to immune and metabolic processes (Figs. 1 and 2).
We mainly validated genes that were up-regulated with RT-qPCR (Fig 5), including FASN (fatty acid synthase, animal type), SOD1 (superoxide dismutase, cu-zn family), LSC2 (succinyl-coa synthetase beta subunit), TYR (tyrosinase), CALM (calmodulin), ACC1 (acetyl-coa carboxylase/biotin carboxylase 1), ADIPOR (adiponectin receptor), CTSL (cathepsin L), ACLy (ATP citrate (pro-s)-lyase), and FATP1 (fatty acid transporter). These genes were primarily involved in biological processes associated with fatty acid biosynthesis, peroxisome, lysosome, ribo avin metabolism, estrogen signaling pathway, AMPK signaling pathway, antigen processing and presentation, fatty acid degradation, and the PPAR signaling pathway. FASN is an important enzyme in lipid metabolism and plays a role in a variety of biological processes. In mammals, FASN is associated with immunity and anti-bacterial activity [19]. Fatty acids can regulate phagocytic activity, the production of in ammatory cytokines and immunoglobulin, as well as proliferation, apoptosis, and migration of immune cells. It has been reported that palmitate, a product of FASN mediated biosynthesis, can promote monocyte adhesion by activating various immune factors in endothelial cells of human adipose tissue in the microvascular system [20,21]. Recent studies have found that when cray sh are infected with White Spot Syndrome Virus (WSSV), FASN expression can be induced through the PI3K/AKT pathway, improving lipid biosynthesis to support the morphogenesis of the virus, which suggests that FASN may play a positive role in WSSV infection (DOI:10.1016/j.dci.2015.06.0010). In order to gain insight into the key role of the FASN gene in the metabolic process, we investigated the speci c involvement of FASN in the Fatty acid biosynthesis signaling pathway, which shows that FASN is involved in almost all aspects of the process. Investigating the speci c regulatory mechanisms of the FASN gene in the growth process will be the focus of our next study. The FATP family is a group of proteins that is predicted to be part of a fatty acid transport pathway. Six distinct subtypes of FATP have been identi ed in mammalian systems, which take part in the importation of exogenous fatty acids or the activation of long-chain fatty acids. It remains controversial whether these proteins function as membrane-bound fatty acid transporters or as acyl-coa synthases, which activate long-chain fatty acids and accompany their transportation [22]. The damage of aging on the body is usually related to certain structural and functional defects in organs [23]. These defects include oxidative modi cation, protein aggregation, and differential expression of speci c genes. Studies in mammals have found that certain oxidative damage is related to aging (DOI: 10.1016/j.exger.2019.110795). SOD1 is an important member of the super-oxidase family, which has been demonstrated to have ROS scavenging ability [24]. During growth, the body produces a large amount of non-speci c oxidation molecules such as proteins, nucleic acids, and lipids, and SOD1 acts as an oxidative scavenger to effectively control the active oxygen content, thereby reducing the accumulation of oxidative damage. Recent studies have found that SOD1 can also take part in metabolism, transcription, and maintenance of redox balance [25]. Our differential expression analysis showed that SOD1 gene expression was signi cantly higher in the LI group than in the SI group. Tyrosinase can be found in animals, plants, and fungi, and has a broad spectrum of substrates for a variety of phenolic compounds (L-tyrosine, L-DOPA, catechol, caffeic acid, tyramine, phenol, p-aminophenol, cresol, p-cresol, dopamine, 4-hydroxyanisole, L-isoproterenol, 4-ethoxyphenpl, 4-butylcatechol, and pyrogallol) [26]. Overall, these DEGs in the hepatopancreas may explain differential development of different P. clarkii, individuals with the same genetic background.
A healthy microbiome is not a constellation of individual species, but rather includes a network of interdependent microbes [27]. The microbiota in the gastrointestinal tract is a symbiotic partner of the host and is crucial for intestinal microvasculature development, metabolism and the development of multiple components of the host immune system [28,29]. The vertebrate intestine hosts 10 14 bacterial cells, including perhaps tens of thousands of bacterial strains [27]. Previous work in Litopenaeus vannamei demonstrated the close relationship between intestinal microbes and intestinal development, with bene cial intestinal ora contributing to overall host health [30]. In addition, intestinal ora can also have a protective effect against pathogens and can stimulate the host immune system through cell signaling in addition to producing antibacterial peptides and competing for adhesion locations (DOI:10.1038/nature11552). Improving the absorption of nutrients during the farming of aquatic animals is one key to increasing yield and quality. Some species, such as Penaeus vannamei, have been reported to select microorganisms for speci c functions according to their needs under certain selective pressures [31].
We found that at the phylum level, Proteobacteria were the most abundant intestinal bacteria in both groups. Proteobacteria has the largest phylogenetic composition and is comprised of 116 validated families. On average, there are 10.1 families per validated phylum (median, 3.0), according to the List of Prokaryotic Names with Standing in Nomenclature (LPSN). Members of the phylum Proteobacteria have greatly variable morphology and diverse physiology, which makes them competitive in surviving in various ecological niches [32]. Proteobacteria have been observed to be ubiquitous in habitats such as soil [32], plants [33], and seawater [34]. Most microorganisms in the animal gut are obligate anaerobes, but Proteus is an interstitial anaerobic bacterium that we speculate may affect host oxygen homeostasis. Microbiomes that are rich in diversity can provide the host with complementary genetic resources, which have an effect on energy extraction, essential vitamin production, intestinal maturation, and immune system development. Furthermore, many pathogenic Vibrio species in this family are potentially harmful to craw sh [35]. In the larger individual group, the proportion of Actinobacteria was second only to Proteobacteria (Table S3). Studies have shown that the larva development of South American white shrimp can alter the structure and function of intestinal bacteria, with the proportion of Actinobacteria increasing at a certain stage of development [36].
In this study, the proportions of Firmicutes and Actinobacteria bacteria in the SI group were relatively high. This correlation suggests that gut bacteria may in uence the speed at which the body develops over the same time. The species richness of Vogesella, Bacillus, Aeromonas, Bacteroides, De uviicoccus, Acinetobacter, Gemmobacter, and Rhodobacter was further compared at the genus level (Fig. 7). Gemmobacter, Acinetobacter, Candidatus-Bacilloplasma, Vogesella, Rhodobacter, Paracoccus, and Bosea had a signi cantly different richness in the LI group compared to the SI group (p<0.05).
The genus Vogesella belongs to the class Betaproteo and includes three recognized heterotrophic species (Fig. 7), including V. indigofera ATCC 19706T, Pergesida perlucida LMG 24214T and V. lacus LMG 24504T [37], which were all isolated from freshwater sites (in sediment or water). Some studies have con rmed the ability of these genera to degrade peptides and several other polymers. This is contrary to the current understanding that members of the genus Vogesella mainly degrade low molecular weight organic compounds, although this genus has not been studied extensively [38]. The different proportion of the taxonomic abundance of Rhodobacter represents an interesting nding since members of this genus can inhibit the growth of Vibrio [39]. Therefore, bacteria of this family are widely used as probiotics and be very e cacious [39]. According to some reports, although there are some Acinetobacter species in the microbiome of healthy human skin and animal intestines, many isolates recovered from clinical specimens show substantial pathogenicity [40]. Studies have shown that the secondary metabolites provided by the gut microbiota take part in the biologically active functions required by the host [41]. Thus, regulation of the host micro ora may represent a way to manage the growth of healthy host species.
By comparing the intestinal microbiota of the LI and SI groups, we found differences in many species associated with responding to environmental stress or preventing infection (Fig. 8). Carbohydrate metabolism function prediction at the second KEGG level revealed enrichment in pathways associated with transport, information-based processing, cellular processes, and signaling. Xenobiotics biodegradation and metabolism, carbohydrate metabolism, transport, genetic information processing, cellular-processes-and-signaling, and transcription processes are critical pathways for cray sh growth. Studies have shown that secondary metabolites provided by the gut ora are involved in essential biological functions of the host and analysis of transcriptomic KEGG data revealed that many genes are involved in the metabolic and immune processes of the organism. Exploring the intrinsic connection between the two sequencing approaches that we employed in the future will undoubtedly yield new ndings related to the relationship between P. clarkii microbiome and growth rate.
The intestinal tract is often subject to microbial parasitism, which requires pathogenic bacteria to adhere to the intestinal mucosa. The genetic diversity of different species of cray sh, despite their similar kinship, is still present, and the differences in transcription and translation of structural and functional genes have implications for growth and development. The physicochemical properties of the intestinal tract and the type of nutrients also determine the structure and function of the adherent symbiotic intestinal microorganisms.
The rate of development of cray sh is in uenced by many factors, including host and non-host factors. The host factors are genetic, developmental, and physiological state speci c, which in turn affect the speci city of the gut ora. The non-host factors (water environment, feed, drugs, and probiotics) may also affect gut microorganisms in the acquired environment. Overall, the complex and variable gut ora is a balanced system of interactions between host, environment, and microorganisms.

Conclusion
In this study, we analyzed the differential expression of the hepatopancreatic transcriptome and compared differences in intestinal micro ora structure between small and large cray sh siblings reared under the same conditions. We found differences in gene expression related to growth and development among different offspring, despite similar genetics. The rapidly growing individuals showed signi cant differences in the expression of genes that are involved in coping with environmental changes, dealing with disease infection, and the processing of metabolism related to growth. The structure of the microbiome may also play a role in these processes since the structure of intestinal micro ora adapts to the environment to promote the development of an organism. In addition, transcriptomic and 16S ribosomal RNA results, including signi cantly enriched GO and KEGG pathways and DEGs involved in organism development, provide further evidence that differences between individuals with different growth rates are accompanied by changes in intestinal ora structure. Our ndings have implications for understanding the roles that the host and non-host environments play in in uencing growth differences among cray sh individuals.

Materials And Methods
Selection and feeding of experimental animals P. clarkii were selected from sexually mature individuals from the Dongping Lake area of Tai'an, Shandong Province, China. Specimens (35-40 g) of male and female cray sh were cultured in circulating water at 24-27°C. Feed was purchased from New Hope Liuhe Group (Linyi, China). After stabilizing the samples for one week, one male and four females were selected from healthy individuals and reared in the same tank (2 m diameter and 30-40 cm depth) for mating. After mating, the fertilized female cray sh were reared individually until the juvenile cray sh hatched successfully. All juveniles 280 in total, were fed under the same conditions until 45 days of age, during which feed, water temperature, density, dissolved oxygen, nitrite, and ammonia levels were kept within the appropriate ranges.

Collection of intestinal contents and hepatopancreas
Animals were selected and collected after 45 days of feeding. Six of the largest and smallest healthy individuals were selected to represent the large and small study groups. The cray sh were sterilized with 75% alcohol, followed by extraction of as much blood as possible from the pericardial cavity with a 1 ml syringe with anticoagulant added inside the syringe in advance until after the experimental animals were incapacitated, followed by dissection of the hepatopancreas and intestinal contents, which were then aseptically placed in a 1.5 ml centrifuge tube and frozen in liquid nitrogen before storage at -80°C, and hemolymph was stored in 1.5ml Eppendorf tubes, subjected to 10 min of centrifugation at 4000 rpm at 4°C. The supernatant (serum) was collected, packaged, and stored at -80 °C. The intestines of experimental animals from the LI and SI groups were collected and analyzed for tissue sections.
RNA extractions, library construction and RNA sequencing RNA was isolated using Trizol Reagent, according to the manufacturer's instructions (Accurate Biotechnology Co. Ltd, Hunan, China). The quality and integrity of these total RNA samples were determined by a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). 1.5% agarose gels were used to determine RNA integrity. We constructed sequencing libraries from gills using a TruseqTM RNA sample prep Kit for the Illumina HiSeq 2000 platform, following the manufacturer's instructions.
Transcriptome assembly and functional annotation of unigenes A total of 1.5 µg RNA per sample was used as input material for RNA sample preparation. Sequencing libraries were generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA), following the manufacturer's recommendations. Samples were then barcoded, followed by clustering and sequencing using Illumina sequencing technology (Novozymes Biotechnology, Beijing, China). In this step, reads containing adapters, reads containing ploy-N, and low-quality reads were removed from the raw data to obtain clean reads. The Q20, Q30, GC-content, and sequence repeatability of the clean data were also calculated.

Gene functional annotation
All downstream analyses were based on high-quality clean reads. BLAST was used to compare unigene sequences against the NCBI non-redundant protein

Analysis of differentially expressed genes (DEGs) and functional enrichment
Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value <0.05 found by DESeq were assigned as differentially expressed DESeq was utilized to identify DEGs across samples which contained biological replicates. The resulting DEGs were then subjected to enrichment analysis of the three major categories of GO (biological process [BP], cellular component [CC], and molecular function [MF]). Next, KEGG analysis was carried out to identify signi cantly enriched metabolic pathways using KOBAS (2.0).

Validation by quantitative Real-Time PCR (RT-qPCR).
RT-qPCR was used to validate selected candidates from the DEGs which were previously identi ed using RNA-Seq. RT-qPCR reactions were carried out using a LightCycler® 2.0 (Roche). A partial sequence of the P. clarkii 18S ribosomal RNA gene (NCBI:AF436001) was used as the reference gene. The PCR was performed using a SYBR Premix Ex Taq kit (Accurate Biotechnology Co., Ltd), and the 20 μl PCR reaction systems consisted of 10 μl SYBR Premix Ex Taq (2×), 0.4 μl of each gene-speci c primer (10 nmol), 2 μl cDNA, and RNase free water up to 20 μl. Cycling conditions were 95°C for 30 s, followed by 40 cycles of 95℃ for 5 s and 60℃ for 30 s. Each sample was run in triplicate.
Sequencing and data processing of intestinal microbial genomic DNA Total genomic DNA from samples was extracted using the CTAB/SDS method. DNA concentration and purity were monitored on 1% agarose gels. According to the concentration, DNA was diluted to 1 ng/µl using sterile water. All PCR reactions were carried out in 30 µl reactions with 15 µl of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 µM of forward and reverse primers and approximately 10 ng template DNA. Thermal cycling consisted of initial denaturation at 98℃ for 1 min, followed by 30 cycles of denaturation at 98℃ for 10 s, annealing at 50℃ for 30 s, and elongation at 72℃ for 30 s. Reactions were nished with a nal 5 min at 72℃. Sequencing libraries were generated using an Ion Plus Fragment Library Kit (Thermo Scienti c), following manufacturer's recommendations. The library quality was assessed using a Qubit 2.0 Fluorometer (Thermo Scienti c). Finally, the library was sequenced on an Ion S5TM XL platform to generate 400/600 bp single-end reads.
To study the species composition of each sample, operational taxonomic unit (OTU) clustering was performed on valid tags with 97% identity for all samples, and then annotated with representative sequences of OTUs. Cutadapt (V1.9.1, http://cutadapt.readthedocs.io/en/stable/) was used to prune low-quality reads.
After initial quality control, chimeric sequences were removed by comparison with the species annotation database. The Read Sequence (https://github.com/torognes/vsearch/) Upanish v7.0.1001 software (http://www.drive5.com/uparse/) was employed for clustering clean reads from all samples. Sequences with 97% identity were clustered, and the sequence with the highest OTU frequency was selected as the representative sequence. OTU sequences were annotated using the Mothur method, and species annotation analysis was performed using the SSUrRNA database of SILVA 132 (http://www.arb-silva.de/), with thresholds set at 0.8 to 1.0. Classi cation information and statistics were generated for each taxonomic level, including kingdom, boundary, order, family, genus, and species.

Intestinal microbial function prediction methods
Functional prediction was carried out using the nearest neighbor method based on the similarity of 16S ribosomal RNA sequences in the R package Tax4Fun [17]. Speci cally, 16S ribosomal RNA sequences from all prokaryote genomes were extracted from the KEGG database and then were compared to the SILVA, SSU Ref, and Nr databases (BLAST bit score >1500) using the BLASTN algorithm to establish a correlation matrix. The KEGG database results, annotated via UProC and PAUDA, were mapped to the SILVA database for functional annotation.

Histological analysis
The intestines of experimental animals in the LI and SI groups were collected and analyzed in tissue sections. After paraformaldehyde xation for 24h, the intestines were embedded in para n blocks and dehydrated with different concentrations of ethanol. The blocks were sectioned to 5μm thickness with a slicer, stained with hematoxylin and eosin, and resin-xed for examination by an optical microscope.

Humoral factor testing
To detect the viability of humoral factors in animals at different group, hemolymph supernatants from each group as described in 2.1 were randomly selected and mixed according to group distribution, with three biological replicates for each group. The activities of plasma phenol oxidase (Ppo), lysozyme (Lyz), glutaminase (Gls), lipase, amylase (AMS), and superoxide dismutase (SOD) of the aforementioned humoral factors were determined using a commercial assay kit (Nanjing Jiancheng Institute, Nanjing, China) according to the manufacturer's instructions.

Statistical analysis
Statistical analysis was performed using SPSS Version 21.0 (IBM Corporation, USA) and GraphPad Prism version 8.0 (GraphPad Software Inc., San Diego, CA, USA) software. Statistical differences were determined via one-way analysis of variance (ANOVA) and Tukey's multiple range tests. Signi cance and high signi cance were set at p-value < 0.05 and p-value <0.01, respectively.

Abbreviations
Not applicable Declarations Ethics approval and consent to participate: This study was conducted following the 2012 International Guidelines for Biomedical Research Involving Animals (Council for International Organizations of Medical Sciences, http://www.cioms.ch), the care and use of laboratory animals were in full compliance with these guidelines. Furthermore, the study species is not listed in the List of Protected Animals in China and experimental research on this species is legal in China no institutional permission was required for the collection of animals in this study.

Consent for publication
Not applicable.
Availability of data and materials section: The datasets generated and/or analyzed during this study are not publicly available due to copyright and the progress of current research, and we will also delve into speci c mechanisms and applications, but are available from the corresponding authors upon reasonable request.

Competing interests
The authors declare that they have no competing interests. LG: Visualization, Investigation.
All authors have read and approved the manuscript.  Total Unigenes 44372 100 Table 3 Enrichment and annotation of the top ten species of intestinal ora at the phylum level.    KEGG metabolic pathway classi cation statistics. The ordinate is the name of the KEGG metabolic pathway, while the abscissa is the number of genes annotated for the pathway and their proportion to the total number of annotated genes. Based on the KEGG metabolic pathways involved, the genes are classi ed into ve branches: a) cellular processes, b) environmental information processing, c) genetic information processing, d) metabolism and e) organic systems.

Figure 5
Expression levels of candidate unigenes were assessed by RT-qPCR, and the expression trend was the same as that of RNA-Seq.

Figure 6
Species abundance clustering map Genus-level species phylogeny Figure 8