Effect of Dietary Forage to Concentrate Ratios on Dynamic Profile Changes and Interactions of Ruminal Microbiota and Metabolites in Holstein Heifers

A better understanding of global ruminal microbiota and metabolites under extensive feeding conditions is a prerequisite for optimizing rumen function and improving ruminant feed efficiency. Furthermore, the gap between the information on the ruminal microbiota and metabolites needs to be bridged. The aim of this study was to investigate the effects of a wide range of forage to concentrate ratios (F:C) on changes and interactions of ruminal microbiota and metabolites. Four diets with different F:C (80:20, 60:40, 40:60, and 20:80) were limit-fed to 24 Holstein heifers, and Illumina MiSeq sequencing and gas chromatography time-of-flight/mass spectrometry were used to investigate the profile changes of the ruminal microbes and metabolites, and the interaction between them. The predominant bacterial phyla in the rumen were Bacteroidetes (57.2 ± 2.6%) and Firmicutes (26.8 ± 1.6%), and the predominant anaerobic fungi were Neocallimastigomycota (64.3 ± 3.8%) and Ascomycota (22.6 ± 2.4%). In total, 44, 9, 25, and 2 genera, respectively, were identified as the core rumen bacteria, ciliate protozoa, anaerobic fungi, and archaea communities across all samples. An increased concentrate level linearly decreased the relative abundance of cellulolytic bacteria and ciliates, namely Fibrobacter, Succinimonas, Polyplastron, and Ostracodinium (q < 0.05), and linearly increased the relative abundance of Entodinium (q = 0.04), which is a non-fibrous carbohydrate degrader. Dietary F:C had no effect on the communities of anaerobic fungi and archaea. Rumen metabolomics analysis revealed that ruminal amino acids, lipids, organic acids, and carbohydrates were altered significantly by altering the dietary F:C. With increasing dietary concentrate levels, the proportions of propionate and butyrate linearly increased in the rumen (P ≤ 0.01). Correlation analysis revealed that there was some utilization relationship or productive association between candidate metabolites and affected microbe groups. This study provides a better understanding of ruminal microbiota and metabolites under a wide range of dietary F:C, which could further reveal integrative information of rumen function and lead to an improvement in ruminant production.


INTRODUCTION
The rumen is a complex microbial ecosystem and is inhabited by a high density of resident microbiota, with bacteria, ciliate protozoa, anaerobic fungi, and archaea involved (Hobson and Stewart, 1997). Rumen microorganisms play important roles in the degradation of feedstuffs and the production of volatile fatty acids (VFAs), lactate, amino acids, lipids, and hydrogen, which are crucial to the maintenance, growth, and production performance of ruminants (Shabat et al., 2016). Due to the complexity of the rumen ecosystem, a better understanding of the ruminal microbiota under extensive feeding conditions is a necessary prerequisite for the manipulation of the rumen microbiota to optimize rumen function and improve feed efficiency (Jami and Mizrahi, 2012;Bannink et al., 2016).
Although the effect of different forage to concentrate ratios (F:C) on rumen microbiota has been investigated, only two or three levels of F:C were used, which does not provide a dataset that is large enough to understand the microbiota variation (Saleem et al., 2012;Mao et al., 2016;Hua et al., 2017). When a high level of concentrate (usually about 50-65% concentrate on a dry matter basis) is used in diets in a sudden change, dairy cows are at a high risk of suffering from subacute ruminal acidosis (SARA), which results in a high risk of chronic inflammation, lameness, mastitis, and laminitis (Zebeli and Metzler-Zebeli, 2012;Mao et al., 2016;Hua et al., 2017). However, limit-feeding dairy heifers with high-concentrate diets was proposed as a potential strategy to improve feed efficiency in dairy farms (Zanton and Heinrichs, 2007 and can serve as an excellent model for examining the effects of a wide range of F:C on the ruminal microbiota. In most cases, limit-feeding in heifers refers to heifers that are fed with high-concentrate diets to meet, but not exceed, energy requirements for an optimal level of average daily gain Heinrichs, 2009, 2016). The importance of the microbiota in the rumen function has been recognized for decades; however, the dynamic changes of the ruminal microbiota and metabolic profiles in animals limit-fed high-concentrate diets remain unclear.
As cows' health, growth, and metabolism are largely dependent on metabolite production in the rumen, a comprehensive analysis of the composition of rumen fluid could offer important insights into the role of rumen-diet interactions on integrated cows' conditions (Saleem et al., 2013).
About 55-60% of rumen fluid metabolites were related to the rumen microbiota (Saleem et al., 2013). Henderson et al. (2015) established a varying rumen microbiome across a wide range of host species and dietary conditions worldwide, showing a more prominent effect of the diet than of the host. However, there is a gap between the information on the rumen microbiome and that on microbial metabolism that needs to be bridged (Bannink et al., 2016). Metabolomics has become an emerging research area that quantitatively measures small molecular metabolites in biological samples (biofluids or tissues) using high-throughput approaches, such as nuclear magnetic resonance (NMR) (Zhao et al., 2014), liquid chromatography tandem-mass spectrometry (LC-MS) (Tian et al., 2016), and gas chromatography-mass spectrometry (GC-MS) (Hua et al., 2017). Gas chromatography time-of-flight/mass spectrometry (GC-TOF/MS) has been widely used in ruminant biofluid detection (Saleem et al., 2013;Sun et al., 2015;Mao et al., 2016) because of its high sensitivity in metabolite analysis. The identification and integrative analysis of these metabolites can enable a comprehensive characterization of metabolism mechanisms at the molecular level.
Nowadays, the more recently developed high-throughput MiSeq sequencing technique has been successfully used to reveal the changes and functions of the ruminal bacteria community (Mao et al., 2015;Jin et al., 2016;Ye et al., 2016). An integrated investigation of the compositions and abundances of ruminal bacteria, ciliates, fungi, and archaea could give a better understanding of the microbiota in heifers limitfed with high-concentrate diets. We hypothesized that a wide range of dietary F:C would affect the ruminal microbiota and metabolites. The objective of this study, using four dietary levels of F:C, was to evaluate both linear and non-linear effects of increasing the proportion of dietary concentrate on ruminal profiles of microbiota (bacteria, protozoa, anaerobic fungi, and archaea) and metabolites in limit-fed Holstein heifers through a combination of the MiSeq sequencing strategy and the GC-TOF/MS technique. Moreover, the relationships between ruminal microbiota abundance and metabolites were also analyzed in the present study.

Animals, Diets, and Experimental Design
Twenty-four half-sib Holstein heifers (8-10 months and 263 ± 30 kg of body weight) with similar body condition from the Sanyuan Dairy Group (Beijing, China) were randomly assigned into four groups and fed diets with four different F:C (80:20, 60:40, 40:60, and 20:80; namely groups C20, C40, C60, and C80, respectively) with corn silage as the sole forage (Table S1). Half-sib heifers were used here to minimize the difference of individual animals. All heifers were fed with an adaptation diet containing 50% forage and 50% concentrate (Table S1) for 4 weeks, and then transferred to treatment diets that were formulated to meet the requirements of average daily gain at 800 g for dairy cattle (NRC, 2001). The amount of feed offered was adjusted weekly based on body weight, and the quantity of diets provided to high-concentrate-fed groups was restricted so that the intake of metabolizable energy (ME) was similar among groups. Heifers were individually fed with total mixed ration twice daily at 12 h intervals (0700 and 1900 h). Dry matter intake (DMI) for each individual heifer was recorded daily. Additionally, diets were formulated to maintain similar crude protein (CP):ME ratios in order to minimize potential effects due to differences in CP intake across diet treatments. Water was available ad libitum. Heifers were housed in individual stalls in a ventilated, environmentally controlled tie-stall barn with rubber mattress bedding, and were allowed access to an exercise lot for 2 h on days when sampling was not occurring.
This study was carried out in accordance with the recommendations of Instructive Notions with Respect to Caring for Experimental Animals, Ministry of Science and Technology of China. The protocol was approved by the Ethical Committee of the College of Animal Science and Technology of China Agricultural University.

Rumen Sample Collection and Measurements
At the end of the feeding period (day 28), rumen content (liquid and small particles; 100 mL/heifer) was collected using an oral stomach tube, approximately 4 h after morning feeding, according to a previously reported procedure (Shen et al., 2012). The device was thoroughly cleaned using fresh warm water between sample collections, and the first fraction of the sample from each cow was always discarded to ensure that sample collected was identified for each individual cow without contaminants from previous animals and her own saliva. Two mL of each rumen sample was immediately frozen in liquid nitrogen and stored at −80 • C to minimize any possible microbial activities for use in later DNA extraction. Another 2 mL of each rumen sample was passed through a disposable 0.22 µm filter, and after freezing in liquid nitrogen, samples were stored at −80 • C for future ruminal metabolomic analysis. Rumen pH was immediately measured after collection by a mobile pH meter (Starter 300; Ohaus Instruments Co. Ltd., Shanghai, China). After pH measurement, samples were passed through four layers of sterile cheesecloth and kept on ice until further processing. Filtered rumen fluid was centrifuged (17,000 g for 30 min at 4 • C) to obtain a clear supernatant, which was further analyzed for NH 3 -N using a phenol-hypochlorite assay (Broderick and Kang, 1980). Freshly prepared metaphosphoric acid (25% w/v; 2 mL) was added to 8 mL of filtered rumen fluid, vortexed, and then centrifuged (17,000 g for 10 min at 4 • C). The supernatant was stored at −20 • C for further measurement of VFAs. The VFA concentrations were determined with a capillary column gas chromatograph (Cao et al., 2008).
The sequences obtained from the MiSeq platform were processed through open-source software pipeline QIIME (Quantitative Insights into Microbial Ecology) version 1.8.0-dev (Caporaso et al., 2010b), with the criteria as described by previous reports (Mao et al., 2015;Ye et al., 2016). Briefly, (1) the reads that had a mean quality score of no <20 and no shorter than 50 bp were retained; (2) discarding occurred for reads that had exact barcode matching, two nucleotide mismatches in primer matching, and 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 (version 7.1, http://drive5.com/uparse/), and chimeric sequences were identified and removed by UCHIME (Edgar, 2010). The most abundant sequence within each OTU from specific libraries (bacteria, protozoa, fungi, and archaea) was designated as the "representative sequence" and aligned against the SILVA bacterial database (version 119) (Pruesse et al., 2007), NCBI-nt protozoa database (Koetschan et al., 2014), Unite fungi ITS database (version 7.0) (Koljalg et al., 2005), and SILVA archaea database using PYNAST (Caporaso et al., 2010a), respectively, with the default parameters set by QIIME. Community richness and diversity, analyzed with measures, such as Good's coverages, observed species, PD whole tree, Pielou, 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).

Preparation of Rumen Sample for GC-TOF/MS, and Identification of Compounds Discovered
The samples of rumen fluid were prepared, detected, and identified using the procedure previously described by Sun et al. (2015), with modifications. Adonitol was used as an internal standard and GC-TOF/MS analysis was performed using an Agilent 7890 gas chromatograph system coupled with a Pegasus HT TOF/MS (LECO, St. Joseph, MI, USA). The system had a DB-5MS capillary column (30 m × 250 µm inner diameter, 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, 270, 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.
The LECO-Fiehn Rtx5 database (Kind et al., 2009) and commercial databases, including KEGG (http://www.genome.jp/ kegg/) and NIST (http://www.nist.gov/index.html), were utilized to further identify and validate different metabolites. The resulted data involving the peak number, sample name, and normalized peak area were fed into the SIMCA software package (V14, Umetrics AB, Umea, Sweden) for principal component analysis (PCA), and orthogonal projections to latent structuresdiscriminant analysis (OPLS-DA). OPLS-DA was used here to obtain maximal covariance between measured data and response variable. The parameters of variable importance projection (VIP) value >1.5, and q (false discovery rate, FDR) < 0.05, were used as criteria to identify different metabolites between every two dietary treatment groups. In order to find the varying patterns of different metabolites with increasing dietary concentrate levels, the linear, quadratic, and cubic effects among treatments were evaluated among all different metabolites from six comparisons.

Data Analysis
Statistical analysis was performed using R software packages. Data of microbial abundance were transformed to log 10 (n + 1) if necessary to ensure normal distribution (Jin et al., 2016). The dietary effects were tested using the Estimability package (https://cran.r-project.org/package=estimability; author, R. V. Lenth; published date, 2015; version, 1.1-1) in R. The lm function of the Estimability package in R was used to evaluate the linear, quadratic, and cubic effects of the dietary concentrate levels on the variables. Standard errors of the mean were reported. Significance was declared at P < 0.05. Also, all P values from the multiple comparison analyses of the ruminal microbial community and the metabolomics data were adjusted by the FDR by the BonEV package (https://cran.r-project.org/package= BonEV; author, D. Li; published date, 2016; version, 1.0) with p. adjust in R. FDR-corrected P values below 0.05 (q < 0.05) were considered significantly different. DMI (5.32,4.97,4.69,and 4.42 kg/day) quadratically decreased (P = 0.02) with increasing dietary concentrate levels. As the level of dietary concentrate increased, rumen pH, ratio of acetate:propionate, and the proportion of acetate linearly decreased (P < 0.01), while the proportions of propionate and butyrate linearly increased (P < 0.01), the concentration of NH 3 -N quadratically increased (P < 0.01) with C80 at the highest value and C40 at the lowest value, the proportion of isobutyrate quadratically increased (P = 0.02) with C40 at the highest value and C80 at the lowest value, and the proportion of isovalerate cubically changed (P = 0.04) with C60 at the highest value and C40 at the lowest value (Table 1).

Ruminal Bacteria Changes within the Four Treatment Groups
In total, 1,888,443 raw reads were obtained for bacterial 16S rRNA genes by the sequencing procedure. After screening with strict criteria (as described in the methods), 1,797,207 valid reads were obtained, accounting for 95.2% of the raw reads. The Good's coverages for all samples were more than 99%. Compared with other groups, the C80 group had the lowest observed species, PD whole tree, and Pielou (P < 0.05; Table 2). Similar to C60 and C40 groups, the C80 group supported significantly less richness (P < 0.05) than that of the C20 group based on Chao1 index ( Figure 1A). Shannon index analysis indicated a similar tendency of diversity ( Figure 1B). Significant differences were found in the microbial communities at the OTU level between C80 and C40 groups (AMOVA, P = 0.03), between C80 and C20 groups (P < 0.01), between C60 and C20 groups (P < 0.01), and between C40 and C20 groups (P = 0.04). There was no significant difference found in the microbial compositions between C80 and C60 groups (P = 0.12), and between C60 and C40 groups (P = 0.46). The results of PCoA with weighted uniFrac distances indicated that the four treatment groups were largely separated from each other at the OTU level (Figure 2A). There were 2,095 OTUs identified in total, among which 1,352 OTUs were found in all four treatment groups and accounted for 64.5% of the total OTUs, revealing the presence of a strong common microbiota ( Figure S1A). OTU2, which was classified into the family Succinivibrionaceae, was the most dominant OTU taking up 6.73% of the total OTUs.

Ruminal Ciliate Protozoa Changes within the Four Treatment Groups
A total of 1,311,850 sequences were derived from the ciliate protozoa 18S rRNA genes by sequencing, with 1,311,783 used for analysis after quality control at 54,658 sequences per sample on average. The Good's coverages for all samples were more than 99%. There were no significant differences found in observed species, PD whole tree, and Pielou among treatments ( Table 2). The Chao1 value was indicative that the C80 group had lower (P < 0.05) richness than that of the C60 and C20 groups, but the richness in the C40 group was similar with that of other three groups ( Figure 1C). However, the Shannon index values were suggestive that ciliate community diversities were similar among the four groups ( Figure 1D).
Significant differences were found in the ciliate communities at the OTU level between C80 and C40 groups (AMOVA, P = 0.02), between C80 and C20 groups (P < 0.01), and between C60 and C20 groups (P < 0.01). There were no significant differences found in the ciliate compositions between C80 and C60 groups (P = 0.10), between C60 and C40 groups (P = 0.11), and between C40 and C20 groups (P = 0.14). The PCoA with weighted uniFrac distances also showed that C80 was separated from the C40 and C20 groups at the OTU level ( Figure 2B). There were 2,782 OTUs identified in total, and 641 OTUs were found in all four treatment groups (23.04%), which represented the common ciliates ( Figure S1B). The ciliate sequences detected at the phylum level largely belonged to Alveolata (>99.9%). There were 34 taxa identified at the genus level through analysis of microbiota composition, and nine genera were present in all samples, which were indicative of the core protozoal community in this study ( Figure 3B). Entodinium were the most abundant group among ciliates detected at the genus level, with a mean relative abundance of 72.3 ± 4.3%, followed by the second most abundant genus group of Ostracodinium (15.0 ± 3.2%) (Figure S3B). With an increasing amount of concentrate in diets, the prevalence of Entodinium linearly increased (q = 0.04), while the relative abundance of Ostracodinium (q = 0.04) and Polyplastron (q = 0.04) linearly decreased ( Table 3).

Ruminal Anaerobic Fungal Changes within the Four Treatment Groups
In total, 1,447,310 valid sequences were obtained for the anaerobic fungal ITS1 genes by sequencing, and 1,231,562 valid reads were obtained after screening, accounting for 85.1% of the raw reads. The Good's coverages for all samples were more than 98%. There were no significant differences found in observed species, PD whole tree, and Pielou among treatments ( Table 2). An increasing amount of dietary concentrate (P > 0.05) did not affect the diversity and richness of the anaerobic fungal communities, shown by the Chao1 and Shannon indices (Figures 1E,F).
There were significant differences found in the microbial communities at the OTU level between C80 and C20 groups (AMOVA, P = 0.02), and between C60 and C20 groups (P < 0.05). The PCoA axis 1 accounted for 52.70% of the variation, and the PCoA axis 2 accounted for 12.21% of the variation (Figure 2C). In total, 7,047 OTUs were identified, and 1,441 OTUs were found in all four treatment groups (20.45%), which represented the common anaerobic fungi in the rumen ( Figure S1C).
There were 335 genera identified at the genus level, and 25 genera were present in all samples, which were indicative of the core fungi community in this study ( Figure 3C). Piromyces, Caecomyces, Neocallimastix, Cyllamyces, and Anaeromyces were the predominant genera, with mean relative abundances of 17.9 ± 2.3%, 11.1 ± 2.1%, 9.3 ± 1.7%, 8.5 ± 1.8%, and 1.4 ± 0.3%, respectively ( Figure S3C). However, no genus was significantly affected by dietary concentrate levels (q > 0.10). Boxes represent the interquartile range (IQR) between the first and third quartiles (25 and 75th percentiles, respectively), and the horizontal line inside the box defines the median. Whiskers represent the lowest and highest values within 1.5 times the IQR from the first and third quartiles, respectively. Boxes with the different letters above their whiskers are significantly different (P < 0.05) among treatments. C20, diet containing 20% of concentrate; C40, diet containing 40% of concentrate; C60, diet containing 60% of concentrate; C80, diet containing 80% of concentrate.

Ruminal Archaea Changes within the Four Treatment Groups
After quality control, a total of 873,225 clean reads were derived from the archaea 16S rRNA genes by sequencing, with a mean length of 391 bp. The Good's coverages for all samples were more than 99%. Compared with other groups, the C80 group had the lowest observed species and PD whole tree (P < 0.05; Table 2). The C80 and C60 groups supported significantly lower richness (P < 0.05) than that of the C20 group, based on Chao1 index ( Figure 1G). No significant difference of archaea community diversities was found among treatments, based on the Shannon index ( Figure 1H).
A significant difference was only found in the microbial communities at the OTU level between C60 and C40 groups (AMOVA, P < 0.05). The results of PCoA with weighted uniFrac distances was indicative that the four treatment groups have a better separation at PCoA axis 1 than at PCoA axis 2 (73.81 vs. 17.07%, Figure 2D). In total, 1,627 OTUs were identified, and 331 OTUs were found in all four treatment groups (20.34%) which represented the common anaerobic fungi in the rumen ( Figure S1D).
The sequences detected at the phylum level largely belonged to Euryarchaeota (>86.3%), with only a few sequences belonging to Thaumarchaeota as the remaining phylum. Unclassified sequences comprised 13.65% of the total reads. There were seven taxa identified at the genus level through analysis of microbiota composition, and two genera were present in all samples, which were indicative of the core archaeal community in this study ( Figure 3D). The most abundant known genera was Methanobrevibacter (78.0 ± 3.7%) ( Figure S3D). Other minor genera, such as Methanosarcina, Methanosphaera, and Methanimicrococcus, accounted for 0.21 ± 0.15%, 0.18 ± 0.02%, and 0.05 ± 0.02%, respectively, with several reads belonging to Nitrososphaera, Methanocorpusculum, and Halostagnicola. The remainder of the sequences (21.60%) was unclassified at the genus level. However, no genus was significantly affected by dietary concentrate levels (q > 0.10).

Identification, Quantification, and Statistical Comparison of GC-TOF/MS Metabolites in the Rumen
As shown in the GC-TOF/MS total ion chromatograms of rumen samples from heifers fed diets with increasing levels of concentrate, there were in total 558 valid peaks identified from the four groups ( Figure S4). Based on the LECO/Fiehn Metabolomics Library, 243 metabolites, including amino acids, fatty acids, sugars, and organic acids, were quantified in the rumen from the four groups (Table S2), and the four groups shared the same metabolite categories.
Data of the GC-TOF/MS spectra among the four groups analyzed using PCA scores plots are listed in Figure S5. The parameters for the assessment of the OPLS-DA model quality in discriminating the four groups could be represented by validation plots (Figures 4A,C,E,G,I,K). The corresponding R 2 Y values of OPLS-DA models in C80 vs. C60, C80 vs. C40, C80 vs. C20, C60 vs. C40, C60 vs. C20, and C40 vs. C20 were 0.937, 0.952,  samples in the score plots were within the 95% of the Hotelling T 2 ellipse.

DISCUSSION
In this study, the effects of wide F:C on ruminal microbiota and metabolites, as well as the interaction between them, were investigated. By definition, limit-feeding always results in decreasing DMI with increasing dietary concentrate levels Heinrichs, 2009, 2016). Even though rumen pH decreased linearly as there was a larger amount of fermentable carbohydrate in high-concentrate diets, the lowest pH value was more than 6.3, which indicated that all the animals had a healthy rumen internal environment (Plaizier et al., 2017). In contrast to previous studies, which focused on the injury caused by highconcentrate diets and low ruminal pH (Saleem et al., 2012;Mao et al., 2016), our study focused on the healthy microbial composition and nutritive metabolites other than bioamine or phenylacetate.
Our findings showed that even with the limit-feeding strategy, the highest concentrate diet (80%) still had a negative effect on FIGURE 4 | Corresponding validation plots and OPLS-DA score plots derived from the GC-TOF/MS metabolite profiles of rumen samples for cows fed increasing levels of concentrate in their diets. Corresponding validation plots and OPLS-DA score plots (respectively) for: (A,B) the C80 group vs. C60 group; (C,D) the C80 group vs. C40 group; (E,F) the C80 group vs. C20 group; (G,H) the C60 group vs. C40 group; (I,J) the C60 group vs. C20 group; (K,L) the C40 group vs. C20 group. C20, diet containing 20% of concentrate; C40, diet containing 40% of concentrate; C60, diet containing 60% of concentrate; C80, diet containing 80% of concentrate. bacterial richness (based on the Chao1 value) and diversity (based on the Shannon index), which was similar with 70% concentrate free-feeding cattle (Mao et al., 2013). Bacteroidetes, Firmicutes, and Proteobacteria were the most predominant bacteria in the rumen, which was recognized by most previous studies (Jami et al., 2013;Henderson et al., 2015;Shabat et al., 2016). Prevotella, Succiniclasticum, Butyrivibrio, Ruminococcus, Fibrobacter, and Selenomonas were both the core bacterial community and main bacteria in this study, which was similar to previous reports (Jami and Mizrahi, 2012;Li et al., 2012;Zhou et al., 2017). The higher abundance of these genera would enable their members to occupy various ecological niches and make full use of nutrients within the rumen (Jami and Mizrahi, 2012;Mao et al., 2015). Among them, Butyrivibrio, Ruminococcus, and Fibrobacter are known to be cellulolytic bacteria (Russell and Baldwin, 1978;Fernando et al., 2010;Zhou et al., 2017), whereas Prevotella is known for the utilization of starch and protein (Jami et al., 2013). Succiniclasticum and Selenomonas are specialized in fermenting succinate and converting it to propionate (van Gylswyk et al., 1997;Fernando et al., 2010;Shabat et al., 2016). The existence of these dominant bacteria was considered to be the result of the large amount of fiber (the mean neutral detergent fiber level was 31.1% dry matter; mostly from corn silage) and carbohydrate (the mean non-fibrous carbohydrate level was 42.8% dry matter; mostly from corn silage and steam-flaked corn) resources coexisting in these diets. Within the bacteria community that has been significantly changed by different dietary concentrate levels, Acidithiobacillus, Christensenella, and RFN54, had low relative abundances (<0.05%), revealing that their affects and functions might be minor.
In line with a number of previous reports (Sylvester et al., 2004;Kittelmann et al., 2013;Mao et al., 2016), Entodinium was the most abundant group in the rumen and linearly increased with increasing amounts of dietary concentrate in the current study. Mao et al. (2016) and Hook et al. (2011) also reported that the percentage of Entodinium that was able to engulf starch increased as the proportion of grain in diets increased. Similarly with a report from Kittelmann et al. (2013) that Ostracodinium was largely dominant in silage-fed ruminants (86% of the sequences obtained from cows fed silage), we found that Ostracodinium was the dominant genera in the rumen and linearly decreased with the increasing proportion of dietary concentrate level.
In the present study, the anaerobic fungi richness (based on Chao1 value) and diversity (based on Shannon index) might have some resistance to high-concentrate diets, which FIGURE 5 | Correlations between the predominant ruminal microbiota affected by treatments and the potential marker compounds. Only the significantly affected ruminal microbiota (at the genera level) that had a mean relative abundance >0.5% were considered in the correlation. Only the ruminal metabolites that had VIP > 1.5, q < 0.05, similarity > 600, and fold change >2 or <0.5 were considered in the correlation. Bacteria are shown by red text color, and protozoa are shown by blue text color. (A) Correlation between the microbiota and metabolites. Cells are colored based on Spearman's correlation coefficient: red represents a positive correlation, and blue represents a negative correlation. (B) Co-occurrence network analysis among the microbiota and metabolites. Each co-occurring pair among microbial populations at the genus level and metabolites has an absolute Spearman rank correlation above 0.55 [gold line, positive correlation (r ≥ 0.55); blue line, negative correlation (r ≤ −0.55)] with a FDR-corrected significance level under 0.05. Microbes are shown by green round nodes, and metabolites are shown by pink diamond nodes.
was similar with the report by Mao et al. (2016), while Kumar et al. (2015) found that high-concentrate (50%)-fed post-calving cows had lower Chao1 value and Shannon index (P < 0.05) than in low-concentrate (20%) fed cows. That difference might be due to transition cows having more sensitive rumen microbiota than heifers. Consistent with previous reports (Liggenstoffer et al., 2010;Koetschan et al., 2014;Mao et al., 2016), Neocallimastigomycota was the most predominant fungal phylum and linearly decreased as the proportion of dietary concentrate increased in this study, which might be explained by the important role of Neocallimastigomycota in the degradation of fibrous plant materials in ruminants (Gruninger et al., 2014;Koetschan et al., 2014). Several ruminal fungi that were rare in previous studies (Kittelmann et al., 2013;Mao et al., 2016), namely Ascomycota, Zygomycota, Basidiomycota, Cercozoa, and Chytridiomycota, were found in the current study. This result might be explained by the primers used for sequencing. In two previous studies (Kittelmann et al., 2013;Mao et al., 2016), the same primers were used and similar fungi taxa were observed; likewise, the present study used the same primers as reported by Kumar et al. (2015), as a consequence, similar fungal taxa were observed in these two studies. As these phyla were mostly detected from marine samples (Jones et al., 2015) and rarely found in herbivores, their metabolic and functional significance in the ruminal ecosystem is unknown. Therefore, the cause of the linear increase of Ascomycota, Basidiomycota, Cercozoa, and Chytridiomycota in high-concentrate-fed animals remains unclear.
The six known and predominant anaerobic fungi genera, namely Anaeromyces, Caecomyces, Cyllamyces, Neocallimastix, Orpinomyces, and Piromyces, were all obtained from this study. This finding was in accordance with most previous studies (Liggenstoffer et al., 2010;Gruninger et al., 2014;Kumar et al., 2015). Anaeromyces, Cyllamyces, Neocallimastix, Orpinomyces, and Piromyces have cellulase and/or xylanase activities (Liu et al., 2008;Kittelmann et al., 2012;Jin et al., 2014). Neocallimastix and Orpinomyces also have carbohydrate active enzymes (CAZYmes) or genes Youssef et al., 2013). Therefore, the dominance of enzymes is likely dependent upon stimulation by the substrates provided. Kumar et al. (2015) reported that Cyllamyces remained unchanged between high (50%) and low (20%) concentration diets in primiparous cows, but the proportion of the genus Cyllamyces increased substantially (45%) in multiparous cows fed high-concentrate diets, while Caecomyces decreased in high-concentrate-diet fed cows, whether they were primiparous or multiparous (Kumar et al., 2015). Fernando et al. (2010) found Orpinomyces was abundant in cows fed a high-protein and hay-based diet. However, the coexistence of high-concentrate and high-forage diets among the four treatments might lead to a balance of the functions of cellulose and carbohydrate degradation, which resulted in the similar anaerobic fungi abundance at the genus level in the current study.
In this study, increasing levels of dietary concentrate decreased ruminal archaeal richness (based on Chao1 value), but had no significant effect on diversity (based on Shannon index), which was similar with results from Mao et al. (2016). In accordance with other studies (Kumar et al., 2015;Mao et al., 2016;Zhou et al., 2017), we also found that Euryarchaeota and Methanobrevibacter were the most abundant phylum and genus in the rumen, respectively. Methanobrevibacter was responsible for using hydrogen and/or formate to produce CH 4 in ruminants (Danielsson et al., 2017). Increasing dietary concentrate levels did not result in a changing archaea community (Kumar et al., 2015;Mao et al., 2016), which is similar to this study. It has been reported that archaeal communities have less variation and diversity across diets and species than other microbial communities (Jeyanathan et al., 2011;Henderson et al., 2015). There might be two reasons for that. On one hand the total density of archaea was quite low (3-4% of total microbes) in the rumen, and on the other hand the archaea community might have the ability to be resilient to dietary changes (Kumar et al., 2015). Therefore, more work is needed to illuminate the abundance and function of archaea in the rumen.
In order to assess the variations of the rumen microbiome and metabolism in response to increasing limit-fed dietary concentrate, the trends or regression lines over the four gradients of concentrate, rather than performing between group comparisons, were investigated. Among the significantly affected metabolites in the present study, most of them were amino acids, such as 3-aminoisobutyric acid, alanine, and proline. Also, some aromatic compounds, like phenylacetic acid (Scott et al., 1964) and hydrocinnamic acid (Moss et al., 1970), were the degradation products of amino acids. This result was indicative that the increasing levels of dietary concentrate had significant effects on ruminal amino acid metabolism.
According to Saleem et al. (2013), approximately 55-60% of the rumen fluid metabolites were related to the rumen microbiota. We investigated the functional correlation between common significantly different metabolites and the changed rumen microbiome in the present study. Acetate, propionate, and butyrate were the main compositions of VFAs and the end-products of numbers of microbiota (Kamke et al., 2016). Among them, propionate was the major VFA contributing to gluconeogenesis (Young, 1977), and butyrate was found as a key factor stimulating the initial growth and differentiation of the rumen to become a mature absorptive organ (Ploger et al., 2012;Kamke et al., 2016). Therefore, the higher percentage of ruminal propionate and butyrate in groups fed higher-concentrate diets could be beneficial to the energy metabolism and development of the rumen in heifers, even though the detection of morphological features of ruminal epithelium cells should be warranted.
In the process of synthesizing acetate, propionate, and butyrate, pyruvate was a competitive substrate for them, and the butyrate level could be raised as well as a higher propionate content occurring (Russell and Rychlik, 2001;Sutton et al., 2003;Ploger et al., 2012). As Anaerosporobacter, Saccharofermentans, and Ostracodinium produced acetate as main end-products (Petzel and Hartman, 1990;Ziemer, 2014;Ziganshina et al., 2014), it could be partly explained that these three bacteria were found to be negatively associated with propionate and butyrate, while positively associated with acetate in the present study ( Figure 5). Entodinium contributed a lot to propionate production in the rumen, which might be due to Entodinium being a starch degrader (Ivan et al., 2000). As acetate also was one of main end-products of Fibrobacter and Succinimonas or their symbiotic microbes (Forano et al., 2008), that could be the reason that acetate was positively related with these two bacteria in this study (r > 0.60, q < 0.01).
As NH 3 -N is one essential nitrogen source for Fibrobacter (Atasoglu et al., 2001), it is easy to understand the negative relationship between Fibrobacter and NH 3 -N that was found in this study (r < −0.71, q < 0.01). Since Fibrobacter and Succinimonas prefer to utilize fiber in high-forage diets and were less prevalent in high-concentrate diets in the present study, NH 3 -N was also negatively associated with Succinimonas ( Figure 5). The NH 3 -N concentration in the C80 group was higher than 5 mg/dL (Table 1), which was the minimum requirement for maximal microbial growth (Satter and Slyter, 1974); this was indicative that a high-concentrate diet could produce more ruminal microbe crude protein for utilization by the animals.
In ruminants, proline, which can be formed by arginine and glutamine/glutamate, was found to be one of the most abundant (approximately 12%) amino acids in casein and played important roles in protein synthesis, structure, metabolism, and nutrition (Bruckental and Alumot, 1984;Wu et al., 2011). When more proline was added into the diet, approximately 20% of arginine was saved (Alumot et al., 1983). Therefore, higher ruminal proline resulting from feeding high-concentrate diets could be beneficial to ruminal amino acid metabolism and whole-body protein synthesis. Also, proline was found to have the ability to enhance the digestion of cellulose by rumen microorganisms (Dehority et al., 1958;Amos et al., 1971). The fact that proline stimulates in vitro rumen microbial cellulose digestion may be explained partially by its conversion to valeric acid (Amos et al., 1971). In this study, proline was negatively associated with Fibrobacter (r ≤ −0.55, q < 0.02), which was suggestive that less proline was degraded in the rumen of heifers fed highconcentrate diets.
Among amino acids, peptides, and analogs, 3-aminoisobutyric acid was an intermediate of valine degradation in the pathways of valine, leucine, and isoleucine degradation (Podebrad et al., 2000). Valine was firstly degraded to isobutyryl-CoA and then to 3-aminoisobutyric acid, while branched chain fatty acid was also formed by isobutyryl-CoA (KEGG, Map 00280). Therefore, 3-aminoisobutyric acid and branched chain fatty acid were competitive in the rumen. Dehority et al. (1958) reported that valine could be oxidatively deaminated and decarboxylated by rumen microorganisms to form isobutyrate, and that isobutyrate had been shown to enhance the digestion of cellulose by rumen microorganisms in vitro. In the current study, 3-aminoisobutyric acid was negatively associated (r < −0.58, q < 0.05) with cellulolytic bacteria (such as Fibrobacter and Anaerosporobacter), as well as isobutyrate quadratically decreasing with increasing dietary concentrate levels (Table 1), which could be explained by most valine being degraded to isobutyryl-CoA then forming 3-aminoisobutyric acid other than isobutyrate in the rumen of high-concentrate-fed heifers.
Even though hydrocinnamate can be derived from phenylalanine (Moss et al., 1970), lactic acid bacteria have the enzymatic ability to convert cinnamate into hydrocinnamate (Broberg et al., 2007). The concentration of hydrocinnamate was found to be substantially higher in silage inoculated with lactic acid bacteria (Broberg et al., 2007). In the present study, the lower amount of corn silage (from 80 to 20% dry matter in diets, Table S1) used in high-concentrate diets could explain the linear decrease of hydrocinnamate in the rumen as the content of dietary concentrate was elevated, which was in line with work by Saleem et al. (2012). It was reported that hydrocinnamate could stimulate the growth of Ruminococcus albus 7, Ruminococcus flavefaciens FD-1, Butyrivibrio fibrisolvens 49, and Lachnospira multipara D-32, while the stimulation of the cellulolytic activity of these bacteria had species specificity (Borneman et al., 1986). This might be indicative that some species of Fibrobacter and Succinimonas could be specifically stimulated by hydrocinnamate in the current study.
Phytanic acid originates from chlorophylls (Dawson et al., 1974). One of the intermediate products, dihydrophytol, was hydrogenated and oxidized by rumen microorganisms to form phytanic acid, which led to a slow accumulation of phytanic acid in the ruminant's body and products (Hansen, 1966;Dawson et al., 1974). Hansen (1966) found that phytanic acid could account for 2.9% of the total fatty acids from the rumen bacteria of dairy cows fed a clover-grass hay diet. The bacteria, such as Anaerosporobacter, Fibrobacter, Saccharofermentans, and Anaeroplasma, that had the ability to degrade cellulose could release more chlorophyll and form phytanic acid in the present study.
In summary, the present study combined microbiome and metabolomics analysis to investigate the effects of limit-feeding four increasing gradient levels of dietary concentrate on ruminal microbial communities and metabolites under healthy conditions. With the increasing amount of dietary concentrate, the cellulolytic related bacteria and ciliate protozoa largely decreased, while other microbiota remained stable or slightly changed. Utilization or productive association was widely found between affected microbiota and certain metabolites. This study has allowed a better understanding of ruminal metabolites and microbial functions, which could further provide integrative information about rumen function and lead to improvements in ruminant production, such as increased digestion and feed efficiency.

AUTHOR CONTRIBUTIONS
YW, SL, ZC, JZ, and HS conceived and designed the experiments. HS, SJ, HZ, and JZ conducted the experiments. JZ, HS, SJ, and HZ performed the statistical analysis of the experimental data. Finally, the paper was written by JZ and HS, and was modified by YW and YH. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
This study was funded by a grant from the National Natural Science Foundation of China (grant no. 31402099), and the National Dairy Industry and Technology System (CARS-37). We are thankful for the help of workers from the Sanyuan Dairy Group.