Multi-omics strategy reveals potential role of antimicrobial resistance and virulence factor genes responsible for Simmental diarrheic calves caused by Escherichia coli

ABSTRACT Escherichia coli (E. coli) is reported to be an important pathogen associated with calf diarrhea. Antibiotic resistance genes (ARGs) and virulence factor genes (VFGs) pose a considerable threat to both animal and human health. However, little is known about the characterization of ARGs and VFGs presented in the gut microbiota of diarrheic calves caused by E. coli. In this study, we used multi-omics strategy to analyze the ARG and VFG profiles of Simmental calves with diarrhea caused by E. coli K99. We found that gut bacterial composition and their microbiome metabolic functions varied greatly in diarrheic calves compared to healthy calves. In total, 175 ARGs were identified, and diarrheal calves showed a significantly higher diversity and abundance of ARGs than healthy calves. Simmental calves with diarrhea showed higher association of VFGs with pili function, curli assembly, and ferrienterobactin transport of E. coli. Co-occurrence patterns based on Pearson correlation analysis revealed that E. coli had a highly significant (P < 0.0001) correlation coefficient (>0.8) with 16 ARGs and 7 VFGs. Metabolomics analysis showed that differentially expressed metabolites in Simmental calves with diarrhea displayed a high correlation with the aforementioned ARGs and VFGs. Phylotype analysis of E. coli genomes showed that the predominant phylogroup B1 in diarrheic Simmental calves was associated with 10 ARGs and 3 VFGs. These findings provide an overview of the diversity and abundance of the gut microbiota in diarrheic calves caused by E. coli and pave the way for further studies on the mechanisms of antibiotic resistance and virulence in the calves affected with diarrhea. IMPORTANCE Simmental is a well-recognized beef cattle breed worldwide. They also suffer significant economic losses due to diarrhea. In this study, fecal metagenomic analysis was applied to characterize the antibiotic resistance gene (ARG) and virulence factor gene (VFG) profiles of diarrheic Simmental calves. We identified key ARGs and VFGs correlated with Escherichia coli isolated from Simmental calves. Additionally, metabolomics analysis showed that differentially expressed metabolites in Simmental calves with diarrhea displayed a high correlation with the aforementioned ARGs and VFGs. Our findings provide an insight into the diversity and abundance of the gut microbiota in diarrheic calves caused by Escherichia coli and pave the way for further studies on the mechanisms of antibiotic resistance and virulence in the diarrheal calves from cattle hosts. Simmental is a well-recognized beef cattle breed worldwide. They also suffer significant economic losses due to diarrhea. In this study, fecal metagenomic analysis was applied to characterize the antibiotic resistance gene (ARG) and virulence factor gene (VFG) profiles of diarrheic Simmental calves. We identified key ARGs and VFGs correlated with Escherichia coli isolated from Simmental calves. Additionally, metabolomics analysis showed that differentially expressed metabolites in Simmental calves with diarrhea displayed a high correlation with the aforementioned ARGs and VFGs. Our findings provide an insight into the diversity and abundance of the gut microbiota in diarrheic calves caused by Escherichia coli and pave the way for further studies on the mechanisms of antibiotic resistance and virulence in the diarrheal calves from cattle hosts.

incidence is strongly associated with gut microbiota (1)(2)(3).Gut microbiota plays vital roles in shaping key aspects of postnatal life, including the development of the gut and immune system in newly born calves (4), influencing the host's energy balance (5), and preventing colonization of foreign pathogens (6).Therefore, understanding the composition and functional properties of gut microbiota is imperative to control the occurrence and adverse effects of neonatal calf diarrhea.
Gut microbiome in ruminants consists of different types of microorganisms that produce thousands of bioactive compounds and coexists harmoniously in the habitat of the host gut (7)(8)(9).Diversity of gut microbiota is an indicator of the effective microbial ecosystem and the host's health.Characterization of the diversity of gut microbiota is imperative for a better understanding of gut ecology and its physiological manifes tations.To date, 16S rRNA gene sequencing and shotgun metagenomic sequencing are the two swordsmen for detecting microorganisms and assessing microbial diver sity.The 16S rRNA gene sequencing has been widely used to identify microbiota and its diversity analysis in different livestock species, such as cattle (10), goats (11), sheep (12), and buffalo (13).In contrast, the shotgun metagenomic sequencing allows researchers to sequence all given genomic DNA from a complex sample, unlike 16S rRNA gene sequencing which only targets 16S rRNA genes.Hence, shotgun metage nomic sequencing can further explore the functional characteristics of microorganisms.Moreover, shotgun metagenomic sequencing has also been widely used for investigat ing microbial diversity in different farm animal species, including cattle (14), goats (15), sheep (16), and buffalo (17).Recently, the use of 16S rRNA gene sequencing of the fecal microbiota of calves revealed that microbial metabolism-related genes, such as starch and sucrose metabolism, sphingolipid metabolism, alanine aspartate, and glutamate metabolism, were significantly altered in diarrheal calves (18).However, the application of shotgun metagenomic sequencing for the identification of microbiota associated with calf diarrhea is still limited.
Numerous studies have characterized microbial diversity using metagenomic analysis; however, little is known about the identification of antibiotic resistance genes (ARGs) and virulence factor genes (VFGs) in the fecal microbiota of animals using the metagenomic approach.The ARGs and VFGs are well known to pose substantial health risks to both animals and humans.While ARGs contribute to resistance to antibiotics, VFGs facilitate infection and survival in adverse environments (19).It is documented that antimicrobial resistance levels vary greatly across countries, regions, and even herds or breeds (20)(21)(22)(23).Therefore, the issues of antibiotic resistance and virulence are attracting wide concern from the scientific community.Identification of the ARGs and VFGs is imperative to better understand pathogenesis and devise strategies to control disease and antibiotic resistance.Although many ARGs and VFGs have been identified for microbes found in the environment, these factors are poorly understood in the gut microbiome of ruminants.
Escherichia coli (E.coli) is an important pathogenic bacterium associated with diarrhea.According to the host or organ tropism, E. coli genomes were categorized into different pathotypes, such as extraintestinal pathogenic E. coli, intestinal pathogenic E. coli, and uropathogenic E. coli (24).Additionally, E. coli genomes were also grouped into the diarrhea-associated hemolytic E. coli and enterohemorrhagic E. coli (EHEC) based on the disease conditions, as well as Shiga toxin-producing E. coli (STEC) based on the virulence factors (25).For instance, several antibacterial resistance genes involved in canine diarrhea have been identified in different strains of E. coli, including EHEC, STEC, and necrotoxic E. coli (26).Eltai et al. (27) reported that diarrhea incidence caused by E. coli (DEC) strains had a high rate of antimicrobial resistance, and the atypical enteropathogenic E.coli and enteroaggregative E. coli were the primary etiological agents of diarrhea in children among DEC pathotypes.It is well documented that E. coli can be commensal or pathogenic because of its high genetic plasticity (28).The phylogroups of E. coli have been described (A, B1, B2, C, D, E, and F) on the basis of multilocus sequence typing (MLST) profiles (29,30).Moreover, E. coli species are currently divided into 183 O-serogroups (lipopolysaccharide) and 53 H-types (flagellar antigen) (31).In epidemiological studies, knowledge of the O-serogroups can help identify pathogenic lineages, classify E. coli for epidemiological studies, and identify outbreaks of diseases (32).Importantly, the relationship between antibiotic resistance and serogroups and virulence factors of E. coli was discovered in different disease conditions, such as children with septicemia (33), patients with urinary tract infection (34), and hospitalized patients (35).Nevertheless, similar studies on calf diarrhea are relatively limited.
In China, calf diarrhea is a common disease in cattle industry that causes huge morbidity and mortality losses.Reducing the occurrence of calf diarrhea is also a key health challenge faced by the Chinese cattle industry.Despite the importance of the gut microbiota in modulating host health, little is known regarding changes in the gut microbiota, ARGs, and VFGs and their potential associations with pathological status in different calf breeds affected by diarrhea.To devise effective strategies to control calf diarrhea, it is essential to better understand the gut microbiome, ARGs, and VFGs.In the present study, samples from calves suffering from diarrhea were taken from a Simmental cattle farm located in Henan province of China.In this regard, the aim of the present study was to characterize the changes in the gut microbiome and resistome of calves affected by diarrhea from Chinese Simmental beef cattle breed.Furthermore, we compared the differences in the relative abundances of ARGs and VFGs between the diarrheic calves due to E. coli and the healthy individuals.We believed that our potential findings would provide crucial insights into the gut microbiome of neonatal calves affected by diarrhea and pave the way for studies aimed at elucidating the antibiotic resistance and virulence mechanisms in diarrhea.

Experimental design
A total of 257 fecal samples were collected from Simmental calves with diarrhea reared at a cattle farm attached to the Institute of Animal Husbandry and Veterinary Medi cine, Henan Academy of Agricultural Sciences (Luoyang, China).Pathogen detection for all fecal samples was performed as previously described (36), and then the IDEXX Rota-Corona-K99 Ag Test was used to confirm the incidence of diarrhea.In briefly, fecal samples were fully resuspended in phosphate-buffered saline (1:10 [wt/vol]) and centrifuged at 10,000 × g for 10 min at 4°C.Next, viral RNA and DNA were extracted from 200 µL of the fecal suspension using the TaKaRa MiniBEST Viral RNA/DNA Extraction Kit Ver.5.0 (TaKaRa Bio Inc.) according to the manufacturer's instructions; another 200 µL of a suspension was used directly for the extraction of bacteria and protozoa nucleic acids with a QIAamp Fast DNA Stool Mini Kit (QIAGEN, Hilden, Germany).Subsequently, the multiplex PCR assay was used in the present study.PCR reaction mixture was composed of 10 µL of 2× Mix (TaKaRa Bio Inc.), 1.0 µL of each primer (10 µM), and 20 ng of template DNA in a total volume of 20 µL.All the PCR reactions were performed using a Thermal Cycler Dice Model TP600 (TaKaRa Bio Inc.).The detected pathogens consisted of the group A rotavirus, group B rotavirus, group C rotavirus, bovine coronavirus, bovine torovirus, bovine norovirus, bovine enteric Nebraska-like calicivirus, bovine nebovirus, bovine viral diarrhea virus, Clostridium perfringens, Salmonella enterica, Salmonella enterica Typhimurium, Escherichia coli, and Eimeria zuernii.Finally, 30 fecal samples, including 15 diarrheal calves with only positive K99 and 15 healthy individuals, were selected for shotgun metagenomic sequencing and non-targeted metabolomic sequencing.While not all calves were kept in the same space, a physical space was created for them to eat and interact with other calves during the experimental period (from calf birth to 18 days old).No routine treatments were performed except navel disinfection.
All selected calves were given the same colostrum at 8% of body weight within 2 h after birth, then antibiotic-free processed milk at 10% of BW was provided twice a day (06:00, 18:00) during their first 7 days of age.No waste milk (milk from cows treated with antibiotics or from cows with clinical mastitis) was added to the good milk fed to all calves.All feeding procedures were carried out by a professional operator using nipple buckets.Calves were gradually removed from bottle feeding and encouraged to drink milk from the buckets.Water was available ad libitum during the study period.All calves were provided ad libitum access to a calf starter (20% CP, Purina) from the seventh day of life until the end of the experiment.
The status and severity of diarrhea in the neonatal calves were assessed through a previously described procedure (37).Briefly, fecal scoring was recorded according to a standard scoring procedure (0 = normal feces; 1 = semi-formed feces; 2 = loose feces; 3 = watery feces).Individual calves were sampled before receiving antibiotics.Fecal samples for each individual were collected by the rectal palpation method while wearing clean disposable latex gloves and then stored at −80°C until DNA extraction.Details regarding sample collection from calves are summarized in Table S1

DNA extraction and shotgun metagenomic sequencing
Total DNA was isolated from 200 mg of fecal samples using the QIAamp PowerFecal DNA Kit (Qiagen, Maryland, USA) following the manufacturer's instructions.DNA concentra tion and purity were measured using a Qubit spectrophotometer (Nanodrop Technolo gies Inc., DE, USA) and 1.5% agarose gel electrophoresis, respectively.Paired-end library construction for each sample was performed by using the TruSeqTM DNA Sample Prep Kit (Illumina Inc., CA, USA) according to the manufacturer's instructions.Overall, genomic DNA samples were first fragmented to an average size of 350 bp using Covaris M220 (Gene Company Limited, Guangzhou, China).Then, the obtained fragments were end-repaired, A-tailed, and further ligated with adapters, as after that those ligated fragments were amplified by PCR, size selected, and purified.Finally, paired-end sequencing was carried out using Illumina NovaSeq System platform (Illumina Inc., CA, USA).

Metagenomics analysis
Raw sequence data for each sample were converted and filtered into clean data by the readfq toolkit ver8 (https://github.com/cjfields/readfq).The filtration criteria applied were as follows: (i) remove reads with low-quality bases that exceed 40 bp; (ii) exclude reads with N bases reaching 10 bp; and (iii) remove reads whose overlaps with adapters exceed 15 bp.In addition, clean data were mapped into the bovine reference genome (ARS-UCD1.2) using the Bowtie ver2.2.4 software (38), and any hits associated with clean data and their mated reads were also removed.
The clean data were used for assembling analysis using the MEGAHIT ver1.0.4 software (39) with default settings.The scaftigs with a length ≥500 bp were used to predict the genes using MetaGeneMark ver2.10 software (40).The CD-HIT ver4.5.8 software (41) with the parameters (-c 0.95 -aS 0.9) was utilized for clustering genes from each sample.Genes with reads >2 in each sample were retained and defined as the final non-redundant gene catalog in this study.Quantification of unique genes in each sample was performed using Salmon ver1.9.0 software (42).The relative abundances of genes was presented with transcripts per kilobase of exon model per million mapped reads (TPM) based on the following equation: where N i is the number of reads mapped to the ith gene, L i is the total length of an exon in the jth gene.Non-redundant gene sequences were taxonomically assigned using DIAMOND ver.0.9.9 software (43) against the NCBI NR database ver.2018.01.02.Taxonomic profiles were determined at the phylum, class, order, family, genus, and species levels, and the relative abundances at those taxonomic ranks were analyzed using MetaPhlAn ver.3.0 software (44).Taxonomic and phylogenetic trees of the gut microbiome from the two studied groups were generated by GraPhlAn ver1.1.3software (45).The amino acid sequences of non-redundant genes were functionally annotated using Diamond ver.0.9.9 software (43) against the Kyoto Encyclopedia of Genes and Genomes (KEGG; version 2018-01-01), carbohydrate-active enzymes (CAZy; version 2018-01-01), and EggNOG (v4.5) databases.

Non-targeted metabolomic analysis of fecal samples
The fecal samples from 30 calves were collected.The sample processing and sequencing were performed as described by our previous study (46).Here, a total of 14 diarrhea calves and 14 healthy individuals met the requirements for non-targeted metabolomic sequencing and were used for further analysis.The Compound Discoverer 3.1 (Thermo Fisher Scientific, Waltham, MA, USA) was used to process the raw data.The normalized data were used to predict the molecular formula based on additive ions, molecular ion peaks, and fragment ions.The qualitative and relative quantitative results of peaks were performed by the mzCloud (https://www.mzcloud.org/),mzVault, and MassList databases.Additionally, three databases including KEGG, Human Metabolome Database, and LIPID MAPS were used for metabolite annotation.Identification of differentially expressed metabolites (DEMs) was performed using univariate data analysis (t-test) with the following thresholds: projection (VIP) > 1, P-value < 0.05, and fold change > 2. KEGG analysis for DEMs was performed using the MetaboAnalyst (V5.0) (47).

Identification of antibiotic resistance genes and virulence factors in metage nomes
Relative abundances of ARGs and VFGs in the gut microbiome were quantified using abricate ver.1.0.0 software (48) against the Comprehensive Antibiotic Resistance Database (CRAD; ver.2.0.1) and virulence factor database (VFDB).To determine the relative abundances of both ARGs and VFGs at the species level, we summed up the TPMs of those ARGs or VFGs under each specific taxonomic group.

Co-occurrence network analysis for ARGs, VFGs, microbial taxa, and metabo lites
To investigate the co-occurrence pattern of microbial communities, resistomes, and metabolites, a correlation matrix was constructed by calculating the pairwise Pearson correlation with the "psych" R package.Both correlation coefficient of >0.8 and adjusted P-value of <0.05 were defined as the threshold levels.Finally, Cytoscape ver.3.10.1 software (49) was applied to convert the resulting correlation matrix into an associated network.

Phylotyping and serotyping for Escherichia coli
To obtain the sequence of E. coli genomes for each individual, the sequence of E. coli was first extracted from the NT database and used as the reference.All the assembled sequences for each sample were mapped into the reference genome using mummer ver.3.23 (50) with the following parameters: nucmer --maxmatch --nosimplify -g 5000 c 200 L 500 show-coords -rlo -L 500.Only the contigs for each sample were defined as the sequence of E. coli genomes with the threshold of similarity > 0.95 and coverage > 0.90.Moreover, we downloaded a total of 10,061 E. coli genomes from NCBI as of 24 November 2023, with the criteria that each assembly has (i) >200 scaffolds/contigs, (ii) GC percentage > 0.47, (iii) genome size > 2 Mb, and (iv) coding DNA sequence > 2,000 bp.The detailed information on the E. coli genomes is listed in Table S2 at https://doi.org/10.7910/DVN/UKE1CC.Sequence types (ST) of E. coli genomes were identified using in silico MLST ver.2.23.0 (51) software based on PubMLST database (https://pubmlst.org).Phylogroup typing of all E. coli genomes was performed using EzClermont ver.0.6.3 (52) software.The in silico prediction of E. coli serotype was identified using ECTyper ver.1.0.0 (53) software.Additionally, the abricate ver.1.0.0 (54) software was applied to predict the ARGs and VFGs for E. coli genomes against the CRAD (identity = 0.98 and coverage = 1.0) and VFDB (identity > 0.98 and coverage = 1.0) database, respectively.

Statistical analysis
One-way analysis of variance was performed using GraphPad Prism ver.10.0.3 (P < 0.05).Differences in relative abundances of taxa, functional genes, metabolic pathways, ARGs, and VFGs between different comparisons were identified using linear discriminant analysis effect size (LEfSe) or Wilcoxon's rank-sum test.LEfSe analysis was performed by using the LEfSe ver.1.1.2software.Significance testing of the principal coordinates analysis (PCoA) was performed based on the Bray-Curtis distances using the "vegan" R package (55).The analysis of similarity (ANOSIM) was applied to determine the difference in abundances of gut microbiome or ARGs between DSM and HSM groups.ANOSIM analysis was performed by using the "vegan" R package.The Benjamini and Hochberg false discovery rate was used to correct the P-value.Simpson, Shannon, and Chao indexes of the gut microbiota, ARGs and VFGs between comparison groups were performed to evaluate their diversity using the "vegan" R package.All plots were mainly visualized by ggplot2 packages (56).

Demographic parameters of the study calves
In this study, 30 Simmental calves were used.Details regarding their characteristics are presented in Table S3 at https://doi.org/10.7910/DVN/UKE1CC.Age, sex, and temperature had no significant differences (P > 0.05, Wilcoxon rank-sum test) between diarrheic and healthy calves.Only the fecal score value showed significant (P < 0.05; Wilcoxon rank-sum test) difference between diarrheal and healthy calves.

Gut microbiome community in diarrheic and healthy calves
Here, a total of 15 diarrheal and 14 healthy fecal samples met the library requirements for shotgun metagenomic sequencing.The results showed that a total of 2.16 Gb of raw reads were generated, of which 0.94 and 1.21 Gb belonged to DSM and HSM groups, respectively (see Table S4 at https://doi.org/10.7910/DVN/UKE1CC).After quality filtering and removing host genome sequences, 2.15 Gb of clean data was retained.Subsequently, a total of 2,699,417 contigs at the scaffold level were generated using the de novo assembly (see Table S5 at https://doi.org/10.7910/DVN/UKE1CC),with a mean N50 length of 7335.89bp.Furthermore, a total of 2,816,110 non-redundant genes were predicted, with a mean protein length of 234.35 and 211.72 bp in the DSM and HSM groups, respectively (see Table S6 at https://doi.org/10.7910/DVN/UKE1CC).
The ANOSIM analysis (R = 0.3039, P = 0.001) revealed significant differences in the relative abundance of the gut microbiota at the genus level between the two calf groups under study (Fig. 1A).The diarrheic calves showed a lower diversity (Shannon index, Simpson's index, and Chao) of gut microbiota at the species level compared to the healthy calves (Fig. 1B, C, and D).

Functional comparison of gut microbiome between diarrheic and healthy calves
We initially predicted the gene functions of the gut microbiome based on the KEGG database.The most abundant genes in the microbiomes of the DSM group were associated with ABC transporters, two-component regulatory system, microbial metabolism in diverse environment, phosphotransferase system (PTS), and biofilm formation-Escherichia coli (Fig. 2A).In contrast, aminoacyl-tRNA biosynthesis was the most abundant pathway in the gut microbiomes of the HSM group (Fig. 2A).The majority of genes enriched in membrane transport were abundant (P < 0.05) in the gut micro biomes of diarrheic calves compared to healthy individuals (Fig. 2B), followed by the (F) LEfSe analysis (P < 0.05 and LDA > 2.5) shows the microbiota compositions at the genus level between groups.(G) LEfSe analysis (P < 0.05 and LDA > 2.5) shows the microbiota compositions at the species level between groups.
signal transduction, cellular community prokaryotes, and xenobiotics biodegradation and metabolism.
Furthermore, we explored CAZymes activity in the gut microbiome based on CAZy database.These CAZymes were clustered into six classes: auxiliary activities (AA), carbohydrate-binding modules, carbohydrate esterases, glycoside hydrolases, glycosyl transferases, and polysaccharide lyases.Total genes of CAZymes showed significant differences (P < 0.05) in abundance between DSM and HSM groups (Fig. 2C).Notably, we found that these genes encoding AAs and cohesin exhibited a significant (P < 0.05) difference between DSM and HSM groups (Fig. 2D).Among them, the AA0, AA4, AA6, AA7, and AA17 in AA family were significantly more abundant in the DSM gut microbiome (Fig. 2E).

Antibiotic resistance genes of gut microbiome in diarrheic calves
CARD database was used to identify the ARGs of gut microbiome in diarrheic calves, generating a total of 175 ARG subtypes.NMDS2 analysis based on the ARGs profile showed that the abundance difference in ARG between DSM and HSM groups was evident (Fig. 3A).ANOSIM analysis also showed a clear separation between the two groups based on the ARGs abundance (R = 0.399, P = 0.001; Fig. 3B).Both Shannon and Simpson's indexes demonstrated that the diarrheic calves had higher (P < 0.05) diversity of ARGs than that of the healthy calves (Fig. 3C and D).
To uncover the abundance differences in ARGs between DSM and HSM groups, the LEfSe analysis revealed that a total of 20 ARG subtypes were significantly more abundant in the guts of diarrheic calves compared to the control group (Fig. 3E).Among them, a total of 13 ARGs were predicted to be associated with the resistance mechanism of antibiotic efflux based on the CARD database, followed by the 4 ARGs for antibiotic inactivation, as depicted in Table S7 at https://doi.org/10.7910/DVN/UKE1CC.

Virulence factor genes of gut microbiome in diarrheic calves
To identify the VFGs of the gut microbiome in diarrheic calves, the microbial genes were analyzed against the VFDB database.A total of 271 VFGs were detected, of which 227 VFGs were common between the DSM and HSM groups (Fig. 4A).PCoA analysis based on the VFGs abundance revealed that the PCoA1 and PCoA2 explained 35.52% and 17.47% of the variability, respectively, demonstrating that the virulence difference between DSM and HSM groups was almost completely separated (Fig. 4B).No significant difference in the diversity of VFGs between DSM and HSM calves was detected by the Shannon and Simpson's index (Fig. 4C and D, respectively).However, significant (P < 0.05) difference in the abundance of VFGs between the two studied groups was observed (Fig. 4E).
To investigate the abundance differences in VFGs between the diarrheic and healthy calves, LEfSe analysis showed that a total of 27 VFGs were significantly different in abundance between DSM and HSM groups (Fig. 4F).Differentially expressed VFGs had differences (P < 0.0001) in abundance between the diarrheic and healthy calves (Fig. 4G).Notably, eight VFGs, including acm, gspG, fepG, fimB, fimE, csgB, csgF, and yagZ/ecpA, were highly expressed abundance in the DSM compared to HSM group (Fig. 4H), indicating that these VFGs might play vital role in the occurrence of diarrhea.

Key metabolites in diarrheic calves caused by E. coli associated with ARGs and VFGs
Evidence revealed that the metabolic changes were associated with the incidence of diarrhea (57).We, therefore, used fecal samples for metabonomic analysis to determine the metabolic changes in diarrheic calves.Here, a total of 11,092 and 7,421 metabolites were identified for the positive and negative ion models, respectively.Furthermore, a total of 15,713 unique metabolites were obtained, and the detail information on these metabolites was listed in Table S8 at https://doi.org/10.7910/DVN/UKE1CC.PCA analysis of these metabolites showed that the studied samples had good quality and system stability (Fig. 6A).Moreover, a total of 2,405 DEMs were identified using the significance analysis with VIP > 1, fold change > 2, and P-values < 0.05, of which 575 metabolites were upregulated and 1,829 downregulated (Fig. 6B).KEGG analysis of DEMs showed that steroid hormone biosynthesis was the most significant pathway for enrichment (Fig. 6C).Moreover, correlation analysis revealed that 28 DEMs in the diarrhea calves caused by E. coli had strong significant correlation coefficients of >0.8 with 16 ARGs and 7 VFGs (Fig. 6D).The 28 DEMs were sodium cholate, Ala-Glu-Lys, Glu-Glu-Lys, Leu-Leu-Glu,  6E).The relative abundance of D-ornithine was found to be higher in the diarrheic calves compared to the healthy ones (Fig. 6F).Additionally, D-ornithine had a strong correlation with two ARGs and one VFGs.The relative abundance analysis showed that abundances of these ARGs (Fig. 6G) and VFGs (Fig. 6H) both were higher in the diarrheic calves compared to the healthy ones.

DISCUSSION
In the present study, the microbiomes and resistomes of purebred Simmental calves suffering from diarrhea were explored based on shotgun metagenomic sequencing.Importantly, we found that several bacterial species varied between the diarrheic and healthy calves.The Pseudomonadota (38.40%),Bacillota (31.57%), and Bacteroidota (22.48%) were the three most dominant phyla accounting for a major proportion of the gut microbiota in the diarrheic calves, which is quite similar to the relative abundances already reported in diarrheic animals (58)(59)(60).Significant differences in the abundance of nine genera, including Escherichia, Shigella, Lactobacillus, Salmonella, Enterocloster, Klebsiella, Cronobacter, Enterobacter, and Agathobaculum, were observed in the guts of diarrheic calves.Considering the species, Escherichia coli and Shigella sonnei were more abundant in the DSM than HSM groups.Previous studies have demonstrated associations between Escherichia (61), Shigella (62), Lactobacillus (63), Salmonella (64), Enterobacter (65), and Klebsiella (66) and diarrhea in animals.In addition, several distinct bacterial taxa were identified at the species level as being different between the healthy calves and those calves suffering from diarrheal disease in our study.It is noteworthy that the Escherichia coli and Shigella sonnei, as pathogens, were significantly (P < 0.05) different in the relative abundance of diarrheic calves compared to healthy calves.It   is well known that both E. coli and Shigella sonnei are the pathogenic bacteria and associated with increased incidence of diarrhea in farm animals (67)(68)(69) and humans (70).These findings showed that E. coli played a vital role in the occurrence of diarrhea.Gut microbiota plays an important role in shaping the health status of animals through alteration of the host's metabolic functions (71).In the present study, we observed differences at the functional gene levels between diarrheic and healthy calves.These microbial genes were mainly involved in the ABC transporters, two-com ponent system, microbial metabolism in diverse environment, PTS, and biofilm formation Escherichia coli.Evidence showed that ABC transporters were increasingly recognized for their ability to modulate the absorption, distribution, metabolism, secretion, and toxicity of xenobiotics (72).Interestingly, our data showed that the abundance of microbial genes in the xenobiotics biodegradation and metabolism pathway was increased in the diarrheal calves compared to healthy individuals.Additionally, the two-compo nent systems were reported to be the predominant way bacteria adapt to changing environments and modulate their fitness in various niches (73)(74)(75).These results indicate that the microbial genes may play a vital role in the adaptability to xenobiotic toxicity.Furthermore, we observed significant differences in the relative abundances of CAZymes between diarrheic and healthy calves, particularly with regard to some genes encoding AAs.It has been reported that AAs are involved in the catabolism of carbohydrates (76).This could imply that the bacterial genes encoding AAs may impact the gut microbiome metabolism in individuals affected by diarrhea.Cheng et al. (77) also highlighted specific carbohydrates that play a crucial role in human diseases by regulating the imbalanced intestinal microbiota.
Widespread antibiotic use causes antibiotic selection pressure, resulting in potential emergence of antibiotic-resistant bacteria and reduction of microbiota diversity, thereby accelerating the rise of AMR levels (78)(79)(80).It is well known that coexisting microbes can produce antibiotics or bacteriocins in response to the competition for nutrients in the intestinal digesta (81)(82)(83).Here, we observed that the diarrheal calves had higher abundance and diversity of ARGs compared to healthy calves.This might be partially explained by the differences in the ability of bacterial species to adapt to changes in nutritional habits, diets, and environmental conditions on farms as well as selection pressures exerted by antibiotics (84).Similar findings have been reported by other studies (85,86).Importantly, significant differences in the richness and the abundance of ARGs between diarrheal and healthy calves were observed.These ARGs contained the floR, APH (6)-Id, TEM, mdtL, emrD, cpxA, CRP, ampC, marA, evgA, H-NS, mphB, gadX, bacA, emrR, gadW, acrS, ramA, dfrA27, and arr-3.According to resistance mechanism analysis of CARD database, a total of 13 ARGs, including the floR, mdtL, emrD, cpxA, CRP, marA, evgA, H-NS, gadX, emrR, gadW, acrS, and ramA, were the antibiotic efflux, followed by 4 ARGs (TEM, ampC, mphB, and arr-3) and APH (6)-Id for antibiotic inactivation, bacA for antibiotic target alteration, and dfrA27 for antibiotic target replacement.These results indicate that the resistance mechanisms of these ARGs in the gut microbiome are mainly involved in the antibiotic efflux.
Co-occurrence pattern showed that a total of 48 bacterial genera displayed significant (P < 0.0001) correlation with 20 ARGs (correlation coefficients >0.8).Importantly, these bacteria have been demonstrated to play essential roles in the drug resistance through the resistance mechanism of antibiotic efflux (87)(88)(89).Notably, our analysis showed that Escherichia coli had highly significant (P < 0.0001) correlation coefficients of >0.8 with 16 ARGs, which were mainly involved in the antibiotic efflux.These findings are in agreement with earlier studies demonstrating that E. coli harbored the largest num ber of ARGs, showing multidrug resistance in a wide variety of environments (90)(91)(92).Moreover, for these E. coli-mediated ARGs, they had higher abundance in the diarrheal calves compared to healthy individuals.These results point out that the development of antimicrobial resistance in diarrhea calves may be related to the use of maternal antibiotics, which are passed through the mother's colostrum.
Virulence factors of a bacterium are molecules secreted, membrane-associated or cytosolic agents that help the bacteria to colonize their hosts (93).Bacterial patho gens use their secretion system to directly incapacitate target cells and enhance their virulence (94).In the present study, we found that 18 VFGs, including the iroC, iroD, iroE, gspG, fepG, fimE, sfaY, sfaE, sfaF, sfaG, iroN, iroB, csgF, csgB, papF, acm, yagZ/ecpA, and fimB, were significantly enriched in the diarrheic calves (P < 0.05 and LDA >2.5) compared to the healthy individuals.In particular, co-occurrence patterns showed that E. coli had strong positive correlation coefficients of >0.8 with the aforementioned VFGs.Previous studies shown that both fimB and fimE were responsible for controlling the expression of type I fimbriae in E. coli (95)(96)(97).Swasthi et al. (98) reported that curli were functional amyloids present on the outer membrane of E. coli, while csgB and csgF were required for curli assembly (98).Four genes (sfaY, sfaE, sfaF, and sfaG) were known to play a vital role in the adhesion of S-fimbriae in the pathogenesis of E. coli meningitis (99).Evidence showed that five iro genes (iroC, iroD, iroE, iroN, and iroB) were required for the virulence of Avian pathogenic E. coli (100).Additionally, it is well documented that fepG was essential for ferrienterobactin transport in E. coli (101), while papF played various functional roles in pilus biogenesis of E. coli (102).These results indicate that these VFGs may be involved in various functions of E. coli pathogenicity.
E. coli isolates are known to be a vast diversity in terms of clonal lineages due to the transition between commensalism, mutualism, and pathogenicity, depending on the selection pressures and niche-specific constraints or adaptations experienced by the host genomes (114).In this study, we demonstrated that the currently investiga ted E. coli genomes constructed from 15 samples consisted of 8 major phylogroups, which is similar to previous studies (115)(116)(117).For the diarrheal Simmental calves, we observed that three individuals had E. coli genomes from the phylogroup B1, while seven individuals had E. coli genomes from unknown phylogroups.Coura et al. (118) found that the phylogroups B1 and E were associated with E. coli isolated from cattle, while phylogroups B2 and D were associated with E. coli collected from water buffalo (118).Recently, Özgen et al. (119) reported that the phylogenetic group B1 has been observed in the case of bovine abortions (119).Notably, Aguirre-Sánchez et al. (120) mentioned that the phylogroup B1 for E. coli isolated from cattle located in Mexico demonstrated a significant relationship with food sources (120).Moreover, we found that 10 ARGs (emrR, CRP, ampC1, mphB, acrS, cpxA, evgA, marA, H-NS, and APH (6)-Id) and 3 VFGs (fimB, fimE, and yagZ/ecpA) showed strong correlations with the E. coli from phylogroup B1.Therefore, it can be inferred that the ARGs are mainly related to antibiotic efflux, while the VFGs are correlated to the pilus function of E. coli.These results imply that these ARGs and VFGs may play a vital role in the occurrence of diarrhea incidences caused by E. coli in the Simmental calves.

Conclusions
In summary, our study concluded that the relative abundances of fecal microbiota in the diarrhea calves displayed drastic changes between diarrhea and healthy calves.Diarrheic Simmental calves exhibited significant differences in abundance of antibiotic resistance and virulence factor genes.Importantly, our findings revealed that E. coli had a strong significant positive correlation with most of the ARGs and VFGs.The ARGs in diarrhea calves were mainly involved in antibiotic efflux, while these VFGs were associated with VFGs with pili function, curli assembly, and ferrienterobactin transport of E. coli.Metabolomics analysis showed that differentially expressed metabolites in Simmental calves with diarrhea displayed a high correlation with the aforementioned ARGs and VFGs.Our study provides crucial insights into the gut microbiome and ARGs of neonatal calves affected with diarrhea.Z.S. and T.D. conceived, designed, and performed the data analysis.Z.S. and T.D. drafted the manuscript.F.-U.H. and H.E.R. participated in the writing of the manuscript.Z.S., Y.L., Y.W., X.Y., X.M., and Z.X. were involved in the animal experiment and sample analysis.H.E.R. deeply reviewed and edited the final manuscript.Z.S. and W.W. supervised and provided continuous guidance for the experiment.All authors contributed to the article and approved the final manuscript.

FIG 1
FIG 1 Gut microbiome community in diarrheic and healthy Simmental calves.(A) ANOSIM analysis displays the difference in the structure of bacterial community between diarrheic and healthy groups.(B) Comparisons of Shannon index of gut microbiota at the species level between diarrheic and healthy groups.(C) Comparisons of Simpson's index of gut microbiota at the species level between diarrheic and healthy groups.(D) Comparisons of Chao index of gut microbiota at the species level between diarrheic and healthy groups.(E) Bar chart representing the microbiota compositions at the genus level between groups.

FIG 2
FIG 2 Gut microbial function analysis between diarrheic and healthy Simmental calves.(A) Difference of KEGG pathway abundance by LEfSe (P < 0.05 and LDA > 2.5).(B) Comparison of main KEGG function.The red color means significantly more abundant in the gut microbiome of diarrheic Simmental calves, and green color indicates significantly more abundant in the gut microbiome of healthy Simmental calves.(C) Total abundance comparison of CAZymes between groups.(D) Abundance comparison of CAZymes between groups.(E) Abundance comparison of the auxiliary activities (AA), polysaccharide lyases (PL), and glycoside hydrolase (GH).*, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.

FIG 3
FIG 3 Antimicrobial resistance gene analysis of gut microbiome between diarrheic and healthy Simmental calves.(A) DMDS analysis displays the resistance difference of bacterial community between diarrhea and healthy groups.(B) ANOSIM analysis of ARGs displays the resistance difference of bacterial community between diarrhea and healthy groups.(C) Comparisons of Shannon index of ARGs between diarrheic and healthy groups.(D) Comparisons of Simpson's index of ARGs between diarrheic and healthy groups.(E) Difference in the antimicrobial resistance of gut microbiome between groups based on LEfSe analysis (P < 0.05 and LDA > 2.5).**, P-value < 0.01; ***, P-value < 0.001.

FIG 4
FIG 4 Virulence factors gene analysis of gut microbiome between diarrheic and healthy Simmental calves.(A) Venn analysis of VFGs between two groups.(B) PCoA analysis for diarrheic and healthy Simmental calves based on VFGs.(C) Comparisons of Shannon index of VFGs between diarrheic and healthy groups.(D) Comparisons of Simpson's index of VFGs between diarrheic and healthy groups.(E) Comparisons of the abundances of all VFGs between two groups.(F) Difference of VFGs abundance between diarrheic and healthy groups by LEfSe (P < 0.05 and LDA > 2.5) analysis.(G) Comparisons of the abundances of 18 VFGs between two groups.(H) Heatmap plot displays the abundance distribution of 18 VFGs across all samples.NS indcates no significant difference.

FIG 5
FIG 5 Correlation analysis between ARGs, VFGs, and microbial taxa.(A) Correlation between ARGs and microbial taxa.(B) Correlation between VFGs and microbial taxa.(C) Correlation between VFGs and ARGs.(D) Co-occurrence patterns among the ARGs, VFGs, and Escherichia coli based on their correlation coefficient of >0.8.

FIG 6
FIG 6 Identification of differentially expressed metabolites related to diarrhea and its correlation with ARGs, VFGs, and microbial taxa.(A) PCA analysis for diarrheic and healthy Simmental calves based on metabolites.(B) Volcano plot displayed the differentially expressed metabolites between diarrheic and healthy Simmental calves.(C) KEGG enrichment analysis of differentially expressed metabolites.(D) Co-occurrence patterns among the DEMs, ARGs, VFGs, and E. coli based on their correlation coefficient of >0.8 and P-value of <0.05.(E) KEGG enrichment analysis of DEMs associated with the ARGs, VFGs, and/or E. coli.(F) Comparisons of the abundances of D-ornithine between two groups.(G) Comparisons of the abundances of mphB and ampC between two groups.(H) Comparisons of the abundances of gspG between two groups.*, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.

FIG 7
FIG 7 Phylogroup analysis of Escherichia coli genome based on their MLST profiles.Note: blue indicates the phylogroup B1, orange indicates the phylogroup A, green indicates phylogroup B2, purple indicates phylogroup D, cyan sky indicates phylogroup E, brown indicates phylogroup C, bright purple indicates phylogroup G, pink indicates phylogroup F, and dark sapphire indicates phylogroup unknown.

TABLE 1
Phylotype, serotype, ARGs, and VFGs of E. coli genomes in diarrheic Simmental calves a,b (Continued) Bold fonts indicate ARGs or VFGs that are highly associated with calf diarrhea.c NA indicates the serotype unknown b