Dietary forage to concentrate ratios impact on yak ruminal microbiota and metabolites.

To improve the rumen fermentation function and growth performance of yaks (Bos grunniens), better understanding of the effect of different dietary forage to concentrate ratios on rumen microbiota and metabolites is needed. In the present study, three diets with different dietary forage to concentrate ratios (50:50, 65:35, and 80:20) were fed to 36 housed male yaks. The changes in the distribution of rumen microorganisms and metabolites and the interactions between them were studied by 16S rRNA gene sequencing and liquid chromatography-mass spectrometry (LC-MS). The diversity and richness of microorganisms in the rumen varied according to diet. The most abundant phyla were Firmicutes and Bacteroidetes. Firmicutes was the most abundant in the C50 group, and the relative abundance of Bacteroidetes was significantly lower in the C65 group than in the C80 group (p < 0.05). The Christensenellaceae_R-7_group, Rikenellaceae_RC9_gut_group, and Methanobrevibacter had the highest relative abundances at the genus level. Among them, Christensenellace_R-7_group had the highest relative abundance in the C50 group. The Rikenellaceae_RC9_gut_group was significantly abundant in the C80 group compared with the C50 group. The Methanobrevibacter content was higher in the C65 group than in the other two groups. Both the concentration and metabolic pathways of rumen metabolites were influenced by the dietary concentrate ratio; lipids, lipid-like molecules, organic acid metabolites, and organic oxide-related metabolites differed between the groups. Significant changes were found for six metabolic pathways, including arginine and proline metabolism; glycine, serine, and threonine metabolism; glyoxylate and dicarboxylate metabolism; arginine biosynthesis; glycerophospholipid metabolism; glycerolipid metabolism; and nitrogen metabolism.

Introduction

Yaks (Bos grunniens), an endemic livestock species living on the Qinghai–Tibetan Plateau (QTP), are an important means of production and livelihood for local herders (Long et al., 2008). In recent years, with the improvement of human living standards, human demand for livestock products is increasing. To meet the increasing demand for livestock products in the market, the production performance of ruminants is generally improved by increasing the proportion of concentrate feed in the diet. However, in actual production, achieving the optimal proportion of concentrate feed in the diet is rare; concentrations that are too high or too low affect the economic efficiency of breeding. Therefore, it is important to choose an appropriate ratio of dietary forage to concentrate to improve the breeding of yak.

The rumen is an important site for the digestion of feed in ruminants and contains a large number of bacteria, fungi, protozoa, and archaea (Qiu et al., 2020). Ration composition and nutrient content are the most important factors influencing rumen microbiota (Yáñez-Ruiz et al., 2015). Therefore, understanding the response of the rumen microbiota of ruminants to dietary changes is essential to improve the growth performance of ruminants. The effect of diet structure on the rumen microbiota of yaks has become a new hot topic in recent years. At present, many studies have investigated the rumen microbes of yaks (Fan et al., 2021; Wu et al., 2021). For example, Ahmad et al. (2020) investigated the effect of different dietary energy levels on the rumen microbiota of yaks. However, there is little published information on the effect of dietary concentrate to forage ratio on the composition of the yak rumen microbial community, therefore, our study focused on this area.

Since the health and growth of ruminants depend heavily on the production of metabolites in the rumen, a comprehensive analysis of rumen fluid provides insight into the interaction between the diet fed to ruminants and the rumen (Saleem et al., 2013). In recent years advanced techniques have been able to identify the metabolites of rumen microbiota, and it has been found that changes in some metabolites may affect microbiota alterations and host activity. Previous studies have shown that approximately 55%–60% of rumen fluid metabolites are associated with rumen microbiota (Saleem et al., 2013). The interaction between microbiota and metabolites facilitates the interpretation of microbiota classification and functional properties of some metabolites (Heinken et al., 2014). Studies on rumen microbiota and metabolites in yaks have revealed that the microbiota composition is significantly different between the concentrate and roughage groups and that the type of feed significantly alters the concentration and metabolism of metabolites involved in protein digestion and absorption (Liu et al., 2019b). Studies on the rumen microbiota and metabolic profile of Holstein heifers found that rumen microbiota were significantly affected and rumen metabolites amino acids, lipids, organic acids, and carbohydrates were significantly changed by changing the dietary forage to concentrate ratios (Zhang et al., 2017).

Although microbiomic and metabolomic studies have been widely used in animal nutrition, there is a lack of reports on the interaction between rumen microbiota and metabolites in yaks at different dietary forage to concentrate ratios, mainly owing to the differences in nutritional and metabolic characteristics between yaks and other breeds of cattle, the specificity of the growing environment, and the lack of technical guidance. Therefore, it is possible that the results of other ruminant studies may not be adapted to the diets of yaks. In this study, the rumen microbiota and metabolites of yaks were analyzed by combining 16S rRNA gene sequencing and LC–MS techniques to illustrate the effects of three different dietary forage to concentrate ratios.

Materials and methods

Animals, diets, and experimental design

This study was conducted from September to December 2019 at the Laozaxi breeding site in Guinan County, Qinghai Province, China. In total, 36 3-year-old healthy male yaks (weight: 164.9 ± 12.9 kg) were selected from grazing pastures and randomly divided into three groups. The yaks in the three groups were fed diets with a 50:50 (C50), 65:35 (C65), or 80:20 (C80) concentrate to forage ratio. According to the survey, the current level of concentrate used in farms is about 50–80%, however, the level of concentrate of high concentrate is usually 70 -80%, and the level of concentrate below 50% does not have the effect of promoting production, so we screened between 50 and 80% and set a three gradient concentrate to forage ration. The diets were prepared according to the Chinese Beef Cattle Feeding Standard (NY/T815-2004), and the composition and nutrient content of the diets are shown in Table 1. All yaks were uniformly numbered and fed alone with free access to water and fed twice a day at 08:00 and 17:00. The pretest period was 15 days, and the normal test period was 90 days.

Table 1

Ingredients and nutritional composition of each diet.

Ingredients (%)	Group
	C50	C65	C80
Oat hay	50.00	35.00	20.00
Corn	22.68	29.75	37.15
Wheat	6.25	8.39	10.42
Wheat bran	6.43	8.56	10.67
Rapeseed meal	6.36	8.55	10.61
Soybean meal	2.17	2.91	3.62
Palm oil powder2	2.11	2.84	3.53
CaHPO4	1.00	1.00	1.00
NaCI	1.00	1.00	1.00
Premix1	2.00	2.00	2.00
Nutrient composition (%)
CP	12.41	13.18	13.87
ME MJ/kg	10.98	11.89	12.80
NDF	37.67	31.36	25.48
ADF	23.32	18.36	13.61
Ca	0.45	0.48	0.50
P	0.48	0.56	0.63

The premix provided the following per kg of the diet: VA 3500 IU, VD 1000 IU, VE 40 IU, Mn 40 mg, Fe 50 mg, Cu 10 mg, Zn 40 mg, Se 0.3 mg;

Others were measured values.

Crude protein (CP), neutral detergent fiber (NDF), acid detergent fiber (ADF), calcium (Ca), and phosphorus (P) content were determined according to the official reference method were determined in each sample in the laboratory, and ME was calculated. Mixed feeds (100 g) were collected and dried in a forced-air oven at 60°C for 48 h and then ground through a 1-mm sieve before analysis. CP、Ca and N contents were determined according to AOAC (2007). The NDF and ADF contents were determined by the method of Vansoest et al. (1991).

Rumen sample collection and measurements

On day 90 of the formal trial, 150 ml of rumen fluid samples were collected using a gastric tube sampler before morning feeding. Ruminal fluid was collected from 36 yaks, and five ruminal fluid samples were randomly selected from each group for measurement and gene sequencing. Rumen fluid samples were filtered through 4 layers of gauze and rumen pH was measured immediately using a benchtop acidity meter (Model HI221, HANNA, Italy). The remaining rumen fluid samples were dispensed into 15 ml centrifuge tubes and immediately frozen in liquid nitrogen, then transferred to a − 80°C refrigerator. The filtered rumen fluid was centrifuged (17,000 × g for 30 min at 4°C) to obtain the supernatant, which was further analyzed for NH3-N using phenol hypochlorite analysis (Broderick and Kang, 1980). Freshly prepared metaphosphoric acid (25% w/v, 2 ml) was added to 8 ml of filtered rumen fluid and then centrifuged (17,000 × g for 10 min at 4°C). Volatile Fatty Acids (VFAs) concentrations were determined using gas chromatography (GC-2014; Shimadzu Corporation, Japan) as described by Cao et al. (2008). The determination conditions were as follows: flame ionization detector (FID), capillary column (FFAP, 30.00 m × 0.32 mm × 0.50 μm); The heating conditions were as follows: the initial temperature was 60°C, and the temperature was increased to 120°C at 10°C/min and retained for 2 min; the temperature was increased to 180°C at 15°C/min and retained for 5 min; the vaporization chamber temperature was 250°C. FID temperature 250°C; Sample size: 1 μl, carrier gas: high purity nitrogen (99.99%), pressure: 0.7 MPa; Hydrogen pressure 0.4 MPa, air pressure 0.4 MPa, capillary column pressure 0.6–0.8 MPa, shunt ratio 40:1.

Microbial DNA extraction, PCR amplification, sequencing, and sequencing data processing

The CTAB method was used to extract microbial DNA from rumen fluid samples. We used 1.0% agarose gel electrophoresis to detect the concentration and purity of DNA. According to the concentration, the DNA was diluted to 1 ng/μl with sterile water. The extracted DNA was amplified by PCR, using barcoded specific primers 515F (5′-GTGCCAGCMGCCGCGG-3′) and 806R (5′-GTGCCAGCMGCCGCGG-3′) to amplify the V3 – V4 hypervariable region of 16S rRNA gene. We used 25 μl amplification system, 5 μmol/l upstream and downstream primers, and ~5 ng template DNA for the PCR reaction. PCR amplification conditions were as follows: 94°C pre-denaturation treatment for 5 min, a denaturation cycle at 94°C for 30 s, annealing at 50°C for 30 s, and extension at 72°C for 60 s for a total of 30 cycles; finally, a 72°C extend for 7 min. The PCR amplification products were detected by 1.0% agarose gel electrophoresis; the recovered products were purified using the MinElute Gel Extraction Kit (Qiagen, Germany) and TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina, United States). The constructed library was quantified by Qubit and Q-PCR, and then sequenced on the Illumina NovaSeq PE250 platform. PCR amplification, PCR product mixing and purification, library construction, and computer sequencing processes were all produced by Novogene Bioinformatics Technology Co., Ltd (Beijing, China) to produce 250 bp paired-end readings.

The sequences obtained from the NovaSeq PE250 platform were processed through open-source software pipeline QIIME (Quantitative Insights into Microbial Ecology) version 1.8.0-dev (Caporaso et al., 2010), with the criteria as described by previous reports (Mao et al., 2015; Ye et al., 2016). Briefly (1) threads that had a mean quality score of no <20 and no shorter than 50 bp were retained; (2) we discarded reads that had exact barcode matching, two nucleotide mismatches in primer matching, and/or ambiguous characters; (3) only sequences that overlapped by more than 10 bp were assembled according to their overlap sequence. Reads that could not be assembled were discarded. Sequences were binned into operational taxonomic units (OTUs) based on 97% identity using UCLUST (version7.1, and chimeric sequences were identified and removed by UCHIME (Edgar, 2010). The most abundant sequence within each OTU from specific libraries (bacteria) was designated as the “representative sequence” and aligned against the (SILVAbacterial) database (version 119; Pruesse et al., 2007), with the default parameters set by QIIME. Community richness and diversity, analyzed with measures such as Good’s coverages, observed species, PD whole tree, Chao1 and Shannon indices, weighted uniFrac distance-based principal coordinate analysis (PCoA), and weighted distance-based analysis of molecular variance (AMOVA), which were used to illustrate significant differences among the samples, were assessed by the program MOTHUR v.1.35.0 (Schloss et al., 2009).

All the raw sequences (16S) after assembling and filtering were submitted to the NCBI Sequence Read Archive, under accession number PRJNA844537.

GC-TOF/MS and identification of compounds

Rumen fluid samples were dried completely in a vacuum concentrator without heating, and then 60 μl of methoxyamination hydrochloride (20 mg/ml in pyridine) was added and incubated for 30 min at 80°C. BSTFA reagent (80 μl, 1% TMCS, v/v) was added to the sample aliquots and incubated for 1.5 h at 70°C. Samples were analyzed by an Agilent7890 gas chromatograph system coupled with a Pegasus HT time-of-flight mass spectrometer (GC-TOF/MS). The system had a DB-5MS capillary column (30 m × 250 μm inner diameters, 0.25 μm film thickness, J&W Scientific, Folsom, CA, USA) coated with 5% diphenyl cross-linked with 95% dimethylpolysiloxane. The temperature was initially kept at 50°C for 1 min, and then increased to 320°C at a rate of 10°C /min; the column was then maintained for 5 min. The temperatures of injection, transfer line, and ion source were 280°C, 270°C, and 220°C, respectively. The mass spectrometry data were acquired in full-scan mode with a mass-to-change ratio (m/z) range of 30–600 at a rate of 20 spectra/s after a solvent delay of 366 s.

Chroma TOF 4.3X software of the LECO Corporation and LECO-Fiehn Rtx5 database (Kind et al., 2009) were used for raw peak exacting, data baselines filtering, calibration of the baseline, peak alignment, deconvolution analysis, peak identification, and integration of the peak area. Both mass spectrum match and retention index match were considered in the metabolite identification.

The resulting peak number, sample name, and normalized peak area data were fed into the SIMCA software package (V14, Umetrics AB, Umea, Sweden) for principal component analysis (PCA), and orthogonal projections to latent structures-discriminant analysis (OPLS-DA). Differentially expressed metabolites (DEMs) were identified, combing Variable Importance in Projection (VIP) obtained from OPLS-DA analysis and t-test (VIP > 1 and p < 0.01). DEMs were further identified and validated by the Bovine Metabolome Database (BMDB) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The variation of different metabolites with diet concentrate ratio was identified and the horizontal, linear, and secondary effects were evaluated for the treatment groups (Zhang et al., 2020). DEM was imported into the MetaboAnalyst web server to view its metabolic pathway distribution and enrichment analysis (Xia et al., 2009).

Correlations between microbial communities and rumen metabolites

Ruminal metabolites with VIP > 1 and p < 0.05, along with significantly affected microbiota were selected for use in R (v3.2.4) to study changes in relevant metabolic processes. For Spearman correlation analysis, p values were calculated using the Psych package (author, W. Revele; publication date, 2016; version, 1.6.9), with absolute Spearman correlations of p > 0.05. These correlations were visualized by Pheatmap (author, R. Kolde; publication date, 2015; version 1.0.8) R and packages in Cytoscape 2.8.2 (Smoot et al., 2011).

Statistical analysis

Statistical analyses were performed using R (v 3.6.1) and SPSS 26.0. According to Zhang et al. (2020), the linearity of the treatments was assessed by the lm function of the estimability package in R.

Results

Rumen fermentation parameters

Ruminal pH decreased linearly (p < 0.01) with increasing dietary concentrate level and had the highest value in the C50 group and the lowest value in the C80 group. Acetate concentration was significantly higher in the C80 group than in the C50 and C65 groups, and was higher in the C50 group than in the C65 group. Ruminal isobutyrate increased quadratically (p = 0.002), and ruminal TVFA, propionate, and butyrate concentrations increased linearly (p < 0.05), all of which were highest in the C80 group and lowest in the C50 group (Table 2).

Table 2

Effect of dietary forage to concentrate ratios on rumen fermentation parameters.

Items	Treatmentsa			SEMb	p-values
	C50	C65	C80		Treatment	Linear	Quadratic
pH	6.49	6.37	5.94	0.077	0.001	0.001	0.411
NH3-N,mg/dL	12.74	12.58	13.24	0.154	0.19	0.192	0.283
MCP	2.50	2.79	2.34	0.114	0.276	0.592	0.161
TVFA, mM	76.31	80.25	90.45	2.044	0.004	0.001	0.654
VFAs, molar of TVFA
Acetate	43.62	42.78	49.55	1.068	0.004	0.018	0.054
Propionate	15.46	17.88	19.52	0.761	0.076	0.022	0.495
Butyrate	12.49	14.13	15.31	0.479	0.038	0.01	0.455
Valerate	1.32	1.52	1.67	0.078	0.174	0.057	0.174
Isobutyrate	0.65	0.83	1.06	0.057	0.002	<0.001	0.002
Isovalerate	2.78	3.15	3.34	0.123	0.141	0.055	0.383
Acetate: Propionate	2.95	2.39	2.55	0.114	0.157	0.157	0.086

C50, diet contained 50% of concentrate; C65, diet contained 65% of concentrate; C80, diet contained 80% of concentrate.

SEM, standard error of the mean.

Bacterial diversity among treatment groups

After screening, a total of 1,157,681 high-quality 16S rRNA gene sequences were obtained from 15 rumen fluid samples. After subdividing each sample into equal sequencing depths (77,178 reads per sample) and clustering, 4,289 OTUs were obtained with 97% concordance. Microbial diversity and richness differed among the three groups according to the Shannon index and Chao1 values, indicating higher diversity and richness in the C50 group and lower diversity and richness in the C65 group, but the differences were not significant (Figures 1A,B). PCoA plots using unweighted Nifrac matrix distances, which show the diversity of the rumen microbiota for bacterial communities clustered by different concentrate to forage ratio diets, clearly showed the structure of the different bacterial communities in the three groups (p < 0.05). Among them, the C65 group was significantly different from the other two groups (Figure 1C), indicating that different dietary concentrate to forage ratios affected the compositions of the bacterial communities.

Figure 1

Bacterial diversity analysis of yak rumen samples for different dietary concentrate ratios. Alpha diversity between different groups’ (A) Chao 1 value and (B) Shannon index; (C) beta diversity principal coordinate analysis (PCoA). Differences based on concentrate to forage ratio of different diets were used.

Bacterial compositions among treatment groups

Across all groups, 59 phyla, 145 classes, 291 orders, 377 families, and 601 genera were identified. The relative abundance of dominant taxa at the phylum and genus levels is shown in Figures 2A,B. At the phylum level, the most abundant phyla were Firmicutes and Bacteroidota, which accounted for 56.56 and 19.57% of the total readings, respectively (Figure 2A). Firmicutes were enriched in the C50 group. The relative abundance of Bacteroidota in the C65 group was significantly lower than that in the C80 group (p < 0.05), and was followed by Actinobacteriota (4.22%), unidentified-Bacteria (6.65%), Euryarchaeota (6.43%), Proteobacteria (2.13%), Chloroflexi (0.83%), Verrucomicrobiota (0.33%), Spirochaetota (0.42%), Acidobacteriota (0.47%. At the genus level, Christensenellaceae_R-7_ group (19.27%), Rikenellaceae_RC9_gut_ group (8.15%), and Methanobrevibacter (6.36%) were the dominant genera. The content of the Christensenellace_R-7_ group was the highest. Compared with the C50 group, the Rikenellaceae_RC9_gut_ group of the C80 group was significantly richer (p < 0.05). The content of Methanobrevibacter in the C65 group was higher than that in the other two groups. Other genera with lower abundance included NK4A214_group (5.14%), Candidatus_Saccharimonas (5.00%), [Eubacterium]_ventriosum_group (3.02%), Olsenella (2.95%), Lachnospiraceae_NK3A20_group (2.57%), Prevotella (2.31%), and Ruminococcus (1.86%). The number of NK4A214_group in the C80 group was significantly higher than the other two groups (p < 0.05; Figures 2B,C).

Figure 2

Bacterial composition of rumen samples from yaks with different dietary concentrate to forage ratio. Bacterial composition at the phylum (A) and genus (B) levels; (C) significantly different bacterial phylum and genus between groups, where different letters indicate significant differences between groups.

Identification, quantification, and statistical comparison of GC-TOF/MS metabolites in rumen

Supplementary Figure 1 shows the Pearson correlation coefficients between QC samples. Metabolites from rumen fluid of yaks fed with different concentrate to forage ratio diets, including QC samples, were analyzed by LC–MS, and PCA after positive and negative mode ionization showed the main unsupervised separation between the three groups (including QC samples; Supplementary Figure 2) To better distinguish the differences between the groups, OPLS-DA positive and negative mode ionization was performed to create scoring plots (Figure 3) to verify the differential metabolites between the two groups. We found that the two groups were clearly separated, with corresponding R2Y values of 0.946, 0.949, 0.96, 0.954, 0.958, and 0.978 for the OPLS DA model in C50 vs. C65, C65 vs. C80, and C50 vs. C80, respectively. All samples in the scoring plots were within 95% of the Hotelling T2 ellipse, indicating that the OPLS-DA model had satisfactory validity. The scoring plots showed clear separation and distinction between the different concentrate to forage ratio diets, indicating that the OPLS-DA model could identify differences among the groups.

Figure 3

Corresponding validation plots and orthogonal projections to latent structures-discriminant analysis (OPLS-DA) score plots derived from gas chromatography time-of-flight mass spectrometer (GC-TOF/MS) metabolite profiles of yak rumen samples with different dietary concentrate ratios. Corresponding validation plots and OPLS-DA score plots for (A,B) group C50 vs. group C65 (C,D) group C50 vs. group C80, and (E,F) group C65 vs. group C80. C50, diets with 50% concentrate; C65, diets with 65% concentrate; C80, diets with 80% concentrate.

Metabolomic profiles in rumen

A total of 999 metabolites were identified. As shown in Supplementary Table 1, 107 differential metabolites were detected in the C50 vs. C65 group, 163 differential metabolites were detected in the C50 vs. C80 group, and 130 differential metabolites were detected in the C65 vs. C80 group (VIP > 1 and p < 0.05). Overall, 55 differential metabolites were identified among the three treatment groups, of which 20 were classified as lipids and lipid-like molecules, 14 were classified as organic acids and derivatives, 4 were classified as organoheterocyclic compounds, 5 were classified as organic oxygen compounds, 4 were classified as benzenes, 3 were classified as nucleosides, nucleotides, and analogues, 2 were classified as organic nitrogen compounds, 2 were classified as phenylpropanoids and polyketides, and 1 was classified as alkaloids and derivatives (Table 3). As shown in Table 3, changes to the dietary concentrate to forage ratio included linear and quadratic increases and decreases.

Table 3

Comparison of significantly changed rumen metabolites in three treatments.

Superclass	Metabolite names	Treatmentsa			SEMb	q-values
		C50	C65	C80		Treatment	Linear	Quadratic
Lipids and lipid-like molecules	LPC 18:1	194691104.7	42708696.68	73549354.12	27916602.05	0.049	0.057	0.056
	(+/−)14(15)-DiHETE	7896496.452	6041559.384	10162068.62	673300.6846	0.028	0.113	0.038
	PE (18:0/18:1)	568223637.8	162,546,434	310528017.1	55846937.5	0.002	0.013	0.002
	Oleanolic acid	10860029.05	25271500.23	31586962.46	3021353.7	0.005	0.002	0.165
	Docosanoic acid	120465591.8	312561086.8	278834439.5	29895781.66	0.007	0.011	0.014
	FAHFA (18:0/20:2)	298961990.2	196,149,233	131557819.5	22268756.13	0.001	<0.001	0.192
	LPS 19:0	38712801.23	26590129.57	17910065.26	2915932.416	0.003	0.001	0.325
	1a,1b-Dihomo prostaglandin E1	6101762.044	10267200.84	15370526.54	1314266.816	0.004	0.001	0.728
	9-KODE	5199482.404	10792441.07	14612730.08	1337414.414	0.004	0.001	0.315
	FAHFA (16:0/19:2)	14935684.43	13776392.47	6942936.672	1166944.149	0.002	0.001	0.271
	8-Isoprostaglandin E2	7185522.319	7089838.662	2292831.814	802492.4605	0.006	0.004	0.175
	5-Phenylvaleric Acid	4951271.466	2543348.639	2308797.869	403568.0036	0.003	0.002	0.029
	Ursolic acid	67370520.13	262907645.3	329005364.6	43450056.48	0.023	0.009	0.209
	PC (17:2/18:5)	10851812.01	3702900.431	3488810.378	1207361.893	0.006	0.004	0.033
	(±)12(13)-DiHOME	1,108,505,784	557158819.6	1,438,674,449	148732595.2	0.035	0.289	0.024
	FAHFA (20:0/22:3)	3800547.937	4999730.565	3051165.382	356261.5261	0.067	0.339	0.044
	FAHFA (20:0/20:0)	2552261.78	3927179.497	2148783.717	298715.7078	0.024	0.499	0.011
	FAHFA (20:1/22:3)	1225762.082	3295555.767	1397935.437	296420.9791	0.001	0.688	0.000
	octadec-9-ynoic acid	142275777.6	119237853.7	358769046.8	35729092.27	0.002	0.002	0.054
	Methyl dihydrojasmonate	39830838.48	25142191.38	55470567.17	4440973.719	0.008	0.069	0.011
Organic acids and derivatives	4-Acetamidobutanoic acid	339558459.3	705539699.7	614883543.4	83776799.16	0.183	0.178	0.145
	Glycohyocholic acid Sodium salt	1121650.667	2160565.954	1293457.42	171128.2522	0.015	0.6	0.004
	Glycochenodeoxycholic acid sodium salt	138001.0444	102694.7378	210195.6813	27122.37049	0.273	0.284	0.290
	N6,N6,N6-Trimethyl-L-lysine	22992800.87	49865733.38	43396522.91	3542298.6	<0.001	0.001	0.025
	3-(1-cyano-1,2-dihydroisoquinolin-2-yl)-3-oxopropyl propionate	9809639.689	5114501.201	3930670.708	847991.9981	0.002	0.001	0.056
	2-Phenylglycine	8903501.364	12633774.31	26737460.27	2435616.441	0.001	<0.001	0.349
	FLK	4400166.146	2468044.319	1353037.499	397453.436	0.001	<0.001	0.133
	L-Cystathionine	1626295.198	1518461.996	935345.7416	109842.2172	0.009	0.004	0.395
	Val-Ser	74101368.05	56742398.55	35016005.23	5274132.354	0.002	<0.001	0.712
	Gly-Tyr-Ala	10585113.74	5431686.462	3484016.941	1071800.745	0.007	0.003	0.157
	Methylmalonic acid	6,854,758,528	10,232,838,938	6,095,994,035	644249580.7	0.008	0.517	0.003
	D-Proline	8790839.633	21238605.66	12225921.38	1605430.561	<0.001	0.122	0.000
	Methyl-2-hydroxyisobutyric acid	5627003.487	11841335.87	6323869.199	993122.1649	0.007	0.696	0.002
	5-S-cysteinyldopa	359632.8043	1151730.365	622027.8277	102360.9832	<0.001	0.082	0.000
Organoheterocyclic compounds	Uracil	2,921,751,171	2,509,346,257	2,244,999,054	152680217.3	0.196	0.079	0.617
	Roquefortine C	121793862.4	69710112.4	44521541.9	12,205,993	0.016	0.006	0.272
	Vitamin B2	9657353.128	8023761.343	5587224.07	626803.3111	0.014	0.004	0.918
	6-Dimethylaminopurine	44812097.53	20080013.45	16138724.04	4571946.438	0.008	0.004	0.072
Organic oxygen compounds	D-(+)-Mannose	1,220,583,273	2,446,884,559	1,865,896,826	184868124.8	0.012	0.08	0.006
	L-Threonic acid	75199858.64	153399194.6	83824327.61	14171561.31	0.032	0.765	0.010
	D-Galactosamine	3924881.77	9228256.303	6480804.621	697657.3705	0.001	0.029	0.000
	L-Kynurenine	7186761.226	8638010.59	9452862.755	354037.1177	0.016	0.005	0.320
	L-(−)-Glyceric acid	241945895.5	324209480.6	188479467.4	21370398.42	0.018	0.211	0.014
Benzenoids	4-Hydroxybenzoic acid	18829863.68	11260717.01	11002717.31	1431899.843	0.025	0.016	0.082
	Isohomovanillic acid	49367075.51	26081257.83	39488821.6	4000030.654	0.044	0.249	0.018
	Homovanillic acid	21776146.29	18347209.78	15121282.09	922797.1842	0.003	0.001	0.432
	2-[(carboxymethyl)(methyl)amino]-5-methoxybenzoic acid	106638622.4	324209480.6	188479467.4	21370398.42	0.085	0.321	0.060
Nucleosides, nucleotides, and analogues	2’-O-Methylguanosine	102175235.9	172895328.6	115928891.7	12413194.3	0.033	0.588	0.010
	1-Methylguanosine	13592217.39	17438884.08	5796441.008	1554608.446	0.001	0.005	0.004
	Pseudouridine	70207322.22	143837465.1	74685038.63	11250678.16	0.002	0.806	0.001
Organic nitrogen compounds	Spermidine	229394701.5	85090134.18	46356478.24	27157907.31	0.004	0.002	0.081
	2,6-di(2-thienylmethylidene)cyclohexan-1-one	92684266.52	47489798.42	47789129.5	6922044.883	0.001	0.001	0.009
Phenylpropanoids and polyketides	LSD-d3	29688333.54	19809762.59	17268311.85	1820801.817	0.003	0.001	0.068
	3-hydroxy-3,4-bis[(4-hydroxy-3-methoxyphenyl)methyl]oxolan-2-one	1912453.091	1004142.123	1065828.76	127195.7192	<0.001	<0.001	0.002
Alkaloids and derivatives	6-Acetylmorphine	12992029.47	7645667.857	5767021.69	1116835.001	0.01	0.004	0.160

C50, diet contained 50% of concentrate; C65, diet contained 65% of concentrate; C80, diet contained 80% of concentrate.

SEM, standard error of the mean. Significance considered at p < 0.05.

Metabolic pathways of differential metabolites

Pathway topology analysis was performed based on the identified metabolites and concentrations. Figure 4 shows seven enriched major metabolic pathways among the three treatment groups, including: arginine and proline metabolism, glycine, serine and threonine metabolism; glyoxylate and dicarboxylate metabolism; arginine biosynthesis; glycerophospholipid metabolism; glycerolipid metabolism; and nitrogen metabolism. The above pathways were identified as significantly different (p < 0.05). Seven metabolites (spermidine, D-proline, 4-acetamidobutanoic acid, phosphatidylcholine, 5-phenylvaleric acid, D-galactosamine, and L-threonic acid) were localized to these important pathways according to the KEGG pathway identification, and thus were identified as characteristic metabolites of 55 metabolic species.

Figure 4

Pathway analysis of differential metabolites in rumen samples from yaks with different dietary concentrate ratios using the Bos taurus Kyoto Encyclopedia of Genes and Genomes (KEGG) database of MetaboAnalyst 4.0. Larger circles indicate richer pathways and darker colors indicate higher pathway impact values. The closer the color is to red, the smaller the p-value.

Correlations between the ruminal metabolomes and microbiomes

Figure 5 represents the correlation analysis of rumen microbiota with metabolites. In total, 24 metabolites (from 55 differential metabolites) and 10 microbial genera were screened, and the results show that the dynamic fluctuations of certain metabolites were closely correlated with the abundance of various microbial groups. Ruminococcus was positively correlated with 3-(1-cyano-1,2-dihydroisoquinolin-2-yl)-3-oxopropyl propionate and LPC 18:1; Christensenellaceae_R-7_group was positively correlated with Homovanillic acid; NK4A214_group was positively correlated with L-cystathionine, Gly-Tyr-Ala, 3-(1-cyano-1,2-dihydroisoquinolin-2-yl)-3-oxopropyl propionate, FAHFA (18:0/20:2), 1-methylguanosine, and uracil; Rikenellaceae_RC9_gut_group was positively correlated with 14, 15-DiHET; 4-acetamido butanoic acid was positively correlated with L-cystathionine, Gly-Tyr-Ala, FAHFA (18,0/20:2), and spermidine; Prevotella was positively correlated with L-kynurenine, but negatively correlated with 3-(1-cyano-1,2-dihydroisoquinolin-2-yl)-3-oxopropyl propionate, Vitamin B2, and 5-phenylvaleric acid.

Figure 5

Correlation matrix between genus level and rumen differential metabolites. Each row in the figure represents a metabolite, each column represents a genus, and each lattice represents the Pearson correlation coefficient between metabolite and genus levels. Red color indicates positive correlation, while blue color indicates negative correlation. *p < 0.05; **p < 0.01.

Discussion

Dietary and nutrient compositions are important factors affecting rumen microbiota in ruminants (Yáñez-Ruiz et al., 2015). An appropriate diet concentrate to forage ratio not only improves rumen microbiota (Shen et al., 2017) but also promotes the growth and development of ruminants. In this study, we investigated the effect of dietary concentrate to forage ratio on rumen microbiota and metabolome in yaks. The results provide a reference for yak growth by regulating the relationship between the dietary concentrate to forage ratio and rumen microbiota.

A stable rumen environment is particularly important for ruminants; pH, NH3-N, and VFAs molar concentrations are important indicators of a stable rumen environment and reflect the status of rumen fermentation (Tomczak et al., 2019). In general, the most suitable pH for microbial growth and reproduction in the rumen is between 6.2 and 7.0, and the ecological environment of rumen microorganisms will be relatively stable within this range. When the rumen pH is below 6.0, the growth and reproduction of fibro-degrading bacteria in the rumen of ruminants is severely inhibited, affecting the digestion of fiber and inducing rumen acidosis (Grilli et al., 2016). In this study, the rumen pH of yaks in the C50 and C65 groups ranged from 6.37 to 6.49, while that of the C80 group was 5.94, indicating that rumen acidosis was induced when the percentage of concentrate increased to a certain level. Ruminal pH decreased linearly with increasing dietary concentrate level, which is consistent with the findings of Yang et al. (2000); this result may be due to the high level of dietary concentrate in the C80 group, which produced a large amount of VFAs under the degradation of rumen microorganisms, resulting in a decrease of ruminal pH.

VFAs is the end product of carbohydrates and proteins in the rumen rations of ruminants. In this study, the rumen TVFA, propionate, and butyrate concentrations increased linearly with increasing dietary concentrate level, similar to the findings of past studies (Klevenhusen et al., 2013; Giger-Reverdin et al., 2014; Polyorach et al., 2014). However, the acetic acid concentration, which was lower in the C65 group than in the C50 group and increased in the C80 group, was not consistent with the results of these past studies. This result may be due to the high level of concentrate in the C80 group, which affected the rumen environment and rumen microbial growth, resulting in abnormal rumen fermentation. Ration concentration to forage ratio increased from 50:50 to 65:35 and had a specific promotion effect on rumen TVFA content, but the ratio of concentration to forage 80:20 ratio had a negative effect instead. This study obtained the best result with a 65:35 forage to the acceptable ratio in the yaks’ diet.

The 16S rRNA gene high-throughput sequencing technology can provide a rapid and comprehensive understanding of the effect of diet concentrate to forage ratio on rumen microbiota of yaks. In this study, sequencing results revealed most of the microbiota in the samples. The Chao1 index measured the abundance of the microbiota, with a higher Chao1 index indicating greater microbiota abundance in the samples. The Shannon index measured the diversity of the microbiota, with a higher Shannon index indicating greater microbiota diversity in the samples. Rumen fluid samples from yaks fed in the C65 group had lower rumen microbiota richness and diversity compared with those fed in the C50 group, which is consistent with the results of other studies (Liu et al., 2019c; Wang et al., 2020; Ramos et al., 2021). However, rumen fluid samples from yaks reared in the C80 group showed higher abundances and diversity of rumen flora compared with yaks reared in the C65 group. It is possible that the rumen microbiota was disturbed due to the feeding of high concentrate (Krause and Oetzel, 2006; Moya et al., 2009), the rumen microbiota underwent value-added, together with the composition of the diet and the special breed of yaks, this resulted in a higher abundance and diversity of rumen microbiota in the C80 group of yaks. PCoA analysis showed that rumen microbiota clustered according to different dietary concentrate ratios and significant differences between rumen microbiota were observed, which is consistent with other studies (Hu et al., 2019; Islam et al., 2021).

In the current study, Firmicutes and Bacteroidota were the dominant phylum in all three groups of yak rumen microorganisms, similar to previous studies (Zhou et al., 2017; Bi et al., 2018). Firmicutes is mainly involved in the catabolism of fibrous material; whereas, Bacteroidota is mainly involved in the degradation of non-fibrous material (Evans et al., 2011; Reigstad and Kashyap, 2013), which is important for the nutritional metabolism of ruminants. The rumen Bacteroidota and Firmicutes were higher in the rumen of yaks in the C50 group than in the C65 group, indicating a decreasing trend of rumen Bacteroidota and Firmicutes with increasing dietary concentrate levels, consistent with the results of past studies (Nagaraja and Titgemeyer, 2007; Chiquette et al., 2008). However, the abundances of both Bacteroidota and Firmicutes in the rumen of yaks in the C80 group were higher than those in the C65 group, reflecting increasing relative abundance after rumen acidosis in yaks in the C80 group. This is consistent with the results of past studies (Mao et al., 2013; Plaizier et al., 2021), which showed relatively large increases in the numbers of bacteria associated with rumen digestion of metabolized starchy carbohydrates. This suggests that it may be necessary for rumen microbiota to maintain survival in an acidic environment by increasing the relative abundance of the phylum Bacteroidota. Combined with the experimental design, it was inferred that Bacteroidota used abundant starch substrates to reproduce rapidly over a short period of time, resulting in an increase in the relative abundance of Bacteroidota. Proteobacteria was the third most abundant phyla (Golder et al., 2014; Metzler-Zebeli et al., 2016). Proteobacteria play an important role in the formation and fermentation of biofilms and the digestion of soluble carbohydrates (Pitta et al., 2016). In this study, the relative abundance of Proteobacteria increased with increasing concentrate level, which is consistent with the findings of Wang et al. (2020). The reason for this may be due to the higher crude protein level and lower fiber content in the concentrate feeds, so it is hypothesized that the relative abundance of the Aspergillus phylum increases with increasing protein level in the feeds and is negatively correlated with the fiber level in the feeds.

At the genus level, Christensenellaceae_R-7_group, Rikenellaceae_RC9_gut_group, and Methanobrevibacter were the dominant microbiota in the rumen. The NK4A214_group, Candidatus_SaccharimonasEubacterium]_ventriosum_group, Olsenella, Lachnospiraceae_NK3A20_group, Prevotella, and Ruminococcus were also detected. The Christensenellaceae_R-7_group is described family in the phylum Firmicutes (Waters and Ley, 2019) that mainly decomposes fibrous material (Evans et al., 2011); in this study, the Christensenellaceae_R-7_group was highest in the C50 group and lowest in the C80 group. The yaks in the C80 group were fed a high level of concentrate, low level of fiber in concentrate feeds, this may account for the lowest relative abundance of the Christensenellaceae_R-7_group in this group, which is consistent with the results of Sha et al. (2021). The Rikenellaceae_RC9_gut_group belongs to the Rikenellaceae, and most are able to ferment unabsorbed polysaccharides in the host intestine, producing short-chain fatty acids (SCFAs), such as acetic acid, propionic acid, and butyric acid (Su et al., 2014). The Rikenellaceae_RC9_gut_group was significantly higher in the C80 group than in the C50 group. Because of the increase in concentrate levels, the contents of acetic acid, propionic acid and butyric acid produced by rumen fermentation in ruminants increased. Therefore that the relative abundance of Rikenellaceae_RC9_gut_group tended to increase with the increased level of concentrate. The NK4A214_group belongs to the family Rumenococcaceae, which is rich in endo-1, 4-beta-xylanas, and cellulase genes, which play an important role in the degradation of cellulose and hemicellulose; as such, their relative abundance is related to the diet concentrate to forage ratio (Koike and Kobayashi, 2009; Biddle et al., 2013). A recent study shows that NK4A214_group is positively correlated with isobutyrate and isovalerate concentrations (Liu et al., 2019b). The relative abundance of NK4A214_group was significantly higher in the C80 group than in the other two groups, similar to the results of a previously study (Chen et al., 2021), that the relative abundance of this genus increased with increasing dietary concentrate level. This can also explain the high isobutyrate and isovalerate concentrations in C80 group.

The results of principal component analysis and OPLS-DA showed significant differences in the metabolome of rumen fluid among the three groups. The study of pathway topology determined that the main metabolic pathway is nitrogen metabolism. The synthesis of rumen microbial proteins is directly related to nitrogen metabolism in the rumen (Li et al., 2019). Appropriate energy and protein levels in the rumen are more conducive to the growth and reproduction of rumen microorganisms and increase the utilization of feed protein by microorganisms (Wang et al., 2016). In the arginine and proline metabolism pathway, arginine is degraded by arginase to ornithine, which is metabolized to produce γ-aminobuytric acid (GABA). The 4-acetamidobutanoic acid is a derivative of GABA, which is an important mediator of arginine and proline metabolism (Gao et al., 2017); GABA acts as an excitatory transmitter in the early stages of life and is involved in brain development (Hasegawa et al., 2018). In our study, 4-acetamidobutanoic acid was positively correlated with Rikenellaceae_RC9_gut_group, which belongs to Rikenellaceae and produces short-chain fatty acids (SCFAs) to provide more energy, which has a positive effect on arginine and proline metabolism. Proline is a cyclic amino acid that enters the peptide chain and undergoes hydroxylation to form 4-hydroxyproline, an important component of animal proteins. In ruminants, arginine and proline are among the amino acids rich in casein, which promotes the digestion of cellulose by rumen microorganisms and play an important role in protein synthesis, metabolism, and nutrition (Amos et al., 1971; Bruckental and Alumot, 1984; Wu et al., 2011). In this study, the highest contents of 4-acetamidobutanoic acid and proline were found in the C65 group. One possible explanation is that the C65 group produced higher concentrations of proline which is more favorable for rumen amino acid metabolism and organism protein synthesis. In addition, Candidatus_Saccharimonas is a potential beneficial bacterium, which is closely related to fiber degradation and can produce cellulase in the rumen to degrade fiber-like substances such as cellobiose, and is a typical class of fiber-degrading bacteria (Liu et al., 2019a), whose increased content can improve fiber digestibility (Patra and Yu, 2015), which can promote the digestion and utilization of rumen nutrients to a certain extent. Correlation analysis showed that proline was positively correlated with Candidatus_Saccharimonas; therefore, the C65 group increased the concentration of proline in the rumen by increasing the relative abundance of Candidatus_Saccharimonas.

Glyoxylate and dicarboxylate metabolism, along with glycine, serine, and threonine metabolism are involved in the regulation of energy metabolism and immune regulation in the body; both metabolic pathways are associated with organic acids and provide important energy metabolic precursors for entry into the citric acid cycle (Yu et al., 2017). The metabolites L-threonic acid and PC (17,2/18:5) are present in this metabolic pathway, and threonine is converted to glycin and serine, which are key metabolites for protein synthesis (Laura et al., 2019). Phosphatidylcholine (PC) is an important component of the cellular lipid bilayer and belongs to a glycerophospholipid. Choline is a major component of PC, and sufficient PC facilitates lipid transport (Shiau and Cho, 2002; Yeh et al., 2013). PC acts as a choline “storage “molecule, and choline can be metabolized to organic osmolytes via the glycine, serine, and threonine metabolism pathway (Athamena et al., 2011; Jiang et al., 2019).

Moreover, this metabolite is also involved in the glycerophospholipid metabolism pathway (Mohamad et al., 2015). Lipid metabolism includes glycerophospholipid metabolism and glycerolipid metabolism. Lipid and lipid-like molecules in gut microbiota are associated with lipid and lipid-like metabolites that may affect lipid metabolism by increasing substrates for energy metabolism in the liver and surrounding tissues (Bäckhed and Crawford, 2010). In the present study, we found that the concentration of L-threonic acid tended to increase and then decrease, with the highest concentration in the C65 group; the concentration of PC (17:2/18:5) increased linearly with increasing concentrate levels. Our results suggest that metabolites related to lipid metabolism play an important role in the growth and rumen health of yaks. Increasing concentrate level increases the concentration of threonine in the rumen of yaks and promotes glycine, serine, and threonine metabolism; in turn, this improves the growth performance of yaks.

Arginine is a conditionally essential amino acid that plays an important role in the physiological function and metabolic regulation of animal organisms. Arginine promotes urea synthesis, regulates ammonia metabolism through the urea cycle, and improves reproductive performance, lactation function, immunity, growth performance, meat quality, and feed utilization (Chacher et al., 2013; Yohan et al., 2016; Zy et al., 2020). Spermidine is a derivative of arginine produced through the ornithine metabolism pathway and is a polyamine with important physiological functions (Minois et al., 2013). In the present study, the concentration of spermidine showed a linear increase with increasing concentration level, and was negatively correlated with Rikenellaceae_RC9_gut_group. This indicates that the concentration of arginine is affected by an increasing level of concentrate; therefore, raising the level of concentrate can promote arginine biosynthesis and contribute to the growth performance of yaks.

The results of this study show that there is a complex relationship between yak rumen microbes and their metabolites, both of which are influenced by diet concentrate to forage ratios. The results showed that different dietary concentrate ratios had major effects on rumen fermentation, microbial composition and metabolic functions in yaks. According to the results of this experiment, we found that the increase of concentrate level had a certain promotion effect on TVFA content in yak rumen and reduced the relative abundance of bacteria related to fiber degradation, increased the concentration of threonine in yak rumen and promoted the metabolism of glycine, serine and threonine, and the growth of yak was improved, the high ratio of forage to concentrate in the diet had negative effects. In summary, this study combined microbiomic and metabolomic analysis of the correlations among microbiota and differential metabolites in the rumen of yaks under different dietary concentrate ratios. The results offer new insight into yak rumen microbiota and metabolites, and provide guidance for optimal yak diet concentrate ratios to improve the growth performance of housed yaks.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement

All procedures in this study were approved by the Institutional Animal Care and Use Committee of Qinghai University.

Author contributions

SW, SC, and KP contributed to conception and design of the study. KP, SW, and XW collected the samples. SW, KP, YY, and SC conducted relevant experiments. KP and SW organized the database. XW and SL performed the statistical analysis. KP wrote the first draft of the manuscript. XW, YY, SL, and SC wrote sections of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by the Qinghai Provincial Science and Technology Department of China (grant no. 2020-ZJ-935Q) and Qinghai Province Thousands of High-end Innovative Talents Plan.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.