The rumen microbiome of yak co-evolves with its host probably adding the adaptation to its harsh environments

Background Rumen microbes play an important role in ruminant energy supply and animal performance. Previous studies showed that yak (Bos grunniens) rumen microbiome and fermentation differ from other ruminants. However, little is understood on the features of the rumen microbiome that make yak adapted to its unique environmental and dietary conditions. Here we investigated the rumen microbiome and metabolome to understand how yak adapts to the coarse forage and harsh environment in the high Qinghai-Tibetan plateau. Result Metataxonomic analysis of the rumen microbiota revealed that yak (Bos grunniens), domesticated cattle (Bos taurus), and dzo (a hybrid between the yak and domestic cattle) have distinct rumen microbiota. Metagenomic analysis displayed a larger gene pool encoding a richer repertoire of carbohydrate-active enzymes (CAZymes) in the rumen microbiome of yak and dzo than cattle. Some of the genes encoding glycoside hydrolases (GH) that mediate the digestion of cellulose and hemicellulose were signicantly enriched in the rumen of yak than cattle, but the cattle rumen microbiome had more genes assigned to GH57 that primarily includes amylases. The rumen fermentation prole differed also, with cattle having a higher molar proportion of acetate but a lower molar proportion of propionate than dzo and yak. Metabolomic analysis showed differences in both rumen microbial metabolic pathways and metabolites, mainly amino acids, carboxylic acids, sugars, and bile acids. Notably, styrene degradation, primary bile acid biosynthesis, glyoxylate, and dicarboxylate metabolism signicantly differed between cattle and dzo; streptomycin biosynthesis was signicantly different between cattle and yak; and the pathways for biotin metabolism and styrene degradation signicantly differed between dzo and yak. Correlation analysis revealed certain microbial species correlated with differential rumen metabolites. Nine differential metabolites showed a positive correlation with seven species belonging to Bacteroides and Alistipes but a negative correlation with ten species belonging most of the cellulase transcripts, most of hemicellulose transcripts in the rumen The high of Ruminococcus, Fibrobacter, and Prevotella and the activities of carboxymethyl cellulase and avicelase in the rumen of yak corroborate its greater ability to digest crude ber. Xue et al. the microbial composition between cows with high and low milk yield and milk protein content, and they found that the rumen species enriched in the high-yielding cows were mainly members of Prevotella, and P. oralis, P. buccae, P. albensis, and R. avefaciens were signicantly higher in the high-yielding cows. We conclude that the above four bacterial species are positively correlated with feed eciency. Interestingly, they were also signicantly more predominant in the rumen of yak than the other two ruminant species. These species, among others, may help enable yak to adapt to its poor diet. In the present study, yak was found to possess more hemicellulases (e.g., GH11, GH5, GH44, GH16, GH17) and rhamnogalacturonan endolyase (e.g., PL11and PL4-4) than cattle. Many CBMs associated with galactose binding (CBM61, CBM62, CBM51) and L-rhamnose binding (CBM67) were also enriched in yak. These results indicate that yak probably has a stronger activity in removing galactose from the main chain of hemicellulose. Rapid hemicellulose digestion can facilitate cellulose digestion in forage because

alteration of rumen microbiota can impact rumen function and energy utilization in the host body, rumen microbiota may also contribute to host adaptation. We hypothesized that besides the adaptative evolution of the respiratory and the circulatory systems, the digestive system, especially the rumen microbiota, of yak had probably also undergone adaptative evolution to ingest and digest the available local feed. However, it remains to be determined if and how the yak rumen microbiota helps its host to adapt to its harsh environment.
Zhang et al. [19] compared the rumen metagenome and rumen epithelial metatranscriptome of cattle and yak living at different altitudes. Their results showed that compared with cattle living at low altitudes, the rumen microbiome of the yak is enriched with VFA fermentation pathways and the yak rumen wall is more effective in absorbing VFA. That study suggests the contributions of the rumen microbiome to the adaptive evolution of ruminants living at high altitudes. However, their VFA results were obtained in vitro, and they only focused on VFA production and methanogenesis. Because diet, environment, and feeding can profoundly affect the rumen microbiome [23,24], these confounding factors should be eliminated in comparative studies of different species or breeds of ruminants. The objective of this study was to elucidate the potential mechanism by which the rumen microbiota contributes to yak adaption. To achieve this goal, we compared the rumen microbiota (composition and structure), fermentation, and function (both the metagenome and metabolome) among yak, cattle, and dzo (designated as a species), all of which were kept at the same high altitude and fed the same diet. This study provided new knowledge of the yak rumen microbiome that might help understand its adaptation to high-altitude environments.

Results
Yak has the highest rumen propionate molar proportion and cellulase activity.
We determined the rumen fermentation characteristics and cellulase activities of fresh rumen samples to compare the feed fermentation and cellulolytic capacities among yak, Qaidam yellow cattle, and dzo. No signi cant difference was observed in rumen pH among the three ruminant species (Table 1), but the rumen concentration of total VFA in dzo was lower (P < 0.05) than in cattle. The molar proportion of propionate and butyrate in yak was higher than in cattle (P < 0.05), while the molar proportion of acetate and acetate : propionate (A:P) ratio showed the opposite trend. The rumen microbiome of yak had the highest (P < 0.05) activity of carboxymethyl cellulase and avicelase (Table 1). Different superscripts in a row designate signi cant difference (P < 0.05).
The rumen microbiota structure and composition are different among the three species.
Using metataxonomic analysis, we compared the rumen microbiota among the three ruminant species. The sequencing depth coverage was > 98% for bacteria and 99% for fungi for all the samples. The three species had different bacterial and fungal rumen microbiota. Yak had a lower (P < 0.05) ACE richness estimate of rumen bacteria than cattle and dzo, but no signi cant difference was observed in Simpson diversity index among the three ruminant species. Simpson evenness was higher (P < 0.05) for yak than for cattle and dzo (Additional le 1: Table S1). With respect to fungi, cattle had a lower (P < 0.05) observed species richness (both ACE and Chao1 estimates) than dzo and yak, but Simpson evenness for cattle was higher (P < 0.05) compared with yak and dzo (Additional le 2: Table S2). The β-diversity of the bacterial and fungal microbiotas was compared using PCoA based on Bray-Curtis dissimilarity, and ANOSIM analysis showed a difference (P = 0.001) in the overall rumen microbiota of both bacteria and fungi among the three ruminant species (Fig. 1a, b). Cattle and dzo had more similar rumen bacterial microbiota compared to yak, while yak and dzo had a more similar fungal microbiota compared to cattle. Overall, the rumen microbiota of bacteria and fungi of cattle clustered separately from that of yak, while that of dzo fell in between.
We further investigated the composition and function of rumen microbiota of the three ruminant species using metagenomics. We selected six samples from each ruminant species based on the metataxonomic results to maximize depth coverage with the budget we had. We obtained a total of > 1.36 billion raw reads (> 75.8 million reads per sample) totaling > 206 Gb of sequence data (> 11.4 Gb per sample). After the removal of host DNA and quality ltering, we obtained > 1.35 billion reads in total and > 75.1 million reads per sample. The average length of the Open reading frames (ORFs) was about 494 bp. PCoA analysis of the metagenomic data also revealed difference (P = 0.001) in the overall bacterial microbiota among the three ruminant species (Fig. 1c). Bacteroidetes, Firmicutes, Proteobacteria, Euryarchaeota were the predominant phyla in all samples, with Bacteroidetes (53.71%) and Firmicutes (24.74%) being the most predominant. The relative abundance of Proteobacteria, Verrucomicrobia, Planctomycetes, Synergistetes, Cyanobacteria, Chloro exi, and Acidobacteria was higher (P < 0.05) in yak than in cattle, while Chlamydiae was more predominant (P < 0.05) in cattle and Actinobacteria more predominant (P < 0.05) in dzo (Additional le 3: Figure S1). At the genus level, Bacteroides, Alistipes, Butyrivibrio, Faecalibacterium, and Pseudobutyrivibrio were more predominant in cattle than in yak (P < 0.05), while the opposite is true (P < 0.05) for Ruminococcus, Selenomonas, Oscillabacter, Treponema, and Fibrobacter (Additional le 4: Figure S2). At the species level, nine species, most of which belong to Bacteroides and Alistipes, were more predominant (P < 0.05) in cattle than in yak or dzo; four species were more predominant (P < 0.05) in dzo; and 18 species (most belonging to Prevotella) were more predominant (P < 0.05) in yak (Fig. 2).
Because the yak rumen had a higher propionate molar proportion (Table 1), we compared the relative abundance of bacterial species that either produce or utilize succinate, a fermentation intermediate of propionate production via the succinate pathway [25]. Of the 12 succinate-producing bacterial species, Prevotella brevis, P. bryantii, P. albensis, Fibrobacter succinogenes, Succinimonas amylolytica, and Mitsuokella jalaludinii were enriched (P < 0.05) in the yak rumen compared to the rumen microbiota of cattle or dzo (Fig. 3a), so was the succinate-utilizing species identi ed, Selenomonas ruminantium (Fig. 3b). The bacterial genera likely using the acrylate pathway for propionate production did not differ (P > 0.05) in relative abundance among the three ruminant species (Additional le 5: Figure S3). The relative abundance of Ruminococcus avefaciens, R. albus, and F. succinogenes, the best-known species of cellulolytic bacteria in the rumen, were enriched (P < 0.05) in yak compared to cattle (Additional le 6: Figure S4).
The yak rumen microbiota has a larger brolytic capacity but a smaller amylolytic capacity.
We further compared the functional potentials of the rumen microbiota among the three ruminant species. Annotation of the metagenomic data using the KEGG database followed with LEfSe analysis showed that replication and repair and endocrine system were enriched in cattle (P < 0.05); dzo had more (P < 0.05) metabolic pathways of other amino acids; whereas metabolism of cofactors and vitamins was represented more abundantly (P < 0.05) in yak (Fig. 4a). At level 3 (Fig. 4b), cattle, yak, and dzo had 24, 18, and 5 metabolic pathways enriched (P < 0.05), respectively, compared with the other two ruminant species. The metabolism of pyrimidine, purine, starch and sucrose, cyanoamino acid, and aminoacyl-tRNA biosynthesis were the top ve for cattle; amino acids biosynthesis, pantothenate and CoA biosynthesis, bio lm formation-Vibrio cholerae, lipopolysaccharide biosynthesis, cysteine and methionine metabolism were the top ve for yak; and carbon metabolism, butanoate metabolism, citrate cycle, sulfur relay system, and selenocompound metabolism were the top ve for dzo. Differential modules based on LEfSe analysis were showed in Additional le 10: Figure S5.
Rumen metabolite pro les are different among the three ruminant species.
We also used metabolomics to comparatively examine the metabolite pro les of the rumen microbiome among the three ruminant species.
Different rumen microbiota may lead to different metabolites.
We assessed the relationship between the relative abundance of bacterial species (shown in Fig. 2) and rumen metabolites (shown in

Discussion
The yak is indigenous livestock on the Qinghai Tibetan Plateau, and it is raised at altitudes between 3,000 and 5,000 m [29]. The high altitude and the environmental conditions are not suitable for domestic cattle. The cattle used in this experiment are Qaidam yellow cattle, which is well adapted to an altitude of 2,800 meters, but its pulmonary artery pressure is still signi cantly higher than that of indigenous yak and Tibetan sheep [29]. Currently, little is known whether the rumen microbiome is co-evolved tighter with ruminants living in highaltitude environments. In this study, we integrated rumen fermentation pro les, metataxonomics, metagenomics, and metabolomics of the rumen microbiome comparatively to explore how the rumen microbiome might help yak to adapt to its living conditions, especially its poor dietary conditions unique to the Qinghai-Tibet high plateau. Compared to the rumen microbiome of cattle, the yak rumen microbiome enriched cellulase and hemicellulase, and PL families, which can help improve ber degradation. On the other hand, yak had decreased amylase-containing GH families and CBM families, which might slow down starch degradation in the rumen and allow more starch to reach the small intestine. Additionally, the increased succinate-producing and utilizing bacterial species might promote propionate production. All the above are consistent with the ability of yak to live on the poor diet available on the high plateau.
Rumen microbiota is primarily responsible for the energy acquisition of ruminants, and many studies showed that differences in rumen microbiota might result in different energy e ciency [6,30]. We rst compared the rumen microbiota among the three ruminant species. Of the 21 phyla identi ed by metagenomics in this study, Bacteroidetes and Firmicutes, which are considered the important microbes in providing the energy required by ruminant animals [31], were the most predominant. Studies have shown that with the increase in dietary energy level and concentrate proportion, the relative abundance of Firmicutes could increase [32,33]. Similarly, the gut microbiota of obese people and mice caecum contained more Firmicutes and less Bacteroidetes [34,35], which agrees with that the relative abundance of Firmicutes can be positively correlated with dietary energy concentration. Compared with the study of Ahmad et al., [33] the relative abundance of Bacteroidetes in the three ruminant species was higher, but that of Firmicute was lower. The discrepancy might be partially attributable to the concentrate proportions (30%, day matter) in the diet of the cited study. We also detected many phyla, such as Proteobacteria and Verrucomicrobia, whose relative abundance was signi cantly higher in yak than in the other two ruminant species. Members of the phylum Proteobacteria are metabolically exible modifying their gene repertoire in response to changes in available substrate and energy sources [36, 37], and their relative abundance was signi cantly higher in the rumen of cows with high milk yield and milk protein content, suggesting that Proteobacteria is positively correlated with feed e ciency [6]. Verrucomicrobia contains species that are highly specialized degraders of fucoidans and other complex polysaccharides [38], and they may play an important role in polysaccharide degradation. In a previous transcriptomic study, Ruminococcus and Fibrobacter were shown to express most of the cellulase transcripts, they, together with Prevotella, probably expressed most of hemicellulose transcripts in the rumen [28]. The high abundance of Ruminococcus, Fibrobacter, and Prevotella and the activities of carboxymethyl cellulase and avicelase in the rumen of yak corroborate its greater ability to digest crude ber. Xue et al.
[6] compared the microbial composition between cows with high and low milk yield and milk protein content, and they found that the rumen species enriched in the high-yielding cows were mainly members of Prevotella, and P. oralis, P. buccae, P. albensis, and R. avefaciens were signi cantly higher in the high-yielding cows. We conclude that the above four bacterial species are positively correlated with feed e ciency. Interestingly, they were also signi cantly more predominant in the rumen of yak than the other two ruminant species. These species, probably among others, may help enable yak to adapt to its poor diet.
In the present study, yak was found to possess more hemicellulases (e.g., GH11, GH5, GH44, GH16, GH17) and rhamnogalacturonan endolyase (e.g., PL11and PL4-4) than cattle. Many CBMs associated with galactose binding (CBM61, CBM62, CBM51) and L-rhamnose binding (CBM67) were also enriched in yak. These results indicate that yak probably has a stronger activity in removing galactose from the main chain of hemicellulose. Rapid hemicellulose digestion can facilitate cellulose digestion in forage because the removal of hemicellulose can expose cellulose for microbial attachment and digestion. Although GH13, which contains the well-known amylases, did not differ among the three ruminant species, its subfamilies GH13-8 (starch branching enzyme) and GH13-12 and GH13-14 (pullulanase) and their associated CBM (CBM48) [39] were more abundant in cattle (Additional le 12: Table S6). Another starch binding CBM, CBM41 [40,41], was also enriched in cattle. The yak rumen microbiome probably has increased its ability to degrade forage not only by enhancing cellulose degradation (as evidenced by enriched GHs involved in cellulose hydrolysis) but also by enhancing hemicellulose and pectin degradation. On the other hand, the rumen microbiome of cattle has more abundance of GHs and CBM involved in starch digestion and binding, respectively. Even though both the yak and cattle used in the present study were fed the same diet, the rumen microbiome of yak probably has evolved to be more adapted to the forage that yak consumes, whereas the rumen microbiome of cattle, which are domesticated and has been fed diets containing starch, likely has evolved to better utilize starch. Yak produces milk of desirable characteristics (e.g., high milk fat, milk protein, milk total solids), but its milk yield is very low (about 1 kg per day) [42,43]. The lack of su cient dietary energy is one limiting factor. It will be interesting to test if the rumen microbiome of yak can evolve and gradually adapt to starch-rich diets over long-term adaptation.
It is well-known that forage typically increases acetate concentration in the rumen than a diet containing starch. In the present study, acetate molar proportion reached 75.1-77.9%, signi cantly higher than that reported in other studies [44][45][46], and the acetate: propionate ratio was also much higher than that reported by others [32,33]. Both were likely attributed to the forage-only diet used in the present study.
One previous study showed that rumen fermentation to propionate could increase energy e ciency [2]. Indeed, glucose fermentation to acetate, propionate, and butyrate can provide 62%, 109%, and 78% of the energy stored in glucose, respectively [47]. The gain of energy in propionate fermentation is because it incorporates hydrogen. This in turn reduces the amount of energy wasted via methane production from hydrogen. Two fermentation pathways are mainly responsible for the propionate production in the rumen: the acrylate pathway and succinic acid pathway [25], with the former using lactate as the input substrate and acryloyl-CoA as an intermediate, while the latter using oxaloacetate as the starting substrate and succinate as an intermediate. In this study, the increased relative abundance of succinateproducing and utilizing bacterial species seen in yak might be attributable to the higher propionate proportion detected in that species. On the other hand, because lactate-producing and utilizing bacterial genera did not differ in relative abundance among the three ruminant species, the acrylate pathway of propionate production might be similar across the three ruminant species. Isobutyrate, isovalerate, and 2methyl-butanoic are considered essential growth factors for most ber-degrading microorganisms [48][49][50]. The increased isobutyrate and isovalerate concentration (data was not shown) noted in the yak might facilitate the growth of some ber-degrading bacteria, which corroborates the increased relative abundance of the three well-known cellulolytic bacteria: R. avefaciens, R. albus, and F. succinogenes, in the yak. Branched-chain volatile fatty acids are used in synthesizing branched-chain amino acids [48]. The higher branched-chain volatile fatty acids concentration in the yak rumen may increase valine, leucine, and isoleucine biosynthesis. It should be noted that although all the three ruminant species had acetate as the major VFA, yak had the highest molar proportion of propionate but the lowest molar proportion of acetate. Thus, we speculated that in the absence of high-quality forage, the yak rumen microbiome is adapted to change the rumen fermentation and enhance energy harvest.
Our metabolomic results showed that cattle enriched styrene degradation and streptomycin biosynthesis, while yak and dzo enriched primary bile acid biosynthesis. Bile acids are important components of the bile, and they play an important role in fat metabolism, mainly in the intestinal and liver circulatory systems where they play a protective role through recirculation [51]. The presence of primary and secondary bile acids in the rumen suggests that the rumen microbiome may play an important role in the conversion of primary to secondary bile acids in ruminants [52]. Further studies are needed to explore the mechanisms by which bile acids circulate in ruminants.
Uridine 5'-monophosphate consists of phosphoric acid, ribose, and uracil, the latter of which is a nucleobase for RNA biosynthesis. The enriched pyrimidine metabolism in the cattle rumen might explain the increased uridine 5'-monophosphate therein. The enhanced nucleotide metabolism pathway in the cattle rumen microbiome may provide more substrates for replication and repair, translation, and transcription. The upregulated starch and sucrose metabolism in the cattle rumen might be attributable to the higher concentration of glucose-6-phosphate therein. Taken together, metabolite alterations were consistent with the metagenome results, and they jointly corroborate the adaption of yak to the high plateau environments.
As demonstrated in the literature [23,24,53], many factors, including feed, age, and management, can affect the gut microbiome composition and its functions in animals. Because yak is not domesticated, and the yak animals used in the present study were not con ned, we could not weigh each of the study animals for the exact bodyweight or track the speci c feed each animal consumed. Differences in body weight and feed consumption might be confounding factors that could have affected the results. Apparently, the yak rumen microbiome features, as determined by metataxonomics, enzymatic assays, metagenomics, and metabolomics, can help the host to adapt to its harsh environments, especially the poor diet available to it, but it remains to be determined how the body of yak affects its rumen microbiome. Future studies are required to explore the adaption relationship between high-altitude animals and the symbiotic microbiome. Such information will be useful to understand the adaptability of high-altitude animals to various environments and their domestication.

Conclusion
Our study revealed a clear difference in the composition, functions, and metabolome of the rumen microbiome among yak, dzo, and domesticated cattle. Although the relative abundance of the genes that code for amylases and their associated CBM families were lower in yak than in cattle, the yak rumen microbiome increased the abundance of cellulolytic bacteria and the genes encoding cellulase and hemicellulase, making yak better adapted to digest its roughage-based diets. Besides, the yak rumen microbiome had more succinateproducing and -utilizing bacterial species, supporting more pyruvate fermentation to propionate. These ndings can help better understand how the yak rumen microbiome aids in yak's adaption to its poor diet in high plateau environments, and provide a foundation for further studies on microbial roles in the physiology and health of agricultural animals.

Experiment design and sample collection
Ten Qaidam yellow cattle, dzo, and yak each with similar body weight (about 200 kg) and age (5-6 years old) were used in this study.
Because of the di culty to weigh each of the grazing animals, one herdsman estimated the bodyweight of each animal visually. All the study animals were raised in the Qinghai-Tibet plateau and had never received any supplementary feed. In summer, they were grazing in the mountain grassland, while they were grazing at the foot of the mountain in winter, where Kobresia myosuroides and Phragmites communis were the predominant pasture species. The sampling site was located in the central part of Qinghai Province, in the eastern part of the Qaidam Basin (36 ° 49 '-37 ° 20' N, 98 ° 87 '-99 ° 27' E). The annual average temperature in this area is 3.5℃, and the average altitude is 4,000 meters. Rumen uid samples were collected from each animal using an oral stomach tube and a pump, both of which were thoroughly cleaned using fresh warm water between sample collections. The rst 10-15 ml of sample from each animal was discarded to avoid contamination from saliva. Approximately 200 mL ruminal uid sample was collected from each animal, and half of it was immediately frozen in liquid nitrogen and then stored at -80℃ until analysis. The rest of each ruminal sample was ltered through four layers of cheesecloth, and the ltrate was used for pH measurement and then preserved at -20℃ for VFA analysis and brolytic enzyme activity assay.

Analysis of VFA and plant cell-degrading enzyme activity
The rumen uid samples were thawed at 4℃, and the solid particles and protein were removed according to the procedures of Li et al. [54]. The VFA concentrations were analyzed using gas chromatography (Agilent Technologies 7820A GC system, Santa Clara, CA) equipped with a 30 m × 0.25 mm × 0.33 μm fused silica column (AE-FFAP, Atech Technologies Co. Ltd., Shanghai, China).
One aliquot of each thawed rumen uid sample was centrifuged at 1,000× g for 10 min at 4℃. The supernatant was immediately collected and analyzed for the activity of carboxymethyl cellulase (CMCase), avicelase, xylanase, and acetylesterase using carboxymethyl cellulose, avicel, birchwood xylan, and p-nitrophenyl acetate, respectively, as the substrates. The enzyme assay reaction was incubated at 39℃ and pH 7.0 for 30 min except for the xylanase assay (for 15 min). The amounts of released reducing sugars were quanti ed using the dinitrosalicyclic acid colorimetry method [55], and the production of p-nitrophenol was determined at 415 nm [56, 57]. One unit of enzyme activity was de ned as the amount of enzyme releasing 1 μmol of reducing sugar (e.g., xylose or glucose equivalent) or p-nitrophenol per min per ml.

Metataxonomic analysis of rumen prokaryotes and fungi
Microbial DNA was extracted from each rumen sample using the E.Z.N.A.® soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer's protocols. The nal DNA concentration and purity were determined using a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scienti c, Wilmington, USA), followed by visual quality checking using agarose gel (1%) electrophoresis. Individual amplicon libraries were prepared for prokaryotes using PCR ampli cation of the V3-V4 hypervariable regions of the 16S rRNA gene with primers 338F The raw sequence reads were demultiplexed, quality-ltered using Trimmomatic, and then merged using FLASH. Operational taxonomic units (OTUs) were clustered (de novo) with a 97% similarity cutoff using UPARSE (version 7.1 http://drive5.com/uparse/) and chimeric sequences were identi ed and removed using UCHIME. The representative 16S rRNA gene sequences of the OTUs were taxonomically classi ed using the RDP Classi er algorithm (http://rdp.cme.msu.edu/) against the Silva 128 database at a con dence threshold of 70%. The ITS sequences were taxonomically assigned using the UNITE7.0 database (https://unite.ut.ee/).

Metagenomic sequencing and analysis
Each of the microbial DNA extracts was fragmented to an average size of about 300 bp using Covaris M220 (Gene Company Limited, China

Metabolomic analysis of rumen uid
One hundred µL of each rumen uid sample was subjected to metabolite extraction using 500 µL methanol: water (4:1, v/v) solution containing 2% L-2 chlorophenylalanine (as internal standard). Then the mixtures were vortexed for 10 s and centrifugation at 12,000 rcf at 4℃ for 20 min. The supernatant was carefully transferred to a glass-derived bottle and vacuum-dried. After 80 µL methoxy amine hydrochloride (15 mg/mL in pyridine) was added, the samples were vortex-mixed for 2 min and incubated at 37℃ for 90 min to carry out oximation reaction. Eighty µL of tri uoroacetamide reagent containing 1% trimethylchlorosilane and 20 µL n-hexane were added to each sample, and then all samples were vortex-mixed for 2 min and incubated at 70℃ for 60 min. The samples were cooled to room temperature and analyzed using gas chromatography and mass spectrometry (GC-MS).
The rumen metabolome was analyzed using gas chromatography (Agilent 7890A, Agilent Technologies, Inc., Santa Clara, CA, USA) coupled to an Agilent 5975C mass selective detector (Agilent, USA) with an inert electron impact ionization (EI) source and ionization voltage at 70 eV. Brie y, metabolites were separated with an HP-5MS capillary column (30 m × 0.25 mm × 0.25 um) using helium (99.999% purity) as the carrier gas at a constant ow rate (1 mL/ min). The GC column temperature was programmed to hold at 60℃ for 0.5 min, rise to 310℃ at a rate of 8℃/min, and then hold at 310ºC for 6 min. A QC sample was prepared by pulling an equal volume of each sample, and the QC sample and rumen uid samples were analyzed in the same manner. To assess the repeatability of the analysis, the QC sample was injected once every 10 rumen uid samples. The GC-MS data were processed using the MassHunter workstation Quantitative Analysis package (version v10.0.707.0) to extract raw peaks, lter and calibrate data baselines, align peaks, deconvolute, identify peaks, and integrate peak areas. The resulting matrixes detected in at least 80% of the samples were retained. After ltering, the missing values of the raw data were lled up by half of the compound minimum. The peak area was normalized in the data analysis. The internal standard was used for data quality control (reproducibility), and the metabolic features whose relative standard deviation (RSD) exceeded that of the QC by >30% were discarded. Mass spectra of these metabolic features were identi ed using the Fiehn database (https:// ehnlab.ucdavis.edu/projects/ ehnlib).
The metabolomic data were analyzed using PCA, and OPLS-DA was used to determine the global metabolic differences among the three ruminant species. Statistical signi cance among species was declared at a VIP value of >1 and a P-value of <0.05. P-values were estimated using paired Student's t-test for single-dimensional statistical analysis. A total of 27 differential metabolites among two of the three ruminant species were mapped into their biochemical pathways through metabolic enrichment and pathway analysis based on search against the KEGG database (http://www. genome.jp/kegg/) [63]. These metabolites were classi ed according to the pathways into which they were mapped or the functions that they could perform. Differential metabolites were cross-listed with the pathways in the KEGG database, and the top differential pathways were identi ed [64]. The relationship between different species and metabolites was visualized as a heat map using the "pheatmap" package in R (www.r-project.org).

Declarations Ethics approval
Animal care and experimental procedures were approved by the Institutional Animal Care and Use Committee (protocol number: NWAFAC1008) of the College of Animal Science and Technology of the Northwest A&F University (Yangling, Shaanxi, China).

Consent for publication
No applicable.

Availability of data and material
All sequencing data are available from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession number PRJNA744001, PRJNA744022 and PRJNA744415.

Competing interests
The authors declare that they have no con ict of interests.  Plots of principal coordinates analysis (PCoA) comparing the overall rumen microbiota among the three ruminant species. A (bacteria) and B (fungi) were based on Bray-Curtis dissimilarity determined by metataxonomics (n=9 per species), and C was based on Bray-Curtis dissimilarity determined by metagenomics (n=6 per species). The ellipses represent 95% of the samples belonging to each group.
Statistical signi cance of difference was tested using ANOSIM with 999 permutations.

Figure 2
Microbial species (identi ed in the rumen metagenome) signi cantly differed in relative abundance among the three ruminant species. The difference was tested using the Kruskal-Wallis H test, with post hoc test done using the Tukey-Kramer test (n=6 per species). Only the microbial species each with a relative abundance > 0.1% in all samples were shown. *P < 0.05, ** P < 0.01.  Differential KEGG pathways at level 2 (A), level 3 (B) among the three ruminant species as determined using LEfSe. Only the pathways that differed signi cantly (P < 0.05) between two of the three ruminant species with an LDA >2 are shown.

Figure 5
Metabolic pathways signi cantly differing in rumen metabolites between cattle and dzo (A), between cattle and yak (B), and between dzo and yak (C). The X-axis and Y-axis represent the pathway impact and pathway enrichment, respectively. The darker the color indicated the smaller the P-value from enrichment analysis, and the larger abscissa indicates greater impact from the pathway topology analysis, respectively (n = 6 per species).

Figure 6
Pearson correlations between signi cantly different rumen microbial species and microbial metabolites. The color and the color intensity of the squares correspond to the direction and strength of the correlation based on the scale to the right.