Abstract

   An intensive feeding system might improve the production cycle of yaks.
   However, how intensive feeding system contributes to yak growth is
   unclear. Here, multi-omics, including rumen metagenomics, rumen and
   plasma metabolomics, were performed to classify the regulatory
   mechanisms of intensive feeding system on yaks. Increased growth
   performance were observed. Rumen metagenomics revealed that
   Clostridium, Methanobrevibacter, Piromyces and Anaeromyces increased in
   the intensively fed yaks, contributing to amino acid and carbohydrate
   metabolism. The grazing yaks had more cellulolytic microbes. These
   microbiomes were correlated with the pathways of “Alanine aspartate and
   glutamate metabolism” and “Pyruvate metabolism”. Intensive feeding
   increased methane degradation functions, while grazing yaks had higher
   methyl metabolites associated with methane production. These rumen
   microbiomes and their metabolites resulted in changes in plasma
   metabolome, finally influencing yaks’ growth. Thus, an intensive
   feeding system altered the rumen microbiome and metabolism as well as
   host metabolism, resulting in improvements of yak growth.

   Subject terms: Metagenomics, Microbiome

Introduction

   The yak (Bos grunniens), inhabiting the Qinghai-Tibetan Plateau (the
   world’s highest plateau) in China, adapts to the harsh environment with
   low partial pressure of oxygen, and a high level of ultraviolet
   radiation^[34]1,[35]2. Traditionally, yaks graze yearly on high
   altitude pastures and experience inadequate nutrient sources in the
   long-cold winter. Thus, some specific anatomical and physiological
   traits in the yaks have been developed, such as a giant heart and lung,
   smaller blood cells, lack of hypoxic pulmonary vasoconstriction,
   increased ability to graze forage and high energy
   metabolism^[36]3,[37]4. However, some common features between yak and
   other ruminants are observed. For example, their rumen has a complex
   microbial community and can ferment fibrous plant materials into
   volatile fatty acids (VFAs) for the host requirement^[38]5. In the
   meantime, an extraordinary ability to degrade lignocellulose, produce
   VFAs and metabolize nitrogen was found in the rumen of yak compared to
   cattle under the same feeding system^[39]6,[40]7. Due to malnutrition
   and the extreme living environment, longer production cycle and poor
   meat quality in yaks are the major problems for yak industry which need
   to be addressed as consumer demands increase^[41]8. Although
   hybridization of yaks and other ruminant breeds has been conducted and
   could improve productivity, improvements of feeding strategy to
   increase purebred yak productivity may be the best approach for the
   long-term developments of the yak industry.

   An intensive feeding system provides the precisely designed diet to
   housed animals to maximize the partitioning of feed nutrients into
   productive outputs and minimize potentially polluting emissions^[42]9.
   Intensive feeding programs have been widely applied to beef and dairy
   cattle. For yaks, researchers have investigated the effects of the
   semi-intensive grassland-based systems and obtained improved results of
   growth performance, rumen fermentation and changes of gut microbiota
   via supplementation of concentrate diet^[43]10–[44]13. Furthermore,
   intensive feeding systems improved nutrient digestibility, body
   metabolism, meat quality and milk yield have been
   reported^[45]14–[46]17. As we know, the rumen as the most important
   digestive organ converts plant materials into VFAs and proteins for the
   host requirement via microbial fermentation^[47]18,[48]19. One study
   found that feedlot fattening changed the rumen fermentation and
   microbial diversity as well as meat quality in yaks^[49]20. Another
   research observed that warm-grazing and cold-indoor feeding regimes
   influenced yak rumen microbiome and greenhouse gas emissions compared
   to traditionally grazing yaks^[50]21. However, these studies use 16S
   rRNA sequencing technology, resulting in less knowledge in low
   abundance bacteria and other microbiota. Changes in the rumen
   microbiome and its functions in response to an intensive feeding system
   and its subsequent influences on the body metabolism in yaks are still
   unclear. Additionally, some questions need to be addressed, such as
   whether an intensive feeding system fits the yaks living in a
   high-altitude area, whether changes in yak rumen microbiome by
   intensive feeding are similar to other ruminants (e.g., beef and dairy
   cattle), and how intensive feeding affect the rumen microbiome
   (composition and functions), microbial metabolites, and the host
   metabolites in yaks.

   Therefore, in this study, we performed rumen metagenomics, rumen
   metabolomics, and plasma metabolomics to identify the impacts of an
   intensive feeding system on yaks compared to the traditional grazing
   system. This study will provide knowledge of the rumen
   microbiome-dependent mechanism to improve the host phenotypes, and
   evidence of the advantages of an intensive feeding system for the
   livestock industry.

Results

Intensive feeding system significantly altered the growth, metabolism and
physiology of yaks

   As designed, the nutrient level of the diet in the intensive feeding
   system was significantly greater than forage on the natural grassland.
   The crude protein and starch of the intensive diet were 16.7% and 36.2%
   compared to 8.5% and 3.0% in forage (Supplementary Table [51]1). In the
   meantime, the NDF and ADF were lower in intensive diet (20.1% and
   31.1%) than that in grazed forage (35.0% and 51.5%). Although the
   forage intake was not measured, significant improvements of the
   bodyweight of yaks consuming an intensive diet were observed as final
   body weight and ADG in the intensive group was higher compared to the
   grazing group (Supplementary Table [52]2).

   Next, serum biochemical parameters associated with the health status of
   animals and nutrient metabolism of diet in the body were measured in
   this study (Supplementary Table [53]3). Alkaline phosphatase (ALP),
   total bilirubin (TBil), direct bilirubin (Dbil), indirect bilirubin
   (Ibil) and glucose in the intensive group were significantly higher
   than those in the grazing group, while concentration of creatinine (Cr)
   and total cholesterol (TC) were significantly greater in the grazing
   yaks (P < 0.05). Other parameters, including the total protein (TP),
   albumin (ALB), globulin (GLB), urea and triglyceride (TG), were not
   different between the two groups (P > 0.05). These results indicated
   that the body metabolism, including carbohydrate and lipid, was
   influenced. In addition, regarding the hash living environment (e.g.,
   high altitude, low oxygen) of Zhongdian yaks, the blood physiological
   parameters were also determined (Supplementary Table [54]4). The
   results showed that red blood cell (RBC) count, hemoglobin (HGB), and
   mean corpuscular hemoglobin concentration (MCHC) were high in grazing
   groups. On the contrary, mean corpuscular volume (MCV), mean cellular
   hemoglobin (MCH), red cell distribution width variance coefficient
   (RDW-CV), platelet count (PLT) and mean platelet volume (MPV) were
   increased by the intensive feeding system.

Rumen metagenomics revealed the diversity and structure of the microbiome
associated with intensive feeding system

   Metagenome sequencing generated a total of 230.5 Gbp data, with
   1,026,889,026 reads (73,349,216 reads per sample). After de novo
   assembly, a total of 2,240,346 contigs were generated (the N50 length
   of 1,310 bp), with 160,024 per sample. The rumen metagenome across all
   samples consisted of 97.56% bacteria, 1.75% eukaryotes, 0.49% archaea,
   and 0.20% viruses (Supplementary Table [55]5). The rumen microbial
   domains, including bacteria and viruses, were increased by the
   intensive feeding system (P < 0.05), while the abundance of eukaryotes
   was greater in the grazing group. The principal coordinate analysis
   (PCoA) showed separations between the two groups based on rumen
   bacterial, eukaryotic and viral species (Fig. [56]1; Supplementary
   Table [57]6). In addition, no differences in the archaeal species were
   observed. At the same time, decreases in richness at the species level
   were observed among bacteria, archaea and eukaryote of yaks in the
   intensive group (Supplementary Figure [58]1).

Fig. 1. Intensive feeding system altered the rumen microbiome structure in
the yaks.

   [59]Fig. 1
   [60]Open in a new tab

   Bacterial (A), archaeal (B), eukaryotic (C) and viral (D) compositional
   structure of the rumen samples from grazing and intensively fed yaks
   based on species visualized using principal-coordinate analysis (PCoA).

Compositional profiles of the rumen microbiome in yaks altered by intensive
feeding system

   The dominant rumen bacterial phyla were Bacteroidetes (52.29%),
   Firmicutes (34.81%) and Proteobacteria (5.41%) across all samples
   (Supplementary Fig. [61]2A). In the meantime, we found that the
   intensive feeding system increased the abundance of Firmicutes (40.85%)
   and decreased Bacteroidetes (47.67%) compared to the grazing group
   (28.76% and 56.90%). At the genus level, dominant bacteria included
   Prevotella (43.96%), Bacteroides (7.94%), Clostridium (4.46%),
   Alistipes (3.53%), Fibrobacter (2.77%) and Ruminococcus (2.44%)
   (Supplementary Fig. [62]3A). At the species level, LEfSe algorithm
   identified 11 Clostridium sp., 3 Prevotella sp., 3 Ruminobacter sp., 2
   Succinatimonas sp., 1 Oscillibacter sp., Succiniclasticum ruminis,
   Ruminococcus bromii, and Faecalibacterium prausnitzii that were
   abundant in the intensive group (Fig. [63]2A). The grazing group had
   greater abundances of 5 Bacteroides sp., 3 Prevotella sp., 3 Alistipes
   sp., 3 Fibrobacter sp. and Pseudomonas fluorescens.

Fig. 2. The rumen microbiome distinguishing grazing and intensive groups.

   [64]Fig. 2
   [65]Open in a new tab

   Heatmap depicting the significantly different bacterial (A), archaeal
   (B), eukaryotic (C) and viral (D) species. E Correlation networks
   showed associations among microbiome species across kingdom. The edge
   colors (red: positive, blue: negative) are proportional to the
   correlation strength. Only strong (Spearman R of > 0.7 or <−0.7) and
   significant (P < 0.05) correlations are displayed. Significant
   differences were tested by linear discriminant analysis effect size
   (LEfSe) analysis, with linear discriminant analysis (LDA) score of > 2
   and P value of < 0.05.

   For archaea, the dominant phyla were Euryarchaeota (92.92%) and
   Candidatus Heimdallarchaeota (2.53%) across all samples (Supplementary
   Fig. [66]2B). The intensive group had a higher abundance of
   Euryarchaeota but lower abundance of Candidatus Heimdallarchaeota,
   Candidatus Bathyarchaeota and Candidatus Woesearchaeota. At the genus
   level, Methanobrevibacter (73.59%), Methanobacterium (3.49%),
   Methanosarcina (3.47%), Methanothrix (2.86%) and Methanomassiliicoccus
   (2.31%) were dominant ((Fig. [67]S3B). Except Methanobrevibacter was
   greater in the rumen of yaks of the intensive group, other dominant
   archaeal genera, including Methanosarcina, Methanothrix,
   Methanomassiliicoccus, Methanosphaera and Methanolobus, were higher in
   the rumen of grazing yaks. The rumen of intensively fed yaks was
   enriched by the species of Methanobrevibacter curvatus,
   Methanobrevibacter filiformis, Methanobrevibacter ccuticularis, and
   Lokiarchaeum sp. GC14_75 (Fig. [68]2B). However, 4 Euryarchaeota
   archaeon sp., 2 Methanolobus sp. (M. profundi and M. psychrophilus), 2
   Candidatus Heimdallarchaeota archaeon sp. and Methanosaeta harundinacea
   were over-represented in the rumen of grazing yaks.

   In addition, for the rumen Eukaryotic genera, we found that
   Neocallimastix, Piromyces and Anaeromyces were top taxa and had greater
   abundances in the intensive group (Supplementary Fig. [69]4A). The
   species, including Piromyces sp. E2, Anaeromyces robustus,
   Neocallimastix californiae, Arthrobotrys oligospora and Cadida auris,
   were greater in the intensive group (Fig. [70]2C). Abundances of
   Rozella allomycis, Rhizophagus irregularis, Basidiobolus meristosporus,
   Rhizoclosmatium gobosum, Allomyces macrogynus, Batrachochytrium
   dendrobatidis, Catenaria anguillulae, Conidiobolus coronatus, and
   Lobosporangium transversale were higher in the grazing group.

   Moreover, the intensive feeding system increased the abundances of
   rumen viral genera, including C5virus, Cp220virus and N4virus (31.83%,
   26.20% and 9.95%), compared to the grazing group (20.58%, 7.30% and
   1.38%) (Supplementary Fig. [71]4B). On the contrary, higher abundances
   of Hokovirus, Alpharetrovirus, Catovirus, Pepy6virus, Bc431virus and
   Cervidpoxvirus were observed in the rumen of grazing yaks. The species,
   including Geobacillus virus E3, 2 Campylobacter virus (CP21 and INN35)
   and Enterobacter phage EcP1, had higher abundances in the intensive
   group (Fig. [72]2D). In grazing yaks, their rumens were abundant with
   Hokovirus HKV1, Avian myeloblastosis virus, Catovirus CTV1, Deerpox
   virus W-848-83 and Klosneuvirus KNV1.

   Additionally, microbiome interactions at the species level were
   determined using network analysis (Fig. [73]2E). Complex correlations
   among bacterial, archaeal, eukaryotic and viral species were observed.
   Bacterial species were the bridge nodes to connect archaea, eukaryotes
   and viruses. For instance, Prevotella sp. (brevis, ruminicola)
   positively correlated with Methanomassiliicoccales archaeon RumEn M1,
   Methanosaeta harundinacea, Candidatus Heimdallarchaeota archaeon AB_125
   and Rhizopus microsporus, but Clostridium sp. CAG-1024 negatively
   correlated with Catovirus CTV1, and Hokovirus HKV1. Archaeal species
   interacted with each other. Methanobrevibacter olleyae,
   Methanobrevibacter ruminantium, Methanobrevibacter wolinii and
   Methanobrevibacter smithii with positive correlation formed a subunit,
   while Methanomassiliicoccus luminyensis, Methanobrevibacter millerae,
   Candidatus Methanomassiliicoccus ntestinalis and Methanobrevibacter sp.
   A27 shaped like another subunit. The dominant eukaryotic and viral
   species had complex interactions.

Intensive feeding system changed functional profiles of the rumen microbiome
in yaks

   The functions of the rumen microbiome were determined by the Kyoto
   Encyclopedia of Genes and Genomes (KEGG) profiles and genes encoding
   CAZymes. For KEGG profiles at the first level category, “Metabolism”
   (52.43%) and “Genetic information processing” (19.50%) were dominant
   (Supplementary Table [74]7). When the KEGG pathways were compared, the
   rumen of intensively fed yaks were significantly enriched with the
   pathway associated with amino acid metabolism, such as “Alanine
   aspartate and glutamate metabolism” (ko00250), “Cysteine and methionine
   metabolism” (ko00270), “Valine, leucine and isoleucine biosynthesis”
   (ko00290), “Arginine biosynthesis” (ko00220), “Histidine metabolism”
   (ko00340), and associated with carbohydrate metabolism of “Pyruvate
   metabolism” (ko00620), “Glycolysis Gluconeogenesis” (ko00010), “Pentose
   phosphate pathway” (ko00030), “Citrate (TCA) cycle” (ko00020),
   “Butanoate metabolism” (ko00650), “Propanoate metabolism” (ko00640)
   (LDA > 2 and P < 0.05; Fig. [75]3). Moreover, other pathways associated
   with energy metabolism (e.g., “Methane metabolism” (ko00680), “Carbon
   fixation in photosynthetic organisms”, “Nitrogen metabolism”) were also
   higher in the intensive group. However, the grazing group was enriched
   with some signal transduction pathways, such as “AMPK signaling
   pathway” (ko04152), “FoxO signaling pathway” (ko04068) and “PI3K-Akt
   signaling pathway” (ko04151). The intensive group was enriched by 11 of
   the KEGG modules related to nucleotide and amino acid metabolism,
   including M00053, M00021, and 21 modules associated with carbohydrate
   metabolism, such as M00173, M00003, M00620, M00004, M00377, M00308,
   M00009, and M00374 (Supplementary Fig. [76]5). Overall, the intensive
   feeding system increased the function of amino acids and carbohydrate
   metabolism (Fig. [77]4).

Fig. 3. Differential KEGG functions of rumen microbiota between the grazing
and intensively fed yaks.

   [78]Fig. 3
   [79]Open in a new tab

   Significant differences were tested by linear discriminant analysis
   effect size (LEfSe) analysis, with linear discriminant analysis (LDA)
   score of > 2 and P value of < 0.05.

Fig. 4. Microbial functions involved in carbohydrate metabolism and energy
metabolism in the rumen of the intensively fed yaks.

   [80]Fig. 4
   [81]Open in a new tab

   Metabolic pathways involved in VFA biosynthesis and methanogenesis were
   displayed.

   For CAZyme profiles, classes of glycoside hydrolases (GHs) (64.68%),
   glycosyltransferase (GTs) (15.60%), carbohydrate-binding modules (CBMs)
   (9.87%), carbohydrate esterases (CEs) (7.77%) were dominant in all
   samples (Supplementary Fig. [82]6). At the family level, 15 GHs, 7 GTs,
   4 CMBs and 2 CEs were identified as signatures for the intensive group,
   while 11 GHs, 3 GTs and 8 CBMs were enriched in the grazing group
   (Supplementary Fig. [83]7). Moreover, the rumen CAZyme activities from
   the intensive group mainly were under GH13 associated with starch
   utilization, such as branching enzyme (EC2.4.1.18),
   cyclomaltodextrinase (EC3.2.1.54), alpha-amylase (EC3.2.1.1),
   isoamylase (EC3.2.1.68), amylosucrase (EC2.4.1.4), starch synthase
   (EC2.4.99.16), and neopullulanase (EC3.2.1.135) (Supplementary Fig.
   [84]8). The rumen of grazing yaks had a greater abundance of GH3
   enzymes to degrade lignocelluloses, such as alpha-L-arabinofuranosidase
   (EC3.2.1.55), beta-glucosidase (EC3.2.1.21) and beta-xylosidase
   (EC3.2.1.37).

Associations between rumen microbiome and its functions

   As protein and carbohydrate (e.g., fiber and starch) contents are two
   significant determining nutrient factors of the diet in the intensive
   feeding and grazing system, we further focused on the functions of
   amino acid and carbohydrate metabolism in the rumen microbiome. Two
   most important pathways involved in amino acid and carbohydrate
   metabolism were “Alanine aspartate and glutamate metabolism” (ko00250)
   and “Pyruvate metabolism” (ko00620) enriched in the rumen of intensive
   group yaks (LDA > 2, P < 0.05; Fig. [85]3). The abundances of genes
   encoding enzymes involved in these two pathways were also enriched in
   the rumen of intensive yaks, such as EC 6.3.5.5, EC 6.3.1.2 associated
   with ko00250 and EC1.2.7.1, EC4.2.1.2 within ko00620 (Supplementary
   Fig. [86]9). In addition, “Methane metabolism” (ko00680) and its
   related enzymes (e.g., EC1.2.7.1, EC 1.1.1.95 and EC 1.1.1.399) were
   higher in the intensive group.

   A network with Spearman’s rank correlation between species and those
   three pathways was then built to estimate the associations between the
   rumen microbiome and its functions (Fig. [87]5). A total of 41 species
   (17 positive and 24 negative correlation) showed significant
   correlation with amino acid and carbohydrate metabolism pathways as
   well as methane metabolism (R > 0.70 and P < 0.05). Bacterial species
   played key roles in both positive and negative correlations. The
   species of Clostridium sp. showed the strongest positive relationships,
   while Bacteroides sp. were prominent members among negative
   associations. Additionally, archaeal species, including
   Methanobrevibacter curvatus, Methanobrevibacter filiformis, and
   Methanobrevibacter ccuticularis, with greater abundance in intensively
   fed yaks’ rumen were positively correlated with these three pathways.
   Correspondingly, signature archaea for the intensive group (e.g.,
   Methanosaeta harundinacea, Candidatus Heimdallarchaeota archaeon
   AB-125) negatively related with “Alanine aspartate and glutamate
   metabolism”, “Pyruvate metabolism” and “Methane metabolism”.
   Interestingly, species of Fibrobacter sp. UWB4 and Methanolobus
   profundi had a negative correlation with “Methane metabolism”.

Fig. 5. Correlation networks showed associations between significantly
different microbiome species and the three most important pathways.

   [88]Fig. 5
   [89]Open in a new tab

   The edge colors (red: positive, blue: negative) are proportional to the
   correlation strength. Only strong (Spearman R of > 0.7 or < − 0.7) and
   significant (P < 0.05) correlations were displayed.

Rumen fermentation parameters and metabolome associated with intensive
feeding system

   The changes in microbial fermentation products affected by the
   intensive feeding system were estimated. We observed that the intensive
   group had significantly higher concentrations of valeric acid and
   isovaleric acid (Supplementary Table [90]8). The higher concentration
   of ammoniacal nitrogen was higher in the grazing group. Thus, the mode
   of rumen fermentation was altered by the feeding system.

   We then analyzed the rumen metabolites using untargeted metabolomics.
   Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA) showed
   remarkable differences between the intensive and grazing groups
   (Supplementary Fig. [91]10A). A total of 256 compounds were identified
   in the rumen metabolomics. After t-test and variable importance in
   projection (VIP) filtering for the relative concentrations of rumen
   metabolites, 42 metabolites were significantly higher in the rumen of
   intensive feeding yaks, such as L-Glutamic acid, 3-Methyladenine,
   Nicotinic acid, Gibberellin A97, Pyroglutamic acid, Valyl-Isoleucine,
   Pyridoxamine and Meta-Tyrosine, while Proline betaine,
   N-Acetyldopamine, 1-nitroaphthalene and N-Acetyldopamine were abundant
   in the grazing group (Supplementary Fig. [92]10B, Fig. [93]6A).
   Additionally, Metaorigins analysis was performed to track where the
   signature metabolites come from the host, microbiome, or both
   (Supplementary Fig. [94]11). We found that the signature metabolites
   differentiating the two groups were from both host and rumen microbiome
   and related to amino acids metabolism, such as “Histidine metabolism”,
   “Alanine, aspartate and glutamate metabolism”, “Arginine and proline
   metabolism”, and “Cysteine and methionine metabolism”. Then, using
   Metabolic pathway analysis (MetPA), a total of 11 pathways of these
   metabolites were identified as significantly different pathways
   (Benjamini-Hochberg false discovery rate (FDR) < 0.01, pathway impact >
   0.1), such as “D-Glutamine and D-glutamate metabolism”, “Alanine,
   aspartate and glutamate metabolism”, “Histidine metabolism”, “Arginine
   and proline metabolism”, “Vitamin B6 metabolism”, “Glutathione
   metabolism” and “Arachidonic acid metabolism” (Fig. [95]6B).

Fig. 6. Rumen metabolome of the grazing and intensively fed yaks.

   [96]Fig. 6
   [97]Open in a new tab

   A The Euclidean distance matrix for the quantitative values of
   differential metabolites was calculated, and a heatmap was selected to
   display the metabolites between the grazing and intensively fed yaks
   using a complete linkage method. B Pathway enrichment analysis
   performed using the significantly different rumen metabolites between
   grazing and intensively fed yaks.

Plasma metabolome influenced by the intensive feeding system

   We also examined plasma metabolomes associated with the intensive
   feeding system. A significant shift of the structure in the plasma
   metabolome between the intensive and grazing groups was observed
   (Supplementary Figure [98]12). The comparison analysis revealed that
   the relative concentrations of 26 metabolites (e.g., D-Proline,
   L-Phenylalanine, Taurine, Graveolinine and Acetone cyanohydrin) were
   significantly higher in the blood of intensively fed yaks, and the
   relative concentrations of 118 metabolites (e.g., L-Carnitine, Proline
   betaine, Creatinine, Betaine, and L-Arginine) were significantly higher
   in the plasma of grazing yaks (P < 0.05, VIP > 1) (Supplementary Fig.
   [99]12, Fig. [100]7A). Using MetPA, the enrichment of 20 pathways was
   significantly different between the two groups, including “Taurine and
   hypotaurine metabolism”, “Phynylalanine, tyrosine and tryptophan
   biosynthesis”, “Phenylalanine metabolism”, “Glycerophospholipid
   metabolism”, “Alanine, aspartate and glutamate metabolism” and “Primary
   bile acid biosynthesis” (FDR < 0.01, pathway impact > 0.1) (Fig.
   [101]7B).

Fig. 7. Plasma metabolome of grazing and intensively fed yaks.

   [102]Fig. 7
   [103]Open in a new tab

   A The Euclidean distance matrix for the quantitative values of
   differential metabolites between grazing and intensive yaks was
   calculated, and the abundances were changed to logscale for significant
   metabolites visualization of the radar chart. B Pathway enrichment
   analysis was performed using the significantly different rumen
   metabolites between Grazing and Intensive groups.

Relationships between the rumen microbiome, rumen metabolome and plasma
metabolome revealing the mechanism of the intensive feeding system to affect
yaks

   A comprehensive analysis of the rumen microbiome, metabolome and plasma
   metabolome could reveal how the intensive feeding system affects the
   microbiome and metabolism of yaks, resulting in an understanding of the
   overview of the mechanism of the intensive feeding system. Network
   analysis was performed to determine their associations (Fig. [104]8).
   We found that rumen bacteria were tightly associated with rumen
   metabolites, although other species also contributed to rumen
   metabolome. The critical microbiotas for each group were associated
   with their specific metabolites, which formed two modules with negative
   correlation (Fig. [105]8, Supplementary Fig. [106]13A). For example,
   Clostridium sp. as the signature microbiota for the intensive group was
   positively correlated with L-Glutamic acid, Pyroglutamic acid, Nitrite
   and N-Acetylputrescine that were abundant in the intensive group.
   Similarly, species, including Alistipes sp. and Bacteroide sp.,
   positively correlated with L-Carnitine and Proline betaine abundant in
   the grazing group. Moreover, archaeal species showed similar patterns,
   such as Methanobrevibacter filiformis, Methanobrevibacter ccuticularis,
   M. profundi and M. psychrophilus. Next, similar relationships between
   rumen microbiome and plasma metabolites were found (Fig. [107]8,
   Supplementary Fig. [108]13B). Eukaryotes and viruses had more
   correlations with other microbial species and rumen metabolites. In
   addition, associations between rumen and plasma metabolites existed
   significant correlations were found, and signature metabolites for the
   intensive and grazing groups shaped two subunits, respectively (Fig.
   [109]8, Supplementary Fig. [110]13C). Notably, shared metabolites
   (L-Carnitine and Proline betaine) with higher abundances in both rumen
   and plasma metabolomics from grazing yaks were positively correlated
   (Supplementary Fig. [111]14). Moreover, the common predicted KEGG
   pathways of rumen and plasma metabolites were observed, including
   “Alanine aspartate and glutamate metabolism”, “Histidine metabolism”,
   “Arginine and proline metabolism”, and “D-Glutamine and D-glutamate
   metabolism” (Supplementary Fig. [112]15). Thus, the rumen microbiome
   and metabolites interacted with each other and influenced the host
   plasma metabolism.

Fig. 8. Interactions between rumen metagenome, metabolome, and plasma
metabolome in the yaks.

   [113]Fig. 8
   [114]Open in a new tab

   Spearman’s rank correlations were calculated, and only strong (Spearman
   R of > 0.7 or < − 0.7) and significant (P < 0.05) correlations were
   shown.

Discussion

   With increasing demands of animal-based products, intensive feeding
   practice for finishing cattle, especially beef cattle, has been widely
   adopted^[115]9,[116]22. However, yaks, the vital livestock species of
   the Qinghai-Tibetan Plateau, are still reared in a wild or semi-wild
   program, resulting in a lower growth rate and longer production cycle.
   To address this question, based on traditional feeding practices
   (grazing pasture yearly), this study house-reared finishing yaks under
   an intensive feeding system via giving a designed high-grain diet. A
   significant increase in growth performance of intensively fed yaks was
   observed as high grain diets with greater protein and starch were
   consumed. Serum biochemical and physiological parameters were also
   altered. Based on the rumen metagenome, the community composition of
   the microbiome and its functions changed by the intensive feeding
   system. In the meantime, we observed that the changes in rumen
   fermentation mode and metabolome were correlated with microbial
   alterations. Moreover, signature plasma metabolites for the intensive
   group were characterized. Additionally, identifying complex
   relationships among rumen microbiome, rumen and plasma metabolites
   provided insight how intensive feeding contributed to the host
   physiology and metabolism, which benefits into yak production.

   The positive effects of the intensive feeding system on the growth
   performance of yaks were not surprising due to their diet composition
   with high crude protein and starch but less fibers. Similar results
   were also reported in yaks^[117]3, beef cattle^[118]22, sheep^[119]23
   and goat^[120]5. The changes in the serum biochemical parameters
   reflect the high nutrient utilization efficiency in the intensive
   group. Notably, for yaks, serum physiological parameters, such as red
   blood cell count and hemoglobin, associated with high altitude and
   lower pressure of oxygen were altered. Red blood and hemoglobin were
   more significant in grazing yaks, which is because grazing yaks
   exercise more during grazing. On the other hand, intensive feeding
   systems changed the oxygen delivery model of yaks. Unfortunately, it is
   unknown whether this condition is suitable for yaks living in
   high-altitude areas or not. Long-term observation of this aspect should
   be investigated in future studies.

   Rumen microbiome plays a critical role in nutrient utilization, methane
   production, host growth and food security^[121]24,[122]25. Most
   previous studies investigated yak rumen microbiota or intensive feeding
   system using 16S rRNA sequencing^[123]26–[124]28, which restricts our
   understanding of low abundance bacteria and other microbial species.
   Based on metagenomics, this study found that the intensive feeding
   system significantly altered the rumen microbiome, including bacteria,
   archaea, eukaryotes, and viruses. This might be mostly due to the
   nutritional factors (proteins and carbohydrate) of the diet of the
   intensive feeding system^[125]29. A previous study found that crude
   protein level in the supplementary diet resulted in the alteration of
   rumen bacteria and fermentation in Jersey-yak^[126]11, which reflected
   the effects of crude protein levels in the current study. Our previous
   research also confirmed that additional high-grain diet supplementation
   could influence the rumen microbiota^[127]5. Bacteria as critical
   players in diet degradation and fermentation was the major kingdom
   influenced by the intensive feeding system, which is similar to
   previous studies^[128]30,[129]31. Clostridium sp. involved in amino
   acid and carbohydrate metabolism were greater in intensive groups,
   indicating their ability to utilize protein and starch^[130]32.
   Prevotella ruminicola and brevis as well as Bacteroides sp. were more
   abundant in the rumen from yak grazing on pastures. This might be due
   to their preferred roles in fiber-based polysaccharides
   digest^[131]33,[132]34. Fibrobacter sp. with cellulolytic properties
   were higher in grazing groups, which reveals the diet effects
   further^[133]35. Furthermore, high abundance of Bacteroidetes was found
   in yak rumen community. A previous study observed that grazing yaks had
   higher abundance of Bacteroidetes than that of feedlot Simmental
   cattle^[134]20. This noted that yak rumen microbiome was different and
   had its specific characteristics compared to other ruminants. In the
   meantime, shared rumen microbial species among ruminants reflected
   their commonalities. Additionally, for archaea, Methanobrevibacter was
   the only dominant genus in intensive group, and its major species
   members included M. curvatus, filiformis, and cuticularis. However, the
   grazing group had too many other archaeal species. The distributions of
   these methanogenic species suggest that high-grain diet may be
   associated with slight decreases in methane metabolism and emissions
   via reducing the number of archaeal species^[135]36. Although the
   abundance and functions of eukaryotes and viruses in the rumen are
   relatively small and largely unknown, the current study did observe
   several species were affected by the intensive feeding, such as
   Neocallimastix californiae and Geobacillus virus E3. Recent studies
   revealed that eukaryotic and viral members are associated with
   carbohydrate enzymes or feed digestibility^[136]37,[137]38, which is
   also found in this study. Another study reported that Neocallimastix
   californiae co-cultured with Clostridium acetobutylicum benefited for
   butyrate production from cellulose^[138]39. All these results of
   microbiome changes corresponded to the consumption of high grain diet.
   Further studies should be conducted to deeply explore functions and
   associations of eukaryotes and viruses with bacteria.

   Corresponding to the changes in the rumen microbiome, we observed that
   the intensive feeding system greatly impacted carbohydrate metabolism.
   Several pathways enriched in the intensive group based on rumen
   metagenomic and metabolomics, including “Glycolysis”, “Pentose
   phosphate pathway”, and “Pyruvate metabolism”, degrade and convert
   starch into pyruvate and VFAs, revealing that high starch intake
   altered the yak rumen fermentation and microbiome. Over-representation
   of starch-related CAZyme (e.g., branching enzyme (EC2.4.1.18),
   alpha-amylase (EC3.2.1.1) and cyclomaltodextrinase (EC3.2.1.54)) also
   confirmed this result. Moreover, high abundance of L-Carnitine from the
   rumen and plasma metabolomics reflected the improvement of carbohydrate
   metabolism in intensively fed yaks since L-Carnitine was correlated
   with fat synthesis and accumulation^[139]40.

   Interestingly, for the grazing yaks, the KEGG pathway of its rumen
   microbiome contained more signal transduction, and more fiber
   degradation-related enzymes were found. This study was conducted in
   winter, and grazing yaks responded to insufficient supply and high
   fiber proportion of withered and yellow forage^[140]41,[141]42. Thus,
   higher VFA concentrations found in yaks under the intensive feeding
   system revealed that dietary starch consumption influenced rumen
   microbiome and metabolism, which finally benefits yak growth.

   Diet protein is another essential and critical nutrient for animals.
   The emission of nitrogenous substances (e.g., ammonia) is related to
   environmental challenges and feed conversion efficiency^[142]43. Thus,
   the rumen microbial metabolic and metabolomic pathways involved in
   amino acid degradation and biosynthesis were deeply investigated in
   this study. For rumen metagenomics, many previous studies focused on
   carbohydrate metabolism and paid less attention to amino acid
   metabolism^[143]11,[144]31. This study observed consistent changes
   associated with amino acid metabolism based on rumen metagenomics and
   metabolomics. For instance, the intensive feeding system increased the
   KEGG pathways associated with amino acid metabolism, such as “Alanine
   aspartate and glutamate metabolism”, “Cysteine and methionine
   metabolism” and “Histidine metabolism”, which indicated that yaks under
   the intensive feeding system might have a higher efficient of protein
   utilization Our results also found that the signature microbiomes for
   the intensive group were correlated with “Alanine aspartate and
   glutamate metabolism”. The microbiome correlated with amino acid
   metabolism was also reported in a previous study^[145]44. Clostridium
   sp. as the intensive group classifier might be hydrogen acceptors in
   amino acid fermentation^[146]45. Moreover, the hydrogenotrophic species
   of Methanobrevibacter with higher abundances in the intensive group may
   branch amino acids into nitrogen via mutualistic cooperation with
   bacteria (e.g., Prevotella)^[147]46. Other eukaryotic species,
   including Piromyces sp. E2 and Anaeromyces robustus, were associated
   with amino acid and carbohydrate degradation^[148]39,[149]47,[150]48.
   Additionally, amino acid metabolism is tightly related to carbohydrates
   as carbon skeletons derived from deamination could yield VFAs^[151]49.
   Similarly, our results found that the microbiomes associated with
   “Alanine aspartate and glutamate metabolism” were also related to
   “Pyruvate metabolism”. However, the cellulolytic microbes (Bacteroides
   sp. and Fibrobacter sp.) that do not produce hydrogen negatively
   interacted with amino acid metabolism, which might be due to
   competition relationships between cellulolytic microorganisms and amino
   acid producers^[152]50. Shared Proline betaine of rumen and plasma
   metabolomics was higher in grazing yaks, indicating that amino acid
   metabolism in the rumen is associated with host metabolism. Therefore,
   hydrogen related microbiotas are critical to rumen microbial amino acid
   biosynthesis and metabolism, and the intensive feeding system in yaks
   could promote their amino acid metabolism in the rumen and plasma.

   Ruminal methane production is defined as energy loss and related to
   greenhouse effect. Ruminant researchers found that diet or feeding
   strategy is associated with reduction of methane emissions. A previous
   report found greater CH[4] production for cattle consuming high-fiber
   diets than for cattle fed high-grain diets^[153]36. In the current
   study, the intensive feeding system might benefit methane emission
   reduction in yaks compared to grazing feeding system since more
   high-grain diet was consumed in yaks. In this study, only three
   Methanobrevibacter species had higher relative abundances in the rumen
   of yaks of the intensive group, and abundant “Methane metabolism” of
   rumen microbiome was also found. However, the rumen of yaks in grazing
   groups had more archaeal species with relative high abundances. Two
   major routes for methanogenesis in the rumen include the
   hydrogenotrophic pathway and methyl groups. Many archaeal, eukaryotic
   and bacterial species producing hydrogen contribute to methane
   production, such as Methanobrevibacter sp. and fibrolytic
   bacteria^[154]51. Moreover, methyl groups contained methylamines
   derived from proline betaine that were abundant in the rumen and plasma
   metabolomics of grazing yaks. Thus, the intensive feeding system for
   yaks might improve its methane emissions via manipulation of the rumen
   hydrogenic microbiome and methyl metabolites. Future investigations
   should measure the methane emission of the yaks to support this
   speculation.

   In ruminants, although previous studies have started to investigate the
   relationship between rumen microbiome, rumen metabolome and blood
   metabolome^[155]28,[156]31,[157]52, the detailed information how yaks
   responded the intensive feeding system is still unclear. A recent study
   investigated the effects of high-grain diet on the 16S rRNA of the
   rumen bacteria and their metabolome in adult yaks and found it altered
   the microbiota and metabolites^[158]28. Another study found that
   house-feeding altered the profile of rumen bacteria based on 16S rRNA
   sequencing and serum metabolome in Jersey cows^[159]27. The
   house-fattening feeding system affecting rumen bacteria and metabolites
   in yaks was estimated^[160]52,[161]53. In our study, rumen metagenome
   was determined to excavate more microbial species and their functions,
   and the associations between multi-omics data and the intensive feeding
   system were more deeply analyzed. Although individual variation was
   observed, our results profoundly illustrated the mechanisms of how
   intensive feeding systems affected the rumen microbiome and metabolites
   and subsequently host metabolites to improve yak growth performance.
   Archaea, eukaryotes, and viruses in yak rumen were incorporated into
   analyses, which unveils their specific roles in nutrient metabolism and
   interactions with rumen bacteria. One of shortness of this study might
   be lack of validation trial of signature microbiotas, which will be
   performed in our future research. Overall, the intensive feeding system
   influenced multiple aspects of yaks in our results and the intensive
   feeding system benefits the yak industry.

   Our study confirmed that the intensive feeding system significantly
   improved growth performance, rumen fermentation and serum parameters
   compared to traditional grazing practices in yaks. Based on multi-omics
   results, we identified the signature features of the rumen microbiome,
   microbial functions, metabolites, and their relationships with plasma
   metabolome which benefits yak growth. High grain diet led to increases
   of rumen Clostridium sp., Methanobrevibacter sp. (M. curvatus, M.
   filiformis, M. ccuticularis) and Piromyces sp. E2 that contributed to
   amino acid and carbohydrate metabolism but decreased methane
   production. The rumen of grazing yaks with high fibrolytic
   microorganisms and methylamines might be the major reason for high
   methane production and lower feed utilization efficiency. The levels of
   small molecules and their end products in the rumen were absorbed and
   transported by the host. In summary, the intensive feeding system may
   improve yak growth and methane emission via rumen microbiome and host
   metabolites. This study provides evidence that the intensive feeding
   system could be adopted for the yak industry with increased production
   efficiency and decreased greenhouse gas emissions.

Methods

Animals and sampling

   The animal trial was conducted at the Tiancheng Lun Zhu Agricultural
   Products Development Co., Ltd. located in the north of Shangri-La
   County, Yunnan Province, China. The farm belonging to the subtropical
   and cold temperate monsoon climate was at an altitude of approximately
   3600 meters.

   A total of 50 healthy Zhongdian yaks aged 3 to 4 years old were
   selected and divided into two treatments. One treatment group was
   traditionally grazed all day on the natural ranch (Grazing group),
   while the other group was reared under an intensive feeding system
   (Intensive group) that housed and fed a total mixed ration diet to
   yaks. The trial lasted 60 days. At the end of the trial, fourteen yaks
   were randomly selected from two treatments to collect samples for
   downstream analysis (e.g., blood parameters, rumen fermentation
   parameters and multi-omics), with seven bulls in each group. The
   ingredients and the chemical composition of the diets are shown in
   Supplementary Table [162]1. A total mixed ration comprised whole plant
   corn silage, beef cattle concentrate and crushed corn. The grassland
   where yaks grazed contained forage including Blysmus sinocompressus,
   Carex tristachya, Potentilla fulgens, Spiranthessinensis, Poa
   pratensis, and Euphorbia jolkinii. The treatment yaks were fed twice at
   9:00 a.m. and 4:00 p.m. every day, and the diet intake was recorded
   daily. The yak’s body weight was measured before early feeding at the
   beginning and end of the trial using fifty yaks, and average daily gain
   (ADG) was calculated. Specifically, yaks grazing on the natural pasture
   were driven to a holding pen at the night for the next day sampling,
   which guaranteed the same sampling time of two treatment groups. All
   animals used in the current study were healthy and did not receive any
   recorded therapeutic or prophylactic antibiotic treatments.

   An oral stomach tube was used for the collection of rumen content
   samples. In brief, at 8 am in the morning, a tube was inserted into the
   rumen, and a vacuum sampler was used to pump rumen fluid. To avoid
   saliva contamination, the first 30 mL of rumen digesta sample was
   discarded. For each animal, 30 ml of rumen fluid was collected, and
   divided into three parts by placing into 10 ml polypropylene tubes.
   Then, all samples were snap-frozen in liquid nitrogen, and brought back
   to the laboratory for long-term storage at −80 °C. The jugular blood
   samples were also collected at the end of the trial. Before feeding in
   the morning, blood samples were collected into three 10 mL tubes. One
   tube without anticoagulant was rested for 1 h and centrifuged (3000 ×
   g, 10 min, 4 °C), and then the aliquots of serum were frozen
   immediately at −20 °C until later chemical parameters analysis. Another
   tube containing heparin sodium anticoagulant was used for blood
   physiological parameter measurement, whereas the remained blood tube
   containing heparin sodium anticoagulant was centrifuged (3000 × g,
   10 min, room temperature), and then the blood plasma was collected and
   stored at −80 °C for metabolomics.

Chemical analysis

   Diet nutrient measurement: The samples of forage and total mixed ration
   were ground to pass through a 1-mm sieve and dried in an oven at 105 °C
   for 2 h to measure the dry matter (DM) content. The ash content,
   neutral detergent fiber (NDF), acid detergent fiber (ADF), calcium, and
   total phosphorus were determined according to previous methods^[163]54.
   Total nitrogen of diet samples was measured by an Auto‐Analyzer after
   Kjeldahl digestion and crude protein (CP) of the diet was calculated as
   6.25 × N. Ether extract (EE) was determined using a reflux system
   (ANKOM XT15, America) with petroleum ether at 90 °C for 1 h. Starch was
   measured by the enzymatic method, involving incubation of the
   autoclaved sample with amyloglucosidase and spectrophotometric
   measurement of released glucose.

   VFAs determination: The rumen fluid samples were thawed at 4 °C and
   then centrifuged at 2500 × g at room temperature. Next, 1 ml of the
   supernatant per sample was transferred into a 1.5 ml centrifuge tube
   which contained 0.2 ml of metaphosphoric acid solution (25% w/v). Then,
   the mixture was centrifuged at 10,000 × g at 4 °C after it was placed
   in an ice water set at 4 °C for 30 min. The collected supernatant was
   transferred to a solution containing an internal standard
   (2-Ethylbutyric acid, 2 EB) and a concentration of 25% metaphosphoric
   acid. The VFA concentration was detected using gas chromatography
   (GC–6800, Beijing Beifen Tianpu Instrument Technology, Co., Ltd.,
   China).

   Serum parameters analyses: The serum analysis was performed at Beijing
   Jinhai Keyu Biotechnology Development Co., Ltd. (Beijing, China). The
   concentration of serum biochemical parameters of alkaline phosphatase
   (ALP), total bilirubin (TBil), Direct Bilirubin (DBil), Indirect
   Bilirubin (IBil), Total Protein (TP), albumin (ALB), globulin (GLB),
   creatinine (Cr); triglyceride (TG) and total cholesterol (TC) was
   determined using A Kehua ZY KHB-1280 automatic biochemical analyzer
   (Shanghai Kehua Biological Engineering Co., Ltd., Shanghai, China)
   according to the manufacturer’s instructions of the standard commercial
   kits. Blood physiological indicators, including Red Blood Cell (RBC)
   Count, Hemoglobin (HGB), hematocrit (HCT), Mean corpuscular volume
   (MCV), mean corpuscular hemoglobin content (MCH), mean corpuscular
   hemoglobin concentration (MCHC), coefficient of variation (RDW),
   platelet number (PLT), mean platelet volume (MPV) and platelet
   distribution width (PDW), were measured using automatic blood analyzer
   (BC-2800Vet, Tokyo, Japan).

Shotgun metagenomic sequencing of rumen microbiome

   Rumen samples were subjected to metagenomic sequencing. All
   experimental procedures were executed in a stringent clean and sterile
   class II B2 biosafety chamber to lower the contamination. Total genomic
   DNA was extracted using a DNeasy PowerSoil Kit (Qiagen, Hilden,
   Germany) according to the manufacturer’s instructions. Negative
   controls during extraction were included. DNA integrity was estimated
   using 1% agarose gel electrophoresis. A Qubit 2.0 Fluorimeter
   (Invitrogen, Carlsbad, California, USA) was employed to measure the DNA
   concentration of samples. The metagenomic library was built using a
   TruSeq DNA PCR-Free Library Preparation Kit (Illumina, San Diego,
   California, USA). Metagenomic sequencing was conducted using an
   Illumina HiSeq 4000 platform with 150 bp paired end reads at the
   Beijing Headquarters of Novogene Biotech Co., Ltd (Beijing, China).

   The metagenomic raw data obtained from the Illumina HiSeq sequencing
   platform were processed using Trimmomatic (version 0.33) to remove
   low-quality base with quality scores < 20, short reads (< 40 bp) and
   “N” base. By aligning to the yak genome (BosGru3.1), host
   contaminations were removed using Bowtie software (version 2.2.4) to
   filter the reads that were of host origin, with the setting of
   –end-to-end, –sensitive, -I 200, -X 400^[164]55. Next, the clean data
   were assembled for each sample using SOAPdenovo software
   (V2.04)^[165]56, and the parameters were as follows: -d 1, -M 3, -R,
   -u, -F, -K 55. MetaGene ([166]http://metagene.cb.k.u-tokyo.ac.jp/) was
   used to predict open reading frames (ORFs) from the assembled contigs
   (length >= 500 bp)^[167]57. Assembled contigs were then pooled to
   remove non-redundancies using CD-HIT with 95% identity and 90% coverage
   for the selection of the representative sequences^[168]58.

Rumen metagenomic bioinformatics

   Regarding taxonomic prediction of rumen microbiome, DIAMOND
   software^[169]59 (V0.9.9) was used to blast the contigs to the
   sequences of bacteria, fungi, archaea and viruses which were all
   extracted from the nr database (Version: 2023-07-28) of NCBI with the
   parameter setting of blastp -e 1e-5. Rumen microbiome taxonomic
   profiles were analyzed at the domain, phylum, genus and species levels
   to identify the changes of microbiome composition. Beta diversity based
   on Bray-Curtis distances was visualized on PCoA plots. Microbial taxa
   with a relative abundance > 0.1% in at least 50% of yaks within each
   group were used for downstream analysis.

   For functional analysis, DIAMOND was used to annotate non-redundant
   genes against the KEGG (Kyoto Encyclopedia of Genes and Genomes)
   database with an E value of 1e-5^[170]60. Abundances of KEGG Orthology
   (KO), pathways, KEGG enzymes, and Modules were normalized for the
   subsequent statistics. Moreover, the Carbohydrate-Active EnZymes
   database (CAZy)^[171]61 were annotated by aligning genes to dbCAN
   database (HMMdb V8) with hmmscan program in HMMER (v3.1b2). The
   abundance values of the CAZymes were assigned by mapping the clean
   reads from each sample to the corresponding unigenes obtained after
   CD-HIT clustering. For all functional annotations, the annotated hit(s)
   with the highest score was used for the subsequent analyses.

Determination and analysis of rumen and plasma metabolome

   The rumen and plasma metabolome were analyzed using LC-MS/MS system
   combined quadrupole Orbitrap mass spectrometer (Q Exactive Orbitrap,
   Thermo Fisher Scientific, USA)^[172]62. In detail, 100 μL of the sample
   was transferred to an EP tube. After the addition of 400 μL of extract
   solution (acetonitrile/methanol = 1; containing isotopically labeled
   internal standard mixture), the tube was vortexed for 30 s, sonicated
   for 5 min in an ice-water bath, and incubated for 1 h at -40 °C to
   precipitate proteins. Then, the sample was centrifuged at 12000 rpm for
   15 min at 4 °C. The resulting supernatant was transferred to a fresh
   glass vial for downstream analysis. The quality control (QC) sample was
   prepared by mixing an equal aliquot of the supernatants from all the
   samples. The mobile phase consisted of 25 mmol/L ammonium acetate and
   25 ammonia hydroxides in water (pH = 9.75) (A) and acetonitrile (B).
   The analysis was carried with elution gradient as follows: 0 ~ 0.5 min,
   95% B; 0.5 ~ 7.0 min, 95% ~ 65% B; 7.0 ~ 8.0 min, 65% ~ 40% B;
   8.0 ~ 9.0 min, 40% B; 9.0 ~ 9.1 min, 40% ~ 95% B; 9.1 ~ 12.0 min, 95%
   B. The column temperature was 25 °C. The auto-sampler temperature was
   4 °C, and the injection volume was 2 μL. The QE HFX mass spectrometer
   was used for its ability to acquire MS/MS spectra on
   information-dependent acquisition (IDA) mode in the control of the
   acquisition software (Xcalibur, Thermo). The acquisition software
   continuously evaluates the full scan MS spectrum in this mode. The ESI
   source conditions were set as follows: sheath gas flow rate as 50 Arb,
   Aux gas flow rate as 10 Arb, capillary temperature 320 °C, full MS
   resolution as 60000, MS/MS resolution as 7500, collision energy as
   10/30/60 in NCE mode, spray Voltage as 3.5 kV (positive) or -3.2 kV
   (negative), respectively.

   The raw data of metabolomics was converted to the mzXML format using
   ProteoWizard and processed with an in-house program, which was
   developed using R and based on XCMS, for peak detection, extraction,
   alignment, and integration. Then, an in-house MS2 database (BiotreeDB)
   was applied to metabolite annotation. The cutoff for annotation was set
   at 0.3. In total, 263 rumen metabolites and 177 serum metabolites were
   identified and transformed to have a zero mean and a unit variance for
   downstream analysis. An online platform, MetaboAnalyst 5.0^[173]63, was
   used for the metabolomic analysis based on targeted metabolites
   utilizing the library of Bos Taurus (cow)^[174]64. Metabolite set
   enrichment analysis (MESA) was performed using MetaboAnalyst 5.0, based
   on the pathway-associated metabolite sets library^[175]65.

Statistics

   Growth performance, serum biochemical and physiological parameters, and
   rumen VFAs concentrations were compared using t test in R (v4.2.0)
   “stats” package. Rumen microbiome at the domain, phylum, and genus
   levels were compared using the Wilcoxon rank-sum test, with the
   Benjamini–Hochberg FDR adjusted P value < 0.05 being considered as
   significantly different. The dissimilarity of bacterial, archaeal,
   eukaryotic and viral composition was tested using Analysis Of
   Similarity (ANOSIM) performed on the R “vegan” package. The linear
   discriminant analysis effect size (LEfSe) algorithm^[176]66 based on an
   online Galaxy server was used to detect the features differentiating
   grazing and intensive groups, such as rumen microbial species, KEGG
   pathways and modules, and CAZyme families. The cutoff for significant
   differences of LEfSe was LDA score > 2 and P value < 0.05.

   MetaboAnalyst 5.0 was employed to perform the statistical analysis for
   both rumen and plasma metabolomic data. The orthogonal projections to
   latent structures- discriminant analysis (OPLS-DA). The t test was
   performed between the two groups, with the FDR adjusted P value of
   Benjamini–Hochberg and the Variable Importance in the Projection
   (VIP) > 1 being considered as significantly different metabolites.
   Additionally, ‘pheatmap’ package in R was used to display the signature
   metabolites in rumen and plasma differentiating two groups. The
   analysis of traceability and enrichment of rumen metabolites was
   performed using MetOrigin^[177]67.

   Spearman’s analysis was performed to calculate the correlations between
   signature microbiota and critical metabolites, using ‘psych’ package in
   R, and the network was visualized using the Cytoscape (version 3.10.1).
   Only significant connections were shown in the network (|r | > 0.7,
   P < 0.05).

Reporting summary

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

Supplementary information

   [179]Supplemental figures and tables^ (2.8MB, pdf)
   [180]Reporting summary^ (2.8MB, pdf)

Acknowledgements