Skip to main content
Download PDF

Metagenomics analysis revealed the distinctive ruminal microbiome and resistive profiles in dairy buffaloes

Animal Microbiome volume 3, Article number: 44 (2021) Cite this article

Abstract

Background

Antimicrobial resistance poses super challenges in both human health and livestock production. Rumen microbiota is a large reservoir of antibiotic resistance genes (ARGs), which show significant varations in different host species and lifestyles. To compare the microbiome and resistome between dairy cows and dairy buffaloes, the microbial composition, functions and harbored ARGs of rumen microbiota were explored between 16 dairy cows (3.93 ± 1.34 years old) and 15 dairy buffaloes (4.80 ± 3.49 years old) using metagenomics.

Results

Dairy buffaloes showed significantly different bacterial species (LDA > 3.5 & P < 0.01), enriched KEGG pathways and CAZymes encoded genes (FDR < 0.01 & Fold Change > 2) in the rumen compared with dairy cows. Distinct resistive profiles were identified between dairy cows and dairy buffaloes. Among the total 505 ARGs discovered in the resistome of dairy cows and dairy buffaloes, 18 ARGs conferring resistance to 16 antibiotic classes were uniquely detected in dairy buffaloes. Gene tcmA (resistance to tetracenomycin C) presented high prevalence and age effect in dairy buffaloes, and was also highly positively correlated with 93 co-expressed ARGs in the rumen (R = 0.98 & P = 5E-11). In addition, 44 bacterial species under Lactobacillus genus were found to be associated with tcmA (R > 0.95 & P < 0.001). L. amylovorus and L. acidophilus showed greatest potential of harboring tcmA based on co-occurrence analysis and tcmA-containing contigs taxonomic alignment.

Conclusions

The current study revealed distinctive microbiome and unique ARGs in dairy buffaloes compared to dairy cattle. Our results provide novel understanding on the microbiome and resistome of dairy buffaloes, the unique ARGs and associated bacteria will help develop strategies to prevent the transmission of ARGs.

Background

Antimicrobial resistance (AMR) is one of the biggest threats to human health and causes several to ten million death annually in the world [1]. The wide spread of antimicrobial resistance genes (ARGs), especially for those that confer resistance to clinical relevance or “last defense line” antibiotics, requires urgent tasks to reveal the harbor sources and transmission routes of diverse ARGs [2]. The microbiome of livestock gastrointestinal tracts serves as important reservoirs for ARGs [3,4,5,6], which is largely attributed to the long-term use of antibiotic in veterinary medicine or sub-therapeutic doses of antibiotic supplied in the feed acting as growth promoter [7]. The microbial ARGs are mainly structured by the bacterial phylogeny [8]. The ruminal bacterial composition and diversity are much higher than fecal microbiome in dairy cows [9], suggesting that rumen microbiome is a larger reservoir of ARGs.

The ARGs harbored in direct-fed microbials or silage inoculants that isolated from rumen microbiome (e.g. lactic acid bacteria or fibrolytic bacteria) [10] could be directly or indirectly transmitted and accumulated in human and other food animals. In addition, the bacterial ARGs in the rumen will easily spread to the environment, such as soil and water, through rumination and saliva, which will lead to drug-selective pressure on environmental and human life related bacterial communities [11]. The high-throughput metagenomics technology has been gradually applied into functional annotation and ARGs identification of large amounts of uncultured bacteria in gut microbiome of different species [12]. Nevertheless, the ruminal resistome information is still limited [13], especially for more diverse ruminant species.

As an important ruminant species, buffaloes are estimated to be more than 200 million worldwide and made great contribution to the nutrition and health of humankind [14]. Compared with dairy cows, dairy buffaloes have some unique physiological patterns, such as higher nutritious milk for cheese and cream production [15], better adaptability to low-quality and less-digestible roughage, hot and humid climate [16], advantages in resistance and susceptibility to common bovine disease [17]. Comparative analysis of rumen microbiome nowadays attracts increasing interest to understand the above distinct physiological patterns in dairy buffaloes [18].

Dairy buffaloes are generally raised in a low-density system with waterlogged conditions [19], which exhibit close interactions with the environment. Therefore, it is hypothesized that more diverse and distinctive microbiome and resistome existed in the rumen of dairy buffaloes compared with dairy cows. Here, we collected rumen samples from 16 dairy cows and 15 dairy buffaloes to perform compositional and functional comparative analysis using metagenomics. Since the average culling time for dairy buffaloes worldwide is generally longer than that of dairy cows (> 10 vs. 5 ~ 6 years old), dairy buffaloes with different ages (from 1 to 10 years old) were selected to reveal the comprehensive rumen microbial profiles. The current study provides comprehensive understanding on microbiome and resistome profiles of dairy buffaloes, and uncovers potential risks to clinical treatment and human health caused by spreading ruminal ARGs from dairy buffaloes.

Methods

Experimental design, animals and sample collection

All experimental procedures were approved by the Animal Care Committee of Zhejiang University (Hangzhou, China). A total of 16 dairy cows and 15 dairy buffaloes were used in this study. Sixteen Holstein dairy cows (3–5 year-old) were selected from a commercial farm in Hangzhou (China). Fifteen Murrah dairy buffaloes with different ages (Y: 12 month-old, n = 5; M: 3–5 year-old, n = 4; E: 6–8 year-old, n = 3; O: > 9-year-old, n = 3) were selected from buffalo farm of Buffalo Research Institute, Chinese Academy of Agricultural Sciences (Nanning, China). The animals were healthy and received no therapeutic or prophylactic antimicrobial treatment during the current lactation cycle or in the past year. All the selected animals were healthy and received no therapeutic or prophylactic antimicrobial treatment based on the veterinary records. The diet of dairy cows and dairy buffaloes were supplemented (Table S1). Fifteen mL of rumen fluid were collected from each animal before morning feeding using an oral stomach tube as described previously [20]. The first 50–100 mL of rumen fluids in each samping were discarded to remove the potential saliva contamination. The rumen content samples were immediately frozen in liquid nitrogen and stored at − 80 °C until further analysis.

DNA extraction, library construction, sequencing and data processing

Total DNA of rumen contents (3 mL) was extracted using the repeat bead-beating plus column method [21]. The concentration and quality of DNA were evaluated using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and 1% agarose gel. The Covaris M220 (Covaris lnc., Woburn, MA, USA) was used to generate 300 bp DNA fragments. The TrueSeq DNA PCR-Free Library Prep Kits (Illumina, San Diego, CA, USA) were applied to perform metagenome libraries construction according to the manufacturer’s protocol. Metagenomic sequencing (PE150) was performed on an Illumina HiSeq 3000 platform (Illumina, San Diego, CA, USA). Totally, 189.5 Gb clean data (from 194.5 Gb raw data) with an average of 61.11 million clean reads in each sample were generated.

The “N” records, 3′- and 5′-ends reads, short reads (< 50p) and low-quality bases (Phred score < 20) were filtered using fastp software [22]. Host genomic and adaptor contaminations were removed by mapping to the Bos taurus and Bubalus bubalis genomes using Burrows-Wheeler alignment (BWA) software [23]. After trimming, a set of high-quality and free of host genomic contaminants reads were obtained for further analysis.

The clean reads were de novo assembled for each sample using Megahit software based on succinct de Bruijn graph method [24]. The assembled contigs with the length greater than 300 bp were predicted to open reading frames using MetaGene software [25]. A total of 23,963,622 contigs with the length from 300 to 415,681 bp were obtained. Assembled contigs were pooled and then used to construct non-redundant (NR) gene catalog using CD-HIT with 90% identity and 90% coverage [26]. Gene abundances were calculated as gene reads/gene length using SOAPaligner software based on 95% identity [27]. And when comparing the taxa difference between two groups, the gene abundance was normalized as relative abundance: reads abundance of taxon A in sample 1/total reads abundances of all taxa in sample 1 × 100, to remove the bias of sequencing depth among different samples.

Taxonomic classification and functional annotation

Taxonomic classification was performed using DIAMOND software [28] to map the sequence of NR gene catalog against the NR database (ftp://ftp.ncbi.nlm.nih.gov/blast/db/). Taxonomic abundance profiles were generated at domain, kingdom, phylum, class, order, family, genus and species levels. The relative abundance of each taxon was calculated and normalized using the proportion of each taxa identified in the same sample. Only the taxa with a relative abundance > 0.1 in more than 50% of animals within each group were selected for downstream analysis. Stacked histogram was used to display the abundance difference of phylum and family levels bacteria.

The sequence of NR gene catalog was mapped to Kyoto Encyclopedia of Genes and Genomes database (KEGG, http://www.genome.jp/kegg/) using Diamond with BLASTP (http://blast.ncbi.nlm.nih.gov/Blast.cgi) type and an e value of 1E-5. Pathway annotation was then conducted using KOBAS 2.0 according the blast results [29]. Carbohydrate-active enzymes were identified based on HMMscan tools against the Carbohydrate-Active enZYmes Database (CAZy) with hmmer type and e-value ≤1E-5 [30]. Abundances of KEGG pathway and CAZymes were normalized into counts per million reads (CPM) for further comparison analysis with the cutoff of CPM > 5 in greater than 50% of animals in at least one group.

ARGs identification and resistome analysis

The sequence of NR gene catalog was mapped to the Comprehensive Antimicrobial Resistance Database (CARD) using Diamond software with BLASTP type and an e-value ≤1E-5 to identify the ARGs [28]. The corresponding AMR mechanisms and target antibiotics of ARGs were also provided in the CARD database [28]. The read counts of ARGs in each sample were normalized into CPM and only the ARGs with CPM value > 1 in more than 50% of animals in at least one group were used for down-stream analysis. Unique ARGs were defined as the genes absent in one group but with CPM values > 1 in all samples in the other group.

Supervised partial least-squares discrimination analysis (PLS-DA) was performed using ropls package [31] in R software to reveal discriminate patterns of ARGs between dairy cow and dairy buffalo groups. Pheatmap package in R was used to analyze the classification of ARGs and samples based on ARGs expression abundance. The co-expressed ARGs modules and their complex correlations with unique ARGs in dairy buffaloes were characterized using weighted gene co-expression network analysis (WGCNA) package in R [32]. To achieve a signed R2 value of scale-free topology model greater than 0.9 [33], the soft-thresholding power value was set as 9. The minimum count of genes assigned in each module was selected at 10. The cut-off of significantly correlated gene modules was defined as module-unique ARGs relationship |R| > 0.7 & P < 0.01. The association between membership and gene significance in the most correlated module was further explored. The correlation between the tcmA and bacterial species in the rumen of dairy buffaloes were calculated by Hmisc package in R. The co-occurence network was visualized by Cytoscape software [34].

Statistical analysis

The Bray-Curtis dissimilarity matrices based principal-coordinate analysis (PCoA) and P values calculation at genus and species level were performed using adae4 (version 1.7.13) [35] and vegan (version 2.5.4) [36] packages in R software. The significantly different features (SDFs) at each taxa level between dairy cows and dairy buffaloes were analyzed using least discriminant analysis (LDA) Effect Size (LEfSe) [37], in which both statistical significance and biological relevance were considered. The SDFs was defined at the thresholds of P < 0.01 & LDA > 3.5. The KEGG pathways and CAZymes between dairy cow and dairy buffalo groups were compared using the Wilcoxon rank-sum test and displayed in the volcano plot using R. The functional memberships with a false discovery rate (FDR) < 0.01 & |fold change (FC)| > 2 were considered as significantly different. The extremely higher CAZymes were defined as FDR < 0.01 & FC > 25. The difference of unique ARGs abundance among Y, M, E and O groups in dairy buffaloes were analyzed using ANOVA test with a P < 0.05 considered as statistical significance.

Results

Different rumen bacterial taxonomic compositions between dairy cows and dairy buffaloes

Bacteria was most predominant in all animals (Table S2) and showed significantly higher relative abundance in dairy buffaloes than in dairy cows (96.56% vs. 94.25%, P < 0.001). At phylum and family levels, dairy buffaloes had similar bacterial composition with dairy cows (Fig. 1A and B). Totally, 15 different bacterial phyla were identified (relative abundance > 0.1% in more than 50% animals), with Bacteroidetes, Firmicutes, Proteobacteria, Actinobacteria and Spirochaetes being the top 5 most abundant (accounting for > 89% of totally bacterial community) phyla in both dairy cows and dairy buffaloes (Fig. 1A). At family level, 45 different bacterial taxa were identified, with Prevotellaceae, Lachnospiraceae, Bacteroidaceae, Ruminococcaceae and Clostridiaceae being the top 5 most abundant (accounting for > 60% of totally bacterial community) families in both dairy cows and dairy buffaloes (Fig. 1B). At genus and species levels, the bacterial community displayed significantly different profiles between dairy cows and dairy buffaloes in PCoA plot (P = 0.001, Fig. 1C and D). Among the 44 SDFs under the cutoff of LDA > 3.5 and P < 0.01, 24 bacterial taxa were significantly higher, while 20 bacterial taxa were significantly lower in dairy buffaloes than dairy cows (Fig. 1E). The key species significantly enriched in the rumen of dairy buffaloes including Bacterium F082, Prevotella ruminicola, Faecalibacterium sp. CAG:74, Lactobacillus amylovorus, Prevotella sp. P6B4, Acetobacter pasteurianus, and Bacterium F083 (Fig. 1E).

Fig. 1
figure 1

Rumen bacterial compositional profiles of dairy cows and dairy buffaloes. A Relative abundance of major bacteria phyla (relative abundance > 0.1% in more than 50% animals) for all individuals. B Relative abundance of major bacteria families (relative abundance > 0.1% in more than 50% animals) for all individuals. C Principal coordinate analysis (PCoA) of bacteria genera based on Bray-Curtis distance. D PCoA of bacteria species based on Bray-Curtis distance. E Cladogram of significantly different taxa identified in the rumen microbiome data sets of dairy cows and dairy buffaloes based on the cut-off of LDA > 3.5 and P < 0.01. Clades significantly enriched in each cohort are highlighted by the colors shown in the legend. DB: dairy buffalo; SDF: significantly different feature; LDA: linear discriminant analysis

Full size image

Functional difference of rumen microbiome between dairy cows and dairy buffaloes

By aligning the NR gene catalog to the KEGG and CAZy database, metabolic pathways and carbohydrate-active enzymes were analyzed and compared between dairy cows and dairy buffaloes. Among the 21 significantly different third-level pathways between two groups (FDR < 0.01 & |FC| > 2), 17 pathways were enriched in the dairy buffaloes and 4 pathway were enriched in the dairy cows (Fig. 2A).

Fig. 2
figure 2

Functional difference of ruminal microbiome between dairy cows and dairy buffaloes. A The volcano plot of significantly different 3rd level KEGG pathways between dairy cows and dairy buffaloes. B The volcano plot of significantly different CAZymes between dairy cows and dairy buffaloes. C The distribution of significantly different CAZymes families between dairy cows and dairy buffaloes. D The relative abundance of extremely higher CAZymes encoded genes in the dairy buffaloes. FC: fold change (dairy buffaloes/dairy cows); GHs: glycoside hydrolases; PLs: polysaccharide lyases; CBMs: carbohydrate-binding modules; GTs: glycosyltransferase; AAs: auxiliary activities; DB: dairy buffalo; DC: dairy cow; CPM: counts per million

Full size image

For the CAZymes profiles, a total of 413 genes encoding family-level CAZymes were identified in the rumen microbiome of dairy cows and dairy buffaloes (Table S3). As shown in Fig. 2B, 59 of CAZymes encoded genes were up-regulated and 36 were down-regulated in dairy buffaloes than in dairy cows (FDR < 0.01 & |FC| > 2). Among the CAZymes encoded genes significantly higher expressed in the dairy buffaloes (Fig. 2C), 56% of them belonged to glycoside hydrolases (GHs), 17% to glycosyltransferase (GTs), 17% to carbohydrate-binding modules (CBMs), 5% to polysaccharide lyases (PLs), and 5% to auxiliary activities (AAs). The extremely higher CAZymes encoded genes in dairy buffaloes were GH71 (FDR = 6.70E-9, FC = 38.72), GH70 (FDR = 3.49E-4, FC = 54.86), GH68 (FDR = 0.0015, FC = 167.7), GT48 (FDR = 0.003, FC = 90.53), AA1 (FDR = 4.89E-05, FC = 26.77) (Fig. 2D).

Resistome profiles of the rumen microbiome and distinct ARGs patterns between dairy cows and dairy buffaloes

Dairy buffaloes showed distinct resistome profiles from dairy cows (Fig. 3). Totally, 505 ARGs were detected in rumen microbiome of dairy cows and dairy buffaloes; and 435 of them were expressed in at least one group (CPM > 1 in more than 50% animals). In the PLS-DA plot, the resistome profiles displayed clearly discrimination between dairy cows and dairy buffaloes, in which the PC1 explained 39.4% of total variations (Fig. 3A). The heatmap revealed that 16 dairy cows and 15 dairy buffaloes were assigned into two separated groups (Fig. 3B). More than 60% of ARGs showed higher expression level in dairy buffaloes than in dairy cows in the heatmap overview (Fig. 3B). The relative abundance of top 10 most abundant ARGs accounted for about 32.10 and 33.13% of total reads in dairy cows and dairy buffaloes, respectively (Table S4). Staphylococcus aureus rpoC conferring resistance to daptomycin was the most abundant ARG in both groups. In terms of AMR phenotypes of ruminal resistome in dairy buffaloes, 33.13% of the sequence reads corresponded to fluoroquinolone (11.87%), lipopeptide antibiotic (7.46%), and aminocoumarin (6.20%) AMR (Table S5).

Fig. 3
figure 3

Rumen resistome profiles in dairy cows and dairy buffaloes. A The principal component analysis of ARGs identified in the rumen of dairy cows and dairy buffaloes. B The heatmap of the relative abundance of all the ARGs in all the individuals

Full size image

Unique ruminal ARGs in dairy buffaloes and their resistance mechanisms

Eighteen ARGs were uniquely detected in the rumen microbiome of dairy buffaloes (CPM > 1 in all samples of dairy buffaloes and CPM = 0 in all samples of dairy cows) (Fig. 4). These ARGs included tcmA (CPM = 399.92), tet43 (CPM = 208.52), tet41 (CPM = 39.44), tet34 (CPM = 70.76), ARO:3003379 (CPM = 75.58), otrC (CPM = 38.84), otrB (CPM = 174.49), opcM (CPM = 232.42), ARO:3003392 (CPM = 216.70), mexV (CPM = 37.12), mexH (CPM = 39.07), mexC (CPM = 64.17), mdtD (CPM = 363.64), sul3 (CPM = 61.65), dfrC (CPM = 9.93), iri (CPM = 158.81), APH (3′)-IIc (CPM = 126.53), and adeA (CPM = 55.19) (Table S6). Unique ruminal ARGs showed large individual variance within dairy buffaloes (Fig. 4). The relative abundance of tcmA (P = 0.002), otrB (P = 0.021), sul3 (P = 0.043), iri (P = 0.021) presented significant difference among different age groups of dairy buffaloes (Fig. 4).

Fig. 4
figure 4

The unique ruminal ARGs detected in the dairy buffaloes. The relative abundance, target antibiotics and resistant mechanisms of 18 unique ARGs identified in each sub-group (Y/M/E/O). ARGs: antibiotic resistance genes

Full size image

Most of the unique ARGs (tcmA, tet43, tet41, ARO:3003379, otrC, otrB, opcM, mexV, mexH, mexC, mdtD, adeA) conferred resistance through the mechanism of “efflux pump conferring antibiotic resistance” (Fig. 4). The resistance mechanisms of tet34, ARO:3003392, sul3 and dfrC, iri and APH (3′)-IIc were “gene conferring antibiotic resistance via molecular bypass”, “ARGs variant or mutant”, “antibiotic target replacement protein”, and “antibiotic inactivation enzyme”, respectively (Fig. 4). These ARGs conferred resistance to tetracenomycin C, tetracycline, glycylcycline, aminoglycoside, rifamycin, diaminopyrimidine, sulfonamide, aminocoumarin, phenicol, fluoroquinolone, penam, macrolide, cephalosporin, acridine dye, isoniazid, and triclosan (Fig. 4).

The relationship between unique ARGs and bacterial species in the rumen of dairy buffaloes

To investigate the roles of unique ARGs in the ruminal resistome of dairy buffaloes, the relationship between co-expressed ARGs modules and unique ARGs was analyzed using WGCNA. Fourteen co-expressed gene modules that contained different number of ARGs (12–93) were identified (Fig. 5A). Totally, 16 significantly correlated relationships between 7 gene modules and 18 unique ARGs were observed (|R| > 0.7 & P < 0.01). Among them, MEturquoise module was positively associated with 6 ARGs, including tcmA (R = 0.98 & P = 5E-11), otrB (R = 0.92 & P = 9E-07), iri (R = 0.90 & P = 4E-06), tet34 (R = 0.74 & P = 0.002), tet43 (R = 0.74 & P = 0.002), and mdtD (R = 0.72 & P = 0.003) (Fig. 5A). The co-expressed membership in MEturquoise module and gene significance for tcmA displayed highest relationship (R = 0.99, P = 3.3E-79) (Fig. 5B).

Fig. 5
figure 5

The relationship between unique ARGs and bacterial species in the rumen of dairy buffaloes. A The relationship between 18 unique ARGs and co-expressed ARGs in the rumen of dairy buffaloes. B The relationship between members in significantly correlated gene module and gene significance for tcmA. C The co-occurrence network of tcmA and its highly associated bacteria species in the rumen of dairy buffaloes

Full size image

The co-occurrence patterns between bacterial species in the rumen of dairy buffaloes and tcmA gene were further explored using network analysis (Fig. 5C). Totally, 65 bacterial species were highly positively correlated (R > 0.95, P < 0.01) with tcmA, with 44 belonging to Lactobacillus genus (Fig. 5C). L. amylovorus (R = 0.993, R = 2.77 E-13) and L. acidophilus (R = 0.992, R = 4.94 E-13) showed strongest positive association with tcmA (Fig. 5C). It is notable that L. amylovorus (FC = 269.9, Padjust < 0.001) and L. acidophilus (FC = 11.84, Padjust < 0.001) were highly enriched in the rumen of dairy buffaloes than in dairy cows.

Discussion

To our knowledge, this study should be one of the first study to compare rumen microbiome between dairy cows and buffaloes using metagenomics. Many factors can affect and shape the rumen microbiome, such as diet, breed, age, etc. It is reported that age plays important roles in rumen microbiome dynamics [38]. Dairy buffaloes with different ages (from 1 to 10 year old) was selected in this study to reveal the comprehensive profiles of rumen microbiome as well as to explore the age effects on resistome. Using q-PCR and 16S rRNA gene sequencing, it has been reported that Bacteroidetes (accounting for 42–72% of total bacteria) and Prevotella (accounting for 22–58% of total bacteria) were the most abundant phylum and genus bacteria in dairy buffaloes [39], which were consistent with our results (accounting for 37–60% and 22–44% of total bacteria). We also found that the top abundant bacterial phyla and families were consistent in dairy cows and dairy buffaloes. Even though dairy buffaloes showed distinct microbial community composition and higher Shannon diversity from dairy cows under the similar diet [18], the major rumen microbial profiles especially at phylum level show strong similarities between cattle and buffaloes [40, 41]. Lin et al. reported that no correlation between rumen microbiota and diet [39] while significant association between rumen microbiota composition and milking performance were observed [42] in dairy buffalo. It is indicated that dairy cows and dairy buffaloes shared some core rumen microbiota regardless of diet and geographical location, but the differed microbiota at higher taxonomic levels may reflect the breed or functional difference. We found the diet-induced SDFs were likely to be included in the comparisons between dairy cows and dairy buffaloes with relative lower strict thresholds, therefore, the SDFs revealed under the cutoff of LDA > 3.5 and Padjust < 0.01 were further investigated.

Prevotella ruminicola was the most abundant bacterial species in our dataset in both dairy cows and buffaloes, consistent with previous reports [43, 44]. Prevotella ruminicola can efficiently utilize ammonia nitrogen and peptides for microbial protein synthesis [45]. The optimal adaptative mechanism of Prevotella ruminicola in the rumen environment makes it maximize nitrogen assimilation [46]. It is reported that the Prevotella ruminicola population increased with the increasing dietary crude protein (CP) levels [47, 48]. However, the dietary CP in dairy buffalo was lower than that in dairy cow (13.8% vs. 16.1%). The opposite results (the Prevotella ruminicola was significantly enriched in dairy buffaloes) suggest that Prevotella ruminicola abundance in the rumen was not only driven by dietary CP but may also be influenced by microbial-microbial interaction and host effect between dairy cows and buffaloes. Another significantly enriched prevotella species (Prevotella sp. P6B4) in dairy buffaloes was also abundant bacteria detected in rumen across the globe [49]. The two species were predicted to have 27 and 26 polysaccharide utilization locus [50], indicating that the rumen bacteria in dairy buffaloes have higher breakdown capacity of complex glycans than that in dairy cows. The higher abundance of Acetobacter pasteurianus and Lactobacillus amylovorus in the rumen contributed to the greater degradation ability of cellulose in dairy buffaloes. The function of Faecalibacterium sp. CAG:74, Bacterium F082 and Bacterium F083 species has not been well characterized, nevertheless, the high metabolic activity of these bacteria was approved [51]. The potential higher levels of metabolic function and carbohydrate degradation in the rumen microbiome of dairy buffaloes were supported by functional analysis.

More than 62% of significantly different metabolic pathways and CAZymes between two groups were enriched in dairy buffaloes, which provides further evidence that dairy buffaloes were enhanced in metabolism, especially for carbohydrates degradation. The enrichment of genes encoding CAZymes involved in different types of enzymes revealed the higher capability of degrading diverse and complex substrates in the rumen of dairy buffaloes. The GHs are common enzymes in rumen that catalyze the hydrolysis of glycosidic bond in complex carbohydrates, which are involved in deconstructing plant biomass such as cellulose, hemicellulose, and starch [52]. The large amount of significantly enriched GHs identified in the rumen microbiome of dairy buffaloes suggests the strong fiber degradation roles in the rumen of dairy buffaloes. The extremely more abundant GHs families (GH68, GH70, GH71) with over 50-fold enrichment revealed higher ability of acting on sucrose, maltodextrins/starch, and energy production in dairy buffaloes [53, 54]. It has been suggested that utilizing microbial enzymes belonging to AA family could help degrade the non-carbohydrate structural components (e.g. lignin) since the AAs contain a group of ligninolytic enzymes or multi-copper oxidases [55]. The AA1 family plays important roles in lignin degradation by acting as laccase [56]. The AA1 family-encoded genes only identified in the rumen of dairy buffaloes with a relative high abundance (CPM = 185) indicates that dairy buffaloes might be more capable or more efficient in degrading the recalcitrant lignin, which contribute to the better adaptability to less-digestible forage.

Liu et al. investigated the fecal resistome of dairy cattle and identified 329 ARGs that putatively conferred resistance to 17 classes of antimicrobials [5]. Recent studies reported that the rumen microbiome carries abundant ARGs, which may be transmitted to humans via environment (i.e. soil and water) through saliva or the flow of rumen microbial biomass to the gut [11]. In this study, more than 50% of ARGs and resistant antimicrobials classes were detected in the rumen of dairy cows and buffaloes, suggesting the necessity to pay attention to the ruminal resistome, especially for those that confer resistance to clinically important antibiotics. Mcr-1 is a plasmid-mediated colistin resistance gene that was firstly discovered in the pork, chicken and human in China, and then identified globally [57], which compromise last-line treatment options for human health especially for the multidrug-resistant Gram-negative pathogen infections [58]. It is suggested that the spread of mcr-1 is from livestock sector to human beings [2]. The high mcr-1 positive detection rate in the rumen of dairy cows (81%) and dairy buffaloes (100%) indicates the mcr-1 widespread in dairy cows and buffaloes in the south of China, thus the use of colistin in dairy industry should be reconsidered. Linezolid is one of the last antimicrobial treatment options in human clinical medicine and has not been approved for applying in livestock industry [6]. However, several linezolid resistant genes (cfrA, cipA, clbA, clbB, clbC) with high relative abundance (average CPM > 600 for each gene) were identified in the rumen samples of dairy cows and dairy buffaloes, and they also confirmed in other studies [59, 60]. The high prevalence of ARGs that conferred resistance to clinical-important antibiotics in the rumen microbiota should be taken into account to address potential concerns for public health. The microbial genomes contain functional integrative conjugative elements that determine the possibility of transferring ARGs into human pathogens or zoonotic pathogens [13]. Further work is needed to investigate these aspects.

It has been reported that diets [5], ages [61], and living environment [12] play important roles in shaping gut resistome in animals. The distinct resistome between dairy cows and buffaloes indentified in this study may be attributed to the combined effects of distinct diets, ages, environmental conditions and host genetics. Among them, the diet of dairy buffaloes hade lower nutritional level (e.g., CP) and less abundant ingredient types. It’s noted that 60% corn silage was utilized in the dairy buffaloes’ diet while only 16.7% used in the dairy cows’ diet (Table S1), which may contribute the more diverse resistome in the rumen of dairy buffaloes. The improved adaptability to less-digestible roughage in dairy buffaloes [16] may harbor varied microbiota with diverse ARGs. In addition, the close interactions between dairy buffaloes and waterlogged conditions [19] increase the ARGs horizontal transfer chances, leading to a more abundant resistome. Interestingly, the 18 unique ARGs were detected in each ruminal sample of dairy buffaloes, but they were absent in any ruminal sample from dairy cows. This is also confirmed in another metagenomic dataset (unpublished) of dairy cows. This may be attributed to the distinct rumen microbiome between two host species since microbial resistome are mainly structured by the bacterial phylogeny [8]. These 18 unique ARGs have not been detected in the ruminal resistome of different beef cattle species with varied diets in a previous study [11]. In another study, Cao et al. collected gut samples from 10 different migratory bird species (total number = 99) and characterized their microbiome and resistome using metagenomics [62]. Four dairy buffalo unique ARGs (adeA, sul3, mexC, tet41) were identified in this dataset, but the detection rate was only 1–14% and none of them was distributed in every sample of one certain species. When compared with the study on the swine farm related samples (human feces, pig feces, soil, sewage, and ventilation dust) [63], three (tet34, sul3, iri) buffalo unique ARGs were also identified. Based on the above comparison with more diverse gut microbiome, most of unique ARGs in dairy buffaloes were conserved, indicating the uniqueness of ruminal resistome and ARGs profiles.

Among the unique ARGs, tcmA was highlighted since it showed greatest positive correlation with a gene module consisting of 93 co-expressed ARGs. The tcmA is predicted to confer resistance to tetracenomycin C via an active efflux system [64]. Tetracenomycin C is a narrow spectrum anthracycline antibiotic, which is active against Gram-positive bacteria and also shows antitumor effect [65]. In this study, we identified more than 40 bacteria under Lactobacillus genus that were highly related to tcmA, suggesting that tcmA might be widely harbored and/or more likely horizontally transferred in Lactobacillus bacteria. Among them, L. amylovorus and L. acidophilus showed strongest association with tcmA, which is supported by the taxonomic alignment using ARGs-containing contigs. The lactic bacteria Lactobacillus amylovorus was firstly isolated from waste-corn fermentations in cattle, and can secret starch-hydrolyzing enzyme (a-amylase) and rapidly digest cornstarch [66]. It has been found that L. amylovorus and L. acidophilus exert beneficial effects on gastrointestinal immunomodulation and health as candidate bacterial therapeutics, such as mitigating the bovine respiratory pathogens in dairy calves [67]. The probiotic properties of L. amylovorus and L. acidophilus have been well studied [68, 69], suggesting the promising cholesterol-lowering properties [70], anti-pathogens effect [71], and prevention of disease in human and dairy calves [72]. The significantly higher relative abundance of these two species especially for L. amylovorus (2% in young dairy buffaloes) in the rumen of dairy buffaloes may provide probiotic resource for fermented feed. The detected tcmA will provide fundamental knowledge and potential strategies to reduce ARG transmission risk when isolated Lactobacillus probiotic from the rumen of dairy buffaloes.

Conclusions

Ruminal bacterial composition varied significantly between dairy cows and dairy buffaloes. The majority of differed KEGG pathways and CAZymes were enriched in the dairy buffaloes with GH68, GH70, GH71, GT48 and AA1were extremely significantly higher than dairy cows. Dairy buffaloes showed distinct resistome profiles and clusters from dairy cows. A total of 18 ARGs conferring resistance to 16 antibiotic classes were uniquely detected in the dairy buffaloes. tcmA plays a central role in the ruminal resistome and is highly associated with L. amylovorus and L. acidophilus. Our results provide novel insights into the microbiome and resistome of dairy buffaloes, the identified unique ARGs and associated bacteria will help develop strategies to prevent the transmission of ARGs.

Availability of data and materials

The metagenomics raw sequences of rumen samples were submitted into the NCBI Sequence Read Archive (SRA) under the accession numbers PRJNA718720, other relevant data are available upon request.

Abbreviations

AAs:

Auxiliary activities

AMR:

Antimicrobial resistance

ARGs:

Antimicrobial resistance genes

BWA:

Burrows-Wheeler alignment

CARD:

Comprehensive Antimicrobial Resistance Database

CAZy:

Carbohydrate-Active enZYmes Database

CBMs:

Carbohydrate-binding modules

CPM:

Counts per million reads

FDR:

False discovery rate

FC:

Fold change

GHs:

Glycoside hydrolases

GTs:

Glycosyltransferase

KEGG:

Kyoto Encyclopedia of Genes and Genomes database

LDA:

Least discriminant analysis

LEfSe:

LDA effect size

NR:

Non-redundant

PCoA:

Principal-coordinate analysis

PLs:

Polysaccharide lyases

PLS-DA:

Partial least-squares discrimination analysis

SDFs:

Significantly different features

WGCNA:

Weighted gene co-expression network analysis

References

  1. Sommer MO, Dantas G, Church GM. Functional characterization of the antibiotic resistance reservoir in the human microflora. Science. 2009;325(5944):1128–31. https://doi.org/10.1126/science.1176950.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Liu Y-Y, Wang Y, Walsh TR, Yi L-X, Zhang R, Spencer J, et al. Emergence of plasmid-mediated colistin resistance mechanism MCR-1 in animals and human beings in China: a microbiological and molecular biological study. Lancet Infect Dis. 2016;16(2):161–8. https://doi.org/10.1016/S1473-3099(15)00424-7.

    Article  CAS  PubMed  Google Scholar 

  3. Munk P, Knudsen BE, Lukjancenko O, Duarte ASR, Van Gompel L, Luiken RE, et al. Abundance and diversity of the faecal resistome in slaughter pigs and broilers in nine European countries. Nat Microbiol. 2018;3(8):898–908. https://doi.org/10.1038/s41564-018-0192-9.

    Article  CAS  PubMed  Google Scholar 

  4. Wang Y, Hu Y, Liu F, Cao J, Lv N, Zhu B, et al. Integrated metagenomic and metatranscriptomic profiling reveals differentially expressed resistomes in human, chicken, and pig gut microbiomes. Environ Int. 2020;138:105649. https://doi.org/10.1016/j.envint.2020.105649.

    Article  CAS  PubMed  Google Scholar 

  5. Liu J, Taft DH, Maldonado-Gomez MX, Johnson D, Treiber ML, Lemay DG, et al. The fecal resistome of dairy cattle is associated with diet during nursing. Nat Commun. 2019;10(1):4406. https://doi.org/10.1038/s41467-019-12111-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Wang Y, Li X, Fu Y, Chen Y, Wang Y, Ye D, et al. Association of florfenicol residues with the abundance of oxazolidinone resistance genes in livestock manures. J Hazard Mater. 2020;399:123059. https://doi.org/10.1016/j.jhazmat.2020.123059.

    Article  CAS  PubMed  Google Scholar 

  7. He T, Wang R, Liu D, Walsh TR, Zhang R, Lv Y, et al. Emergence of plasmid-mediated high-level tigecycline resistance genes in animals and humans. Nat Microbiol. 2019;4(9):1450–6. https://doi.org/10.1038/s41564-019-0445-2.

    Article  CAS  PubMed  Google Scholar 

  8. Forsberg KJ, Patel S, Gibson MK, Lauber CL, Knight R, Fierer N, et al. Bacterial phylogeny structures soil resistomes across habitats. Nature. 2014;509(7502):612–6. https://doi.org/10.1038/nature13377.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Liu JH, Zhang ML, Zhang RY, Zhu WY, Mao SY. Comparative studies of the composition of bacterial microbiota associated with the ruminal content, ruminal epithelium and in the faeces of lactating dairy cows. Microbial Biotechnol. 2016;9(2):257–68. https://doi.org/10.1111/1751-7915.12345.

    Article  CAS  Google Scholar 

  10. Guo L, Yao D, Li D, Lin Y, Bureenok S, Ni K, et al. Effects of lactic acid Bacteria isolated from rumen fluid and feces of dairy cows on fermentation quality, microbial community, and in vitro digestibility of alfalfa silage. Front Microbiol. 2020;10:2998. https://doi.org/10.3389/fmicb.2019.02998.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Auffret MD, Dewhurst RJ, Duthie CA, Rooke JA, John Wallace R, Freeman TC, et al. The rumen microbiome as a reservoir of antimicrobial resistance and pathogenicity genes is directly affected by diet in beef cattle. Microbiome. 2017;5(1):159. https://doi.org/10.1186/s40168-017-0378-z.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Campbell TP, Sun X, Patel VH, Sanz C, Morgan D, Dantas G. The microbiome and resistome of chimpanzees, gorillas, and humans across host lifestyle and geography. ISME J. 2020;14(6):1584–99. https://doi.org/10.1038/s41396-020-0634-2.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Sabino YNV, Santana MF, Oyama LB, Santos FG, Moreira AJS, Huws SA, et al. Characterization of antibiotic resistance genes in the species of the rumen microbiota. Nat Commun. 2019;10(1):5252. https://doi.org/10.1038/s41467-019-13118-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Singh B, Mal G, Kues WA, Yadav PS. The domesticated buffalo-an emerging model for experimental and therapeutic use of extraembryonic tissues. Theriogenol. 2020;151:95102.

    Article  Google Scholar 

  15. Aspilcueta-Borquis R, Neto FA, Baldi F, Bignardi A, Albuquerque LG, Tonhati H. Genetic parameters for buffalo milk yield and milk quality traits using Bayesian inference. J Dairy Sci. 2010;93(5):2195–201. https://doi.org/10.3168/jds.2009-2621.

    Article  CAS  PubMed  Google Scholar 

  16. Hamid M, Siddiky M, Rahman M, Hossain K. Scopes and opportunities of buffalo farming in Bangladesh: a review. SAARC J Agri. 2016;14(2):63–77.

    Article  Google Scholar 

  17. El Debaky HA, Kutchy NA, Ul-Husna A, Indriastuti R, Akhter S, Purwantara B, et al. Potential of water buffalo in world agriculture: challenges and opportunities. Appl Anim Sci. 2019;35(2):255–68. https://doi.org/10.15232/aas.2018-01810.

    Article  Google Scholar 

  18. Iqbal MW, Zhang Q, Yang Y, Li L, Zou C, Huang C, et al. Comparative study of rumen fermentation and microbial community differences between water buffalo and Jersey cows under similar feeding conditions. J Appl Anim Res. 2018;46(1):740–8. https://doi.org/10.1080/09712119.2017.1394859.

    Article  Google Scholar 

  19. De Rosa G, Grasso F, Braghieri A, Bilancione A, Di Francia A, Napolitano F. Behavior and milk production of buffalo cows as affected by housing system. J Dairy Sci. 2009;92(3):907–12. https://doi.org/10.3168/jds.2008-1157.

    Article  CAS  PubMed  Google Scholar 

  20. Shen JS, Chai Z, Song LJ, Liu JX, Wu YM. Insertion depth of oral stomach tubes may affect the fermentation parameters of ruminal fluid collected in dairy cows. J Dairy Sci. 2012;95(10):5978–84. https://doi.org/10.3168/jds.2012-5499.

    Article  CAS  PubMed  Google Scholar 

  21. Yu Z, Morrison M. Improved extraction of PCR-quality community DNA from digesta and fecal samples. BioTechniques. 2004;36(5):808–12. https://doi.org/10.2144/04365ST04.

    Article  CAS  PubMed  Google Scholar 

  22. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90. https://doi.org/10.1093/bioinformatics/bty560.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Li H, Durbin R. Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doi.org/10.1093/bioinformatics/btp324.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Liu C-M, Li D, Sadakane K, Luo R, Lam T-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31(10):1674–6. https://doi.org/10.1093/bioinformatics/btv033.

    Article  CAS  PubMed  Google Scholar 

  25. Noguchi H, Park J, Takagi T. MetaGene: prokaryotic gene finding from environmental genome shotgun sequences. Nucleic Acids Res. 2006;34(19):5623–30. https://doi.org/10.1093/nar/gkl723.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Niu B, Fu L, Wu S, Li W, Zhu Z. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012;28(23):3150–2.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Yu C, Wang J, Kristiansen K, Li R, Yiu S-M, Lam T-W, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25(15):1966–7.

    Article  PubMed  Google Scholar 

  28. Alcock BP, Raphenya AR, Lau TTY, Tsang KK, Bouchard M, Edalatmand A, et al. CARD 2020: antibiotic resistome surveillance with the comprehensive antibiotic resistance database. Nucleic Acids Res. 2019;48(D1):D517–25.

    PubMed Central  Google Scholar 

  29. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(suppl_2):W316–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Lombard V, Golaconda Ramulu H, Drula E, Coutinho PM, Henrissat B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 2014;42(D1):D490–5. https://doi.org/10.1093/nar/gkt1178.

    Article  CAS  PubMed  Google Scholar 

  31. Thévenot EA, Roux A, Xu Y, Ezan E, Junot C. Analysis of the human adult urinary metabolome variations with age, body mass index, and gender by implementing a comprehensive workflow for univariate and OPLS statistical analyses. J Proteome Res. 2015;14(8):3322–35. https://doi.org/10.1021/acs.jproteome.5b00354.

    Article  CAS  PubMed  Google Scholar 

  32. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559. https://doi.org/10.1186/1471-2105-9-559.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Sun H-Z, Zhou M, Wang O, Chen Y, Liu J-X, Guan LL. Multi-omics reveals functional genomic and metabolic mechanisms of milk production and quality in dairy cows. Bioinformatics. 2020;36(8):2530–7. https://doi.org/10.1093/bioinformatics/btz951.

    Article  CAS  PubMed  Google Scholar 

  34. Legeay M, Doncheva NT, Morris JH, Jensen LJ. Cytoscape app [version 1; peer review: 1 approved]; 2020.

    Google Scholar 

  35. Zapala MA, Schork NJ. Multivariate regression analysis of distance matrices for testing associations between gene expression patterns and related variables. P Natl Acad Sci USA. 2006;103(51):19430–5. https://doi.org/10.1073/pnas.0609333103.

    Article  CAS  Google Scholar 

  36. Oksanen J, Blanchet FG, Kindt R, Legendre P, O’hara R, Simpson GL, et al. Vegan: community ecology package. R package version 1; 2010. p. 17–4. URL http://CRAN.R-project.org/package=vegan.

    Google Scholar 

  37. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genom Biol. 2011;12(6):R60. https://doi.org/10.1186/gb-2011-12-6-r60.

    Article  Google Scholar 

  38. Furman O, Shenhav L, Sasson G, Kokou F, Honig H, Jacoby S, et al. Stochasticity constrained by deterministic effects of diet and age drive rumen microbiome assembly dynamics. Nat Commun. 2020;11(1):1904. https://doi.org/10.1038/s41467-020-15652-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Lin B, Henderson G, Zou C, Cox F, Liang X, Janssen PH, et al. Characterization of the rumen microbial community composition of buffalo breeds consuming diets typical of dairy production systems in southern China. Anim Feed Sci Tech. 2015;207:75–84. https://doi.org/10.1016/j.anifeedsci.2015.06.013.

    Article  CAS  Google Scholar 

  40. Imai S, Abdullah N, Ho Y, Jalaludin S, Hussain H, Onodera R, et al. Comparative study on the rumen ciliate populations in small experimental herds of water buffalo and Kedah Kelantan cattle in Malaysia. Anim Feed Sci Tech. 1995;52(3–4):345–51. https://doi.org/10.1016/0377-8401(94)00726-P.

    Article  Google Scholar 

  41. Lwin K-O, Matsui H, Ban-Tokuda T, Kondo M, Lapitan RM, Herrera JRV, et al. Comparative analysis of methanogen diversity in the rumen of crossbred buffalo and cattle in the Philippines by using the functional gene mcrA. Mol Biol Rep. 2012;39(12):10769–74. https://doi.org/10.1007/s11033-012-1969-1.

    Article  CAS  PubMed  Google Scholar 

  42. Zou C, Gu Q, Zhou X, Xia Z, Muhammad WI, Tang Q, et al. Ruminal microbiota composition associated with ruminal fermentation parameters and milk yield in lactating buffalo in Guangxi, China: a preliminary study. J Anim Physiol Anim Nutr. 2019;103(5):1374–9. https://doi.org/10.1111/jpn.13154.

    Article  CAS  Google Scholar 

  43. Costa-Roura S, Balcells Terés J, de la Fuente OG, Mora-Gil J, Llanes N, Villalba Mata D. Effects of protein restriction on performance, ruminal fermentation and microbial community in Holstein bulls fed high-concentrate diets. Anim Feed Sci Tech. 2020;264:114479. https://doi.org/10.1016/j.anifeedsci.2020.114479.

    Article  CAS  Google Scholar 

  44. Koringa PG, Thakkar JR, Pandit RJ, Hinsu AT, Parekh MJ, Shah RK, et al. Metagenomic characterisation of ruminal bacterial diversity in buffaloes from birth to adulthood using 16S rRNA gene amplicon sequencing. Funct Integr Genomic. 2019;19(2):237–47. https://doi.org/10.1007/s10142-018-0640-x.

    Article  CAS  Google Scholar 

  45. Purushe J, Fouts DE, Morrison M, White BA, Mackie RI, Coutinho PM, et al. Bacteria NACfR: comparative genome analysis of Prevotella ruminicola and Prevotella bryantii: insights into their environmental niche. Microb Ecol. 2010;60(4):721–9. https://doi.org/10.1007/s00248-010-9692-8.

    Article  PubMed  Google Scholar 

  46. Kim JN, Méndez-García C, Geier RR, Iakiviak M, Chang J, Cann I, et al. Metabolic networks for nitrogen utilization in Prevotella ruminicola 23. Sci Rep. 2017;7(1):7851. https://doi.org/10.1038/s41598-017-08463-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Wang C, Liu Q, Guo G, Huo W, Liang Y, Pei C, et al. Effects of different dietary protein levels and rumen-protected folic acid on ruminal fermentation, degradability, bacterial populations and urinary excretion of purine derivatives in beef steers. J Agri Sci. 2017;155(9):1477–86. https://doi.org/10.1017/S0021859617000533.

    Article  CAS  Google Scholar 

  48. Liu Q, Wang C, Li H, Guo G, Huo W, Pei C, et al. Effects of dietary protein levels and rumen-protected pantothenate on ruminal fermentation, microbial enzyme activity and bacteria population in blonde d'Aquitaine× Simmental beef steers. Anim Feed Sci Tech. 2017;232:31–9. https://doi.org/10.1016/j.anifeedsci.2017.07.014.

    Article  CAS  Google Scholar 

  49. Henderson G, Cox F, Ganesh S, Jonker A, Young W, Abecia L, et al. Rumen microbial community composition varies with diet and host, but a core microbiome is found across a wide geographical range. Sci Rep. 2015;5(1):14567. https://doi.org/10.1038/srep14567.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Terrapon N, Lombard V, Drula E, Lapébie P, Al-Masaudi S, Gilbert HJ, et al. PULDB: the expanded database of polysaccharide utilization loci. Nucleic Acids Res. 2018;46(D1):D677–83. https://doi.org/10.1093/nar/gkx1022.

    Article  CAS  PubMed  Google Scholar 

  51. Wirth R, Kádár G, Kakuk B, Maróti G, Bagi Z, Szilágyi Á, et al. The planktonic core microbiome and core functions in the cattle rumen by next generation sequencing. Front Microbiol. 2018;9:2285. https://doi.org/10.3389/fmicb.2018.02285.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Brulc JM, Antonopoulos DA, Miller MEB, Wilson MK, Yannarell AC, Dinsdale EA, et al. Gene-centric metagenomics of the fiber-adherent bovine rumen microbiome reveals forage specific glycoside hydrolases. P Natl Acad Sci USA. 2009;106(6):1948–53. https://doi.org/10.1073/pnas.0806191105.

    Article  Google Scholar 

  53. Gangoiti J, Van Leeuwen SS, Gerwig GJ, Duboux S, Vafiadi C, Pijning T, et al. 4, 3-α-Glucanotransferase, a novel reaction specificity in glycoside hydrolase family 70 and clan GH-H. Sci Rep. 2017;7(1):39761. https://doi.org/10.1038/srep39761.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Cantarel BL, Lombard V, Henrissat B. Complex carbohydrate utilization by the healthy human microbiome. PLoS One. 2012;7(6):e28742. https://doi.org/10.1371/journal.pone.0028742.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. He B, Jin S, Cao J, Mi L, Wang J. Metatranscriptomics of the Hu sheep rumen microbiome reveals novel cellulases. Biotechnol Biofuels. 2019;12(1):153. https://doi.org/10.1186/s13068-019-1498-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Park Y-J, Kong W-S. Genome-wide comparison of carbohydrate-active enzymes (CAZymes) repertoire of Flammulina ononidis. Mycobiol. 2018;46(4):349–60. https://doi.org/10.1080/12298093.2018.1537585.

    Article  Google Scholar 

  57. Shen Z, Wang Y, Shen Y, Shen J, Wu C. Early emergence of mcr-1 in Escherichia coli from food-producing animals. Lancet Infect Dis. 2016;16(3):293. https://doi.org/10.1016/S1473-3099(16)00061-X.

    Article  PubMed  Google Scholar 

  58. Wang Y, Tian G-B, Zhang R, Shen Y, Tyrrell JM, Huang X, et al. Prevalence, risk factors, outcomes, and molecular epidemiology of mcr-1-positive Enterobacteriaceae in patients and healthy adults from China: an epidemiological and clinical study. Lancet Infect Dis. 2017;17(4):390–9. https://doi.org/10.1016/S1473-3099(16)30527-8.

    Article  CAS  PubMed  Google Scholar 

  59. Zhao Q, Wang Y, Wang S, Wang Z, Du X-d, Jiang H, et al. Prevalence and abundance of florfenicol and linezolid resistance genes in soils adjacent to swine feedlots. Sci Rep. 2016;6:32192.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Li J, Shao B, Shen J, Wang S, Wu Y. Occurrence of chloramphenicol-resistance genes as environmental pollutants from swine feedlots. Environ Sci Technol. 2013;47(6):2892–7. https://doi.org/10.1021/es304616c.

    Article  CAS  PubMed  Google Scholar 

  61. Gasparrini AJ, Wang B, Sun X, Kennedy EA, Hernandez-Leyva A, Ndao IM, et al. Persistent metagenomic signatures of early-life hospitalization and antibiotic treatment in the infant gut microbiota and resistome. Nat Microbiol. 2019;4(12):2285–97. https://doi.org/10.1038/s41564-019-0550-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Cao J, Hu Y, Liu F, Wang Y, Bi Y, Lv N, et al. Metagenomic analysis reveals the microbiome and resistome in migratory birds. Microbiome. 2020;8(1):26. https://doi.org/10.1186/s40168-019-0781-8.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Sun J, Liao X-P, D’Souza AW, Boolchandani M, Li S-H, Cheng K, et al. Environmental remodeling of human gut microbiota and antibiotic resistome in livestock farms. Nat Commun. 2020;11(1):1427. https://doi.org/10.1038/s41467-020-15222-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Guilfoile P, Hutchinson C. Sequence and transcriptional analysis of the Streptomyces glaucescens tcmAR tetracenomycin C resistance and repressor gene loci. J Bacteriol. 1992;174(11):3651–8. https://doi.org/10.1128/jb.174.11.3651-3658.1992.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Blanco G, Patallo EP, Braña AF, Trefzer A, Bechthold A, Rohr J, et al. Identification of a sugar flexible glycosyltransferase from Streptomyces olivaceus, the producer of the antitumor polyketide elloramycin. Chem Biol. 2001;8(3):253–63. https://doi.org/10.1016/S1074-5521(01)00010-2.

    Article  CAS  PubMed  Google Scholar 

  66. Nakamura L. Lactobacillus amylovorus, a new starch-hydrolyzing species from cattle waste-corn fermentations. Int J Syst Evol Microbiol. 1981;31(1):56–63.

    CAS  Google Scholar 

  67. Amat S, Alexander TW, Holman DB, Schwinghamer T, Timsit E. Intranasal bacterial therapeutics reduce colonization by the respiratory pathogen Mannheimia haemolytica in dairy calves. Msystems. 2020;5(2):e00629–19.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Aristimuño Ficoseco C, Mansilla FI, Maldonado NC, Miranda H, Fátima Nader-Macias ME, Vignolo GM. Safety and growth optimization of lactic acid bacteria isolated from feedlot cattle for probiotic formula design. Front Microbiol. 2018;9:2220. https://doi.org/10.3389/fmicb.2018.02220.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Kos B, Šušković J, Vuković S, Šimpraga M, Frece J, Matošić S. Adhesion and aggregation ability of probiotic strain Lactobacillus acidophilus M92. J Appl Microbiol. 2003;94(6):981–7. https://doi.org/10.1046/j.1365-2672.2003.01915.x.

    Article  CAS  PubMed  Google Scholar 

  70. Omar JM, Chan Y-M, Jones ML, Prakash S, Jones PJ. Lactobacillus fermentum and Lactobacillus amylovorus as probiotics alter body adiposity and gut microflora in healthy persons. J Funct Foods. 2013;5(1):116–23. https://doi.org/10.1016/j.jff.2012.09.001.

    Article  Google Scholar 

  71. Adetoye A, Pinloche E, Adeniyi BA, Ayeni FA. Characterization and anti-salmonella activities of lactic acid bacteria isolated from cattle faeces. BMC Microbiol. 2018;18(1):96. https://doi.org/10.1186/s12866-018-1248-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Fernández S, Fraga M, Silveyra E, Trombert A, Rabaza A, Pla M, et al. Probiotic properties of native Lactobacillus spp. strains for dairy calves. Benef Microbes. 2018;9(4):613–24. https://doi.org/10.3920/BM2017.0131.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We acknowledge Dr. Bo Lin and his group in Guangxi University (Nanning, China) for their assistance in the sample collection.

Funding

This research was supported by grants from the National Natural Science Foundation of China (32002207, 31872380), Fundamental Research Funds for the Zhejiang Provincial Universities (2021XZZX027), and the “Hundred Talents Program” Research Professor Start-up Fund of Zhejiang University.

Author information

Authors and Affiliations

  1. Institute of Dairy Science, Ministry of Education Key Laboratory of Molecular Animal Nutrition, College of Animal Sciences, Zhejiang University, Hangzhou, 310058, China

    Hui-Zeng Sun, Ke-Lan Peng, Ming-Yuan Xue & Jian-Xin Liu

Authors
  1. Hui-Zeng Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Ke-Lan Peng
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ming-Yuan Xue
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jian-Xin Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JL, HS and MX designed the study. KP collected the rumen samples and extracted DNA. HZ performed the bioinformatics analysis and wrote the manuscript. HZ and JL revised the manuscript. All authors read, revised and approved the final manuscript.

Corresponding author

Correspondence to Jian-Xin Liu.

Ethics declarations

Ethics approval and consent to participate

Animal care and experimental procedures were approved by the Animal Care Committee of Zhejiang University (Hangzhou, China) and were under the university’s guidelines for animal research.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

The diet ingredents of dairy cows and dairy buffaloes. Table S2. The relative abundance of bacteria at domain level in individuals. Table S3. The significantly different CAZymes encoded genes identified in the rumen microbiota between dairy cows and buffaloes. Table S4. ARGs and their abundance in dairy cows and dairy buffaloes. Table S5. ARGs and their AMR phenotypes in the dairy buffaloes. Table S6. The relative abundance of unique ARGs in different age groups of dairy buffaloes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, HZ., Peng, KL., Xue, MY. et al. Metagenomics analysis revealed the distinctive ruminal microbiome and resistive profiles in dairy buffaloes. anim microbiome 3, 44 (2021). https://doi.org/10.1186/s42523-021-00103-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42523-021-00103-6

Keywords

  • Dairy buffaloes
  • Dairy cattle
  • Rumen
  • Antimicrobial resistant genes
  • Microbiome

Animal Microbiome

ISSN: 2524-4671