Method for analyzing intestinal health of crisp tilapia mossambica

Abstract

The invention relates to the technical field of genomics and metabonomics, in particular to a method for analyzing intestinal health of crisp tilapia mossambica. The tilapia mossambica intestinal health judgment model is established by dynamically monitoring tilapia mossambica muscle texture and intestinal morphology through high time resolution, analyzing intestinal health three-dimensional indexes (tissue morphology-flora structure-host gene expression) and serum metabolism response after embrittlement succeeds, and combining cross-omics data association analysis. The change condition of the muscle texture can be known more timely, a more accurate basis is provided for adjustment of a feeding strategy, and the stability of the fish meat quality can be improved. The intestinal health is evaluated from three dimensions of histomorphology, flora structure and host gene expression, so that the health condition of the intestinal tract can be known more comprehensively, and more comprehensive reference is provided for optimization of a feeding method. The correlation between different omics data can be deeply analyzed, the physiological and metabolic mechanisms of the fish body are disclosed, and deeper theoretical support is provided for optimization of a feeding method.

Description

Technical Field

[0001] This invention patent relates to the technical fields of genomics and metabolomics, and specifically relates to a method for analyzing the intestinal health of crispy tilapia. Background Art

[0002] Aquaculture, as an important source of global animal protein supply, faces the dual challenges of optimizing the farming cycle and managing animal health. During the intensive farming of fish, intestinal health is a key indicator determining nutrient absorption efficiency, immune function, and stress resistance. In recent years, the integrated application of multi-omics technologies (metabolomics, transcriptomics, microbiomics) has provided a new paradigm for comprehensively evaluating the physiological state of farmed organisms. Research has shown that the intestinal microbiota structure of fish is significantly correlated with the host individual (Meng Xiaolin, 2019, Acta Hydrobiologica Sinica), and serum metabolome can reflect the overall physiological state of the host.

[0003] Tilapia is the second most farmed fish globally. Its crispy variety forms a unique muscle texture through specific feed regulation, but the relationship between the farming cycle and fish health has not been clarified. Existing research mainly focuses on the impact of the crisping process on muscle quality (such as: Qin Zhiqing, 2023, Feed Research; Xie Xi, 2022, Acta Hydrobiologica Sinica; Liang Ping, 2022, Fujian Agricultural Science and Technology; Qingqing Li et al., 2023, Frontiers in Nutrition), and the research on crispy feed (such as: Peng Kai, 2023, Feed Industry). There are obvious gaps in the research on the intestinal health system. The determination of the traditional farming cycle mainly relies on growth texture or muscle section, etc., and the research on intestinal health only focuses on intestinal morphology and microbiota, etc. (such as: Fu Bing, 2023, Chinese Journal of Animal Nutrition; Xiaogang He, 2023, Animals), lacking a multi-dimensional health evaluation system.

[0004] Existing patents mainly focus on the feeding methods and feed formulations of crispy tilapia, and pay more attention to the effect of improving the meat quality of crispy tilapia. Regarding the intestine, it mainly involves reducing the occurrence of intestinal inflammation and improving digestion and absorption, lacking patents on in-depth research of the intestinal health system of crispy tilapia.

[0005] For example, the patent application with the publication number CN113974025A discloses a method for raising crispy tilapia. By adding black soldier fly larvae powder to the feed for crispy tilapia and feeding feeds with different compositions to crispy tilapia at different stages, the results show that this method can reduce the breeding cost of crispy tilapia, improve the body immunity of crispy tilapia, and reduce the occurrence probability of intestinal inflammation in crispy tilapia. Another example is that the patent application with the publication number CN113974025A discloses a feed for crispy tilapia and its preparation method. The tilapia feed uses mixed bacteria fermented soybean meal to replace the protein source, and the fermentation effectively removes various anti-nutritional factors such as trypsin inhibitors and goitrogens in soybean meal. Using broad bean, wheat, soybean meal, wheat bran, peanut bran, and secondary powder as the main raw materials, it improves the crude fat content and collagen content in the muscle of tilapia. Microwave ripening is used to increase the fragrance of the feed, and low-temperature grinding and ultrafine grinding are carried out, which is convenient for puffing processing, preserves the nutrients in the feed, and is convenient for young fish to digest. The above patents all relate to the healthy breeding of tilapia, but do not deeply study the specific mechanism and influencing factors of intestinal inflammation.

[0006] The prior art often has the following limitations:

[0007] (1) Single index: Most existing patents and literatures use growth performance (such as specific growth rate) or single tissue indexes (such as muscle texture, intestinal enzyme activity, liver enzyme activity, and liver metabolism, etc.) as the judgment basis for the breeding cycle, lacking systematic health assessment.

[0008] (2) Weak mechanism interpretation: The existing solutions fail to establish a linkage analysis framework of "microbiota-host gene expression-metabolic response", and cannot reveal the molecular mechanism of the impact of crispy feed on the intestinal health of tilapia.

[0009] (3) Lack of breed specificity: The existing research on crispy fish focuses on grass carp, lacking research on the specific physiological characteristics of the intestinal health of tilapia after crispy treatment. Summary of the Invention

[0010] In order to solve the above technical problems, the purpose of the present invention is to provide an analysis method for the intestinal health of crispy tilapia. By dynamically monitoring the muscle texture and intestinal morphology of tilapia with high time resolution, after successful crisping, analyze the three-dimensional indexes of intestinal health (tissue morphology-microbiota structure-host gene expression) and serum metabolic response, and combine cross-omics (serum metabolome-intestinal transcriptome, serum metabolome-intestinal content microbiome) data correlation analysis to establish a method for judging the intestinal health of tilapia.

[0011] The purpose of the present invention can be achieved by the following technical solutions:

[0012] The present invention provides an analysis method for the intestinal health of crispy tilapia, including the following steps:

[0013] (1) After domesticating tilapia for half a month, it was divided into a control group and an experimental group. The control group was fed with commercial tilapia feed, and the experimental group was fed with brittle tilapia feed.

[0014] (2) On the 30th day, 60th day, 90th day, and 120th day respectively, samples were selected from the brittle group and the control group for muscle tissue sectioning, texture analysis, and intestinal tissue sectioning.

[0015] (3) On the 120th day, samples were selected from the brittle group and the control group respectively for detecting the activity of intestinal antioxidant enzymes, Q-PCR detection of intestinal anti-inflammatory factor genes, ELISA protein detection of the intestine, 16S rRNA sequencing of intestinal contents, serum metabolome detection, and intestinal epithelial transcriptome sequencing.

[0016] (4) Dynamic brittle degree analysis: Understand the muscle quality through muscle tissue sectioning and texture analysis, and preliminarily understand the degree of intestinal damage through intestinal tissue sectioning.

[0017] (5) Three-dimensional evaluation of intestinal health: Conduct intestinal morphological evaluation through the villus height and crypt depth of intestinal sections, analyze the intestinal flora structure through α-diversity, β-diversity, and functional prediction in 16S rRNA sequencing of intestinal contents, and conduct host gene analysis through screening pathway enrichment analysis of differentially expressed genes in intestinal epithelial transcriptome sequencing.

[0018] (6) Cross-omics integration analysis: Conduct joint analysis through serum metabolome detection - host gene analysis of intestinal epithelium and serum metabolome detection - intestinal flora structure analysis, integrate the relationship between intestinal genes and microorganisms and the overall metabolic level of fish, and conduct cross-omics analysis of the intestinal health of crispy tilapia.

[0019] Further, in step (1), the domestication refers to feeding twice a day (8:00 and 17:00) with apparent satiation using commercial tilapia feed until the fish body is stable and adapts to the experimental environmental conditions.

[0020] Further, in step (1), the weight of the tilapia is 350 - 450 g; three parallels are set for each of the control group and the experimental group, with 30 tilapia in each parallel; the number of samples is 9 tilapia for each of the control group and the experimental group.

[0021] Further, in step (2), the operation and analysis method of the muscle tissue sectioning include the following steps: material sampling, dehydration, clearing, wax infiltration, embedding, sectioning, baking, and HE staining; the texture analysis is measured using a Universal TA texture analyzer; the operation and analysis method of the intestinal tissue sectioning include the following steps: material sampling, dehydration, clearing, wax infiltration, embedding, sectioning, baking, and HE staining.

[0022] Further, in step (3), the detection of intestinal antioxidant enzyme activity is determined using a kit from Nanjing Jiancheng Bioengineering Institute.

[0023] Further, in step (3), the method for Q-PCR detection of intestinal anti-inflammatory factor genes includes the following steps: extracting total RNA from tilapia intestines, quantifying the RNA concentration using a spectrophotometer, removing genomic DNA, reverse transcribing it into cDNA using a reverse transcription kit, and finally performing real-time fluorescence quantitative PCR detection.

[0024] Further, in step (3), the intestinal ELISA protein detection is determined using a kit from Jiangsu Enzyme Immunoassay Industry Co., Ltd.

[0025] Further, in step (3), the method for 16S rRNA sequencing of intestinal contents includes the following steps: extracting total microbial group DNA, PCR amplification of target fragments, magnetic bead purification and recovery of amplification products, fluorescence quantification of amplification products, preparation of sequencing libraries, and high-throughput sequencing on a machine.

[0026] Further, in step (3), the serum metabolome detection is performed by Suzhou Panomic Biopharmaceutical Technology Co., Ltd. for metabolome determination.

[0027] Further, in step (3), the method for intestinal epithelial transcriptome sequencing includes the following steps: extraction of total RNA, quality detection of total RNA, purification of mRNA, fragmentation of mRNA, cDNA synthesis, PCR enrichment of library fragments, library quality inspection, and sequencing on an Illumina platform.

[0028] The beneficial effects that this application can produce are as follows:

[0029] (1) A real-time monitoring system for dynamic muscle texture changes with high temporal resolution (≤30-day interval) can understand the changes in muscle texture more timely. Compared with the traditional sampling interval of ≥60 days, it can capture the nodes of muscle texture changes and intestinal morphological changes more accurately, providing a more accurate basis for adjusting feeding strategies and helping to improve the stability of fish meat quality; through high-temporal-resolution monitoring, the key nodes of muscle texture changes can be grasped more precisely, providing more detailed data support for scientific research and production practice.

[0030] (2) A three-dimensional evaluation system for intestinal health (histomorphology - microbial community structure - host gene expression) comprehensively evaluates intestinal health from three dimensions of histomorphology, microbial community structure, and host gene expression, can understand the health status of the intestine more comprehensively, and provides a more comprehensive reference for optimizing feeding methods; based on scientific research methods and indicators, it can evaluate intestinal health more scientifically and provide a more reliable basis for feeding management.

[0031] (3) The ability of cross-omics data correlation analysis (metabolomics-transcriptomics-microbiome) can deeply analyze the correlations between different omics data, reveal the physiological and metabolic mechanisms of fish bodies, and provide more in-depth theoretical support for the optimization of feeding methods; by integrating multiple omics data, it is possible to understand the health status and growth performance of fish bodies from a systematic level and provide more systematic guidance for feeding management. Compared with only using single omics indicators (such as only gut microbiota or serum metabolomics), it can better explain the "microbiota-host" interaction mechanism and avoid the problem that single indicators are easily interfered by accidental factors. Description of the Drawings

[0032] Figure 1 It is a histological section diagram of the intestinal tract of the control group and the brittle group on the 30th day, where the left side is the control group and the right side is the experimental group. From top to bottom, they are the anterior intestine, middle intestine, and posterior intestine.

[0033] Figure 2 It is a histological section diagram of the muscle tissue of the control group and the experimental group on the 30th day, where the left side is the control group and the right side is the brittle group.

[0034] Figure 3 It is a histological section diagram of the intestinal tract of the control group and the brittle group on the 60th day, where the left side is the control group and the right side is the experimental group. From top to bottom, they are the anterior intestine, middle intestine, and posterior intestine.

[0035] Figure 4 It is a histological section diagram of the muscle tissue of the control group and the experimental group on the 60th day, where the left side is the control group and the right side is the brittle group.

[0036] Figure 5 It is a histological section diagram of the intestinal tract of the control group and the brittle group on the 90th day, where the left side is the control group and the right side is the experimental group. From top to bottom, they are the anterior intestine, middle intestine, and posterior intestine.

[0037] Figure 6 It is a histological section diagram of the muscle tissue of the control group and the experimental group on the 90th day, where the left side is the control group and the right side is the brittle group.

[0038] Figure 7 It is a histological section diagram of the intestinal tract of the control group and the brittle group on the 120th day, where the left side is the control group and the right side is the experimental group. From top to bottom, they are the anterior intestine, middle intestine, and posterior intestine.

[0039] Figure 8 It is a histological section diagram of the muscle tissue of the control group and the experimental group on the 120th day, where the left side is the control group and the right side is the brittle group.

[0040] Figure 9 It is a diagram of the change in muscle texture of the control group and the experimental group, where b represents the experimental group and c represents the control group.

[0041] Figure 10 It is a comparison chart of growth parameters for the control group and the experimental group, where b represents the experimental group and c represents the control group.

[0042] Figure 11 It is a comparison chart of intestinal antioxidant enzyme activities for the control group and the experimental group, where b represents the experimental group and c represents the control group.

[0043] Figure 12 It is a flow chart for serum metabolome analysis.

[0044] Figure 13 It is a PLS-DA score chart obtained from serum metabolome analysis, where b represents the experimental group and c represents the control group.

[0045] Figure 14 It is a statistical bar chart of primary differential metabolites obtained from serum metabolome analysis.

[0046] Figure 15 It is a volcano plot of primary differential metabolites obtained from serum metabolite analysis.

[0047] Figure 16 It is a statistical bar chart of secondary differential metabolites obtained from serum metabolome analysis.

[0048] Figure 17 It is a heat map of secondary differential metabolites obtained from serum metabolome analysis.

[0049] Figure 18 It is a volcano plot of secondary differential metabolites obtained from serum metabolome analysis.

[0050] Figure 19 It is a scatter plot of m / z ratio and P value of secondary differential metabolites obtained from serum metabolome analysis.

[0051] Figure 20 It is a bubble plot of metabolic pathway impact factors obtained from serum metabolome analysis.

[0052] Figure 21 It is a flow chart of the experiment for intestinal epithelial transcriptome sequencing detection.

[0053] Figure 22 It is a flow chart of the analysis for intestinal epithelial transcriptome sequencing detection.

[0054] Figure 23 It is a statistical bar chart of differential expression detection between two samples obtained from the analysis of intestinal epithelial transcriptome sequencing detection.

[0055] Figure 24 It is a volcano plot of differential expression detection between two samples obtained from the analysis of intestinal epithelial transcriptome sequencing detection.

[0056] Figure 25 It is a clustering analysis chart obtained from the analysis of intestinal epithelial transcriptome sequencing detection.

[0057] Figure 26 It is a bar chart of GO enrichment analysis obtained from the detection and analysis of intestinal epithelial transcriptome sequencing.

[0058] Figure 27 It is a bubble chart of GO enrichment analysis obtained from the detection and analysis of intestinal epithelial transcriptome sequencing.

[0059] Figure 28 It is a bar chart of KEGG enrichment analysis obtained from the detection and analysis of intestinal epithelial transcriptome sequencing.

[0060] Figure 29 It is a bubble chart of enrichment analysis obtained from the detection and analysis of intestinal epithelial transcriptome sequencing.

[0061] Figure 30 It is a flow chart of the analysis of 16S rRNA sequencing detection of intestinal contents.

[0062] Figure 31 It is a taxonomic composition analysis chart obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents.

[0063] Figure 32 It is an Alpha diversity index chart obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents.

[0064] Figure 33 It is a Beta diversity analysis (PCoA analysis, weighted UniFrac distance) chart obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents.

[0065] Figure 34 It is a Venn diagram of species difference analysis and marker species ASV / OTU obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents.

[0066] Figure 35 It is a heat map of species composition obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents.

[0067] Figure 36 It is a heat map of species composition (between groups) obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents.

[0068] Figure 37 It is a LEfSe analysis (bar chart of the LDA effect value of marker species) obtained from the detection and analysis of 16S rRNA sequencing of intestinal contents. Detailed implementation methods

[0069] Next, the technical solutions in the embodiments of the present application will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present application without creative efforts shall fall within the protection scope of the present application.

[0070] Example

[0071] In this example, the ingredients of tilapia commercial feed and tilapia crisping feed are shown in Table 1 below.

[0072] Table 1

[0073]

[0074] 180 tilapia weighing about 400 g were fed twice a day (8:00 and 17:00) with tilapia commercial feed until the fish bodies were stable and adapted to the experimental environmental conditions. The domestication time was about 2 weeks. Then, they were divided into a control group of 90 tails and an experimental group of 90 tails. The control group was fed with tilapia commercial feed, and the experimental group was fed with tilapia crisping feed. Both the control group and the experimental group had 3 replicates, with 30 tails in each replicate.

[0075] On the 30th day, 60th day, 90th day, and 120th day respectively, 9 tails were selected as samples from the crisping group and the control group. After being anesthetized with anesthetic MS-222 (60 mg, Sigma, GYT0202813), the samples on the 120th day were bled from the caudal vein with a sterile blood collection needle into a sterile blood collection negative pressure tube and centrifuged at 3000 r / min for 10 min to obtain the supernatant, which was stored in an ultra-low temperature refrigerator at -80 °C. On the 30th day, 60th day, 90th day, and 120th day, muscle and intestinal tissues were collected under ice bath and sterile conditions, and intestinal contents were collected from the samples on the 120th day. Among them, the muscle was divided into three parts: one part was put into a general tissue fixing solution (biosharp, White Shark Biotechnology Co., Ltd.) for section analysis, one part was cut into a specified size (three pieces were taken from each fish) and immediately subjected to texture analysis, and the last part was quickly frozen in liquid nitrogen and stored in a refrigerator at -80 °C for subsequent experimental analysis; the intestine was divided into two parts: after being rinsed clean with sterile normal saline and slightly drained, one part was put into a general tissue fixing solution for section analysis, and the other part was quickly frozen in liquid nitrogen and stored in a refrigerator at -80 °C for subsequent experimental analysis; the intestinal contents were quickly frozen in liquid nitrogen and stored in a refrigerator at -80 °C for subsequent experimental analysis.

[0076] The methods for analyzing muscle tissue sections and intestinal tissue sections are as follows:

[0077] Place the tissue in a labeled embedding cassette, and then dehydrate, clear, and infiltrate the tissue with 70% alcohol for 1.5 h, 80% alcohol for 30 min, 90% alcohol for 30 min, 95% alcohol I for 30 min, 95% alcohol II for 30 min, 100% alcohol I for 30 min, 100% alcohol II for 30 min, xylene-absolute ethanol mixture (1:1) for 10 min, xylene I for 5 min, xylene II for 5 min, paraffin I at 58 - 60 °C for 30 min, and paraffin II at 58 - 60 °C for 60 min. Then open the lid of the embedding cassette, break it off and discard it. Use forceps to pick out the tissue and place it on the wax table. Place the bottom of the embedding cassette on one side of the wax table. Subsequently, use forceps to take out the metal mold from the wax bath and place it on the wax table. Add a layer of paraffin liquid to the metal mold. When the surface of the liquid turns white, use forceps to pick up the tissue and put it into the paraffin liquid. After waiting for 3 s, place the white embedding cassette on the metal mold. Then add fresh paraffin to the mold to completely cover the tissue. Stop adding paraffin when the surface of the liquid is parallel to the edge of the embedding cassette. Translate the metal mold with the dropped paraffin to the cooling table in the right hand. The wax block on the cooling table needs to be placed for 30 min. Place the wax block with the metal mold in a -20 °C refrigerator for 10 min. Separate the metal mold from the wax block. Put the embedded wax block into a self-sealing bag and place it at 4 °C overnight. Take the wax block out of the 4 °C refrigerator and pre-cool it in a -20 °C refrigerator for 30 min. During the pre-cooling process, prepare the glass slides and make corresponding marks on the bottom area of the glass slides with a pencil. Fix the blade on the microtome and confirm that the microtome is in the "locked" state. At the same time, turn on the oven and add 3 / 4 of ultrapure water to the basin. When the water temperature reaches 41 °C, fix the pre-cooled wax block on the microtome. First, adjust the trimming thickness to 10 μm. When it is found that the entire morphology of the tissue can be cut on the same slide, adjust the trimming thickness to 4 μm. Use forceps to pick up a corner of the cut slide and move it to the water basin of the oven. After the slide is flattened, use the marked glass slide to fish out the slide. Observe under the microscope whether the slide is fully flattened. After confirmation, place the slide on the oven. Set the temperature of the oven to 60 °C. Bake the section in a 60 °C oven for 4 h. Then dewax with xylene I for 15 min, xylene II for 15 min, 100% ethanol for 5 min, 95% ethanol for 5 min, 80% ethanol for 5 min, 70% ethanol for 5 min, ultrapure water I for 5 min, and ultrapure water II for 5 min. Stain with hematoxylin stain for 1 min, rinse with ultrapure water for 3 times, blue in tap water for 10 min, stain with eosin stain for 2 min, rinse with tap water for 5 times. Place the section in a 60 °C oven for drying and seal the section with neutral balsam. The intestinal tissue section diagrams of the control group and the embrittlement group on the 30th day are as follows Figure 1 shown, where the left side is the control group and the right side is the experimental group. From top to bottom are the anterior intestine, middle intestine, and posterior intestine. The muscle tissue section diagrams of the control group and the experimental group on the 30th day are as follows Figure 2As shown, the intestinal tissue section diagrams of the control group and the embrittlement group on the 60th day are as follows Figure 3 As shown, with the control group on the left and the experimental group on the right, from top to bottom are the anterior intestine, middle intestine, and posterior intestine. The muscle tissue section diagrams of the control group and the experimental group on the 60th day are as follows Figure 4 As shown, the intestinal tissue section diagrams of the control group and the embrittlement group on the 90th day are as follows Figure 5 As shown, with the control group on the left and the experimental group on the right, from top to bottom are the anterior intestine, middle intestine, and posterior intestine. The muscle tissue section diagrams of the control group and the experimental group on the 90th day are as follows Figure 6 As shown, the intestinal tissue section diagrams of the control group and the embrittlement group on the 120th day are as follows Figure 7 As shown, with the control group on the left and the experimental group on the right, from top to bottom are the anterior intestine, middle intestine, and posterior intestine. The muscle tissue section diagrams of the control group and the experimental group on the 120th day are as follows Figure 8 As shown.

[0078] The texture analysis method is as follows:

[0079] Determination was carried out using a Universal TA texture analyzer (Shanghai Tengba Instrument Technology Co., Ltd.). A p / 36R cylindrical probe was used on the analyzer. The pre-test speed was 2 mm / s, the post-test speed was 5 mm / s, and the compression speed test was carried out at a test speed of 1 mm / s. The compression time interval was 2 s, and the compression ratio was 25%. The muscle texture change diagrams of the control group and the experimental group obtained are as follows Figure 9 As shown, it can be seen that there are significant differences in muscle hardness, elasticity, adhesiveness, chewiness, cohesiveness, and resilience between the control group and the experimental group.

[0080] The growth parameter analysis method is as follows:

[0081] During sampling, body weight and total food intake were recorded, and the specific growth rate, feed coefficient, feed efficiency, and condition factor were calculated. The growth parameter comparison diagrams of the control group and the experimental group obtained are as follows Figure 10 As shown.

[0082] On the 120th day, samples were selected from the embrittlement group and the control group respectively for intestinal antioxidant enzyme activity detection, intestinal anti-inflammatory factor gene Q-PCR detection, intestinal ELISA protein detection, intestinal content 16S rRNA sequencing, serum metabolome detection, and intestinal epithelial transcriptome sequencing;

[0083] The intestinal antioxidant enzyme activity was measured using a kit from Nanjing Jiancheng Bioengineering Institute. The intestinal antioxidant enzyme activity comparison diagrams of the control group and the experimental group obtained are as follows Figure 11 As shown.

[0084] For the Q-PCR detection of intestinal anti-inflammatory factor genes, total RNA of tilapia intestine was extracted using the Bori MagaBio plus Total RNA Purification Kit (Hangzhou Bori Technology Co., Ltd.). After quantifying the RNA concentration using a spectrophotometer from TIANGEN (Tiangen Biochemical Technology Co., Ltd., Beijing), it was stored at -80 °C for later use. Subsequently, genomic DNA was removed, and cDNA was reverse transcribed using a reverse transcription kit (TaKaRa, Japan) in a MiniAmp TM Plus (USA) thermal cycler and stored at -80 °C for later use. Real-time fluorescence quantitative PCR (qRT-PCR) was performed using a kit (TaKaRa, Japan) on a Light 96-Time PCR Detection System (Roche, Switzerland).

[0085] For the ELISA protein detection of the intestine, a kit from Jiangsu Enzyme Immunoassay Industry Co., Ltd. was used for determination.

[0086] The experimental data of texture, growth parameters, intestinal antioxidant enzymes, intestinal anti-inflammatory factor genes, and intestinal ELISA proteins were statistically analyzed using SPSS 17.0 software. The experimental results were expressed as mean ± standard error (Mean ± SEM, n = 6). One-way analysis of variance (one-way ANOVA) was used for the analysis method. When there were significant differences between experimental groups, Duncan's multiple comparison was used, and statistical differences were indicated by * with P < 0.05. After the data analysis, graphs were plotted using Graphpad Prism 10.1.2 software.

[0087] The method for serum metabolome detection is as follows:

[0088] Serum samples stored at -80 °C during sampling were transported on dry ice to Suzhou Panomic Biomedical Technology Co., Ltd. for determination using a Thermo Orbitrap Exploris 120 with positive and negative ion mode switching. The parameters for the positive ion mode were: comparison: b vs c, pre: 3, R2X: 0.468 cum, R2Y: 0.998 cum, Q2: 0.891 cum; the parameters for the negative ion mode were: comparison: b vs c, pre: 2, R2X: 0.316 cum, R2Y: 0.99 cum, Q2: 0.771 cum; and then through such as Figure 12The analysis process shown uses software such as the R packages pheatmap, ropls, dendextend, cor, Wilcox.test, etc. for analysis. Through partial least squares discrimination analysis (PLS-DA), a regression model between metabolite expression levels and sample categories is established using partial least squares regression to achieve the prediction of sample categories. PLS-DA models for each comparison group are established, and the model evaluation parameters R2 (model interpretability) and Q2 (model predictability) are obtained through cross-validation. If R2 and Q2 are closer to 1, it indicates that the model is more stable and reliable. pre, number of principal components; R2X, model interpretability (for the X variable dataset); R2Y, model interpretability (for the Y variable dataset); Q2, model predictability. After analysis using the Ropls package in R language and plotting, the obtained PLS-DA score plot is as shown in Figure 13 shown, where the abscissa PC1 represents the first principal component score value, the ordinate PC2 represents the second principal component score value, the points represent samples, represents the 95% confidence interval, and the colors represent different groupings. Differentially expressed metabolites are searched from the list of primary metabolites, differential analysis is performed using the set statistical test method, and significant difference screening is carried out through P value and VIP. P value refers to the Student t-test, a statistically significant difference, and its screening threshold is P value less than 0.05. VIP refers to the variable importance in projection of the first principal component of OPLS-DA, and its screening threshold is VIP greater than 1.0. The statistical bar chart of primary differentially expressed metabolites plotted is as shown in Figure 14 shown, where the X-axis represents the number of differentially expressed metabolites, the Y-axis represents the comparison groups, red represents the up-regulated number, blue represents the down-regulated number, Statistic of Differently Expressed Metabolite: Statistics of differentially expressed metabolites, Metabolitecount: Metabolite count, Regulation: Up and down regulation relationship. The volcano plot of primary differentially expressed metabolites obtained is as shown in Figure 15As shown, the volcano plot can visually display the distribution and change trends of differential metabolites in two groups of samples. Usually, the abscissa is represented by log2(FC), and the ordinate is represented by -log10(P value), which is the negative logarithm of the statistical significance P value. The greater the change in the quantitative value, the more significant the difference, and the metabolites are distributed at both ends. Using the differential metabolite screening conditions such as the preset FC value, P value, and VIP, a metabolite volcano plot is drawn. The abscissa represents the logarithm of the Log2 of the difference multiple of the quantitative value of a metabolite in the two samples, and the ordinate represents the -log10 of the P value. Each point in the figure represents a metabolite. The larger the absolute value of the abscissa, the greater the difference in the expression multiple of a metabolite between the two samples. The larger the ordinate value, the more significant the differential expression, and the more reliable the screened differential expression metabolites. The size of the point represents the size of the VIP value. Red points represent upregulation differences, blue points represent downregulation differences, and gray points represent metabolites that do not meet the differential screening conditions. By default, the mz values of the top 5 metabolites with the smallest P value are displayed. P value: statistical p value, the smaller it is, the more significant the difference; -log10(P value): the -log10 value of the statistical P value; VIP: the variable importance in projection value of the first principal component of OPLS-DA; mz: mass-to-charge ratio; FC: the difference multiple value of a metabolite in different groups; log2(FC): the log2 value of the difference multiple of a metabolite in different groups. The difference multiple FC = average value (experimental group) / average value (control group). Substance identification is performed by retrieving and comparing (searching the library) using spectral databases such as HMDB, massbank, LipidMaps, mzcloud, KEGG, and the self-built metabolite standard database of Nomi Metabolomics. The metabolites containing secondary spectra in the quantitative list are compared and matched with the fragment ion and other information of each secondary spectrum in the database to achieve the secondary qualitative identification of metabolites. The secondary identified metabolites are selected from the identification results and screened using the preset P value and VIP thresholds in the statistical test of the primary differential analysis results to obtain secondary differential metabolites. The statistical bar chart of the obtained secondary differential metabolites is as shown in Figure 16 As shown, where the X-axis represents the number of differential metabolites, the Y-axis represents the comparison groups, red represents the number of upregulations, and blue represents the number of downregulations. The R Pheatmap package is used to scale the matrix data of secondary differential metabolites (Scale), and two-way clustering is performed on the samples and differential metabolites to draw a clustering heatmap. The heatmap of the obtained secondary differential metabolites is as shown in Figure 17 As shown, where the columns represent samples, the rows represent metabolites, the clustering tree on the left is the differential metabolite clustering tree, and the top is the sample clustering tree. The gradient color represents the size of the quantitative value. The redder the color, the higher the expression level, and the bluer the color, the lower the expression level. When the number of metabolites exceeds 150, the metabolite names are not displayed. The volcano plot of the obtained secondary differential metabolites is as shown inFigure 18 As shown, the abscissa represents the logarithm of the Log2 of the difference multiple of the quantitative values of a certain metabolite in two samples, and the ordinate represents the -log10 of the P value. Each point in the figure represents a metabolite. The larger the absolute value of the abscissa, the greater the difference in the expression level multiple of a certain metabolite between the two samples. The larger the ordinate value, the more significant the differential expression, and the more reliable the differentially expressed metabolites screened. The size of the point represents the VIP value. Red dots represent upregulated differences, blue dots represent downregulated differences, and gray dots represent metabolites that do not meet the differential screening conditions. By default, the names of the top 5 metabolites with the smallest P value are displayed. According to the mass-to-charge ratio of the metabolite and the P value, a differential scatter plot is drawn, which can display the differential characteristics under the mass-to-charge ratio size distribution of the metabolite. The scatter plot of the mass-to-charge ratio and P value of the secondary differential metabolites obtained is as Figure 19 shown. Among them, the abscissa represents the mass-to-charge ratio of a certain metabolite, and the ordinate represents the -log10 of the P value. Each point in the figure represents a metabolite. Red dots represent upregulated differences, blue dots represent downregulated differences, and gray dots represent metabolites that do not meet the differential screening conditions. The size of the point represents the VIP value. Red dots are used to represent differential metabolites. By default, the names of the top 5 metabolites with the smallest P value are displayed. Perform KEGG pathway enrichment analysis on the differential metabolite list. The enrichment method is based on the hypergeometric distribution test, and the topological analysis uses betweenness. The topological analysis aims to evaluate whether a given gene or metabolite plays an important role in biological reactions based on its position in the pathway. The bubble plot of the metabolic pathway impact factor obtained is as Figure 20 shown. Among them, the abscissa is the Impac value enriched in different metabolic pathways, the ordinate is the enriched pathway, the size of the point represents the number of corresponding metabolites on the pathway, and the color is related to the P value. The redder the color, the smaller the P value, and the bluer the color, the larger the P value. Hits: The total number of differential metabolites in the target metabolic pathway, Pvalue: The P value of the hypergeometric distribution test. The smaller the P value, the more significant the impact of the detected differential metabolites on this pathway. Impact: The metabolic pathway impact value. The larger it is, the greater the impact of the differential metabolites detected this time on the target pathway.

[0089] The detection method for intestinal epithelial transcriptome sequencing is as follows:

[0090] The intestinal tissue was sampled at -80℃ for the experiment. The experimental process was to enrich the mRNA with polyA structure in the total RNA by Oligo (dT) magnetic beads, and the RNA was broken into fragments of about 300bp in length by ion shearing. The fragments of 300bp in length were selected because the length of the linker is fixed. If the length of the broken fragment is shorter, the proportion of the linker sequence will be higher, thereby reducing the proportion of valid data; if the length of the broken fragment is longer, it will be unfavorable for the generation of clusters during the sequencing process. RNA was used as a template, 6-base random primers and reverse transcriptase were used to synthesize the first-chain cDNA, and the first-chain cDNA was used as a template for the synthesis of the second-chain cDNA, and the library construction was completed. After that, PCR amplification was used to enrich the library fragments, and then the library was selected according to the fragment size. The library size was 450bp. Then, the library was quality checked by Agilent2100Bioanalyzer, and the total concentration and effective concentration of the library were tested. Then, according to the effective concentration of the library and the amount of data required for the library, the libraries containing different Index sequences (each sample was added with a different Index, and finally the offline data of each sample was distinguished according to the Index) were mixed in proportion. The mixed library was uniformly diluted to 2nM and formed into a single-stranded library through alkaline denaturation. After RNA extraction, purification and library construction, the samples were subjected to paired-end (PE) sequencing using the second-generation sequencing technology (Next-Generation Sequencing, NGS) based on the Illumina sequencing platform. The experimental flow chart is shown below. Figure 21 The analysis process includes:

[0091] 1. Data collation: After the samples are sequenced on the machine, image files are obtained, which are converted by the software provided by the sequencing platform to generate the raw data (Raw Data) of FASTQ, i.e. the offline data. Then, the offline data (Raw Data) of each sample is statistically analyzed, including the sample name, Q30, the percentage of ambiguous bases, and Q20 (%) and Q30 (%). Q30 (bp): the total number of bases with a base recognition accuracy of more than 99.9%; N (%): the percentage of ambiguous bases; Q20 (%): the percentage of bases with a base recognition accuracy of more than 99%; Q30 (%): the percentage of bases with a base recognition accuracy of more than 99%;

[0092] 2. Data filtering: Sequencing data contains some reads with adapters and low quality. These sequences will cause great interference to subsequent information analysis. Therefore, the sequencing data needs to be further filtered. The data filtering criteria mainly include: 1) using Fastp to remove sequences with adapters at the 3' end; 2) removing reads with an average quality score lower than Q20;

[0093] 3. Data quality assessment: including base quality distribution, base content distribution, and average quality distribution of Reads, which is used to evaluate data quality;

[0094] 4. Comparative analysis:

[0095] (1) Basic statistics of alignment results: Use the upgraded HISAT2 (http: / / ccb.jhu.edu / software / hisat2 / index.shtml) of TopHat2 to align the filtered Reads to the reference genome. HISAT2 uses an improved BWT algorithm (Sirén et al., 2014) with faster speed and less resource consumption. When aligning with HISAT2, for non-strand-specific libraries, default parameters are used. For strand-specific libraries, the library type needs to be specified (i.e., --rna-strandness RF for first and --rna-strandness FR for second). If the reference genome is selected appropriately and there is no contamination in the relevant experiment, the mapping ratio of the sequencing sequences generated by the experiment is generally higher than 70%;

[0096] (2) Statistical analysis of alignment region distribution: Statistically analyze the distribution of Reads aligned to the genome. The positioning regions are divided into CDS (coding region), Intron (intron), Intergenic (intergenic region), and UTR (5' and 3' untranslated regions). In species with relatively complete genome annotation, usually, the content of Reads aligned to CDS (coding region) is the highest. The Reads aligned to the Intron (intron) region come from the remnants of pre-mRNA or are caused by intron retention events occurring during alternative splicing. The Reads aligned to the Intergenic (intergenic region) may be transcribed from new genes or new non-coding RNAs;

[0097] (3) Gene coverage uniformity: The distribution of the coverage of sequencing Reads on genes, which shows the sequence coverage of all genes of each sample from the 5' to 3' regions, and is used to evaluate the uniformity (or bias) of the sequencing results. Under ideal conditions, the distribution of Reads on all expressed genes should show a uniform distribution; its analysis flow chart is as Figure 22 shown.

[0098] Use the DESeq software to perform differential analysis of gene expression. The conditions for screening differentially expressed genes are: the fold change of expression |log2FoldChange| > 1, and the significance P-value < 0.05. The obtained bar chart of differential expression detection for the two samples is as Figure 23As shown, where the abscissa represents the comparison groups for differential analysis, the ordinate represents the number of differential genes, and in the color, red represents up-regulated genes and green represents down-regulated genes. The volcano plot of differentially expressed genes was drawn using the ggplots2 package in R language. The volcano plot shows the gene distribution, the fold change difference in gene expression, and the significance results. Under normal conditions, the distribution of differential genes on the left and right of this figure should be roughly symmetric. The left side shows the down-regulated genes in Case compared to Control, and the right side shows the up-regulated genes in Case compared to Control. The volcano plot of differential expression detection for the two samples obtained is as Figure 24 shown, where the abscissa is the value of the logarithm to the base 2 of the product of the gene expression levels of the two samples, that is, log2(A*B), where A and B represent the expression levels of the gene in the two samples respectively, and the ordinate is the value of the logarithm to the base 2 of the quotient of the expression levels, that is, log2(A / B). Red dots represent the up-regulated genes in this group, blue dots represent the down-regulated genes in this group, and gray dots represent genes with non-significant differential expression. Control: control group, Treat: experimental group, log2FoldChange: the value of the logarithm to the base 2 of the fold change in expression difference. Cluster analysis is used to judge the expression patterns of differentially expressed genes under different experimental conditions; genes with high expression correlation between samples are grouped together. Usually, these genes have actual connections in certain biological processes, or a certain metabolic or signaling pathway. Therefore, through expression clustering, we can discover unknown biological connections between genes. The Pheatmap package in R language was used to perform two-way cluster analysis on the union of differential genes and samples for all comparison groups. Clustering was performed according to the expression levels of the same gene in different samples and the expression patterns of different genes in the same sample. The Euclidean method was used to calculate the distance, and the hierarchical clustering longest distance method (Complete Linkage) was used for clustering. The obtained cluster analysis figure is as Figure 25 shown, where the horizontal direction represents genes, each column represents a sample, red represents highly expressed genes, and green represents lowly expressed genes. The topGO software was used for GO enrichment analysis. During the analysis, the differential genes annotated with GO term were used to calculate the gene list and the number of genes for each term, and then the P-value was calculated by the hypergeometric distribution method (the standard for significant enrichment is P-value < 0.05). The GO terms in which the differential genes are significantly enriched compared to the whole genome background were found, so as to determine the main biological functions performed by the differential genes. The GO enrichment analysis results of the differentially expressed genes were classified according to molecular function MF, biological process BP, and cellular component CC. The top 10 GO term entries with the smallest p-value, that is, the most significantly enriched, were selected for each GO classification for display. The obtained GO enrichment analysis bar chart is as Figure 26As shown, where the abscissa is the term of Go level2, and the ordinate is the -log10(p-value) enriched for each term. GO_Term: the enriched GO entry, Category: the classification where the enriched GO Term is located, GO_Term: the enriched GO entry. According to the GO enrichment results, the degree of enrichment is measured by the Rich factor, FDR value, and the number of genes enriched on this GO Term. The obtained GO enrichment analysis bubble chart is as Figure 27 shown, where the Rich factor refers to the ratio of the number of differentially expressed genes enriched in this GO Term to the number of differentially expressed genes annotated. The larger the Rich factor, the greater the degree of enrichment. The FDR generally ranges from 0 to 1. The closer it is to zero, the more significant the enrichment. The top 20 GO Term entries with the smallest FDR value, that is, the most significantly enriched, are selected for display. FDR: the corrected P value, that is, a more stringent P value. According to the KEGG enrichment analysis results of the differentially expressed genes, the top 30 Pathways with the smallest p-value, that is, the most significantly enriched, are selected for display. The obtained KEGG enrichment analysis bar chart is as Figure 28 shown. According to the KEGG enrichment results, the degree of enrichment is measured by the Rich factor, FDR value, and the number of genes enriched on this pathway. The obtained enrichment analysis bubble chart is as Figure 29 shown, where the Rich factor refers to the ratio of the number of differentially expressed genes enriched in this pathway to the number of differentially expressed genes annotated. The larger the Rich factor, the greater the degree of enrichment. The FDR generally ranges from 0 to 1. The closer it is to zero, the more significant the enrichment. The top 20 KEGG pathways with the smallest FDR value, that is, the most significantly enriched, are selected for display.

[0099] The detection method for 16S rRNA sequencing of intestinal contents is as follows:

[0100] Take the intestinal contents stored at -80°C during sampling for experiments. The specific experimental procedures are as follows:

[0101] 1. Extraction of total DNA from the microbiome: For microbiome samples from various sources, based on past project experience, select the most suitable method for extracting total DNA. At the same time, use Nanodrop to quantify the DNA and detect the quality of DNA extraction by 1.2% agarose gel electrophoresis;

[0102] 2. PCR amplification of target fragments: Usually, target sequences such as microbial ribosomal RNA or specific gene fragments that can reflect the composition and diversity of the microbial community are used as targets. Corresponding primers are designed based on the conserved regions in the sequences, and sample-specific Barcode sequences are added. Then, the variable regions (single or multiple consecutive ones) of the rRNA gene or specific gene fragments are amplified by PCR. The Pfu high-fidelity DNA polymerase from TransGen Biotech is used for PCR amplification, and the number of amplification cycles is strictly controlled to keep it as low as possible while ensuring the same amplification conditions for the same batch of samples. At the same time, negative controls are set up. The negative controls can detect microbial contamination in the environment, reagents, etc. Any sample group with bands in the negative control amplification cannot be used for subsequent experiments;

[0103] 3. Magnetic bead purification and recovery of amplification products:

[0104] (1) Add magnetic beads (Vazyme VAHTSTM DNA CleanBeads) with a volume 0.8 times that of the 25 μl PCR product. After shaking well to suspend, adsorb on the magnetic stand for 5 min, and carefully aspirate the supernatant with a pipette;

[0105] (2) Add 20 μl of 0.8-fold magnetic bead washing solution. After shaking well to suspend, place on the magnetic stand and adsorb for 5 min, then carefully aspirate the supernatant;

[0106] (3) Add 200 μl of 80% ethanol, place it upside down on the magnetic stand to adsorb the magnetic beads to the other side of the PCR tube. After full adsorption, aspirate the supernatant;

[0107] (4) Let it stand at room temperature for 5 min until the alcohol has completely evaporated and the magnetic beads show cracks;

[0108] (5) Add 25 μl of Elution Buffer for elution;

[0109] (6) Place the PCR tube on the adsorption rack for 5 min for full adsorption, then transfer the supernatant to a clean 1.5 ml centrifuge tube for storage;

[0110] 4. Fluorescent quantification of amplification products: Fluorescent quantification is performed on the PCR amplification and recovery products. The fluorescent reagent is Quant-iT PicoGreen dsDNA Assay Kit, and the quantification instrument is a Microplate reader (BioTek, FLx800). According to the fluorescent quantification results, each sample is mixed in the corresponding proportion according to the sequencing amount requirement of each sample;

[0111] 5. Preparation of sequencing library: The TruSeq Nano DNA LT Library Prep Kit is used to prepare the sequencing library;

[0112] (1) First, perform end repair on the above amplification products. Use End Repair Mix 2 in the kit to excise the protruding bases at the 5' end of the DNA sequence, add a phosphate group simultaneously, and fill in the missing bases at the 3' end;

[0113] (2) Add an A base to the 3' end of the DNA sequence to prevent self-ligation of DNA fragments and ensure that the target sequence can be ligated to the sequencing adapter (there is a protruding T base at the 3' end of the sequencing adapter);

[0114] (3) Add a sequencing adapter containing a library-specific tag (i.e., Index sequence) to the 5' end of the sequence so that the DNA molecule can be immobilized on the Flow Cell;

[0115] (4) Use BECKMAN AMPure XP Beads to remove adapter self-ligated fragments through magnetic bead screening and purify the library system after adding the adapter;

[0116] (5) Perform PCR amplification on the DNA fragments ligated with the adapter above to enrich the sequencing library template, and use BECKMAN AMPure XP Beads to purify the library enrichment product again;

[0117] (6) Perform final fragment selection and purification on the library by 2% agarose gel electrophoresis;

[0118] 6. Perform high-throughput sequencing on the machine:

[0119] (1) Before sequencing on the machine, it is necessary to first perform quality inspection on the library on the Agilent Bioanalyzer using the Agilent High Sensitivity DNA Kit. A qualified library has only a single peak and no adapter;

[0120] (2) Then, use the Quant-iT PicoGreen dsDNA Assay Kit to quantify the library on the Promega QuantiFluor fluorescence quantification system. The concentration of a qualified library should be above 2 nM;

[0121] (3) After gradient dilution of the qualified sequencing libraries for each run (the Index sequences cannot be repeated), mix them in the corresponding proportion according to the required sequencing amount, and denature them with NaOH into single strands for sequencing on the machine;

[0122] (4) If paired-end sequencing is performed using the MiSeq sequencer, the corresponding reagent is the MiSeq Reagent Kit V3 (600 cycles); if paired-end sequencing is performed using the NovaSeq sequencer, the corresponding reagent is the NovaSeq 6000SP Reagent Kit (500 cycles); due to the short read length of the MiSeq, and to ensure sequencing quality, it is recommended that the optimal sequencing length of the target fragment be 200 - 450 bp. Ribosomal RNA contains multiple conserved regions and highly variable regions. Usually, we design primers based on the conserved regions to amplify single or multiple variable regions of the rRNA gene, and then sequence and analyze microbial diversity. Due to the limitation of the MiSeq read length and to ensure sequencing quality, the optimal range of the inserted fragment for sequencing is 200 - 450 bp.

[0123] According to Figure 30 Perform the following analysis process: First, preliminarily screen the original off-machine data of high-throughput sequencing according to sequence quality; re-sequence and supplement the problem samples; divide the libraries and samples according to the index and Barcode information for the original sequences that pass the quality preliminary screening, and remove the barcode sequences; perform sequence denoising or OTU clustering according to the QIIME2 dada2 analysis process or the analysis process of the Vsearch software; display the specific composition of each sample (group) at different species taxonomic levels to understand the overall situation; evaluate the Alpha diversity level of each sample based on the distribution of ASV / OTU in different samples, and reflect whether the sequencing depth is appropriate through the rarefaction curve; at the ASV / OTU level, calculate the distance matrix of each sample, and through various unsupervised sorting and clustering methods, combined with the corresponding statistical test methods, measure the beta diversity difference and difference significance between different samples (groups); at the species taxonomic composition level, through various unsupervised, supervised sorting, clustering and modeling methods, combined with the corresponding statistical test methods, further measure the difference in species abundance composition between different samples (groups), and try to find biomarker species; construct an association network based on the composition distribution of species in each sample, calculate the topological index, and try to find key species; based on the sequencing results of 16S rRNA, 18S rRNA and ITS genes, the microbial metabolic functions of the samples can also be predicted, differential pathways can be found, and the species composition of specific pathways can be obtained; based on the above results, draw charts at the level suitable for paper publication and perform statistical test analysis.

[0124] The bioinformatics analysis methods are as follows:

[0125] 1. DADA2 sequence denoising

[0126] Analysis software: QIIME2 (2019.4);

[0127] Analysis steps: First, call qiime cutadapt trim-paired to excise the primer fragments of the sequences and discard the sequences that do not match the primers; then call DADA2 through qiime dada2 denoise-paired for quality control, denoising, splicing, and chimera removal. The above steps are analyzed separately for each library. After denoising all libraries, merge the ASVs feature sequences and the ASV table, and remove the singletons ASVs (i.e., ASVs with a total sequence count of only 1 in all samples, which is the default operation).

[0128] 2. Statistical analysis of sequence length distribution

[0129] Use R language scripts to statistically analyze the length distribution of high-quality sequences contained in all samples.

[0130] 3. Taxonomic annotation of species

[0131] Analysis software: QIIME2 (2019.4);

[0132] Database: 1) For the 16S rRNA gene of bacteria or archaea, the Greengenes database (Release 13.8, http: / / greengenes.secondgenome.com / ) (DeSantis et al., 2006) is defaultly selected, and the Silva database (Release 132, http: / / www.arb-silva.de) (Quast et al., 2013) can also be selected; 2) For the 18S rRNA gene of eukaryotic microorganisms, the Silva database (Release 132) is defaultly selected; 3) For fungal ITS sequences, the UNITE database (Release 8.0, https: / / unite.ut.ee / ) (Koljalg et al., 2013) is defaultly selected; 4) For functional genes or other requirements, we use the localized nt (downloaded in August 2019, ftp: / / ftp.ncbi.nih.gov / blast / db / ) database for annotation.

[0133] Analysis steps: (1) For the first three types of databases, the classify-sklearn algorithm of QIIME2 (Bokulich et al., 2018) (https: / / github.com / QIIME2 / q2-feature-classifier) is used: For the feature sequences of each ASV or the representative sequences of each OTU, in the QIIME2 software with default parameters, a pre-trained Naive Bayes classifier is used for species annotation. (2) For the nt database, the BROCC algorithm (Nilsson et al., 2006) is used: First, blastn is used to align it with the nt database (or specific sequences screened from it); then the brocc.py script is called to obtain annotation information according to the recommended parameters.

[0134] 4. Construct a phylogenetic tree

[0135] Analysis software: QIIME2 (2019.4);

[0136] Analysis steps: Use the analysis process of "qiime phylogeny align-to-tree-mafft-fasttree" to call mafft (Katoh, 2002) for multiple sequence alignment, mask the parts without phylogenetic information, and then call FastTree (Price et al., 2009) to construct a phylogenetic tree and generate a tree file.

[0137] 5. Flatten the ASV / OTU table

[0138] In the previous analysis steps, an abundance table of ASV / OTU has been generated, and some subsequent analysis steps need to be carried out at the same sequencing depth level for each sample. Therefore, certain transformation processing is required for this table. The rarefaction method can be used. It randomly extracts a certain number of sequences from each sample to reach a unified depth, so as to predict the ASVs or OTUs and their relative abundances that can be observed in each sample at this sequencing depth (Heck et al., 1975; Kemp and Aller, 2004). Therefore, this process is also called flattening.

[0139] Analysis software: QIIME2 (2019.4);

[0140] Analysis steps: Use the qiime feature-table rarefy function, and set the flattening depth to 95% of the lowest sample sequence volume.

[0141] Analysis software used: QIIME2 (2019.4); ggplot2 package in R language, etc. By calling the "qiime taxa barplot" command and using the feature table after removing singletons, the compositional distribution of each sample at six taxonomic levels of phylum, class, order, family, genus, and species was visualized, and the analysis results were presented as a bar chart, obtaining a taxonomic composition analysis chart as shown in Figure 31 where phylum represents the phylum level, genus represents the genus level, and Relative Abundance represents the relative abundance.

[0142] Alpha diversity refers to the indicators of species richness, diversity, and evenness in a local homogeneous habitat, and is also known as within-habitat diversity. To comprehensively evaluate the alpha diversity of the microbial community, in this process, the Chao1 (Chao, 1984) and Observed species indices were used to characterize richness, the Shannon (Shannon, 1948a, b) and Simpson (Simpson, 1949) indices were used to characterize diversity, the Faith’s PD (Faith, 1992) index was used to characterize the diversity based on evolution, the Pielou’s evenness (Pielou, 1966) index was used to characterize evenness, and the Good’s coverage (Good, 1953) index was used to characterize coverage. Analysis software: QIIME2 (2019.4) was used. Using the unrarefied ASV / OTU table, the "qiime diversity alpha-rarefaction" command was called, and the parameters "--p-steps 10 --p-min-depth 10 --p-iterations 10" were set, that is, the minimum rarefaction depth was 10, and the parameter "--p-max-depth" was set to 95% of the sequence volume of the sample with the lowest sequencing depth among all samples. Then, 10 depth values were evenly selected between this depth and the minimum depth, and each depth value was rarefied 10 times to calculate the above indices. By default, the average score at the maximum rarefaction depth was selected as the alpha diversity index, and the obtained Alpha diversity index chart is as shown in Figure 32As shown, where each panel corresponds to an alpha diversity index, which is identified in the gray area at its top. In each panel, the abscissa is the grouping label, and the ordinate is the value of the corresponding alpha diversity index. In the box plot, the meanings of each symbol are as follows: the upper and lower end lines of the box are the upper and lower quartiles (Interquartile range, IQR); the median line is the median; the upper and lower edges are the maximum and minimum inner fence values (1.5 times the IQR); the points outside the upper and lower edges represent outliers. Panel: group, Chao1 richness estimator index: first proposed by Chao, which estimates the actual number of species in a community by calculating the number of ASVs / OTUs detected only once and twice in the community (i.e., "Singleton" and "Doubleton"), Shannon diversity index: comprehensively considers the richness and evenness of the community.

[0143] Beta diversity refers to the dissimilarity in species composition between different communities along an environmental gradient or the turnover rate of species along an environmental gradient. Therefore, it is also known as between-habitat diversity. Principal coordinates analysis (PCoA) is one of the most classic non-constrained ordination (Classical Multidimensional Scaling, cMDScale) analysis methods (Ramette, 2007). It unfolds the sample distance matrix after projection in a low-dimensional space and maximally preserves the distance relationship of the original samples. PCoA considers the sample distance as a whole and is more in line with the characteristics of ecological data compared to principal components analysis (PCA). Therefore, as an ordination analysis method, it is more recommended. Analysis software used: QIIME2 (2019.4). Using the rarefied ASV / OTU table, call the "qiime diversity core-metrics-phylogenetic" command to calculate four dissimilarity distance matrices and perform PCoA analysis on these distance matrices, output the QZV file. The resulting Beta diversity analysis (PCoA analysis, weighted UniFrac distance) figure is as Figure 33As shown, each point in the figure represents a sample, and points of different colors indicate different samples (groups). The percentages in the parentheses of the coordinate axes represent the proportion of the sample difference data (distance matrix) that the corresponding coordinate axes can explain. An oval dashed circle is provided in the figure, which is a 95% confidence ellipse (that is, 95 out of 100 samples in this sample group will fall within it). Unweighted UniFrac: It is an analytical method used to compare the differences in microbial communities of environmental samples. It utilizes the evolutionary information between sample sequences to calculate the Unifrac distance matrix between samples for evaluating beta diversity. Unweighted UniFrac measures the differences between different environmental samples based on the lengths of the constructed evolutionary branches, and the differences are represented by 0-1 distance values. The distance between the earliest-differentiated branches on the evolutionary tree is 1, indicating the greatest difference. Weighted UniFrac based PCoA: Based on unweighted UniFrac, it also considers the differences in species abundances. It adds abundance weights to each evolutionary branch, thereby more comprehensively reflecting the degree of difference between community samples in a quantitative manner.

[0144] To study which species are common and which are unique among different samples (groups), a Venn diagram (Venn; https: / / en.wikipedia.org / wiki / Venn_diagram) is used for community analysis. The Venn diagram is made using the ASV / OTU abundance table. The number of members in each set is counted according to their presence or absence among different samples (groups), that is, the number of ASV / OTUs unique to each sample (group) and the number of ASV / OTUs common among samples (note that it is not the abundance value). The analysis software used: R script, VennDiagram package or plotrix package. The Venn diagram of species difference analysis and marker species ASV / OTU is as Figure 34 shown, where each ellipse represents a sample (group), and the overlapping area between ellipses indicates the common ASV / OTUs among samples (groups). The number in each block indicates the number of ASV / OTUs contained in that block.

[0145] To further compare the differences in species composition between samples and display the distribution trend of species abundances in each sample, a heatmap can be used for species composition analysis. In the study, the taxonomic unit composition at the genus level is generally used as the analysis object, so the abundance data of the top 50 genera in terms of average abundance are default used to draw the heatmap. The horizontal and vertical coordinates of the heatmap can be arranged in a specific order, such as sorted by the average abundance of taxa, the sampling time of samples, etc.; or a clustering tree can be drawn for sorting based on the correlation between taxa or samples, that is, a clustered heatmap is drawn. The analysis software used: R script, pheatmap package. The species composition heatmap obtained by analysis is as Figure 35As shown, in the graph of sample clustering, the samples are clustered by UPGMA according to the Euclidean distance of species composition data and arranged according to the clustering results; otherwise, they are arranged according to the grouping or default order of the samples. In the graph of species clustering, by default, the species are clustered by UPGMA according to the Pearson correlation coefficient matrix of their composition data and arranged according to the clustering results; otherwise, they are sorted according to the average abundance of the species in the samples. In the graph, the red color block represents that the abundance of the genus in this sample is higher than that in other samples, and the blue color block represents that the abundance of the genus in this sample is lower than that in other samples. Pearson correlation coefficient: In statistics, the Pearson correlation coefficient is used to measure the correlation (linear correlation) between two variables X and Y.

[0146] The obtained heatmap of species composition (between groups) is as Figure 36 shown, where the red color block represents that the abundance of the genus in this group is higher than that in other groups, and the blue color block represents that the abundance of the genus in this group is lower than that in other groups.

[0147] LEfSe (LDA Effect Size) analysis is an analytical method that combines non-parametric Kruskal-Wallis and Wilcoxon rank-sum tests with the effect size of linear discriminant analysis (LDA). The analysis results of LEfSe include three parts, namely, the bar chart of the LDA value distribution of significantly different species, which is used to show the significantly enriched species (note that significantly down-regulated ones are not shown) in each group and their degree of importance; the cladogram of species taxonomy, which is used to show the taxonomic hierarchical distribution of the marker species in each group of samples; the bar chart of abundance between groups, which is used to show the specific distribution of the marker species in different grouped samples; these displays are more suitable for samples with relatively simple species composition and appropriate number of different species, such as intestinal samples of animals and humans, etc. The analysis software used: Python LEfse package, Rggtree, ggplot2 package, etc., or the Galaxy online analysis platform (http: / / huttenhower.sph.harvard.edu / galaxy / ) is used for analysis. The obtained LEfSe analysis (bar chart of the LDA effect value of the marker species) is as Figure 37 shown, where the vertical axis is the taxonomic units with significant differences between groups, and the horizontal axis intuitively shows the logarithmic score values of the LDA analysis of each taxonomic unit in the form of a bar chart; the taxonomic units are sorted according to the score values to describe their specificity in sample grouping; the longer the length, the more significant the difference of the taxonomic unit, and the color of the bar chart indicates the sample group with the highest abundance corresponding to the taxonomic unit.

[0148] The foregoing has described the preferred embodiments of the present invention in detail. However, the present invention is not limited to the described embodiments. Those skilled in the art can make various equivalent deformations or substitutions without departing from the spirit of the present invention, and these equivalent deformations or substitutions are all included within the scope defined by the claims of this application.

Claims

1. A method for analyzing the intestinal health of crispy tilapia, characterized in that, It includes the following steps: (1) After domesticating tilapia for half a month, it is divided into a control group and an experimental group. The control group is fed with commercial tilapia feed, and the experimental group is fed with crisp tilapia feed; (2) At the 30th day, 60th day, 90th day, and 120th day respectively, samples are selected from the crisp group and the control group for muscle tissue sectioning, texture analysis, and intestinal tissue sectioning; (3) At the 120th day, samples are selected from the crisp group and the control group respectively for detecting intestinal antioxidant enzyme activity, detecting intestinal anti-inflammatory factor gene Q-PCR, detecting intestinal ELISA protein, performing 16S rRNA sequencing on intestinal contents, detecting serum metabolome, and performing intestinal epithelial transcriptome sequencing; (4) Dynamic crispness degree analysis: Understand the muscle quality through muscle tissue sectioning and texture analysis, and preliminarily understand the degree of intestinal damage through intestinal tissue sectioning; (5) Three-dimensional evaluation of intestinal health: Conduct intestinal morphological evaluation through the villus height and crypt depth of intestinal sections, analyze the intestinal flora structure through α-diversity, β-diversity, and function prediction in 16S rRNA sequencing of intestinal contents, and conduct host gene analysis through screening pathway enrichment analysis of differentially expressed genes in intestinal epithelial transcriptome sequencing; (6) Cross-omics integration analysis: Conduct joint analysis through serum metabolome detection - host gene analysis of intestinal epithelium, serum metabolome detection - intestinal flora structure analysis, integrate the relationship between intestinal genes and microorganisms and the overall metabolic level of fish, and conduct cross-omics analysis of the intestinal health of crisp tilapia.

2. The method for analyzing the intestinal health of crispy tilapia according to claim 1, wherein In step (1), the weight of the tilapia is 350 - 450 g; each of the control group and the experimental group has 3 parallels, with 30 tilapia in each parallel; the number of samples is 9 tilapia in each of the control group and the experimental group.

3. The method for analyzing the intestinal health of crispy tilapia according to claim 1, characterized in that, In step (2), the operation and analysis method of the muscle tissue sectioning include the following steps: material sampling, dehydration, clearing, wax infiltration, embedding, sectioning, baking, and HE staining; the texture analysis is measured using a Universal TA texture analyzer; the operation and analysis method of the intestinal tissue sectioning include the following steps: material sampling, dehydration, clearing, wax infiltration, embedding, sectioning, baking, and HE staining.

4. The analysis method for the intestinal health of crispy tilapia according to claim 1, characterized in that, In step (3), the detection of intestinal antioxidant enzyme activity is measured using a kit from Nanjing Jiancheng Bioengineering Institute.

5. The analysis method for the intestinal health of crispy tilapia according to claim 1, wherein, In step (3), the method for detecting intestinal anti-inflammatory factor gene Q-PCR includes the following steps: extracting total RNA from tilapia intestine, quantifying the RNA concentration using a spectrophotometer, removing genomic DNA, reverse transcribing it into cDNA using a reverse transcription kit, and finally performing real-time fluorescence quantitative PCR detection.

6. The method for analyzing the intestinal health of crispy tilapia according to claim 1, characterized in that In step (3), the detection of intestinal ELISA protein is measured using a kit from Jiangsu Enzyme Immunoassay Industry Co., Ltd.

7. The analysis method for the intestinal health of crispy tilapia according to claim 1, wherein, In step (3), the method for 16S rRNA sequencing of intestinal contents includes the following steps: extracting total DNA of the microbiome, PCR amplification of the target fragment, magnetic bead purification and recovery of the amplification product, fluorescence quantification of the amplification product, preparation of the sequencing library, and high-throughput sequencing on the machine.

8. The analysis method for the intestinal health of crispy tilapia according to claim 1, characterized in that In step (3), the serum metabolome detection is performed by Suzhou Panomics Biomedical Technology Co., Ltd. for metabolome determination.

9. The method for analyzing the intestinal health of crispy tilapia according to claim 1, wherein In step (3), the method for intestinal epithelial transcriptome sequencing includes the following steps: extraction of total RNA, detection of total RNA quality, purification of mRNA, fragmentation of mRNA, cDNA synthesis, PCR enrichment of library fragments, library quality inspection, and Illumina platform-based on-machine sequencing.

Images