Introduction

The Tibetan plateau, called ‘the third pole’, forms the largest and highest year-round grazing area in the world. The region is characterized by a harsh climate of extreme cold and aridity, and high ultraviolet radiation and hypoxia, which challenge the survival of humans and other mammals1. The yak (Bos grunniens), an iconic symbol of high altitude and a mainstay for Tibetan people2, has anatomical and physiological adaptations and a genetic basis for mammalian adaptations3, as well as a co-evolved microbiome4, that equip the animal for the high altitude and extreme environment. Yet gut microbiota of the yak and their relationship to seasonal dietary shifts in their natural habitat remain largely unknown, although this information could contribute to the understanding of adaptation to the high-altitude Tibetan plateau.

Gut microbiota are complex and dynamic5, being sensitive to perturbations, such as dietary changes, environmental factors6 and enteric pathogens. They play an integral role in nutrient intake, behavior, metabolism, immune function, and development of the host7,8. Substantial changes in mammalian microbiota composition have been observed in response to seasonal diet availability among and within individuals, as evidenced by longitudinal analyses of gut microbiota in Hadza hunter-gatherers9, wild wood mice (Apodemus sylvaticus)10, red squirrels (Tamiasciurus hudsonicus)6, giant pandas (Ailuropoda melanoleuca)11, wild great apes12,13 and North American bison (Bison bison)14. However, most studies on mammalian gut microbiota dynamics were done in non-stressful environments and without quantitative dietary information related to habitats. Recent reports on gut microbiota composition from large herbivores in the semi-arid East African savanna revealed a greater seasonal turnover and diet–microbiota association in domesticated than in wild species15. These studies provide a better understanding of intra-specific and inter-specific diet–microbiota associations in wild and domesticated species. Yet studies on seasonal diet and microbiota relations are lacking in high-altitude mammals.

To examine the fine-scale relationship between quantitative dietary consumption and gut microbiota, we conducted a spatio-temporal study of the impact of seasonal diet on gut microbiota from 302 individual yaks. Using DNA metabarcoding and 16S rRNA gene analyses, measurements were made in free-grazing yaks across four seasons on the eastern part of the Tibetan plateau, in which either a transhumance (TH) or open-continuous grazing (OCG) regime was followed. These in-depth and longitudinal analyses could contribute to the understanding of the adaptations of mammals to the harsh, high altitude environment.

Results

Established reference database for DNA metabarcoding analysis

According to the National Center for Biotechnology Information (NCBI) and Bold system databases, there were 199 plant species in the study areas. Combined with plant taxonomic identification, this study generated an additional 212 local plant species DNA barcode reference library using P6 loop of the chloroplast trnL (UAA) intron marker genes16 (see the “Methods” section). This revised library was used for DNA metabarcoding analysis. In total, our library comprised 411 plant species (included all species from the sampling region; local plant list recorded 386 species17) throughout the alpine grassland area of the study in the TH and OCG regimes.

Diet diversity and composition across seasons

The diet data included 30,534,414 high-quality sequences, 81,530 unique sequences after removal of singletons and 2010 operational taxonomic units (OTUs) (Supplementary Table 3). All dietary sequences represented 41 plant families, 83 genera, and 80 species. Dietary composition was more diverse in winter than in summer in both grazing regimes and displayed an evident separation seasonally (Fig. 1a, b and Supplementary Figs. 4, 5b, e). Gramineae and Rosaceae were consumed by yaks frequently throughout the year in both TH and OCG grasslands (Fig. 2a, c). The highest relative abundances were Polygonaceae, Rosaceae, and Gramineae in spring and summer, and Gramineae, Rosaceae, and Compositae in winter in both TH and OCG grasslands. In autumn, Gramineae and Rosaceae were highest in both TH and OCG grasslands, and Salicaceae was also high in TH. As indicator plant species, Polygonaceae was identified in spring, Scrophulariaceae and Compositae in winter in both TH and OCG regimes, Polygonaceae in TH and Rosaceae in OCG in summer and Salicaceae in TH and Gramineae in OCG in autumn (Fig. 2b, d; Supplementary Fig. 7a, b). The distinct seasonal diets provided a sound basis to identify the effects of seasonal diet patterns on yak gut microbiota.

Fig. 1: Seasonal changes in both diet and gut microbiota community structures of yaks in transhumance and open-continuous grazing regimes.
figure 1

Within and among seasons, Bray–Curtis dissimilarity in diet and microbiota are presented in Supplementary Table 2. Rows show the same ordinations for diet (a and b) and microbiota (c and d) compositions. Diet composition and gut microbiota represent transhumance (a and c) and open-continuous grazing (b and d) regimes. Individual yak diet compositions from samples collected in (a) spring (n = 32), summer (n = 33), autumn (n = 37), and winter (n = 45) in transhumance grassland (anosim analysis: R = 0.94, p = 0.0001; adonis analysis: R2 = 0.78, p = 0.0001), (b) spring (n = 31), summer (n = 39), autumn (n = 38), and winter (n = 47) in open-continuous grazing grassland (anosim analysis: R = 0.88, p = 0.0001; adonis analysis: R2 = 0.67, p < 0.0001), and gut microbiota compositions in (c) spring (n = 31), summer (n = 31), autumn (n = 37) and winter (n = 48) in transhumance grassland (anosim analysis: R = 0.50, p < 0.0001; adonis analysis: R2 = 0.16, p < 0.0001) and (d) spring (n = 31), summer (n = 37), autumn (n = 38), winter (n = 47) in open-continuous grazing grassland (anosim analysis: R = 0.47, p < 0.0001; adonis analysis: R2 = 0.16, p < 0.0001) plotted on nonmetric multidimensional scaling (NMDS) according to the Bray–Curtis dissimilarity. Analysis of similarities (ANOSIM), adonis analysis and permutational multivariate analysis of variance (PERMANOVA) were used for statistical testing of treatment similarities. The dotted ellipse borders represent the 95% confidence interval.

Full size image
Fig. 2: Seasonal changes of dietary compositions of yaks in transhumance and open-continuous grazing regimes.
figure 2

Stream-graph displays the relative abundance of plant family-level taxa in spring, summer, autumn, and winter in transhumance (a) and open-continuous grazing (c) regimes. Low abundance taxa (<5%) are grouped together as “others”. Indicator families that are related to each season are tracked using Sankey plots in transhumance (b) and open-continuous grazing (d) regimes. Lines represent associations between indicator families and seasons, which are colored by plant family. Line width is scaled to reflect indicator value (higher indicator value of family is more strongly associated with season). Indicator values are presented in Supplementary Fig. 7. The statistical p values mean the family associated with seasons, *p < 0.05, **p < 0.01, ***p < 0.001.

Full size image

In the OCG regime, above-ground biomass (AGB) was highest in the warm season (from June to September), with 2891 ± 148 kg DM/ha in September, and was lowest in the cold season (from October to May) with 99 ± 10.8 kg DM/ha in May (Fig. 3a). Neutral detergent fiber (NDF) and acid detergent fiber (ADF) were higher by 10% and 16%, respectively, while crude protein (CP) and ether extract (EE) were lower by 53% and 25%, respectively, in the cold than the warm season (Fig. 3b). Specifically, the highest CP (166 ± 10.2 g/kg DM) occurred in June, and the lowest (39 ± 5.5 g/kg DM) in April. Similar trends were observed in the TH regime (Supplementary Fig. 3a, b). Consequently, yaks were forced to cope with sparse forage of low-nitrogen content in the cold season.

Fig. 3: Seasonal dynamics in above-ground biomass (AGB) and chemical composition (dry matter basis) of the diets year-round in the open-continuous grazing grassland.
figure 3

a Line chart represents the AGB (kg DM/ha). The dashed line is the mean AGB year-round. b Line chart represents chemical composition (crude protein (CP), ether extract (EE), acid detergent fiber (NDF), and neutral detergent fiber (ADF)) of diet year-round. The dashed lines are the mean of each chemical composition year-round and are colored by each chemical composition. Values and error bars are shown as means ± SE. Average contents of CP are correlated negatively with those of NDF (R2 = 0.59, p < 0.001) (c) and ADF (R2 = 0.78, p < 0.001) (d). Statistic tests are performed by the t-test with FDR (false discovery rate) corrected p-value, *p < 0.05, **p < 0.01, ***p < 0.001.

Full size image

Gut microbiota diversity and composition across seasons

The gut microbiota included 25,375,771 high-quality sequences, 1,125,235 unique sequences after removal of singletons and 14,239 OTUs (Supplementary Table 3). Seasonal shifts in gut microbiota were evident in both the TH and OCG regimes (Fig. 1c, d). In total, 18 gut bacteria phyla were identified, with Firmicutes and Bacteroidetes the most abundant, regardless of season and grazing regime (Supplementary Fig. 6, Supplementary Fig. 7c, d). Distinct seasonal dynamics were exhibited at the genus-level in both TH and OCG grasslands, with Ruminococcaceae_UCG-010, Ruminococcaceae_UCG-005, Rikenellaceae_RC9_gut_group, unclassified Ruminococcaceae and unclassified Bacteroidales the most abundant across seasons (Fig. 4a, c). Moreover, the relative abundance of Prevotellaceae_UCG-004 increased during the summer. Indicators at the genus-level displayed seasonal and spatial differences. In the TH regime, in spring, autumn, and winter, the indicator species were Ruminococcaceae_UCG-010 and Clostridiales_vadinBB60_group; whereas, in summer, they were Ruminococcaceae_UCG-005 and Prevotellaceae_UCG-004. In the OCG regime, in spring and summer, the indicator species were Ruminococcaceae_UCG-010 and Clostridiales_vadinBB60_group; in autumn, they were Flavonifractor, Ruminococcaceae_UCG-005, and Prevotellaceae_UCG-004; whereas, in winter it was the unclassified Paludibacteraceae (Fig. 4b, d; Supplementary Fig. 7e, f).

Fig. 4: Seasonal variations of yak gut microbiota at the genus level in transhumance and open-continuous grazing regimes.
figure 4

Relative abundance of the 20 most abundant genera over seasons (spring, summer, autumn, and winter) are aggregated and colored on a stream-graph in (a) transhumance and (c) open-continuous grazing regimes. Low abundance taxa (except for 20 most abundant genera) are grouped together as “others”. Indicator genera that are related to each season are tracked using Sankey plots in transhumance (b) and open-continuous grazing (d) regimes. Lines represent associations between indicator genera and seasons, which are colored by genus level. Line width is scaled to reflect indicator value (higher indicator value of genus is more strongly associated with the season). Indicator value are shown in Supplementary Fig. 7. The statistical p values mean the genus associated with seasons. *p < 0.05, **p < 0.01, ***p < 0.001.

Full size image

Diet associated with overall microbiota composition across seasons

We applied Procrustes analysis to test diet and microbiota variations across seasons. When analyzed using Bray–Curtis (BC) dissimilarities, seasonal diet composition was associated with microbiota composition consistently in both TH (p = 0.0001, Fig. 5a) and OCG regimes (p = 0.0001, Fig. 5b). However, microbiota richness was not correlated with dietary richness (Supplementary Fig. 8).

Fig. 5: Seasonal clustering patterns in diet–microbiota lineages of yak in transhumance and open-continuous grazing regimes.
figure 5

The diet and microbiota compositions are determined with Bray–Curtis distance in (a) transhumance and (b) open-continuous grazing regimes, which are colored by season. Procrustes rotates the results of separate principal coordinates of diet composition (circle symbols) and gut microbiota composition (triangle symbols). Significant dissimilarities occurred between diet and microbiota ordinations by p-value.

Full size image

Diet and microbiota dissimilarities were consistently higher among than within seasons (Fig. 6 and Supplementary Table 2). Diet separations across seasons were smallest between spring and summer (0.482) in TH and between autumn and winter (0.566) in OCG (Supplementary Table 2), and were greatest between summer and winter (0.779) in TH and between spring and autumn (0.734) in OCG, suggesting that consumption of plant species occurred according to availability (Supplementary Table 2). Most overlapping microbiota occurred between autumn and winter while most non-overlapping microbiota occurred between spring and summer in both TH and OCG, which indicated that microbiota remained relatively stable (Supplementary Table 2). Within season, there was a larger difference in diet dissimilarity (0.117–0.385) than in microbiota dissimilarity (0.432–0.486) in the two grazing regimes, suggesting that stable gut microbiota communities may contribute to host adaptation to the extreme environment (Fig. 6 and Supplementary Table 2).

Fig. 6: Bray–Curtis dissimilarities within seasons in dietary diversity and gut microbial diversity across seasons in transhumance and open-continuous grazing regimes.
figure 6

Rows show the dissimilarities for diet (a and b) and gut microbiota (c and d) within seasons. Each column organizes the data so that diet and gut microbiota Bray–Curtis dissimilarities represent transhumance (a, c) and open-continuous grazing (b, d) regimes. All boxplot distributions are tested by non-parametric Kruskal–Wallis and Wilcoxon with FDR (false discovery rate) corrected p-value, center values indicate the median and error bars. *p < 0.05, **p < 0.01, ***p < 0.001.

Full size image

Gut enterotypes and functional context represented by Akkermansia and uncultured Eubacterium WCHB1-41 for underlying cold adaptation

Based on the report that enterotypes exhibit functional differences18, we examined whether yak gut microbiota partitioned into clusters that differ in functional properties according to seasonal dietary intake. Principal component analysis (PCA) revealed that the samples formed three distinct enterotype clusters based on BC dissimilarities. Each cluster was driven by the variation of its representative genera level: Akkermansia and uncultured Eubacterium WCHB1-41 in Enterotype 1, Ruminococcaceae_UCG-005 in Enterotype 2 and Ruminococcaceae_UCG-010 in Enterotype 3 (Fig. 7a, d–g; Supplementary Figs. 9 and 10). In TH and OCG, there was a change in enterotype assignment across seasons. Enterotype 1 occurred predominantly in the cold season (spring, autumn, and winter), Enterotype 2 in the warm season (summer) (Fig. 7b, c; p < 0.05, Fisher’s exact test) and Enterotype 3 was prevalent throughout the year (Fig. 7b, c). This study identifies the distribution of different enterotypes across seasons in high altitude yaks.

Fig. 7: Enterotype distributions of yak gut microbiota associated with season using Bray–Curtis dissimilarity.
figure 7

Identification of yak enterotypes is presented in Supplementary Fig. 9. a Visualizations of enterotypes, as identified by PAM (partitioning around medoid) clustering. Genera corresponding to each enterotype are identified by their relative abundance (see Supplementary Fig. 10). b and c Proportion of samples for each enterotype in spring, summer, autumn, and winter in (b) transhumance and (c) open-continuous grazing regimes. d–g Relative abundance of bacterial taxa characteristic of each enterotype. Ten genera were chosen based on their average contribution to overall Bray–Curtis dissimilarity. All six bacterial genera are presented in Supplementary Fig. 10. Colors correspond to enterotype clusters. All bar distributions are tested by Fisher’s exact test with FDR (false discovery rate) corrected two-tailed p-values (b and c). All boxplot distributions are tested by non-parametric Kruskal–Wallis and Wilcoxon with FDR-corrected p-value, center values indicated the median and error bars (d–g). *p < 0.05, **p < 0.01, ***p < 0.001.

Full size image

These seasonal distributions of gut enterotypes led us to hypothesize that the fixed gut enterotype, represented by Akkermansia and uncultured Eubacterium WCHB1-41, plays a vital role in regulating nutritional requirements in the cold season with sparse forage. To test this hypothesis, we examined the functional relevance of Akkermansia and Kiritimatiellaeota based on genus-level pan-genomes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Both Akkermansia and Kiritimatiellaeota showed convergent enrichment of enzymes that are involved in arginine and fatty acid biosynthesis pathways (map00220 and map00061) (Fig. 8). Pyruvate generated acetyl-CoA, mediated by both Akkermansia and Kiritimatiellaeota, enters the tri-carboxylic acid (TCA) cycle and involves arginine biosynthesis, and also regulates fatty acid synthesis. Notably, we observed 12 enzymes that play key roles in arginine biosynthesis, but no enzyme is involved in urea synthesis to reduce nitrogen loss in urine under conditions of low-nitrogen stress. Six of a total of seven enzymes participated directly in fatty acids biosynthesis and contributed to energy deposition. These results indicated that both arginine and fatty acid biosynthesis pathways evolved in high altitude mammals for efficient nitrogen utilization and energy deposition in the cold season.

Fig. 8: Metabolic pathway associated with Akkermansis and uncultured Eubacterium WCHB1-41_ge on host metabolism.
figure 8

The compositions of Akkermansis and uncultured Eubacterium WCHB1-41_ge are the key factors that influence arginine biosynthesis and fatty acid biosynthesis and were significantly higher during the cold season. These microbiota may contribute to energy production during cold season, which is the period of insufficient dietary and protein intake. All enzymes and EC (Enzyme Nomenclature) numbers were obtained from Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Full size image

Discussion

We determined seasonal shifts in diet and gut microbiota in a high-altitude large herbivore. Yaks graze the grasslands of the Tibetan plateau all year, and face severe challenges, in particular extreme cold and limited food availability during the cold season. To survive under such conditions, yaks have evolved anatomical and physiological adaptations. Furthermore, convergent evolution was reported for yaks and their rumen microbiota. The rumen microbiota follow unique maturation strategies19 and differ in composition from microbiota of lowland ruminants4.

In this study, seasons constrained both diet and gut microbiota, but less so for microbiota than diet. Such effect of season on microbiota composition in response to food availability and chemical composition might be due to a re-assembly of the community structure of the gut microbiota resulting from substrate availability. Surprisingly, among seasons, the changes in gut microbiota composition were relatively stable compared with changes in dietary composition, which suggests that high altitude mammals evolved stable patterns of gut microbiota composition across seasons. By contrast, gut microbiota in humans and other mammals differed markedly with seasonal dietary shifts9,10. In response to cyclic seasonal dietary fluctuations of the host, gut microbiota could alter their metabolic rate and energy extraction from complex carbohydrates, and, ultimately, promote the co-evolution of hosts and microbes6. Furthermore, our findings support the present grassland polices of TH of moving yaks to different landscapes and searching for favorable diet. This practice mitigates grazing pressure, thereby enhancing animal performance. As demonstrated by our diet and microbiota analyses, a better understanding of relations between yaks and plants can yield key insights in supporting TH with a proper stocking rate and improving grassland management strategies. The key microbiota associated with feed efficiency could provide an alternative solution for strengthening microbiota-led breeding programs to increase yak performance, particularly in the cold season.

The lack of precise methods to assess and identify the complex diet consumed by grazing herbivorous mammals has been a long-standing obstacle in field studies. Conventional methods include field observations, microhistology and natural n-alkanes from plant cuticular waxes20,21. However, these methods are not suitable for diet analyses of herbivorous mammals feeding on pasture with high plant diversity. Field observations require high visibility and are prone to omission of plant species22, while histology is very tedious and often inaccurate23. The n-alkane method is limited when the animal consumes a large number of plant species24. When the n-alkane method was used to determine diet composition of grazing yaks on the Tibetan plateau, only the two or three main dietary components could be detected20. Studies of DNA metabarcoding identified the level of forb dominance in megafaunal diets from permafrost sediment samples25 and allowed fine-grained niche separation from seven large mammalian herbivores26. These studies indicated that ingested plants can be effectively identified and quantified using the DNA metabarcoding method27.

This present study provides biological insights into the clusters and functionality of gut enterotypes. In particular, we describe gut enterotypes as well as the functional genomic information following seasonal dietary changes in free grazing yaks. The Ruminococcaceae_UCG-005 enterotype was associated with high protein and low fiber diets in the warm season and the Akkermansia and uncultured Eubacterium WCHB1-41 enterotype with low protein and high fiber diets in the cold season. It was suggested that the proportion of protein and carbohydrate contents in the diet mediates the host’s enterotype shift, at least in baboons28,29. Seasonal variations in dietary protein and carbohydrate contents could provide an attractive explanation in terms of enterotype dynamics, and also contribute to determine the enterotype for high altitude herbivores. Gut enterotype in humans remained stable when on a 10-day low-fat/high-fiber diet intervention30. However, the most compelling evidence for seasonal diet shifts in enterotypes in the present study was that enterotypes 1 and 2 shifted consistently between seasons and Ruminococcaceae_UCG-01 enterotype 3 remained stable throughout the year. This might be the result of long-term co-evolution between the host and environment and suggests that the fixed seasonal enterotype dynamics play a key role in species formation and adaptation of mammals to the extreme environment at high altitude.

Akkermansia resides in the mucus layer of the gut and is known for its degradation of mucin9,31. It was reported that A. muciniphila plays a pivotal role in the prevention of obesity and type 1 and type 2 diabetes and, hence, promotes health in humans32,33. In addition, there is evidence that A. muciniphila functions as an energy sensor and regulates energy homeostasis for host microbial mutualism34,35. The present study demonstrated a substantial increase in Akkermansia and uncultured Eubacterium WCHB1-41 enterotype 1 with low nitrogen and low energy intake in the cold season, suggesting that this enterotype responds to a high fiber and low protein diet. Previous studies reported that a high fat diet decreased mucus thickness32, whereas, a fiber-rich diet increased mucus thickness, which improved the gut positive barrier function36. Akkermansia and uncultured Eubacterium WCHB1-41 degraded mucins and converted them into short-chain fatty acids (SCFAs), which provided nutrients for other resident bacteria and cells37.

Genetic studies in yaks demonstrated enrichments in the amino-acid metabolism gene (Whsc1, Glul), fatty acid biosynthesis and metabolism gene (Hsd17b12)3, and up-regulation in 36 genes that were related to volatile fatty acids transport and absorption in the ruminal epithelium of high-altitude ruminants4. These studies provide additional adaptive responses to insufficient energy intake that yaks experience. In addition, it was reported that yaks have a relatively low N requirement and utilize dietary N more efficiently than cattle in the cold season38,39,40,41. The current study also provides a mechanistic explanation as Akkermansia and Kiritimatiellaeota are involved in arginine and fatty acid biosynthesis pathways to utilize dietary N and generate energy in the cold season. In the context of the increased metabolic disease prevalence, this study provides an insight on the influence of Akkermansia and uncultured Eubacterium WCHB1-41 enterotype in nitrogen and energy utilization in the host, and implies that this enterotype plays a key role in mediating the nutritional homeostasis in high altitude animals with large therapeutic potential of metabolic diseases in humans. Further analyses of genomic functions of gut microbiota of yaks would be beneficial for a better understanding of mammals living in an extreme environment. Findings in the present study demonstrate that gut microbiota respond to seasonal dietary shifts, which allows yaks to better utilize poor forage of low protein content. The understanding of diet–gut microbiota interaction improves our understanding of how yaks adapt to extreme environments.

Methods

Fecal sample collection and processing

The Tibetan plateau forms the high-altitude core of Asia and presents a challenge for mammals to survive1. In this study, yaks grazed on the Tibetan plateau in both TH and OCG regimes (Supplementary Fig. 1). In OCG, the yaks grazed freely on the same area at 3010–3300 m above sea level (2000 ha; 37°10′–37°12′N, 102°44′–102°51′E) year-round. Dominant and associated plant species included Kobesia humilis, K. capillifolia, Polygonum viviparum, Stipa capillata, Elymus nutans, Thalictrum alpinum, Medicago ruthenica, and Artemisia smithii. In the TH regime, yaks grazed pasture at 2930–3000 m above sea level (13.3 ha; 37°12′N, 102°46′E) in winter-spring; at 3130–3300 m above sea level (666.7 ha; 37°10′N, 102°44′E) in summer; and at 3015–3100 m above sea level (11.3 ha; 37°11′N, 102°44′E) in autumn. The dominant and associated species in winter-spring pasture included E. nutans, K. humilis, K. capillifolia, S. capillata, Leymus secalinus, P. viviparum, and A. smithii; in summer included P. fruticose, P. viviparum, Juncus himalensis, Deschampsia cespitosa, Festuca ovina, Saussurea amara, Carex atrofusca, and C. moorcroftii; in autumn included P. fruticose, Salix oritrepha, P. viviparum, E. nutans, D. cespitosa, K. tibetica, K. humilis, C. kansuensis, J. castaneus, P. anserine, T. alpinum, and S. katochaete. The annual mean temperature was −0.1 °C in the TH and 0 °C in the OCG with a peak in summer (June–August) and a trough in winter (December–February). The grassland was alpine meadow and alpine shrub meadow and the plant growing season was 90–120 days.

According to previous studies42,43, and the average monthly temperature and precipitation recoded by a nearby meteorological station in the study area, four seasons can be identified as spring (April–May), summer (June–August), autumn (September–October), and winter (November–March the following year) (Supplementary Fig. 2). In this study, sampling periods spanned the four seasons in 2017, namely, spring (May), summer (August), autumn (October), and winter (December) (for TH: n = 32, spring; n = 33, summer; n = 37, autumn; n = 45, winter and for OCG: n = 31, spring; n = 39, summer; n = 38, autumn; n = 47, winter). Fresh yak feces were collected from the TH and OCG grasslands, mixed thoroughly in an unused freezing tube, placed immediately into liquid nitrogen containers in the field, and transported to Lanzhou University until further processing. An amount of 0.2 g of fresh feces was used for DNA extraction with QIAamp® Fast DNA Stool Mini Kit (50, QIAgen GmbH) and an extraction blank was processed to monitor for cross-contamination. DNA was quantified using the NanoDrop-2000 UV–Vis Spectrophotometer (Thermo Scientific, Wilmingtoo, DE, USA). The DNA samples were used for both diet (n = 302) and microbiota (n = 300) analyses. The studies and all procedures involving the animals were approved by experimental field management protocols (EAF2021012) of Lanzhou University.

AGB and chemical composition analysis

Three plots (100 m × 100 m) were selected randomly in each of the TH and OCG grasslands in each of the four seasons. Within each plot, five quadrats (50 cm × 50 cm) were selected randomly. AGB was harvested in each quadrat, oven-dried at 65 °C to a constant weight, ground to pass through a 1-mm screen and stored at room temperature for chemical composition analysis. Dry matter (DM) was determined by oven-drying at 105 °C for 24 h, EE by extraction with petroleum ether and N content by the Kjeldahl method44. NDF and ADF were measured according to Van Soest et al.45, with heat stable alpha amylase and sodium sulfite used in the NDF procedure.

Diet DNA metabarcoding

The P6 loop of the chloroplast trnL (UAA) region was used for DNA metabarcoding with primers trnL (UAA) g and trnL (UAA) h16,25,26,46 (Supplementary Table 1). For the PCR assays, 10 μL reactions of each of 0.3 μL primers, 0.2 μL KOD FX Neo, 2 μL dNTP, 5 μL KOD FX Neo buffer, and 50 ng of DNA template were mixed. Thermocycling followed a program of initial denaturing at 95 °C for 4 min, followed by 35 cycles of 94 °C for 30 s, 50 °C for 30 s, and 72 °C for 1 min, with a 5-min final extension at 72 °C. All PCRs were conducted with a no-template negative control and a positive control (consisting of DNA extracted from plant species from our local DNA reference library). The 5′ end of each primer was tagged by a 16-nt multiplex identification tag that differed by 8-nt from the other tag, allowing uniquely tagged PCR products. The sequence was carried out on Illumina HiSeq 2500 platform.

16S rRNA gene Illumina sequencing

The V3–V4 region of the 16S rRNA gene was sequenced on Illumina MiSeq 2500 platform with primers (341F/806R)47 (Supplementary Table 1). For the PCR assays, 50 μL of each of the 30 ng DNA template, fusion primer, and PCR master mix were mixed. The PCR cycles started with a 3 min denaturation at 94 °C, followed by 30 cycles each consisting of 94 °C for 30 s, 56 °C for 45 s, 72 °C for 45 s, and followed by a final step of 72 °C for 10 min. PCR products were purified with AmpureXP beads and eluted in elution buffer. Libraries were qualified by the Agilent 2100 bioanalyzer (Agilent, USA). The amplicons were sequenced on Illumina MiSeq 2500 and generated 2 × 300 bp paired-end reads.

Reference plant DNA libraries

To identify diet plant sequences from fecal samples, we established an extensive DNA reference database from plant species throughout the alpine grassland area of the study, Yongfengtan and Wushaoling. The collection included 212 species that were most abundant in the study area. All plants were identified to species-level by expert botanists.

Reference plant DNA was extracted with DNeasy Plant Mini Kit (50, QIAgen GmbH) using 0.2 g of leaves, and sequenced chloroplast trnL-P6 (UAA) using established primers and protocols (Supplementary Table 1)26,48. DNA was quantified using the NanoDrop-2000 UV–Vis Spectrophotometer (Thermo Scientific, Wilmingtoo, DE, USA). trnL (UAA) was sequenced in 25 μL PCR reaction that included 2.5 μL MgCl2, 4 μL dNTP, 0.5 μL of each primer (trnL(UAA)c/trnL(UAA)d), 12.5 μL 2 × GC Buffer I, 0.25 μL TaKaRa LA Taq®, and 0.5 μL DNA template. Thermocycling for trnL (UAA) proceeded at 94 °C for 1 min, 35 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 1 min, with a 5 min extension at 72 °C. Programs Geneious and MEGA7.0.14 were used for sequence alignment and analysis.

DNA metabarcoding sequence analysis

Sequence demultiplexing, quality, and preliminary identifications were conducted by QIIME 1.9.1. Demultiplexing used split_libraries_fastq.py49. Sequences shorter than 10 bp and mean Illumina fastq quality scores <20 were not considered. Paired reads were merged using USEARCH1150 and then all merged sequences of each sample were pooled. Quality filtering was performed on the pooled sequence with more than 0.5 expected error using fastq_filter command in USEARCH11 and only sequence lengths ≥ 10 were retained. Fastx_uniques command in USEARCH11 was used to find a set of unique sequences from filtered sequences, and singletons (sequence abundance = 1 across all samples) were removed. The remaining sequences were denoised (cluster at 100% similarity) using UNOISE algorithm51, during which OTU representative sequences were generated and potential chimeras were excluded. All the initially pooled sequences were mapped into the denoised sequence (zOTU) to generate an OTU table using otutab command implemented in USEARCH11. Strictly identical sequences were merged and assigned plant species based on their unique sequences to DNA metabarcode sequences with exact matches (100% identity) to reference sequences. Only unique sequences with an identity of 100% to reference sequences were kept for further analysis. When 100% identities were acquired from the local reference and NCBI libraries; preference was given to the local reference library. When a diet sequence matched multiple reference sequences exactly, assignments were revised to the finest taxonomic level by blasting with NCBI. We used the summarize_taxa command to group identical sequences, tally them within samples and then quantify the relative read abundance of each sequence, which is widely used to quantify the proportional foods consumed by animals46,52,53,54, and has been confirmed in several studies26,55,56. The resulting OUT counts per sample were rarefied to 40,000.

Microbiota community analysis

Quality control, merging of pair ends, OTU clustering, and taxonomic assignation were performed using the QIIME 1.9.1. Illumina fastq quality scores < 20, ambiguous nucleotides and chimeras were discarded. The reads were assigned to OTUs using UNOISE351 with a threshold of 100% identity and seeded with SILVA rRNA gene databases57. Sequences identified as archaea, mitochondria, and chloroplast were removed. After filtering and identification, the bacterial 16S rRNA gene data included sequences across 300 samples. The resulting OTU counts per sample were rarefied to 10,000.

Application of enterotype clustering methodology

We applied methods described in humans18,29,30,58 to test for the presence of enterotypes in high-altitude yak. The genus-level relative abundance profiles of samples were clustered using Jensen–Shannon divergence (JSD) and BC dissimilarity and partitioning around medoid (PAM) clustering in R. The robustness of clusters was assessed by the Calinski–Harabasz (CH) index and silhouette score59. We applied the PAM, CH index, and silhouette score to clustering using BC and JSD methods, for which results did not differ (Supplementary Fig. 9). Furthermore, previous studies suggested that BC is related strongly to JSD29,58 and based on abundance method and suitable for revealing variations in abundance taxa, especially those with enterotypes. Thus, BC was implemented in genus-level abundance. To identify genus taxa contributing to enterotype groups based on BC, we applied the SIMPER method29,60, which identifies genus taxa contributing to similarity within- and dissimilarity between enterotypes and ranks their contribution.

Statistical analysis

Standard R commands were performed to generate variations in relative abundance across seasons, and the Wilcoxon test (two-sample comparisons) or the Kruskal–Wallis test (multiple groups) was used to measure significance in non-parametric relative abundance profiles. The t-test was used to measure significance in AGB and chemical composition of the diet profiles compared with the average AGB or chemical composition across seasons. R was performed to visualize the seasonal dynamics in AGB and chemical composition of the diet as line chart and relative abundance of diet and microbiota taxa across seasons as streamplots or boxplots. We used BC dissimilarity and nonmetric multidimensional scaling (NMDS) in vegan61. Pairwise differences within and across seasonal diets and microbiota variation were permutested with 9999 permutations and false discovery rate (FDR) correction. To assess which plant and microbiota taxa were most responsible for seasonal differences in diet and microbiota variations, we performed indicator species analysis62 with indicspecies in R. We used the multipatt function with 9999 permutations to the list of species that were related with a group of samples and r.g. function determined the correlation between two binary vectors. Within seasons, the relationship between diet and microbiota richness in each sample was assessed by linear regression using the data across seasons in TH and OCG regimes. The correlations between diet and microbiota composition were computed based on Procrustes analysis, a program that compares the relative positions of points in two multivariate datasets63, and conducted in R using the vegan package. Monte Carlo p-values for rotational agreement significance testing were determined from 9999 permutations. For enterotype comparisons, samples were pooled into bins (spring, summer, autumn, and winter), and the significance among seasons were identified using Fisher’s exact test with FDR correction of p-values. FDR was applied at a level of 0.05 per tested correlation and significance for multiple comparison.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.