Three study regions, including Zhongba County, Xigatse City (84°03′ E, 29°77′ N; about 4,800 m altitude, HAL), Nagqu City (92°07′ E, 31°48′ N; about 4,500 m altitude, MAL), and Dangxiong County, Lhasa City (91°05′ E, 30°51′ N; about 3,800 m altitude, LAL), were selected. The three regions are key breeding grounds for yaks; their average annual air temperature is −4.0°C, −0.5°C, and 2.0°C, while precipitation is 280, 433, and 485 mm, respectively.

In September 2019, ten healthy female yaks (average 280 kg of BW) in 5–8 years of life were randomly selected from every study region. At the study regions, yaks are allowed to live free-range, follow their dams, and are not given any artificial feed, relying instead on typical alpine meadows as the main type of vegetation.

Before morning grazing, ruminal chyme (30 samples in total) was collected through the mouth of yaks using an esophageal tube as described previously (Ramos-Morales et al., 2014). Next, the chyme was filtered through a four-layer woven gauze to collect the ruminal fluid, which was immediately frozen in liquid nitrogen (−80°C) for subsequent microbiota and fermentation characterization. The pH value was detected using a portable pH meter (206-pH1; Testo). For VFA analysis, three tubes containing 1 ml of rumen fluid were taken from each sample after thawing and centrifugation at 20,000 × g and 4°C for 15 min. Then, 0.25 ml metaphosphoric acid (25 g/100 ml) was added per tube and the sample was analyzed by gas chromatography equipping a megabore HP-MOLSIV column (film thickness: 30 m × 0.53 mm × 25 μm) (SP-3420A; Beifenrili Analyzer Associates, Beijing, China) (Yue et al., 2009).

Microbial DNA was extracted from rumen fluid samples (30 samples in total) using PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA) following the manufacturer’s guidelines. The concentration and purity of the extracted DNA were measured by a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). The V3–V4 region of the bacterial 16S rRNA gene was amplified using primer 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Takahashi et al., 2014). PCR amplification proceeded through an initial denaturing step at 94°C for 5 min; 28 cycles at 94°C for 30 s, 55°C for 30 s, and 72°C for 60 s; and a final extension at 72°C for 7 min using 25-μl reaction volumes, containing 12.5 μl 2 × Taq PCR MasterMix, 3 μl BSA (2 ng/μl), 1 μl forward primer (5 μM), 1 μl reverse primer (5 μM), 2 μl template DNA, and 5.5 μl ddH2O. The PCR products were detected by 2% agarose gel electrophoresis, purified with the Agencourt AMPure XP kit (Beckman Coulter, La Brea, CA, United States) according to the manufacturer’s instructions, and quantified by QuantiFluor-ST (Promega, Madison, WI, United States). Purified amplicons were sequenced on an Illumina MiSeq-PE300 platform (Illumina Inc., San Diego, CA, United States), generating 2 × 300 bp paired-end reads.

The obtained paired-end reads in the original DNA fragments were merged using Flash version 1.20 (Mago and Salzberg, 2011), and then each sample was separated according to a unique barcode. After removing barcodes, primers, and splice variants, raw reads were obtained. To generate high-quality reads, a specific sliding window strategy was adopted for Trimmomatic version 0.36, the window size was set to 50 bp, the average quality value was 20, and the minimum reserved sequence length was 120 bp; in addition, sequences containing N were removed by Pear version 0.9.6 (Bolger et al., 2014). Afterward, Flash and Pear were used to merge the sequences at both ends according to PE overlap correlation, the minimum overlap was set to 10 bp, and the mismatch rate was 0.1, allowing for fasta sequences to be obtained. Finally, the chimera containing fasta sequences was eliminated by comparison with the Gold database using the UCHIME algorithm, whereas the unknown database was removed by de novo means (Knight, 2011). Simultaneously, the short sequences, failing to meet the requirements, were removed, while high-quality sequences of clean reads were obtained. Subsequently, these clean tags were clustered into operational taxonomic units (OTUs) based on a 97% sequence similarity threshold using the UPARSE algorithm in Vsearch version 2.7.1. In each OTU, the richest sequence was filtered as the representative sequence (Edgar, 2013). To derive species classification information corresponding to each OTU, the representative sequences were compared and analyzed by the RDP Classifier algorithm version 2.2 (Wang, 2007) and Silva database¹ (Quast et al., 2012), allowing community annotation at kingdom, phylum, class, order, family, and genus level. After the sample number of the lowest sequence was flattened, alpha diversity (Chao1, Shannon, Simpson, and observed species) was calculated by QIIME 2, while intergroup alpha index variability was demonstrated by the Kruskal–Wallis test in R version 4.0.2 (Miller et al., 2016). Principal coordinate analysis (PCoA) based on Bray–Curtis dissimilarity matrices and QIIME 2 was used to calculate differences in bacterial communities among groups. The linear discriminant analysis effect size (LEfSe, LDA > 3) was used to identify significant bacteria among the three groups (Miller et al., 2016). PICRUSt2 software was used to predict microbiota function and explore differences among the three groups (Douglas et al., 2019).

Rumen fermentation parameters were analyzed using one-way ANOVA and Duncan’s multiple comparison in SPSS 22.0 software (IBM, Chicago, IL, United States). Differences in alpha diversity, relative abundance at phylum, family, and genus level, as well as microbiota function among the three groups were tested using the Kruskal–Wallis method in R version 4.0.3. PCoA, Venn, column, and LEfSe results were visualized using the “ape” and “ggplot2,” “limma,” “ggplot2,” and “ggtree” packages in R, respectively. The “corrplot” package in R was used to analyze the correlation between the top 20 genera of all samples, rumen fermentation parameters, and altitude (based on Spearman’s coefficient). The network containing the top 20 genera was visualized using the “igraph” package in R version 4.0.3. All data are reported as means, and intergroup differences of p < 0.05 were considered significant.

As reported in Table 1, rumen fermentation parameters differed significantly between the three groups (p < 0.05), whereas the pH value had no difference between the three groups (p > 0.05). Also, the MAL group exhibited the highest (p < 0.05) total VFA content and molar proportions of propionic and butyric acids, whereas the molar proportion of acetic acid and A/P (acetic acid/propionic acid) were lower (p < 0.05) than in other groups. Interestingly, rumen fermentation parameters were similar between HAL and LAL groups (p > 0.05).

Gradual stabilization of the OTU rank curve chart indicates high coverage of test samples (Supplementary Figure 1). The rumen microbiota consisted of 5,491,884 high-quality reads, with an average of 183,063 reads per sample. A total of 7,977 OTUs were obtained based on 97% nucleotide sequence identification: 6,895 OTUs in HAL, 5,544 in MAL, and 6,109 in LAL groups (Figure 1A). Of these, 4,219 OTUs (52.89% of the total) were shared among samples from different groups. Alpha diversity calculations (Table 2) revealed significantly higher indices of Shannon and community richness (Chao1 and observed species) in HAL and LAL groups compared with the MAL group (p < 0.05), whereas no significant difference was detected between HAL and LAL groups (p > 0.05). Interestingly, we found Simpson index in HAL significantly higher (p < 0.05) than MAL group, whereas no significant difference was detected between LAL and other groups (p > 0.05). PCoA plots of bacterial structure profiles (Figure 1B) revealed distinct clustering of the three regions. The segregation and dissimilarities observed at OTU level were revealed using Bray–Curtis matrices (PERMANOVA, p < 0.001).

Next, we proceeded with bacterial identification of the 30 samples at different levels. At the phylum level, the top five ruminal bacteria of yaks across all samples were Bacteroidetes (50.69%), Firmicutes (41.11%), Proteobacteria (1.45%), Tenericutes (1.40%), and Actinobacteria (1.05%) (Figure 2A and Supplementary Table 1). The MAL group presented the lowest (p < 0.05) relative abundances of Bacteroidetes and Verrucomicrobia, but the highest (p < 0.05) relative abundance of Firmicutes compared with the HAL and LAL groups. At the family level, Prevotellaceae (24.67%) was the most abundant taxon, followed by Ruminococcaceae (14.76%) and Rikenellaceae (11.67%) across all samples; other families included Lachnospiraceae (9.15%), Christensenellaceae (8.42%), Bacteroidales_BS11_gut_group (7.20%), Veillonellaceae (4.05%), and Bacteroidales_S24-7_group (3.44%) (Figure 2B and Supplementary Table 2). Notably, the relative abundances of Lachnospiraceae, Christensenellaceae, and Veillonellaceae were higher (p < 0.05) in the MAL compared with other groups, whereas the relative abundance of Rikenellaceae increased significantly (p < 0.05) according to the following order: MAL < LAL < HAL. At the genus level, the dominant genera across all groups were Prevotella_1 (16.82%), Rikenellaceae_RC9_gut_group (11.26%), and Christensenellaceae_R-7_group (8.06%) (Figure 2C and Supplementary Table 3). Specifically, the MAL group displayed the highest (p < 0.05) relative abundances of Christensenellaceae_R-7_group, Ruminococcaceae_NK4A214 _group, Quinella, Ruminococcaceae_UCG-005, Olsenella, Butyrivibrio_2, Lachnospiraceae_NK3A20_group, and Acetitomaculum, but the lowest (p < 0.05) relative abundances of Prevotellaceae_UCG-003, Rikenellaceae_RC9_gut_group, Eubacterium_coprostanoligenes_group, Saccharofermentans, and Prevotellaceae_NK3B31_group compared with the LAL and HAL groups. The relative abundance of Ruminococcaceae_UCG-010 was significantly (p < 0.05) higher in the HAL group compared with the MAL and LAL groups, whereas the relative abundance of Papillibacter was higher (p < 0.05) in the LAL group than in the other groups.

To better understand the dominance of specific bacteria within the three groups, we used the LEfSe method (Figure 3). The Firmicutes and Saccharibacteria phyla, including genera Quinella, Christensenellaceae_R_7_group, Lachnospiraceae _NK3A20_group, Acetitomaculum, Olsenella, Ruminococcaceae _UCG_005, Butyrivibrio_2, and Ruminococcus_2, were abundant in the MAL group. The Bacteroidetes, Verrucomicrobia, and Synergistetes phyla, including Rikenellaceae_RC9_gut _group, Saccharofermentans, Lachnospiraceae_AC2044_group, Ruminococcaceae_UCG_010, Anaerovorax, and Family_XIII, were enriched in the HAL group. Finally, the Papillibacter, Tenericutes, Cyanobacteria, and Bacteroidales_S24_7_group phyla, including Prevotellaceae_UCG_003, Prevotellaceae_NK3 B31_group, Papillibacter, Anaeroplasma, and Lachnospiraceae _UCG_006, were over-represented in the LAL group.

The correlation between relative abundance of the top 20 bacterial genera, fermentation parameters, and altitude (Figure 4A) revealed that the relative abundances of Rikenellaceae_RC9_gut_group, Quinella, Prevotellaceae_UCG-003, Lachnospiraceae_NK3A20_group, Papillibacter, Ruminococcaceae_UCG-010, Prevotellaceae_NK3B31_group, and Ruminococcaceae_UCG-005 correlated with altitude and rumen fermentation (p < 0.05). For example, the relative abundance of Rikenellaceae_RC9_gut_group correlated positively with altitude (r = 0.84) and A/P (r = 0.47), whereas the relative abundance of Quinella correlated negatively with altitude (r = −0.49) and butyric acid (r = −0.43). The relative abundance of Papillibacter correlated positively with altitude (r = 0.41), total VFA (r = 0.62), propionic acid (r = 0.45), and butyric acid (r = 0.62), but negatively with acetic acid (r = −0.52) and A/P (r = −0.48). Interestingly, we found a significant negative correlation between Papillibacter and Ruminococcaceae_UCG-005 (|r|> 0.6, p < 0.05; Figure 4B).

The function of rumen microbial communities in yaks from three different regions was predicted by PICRUSt2 software and related to KEGG pathways (Table 3). Notably, “Biosynthesis of ansamycins” (4.47%) was the most abundant pathway in all three groups, followed by “Biosynthesis of vancomycin group antibiotics” (2.47%), “Valine, leucine and isoleucine biosynthesis” (2.25%), and “D-Glutamine and D-glutamate metabolism” (2.19%). The Kruskal–Wallis test confirmed that a total of 25 pathways showed significant (p < 0.05) differences among the three groups (Table 3). Specifically, “Amino acid metabolism,” “D-Alanine metabolism,” and “Thiamine metabolism” pathways, as well as “Carbohydrate metabolism,” “Cell motility,” and “Folding, sorting and degradation” gene categories were significantly (p < 0.05) more abundant in the MAL than in the other two groups.