Skip to main content

Survey of rumen microbiota of domestic grazing yak during different growth stages revealed novel maturation patterns of four key microbial groups and their dynamic interactions

Abstract

Background

The development and maturation of rumen microbiota across the lifetime of grazing yaks remain unexplored due to the varied lifestyles and feed types of yaks as well as the challenges of obtaining samples. In addition, the interactions among four different rumen microbial groups (bacteria, archaea, fungi and protozoa) in the rumen of yak are not well defined. In this study, the rumen microbiota of full-grazing yaks aged 7 days to 12 years old was assessed to determine the maturation patterns of these four microbial groups and the dynamic interactions among them during different growth stages.

Results

The rumen microbial groups (bacteria, archaea, protozoa and fungi) varied through the growth of yaks from neonatal (7 days) to adult (12 years), and the bacterial and archaeal groups were more sensitive to changes in growth stages compared to the two eukaryotic microbial groups. The age-discriminatory taxa within each microbial group were identified with the random forest model. Among them, Olsenella (bacteria), Group 10 sp., belonging to the family Methanomassiliicoccaceae (archaea), Orpinomyces (fungi), and Dasytricha (protozoa) contributed the most to discriminating the age of the rumen microbiota. Moreover, we found that the rumen archaea reached full maturation at 5 approximately years of age, and the other microbial groups matured between 5 and 8 years of age. The intra-interactions patterns and keystone species within each microbial group were identified by network analysis, and the inter-interactions among the four microbial groups changed with growth stage. Regarding the inter-interactions among the four microbial groups, taxa from bacteria and protozoa, including Christensenellaceae R-7 group, Prevotella 1, Trichostomatia, Ruminococcaceae UCG-014 and Lachnospiraceae, were the keystone species in the network based on betweenness centrality scores.

Conclusions

This study depicted a comprehensive view of rumen microbiota changes in different growth stages of grazing yaks. The results revealed the unique microbiota maturation trajectory and the intra- and inter-interactions among bacteria, archaea, fungi and protozoa in the rumen of grazing yaks across the lifetime of yaks. The information obtained in this study is vital for the future development of strategies to manipulate rumen microbiota in grazing yaks for better growth and performance in the harsh Qinghai-Tibetan Plateau ecosystem.

Background

The yak (Bos grunniens), inhabiting the Qinghai-Tibetan Plateau in China, diverged from other ruminant animals millions of years ago [1]. The Qinghai-Tibetan Plateau is the world’s highest plateau with altitudes ranging from 4000 to 5500 m [2], featuring an extreme environment with low ambient temperature and partial pressure of oxygen, as well as a high level of ultraviolet radiation [3]. Yak has developed many anatomical and physiological traits to adapt to this extreme living habitat, including large lungs and hearts, increased foraging ability, high energy metabolism, and lack of hypoxic pulmonary vasoconstriction [4, 5]. Traditionally, yaks graze yearly on native pastures with coarse grasses as their only food source [6]. Although they suffer from inadequate feed and malnutrition in the long and cold seasons [7], they degrade lignocellulose better and have more efficient energy (producing more short chain fatty acids in the rumen) and nitrogen metabolism (higher nitrogen retention) compared to cattle under the same conditions [8, 9].

Similar to other ruminants, the rumen of yak contains complex symbiotic microorganisms that ferment fibrous plant materials, providing the host with usable nutrients such as volatile fatty acids (VFAs) and microbial proteins [10]. A recent study revealed that the proportion of uncultured microbial species was higher in the rumen of naturally grazing yaks than in the rumen of house-farmed Jinnan cattle fed diets with similar energy densities or with high-energy diets [11], and yaks had a unique rumen archaeal community in comparison to domestic cattle under the same grazing conditions [12]. Moreover, it has been reported that individuality [13], feed type [14] and feeding regimes [15, 16] can affect bacterial and archaeal composition in the rumen of yak. Recent research based on rumen metagenomes and host transcriptomes revealed that the rumen of yak contained fewer microbial genes involved in methane emission and a higher abundance of genes involved in host VFA absorption compared with its low-altitude relatives (cattle), suggesting that rumen microbiota co-evolved with the host genome to adapt to extreme environmental conditions [17]. However, most studies mentioned above have only focused on bacterial and archaeal domains and one sampling timepoint. It is known that eukaryotic microbes in the rumen also play important roles in rumen fermentation. Rumen fungi play a significant role in lignocellulosic material degradation, providing a source of fermentable sugars for other microbes and the host [18]. Indeed, it has been reported that the abundance of rumen protozoa was positively associated with total VFA production and the molar proportion of butyrate in the rumen of lambs [19]. Thus, we comprehensively investigated rumen microbial composition (including bacteria, archaea, protozoa and fungi) and their interactions with each other, which is essential for developing appropriate manipulation strategies that promote yaks to better adapt to the extreme grazing environment.

Increasing evidence has highlighted the importance of early-life microbial colonization and its impact on lifelong animal productivity and health [20,21,22,23,24,25,26,27]. It has been shown that the rumen bacterial community is established before the intake of solid food and its composition changes with age [20, 21], and the rumen microbiota (bacteria, archaea and fungi) exhibits an adult-like microbiota between weaning and 1 year of age in dairy cows [22]. In addition, studies have reported that rumen bacterial and archaeal communities continue to evolve and mature with aging even after adulthood [23, 24]. Studies on early-life interventions of rumen microbial community colonization have found that different diets affect microbiota (bacteria, archaea and fungi) development in calves, and some calf-diet-driven differences are apparent in the microbiota of adult cows [25, 26]. For example, feeding a methanogen inhibitor (bromochloromethane) to goat kids affected rumen archaeal community colonization, and the effect of this inhibitor on some less abundant archaeal groups persisted 4 months after exposure [27]. These findings suggest that early-life interventions could lead to assembly of a specific composition or promote the development of rumen that potentially persists later in life, affecting health and productivity [28]. However, there is a paucity of knowledge of how rumen microbiota develops across different growth stages and when it becomes fully maturated over the lifetime of grazing yaks. Therefore, in the present study, we assessed the rumen microbiota of grazing yaks as a whole, including bacteria, archaea, protozoa and fungi, from 7 days to 12 years of age using amplicon sequencing, with the aim of addressing the above knowledge gaps. In addition, the dynamic interactions within and among different microbial groups under different growth stages were also explored to identify whether the compositional changes at each microbial group level could affect the microbial-microbial interactions. Understanding the phylogenetic composition, interactions and maturation of the whole rumen microbiota across the lifetime of yaks will provide the fundamental knowledge needed for developing strategies to improve rumen function in the extreme Qinghai-Tibetan Plateau environment.

Results

Diversity of rumen microbiota of grazing yaks in different growth stages

Amplicon sequencing of bacteria, archaea, protozoa, and fungi in 100 rumen samples collected from yaks aged 7 days to 12 years old across 14 time points (Fig. 1) was performed; however, high-quality sequence data were obtained for 80, 75, 71, and 65 samples, for bacteria, archaea, protozoa, and fungi, respectively. Detailed information on the number of samples is listed in Table 1. Here, 7 days to 1 year of age was considered the pre-weaned stage; 2 to 3 years of age was considered the youth stage (puberty); 5 to 12 years of age was considered the adult stage based on a previous description [29].

Fig. 1
figure 1

Experimental design for rumen sample collection

Table 1 Sample information and sequencing statistics obtained using the DADA2 algorithm in QIIME2

Bacteria were detected at 7 days of age, archaea were detected at 14 days of age, and fungi and protozoa were both detected at 1 month of age (Table 1). For bacteria, 1,296,989 quality-controlled reads were generated, with a mean ± SD of 16,212 (± 7253) sequences per sample, and 26,738 unique exact sequence variants (ESVs) were identified (ranging from 123 to 1437 per sample) (Table 1). Similarly, a total of 581,950 (7759 ± 2364 per sample), 2,824,701 (43,457 ± 10,633 per sample), and 551,668 (7770 ± 4054 per sample) high-quality reads were obtained for archaea, fungi and protozoa, yielding 359 (ranging from 18 to 68), 3230 (ranging from 13 to 115), and 387 unique ESVs (ranging from 7 to 124), respectively. Based on the sampling depth at which rarefaction curves tend to plateau, we rarefied each sample to 7158 (bacteria), 3794 (archaea), 4761 (fungi) and 1818 (protozoa) reads (Additional file 1: Figure S1). The Good’s coverage index was 98.7% (± 0.7%) for bacteria, 99.9% (± 0.1%) for archaea, 99.9% (± 0.03%) for fungi, and 99.7% (± 0.2%) for protozoa (Table 1), indicating that the sequencing depth was adequate to represent each rumen microbial community. To gain insight into the diversity of the rumen microbiota, we compared the Chao1 and Faith’s phylogenetic diversity (PD) indices across age groups. The Chao1 index of the rumen bacteria increased with age up to 2 years old and then remained stable (Fig. 2a), while those of rumen archaea and protozoa were similar up to 2 months, generally increased until the age of 2 years and remained stable afterwards (Fig. 2b and d). However, the Chao1 index of fungi increased from 1 month to 3 months, decreased by the age of 2 years, and further increased with age (Fig. 2c). The PD index of the rumen bacteria and fungi (P < 0.05) significantly increased with age, whereas that of the rumen archaea and protozoa decreased with age, although there were some fluctuations (Fig. 2, Additional file 2: Table S1).

Fig. 2
figure 2

Chao1 index and Faith’s phylogenetic diversity of rumen bacteria (a), archaea (b), fungi (c) and protozoa (d). The table on the right shows the statistical significance

Microbial profiles of rumen microbiota in grazing yaks

Next, principal coordinates analysis (PCoA) based on unweighted and weighted UniFrac distance matrices was performed to determine whether the microbial community structure changed with increasing age. The PCoA plot showed clear age-based separation of rumen bacteria between 7 days to 1 month, 2–4 months, and 6 months to 12 years of age (Fig. 3a). Similarly, three clusters were formed for rumen archaea, where 14 days to 1 month formed one cluster, 2 months to 1 year formed another cluster, and the remaining age groups (2–12 years) formed another cluster (Fig. 3b). For fungi, we observed clear separations between 1 month to 1 year and 2 to 12 years of age (Fig. 3c). However, no clear separation was found for rumen protozoa (Fig. 3d).

Fig. 3
figure 3

Principal coordinate analysis (PCoA) based on weighted and unweighted UniFrac distances. Changes in the rumen bacterial (a), archaeal (b), fungal (c) and protozoal (d) community structure with growth stage

In addition, the permutational analysis of variance (PERMANOVA) showed that profiles of bacteria, archaea and fungi significantly differed among different age groups (Benjamini-Hochberg corrected P < 0.05, Additional file 2: Table S2). For protozoa, it was especially significant between young age groups (1 to 4 months of age) and old age groups (3 to 12 years of age) (Benjamini-Hochberg corrected P < 0.01). The UniFrac dissimilarity of the bacterial community decreased with age (P < 0.01), as the distance between samples decreased with increasing age (Additional file 3: Figure S2A). The highest dissimilarity of the archaeal community was observed at 1 month of age and then decreased with age (P < 0.01, Additional file 3: Figure S2B). However, the decrease in dissimilarity with the increase in age was less apparent for rumen fungi and protozoa (P < 0.05, Additional file 3: Figure S2C and 2D), except for the weighted UniFrac dissimilarity in the fungal community (Additional file 3: Figure S2C).

Shifts in rumen microbiota at different taxonomic levels

Changes in the rumen bacterial composition with yak growth

For bacteria, sequence variants were classified into 368 genera belonging to 26 phyla. At the phylum level, 16 phyla were referred to as the detected bacterial phyla (relative abundance > 0.1% and present in more than half of the total animal populations per age group). The most abundant phylum was Bacteroidetes (41.7 to 74.1%), followed by Firmicutes (18.2 to 41.7%), Proteobacteria (0.6 to 11.3%), Patescibacteria (0.03 to 8.9%), Actinobacteria (0.2 to 2.9%) and Planctomycetes (0.01 to 3.1%) (Additional file 4: Figure S3A). The relative abundances of five phyla, Proteobacteria, Spirochaetes, Kiritimatiellaeota, Actinobacteria, and Patescibacteria, were significantly different among age groups based on DESeq2 analysis (P < 0.05, Additional file 2: Table S3). At the genus level, 64 genera were considered detectable using the same cut-offs mentioned above (Additional file 5: Figure S4A). Among them, Prevotella 1 (13.1 ± 0.1%) and Rikenellaceae RC9 gut group (10.8 ± 0.1%) were predominant regardless of age. Further differential abundance analysis using “DESq2” [30] identified 38 genera that were significantly different among age groups (Additional file 2: Table S3).

Changes in the rumen archaeal composition with yak growth

Taxonomic analysis revealed that most ESVs, from 94.4 to 99.8%, belonged to the phylum Euryarchaeota (Additional file 4: Figure S3B). At the species level, 35 species were identified, and 10 of them were considered detectable (relative abundance > 0.1% and present in more than half of the total animal populations per age group, Additional file 5: Figure S4B). Among these species, the Methanobrevibacter ruminantium clade (37.8 ± 10.2%), the Methanobrevibacter gottschalkii clade (34.8 ± 0.1%), Methanobacterium alkaliphilum (10.2 ± 0.1%), and Group 10 sp. belonged to the family Methanomassiliicoccaceae (5.9 ± 0.04%); Methanosphaera sp. ISO3-F5 (5.0 ± 0.02%) was the top five most abundant species, accounting for 93.7 ± 0.1% of the total number of sequences. Differential abundance analysis using DESeq2 showed that the relative abundances of the Methanobrevibacter ruminantium clade (55.8 ± 0.2%), Group 10 sp. (1.2 ± 0.01%) and Group 12 sp. ISO4 − H5 (0.1 ± 0.0%) were significantly different among age groups (P < 0.05, Additional file 2: Table S3). Specifically, ESVs belonging to the species of the Methanobrevibacter ruminantium clade were predominant at 1 month of age, while ESVs belonging to species Group 10 sp. belonged to the family Methanomassiliicoccaceae, and Methanosphaera sp. ISO3-F5 was less abundant at 1 year and 5 years of age, respectively.

Changes in the rumen fungal composition with yak growth

The fungal phylum Neocallimastigomycota (99.9 ± 0.3%, Additional file 4: Figure S3C) was the most dominant irrespective of age. In total, 7 fungal taxa were identified at the genus level, and 6 were considered detectable according to the selected criteria mentioned above (Additional file 5: Figure S4C). Among them, the top three genera were Caecomyces (35.3 ± 0.3%), Orpinomyces (29.1 ± 0.3%) and Neocallimastix (4.8 ± 0.1%). No ESVs were differentially abundant among age groups.

Changes in the rumen protozoal composition with yak growth

Four protozoal phyla were identified, of which the most abundant phylum was SAR, with a relative abundance of 99.4% (± 0.02) (Additional file 4: Figure S3D). In this phylum, the most abundant ESVs were classified as subclass Trichostomatia (40.4 ± 0.3%). At the genus level, a total of 11 genera were identified, and 7 genera were detectable, with Entodinium (19.2 ± 0.2%), Dasytricha (17.2 ± 0.2%) and Diplodinium (9.2 ± 0.2%) being the most abundant (Additional file 5: Figure S4D). There were no ESVs that differed significantly among age groups.

Maturation patterns differed among different microbial groups

The random forest regression model [31] was used to identify the age-discriminatory taxa of the rumen microbiota. In total, 15 bacterial (Fig. 4a), 9 archaeal (Fig. 4d), 2 protozoal (Fig. 5d) and 2 fungal (Fig. 5a) age-discriminatory taxa were identified according to the feature importance scores and 10-fold cross-validation (Additional file 2: Table S4). Specifically, the bacterial genus Olsenella (0.8 ± 0.01%) contributed the most to discriminating yak rumen microbiota differences according to age (Fig. 4a), while for the other three microbial groups, Group 10 sp. (6.1 ± 0.1%, archaea), Orpinomyces (22.9 ± 0.3%, fungi) and Dasytricha (17.2 ± 0.2%, protozoa) contributed the most to discriminating yak rumen microbiota differences according to age (Fig. 4d, Fig. 5a and d).

Fig. 4
figure 4

Bacterial and archaeal taxonomic biomarkers for defining rumen microbiota maturation in grazing yaks. a, d Heatmap of the mean relative abundance of the age-discriminatory taxa against time of sampling for each age group. b, e Age-discriminatory taxa were identified by applying random forest regression of the relative abundances of taxa in rumen fluid against chronological age of grazing yaks. Importance was determined based on the percentage increase in the mean-squared error of the predicted microbiota age when the relative abundance values of each taxon were randomly permuted (mean importance ± sem., n = 100 replicates). The inset shows 10-fold cross-validation error as a function of the number of input taxa used to regress against the age of yaks in the training set, in order of variable importance. c, f Microbiota age predictions plotted against chronological age. The curve is a smoothed spline fit between microbiota age and chronological age

Fig. 5
figure 5

Fungal and protozoal taxonomic biomarkers for defining rumen microbiota maturation in grazing yaks. a, d Age-discriminatory taxa were identified by applying random forest regression of the relative abundances of taxa in rumen fluid against chronological age of grazing yaks. Importance was determined based on the percentage increase in the mean-squared error of the predicted microbiota age when the relative abundance values of each taxon were randomly permuted (mean importance ± sem., n = 100 replicates). b, e The 10-fold cross-validation error as a function of the number of input taxa used to regress against the age of yaks in the training set, in order of variable importance. c, f Microbiota age predictions plotted against chronological age. The curve is a smoothed spline fit between microbiota age and chronological age

The regression model explained 70–72%, 45–54%, 23–25% and 48–50% of the variance in rumen bacteria, archaea, fungi and protozoa, respectively, based on the above age-discriminatory taxa. There was little improvement in the predictive performance (estimated by 10-fold cross validation) when any additional taxa beyond the age-discriminatory taxa were added for all four microbial groups (Fig. 4b and e; Fig. 5b and e). To quantitatively identify when the rumen microbiota mature, a smoothing spline function was fit between microbiota age and the corresponding chronological age of the animals. As a result, we found that the rumen bacteria reached full maturation between 5 and 8 years of age and archaea reached full maturation at approximately 5 years of age as the smoothed spline fit curve was saturated during this period (Fig. 4c and f). Similar to bacteria, fungi and protozoa reached full maturation between 5 and 8 years of age (Fig. 5c and f).

Intra- and inter-group microbial interactions in the rumen of grazing yaks

We first employed a co-occurrence network based on correlation relationships and P-values adjusted with FDR (false discovery rate) to explore the intragroup interactions of core taxa identified in all age groups (> 50% of the population of each age group). For bacteria, there were 17 nodes and 38 edges, with most nodes belonging to the phyla Bacteroidetes and Firmicutes (Fig. 6a). The topological properties were calculated to describe the complexity of the network (Additional file 2: Table S5). The top three genera identified as keystone taxa were Prevotella 1, Ruminococcaceae NK4A214 group and Ruminococcus 1 based on betweenness centrality scores (Fig. 6a), which measures the number of shortest paths going through a given node, as a proxy for the location of this node in relation to other nodes [32]. For archaea, the resulting co-occurrence network consisted of 5 nodes and 12 edges (Fig. 6b), and the Methanobrevibacter gottschalkii clade was identified as a keystone taxon based on its betweenness centrality score in the network (Additional file 2: Table S5). For fungi and protozoa, the network was less complicated. In the fungal network, Caecomyces was negatively correlated with the rest of the nodes (Neocallimastix and Neocallimastigaceae) that had a positive correlation (Fig. 6c). In the protozoal network, Dasytricha had a positive correlation with the other nodes (Entodinium and Trichostomatia) that had a positive correlation (Fig. 6d).

Fig. 6
figure 6

Co-occurring network analysis of intra-interactions of core taxa in each rumen microbial kingdom. The association patterns of bacteria (a), archaea (b), fungi (c) and protozoa (d) from birth to adulthood. The size of each node is proportional to the relative abundance. The lines in red and blue denote positive and negative correlations, respectively. The nodes are referred to as keystone species and are highlighted in light purple, and the taxonomic names of keystone species are indicated in the networks. The taxonomic information of the nodes is shown at the bottom of the figures

Next, we performed a meta-community co-occurrence network to explore the inter-group interactions of core taxa (please see methods). This generated a meta-community network with 46 links from 24 nodes (Additional file 6: Figure S5), including 17 bacterial nodes, 3 archaeal nodes, 1 fungal node, and 3 protozoal nodes. The relative importance of individual nodes within the network was computed via the topological features (Additional file 2: Table S6). Based on betweenness centrality scores, the top five taxa regarded as keystone species (Additional file 6: Figure S5) were Christensenellaceae R-7 group (3.9 ± 0.04%), Prevotella 1 (13.3 ± 0.1%), Trichostomatia (40.4 ± 0.3%), Ruminococcaceae UCG-014 (0.5 ± 0.02%) and Lachnospiraceae (1.2 ± 0.01%).

The intra- and inter-group interactions of core taxa under each age group (age-specific taxa) were further explored using co-occurrence network analysis. We first investigated the intra-interactions of age-specific taxa and found that the network complexity (as determined by the clustering coefficient and average degree scores) in rumen archaea, fungi and protozoa varied greatly with age, while that in bacteria changed at the first growth stages and remained stable after 3 years of age (Table 2). In addition, we identified 5 and 3 keystone species for bacteria and archaea, respectively, based on the betweenness centrality score. For bacteria, the keystone species varied with age, and only a few taxa were shared by some age groups. For example, Bacteroidales RF16 group uncultured bacterium (shared in 3 months, 6 months, 2 years and 10 years of age) and Streptococcus (shared in month 1, month 3 and 5 years of age) (Additional file 2: Table S7). For archaea, although the keystone species (top 3) differed among each age group, some taxa were shared by multiple age groups, including the Methanobrevibacter ruminantium clade and Methanobrevibacter gottschalkii clade (Additional file 2: Table S7). Similar to bacteria and archaea, the keystone species of fungi and protozoa changed with age, and some were shared by many age groups (Additional file 2: Table S7). Specifically, the fungi Caecomyces, Anaeromyces, Orpinomyces, and Neocallimastigaceae were shared by at least three age groups, while the protozoa Entodinium, Trichostomatia and Trichostomatia uncultured were shared by at least four age groups (Additional file 2: Table S7). Next, we examined the inter-kingdom interactions between age-specific taxa and found that the complexity of networks fluctuated greatly at the early growth stages (1 month to 2 years); the majority of keystone species (top 5) were from bacteria in all age groups, and unique keystone species were found in each age group (Table 3).

Table 2 Global topological properties of co-occurring networks of intra-interactions in each microbial kingdom
Table 3 Topological properties and keystone species of co-occurring networks of inter-kingdom interactions

Discussion

This study assessed the microbial groups (bacteria, archaea, fungi and protozoa) in the rumen and identified the age-discriminatory taxa and full maturation of rumen microbiota, keystone species of dynamic intra- and inter- microbial group interactions of natural grazing yaks at varied growth stages from birth to 12 years of age.

First, this study determined that rumen prokaryotes (bacteria and archaea) and eukaryotes (fungi and protozoa) colonized the rumen of natural grazing yaks at different life stages and matured differently. Among the four groups, the bacterial community colonized the rumen before day 7, which is similar to studies on dairy calves. For example, previous studies revealed that bacteria were detected in the rumen of dairy calves as early as birth and for the first week (1–7 days) of age [21, 33, 34], suggesting that rumen bacterial colonization occurs before the intake of solid food. The predominant bacterial phylum, Bacteroidetes, was found across all age groups, and the relative abundance of this phylum reached a maximum at 6 months of age, which is similar to the changes in the rumen of dairy cows from birth to 2 years of age [21]. Prevotella 1 was the dominant genus (21.7%), and its abundance remained stable after 3 years of age (fed solely with grass). At the early growth stage (days 7 and 14, the main diet was colostrum), the relative abundance of Bacteroides was higher than that of Prevotella 1 (0.53% vs 0.46 and 4.1% vs 1.9%, respectively). A previous study reported that the genus Prevotella (72%) was prominent in older animals (6 months and 2 years of age, fed solely with concentrate), while the genus Bacteroides dominated in the rumen of newborns (1–3 days of age, fed solely with colostrum) [21]. In addition, as the diet changed from milk to concentrate, the relative abundance of the genus Bacteroides in the rumen of dairy calves decreased from 16.9 to 7.1%, and the relative abundance of the genus Prevotella increased to 41.5% with increasing intake of concentrate [20]. However, in the present study, the relative abundance of the genus Prevotella 1 varied after grazing solely on natural grass (between 2 and 3 years of age), suggesting that the external living environment of grazing yaks is more variable than that of intensively farmed dairy cows. The dominance of the genus Prevotella 1 has also been reported in the rumen of 4-year-old captive yaks fed different doses of slow-release urea [35] as well as yaks (45 ± 5 months) under different feeding regimes [15]. This finding suggests that Prevotella 1 may be a conserved ‘core microbiota’ member in the rumen of grazing yaks. Previous studies have determined that members of this genus are involved in hemicellulose, protein [36] and starch degradation [37] but are not regarded as highly cellulolytic bacteria [38]. Future studies are needed to identify whether and how this predominant taxon functions in the rumen of natural grazing yaks raised in the harsh Qinghai-Tibetan Plateau environment.

Archaea were only detected after 14 days of age, which is different from dairy cows. A few studies have reported that archaea colonize the rumen of dairy calves during the first week of life [33, 38, 39], and one study even reported colonization from birth using a qPCR-based approach [40]. These findings suggest that the colonization of archaea in the rumen of natural grazing yaks differs from that of conventionally reared dairy cattle. At the species level, the Methanobrevibacter ruminantium clade and Methanobrevibacter gottschalkii clade were dominant in the rumen of yaks irrespective of age, suggesting that methanogenesis potentially occurs during the early growth stages as these methanogens are responsible for methane production [39]. In the present study, the relative abundance of the Methanobrevibacter gottschalkii clade peaked at 6 months of age and then gradually decreased with age, whereas that of the Methanobrevibacter ruminantium clade reached a minimum at 6 months of age and then increased with age, suggesting that these two species may compete in the rumen of grazing yaks during the early growth stages. In addition, these two species were dominant out of 32 ruminant species [41], indicating that these two species may be core methanogens in the rumen of all ruminants.

Fungi were only detected after 1 month of age, which is different from cattle. A previous study reported that fungi were detected in the rumen of dairy calves [33] and lambs [40] at 7 days of age, which is earlier than what we found in yaks. It is noticeable that we only collected rumen samples at 7 days of age, which may not be early enough to determine the exact time when fungi started to inhabit the rumen of grazing yaks. Future studies using the same method to analyze rumen samples collected at birth and confirm the time of initial colonization fungi in the rumen of grazing yaks are necessary. The phylum Neocallimastigomycota was predominant in the rumen of grazing yaks irrespective of age, and this result is consistent with findings of the rumen of dairy calves (from 7 days to 63 days) [33] and dairy cows [24]. Similar to a previous finding in dairy calves [33], unidentified Neocallimastigaceae, Caecomyces, and Orpinomyces were dominant in the rumen of grazing yaks irrespective of growth stage, indicating that these taxa may play important roles in rumen development during the early growth stages of ruminants.

Similar to fungi, protozoa were detected in the rumen of yaks after 1 month of age. A previous study detected protozoa in the rumen of lambs at approximately 21 days of age [42], which is earlier than that of yaks. In the present study, grazing yaks drank water from rivers or lakes instead of water troughs, which may preclude the colonization of rumen protozoa, as drinking water has been identified as a main source for rumen protozoa colonization [43]. Furthermore, rumen samples of newborn animals were collected in the winter, and rumen protozoa colonization may also be affected, as environmental temperature has been reported to affect the population and diversity of protozoa from soil and water [44, 45]. At the genus level, Entodinium and Dasytricha were predominant in the rumen of grazing yaks, which is consistent with previous findings in Yellow [46] and Hanwoo [47] cattle. In addition, a recent study found that the abundance of Dasytricha spp. increased, whereas the abundance of Entodinium decreased with increasing dietary fiber content [48]. In the present study, the genera Isotricha and Dasytricha were detected after 3 months of age when the diet of grazing yaks naturally changed from mother milk to native grass. Isotricha and Dasytricha play very important roles in utilizing soluble sugars and controlling the rate of carbohydrate fermentation [49]. This observation indicates that the selection of protozoa in the rumen of grazing yaks may be driven by dietary fiber contents.

In addition to identifying the microbiota that colonize rumen, it is also important to know when such microbiota becomes established. This provides information for the potential window to manipulate microbiota, as it is more difficult to alter fully established microbiota [50]. We characterized the rumen microbiota in grazing yaks fully matured between 5 and 8 years of age. Specifically, rumen archaea became fully maturated at approximately 5 years old, while the other microbial groups (bacteria, fungi and protozoa) became fully maturated between 5 and 8 years of age, indicating that the time for rumen microbiota to become fully maturated differs among these four microbial groups. We also observed that this maturation pattern differs based on the changes in alpha diversity. Although the alpha diversity indices became stable at 2 years of age, this does not mean that the rumen microbiota become fully mature as the function of some taxa may still change with host physiological development [51]. In addition, it has been reported that malnutrition delays gut microbiota maturation [52]. The grazing yaks in this study are malnourished conditions, as in the long cold season, the Qinghai-Tibetan Plateau provides limited food resources. Thus, rumen microbial maturation in grazing yaks may be delayed. Yak is traditionally considered an adult between 4 and 7 years old [53, 54], which is similar to the predicted maturation time of rumen microbiota in this study. These findings suggest that the stability of rumen microbiota diversity may not be an indicator of microbiota maturation in the rumen of grazing yaks. Moreover, the age-discriminatory taxa within each microbial group were identified, and these taxa have the power to discriminate the age of the rumen microbiota. Among these, Olsenella (bacteria), Group 10 sp., belonging to the family Methanomassiliicoccaceae (archaea), Orpinomyces (fungi) and Dasytricha (protozoa) contributed the most to discriminating the age of the rumen microbiota, as their abundance tended to be stable across growth stages, and they were present in most age groups [55]. Therefore, future studies are necessary to identify to what extent the age-discriminatory taxa are indicative of the normal development of the rumen microbiota as well as to determine when biomarkers of rumen microbiota are established.

Similar to previous reports in dairy calves from birth to adulthood [21] and pre-weaning dairy calves [33], we also found that inter-individual animal variation in the bacteria and archaea of rumen decreased with age in grazing yaks. This indicates that high fluctuations of rumen bacterial and archaeal communities may be prevalent in the neonates of ruminant species. Based on PCoA analysis, the community of rumen bacteria and archaea was clearly separated between the young age groups (7 days, 14 days and 1 month of age) and the older age groups (2–12 years of age), which also was corresponded to two feeding regimes: mother’s milk and natural grass. This suggests that both age and diet type contribute to the colonization of bacteria and archaea in the rumen of grazing yaks, as previously reported in dairy calves [25]. However, previous studies on dairy calves [33] and goats [56] have found that rumen archaeal communities are less sensitive to age change, and further studies should be performed to determine which factor contributes more to determining the development of archaea in the rumen of grazing yaks. No clear separation was observed in the protozoal and fungal profiles. The rumen fungal community was also reported to be resistant to changes in age between primiparous and multiparous cows [24], and the rumen protozoal diversity changes little throughout life, although the relative abundance of protozoal species in rumen fluctuates with diet [28]. These results suggest that fungal and protozoal communities may be less sensitive to changes in age and diet type than other microbial groups.

Rumen microbiota work synergistically to perform various metabolic functions that ferment fibrous plant materials. Based on the intra-interactions of rumen bacteria, archaea, fungi and protozoa across the lifetime of yaks, we identified different interactions within each microbial group. In the bacterial network, few strong negative correlations were found; however, few strong positive correlations were found within bacteria in the rumen of other ruminant species [41]. Previous studies have reported that the dynamic interactions of rumen microbiota are specific to diet in dairy cows [57] and can be affected by the host [41]. Furthermore, these interactions differed in the rumen of multiparous cows compared to those in primiparous cows [24], suggesting that the interactions of rumen bacteria could also be affected by age and host physiological state. For the archaeal network, the Methanobrevibacter gottschalkii clade and Methanobrevibacter ruminantium clade had a negative relationship, which is similar to previous findings in the rumen of adult sheep, cattle and deer fed different diets [58] indicating that these two species compete for H2 in the rumen [59]. For the fungal network, the genus Caecomyces and the unclassified family Neocallimastigaceae had a strong negative correlation, which is different from what has been reported in the rumen of adult cows fed different diets (80% forage and 20% concentrate) [24], suggesting that the interaction between rumen fungi may be affected by host and diet. A significant negative correlation was observed between subclass Trichostomatia and Entodinium. Members of the subclass Trichostomatia have greater endoglucanase and xylanase activity, while members of the genus Entodinium have only weak endoglucanase and xylanase activity [60], suggesting that the lack of co-occurrence may be due to the exploitation of different opportunities. A previous study revealed that the complexity of the relationships of rumen microbes shifts from nongrazing to grazing and that these microbes work together to adapt to the dietary shift in sheep [61], which is similar to our results. Based on these findings, it is suggested that specific associations within each microbial group may exist across diets and hosts, with the ability to adapt to the specific life environment and potentially utilize the available substrates.

For the identified inter-interactions among bacteria, archaea, fungi and protozoa under different growth stages; fungi were (Caecomyces) positively correlated with bacteria (Prevotellaceae), reflecting a synergistic relationship, and the degradation of plant fiber by fungi appears to facilitate more rapid breakdown of forage by fibrolytic bacteria [57]. The lack of strong associations between archaea and protozoa in this study is similar to previous observations in the rumen of multiparous Nordic Red dairy cows [57] and many different ruminant species from different locations [41], suggesting that the associative patterns between archaea and protozoa are less specific and more random in vivo. In addition, no strong associations were found between bacteria and archaea except for the positive association between Methanosphaera sp. ISO3-F5 and Ruminococcaceae NK4A214 group. This association is unexpected given that rumen bacteria breakdown complex compounds and produce substrates that methanogens use for growth, mainly, hydrogen and methyl-containing compounds [62]. This indicates that the interactions between bacteria and archaea in the rumen are probably host-specific and warrant future studies.

Last, co-occurrence analysis revealed distinct associative patterns of age-specific taxa in different microbial groups. Specifically, the complexity of the bacterial network of the intra-interactions of age-specific taxa fluctuated greatly at early stages of growth (7 days to 2 years old) but remained stable afterwards, and unique keystone species were identified in almost every age group, indicating that the dynamic association patterns between rumen microbial groups changed with growth stages or diet. Regarding the inter-interactions of age-specific taxa, the change in network complexity is similar to that in the intra-interaction, and the majority of keystone species under different growth stages belong to bacteria, indicating that bacteria may play a more central role in rumen biological networks. In addition, the keystone species identified under each age included age specific and core taxa, suggesting they both play important roles in the rumen of yaks at different growth stages.

Conclusions

This study comprehensively explored the rumen microbiota (including bacteria, archaea, protozoa and fungi) in grazing yaks for the first time. Our results revealed that the alpha diversity of rumen microbiota increases with age and generally remains stable after 2 years of age. The rumen microbial communities underwent multiple changes with growth stage, where bacteria and archaea were more sensitive to changes in age compared to other microbial groups. In addition, bacteria and archaea were observed before 7 days, while protozoa and fungi were detected at approximately 1 month of age, suggesting that rumen prokaryotes generally appeared earlier than eukaryotes in the rumen of grazing yaks. Furthermore, four rumen microbial groups had their own maturation trajectory over the lifetime of yaks, where rumen archaea fully maturated at approximately 5 years of age and the other microbial groups maturated between 5 and 8 years of age. Distinct associative patterns among ruminal microbial groups were observed in each age group, and the dynamic intra- and inter-interactions among four microbial groups changed across the lifetime of yaks. Rumen bacteria play a central role in the rumen biological networks, as most keystone species in all age groups belonged to bacteria (including Christensenellaceae R-7 group, Prevotella 1, Ruminococcaceae UCG-014 and Lachnospiraceae). The interactive analysis of this study provides novel insight into elucidating the dynamic intra- and inter- interactions of rumen microbial groups in grazing yaks across a lifetime, providing a solid basis for the manipulation of rumen at different growth stages to improve performance in the harsh Qinghai-Tibetan Plateau ecosystem. One of the limitations of the current study was the lack of rumen fermentation parameters (such as pH, VFA and ammonia) and quantification of microbial abundance by qPCR due to the challenges of collecting enough rumen samples from yaks, especially from the young age groups (from 7 days to 3 months of age). These data, if available, could provide evidence from a functional perspective to reflect rumen development. Regardless, the findings of this study provide fundamental knowledge of the rumen microbiota of grazing yaks during different growth stages. Future studies on the function of the rumen microbiota and its relationship with host development should be explored to identify whether the rumen microbiota and its maturation and interaction can contribute to host growth and productivity.

Methods

Animals and sampling

The yaks enrolled in this study were weaned naturally (approximately 12–18 months old) and grazed naturally together year round (without concentrate supplementation) on native pastures from Wushaoling of Tianzhu Autonomous County, Gansu Province (37°12.4′N, 102°51.7′E; altitude of 3154 m). Specifically, all the animals used in this study were from the same herd and they all grazed together in an alpine meadow on the Qinghai-Tibetan Plateau, and they drank water from the local river or the snow meltwater. In addition, animals were grazing on the grassland except at sampling time (animals were rounded up at night before the day of sampling). Rumen fluid samples were collected from yaks ranging from 7 days to 12 years old, and the experimental design is shown in Fig. 1. Briefly, rumen fluid samples were repeatedly sampled from 8 yaks from 7 days to 1 year old. Rumen fluid samples from another 36 yaks aged 2 years (n = 6), 3 years (n = 6), 5 years (n = 6), 8 years (n = 6), 10 years (n = 6) and 12 years (n = 6) were collected in 1 day within 3 h. Specifically, rumen samples from yaks aged 6 months to 12 years old were collected via a perforated stainless-steel stomach tube connected to a suction pump prior to morning grazing. The procedure for sample collection from calves (7 days to 4 months, isolated from their mothers on the day before sample collection) was similar except the stainless-steel stomach tube was replaced by a plastic flexible stomach tube. The first 5 ml (from 7 days to 4 months of age) and 20 ml (from 6 months to 12 years of age) rumen fluid were discarded to avoid any contamination with saliva. The rumen fluid (approximately 20–100 ml for yaks aged 7 days to 4 months, 250 ml for yaks aged 6 months to 12 years) was collected and divided into aliquots in 10 ml polypropylene tubes. The rumen samples were immediately deep frozen using liquid nitrogen, transported to the laboratory and stored at − 80 °C prior to DNA extraction.

DNA extraction

Total genomic DNA was extracted from freeze-dried samples using a PowerSoil DNA Isolation Kit (MoBIO Laboratories, Inc., Carlsbad, CA, USA) according to the manufacturer’s instructions. The concentration and quality of DNA were assessed using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Scoresby, Australia) and agarose gel (1.0%) electrophoresis, respectively. The isolated DNA was stored at − 20 °C until downstream analysis.

Amplification of target genes for rumen microbiota profiling

For rumen microbial profiling, primers Ba9f (GAGTTTGATCMTGGCTCAG) and Ba515Rmod1 (CCGCGGCKGCTGGCAC) targeting the bacterial partial 16S rRNA gene, primers Ar915aF (AGGAATTGGCGGGGGAGCAC) and Ar1386R (GCGGTGTGTGCAAGGAGC) targeting the archaeal partial 16S rRNA gene [41], primers Reg841F (GACTAGGGATTGGAGTGG and Reg1302R (AATTGCAAAGATCTATCCC) targeting the ciliate protozoal 18S rRNA gene [58], and primers MN100F (TCCTACCCTTTGTGAATTTG) and MNGM2 (CTGCGTTCTTCATCGTTGCG) targeting the internal transcribed spacer region of anaerobic fungi [63] were used to generate amplicons for each rumen microbial group. The PCR amplification products were verified using agarose gel (2%) electrophoresis and purified with a Qiagen Gel Extraction Kit (Qiagen, Germany). Mixed samples were prepared by pooling equal amounts of PCR amplicons from each sample and then sequenced on an Illumina PE MiSeq 300 platform to generate 300-bp paired end reads.

Sequencing data analysis

Among the samples, only 80, 75, 65 and 71 samples were successfully amplified for bacteria, archaea, fungi and protozoa, respectively. Detailed sample numbers for each age group are listed in Table 1. The raw sequence data were assigned to each sample according to the corresponding barcode and were processed using QIIME2 [64]. The DADA2 algorithm [65] as a QIIME2 plugin was used to preprocess the demultiplexed paired-end sequence reads, including quality filtering, denoising, joining paired ends, and removing chimeric sequences. Sequences were clustered into exact sequence variants (ESVs). Next, taxonomy was assigned to ESVs via the “qiime feature-classifier” command using the “classify-sklearn” option [66] against the SILVA 132 database for bacteria and protozoa, the RIM-DB database for archaea, and the UNITE database for fungi. For diversity analyses, sequences were aligned using Mafft [67], and noninformative positions in the alignment were filtered with the ‘mask’ command. Next, a midpoint-rooted phylogenetic tree followed by a phylogenetic tree was performed using the FastTree plugin [68]. Finally, each sample was rarefied to 7158 (bacteria), 3794 (archaea), 4761 (fungi) and 1818 (protozoa) sequences respectively prior to calculating alpha- and beta- diversity metrics; including Faith’s phylogenetic diversity (PD), Chao1 index, and observed OTUs for alpha diversity, and unweighted UniFrac and weighted UniFrac distance for beta diversity.

Definition of rumen microbiota maturation in grazing yaks using random forests

To identify the characteristics of ESVs of rumen microbiota maturation, random forest regression [31] was used, and the frequency of ESVs in the temporal profiles of rumen microbiota against chronological age with default parameter was regressed in R studio (3.5.3). The model was randomly rebuilt 100 times, and the feature importance scores were averaged across the 100 models [69]. To estimate the minimal number of taxa that generated the lowest cross-validation error needed to predict the final model, the “rfcv” function implemented in the “randomForest” package was applied 100 times. A sparse model was built on a subset of predictive variables that were determined according to their feature importance scores and the number of taxa to be included, which was used to predict chronological age, described as “microbiota age”. A smoothing spline function in R studio (3.5.3) was fit between microbiota age and the corresponding chronological age of the animals (at the time of rumen sample collection). When the curve reached a plateau, the microbiota reached full maturity [69].

Co-occurrence network analysis

ESVs with relative abundances > 0.01% were used for co-occurrence analysis to explore the co-occurrence patterns within and between bacteria, archaea, fungi and protozoa. Taxa that occurred in > 50% of the population from each age group and were present in all age groups were termed “core taxa” [70, 71]; these taxa were used to explore interactions of rumen microbes from neonates to adults. Meanwhile, the taxa that were only detected in at least 50% of the samples from each age group were regarded as “age-specific microbiota”; these taxa were applied to elucidate hub taxa in different age groups. The co-occurrence network was built based on the Spearman correlation matrix calculated with the WGCNA package [72], and P-values for multiple testing were adjusted using the Benjamini and Hochberg false discovery rate (FDR) controlling procedure [73]. The nodes in the network represent rumen taxa, whereas the edges (connections) correspond to the correlation between nodes. The network images were generated using Cytoscape 3.7.1 [74]. When exploring the interactions of core taxa, only correlations with an R-corr absolute value greater than 0.3 and adjusted P-value less than 0.05 were plotted. Regarding the threshold of age-specific taxa, to better reflect the entire microbial relationships, correlations with an R-corr absolute value ≥0.5 were used for further analysis without filtering P values. In addition to network topology parameters, betweenness centrality was used to measure the extent to which a node lied on the shortest path between other nodes in the network. Other topological features, such as the degree centrality, transitivity and closeness centrality, were also calculated using the igraph package to describe the complexity of this network. Nodes with the highest betweenness centrality scores were considered keystone species in co-occurrence networks [75].

Statistical analysis

The alpha diversity indices among different age groups were compared using the nonparametric Kruskal-Wallis testing method to determine if there were significant differences using RStudio (3.5.3). Principal coordinates analysis (PCoA) was performed using the qiime2R package [76] to visualize dissimilarity based on weighted and unweighted UniFrac distance metrics. A pairwise PERMANOVA test was carried out to test the significant differences in microbiome beta diversity among age groups, and default permutations were used to calculate the P-value in QIIME2 [77]. In addition, we employed DESeq2 [30] implemented in RStudio to investigate the differentially abundant taxa in different age groups, and this method could handle the uneven sample size during pairwise comparison between groups. The Mann-Whitney U test, which can handle uneven samples [78], was selected to test the alpha diversity indices to compare each group. DESeq2 uses differential expression statistical Wald tests and is adjusted by applying the Benjamini–Hochberg method to correct for multiple hypothesis testing. The FDR cutoff was set at 0.05. Most of the figures presented in this study were generated by RStudio using the corresponding packages unless otherwise noted, and d, m and y indicate day, month and year, respectively, in all tables and figures throughout the manuscript.

Availability of data and materials

Raw sequencing data have been deposited at the NCBI Sequence Read Archive (SRA) under BioProject ID PRJNA566442.

Abbreviations

ESVs:

Exact sequence variants

16S rRNA:

16 Svedberg ribosomal ribonucleic acid

18S rRNA:

18 Svedberg ribosomal ribonucleic acid

PCR:

Polymerase chain reaction

PCoA:

Principal coordinate analysis

FDR:

False discovery rate

References

  1. Gu Z, Zhao X, Li N, Wu C. Complete sequence of the yak (Bos grunniens) mitochondrial genome and its evolutionary relationship with other ruminants. Mol Phylogenet Evol. 2007;42:248–55.

    CAS  PubMed  Google Scholar 

  2. Qiu Q, Wang L, Wang K, Yang Y, Ma T, Wang Z, et al. Yak whole-genome resequencing reveals domestication signatures and prehistoric population expansions. Nat Commun. 2015;6:1–7.

    CAS  Google Scholar 

  3. Ge RL, Cai Q, Shen YY, San A, Ma L, Zhang Y, et al. Draft genome sequence of the Tibetan antelope. Nat Commun. 2013;4:1–7.

    CAS  Google Scholar 

  4. Dolt KS, Mishra MK, Karar J, Baig MA, Ahmed Z, Pasha MAQ. cDNA cloning, gene organization and variant specific expression of HIF-1α in high altitude yak (Bos grunniens). Gene. 2007;386:73–80.

    CAS  PubMed  Google Scholar 

  5. Shao B, Long R, Ding Y, Wang J, Ding L, Wang H. Morphological adaptations of yak (Bos grunniens) tongue to the foraging environment of the Qinghai-Tibetan plateau. J Anim Sci. 2010;88:2594–603.

    CAS  PubMed  Google Scholar 

  6. Long RJ, Ding LM, Shang ZH, Guo XH. The yak grazing system on the Qinghai-Tibetan plateau and its status. Rangel J. 2008;30:241–6.

    Google Scholar 

  7. Dong QM, Zhao XQ, Ma YS, Xu SX, Li QY. Live-weight gain, apparent digestibility, and economic benefits of yaks fed different diets during winter on the Tibetan plateau. Livest Sci. 2006;101:199–207.

    Google Scholar 

  8. Zhou JW, Liu H, Zhong CL, Degen AA, Yang G, Zhang Y, et al. Apparent digestibility, rumen fermentation, digestive enzymes and urinary purine derivatives in yaks and Qaidam cattle offered forage-concentrate diets differing in nitrogen concentration. Livest Sci. 2018;208:14–21.

    Google Scholar 

  9. Zhou JW, Zhong CL, Liu H, Degen AA, Titgemeyer EC, Ding LM, et al. Comparison of nitrogen utilization and urea kinetics between yaks (Bos grunniens) and indigenous cattle (Bos taurus). J Anim Sci. 2017;95:4600–12.

    CAS  PubMed  Google Scholar 

  10. Ishaq SL, Wright ADG. Insight into the bacterial gut microbiome of the north American moose (Alces alces). BMC Microbiol. 2012;12:212.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. An D, Dong X, Dong Z. Prokaryote diversity in the rumen of yak (Bos grunniens) and Jinnan cattle (Bos taurus) estimated by 16S rDNA homology analyses. Anaerobe. 2005;11:207–15.

    CAS  PubMed  Google Scholar 

  12. Huang XD, Tan HY, Long R, Liang JB, Wright ADG. Comparison of methanogen diversity of yak (Bos grunniens) and cattle (Bos taurus) from the Qinghai-Tibetan plateau, China. BMC Microbiol. 2012;12:237.

    PubMed  PubMed Central  Google Scholar 

  13. Guo W, Li Y, Wang L, Wang J, Xu Q, Yan T, et al. Evaluation of composition and individual variability of rumen microbiota in yaks by 16S rRNA high-throughput sequencing technology. Anaerobe. 2015;34:74–9.

    CAS  PubMed  Google Scholar 

  14. Liu C, Wu H, Liu S, Chai S, Meng Q, Zhou Z. Dynamic alterations in yak rumen bacteria community and metabolome characteristics in response to feed type. Front Microbiol. 2019;10:1–19.

    PubMed  PubMed Central  Google Scholar 

  15. Zhou Z, Fang L, Meng Q, Li S, Chai S, Liu S, et al. Assessment of ruminal bacterial and archaeal community structure in yak (Bos grunniens). Front Microbiol. 2017;8:1–10.

    Google Scholar 

  16. Xue D, Chen H, Zhao X, Xu S, Hu L, Xu T, et al. Rumen prokaryotic communities of ruminants under different feeding paradigms on the Qinghai-Tibetan plateau. Syst Appl Microbiol. 2017;40:227–36.

    PubMed  Google Scholar 

  17. Zhang Z, Xu D, Wang L, Hao J, Wang J, Zhou X, et al. Convergent evolution of rumen microbiomes in high-altitude mammals. Curr Biol. 2016;26:1873–9.

    CAS  PubMed  Google Scholar 

  18. Rabee AE, Forster RJ, Elekwachi CO, Kewan KZ, Sabra EA, Shawket SM, et al. Community structure and fibrolytic activities of anaerobic rumen fungi in dromedary camels. J Basic Microbiol. 2019;59:101–10.

    CAS  PubMed  Google Scholar 

  19. Belanche A, Yáñez-Ruiz DR, Detheridge AP, Griffith GW, Kingston-Smith AH, Newbold CJ. Maternal versus artificial rearing shapes the rumen microbiome having minor long-term physiological implications. Environ Microbiol. 2019;21:4360–77.

    PubMed  PubMed Central  Google Scholar 

  20. Rey M, Enjalbert F, Combes S, Cauquil L, Bouchez O, Monteils V. Establishment of ruminal bacterial community in dairy calves from birth to weaning is sequential. J Appl Microbiol. 2014;116:245–57.

    CAS  PubMed  Google Scholar 

  21. Jami E, Israel A, Kotser A, Mizrahi I. Exploring the bovine rumen bacterial community from birth to adulthood. ISME J. 2013;7:1069–79.

    PubMed  PubMed Central  Google Scholar 

  22. Dill-Mcfarland KA, Breaker JD, Suen G. Microbial succession in the gastrointestinal tract of dairy cows from 2 weeks to first lactation. Sci Rep. 2017;7:1–12.

    Google Scholar 

  23. Liu C, Meng Q, Chen Y, Xu M, Shen M, Gao R, et al. Role of age-related shifts in rumen bacteria and methanogens in methane production in cattle. Front Microbiol. 2017;8:1–14.

    Google Scholar 

  24. Kumar S, Indugu N, Vecchiarelli B, Pitta DW. Associative patterns among anaerobic fungi, methanogenic archaea, and bacterial communities in response to changes in diet and age in the rumen of dairy cows. Front Microbiol. 2015;6:1–10.

    Google Scholar 

  25. Dill-McFarland KA, Weimer PJ, Breaker JD, Suen G. Diet influences early microbiota development in dairy calves without Long-term impacts on Milk production. Appl Environ Microbiol. 2019;85:1–12.

    Google Scholar 

  26. Yáñez-Ruiz DR, Abecia L, Newbold CJ. Manipulating rumen microbiome and fermentation through interventions during early life: a review. Front Microbiol. 2015;6:1–12.

    Google Scholar 

  27. Abecia L, Martínez-Fernandez G, Martín-García AI, Ramos-Morales E, Yáñez-Ruiz DR, Waddams KE, et al. An antimethanogenic nutritional intervention in early life of ruminants modifies ruminal colonization by archaea. Archaea. 2014;2014:1–12.

  28. Huws SA, Creevey CJ, Oyama LB, Mizrahi I, Denman SE, Popova M, et al. Addressing global ruminant agricultural challenges through understanding the rumen microbiome: past, present, and future. Front Microbiol. 2018;9:1–33.

    Google Scholar 

  29. Wiener GH, Jianlin H, Ruijun L. Origins, domestication and distribution of yak. The yak. 2nd ed: RAP Publication; 2003.

  30. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.

    Google Scholar 

  31. Ruggles KV, Wang J, Volkova A, Contreras M, Noya-Alarcon O, Lander O, et al. Changes in the gut microbiota of urban subjects during an immersion in the traditional diet and lifestyle of a rainforest village. mSphere. 2018;3:1–8.

    Google Scholar 

  32. Ma B, Wang H, Dsouza M, Lou J, He Y, Dai Z, et al. Geographic patterns of co-occurrence network topological features for soil microbiota at continental scale in eastern China. ISME J. 2016;10:1891–901.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Dias J, Marcondes MI, Noronha MF, Resende RT, Machado FS, Mantovani HC, et al. Effect of pre-weaning diet on the ruminal archaeal, bacterial, and fungal communities of dairy calves. Front Microbiol. 2017;8:1553.

    PubMed  PubMed Central  Google Scholar 

  34. Malmuthuge N, Liang G, Guan LL. Regulation of rumen development in neonatal ruminants through microbial metagenomes and host transcriptomes. Genome Biol. 2019;20:1–6.

    CAS  Google Scholar 

  35. Yan XT, Yan BY, Ren QM, Dou JJ, Wang WW, Zhang JJ, et al. Effect of slow-release urea on the composition of ruminal bacteria and fungi communities in yak. Anim Feed Sci Technol. 2018;244:18–27.

    CAS  Google Scholar 

  36. Thoetkiattikul H, Mhuantong W, Laothanachareon T, Tangphatsornruang S, Pattarajinda V, Eurwilaichitr L, et al. Comparative analysis of microbial profiles in cow rumen fed with different dietary fiber by tagged 16S rRNA gene pyrosequencing. Curr Microbiol. 2013;67:130–7.

    CAS  PubMed  Google Scholar 

  37. Dias J, Marcondes MI, de Souza SM, da Mata e Silva BC, Noronha MF, Resende RT, et al. Bacterial community dynamics across the gastrointestinal tracts of dairy calves during preweaning development. Appl Environ Microbiol. 2018;84:e02675–17.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Anderson KL, Nagaraja TG, Morrill JL, Avery TB, Galitzer SJ, Boyer JE. Ruminal microbial development in conventionally or early-weaned calves. J Anim Sci. 1987;64:1215–26.

    CAS  PubMed  Google Scholar 

  39. Wang Z, Elekwachi CO, Jiao J, Wang M, Tang S, Zhou C, et al. Investigation and manipulation of metabolically active methanogen community composition during rumen development in black goats. Sci Rep. 2017;7:422.

    PubMed  PubMed Central  Google Scholar 

  40. Jiao J, Li X, Beauchemin KA, Tan Z, Tang S, Zhou C. Rumen development process in goats as affected by supplemental feeding v. grazing: age-related anatomic development, functional achievement and microbial colonisation. Br J Nutr. 2015;113:888–900.

    CAS  PubMed  Google Scholar 

  41. Henderson G, Cox F, Ganesh S, Jonker A, Young W, Janssen PH, 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:14567.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. Chaucheyras-Durand F, Ameilbonne A, Auffret P, Bernard M, Mialon MM, Dunière L, et al. Supplementation of live yeast based feed additive in early life promotes rumen microbial colonization and fibrolytic potential in lambs. Sci Rep. 2019;16:1–6.

    Google Scholar 

  43. Bird SH, Hegarty RS, Woodgate R. Modes of transmission of rumen protozoa between mature sheep. Anim Prod Sci. 2010;50:414–7.

    Google Scholar 

  44. Opperman MH, Wood M, Harris PJ. Changes in microbial populations following the application of cattle slurry to soil at two temperatures. Soil Biol Biochem. 1989;21:263–8.

    Google Scholar 

  45. Canals O, Serrano-Suárez A, Salvadó H, Méndez J, Cervero-Aragó S, De Porras VR, et al. Effect of chlorine and temperature on free-living protozoa in operational man-made water systems (cooling towers and hot sanitary water systems) in Catalonia. Environ Sci Pollut R. 2015;22:6610–8.

    CAS  Google Scholar 

  46. Leng J, Zhong X, Zhu RJ, Yang SL, Gou X, Mao HM. Assessment of protozoa in Yunnan yellow cattle rumen based on the 18S rRNA sequences. Mol Biol Rep. 2011;38:577–85.

    CAS  PubMed  Google Scholar 

  47. Shin EC, Cho KM, Lim WJ, Hong SY, An CL, Kim EJ, et al. Phylogenetic analysis of protozoa in the rumen contents of cow based on the 18S rDNA sequences. J Appl Microbiol. 2004;97:378–83.

    CAS  PubMed  Google Scholar 

  48. Mao SY, Huo WJ, Zhu WY. Microbiome-metabolome analysis reveals unhealthy alterations in the composition and metabolism of ruminal microbiota with increasing dietary grain in a goat model. Environ Microbiol. 2016;18:525–41.

    CAS  PubMed  Google Scholar 

  49. Wright ADG. Rumen protozoa. In: Rumen microbiology: from evolution to revolution; 2015.

    Google Scholar 

  50. Zhou M, Peng YJ, Chen Y, Klinger CM, Oba M, Liu JX, et al. Assessment of microbiome changes after rumen transfaunation: implications on improving feed efficiency in beef cattle. Microbiome. 2018;6:62.

    PubMed  PubMed Central  Google Scholar 

  51. Friedman N, Jami E, Mizrahi I. Compositional and functional dynamics of the bovine rumen methanogenic community across different developmental stages. Environ Microbiol. 2017;19:3365–73.

    CAS  PubMed  PubMed Central  Google Scholar 

  52. Subramanian S, Huq S, Yatsunenko T, Haque R, Mahfuz M, Alam MA, et al. Persistent gut microbiota immaturity in malnourished Bangladeshi children. Nature. 2014;510:417–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Long RJ, Zhang DG, Wang X, Hu ZZ, Dong SK. Effect of strategic feed supplementation on productive and reproductive performance in yak cows. Prev Vet Med. 1999;38:195–206.

    CAS  PubMed  Google Scholar 

  54. Zi XD, Zhong GH, Wen YL, Zhong JC, Liu CL, Ni YA, et al. Growth performance, carcass composition and meat quality of Jiulong-yak (Bos grunniens). Asian-Austr J Anim Sci. 2004;17:410–4.

    Google Scholar 

  55. Oh J, Byrd AL, Park M, Kong HH, Segre JA. Temporal stability of the human skin microbiome. Cell. 2016;165:854–66.

    CAS  PubMed  PubMed Central  Google Scholar 

  56. Wang L, Xu Q, Kong F, Yang Y, Wu D, Mishra S, et al. Exploring the goat rumen microbiome from seven days to two years. PLoS One. 2016;11:e0154354.

    PubMed  PubMed Central  Google Scholar 

  57. Tapio I, Fischer D, Blasco L, Tapio M, Wallace RJ, Bayat AR, et al. Taxon abundance, diversity, co-occurrence and network analysis of the ruminal microbiota in response to dietary changes in dairy cows. PLoS One. 2017;12:e0180260.

    PubMed  PubMed Central  Google Scholar 

  58. Kittelmann S, Seedorf H, Walters WA, Clemente JC, Knight R, Gordon JI, et al. Simultaneous amplicon sequencing to explore co-occurrence patterns of bacterial, archaeal and eukaryotic microorganisms in rumen microbial communities. PLoS One. 2013;8:e47879.

    CAS  PubMed  PubMed Central  Google Scholar 

  59. Janssen PH, Kirs M. Structure of the archaeal community of the rumen. Appl Environ Microbiol. 2008;74:3619–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  60. Newbold CJ, De la Fuente G, Belanche A, Ramos-Morales E, McEwan NR. The role of ciliate protozoa in the rumen. Front Microbiol. 2015;6:1–14.

    Google Scholar 

  61. Belanche A, Kingston-Smith AH, Griffith GW, Newbold CJ. A multi-kingdom study reveals the plasticity of the rumen microbiota in response to a shift from non-grazing to grazing diets in sheep. Front Microbiol. 2019;10:122.

    PubMed  PubMed Central  Google Scholar 

  62. Danielsson R, Dicksved J, Sun L, Gonda H, Müller B, Schnürer A, et al. Methane production in dairy cows correlates with rumen methanogenic and bacterial community structure. Front Microbiol. 2017;8:1–15.

    Google Scholar 

  63. Li Z, Wright ADG, Si H, Wang X, Qian W, Zhang Z, et al. Changes in the rumen microbiome and metabolites reveal the effect of host genetics on hybrid crosses. Environ Microbiol Rep. 2016;8:1016–23.

    CAS  PubMed  Google Scholar 

  64. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.

    CAS  PubMed  PubMed Central  Google Scholar 

  66. Hakim J, Schram J, Galloway A, Morrow C, Crowley M, Watts S, et al. The Purple Sea urchin Strongylocentrotus purpuratus demonstrates a compartmentalization of gut bacterial microbiota, predictive functional attributes, and taxonomic co-occurrence. Microorganisms. 2019;7:35.

    CAS  PubMed Central  Google Scholar 

  67. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Price MN, Dehal PS, Arkin AP. FastTree 2 - approximately maximum-likelihood trees for large alignments. PLoS One. 2010;5:e9490.

    PubMed  PubMed Central  Google Scholar 

  69. Gao P, Ma C, Sun Z, Wang L, Huang S, Su X, et al. Feed-additive probiotics accelerate yet antibiotics delay intestinal microbiota maturation in broiler chicken. Microbiome. 2017;5:91.

    PubMed  PubMed Central  Google Scholar 

  70. Cernava T, Erlacher A, Soh J, Sensen CW, Grube M, Berg G. Enterobacteriaceae dominate the core microbiome and contribute to the resistome of arugula (Eruca sativa mill.). Microbiome. 2019;7:13.

    PubMed  PubMed Central  Google Scholar 

  71. Pérez-Losada M, Alamri L, Crandall KA, Freishtat RJ. Nasopharyngeal microbiome diversity changes over time in children with asthma. PLoS One. 2017;12:e0170543.

    PubMed  PubMed Central  Google Scholar 

  72. Manuscript A, Structures T. Fast R functions for robust correlations and hierarchical clustering. NIH Public Access, JStat Softw. 2009;6:247–53.

    Google Scholar 

  73. Krieger AM, Yekutieli D. Adaptive linear step-up Proceudres that control the false discovery rate. Biometrika. 2006;93:491–507.

    Google Scholar 

  74. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    CAS  PubMed  PubMed Central  Google Scholar 

  75. Jiao S, Liu Z, Lin Y, Yang J, Chen W, Wei G. Bacterial communities in oil contaminated soils: biogeography and co-occurrence patterns. Soil Biol Biochem. 2016;98:64–73.

    CAS  Google Scholar 

  76. Kunz IGZ, Reed KJ, Metcalf JL, Hassel DM, Coleman RJ, Hess TM, et al. Equine fecal microbiota changes associated with anthelmintic administration. J Equine Vet Sci. 2019;77:98–106.

    PubMed  Google Scholar 

  77. Anderson MJ. PERMANOVA Permutational multivariate analysis of variance. Austral Ecol. 2005;26:32–46.

    Google Scholar 

  78. Zimmerman DW. Comparative power of student t test and Mann-Whitney U test for unequal sample sizes and variances. J Exp Educ. 1987;55:171–4.

    Google Scholar 

Download references

Acknowledgements

We would like to thank D. Liu and Q. Yan for sample collection.

Funding

This study was funded by the National Natural Science Foundation of China (No. 31672453), and W. Guo and T. Ma are supported by the China Scholarship Council (CSC). This study was also supported by an NSERC Discovery Grant awarded to L.L. Guan.

Author information

Authors and Affiliations

Authors

Contributions

WG and RL designed the study. WG, WW, and SB performed the experiments. YZ and DH provided suggestions on the study design. WG performed all bioinformatics and statistical analyses. MZ and TM guided the data analysis and revised the manuscript. WG and LLG interpreted the data and wrote the manuscript. The authors read and approved the final manuscript.

Corresponding authors

Correspondence to Le Luo Guan or Ruijun Long.

Ethics declarations

Ethics approval

All experiments were conducted in accordance with the Animal Ethics Committee of the Chinese Academy of Lanzhou University (Protocol No., SCXK Gan 20,140,215).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

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: Figure S1

Alpha diversity rarefaction curves of rumen bacteria, archaea, fungi and protozoa.

Additional file 2 Table S1.

The statistical analysis of alpha diversity from 7 days to 1 year based on repeated methods ANOVA (mixed-effects model). Table S2. The results of the pairwise PERMANOVA test. Table S3. The results of DESeq2 analysis. Table S4. The results of age discriminatory feature importance score. Table S5. The topological features of intra-interactions of core taxa in each rumen microbial kingdom. Table S6. The topological features of inter-interactions of core taxa among rumen microbial kingdoms. Table S7. The keystone species in the network within each age group.

Additional file 3 Figure S2A, S2B, S2C, S2D

Microbial UniFrac dissimilarity in each rumen microbial group in grazing yaks. Box plot showing within-group similarity, which was calculated using weighted and unweighted UniFrac metrics based on the average pairwise dissimilarity between each paired sample within different groups. A-D indicate bacteria, archaea, fungi and protozoa, respectively.

Additional file 4 Figure S3A., S3B, S3C, S3D

Taxonomic composition of each ruminal microbial kingdom at the phylum level. Bacteria (A), archaea (B), fungi (C), and protozoa (D).

Additional file 5 Figure S4A., S4B, S4C, S4D

Heatmap analysis of the relative abundance of rumen microbiota. The relative abundance was log10 transformed, and A-D indicate bacteria, archaea, fungi and protozoa, respectively.

Additional file 6 Figure S5

Co-occurrence network analysis of inter-interactions between core taxa among microbial kingdoms of rumen. Blue edges correspond to negative correlations, and red edges correspond to positive correlations. The size of the nodes is related to the relative abundance of the taxa. The lines in red and blue denote positive and negative correlations, respectively. The nodes are referred to as keystone species and are highlighted in light purple, and the taxonomic names of keystone species are indicated in the network figures. The taxonomic information of the nodes is shown at the bottom of the figure.

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

Guo, W., Zhou, M., Ma, T. et al. Survey of rumen microbiota of domestic grazing yak during different growth stages revealed novel maturation patterns of four key microbial groups and their dynamic interactions. anim microbiome 2, 23 (2020). https://doi.org/10.1186/s42523-020-00042-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42523-020-00042-8

Keywords