Changes in Intestinal Microflora Diversity of Diarrhea Calves among Different Seasons using Metagenomic Sequencing Analysisor

The gut microbiota of calves is diverse and abundant, making it a complex and relatively balanced microecosystem (Lee et al., 2019). The microbial community has an impact on the development, nutrient digestion and absorption and immunity of the body’s intestines. When calves are affected by the external environment, their intestinal microbiota balance may be disrupted, leading to changes in the structure and quantity of microbiota. In addition, the morphology, structure and function of the intestinal mucosa in calves are not sound and calves with poor physical fitness are prone to diarrhea (Brooks et al., 1998). In terms of infection, it can occur in any season. It is common in winter and spring. Its incidence rate is different, but its mortality is high. The disease is currently occurring globally, widely distributed in countries such as the United States and the United Kingdom. With the continuous expansion and development of the cattle farming industry, the incidence of the disease is also on the rise, and the serious harm to the cattle farming industry is increasing (Kim et al., 2021).  Therefore, conducting research and analysis on gut microbiota is of great significance for the prevention and control of new calf infectious diseases and the risk assessment of potential zoonotic diseases.
       
With the development of high-throughput sequencing technology, significant breakthroughs have been made in the study of gut microbiota, including nutrition and metabolism between gut microbiota and hosts; Immune regulation is closely related to growth and development, as well as the health and diseases (Handelsman et al., 1998). It avoids the separation and cultivation of microorganisms in the sample, providing a way to study microorganisms that cannot be separated and cultivated and more realistically reflects the composition and interaction of microorganisms in the sample, while studying their metabolic pathways and gene functions at the molecular level (Tringe et al., 2005).
       
In recent years, some studies have shown that, the composition and diversity of animal gut microbiota vary with seasons (Elle et al., 2020; Yuanqiu Dong et al., 2019). By using 16S rRNA sequencing, metagenomic sequencing and bioinformatics analysis methods, the diversity, structure and function of gut microbiota in diarrhoeic calves from different seasons were analyzed (Redding et al., 2021). The effects of seasonal changes, gender differences, nutritional changes, and antibiotic changes on the composition and abundance of gut microbiota in diarrhea calves have been explored and the response mechanism of gut microbiota in diarrhea calves under seasonal changes through gene functional annotation have been analyzed. This study conducted diversity analysis on the gut microbiota of diarrhea calves using high-throughput sequencing technology, preliminarily exploring the impact of diarrhea on the changes of gut microbiota in calves, comparing the changes in gut microbiota structure in different seasons, and predicting the high incidence period of potential diseases. This provides a theoretical basis for improving the prevention and control strategies for calf diseases and the level of feeding management.  

Sample collection and preservation
 
From January to December 2022, 788 fecal samples of diarrhoeic calves within 30 days of age were collected from 38 different scale dairy farms in 13 counties of Tangshan City, Hebei Province. Marking of the calf number, sampling date and farm name on the corresponding sterile centrifuge tube was done, 30-200 g of fecal samples were collected and stored in liquid nitrogen for future use. After all samples were collected, they were divided into four groups, namely the spring group (F1), summer group (C1), autumn group (C2) and winter group (C3). They were stored in liquid nitrogen and sent to Beijing Nuohe Zhiyuan Biological Information Technology Co., Ltd. for sequencing. The specific grouping of sampling information is shown in Table 1.
 

 
 PCR amplification and product mixing and purification
 
The V4 region was amplified using universal primers 515F (5' GTGCCAGCMGCCGCGGTAA-3') and 806R (5'-GGACT ACVSGGGTATCTAAT-3') targeting the 16S RNA gene of bacteria. Using specific primers with Barcode, Phusion from New England Biolabs® High Fidelity PCR Master Mix with GC Buffer and efficient high fidelity enzymes for PCR. The reaction system is 30 µL: Phusion Master Mix (2X) 15 µL, ultrapure water 2 µL, upstream and downstream primers 1.5 µL each and template DNA 10 µL. PCR amplification conditions: 98°C 3 min, (98°C 10 s, 50°C 30 s, 72°C 30 s), 30 cycles, 72°C 5 min.Mix the same amount of samples according to the concentration of  PCR products. After fully mixing, 2% agarose gel electrophoresis was used to detect the PCR products. For the target strip, use the gel recovery kit provided by Qiagen Company to recover the products.
 
High throughput sequencing and data analysis
 
After the sequencing of all samples is completed, the original data will be obtained by splicing the double ended data using Illumina PE150 technology for quality control and chimeric filtering, filtering out sequences with abundance less than 5 (Li et al., 2020; Callahan et al., 2017) and obtaining the final ASVs (Amplicon sequence variations). Use the concept of ASVs to construct a class OTU table, use the RDP classifier Bayesian algorithm to perform taxonomic analysis on OTU representative sequences and calculate the community composition of each sample at the level of phylum, family, genus and species. Obtain an OTU distribution table based on the number of sequences in each OTU, analyze species abundance (Avershina et al., 2013; Rivas et al., 2013; Jiang et al., 2013; Edgar, 2004; Edgar, 2013) and obtain the final feature table and sequence Callahan et al., (2017). Diversity analysis, species classification annotation, difference analysis, functional analysis, etc. will be conducted by the biological analysis platform of Beijing Nuohe Zhiyuan Biological Information Technology Co., Ltd.
 
Data processing
 
Use statistical analysis methods such as T-test, MetaStat, and LEfSe to test the significance of differences in species composition and community structure of grouped samples. The data is presented in the form of mean ± standard error (Mean ± SD), with P<0.05 as the significant difference and P<0.01 as the extremely significant difference.

Quality evaluation of microbial sequencing in calf diarrhea feces
 
Use SOAP denovo and CD-HIT software to assemble and de redundant data, use MetaGeneMark for ORF (Open Reading Frame) prediction and filter out information with a length of less than 100 nt based on the prediction results. After sequencing, effective sequences were obtained through quality control, chimerism removal and splicing. A total of 698,029 original sequences were predicted, with an average of 174,507 sequences per sample. After filtering through relevant software, 531,393 valid sequences were ultimately obtained. The average length is 693.82 bp and the GC content is 46.69%. Among them, there are 268,345 complete genes, accounting for 50.5% of the total number of non redundant genes. In addition, the effective data rate for quality control of calf diarrhea fecal samples is 99.65% (Table 2). The number of C1 characteristic sequences in calf diarrhea fecal samples is the highest and the F1 group is the lowest. The highest overlap is 93,221 feature sequences and the lowest is 15,346 feature sequences. The highest order of the four overlapping parts is 22,097 feature sequences and the lowest is 14,783 feature sequences (Fig 1).
 

 

       
The original data volumes for F1, C1, C2 and C3 are 7,802.78, 7,063.78, 7,356.44 and 6212.20, respectively. The final valid data obtained by filtering F1, C1, C2 and C3 through relevant software are 7, 782.27, 7045.09, 7313.78 and 6195.05, respectively. The Q30 of F1, C1, C2 and C3 filtered data were 92.36, 92.84, 93.19 and 93.00, respectively. The GC contents (%) of F1, C1, C2 and C3 filtered data were 45.67, 45.61, 46.35 and 44.98, respectively. The percentages of effective data for F1, C1, C2, and C3 compared to the original data are 99.737%, 99.735%, 99.420% and 99.724%, respectively, indicating that the sequencing data is reliable and of good quality (Table 3).
 

 
Beta diversity analysis
 
Perform PCoA (Principal Co coordinates Analysis) analysis based on Bray Curtis distance and select the principal coordinate combination with the highest contribution rate for graphical display. Observing the differences between sample groups through PCoA, different colors in the results represent different groups. The closer the sample is, the more similar the microbial composition and structure between the samples are. Conversely, the greater the difference. The distance between the C2 and C3 treatment groups is relatively close, while the distance between the C1 group and the other three groups is relatively far, indicating significant differences in microbial community types (Fig 2).
 

       
In order to study the similarity of different samples, cluster analysis can also be conducted on the samples to construct a cluster tree of the samples. Non metric multidimensional scale (NMDS) analysis based on species abundance at the genus level can also fully support the above analysis results. If the species composition of the samples is more similar, the distance between them in the NMDS diagram is closer (Fig 3).
 

 
Species composition analysis
 
The number of ORFs that can be annotated into the NR database is 465,393 (87.58%) based on 531,393 predicted genes after original de redundancy. Based on the annotation results of feature sequences and the characteristic tables of each sample, a species abundance table at the level of kingdom, phylum, class, order, family, genus and species was obtained. The species composition and inter group differences of the gut microbiota of four groups of samples were analyzed for different levels of species abundance tables. Among them, the ORFs that can be annotated into the NR database have a proportion of 86.45% at the boundary level, 83.53% at the phylum level, 78.59% at the class level, 77.95% at the order level, 62.78% at the family level, 57.58% at the genus level and 38.18% at the species level. The dominant phylum includes Firmicutes, Bacteroidetes, and Actinobacteria, among others. This study focuses on changes at the phylum, genus and species levels (Table 4).
 

 
Microbial composition at the phylum level
 
At the phylum level, species composition analysis yielded a total of 10 phyla. Among them, Firmicutes, Actinobacteria, and Bacteroidetes are the dominant phyla in F1, C2 and C3. The dominant phyla in C1 are Firmicutes, Actinobacteria, and Proteobacteria (Fig 4). Whether in summer or winter, Firmicutes and Actinobacteria have always been significant dominant bacteria, which is similar to other studies and consistent with the results of studies on the composition of gut microbiota in herbivorous animals (Fountain et al., 2020; Oikonomou et al., 2013). These two phyla are the main dominant phyla of ruminants, with the highest relative abundance. They participate in important processes such as food digestion, nutrient regulation and absorption, energy metabolism, and host gut defense against invasion of foreign pathogens (Wang et al., 2017). However, their relative abundance may vary in different seasons. Proteobacteria is widely present in nature and is a common opportunistic pathogen that can colonize the body’s skin, respiratory tract, gastrointestinal tract, etc. This study found that the abundance of Proteobacteria in the intestinal microbiota of winter diarrhea calves is higher than other seasons, indicating that the immune system of calves is more likely to decrease in winter, leading to diseases .  
 

 
 Microbial composition at the genus level
 
At the genus level, species composition analysis yielded a total of 10 genera. The main genera in F1 are Blautia (17.269%), Collins ella (13.88%), and Bacteroides (3.996%). The main bacterial genera in C1 are Bacteroides (13.804%), Prevotella (8.035%) and Clostridium (2.984%). The main genera of bacteria in C2 are Bacteroides (15.535%), Lactobacillus (10.581%) and Faecalibacterium (8.467%). The main genera of bacteria in C3 are Lactobacillus (18.939%), Faecalibacterium (6.679%) and Bacteroides (4.858%) (Fig 5).
 

 
Microbial composition at the species level
 
At the species level, a total of 10 species were obtained through species composition analysis. The main strains in F1 are Blautia sp. CAG:257 (11.709%), Clostridium hiranonis (7.248%) and Firmicutes bacterium CAG:424 (5.044%). The main strains in C1 are Firmicutes bacterium CAG:110 (3.647%), Prevotella sp. CAG:485 (2.787%) and Firmicutes bacterium CAG:424 (0.214%). The main strains in C2 are Faecalibacterium prausnitzii (3.723%), Lactobacillus reuteri (2.224%) and Firmicutes bacterium CAG:424 (0.126%). The main strains in C3 are Lactobacillus reuteri (4.901%), Firmicutes bacterium CAG: 424 (3.481%)and Faecalibacterium prausnitzii (3.084%) (Fig 6).
 

 
PICRUS function prediction analysis
 
Out of 531,393 predicted genes, 380,152 (71.54%) genes can be compared to the KEGG database. Among them, 206,460 (38.85%) genes can be compared to 4,818 KEGG ortholog groups in the database; 378,220 (71.18%) genes can be compared to the egg NOG database; There are 19,712 (3.71%) genes that can be compared to the CAZy database (Table 5).
 

 
KEGG function prediction analysis
 
Based on the functional annotations and abundance information of all samples in the KEEG database, select the top 35 abundant functions and their abundance information in each sample to draw a bar chart. Use the KO based relative abundance bar chart as an example to display (Fig 7). The KEEG database annotation results statistical table shows that amino acid transport metabolism, energy metabolism and other metabolic functions are abundant, as well as genetic functions such as protein modification, folding, and translation are abundant (Table 6). The results indicate that the microbial metabolism and genetic information processing functions of calf diarrhea feces are very rich.
 



 
CARD function prediction analysis
 
Out of 531,393 predicted genes after de redundancy, 307 genes were able to be compared to the CARD database, containing a total of 228 ARO (the Antibiotic Resistance Ontology). Among them, F1 is resistant to 117 drugs, and the top three resistance genes in abundance are tetO, tetW/N/W and mefA; C1 is resistant to 148 drugs, and the top three resistance genes in abundance are tetW/N/W, lnuC and ErmF; C2 is resistant to 170 drugs, and the top three resistance genes in abundance are APH3-Ib, tetW/N/W and OXA-85. C3 is resistant to 163 drugs and the top three resistance genes in abundance are tetO, tetW/N/W and AAC6-Ie-APH2-Ia (Fig 8).
 

       
Seasonal changes are often considered a factor that significantly affects microbial communities and antibiotic resistance genes in different environments (He et al., 2020). The seasonal changes in ruminant populations have a significant impact on the diffusion of antibiotic resistance genes. In this study, the abundance of tetW/N/W was relatively abundant in all four seasons.

The change of gut microbiota in newborn calves during diarrhea is a stable and slow process, and the experiment investigates the differences in gut microbiota in different seasons throughout the entire process. By annotating the gut microbiota species of diarrhoeic calves, it was found that the dominant phyla is Firmicutes and Actinobacteria, and the relative abundance of the phyla varies according to seasonal changes. Different seasons can have an impact on the diversity and function of the gut microbiota in diarrhea. The gut microbiota also plays a crucial role in the nutrient absorption, energy metabolism and immune system of calves.

We are grateful to all lab members for helpful discussions and suggestions in this research paper writing.

This research was supported by projects Tanshan Science and Technology Plan project (23150202 A); Scientific Research Fund of Tangshan Normal University (2023B24); Establishment of Rapid Detection Method for Pathogens of Important Epidemic Diseases in Dairy Cattle and Development of Reagent Kit (HX202207); Tangshan Experimental Station of the Innovation Team for Milk and Beef Cattle in the Modern Agricultural Industry Technology System of Hebei Province (HBCT2023180403).

Not applicable” for studies not involving humans or animals.

Not applicable.

We didn’t create the data.

The authors declare no conflict of interest.