Stability Assessment of the Rumen Bacterial and Archaeal Communities in Dairy Cows Within a Single Lactation and Its Association With Host Phenotype

Better characterization of changes in the rumen microbiota in dairy cows over the lactation period is crucial for understanding how microbial factors may potentially be interacting with host phenotypes. In the present study, we characterized the rumen bacterial and archaeal community composition of 60 lactating Holstein dairy cows (33 multiparous and 27 primiparous), sampled twice within the same lactation with a 122 days interval. Firmicutes and Bacteroidetes dominated the rumen bacterial community and showed no difference in relative abundance between samplings. Two less abundant bacterial phyla (SR1 and Proteobacteria) and an archaeal order (Methanosarcinales), on the other hand, decreased significantly from the mid-lactation to the late-lactation period. Moreover, between-sampling stability assessment of individual operational taxonomic units (OTUs), evaluated by concordance correlation coefficient (C-value) analysis, revealed the majority of the bacterial OTUs (6,187 out of 6,363) and all the 79 archaeal OTUs to be unstable over the investigated lactation period. The remaining 176 stable bacterial OTUs were mainly assigned to Prevotella, unclassified Prevotellaceae, and unclassified Bacteroidales. Milk phenotype-based screening analysis detected 32 bacterial OTUs, mainly assigned to unclassified Bacteroidetes and Lachnospiraceae, associated with milk fat percentage, and 6 OTUs, assigned to Ruminococcus and unclassified Ruminococcaceae, associated with milk protein percentage. These OTUs were only observed in the multiparous cows. None of the archaeal OTUs was observed to be associated with the investigated phenotypic parameters, including methane production. Co-occurrence analysis of the rumen bacterial and archaeal communities revealed Fibrobacter to be positively correlated with the archaeal genus vadinCA11 (Pearson r = 0.76) and unclassified Methanomassiliicoccaceae (Pearson r = 0.64); vadinCA11, on the other hand, was negatively correlated with Methanobrevibacter (Pearson r = –0.56). In conclusion, the rumen bacterial and archaeal communities of dairy cows displayed distinct stability at different taxonomic levels. Moreover, specific members of the rumen bacterial community were observed to be associated with milk phenotype parameters, however, only in multiparous cows, indicating that dairy cow parity could be one of the driving factors for host–microbe interactions.

Introduction

Ruminants live in a symbiotic relationship with a complex consortium of microorganisms residing in the rumen, responsible for the breakdown of feed ingredients, especially plant material rich in oligo- and polysaccharides (Flint et al., 2008). This functional feature results in the production of a wide range of microbial metabolites in the rumen, including short-chain fatty acids (SCFA), amino acids, vitamins, and methane. Some of these metabolites, in particular SCFA, represent a crucial source of nutrients to the host, contributing to the complex interaction between the rumen microbiota and host phenotypes, like milk production and composition, and feed efficiency (Vinet and Zhedanov, 2010; Jami et al., 2014; Lima et al., 2015).

Using culture-independent molecular approaches, the rumen microbiota composition has been described in a range of studies, with Firmicutes and Bacteroidetes as the dominant phyla (Pitta et al., 2014a, 2016; Jewell et al., 2015; Lima et al., 2015). Determining the driving factors for changes in the rumen microbiota composition of dairy cows has been an intense area of research, with diet as the most extensively investigated single factor (Henderson et al., 2015; Kumar et al., 2015; Nathani et al., 2015). Recent studies on large dairy cow and beef cattle cohorts have, however, identified heritable taxa of the rumen microbiota, suggesting that host genetics is also a determinant factor for rumen microbiota composition (Difford et al., 2018; John Wallace et al., 2019; Li et al., 2019). Findings like this have been a driver for hypothesizing the rumen microbiota to be composed of two key components; (i) a stable “core microbiota,” with conserved house-keeping features, crucial for functional stability of the rumen microbial community as well as the host and (ii) a “dynamic microbiota,” with an inherent plasticity allowing for adaption to swift changes in environmental stimuli, provided e.g., via the diet.

Furthermore, research characterizing the links between the variation in rumen microbiota and variation in host phenotypes, such as milk production and composition, feed efficiency, and methane emission has increased drastically. The Firmicutes to Bacteroidetes ratio has been reported to be strongly correlated with daily milk fat yield (Jami et al., 2014), and multiple bacterial taxa, including Micrococcaceae, Deltaproteobacteria, Erysipelotrichaceae, and Prevotellaceae, have been found to be correlated with milk production, and their presence to depend on cow parity and stage of lactation (Lima et al., 2015). Moreover, high feed-efficient cows have been reported to harbor a rumen microbiota with a higher relative abundances of Bacteroidetes, like Prevotella, and a lower relative abundance of Firmicutes (Delgado et al., 2019). Comparing the rumen microbiota composition of beef cattle varying in methane emission, revealed a higher Succinivibrionaceae abundance in low methane emitters (Wallace et al., 2015). Moreover, when testing different dosages of chloroform (an anti-methanogenic compound) to reduce methane emission, concurrent increase in the Bacteroidetes to Firmicutes ratio and a decrease in Archaea were observed (Martinez-Fernandez et al., 2016). These findings illustrate the need of cognizance of the rumen microbiome when evaluating different methane mitigation strategies on the host. Other parameters, like temporal rumen microbiota dynamics have received some research interest. For example, stage of lactation (Jewell et al., 2015; Zhu et al., 2018) and parity (Lima et al., 2015; Difford et al., 2018) have been reported to affect the rumen microbiota composition. However, temporal dynamics in rumen microbial composition e.g., caused by dietary interventions, and potential effects on associations between rumen microbiota and host phenotypes remains to be further investigated. For instance, shifts in the dairy cow rumen bacterial community over the transition period have been reported in a number of studies (Lima et al., 2015; Dieho et al., 2017; Zhu et al., 2018). In these studies, dietary effects, like switching from a high-forage prepartum to a high-concentrate postpartum diet, often dominates over host physiological status, thereby confounding dietary and host physiological changes with rumen microbial dynamics. In comparison, a relatively consistent dietary formula is typically provided to dairy cows during the lactation period, rendering a unique opportunity for obtaining insights into the temporal rumen microbial dynamics free of the dietary effects.

In the present study, we included a cohort of 60 lactating mixed-parity dairy cows, fed a standard total mixed ration, and investigated the rumen bacterial and archaeal dynamics over the lactation period by collecting rumen samples twice within a single lactation with an interval of 122 days. Rumen microbiome data were generated by 16S rRNA gene amplicon sequencing, and host phenotype data, such as milk yield, milk composition, and methane production were registered during the sampling period. We performed between-sampling stability assessment of individual bacterial and archaeal operational taxonomic units (OTUs) and investigated the associative patterns between the rumen microbial community members and selected host phenotype characteristics.

Materials and Methods

The animal experimental procedure followed a protocol approved by The Animal Experiments Inspectorate, Danish Veterinary and Food Administration, Ministry of Environment and Food of Denmark (Approval number 2016-15-0201-00959). Sixty lactating Holstein dairy cows (27 primiparous and 33 multiparous), between 13 and 320 days in milk (DIM), were housed in the same barn at a Danish commercial farm. The farm was equipped with an automated milking robot (Lely International N.V., Maassluis, Netherlands), and the cows were rewarded with up to 3 kg of concentrate at each milking. Throughout the lactation cycle, all the cows were fed ad libitum with the same total mixed ration, mainly constituted of grass silage (8.6 kg/21.5 kg), corn silage (6.8 kg/21.5 kg) and crashed barley (1.3 kg), and had free access to drinking water. All individuals were sampled twice, where the first and second samplings were carried out primo May and primo September, with a time interval of 122 days. A documented rumen sampling strategy was applied by inserting a rumen flora scoop via the mouth and esophagus into the rumen (Geishauser et al., 2012). Subsequently, approximately 40 mL rumen liquid mixed with fine particles were retrieved from the rumen, poured into a 50 mL centrifuge tube, immediately put into an icebox and transferred to the lab for further handling. Aliquots (five per sample) of 1.2 mL vigorously mixed rumen liquid were transferred into1.5 mL vials, snap frozen in liquid nitrogen, and stored at –80°C until further analysis. The remaining volume of the rumen samples was stored at –20°C for SCFAs analysis.

Milk Components Analysis and Methane Measurement

Milk yield of individual cows were recorded on the milking robot (Lely International N.V., Maassluis, Netherlands) upon each visit. On each rumen sampling day, samples of milk were withdrawn directly from the milking robot, and milk composition (protein and fat percentages) was determined using a CombiFoss (Foss, Hillerød, Denmark) at Eurofins (Vejen, Denmark). Methane production of the individual cows was measured with an on-farm methane measurement technique Gasmet DX-4000 (Gasmet Technologies, Helsinki, Finland) as described in Madsen et al. (2010) and Lassen and Løvendahl (2016). The raw data is provided in the Supplementary Material.

Short-Chain Fatty Acids Analysis

Rumen samples, stored at –20°C, were thawed at room temperature. After vigorously mixing for several minutes, 1 mL of sub-sample was adapted to SCFAs analysis according to a previously described protocol (Canibe et al., 2007). Determination of SCFAs was performed by gas chromatography, using a Hewlett Packard gas chromatograph (model 6890; Hewlett Packard, Agilent Technologies, Nærum, Denmark) equipped with a flame ionization detector and a 30-m ZB-5 column with an internal diameter of 0.32 mm and coated with 5%-phenyl 95%-di-methylpolysiloxane with a film thickness of 0.25 μm. The raw data is provided in the Supplementary Material.

Illumina Amplicon Sequencing and Data Analysis

All rumen samples stored at –80°C were handled in parallel and shipped on dry ice to a commercial sequencing company (GATC Biotech, Konstanz, Germany) for analysis. In order to minimize the bias in run difference, rumen samples were cross placed on 96 well plates and adapted to DNA extraction, PCR amplification and sequencing library construction. The universal bacterial 16S rRNA gene primer pairs (27 F: 5′-AGAGTTTGATCCTGGCTCAG-3′ and 534R or 518R: 5′-ATTACCGCGGCTGCTGG-3′) targeting the V1–V3 regions of 16S rRNA gene were used for studying the ruminal bacterial community. In a comparative study of different primer pairs, targeting different hypervariable regions of the bacterial 16S rRNA, this primer pair outperformed the others and provided higher species richness and diversity of the bacterial populations in the rumen (Pitta et al., 2014b). Likewise, the archaeal 16S rRNA gene primer pairs (519F: 5′-CAGCMGCCGCGGTAA-3′ and 1017R: 5′-GGCCATGCACCWCCTCTC-3′) targeting the V4–V6 regions of archaeal 16S rRNA gene were used for studying the ruminal archaeal community. One hundred twenty rumen samples from 60 cows were sequenced on the Illumina platform: for the bacterial 16S rRNA gene analysis, 36 samples were processed by Illumina HiSeq and 82 samples by Illumina MiSeq, while for the archaeal 16S rRNA gene analysis, 50 samples were processed by Illumina HiSeq and 70 samples by Illumina MiSeq. Two samples failed the bacterial 16S rRNA gene profiling. There was an inherent difference in sequencing output between the two Illumina sequencing platforms, with up to 700,000 reads produced from the Illumina HiSeq platform and up to 200,000 reads produced from the Illumina MiSeq platform. A normalization step was adopted to minimize the influence of two different sequencing platforms, and this was performed after the sequencing data analysis. The sequencing data analysis was performed with LotuS pipeline (Less OTUs scripts), developed by a research team in Belgium (Hildebrand et al., 2014). The sequence quality filtering, de-multiplexing and preparation was carried out with sdm options file. With the de-replication step of sequencing data, a minimum of 50 reads with exact 100% identity was applied for the bacterial sequencing data and a minimum of 10 reads with exact 100% identity was applied for the archaeal sequencing data. The total sequence length was truncated to 230 bp after barcode, adaptor, and primer removals. Sequences with homo-nucleotide over 8 were discarded. Middle quality sequences with minimum sequence length of 230 bp and minimum average quality score of 20 were used for estimating OTU abundance, while high quality sequences with minimum sequence length of 230 bp and minimum average quality score of 27 were used for OTU clustering. Filtered reads were clustered into OTUs using the UPARSE pipeline with 97% similarity (Edgar, 2013) and the OTU sequence was assigned taxonomy using RDP classifier (Wang et al., 2007). The Greengenes database (gg_13_5)¹ was used as a reference database for taxonomy assignment. A phylogenetic tree was constructed from sequences aligned with FastTree2 (Price et al., 2009). The OTU table and tree file from the Lotus pipeline were used as input for subsequent analysis in the QIIME-1.8.0 software package. The OTU table, generated from the bacterial 16S rRNA gene profile, was normalized to 76,000 sequences for each sample and the OTU table, generated from archaeal 16S rRNA gene profile, was normalized to 70,000 sequences for each sample with the code (single_rarefaction.py) prior to the calculation of alpha diversity index (chao1, observed species, and PD whole tree) and beta diversity metrics (weighted UniFrac distance matrix). Consequently, six samples were excluded from the large dataset, and the remaining 114 samples from 57 cows with successful characterization of both bacterial and archaeal 16S rRNA gene profile were kept for downstream analysis. All the sequencing data were deposited in the European Nucleotide Archive (ENA) under accession number PRJEB28065.²

Statistical Analysis

Based on the corresponding DIM of the cows, the 114 samples were grouped into three categories as follows: “early lactation (13–50 DIM),” “mid-lactation (50–100 DIM),” or “late-lactation (100–320 DIM).” A single-step multiple comparison procedure (Tukey’s honest significance test) was then applied to compare the relative abundance of bacterial phyla and archaeal orders across the lactation period. Significance was claimed at a P value < 0.05.

Lin’s concordance correlation coefficient (C value) is a scaled agreement index which is commonly used to assess the agreement between two measurement variables (NCSS, 2012). This index was adopted in the present study to evaluate the stability of ruminal bacterial and archaeal communities. Individual OTUs were assessed across all individual cows and the two sampling events. The calculation of concordance correlation coefficient was based on the following equation:

$$C=\frac{2\mathrm{\rho }{\mathrm{\sigma }}_{t1}{\mathrm{\sigma }}_{t2}}{{\mathrm{\sigma }}_{t1}^2+{\mathrm{\sigma }}_{t2}^2+{\left({\mathrm{\mu }}_{t1}-{\mathrm{\mu }}_{t2}\right)}^2}$$

where C, concordance correlation coefficient (between 0 and 1); ρ, Pearson correlation coefficient; σt₁², variance across all animals at the first sampling; σt₂², variance across all animals at the second sampling; μt₁, mean value across all animals at the first sampling; and μt₂, mean value across all animals at the second sampling. The concordance correlation coefficient (C value) was introduced as a parameter to evaluate the stability of individual OTUs, where C value >0.5 is denoted as “stable,” otherwise “unstable.”

Since replicate measures were taken per subject, a linear mixed model controlling for the random effect of subject was employed to test the significance of each bacterial and archaeal OTU fitted separately by means of the MIXED procedure in SAS (SAS 9.3, SAS Institute Inc.) as follows:

$$y_{ijkl}=\mathrm{\mu }+b_i(dim)+b_j\left(OTU\right)+P_k+a_l+e_{ijkl}(1)$$

where y_ijkl is the response variable (milk production, milk fat percentage, milk protein percentage, energy corrected milk, or methane production). μ is the model intercept. The term b_i is the fixed regression coefficient of the linear covariate days in milk (dim). Term b_j is the fixed regression coefficient of the natural logarithm of rarified OTU counts fitted for each bacterial and archaeal OTU separately. P_k is the fixed effect of the kth parity (k = 1–3). Term a_l is the random effect parameters for each animal a_l ∼ND (0, σ²_B) and e_ijkl is the random residual error terms with distributional properties of e_ijkl ∼ND (0, Iσ²_e). The significance values for term b_j where retained and negative log transformed for post hoc testing using a Bonferroni correction to control for false discovery rates, the results of which can be found in Figures 4–6.

Pearson’s product moment correlations between VFA components, milk traits and bacterial taxa, were calculated by means of CORR procedure in SAS (SAS 9.3, SAS Institute Inc.). With the same program, the correlation between archaeal taxa and methane yield was analyzed. Correlations between bacterial and archaeal taxa were inferred using SparCC program (Friedman and Alm, 2012) in mothur 1.35.1 (Schloss et al., 2009), where the permutations parameter was set as 1000. Co-occurrence network was constructed from those taxa, with pearson correlation coefficient (Pearson r) >0.5 or <–0.5 and with P value <0.05, by Cytoscape 3.5.1 (Shannon et al., 2003). Furthermore, a Paired t-test was applied to compare the alpha diversity index of ruminal bacterial and archaeal communities in primiparous and multiparous cows using SigmaPlot 11.0 (SYSTAT Software, Inc).

Results

Composition of the Ruminal Bacterial and Archaeal Communities Over the Lactation Period

The rumen bacterial and archaeal communities were represented by 6,363 and 79 OTUs, respectively, as characterized by 16S rRNA gene amplicon sequencing. The rumen bacterial community was dominated (values in parentheses represent mean relative abundances for the two samplings) by Bacteroidetes (74.9%), Firmicutes (14.5%) and Fibrobacteres (1.5%), accounting for more than 90% of the analyzed amplicons (Figure 1) and the relative abundances of these dominant phyla did not change significantly (P > 0.05) over the lactation period. The ratio between Bacteroidetes and Firmicutes (B/F ratio) differed between primiparous and multiparous cows over the lactation period, with a significant decrease in B/F ratio in primiparous cows (Slope –0.0104, P-value < 0.05, R² = 0.25) and a relatively constant B/F ratio in multiparous cows (Slope 0.0023, P-value > 0.05, R² = 0.0048) over the lactation period (Supplementary Material).

FIGURE 1

Figure 1. Phylum level composition of the rumen bacterial community. Bacterial 16S rRNA gene profile is shown in the area of the graph and distinct phyla are indicated by different colors. The x-axis represents days in milking (DIM) of cows ranging from day 13 to day 320 and the y-axis represents the relative abundance of bacterial taxa.

The relative abundance of the two least abundant phyla, SR1 (1.1%) and Proteobacteria (0.8%), on the other hand, decreased significantly (P < 0.05) from 1.2 and 1.3% in the mid-lactation to 0.9 and 0.6% in the late-lactation period, respectively. At the genus level, the rumen bacterial community was dominated by Prevotella (43.0%), Succiniclasticum (2.3%), Fibrobacter (1.5%), SR1_genera_incertae_sedis (1.1%) as well as a range of genera with a mean relative abundance below 1%. In comparison, the rumen archaeal community was less diverse, with Euryarchaeota being the most abundant phylum; within this phylum, five dominant orders, including Methanobacteriales (53.07%), Methanomassiliicoccales (42.73%), Methanosarcinales (1.69%), Methanomicrobiales, and Halobacteriales (both less than 1%), were identified; the least abundant two orders were only detected in a few rumen samples (Figure 2). The relative abundance of the Methanosarcinales order decreased significantly from 2.5% in the mid-lactation to 1.2% in the late-lactation period (P < 0.05). At the genus level, five dominant genera of the rumen archaeal community were shared by all the rumen samples, among which Methanobrevibacter (46.02%) was the most abundant genus, followed by unclassified Methanomassiliicoccales (27.78%), vadinCA11 (14.95%), Methanosphaera (6.28%), and unclassified Methanosarcinaceae (1.69%).

FIGURE 2

Figure 2. Order level composition of the rumen archaeal community. Archaeal 16S rRNA gene profile is shown in area graph and distinct orders are indicated by different colors. The x-axis represents days in milking (DIM) of cows ranging from day 13 to day 320 and the y-axis represents the relative abundance of bacterial taxa.

Diversity of the Rumen Bacterial and Archaeal Communities Over the Lactation Period

Diversity of the rumen bacterial and archaeal communities over the lactation period was assessed in multiparous and primiparous cows, separately. The diversity of the rumen bacterial and archaeal communities in primiparous cows as well as the archaeal community of multiparous cows showed no significant changes between the two samplings, however, the diversity of the rumen bacterial community in multiparous cows decreased significantly (P < 0.001) between samplings, as indicated by Chao1 and observed species indexes (Table 1).

TABLE 1

Table 1. Comparison of alpha diversity index in primiparous and multiparous cows between samplings.

OTU-Based Stability Assessment of the Rumen Bacterial and Archaeal Communities

Stability assessment of individual OTUs was evaluated by analyzing the concordance correlation coefficient (C-value), ranging from 0 to 1, with a C-value of 1 indicating complete stability; a C-value of 0.5 was chosen as the threshold between “stable” and “unstable.” The relative frequency of the 6363 bacterial OTUs with a C-value above 0.5 was 2.7 % (Figure 3A and Supplementary Figure 1). In other words, this assessment resulted in 176 bacterial OTUs (~2.7% of a total number of 6,363 bacterial OTUs) being stable over the course of the two samplings, whereas the remaining 6,187 bacterial OTUs were unstable. Of those 176 OTUs, 116 OTUs were detected in more than a half of the samples, with an average number of reads between 1 and 1,127 across all the samples, while the remaining 60 OTUs were detected in less than a half of the samples. The 176 stable OTUs were taxonomically assigned to Prevotella (36 OTUs), unclassified Prevotellaceae (27 OTUs), unclassified Bacteroidales (25 OTUs), unclassified Bacteroidetes (16 OTUs), unclassified Lachnospiraceae (15 OTUs), unclassified Clostridiales (9 OTUs), unclassified Alphaproteobacteria (4 OTUs), Ruminococcus (2 OTUs), and others (42 OTUs). In contrast, none of the 79 archaeal OTUs was stable over the course of the two samplings; all had C-values below 0.5 (Figure 3B).

FIGURE 3

Figure 3. Stability assessment of the rumen bacterial and archaeal communities. The calculation of concordance correlation coefficient (C value) is described in the statistical analysis section and percentages of OTUs with similar C value are displayed. An OTU with C value greater than 0.5 and P value of Pearson correlation coefficient smaller than 0.05 is denoted as “stable,” as indicated by red dash line, otherwise “unstable.” Consequently, the C value distribution of 6,363 bacterial OTUs (A) and 79 archaeal OTUs (B) are shown in bar graph.

Milk Phenotype-Based Screening Analysis of Bacterial OTUs

Linear mixed model-based screening analysis indicated that none of the 6,363 bacterial OTUs was significantly associated with milk production and energy corrected milk yield (Figures 4A,B); however, 32 OTUs were significantly associated with milk fat percentage (Figure 4C), and 6 OTUs were significantly associated with milk protein percentage (Figure 4D), all with a relative abundance below 0.5%. The 32 milk fat percentage-associated OTUs were taxonomically assigned to Lachnospiraceae (11 OTUs), Ruminococcaceae (2 OTUs), Porphyromonadaceae (2 OTUs), Prevotellaceae (1 OTU), unclassified Clostridiales (3 OTUs) and unclassified Bacteroidales (1 OTU), while the remaining OTUs were assigned to unclassified Bacteroidetes (10 OTUs) and unclassified Bacteria (2 OTUs). At genus level, 3 OTUs were further assigned to Butyrivibrio (1 OTUs) and Porphyromonas (2 OTUs). The six-milk protein percentage-associated OTUs were taxonomically assigned to Ruminococcus (1 OTU), unclassified Ruminococcaceae (1 OTU), unclassified Lachnospiraceae (1 OTU) and unclassified Bacteria (3 OTUs).

FIGURE 4

Figure 4. Scatter plot illustrating OTUs significantly associated with milk traits. Rarified sequence number of individual OTUs was used as input for the linear mixed model-based significance detection. Negative log base 10 transformed P value of individual OTUs (y-axis) was plotted against their corresponding OTU identities (x-axis). Bonferroni corrected false discovery threshold (0.05/6,363) is used as baseline for significance control as indicated by red line and individual OTUs are represented by open blue circle. The OTUs above the red baseline are considered as OTUs with significant effect on milk traits. (A) Milk production, (B) energy corrected milk, (C) milk fat percentage, and (D) milk protein percentage are shown in Manhattan plots.

Given the differences in milk quantity and quality between multiparous and primiparous cows, we separated their OTU profiles from each other, rendering a slightly higher number of bacterial OTUs in multiparous cows (6,363 OTUs) than in primiparous cows (6,355 OTUs). The need to perform the separated screening analysis was emphasized by the finding of 51 out of 6,363 bacterial OTUs, significantly associated with milk fat percentage (Figures 5A,B), and 3 out of 6,363 bacterial OTUs, significantly associated with milk protein percentage, in multiparous cows, whereas none of the 6,355 identified bacterial OTUs in primiparous cows showed any significant associations with milk phenotypes (Figures 5C,D). The 51-milk fat percentage-associated OTUs were taxonomically assigned to Lachnospiraceae (16 OTUs), Porphyromonadaceae (5 OTUs), Ruminococcaceae (2 OTUs), Prevotellaceae (2 OTUs), while the remaining 26 OTUs were mainly assigned to unclassified Bacteroidetes (14 OTUs) and others (12 OTUs). Among those OTUs, 7 OTUs were further classified as Porphyromonas (4 OTUs), Butyrivibrio (2 OTUs), and Prevotella (1 OUT). The three-milk protein percentage-associated OTUs were assigned to unclassified Ruminococcaceae, unclassified Alphaproteobacteria and unclassified Bacteria.

FIGURE 5

Figure 5. Comparison between multiparous and primiparous cows in OTUs significantly associated with milk traits. Rarified sequence number of individual OTUs was used as input for the linear mixed model-based significance detection. Negative log base 10 transformed P value of individual OTUs (y-axis) was plotted against their corresponding OTU identities (x-axis). Bonferroni correlation (0.05/6,363) is used as baseline for significance control as indicated by red line and individual OTUs are represented by open blue circle. The OTUs above the red baseline are considered as OTUs with significant effect on milk traits. Milk fat percentage of (A) multiparous cows and (B) primiparous cows, milk protein percentage of (C) multiparous cows and (D) primiparous cows are shown in scatter plots.

Methane Production-Based Screening Analysis of Archaeal OTUs

Linear mixed model-based screening analysis did not reveal any significant associations between the 79 archaeal OTUs and methane production (Figure 6).

FIGURE 6

Figure 6. Scatter plot illustrating bacterial and archaeal OTUs significantly affecting methane production (g/day). Rarified sequence number of individual OTUs was used as input for the linear mixed model-based significance detection. Negative log base 10 transformed P value of individual OTUs are shown in y-axis and their corresponding identities are shown in x-axis. Bonferroni correction is used as baseline for significance control, bacterial OTU (0.5/6363) (A) and archaeal OTU (0.5/79) (B), as indicated by red line and individual OTUs are represented by green dot.

Correlation Analysis Between SCFAs, Milk Traits and Genus Abundance

Correlation analysis was performed to evaluate the relationship between rumen microbial taxa, rumen SCFA profiles, and milk traits, and only those microbial taxa with a relative abundance of more/higher than 0.1% were included in the analysis. Furthermore, the taxa not showing any significant correlation to either of the tested parameters were excluded from the results (Figure 7). Thus, Saccharibacteria showed the most positive correlation (Pearson r = 0.38), whereas unclassified Ruminococcaceae showed the most negative correction to total SCFA concentration (Pearson r = –0.33), with a similar correlation pattern for the two genera to acetate concentration (Pearson r = 0.39 and -0.40, respectively). In addition, Saccharibacteria showed the strongest correlation to both acetate (Pearson r = 0.39) and propionate concentration (Pearson r = 0.43), and negative correlation to butyrate concentration (Pearson r = –0.66). Unclassified Ruminococcaceae and unclassified Veillonellaceae showed a positive correlation to butyrate concentration (Pearson r = 0.37 and 0.31, respectively). The remaining microbial taxa (e.g., unclassified Marinilabiliaceae, unclassified Porphyromonadaceae, Methanosphaera and Lachnospira) showed moderate correlations with individual SCFA proportions (Pearson r between –0.25 and 0.26). Similarly, several taxa (e.g., Hallella, Ruminococcus, Saccharibacteria, Saccharofermentans and unclassified Porphyromonadaceae) showed positive correlations to A/P ratio (Pearson r ranging from 0.19 to 0.28). Among all taxa, unclassified Veillonellaceae showed the most negative correlation to the A:P ratio (Pearson r = –0.43).

FIGURE 7

Figure 7. Correlation between SCFAs, milk traits and genus abundance. Pearson correlation analysis between dominant bacterial genera (with relative abundance greater than 0.1%) and concentrations (mmol/kg wet sample) and proportions (%) of SCFAs, milk components (milk yield, energy-corrected milk yield (ECM), milk fat percentage, and milk fat percentage), and methane production (g/day) across the two samplings was performed and the results are shown in Heat map. Only those bacterial taxa showing significant correlations (P < 0.05), positive (red) and negative (blue), with either of the listed items on the left side are shown.

Unclassified Prevotellaceae and Hallella (both Bacteroidetes) were positively correlated with milk yield (Pearson r = 0.21 and 0.33, respectively) and energy-corrected milk yield (Pearson r = 0.19–0.21, respectively), whereas Lachnospiraceae_incertae_sedis and unclassified Ruminococcaceae (both Firmicutes) were negatively correlated with milk yield (Pearson r = –0.22 and –0.25, respectively). Butyrivibrio showed the most positive correlation to milk fat percentage (Pearson r = 0.44), whereas Prevotella showed the most negative correlation to milk fat percentage (Pearson r = –0.22). Three Firmicutes genera (Butyrivibrio, Lachnospiraceae_incertae_sedis, and Pseudobutyrivibrio) showed positive correlations with milk protein percentage (Pearson r ranging from 0.23 to 0.30).

Co-occurrence Analysis of Bacterial and Archaeal Taxa

Complex interactions both between and within the rumen bacterial and archaeal communities were revealed by co-occurrence analysis (Figure 8). For example, Fibrobacter was positively correlated with vadinCA11 (Pearson r = 0.76); Ruminococcus, along with Lachnospiraceae_incertae_sedis and Fibrobacter, were positively correlated with unclassified Methanomassiliicoccaceae (Pearson r = 0.53, 0.63 and 0.64, respectively). Within the rumen archaeal community, both Methanosphaera and Methanobrevibacter were negatively correlated with vadinCA11 (Pearson r = –0.53 and –0.56); Methanobrevibacter was also negatively correlated with unclassified Methanomassiliicoccaceae (Pearson r = –0.54). In addition, betweenness centrality score of the co-occurrence network was used to identify keystone species; a high score of a taxon reflects the importance of control over the interactions of other taxons in the rumen microbial network. As a result, Saccharibacteria and Succiniclastium of the rumen bacterial community had the highest centrality score (1.0 and 0.50, respectively); similarly, unclassified Methanomassiliicoccaceae and Methanosphaera of the rumen archaeal community had the highest centrality score (0.27 and 0.12, respectively).

FIGURE 8

Figure 8. Co-occurrence network analysis between rumen bacterial taxa and archaeal taxa. Only correlations with coefficient >0.5 or <–0.5 and with P value <0.05 were displayed. The width of circle corresponds to the relative abundance of individual taxa. Different line colors are used to indicate a positive (red) or negative correlation (blue) between two taxa and the strength of correlation is indicated by the thickness of the lines.

Discussion

Dynamic Shift of the Rumen Bacterial and Archaeal Communities Over the Lactation Period

At the phylum level, we observed Bacteroidetes and Firmicutes to be dominating the rumen bacterial community and stable in their relative abundance across individual cows and irrespective of parity, DIM, and sampling time, in accordance with previous studies (Pitta et al., 2014a; Jewell et al., 2015; Lima et al., 2015). Likewise, Methanobacteriales and Methanomassiliicoccales were dominant orders of the rumen archaeal community, as also previously reported (Seedorf et al., 2015; Zhu et al., 2017). However, the significant decrease in relative abundances of two bacterial phyla (Proteobacteria and SR1) and one archaeal order (Methanosarcinales) between the two samplings, suggest that the rumen bacterial and archaeal communities (as evaluated on the phylum level) cannot simply be assumed to be completely stable over the lactation period.

The repeated sampling within a single lactation period allowed us to assess the stability of individual OTUs. The choice of the C-value threshold of 0.5 is arbitrary and motivated by being precisely in between perfect concordance (C-value = 1) and complete lack of concordance at (C-value = 0). Resultantly, our observation of unstable bacterial OTUs, along with decreases in the relative abundance of the two phyla, are in agreement with previous findings that the rumen bacterial communities of dairy cows, under the same management system, differed in their composition (Jewell et al., 2015; Bainbridge et al., 2016). Lactation stage has been considered as an important factor responsible for shifts in the rumen bacterial communities of lactating dairy cows (Jewell et al., 2015; Bainbridge et al., 2016) and transition cows (Pitta et al., 2014a; Lima et al., 2015; Dieho et al., 2017; Zhu et al., 2018), and likewise for parity (Pitta et al., 2014a; Jewell et al., 2015; Lima et al., 2015). Accordingly, the fact that the dairy cows involved in the present study varied in lactation stage and parity may account for the observed differences in the rumen bacterial community composition.

Likewise, individual archaeal OTUs-based stability assessment, along with changes in the relative abundance of the order Methanosarcinales, demonstrated an unstable methanogenic archaeal community over the lactation period. The rumen methanogen community is in general resilient to changes in dietary composition and all the animals were in our case maintained on the same diet over time. Diets with different forage-to-concentrate ratios (80:20, 60:40, 40:60, and 20:80) have been observed to have no significant effects on the rumen archaeal community structure at genus level, but significant differences were observed at OTU level between dietary treatments (60:40 versus 40:60) (Zhang et al., 2017). This information suggests that individual OTU-based analysis can provide different insights into the rumen methanogenic archaeal community. It has been reported that the rumen methanogen community was more influenced by DIM than breed differences (Cersosimo et al., 2016). More specifically, the mean relative abundance of Methanosphaera was highest around mid-lactation and lowest around early- and late-lactation, further iterating that the difference in DIM (lactation stage) of the host plays a critical role in the changes of the rumen methanogen community (Cersosimo et al., 2016). In addition, age-related compositional and functional investigations of the rumen methanogenic community revealed that methylotrophic methanogenesis was more active in the young animal rumen, whereas hydrogenotrophic methanogenesis was more active in the mature animal rumen (Friedman et al., 2017). Potential age- or parity-related changes in the rumen methanogen community should thus be taken into consideration. Again, we presented here an unstable rumen methanogen community of dairy cows over the lactation period.

Stable Feature of the Rumen Bacterial Community

Despite the shifts of the rumen bacterial community composition at the phylum and OTU level, we identified 176 stable bacterial OTUs, among which 36 OTUs were taxonomically assigned to Prevotella, representing the most prevalent genus in the rumen. Functionally, members of Prevotella, like Prevotella ruminicola (former Bacteroides ruminicola), have demonstrated abilities for xylan, starch, pectin, and peptide degradation (Dehority, 1966, 1969; Cotta, 1988; Kim et al., 2017). A recent multi-omics study revealed positive correlations between mainly Prevotella species and ruminal microbial metabolites related to amino acid and carbohydrate metabolism (Xue et al., 2020). The functional flexibility may allow them to adjust to the challenging rumen conditions. Additionally, it has been reported that the rumen microbiota can display substantial redundancy (overlap of function among multiple species/taxa) and resilience (resistant to perturbation) (Weimer, 2015). These features may be related to the genetic diversity of the Prevotella group, potentially rendering members of the genus functionally flexible (Avgustin et al., 1994). Further efforts are still required to isolate and characterize members of this taxon. Other stable OTUs were assigned to unclassified Prevotellaceae (27 OTUs), unclassified Bacteroidales (25 OTUs), and unclassified Lachnospiraceae (15 OTUs); these bacterial taxa have been suggested as part of the “core bacterial microbiome” in ruminants (Henderson et al., 2015). Multiple OTUs belonging to the Prevotellaceae, Lachnospiraceae, and Bacteroidales, were correlated with gross feed efficiency of dairy cows (Jewell et al., 2015), suggesting that the stable OTUs (bacterial species) we observed may possess important and conserved ecological and functional implications to the ruminant host.

Association Between Microbial Taxa, SCFAs and Milk Components

Our linear mixed model-based screening analysis revealed that the OTUs (32 in total), significantly associated with milk fat percentage were taxonomically assigned to Bacteroidetes [e.g., unclassified Bacteroidetes (10 OTUs)], as well as Firmicutes [e.g., unclassified Lachnospiraceae (9 OTUs) and Butyrivibrio (2 OTUs)]. This observation is in alignment with the results from correlation analysis between bacterial taxa abundance and milk fat percentage, such as Lachnospiraceae_incertae_sedis and milk fat percentage (Pearson r = 0.20), Butyrivibrio and milk fat percentage (Pearson r = 0.44), and Pseudobutyrivibrio and milk fat percentage (Pearson r = 0.28). These observations may be attributable to their abilities to produce specific SCFAs, as indicated by the positive correlations between Lachnospiraceae_incertae_sedis abundance and concentrations of acetate (Pearson r = 0.21) and propionate (Pearson r = 0.20). Acetate can be incorporated into a wide range of microbial fatty acids, and propionate and butyrate can serve as substrate for lipogenesis to some extent (Emmanuel et al., 1974). However, despite the dominance of Prevotella in the rumen microbiota and its positive correlations with total and individual SCFA concentrations, this genus showed the most negative correlation with milk fat percentage (Pearson r = –0.22). This observation is comparable to a previous finding that the abundance of Prevotella, was negatively correlated with milk-fat yield (Pearson r = –0.69) (Jami et al., 2014). On the contrary, other studies reported a positive correlation between Prevotella and milk fat percentage (Chiquette et al., 2008; Indugu et al., 2017). One explanation of this discrepancy could be due to the dilution effect of milk. We observed a positive correlation between Prevotella and milk yield, as also reported by another study (Indugu et al., 2017). A negative correlation between milk yield and milk fat percentage was found when we examined the relationship between these two parameters at host level (data not shown). This information suggests that an increase in the total quantity of milk may lead to a decrease in milk fat percentage. Another explanation of this discrepancy could be related to the genetic diversity of Prevotella, with the inclusion of different functional species occupying different metabolic niches (Bekele et al., 2010). Varied responses of individual genus members to dietary components may contribute to the negative correlation between Prevotella and milk fat percentage. The genus Pseudobutyrivibrio harbors butyrate-producing species (e.g. Pseudobutyrivibrio xylanivorans (Kopečný et al., 2003), and within the genus Butyrivibrio, Butyrivibrio fibrisolvens even plays a major role in regulation of rumen lipid metabolism (Kepler et al., 1966; Maia et al., 2010).

Linear mixed model-based analysis revealed that 6 OTUs, related to Ruminococcus, unclassified Ruminococcaceae, unclassified Lachnospiraceae and unclassified Bacteria, were significantly associated with milk protein percentage. Ruminococcus and unclassified Lachnospiraceae are not only abundant members of the rumen microbial community, but also play a very important role in fermenting diverse plant polysaccharides. For instance, two important cellulolytic members of the Ruminococcus genus, Ruminococcus flavefaciens and Ruminococcus albus, can degrade complex plant polysaccharides into oligosaccharides (Krause et al., 2003); this fundamental step in the rumen fermentation scenario may affect the growth and activity of other members of the rumen microbial ecosystem, thereby regulating the ruminal microbial biomass.

Interestingly, only two OTUs were shared between the milk phenotype-associated OTUs and the stable OTUs (Supplementary Figure 1), suggesting these OTUs (bacterial species) may fill different functional roles in the rumen microbial ecosystem. Most of the stable OTUs were taxonomically assigned to Prevotella and unclassified Prevotellaceae, and a multitude of correlations were detected between Prevotella, along with Lachnospira, Succiniclasticum, and Saccharofermentans, SCFAs, and milk traits. It is possible that these stable OTUs are fundamental for the rumen anaerobic food chain, while milk-phenotype associated OTUs may be functionally active in the rumen. A stable profile of individual SCFA proportions was observed across individual animals, parity and lactation stage (Supplementary Figure 2), indicating that overall, the rumen microbial activity was not affected by shifts in the rumen microbial community structure over the lactation period. These observations may be closely related to the capacities of the stable Prevotella OTUs (species) to degrade a wide variety of feed ingredients into SCFAs (Dehority, 1966; Cotta, 1988). On the other hand, pectin degradation by Lachnospira species (Rode et al., 1981) and glucose degradation by Saccharofermentans species (Chen et al., 2010) were also enhanced. Perhaps Prevotella species may act as primary degraders of feed ingredients, providing substrates for other functionally active bacteria species which may further contribute to milk phenotype. However, due to lack of cultured isolates, the understanding of the metabolic functions and networks of these, as well as a range of unclassified groups is still in its infancy. Therefore, taking advantage of type species or strains obtained from existing culture collections or novel culturomics approaches to gain mechanistic insights into the host-microbe interactions through more advanced computational approaches such as genome-scale metabolic modeling is highly encouraged (Perez-Garcia et al., 2016).

A large number of the milk phenotype-associated OTUs, identified from the entire dataset, were only found in the multiparous cows and not the primiparous cows, suggesting parity to be a driver of the interactions between the rumen microbiome and milk phenotype. Indeed, previous studies have reported significant differences in milk production (Wathes et al., 2007; Lang et al., 2012) as well as rumen bacterial community composition and function (Pitta et al., 2016) between primiparous and multiparous cows. Most of the milk phenotype-associated OTUs were taxonomically assigned to the phyla Bacteroidetes and Firmicutes. Bacteroidetes was reported to contribute to a majority of metabolic functions in primiparous cows while contributions of Firmicutes and Proteobacteria were more pronounced in multiparous cows (Pitta et al., 2016). This information can be supportive of our observation that B/F ratio in primiparous cows showed a decreased trend over the lactation period. Perhaps an increased proportion of Firmicutes as the primiparous cows advanced through the lactation is required to maximize the contributions of the rumen microbiota to milk phenotype. Difference in the lactation stage of the involved cows in our case is indeed a critical factor for determining stability of the rumen bacterial community. Moreover, the absence of milk phenotype-associated OTUs in primiparous cows is in accordance with previous findings, where more complex interactions, both within and between the microbial domains (as represented by anaerobic fungi, methanogens and bacteria), were observed in multiparous cows than in primiparous cows (Kumar et al., 2015). Therefore, future efforts are required to investigate the stability of the entire rumen microbiome of dairy cows over the lactation period.

Although methanogenic archaea are the only direct methane producers, we did not observe any significant correlation between relative abundance of the archaeal taxa and methane production. It has been reported that high and low methane-emitting sheep could be differentiated by rumen methanogen mRNA rather than the DNA level (Shi et al., 2014). Thus, methane production may be linked to the activity rather than the relative proportion of specific methanogens, advocating the involvement of transcriptomic-based techniques for future investigations of the rumen archaeal community (Poulsen et al., 2013).

Interaction Network Between Bacterial and Archaeal Taxa

Co-occurrence analysis revealed distinctive associative patterns between bacterial (Fibrobacter, Ruminococcus, and Lachnospiraceae) and archaeal taxa (unclassified Methanomassiliicoccaceae and vadinCA11), indicating their inter-dependence in the rumen microbial ecosystem. Fibrobacter succinogenes is important in rumen fiber digestion (Kobayashi et al., 2008), and this fundamental step may provide substrates to other bacteria that are capable of producing H₂ or methyl-containing compound. Other well-characterized bacterial members, such as Ruminococcus albus 7 (Christopherson et al., 2014) and Lachnospira multiparus (Rode et al., 1981), have the abilities to produce methyl-containing compounds and/or H₂ from fibrous substrates, thus favoring the growth of methylotrophic methanogens, as represented by members of the Methanomassiliicoccales order. Furthermore, network analysis revealed Succiniclasticum as a keystone member of the rumen microbiome. A known isolate from this group, Succiniclasticum ruminis (Van Gylswyk, 1995), is capable of converting succinate to propionate as the sole energy-yielding strategy. This step is of critical importance in the rumen fermentation network, and the resulting product propionate is an important precursor for gluconeogenesis in ruminants. Similarly, our observation that Saccharibacteria had the highest centrality score may be related to its positive correlations with total SCFAs, acetate, and propionate concentrations.

In conclusion, the results from the present study demonstrated the rumen bacterial and archaeal communities of dairy cows to be largely unstable over the lactation period, and unstable OTUs were most often significantly associated with host phenotypic variation. However, we identified 176 stable bacterial OTUs over the lactation period and these were mostly unrelated to host phenotypes, but may possess important conserved ecological and/or metabolic functions. Associations between milk phenotypes and OTUs were detected in the multiparous cows but absent in the primiparous cows as well as higher parity cows had reduced diversity, suggesting that increased milk production in higher parity cows could be associated with specialization of the rumen microbiome. Further research into the stable and unstable rumen microbiome features is needed to characterize associations between host phenotypic variation and changes in the rumen microbiota. Technically, transcriptomic-based approaches should certainly be included to elucidate functional stability or dynamics of the rumen microbiota, e.g., over the lactation period.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: European Molecular Biology Laboratory’s European Bioinformatics Institute (EMBL-EBI) European Nucleotide Archive (ENA), https://www.ebi.ac.uk/ena/browser/view/PRJEB28065.

Ethics Statement

The animal study was reviewed and approved by The Animal Experiments Inspectorate, Danish Veterinary and Food Administration, Ministry of Environment and Food of Denmark (Approval number 2016-15-0201-00959).

Author Contributions

PL, OH, and ZZ: conceptualization. PL, OH, and SN: supervision. JL and ZZ: rumen sample collection. PL, JL, and GD: host phenotype data creation. ZZ, GD, and SN: rumen microbiome analysis. ZZ: writing – original draft. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by Innovation Fund Denmark (formerly The Danish Council for Strategic Research) grant number 12-132447.

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.

Acknowledgments

We thank Bernt Guldbrandtsen for his kind help regarding the python scripts used for bioinformatics analysis of next generation sequencing data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.636223/full#supplementary-material

Supplementary Figure 1 | Venn diagram showing OTU distribution. Shared OTUs between milk phenotype-associated (MPA) OTUs (38 OTUs), stable OTUs (176 OTUs), and a total number of 6363 bacterial OTUs, are illustrated in Venn diagram, indicating that only 2 OTUs, assigned to unclassified Firmicutes and unclassified Bacteroidetes, are shared between the MPA and stable OTUs.

Supplementary Figure 2 | SCFA proportions over the lactation period. The relative proportions of individual SCFAs (A: Acetate; P: Propionate; B: Butyrate) were plotted against days in milking ranging from day 13 to day 320. Major SCFA components are designated as open circles with distinct colors. Repeatability of individual SCFA proportions (acetate, propionate and butyrate) are 0.19, 0.11, and 0.21, respectively.

Footnotes

1. ^ https://greengenes.secondgenome.com
2. ^ https://www.ebi.ac.uk/ena/data/view/PRJEB28065

References

Avgustin, G., Wright, F., and Flint, H. J. (1994). Genetic diversity and phylogenetic relationships among strains of Prevotella (Bacteroides) ruminicola from the rumen. Int. J. Syst. Bacteriol. 44, 246–255. doi: 10.1099/00207713-44-2-246

PubMed Abstract | CrossRef Full Text | Google Scholar

Bainbridge, M. L., Cersosimo, L. M., Wright, A. D. G., and Kraft, J. (2016). Rumen bacterial communities shift across a lactation in Holstein, Jersey and Holstein × Jersey dairy cows and correlate to rumen function, bacterial fatty acid composition and production parameters. FEMS Microbiol. Ecol. 92:fiw059. doi: 10.1093/femsec/fiw059

PubMed Abstract | CrossRef Full Text | Google Scholar

Bekele, A. Z., Koike, S., and Kobayashi, Y. (2010). Genetic diversity and diet specificity of ruminal Prevotella revealed by 16S rRNA gene-based analysis. FEMS Microbiol. Lett. 305, 49–57. doi: 10.1111/j.1574-6968.2010.01911.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Canibe, N., Højberg, O., Badsberg, J. H., and Jensen, B. B. (2007). Effect of feeding fermented liquid feed and fermented grain on gastrointestinal ecology and growth performance in piglets. J. Anim. Sci. 85, 2959–2971. doi: 10.2527/jas.2006-744

PubMed Abstract | CrossRef Full Text | Google Scholar

Cersosimo, L. M., Bainbridge, M. L., Kraft, J., and Wright, A.-D. G. (2016). Influence of periparturient and postpartum diets on rumen methanogen communities in three breeds of primiparous dairy cows. BMC Microbiol. 16:78. doi: 10.1186/s12866-016-0694-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Niu, L., and Zhang, Y. (2010). Saccharofermentans acetigenes gen. nov., sp. nov., an anaerobic bacterium isolated from sludge treating brewery wastewater. Int. J. Syst. Evol. Microbiol. 60, 2735–2738. doi: 10.1099/ijs.0.017590-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiquette, J., Allison, M. J., and Rasmussen, M. A. (2008). Prevotella bryantii 25A used as a probiotic in early-lactation dairy cows: effect on ruminal fermentation characteristics, milk production, and milk composition. J. Dairy Sci. 91, 3536–3543. doi: 10.3168/jds.2007-0849

PubMed Abstract | CrossRef Full Text | Google Scholar

Christopherson, M. R., Dawson, J. A., Stevenson, D. M., Cunningham, A. C., Bramhacharya, S., Weimer, P. J., et al. (2014). Unique aspects of fiber degradation by the ruminal ethanologen Ruminococcus albus 7 revealed by physiological and transcriptomic analysis. BMC Genomics 15:1066. doi: 10.1186/1471-2164-15-1066

PubMed Abstract | CrossRef Full Text | Google Scholar

Cotta, M. A. (1988). Amylolytic activity of selected species of ruminal bacteria. Appl. Environ. Microbiol. 54, 772–776.

Google Scholar

Dehority, B. A. (1966). Characterization of several bovine rumen bacteria isolated with a xylan medium. J. Bacteriol. 91, 1724–1729. doi: 10.1128/jb.91.5.1724-1729.1966

PubMed Abstract | CrossRef Full Text | Google Scholar

Dehority, B. A. (1969). Pectin-fermenting bacteria isolated from the bovine rumen. J. Bacteriol. 99, 189–196. doi: 10.1128/jb.99.1.189-196.1969

PubMed Abstract | CrossRef Full Text | Google Scholar

Delgado, B., Bach, A., Guasch, I., González, C., Elcoso, G., Pryce, J. E., et al. (2019). Whole rumen metagenome sequencing allows classifying and predicting feed efficiency and intake levels in cattle. Sci. Rep. 9:11. doi: 10.1038/s41598-018-36673-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Dieho, K., van den Bogert, B., Henderson, G., Bannink, A., Ramiro-Garcia, J., Smidt, H., et al. (2017). Changes in rumen microbiota composition and in situ degradation kinetics during the dry period and early lactation as affected by rate of increase of concentrate allowance. J. Dairy Sci. 100, 2695–2710. doi: 10.3168/jds.2016-11982

PubMed Abstract | CrossRef Full Text | Google Scholar

Difford, G. F., Plichta, D. R., Løvendahl, P., Lassen, J., Noel, S. J., Højberg, O., et al. (2018). Host genetics and the rumen microbiome jointly associate with methane emissions in dairy cows. PLoS Genet. 14:e1007580. doi: 10.1371/journal.pgen.1007580

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998. doi: 10.1038/nmeth.2604

PubMed Abstract | CrossRef Full Text | Google Scholar

Emmanuel, B., Milligan, L. P., and Turner, B. V. (1974). The metabolism of acetate by rumen microorganisms. Can. J. Microbiol. 20, 183–185. doi: 10.1139/m74-028

PubMed Abstract | CrossRef Full Text | Google Scholar

Flint, H. J., Bayer, E. A., Rincon, M. T., Lamed, R., and White, B. A. (2008). Polysaccharide utilization by gut bacteria: potential for new insights from genomic analysis. Nat. Rev. Microbiol. 6, 121–131. doi: 10.1038/nrmicro1817

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J., and Alm, E. J. (2012). Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8:e1002687. doi: 10.1371/journal.pcbi.1002687

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, N., Jami, E., and Mizrahi, I. (2017). Compositional and functional dynamics of the bovine rumen methanogenic community across different developmental stages. Environ. Microbiol. 19, 3365–3373. doi: 10.1111/1462-2920.13846

PubMed Abstract | CrossRef Full Text | Google Scholar

Geishauser, T., Linhart, N., Neidl, A., and Reimann, A. (2012). Factors associated with ruminal pH at herd level. J. Dairy Sci. 95, 4556–4567. doi: 10.3168/jds.2012-5380

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, G., Cox, F., Ganesh, S., Jonker, A., Young, W., Janssen, P. H., et al. (2015). Rumen microbial community composition varies with diet and host, but a core microbiome is found across a wide geographical range. Sci. Rep. 5:14567. doi: 10.1038/srep14567

PubMed Abstract | CrossRef Full Text | Google Scholar

Hildebrand, F., Tadeo, R., Voigt, A. Y., Bork, P., and Raes, J. (2014). LotuS: an efficient and user-friendly OTU processing pipeline. Microbiome 2:30. doi: 10.1186/2049-2618-2-30

PubMed Abstract | CrossRef Full Text | Google Scholar

Indugu, N., Vecchiarelli, B., Baker, L. D., Ferguson, J. D., Vanamala, J. K. P., and Pitta, D. W. (2017). Comparison of rumen bacterial communities in dairy herds of different production. BMC Microbiol. 17:190. doi: 10.1186/s12866-017-1098-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Jami, E., White, B. A., and Mizrahi, I. (2014). Potential role of the bovine rumen microbiome in modulating milk composition and feed efficiency. PLoS One 9:e85423. doi: 10.1371/journal.pone.0085423

PubMed Abstract | CrossRef Full Text | Google Scholar

Jewell, K. A., McCormick, C. A., Odt, C. L., Weimer, P. J., and Suen, G. (2015). Ruminal bacterial community composition in dairy cows is dynamic over the course of two lactations and correlates with feed efficiency. Appl. Environ. Microbiol. 81, 4697–4710. doi: 10.1128/AEM.00720-15

PubMed Abstract | CrossRef Full Text | Google Scholar

John Wallace, R., Sasson, G., Garnsworthy, P. C., Tapio, I., Gregson, E., Bani, P., et al. (2019). A heritable subset of the core rumen microbiome dictates dairy cow productivity and emissions. Sci. Adv. 5:eaav8391. doi: 10.1126/sciadv.aav8391

PubMed Abstract | CrossRef Full Text | Google Scholar

Kepler, C. R., Hirons, K. P., McNeill, J. J., and Tove, S. B. (1966). Intermediates and products of the biohydrogenation of linoleic acid by Butyrinvibrio fibrisolvens. J. Biol. Chem. 241, 1350–1354.

Google Scholar

Kim, J. N., Méndez–García, C., Geier, R. R., Iakiviak, M., Chang, J., Cann, I., et al. (2017). Metabolic networks for nitrogen utilization in Prevotella ruminicola 23. Sci. Rep. 7:7851. doi: 10.1038/s41598-017-08463-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, Y., Shinkai, T., and Koike, S. (2008). Ecological and physiological characterization shows that Fibrobacter succinogenes is important in rumen fiber digestion — Review. Folia Microbiol. (Praha) 53, 195–200. doi: 10.1007/s12223-008-0024-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopečný, J., Zorec, M., Mrázek, J., Kobayashi, Y., and Marinšek-Logar, R. (2003). Butyrivibrio hungatei sp. nov. and Pseudobutyrivibrio xylanivorans sp. nov., butyrate-producing bacteria from the rumen. Int. J. Syst. Evol. Microbiol. 53, 201–209. doi: 10.1099/ijs.0.02345-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Krause, D. O., Denman, S. E., Mackie, R. I., Morrison, M., Rae, A. L., Attwood, G. T., et al. (2003). Opportunities to improve fiber degradation in the rumen: microbiology, ecology, and genomics. FEMS Microbiol. Rev. 27, 663–693. doi: 10.1016/S0168-6445(03)00072-X

CrossRef Full Text | Google Scholar

Kumar, S., Indugu, N., Vecchiarelli, B., and Pitta, D. W. (2015). Associative patterns among anaerobic fungi, methanogenic archaea, and bacterial communities in response to changes in diet and age in the rumen of dairy cows. Front. Microbiol. 6:781. doi: 10.3389/fmicb.2015.00781

PubMed Abstract | CrossRef Full Text | Google Scholar

Lang, S. L. C., Iverson, S. J., and Bowen, W. D. (2012). Primiparous and multiparous females differ in mammary gland alveolar development: implications for milk production. J. Exp. Biol. 215, 2904–2911. doi: 10.1242/jeb.067058

PubMed Abstract | CrossRef Full Text | Google Scholar

Lassen, J., and Løvendahl, P. (2016). Heritability estimates for enteric methane emissions from Holstein cattle measured using noninvasive methods. J. Dairy Sci. 99, 1959–1967. doi: 10.3168/jds.2015-10012

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Li, C., Chen, Y., Liu, J., Zhang, C., Irving, B., et al. (2019). Host genetics influence the rumen microbiota and heritable rumen microbial features associate with feed efficiency in cattle. Microbiome 7:92.

Google Scholar

Lima, F. S., Oikonomou, G., Lima, S. F., Bicalho, M. L. S., Ganda, E. K., de Oliveira Filho, J. C., et al. (2015). Prepartum and postpartum rumen fluid microbiomes: characterization and correlation with production traits in dairy cows. Appl. Environ. Microbiol. 81, 1327–1337. doi: 10.1128/AEM.03138-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Madsen, J., Bjerg, B. S., Hvelplund, T., Weisbjerg, M. R., and Lund, P. (2010). Methane and carbon dioxide ratio in excreted air for quantification of the methane production from ruminants. Livest. Sci. 129, 223–227. doi: 10.1016/j.livsci.2010.01.001

CrossRef Full Text | Google Scholar

Maia, M. R., Chaudhary, L. C., Bestwick, C. S., Richardson, A. J., McKain, N., Larson, T. R., et al. (2010). Toxicity of unsaturated fatty acids to the biohydrogenating ruminal bacterium, Butyrivibrio fibrisolvens. BMC Microbiol. 10:52. doi: 10.1186/1471-2180-10-52

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Fernandez, G., Denman, S. E., Yang, C., Cheung, J., Mitsumori, M., and McSweeney, C. S. (2016). Methane inhibition alters the microbial community, hydrogen flow, and fermentation response in the rumen of cattle. Front. Microbiol. 7:1122. doi: 10.3389/fmicb.2016.01122

PubMed Abstract | CrossRef Full Text | Google Scholar

Nathani, N. M., Patel, A. K., Mootapally, C. S., Reddy, B., Shah, S. V., Lunagaria, P. M., et al. (2015). Effect of roughage on rumen microbiota composition in the efficient feed converter and sturdy Indian Jaffrabadi buffalo (Bubalus bubalis). BMC Genomics 16:1116. doi: 10.1186/s12864-015-2340-4

PubMed Abstract | CrossRef Full Text | Google Scholar

NCSS (2012). “Lin’s Concordance Correlation Coefficient,” in PASS Sample Size Software NCSS.Com, chap. 812, 1–9. doi: 10.2307/2532051

CrossRef Full Text | Google Scholar

Perez-Garcia, O., Lear, G., and Singhal, N. (2016). Metabolic network modeling of microbial interactions in natural and engineered environmental systems. Front. Microbiol. 7:673. doi: 10.3389/fmicb.2016.00673

PubMed Abstract | CrossRef Full Text | Google Scholar

Pitta, D. W., Indugu, N., Kumar, S., Vecchiarelli, B., Sinha, R., Baker, L. D., et al. (2016). Metagenomic assessment of the functional potential of the rumen microbiome in Holstein dairy cows. Anaerobe 38, 50–60. doi: 10.1016/j.anaerobe.2015.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Pitta, D. W., Kumar, S., Vecchiarelli, B., Baker, L. D., Ferguson, J. D., Thomsen, N., et al. (2014a). Temporal dynamics in the ruminal microbiome of dairy cows during the transition period. J. Anim. Sci. 92, 4014–4022. doi: 10.2527/jas.2014-7621

PubMed Abstract | CrossRef Full Text | Google Scholar

Pitta, D. W., Parmar, N., Patel, A. K., Indugu, N., Kumar, S., Prajapathi, K. B., et al. (2014b). Bacterial diversity dynamics associated with different diets and different primer pairs in the rumen of Kankrej cattle. PLoS One 9:e111710. doi: 10.1371/journal.pone.0111710

PubMed Abstract | CrossRef Full Text | Google Scholar

Poulsen, M., Schwab, C., Borg Jensen, B., Engberg, R. M., Spang, A., Canibe, N., et al. (2013). Methylotrophic methanogenic Thermoplasmata implicated in reduced methane emissions from bovine rumen. Nat. Commun. 4:1428. doi: 10.1038/ncomms2432

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, M. N., Dehal, P. S., and Arkin, A. P. (2009). Fasttree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol. Biol. Evol. 26, 1641–1650. doi: 10.1093/molbev/msp077

PubMed Abstract | CrossRef Full Text | Google Scholar

Rode, L. M., Genthner, B. R. S., and Bryant, M. P. (1981). Syntrophic association by cocultures of the methanol- and CO2-H2-utilizing species Eubacterium limosum and pectin-fermenting Lachnospira multiparus during growth in a pectin medium. Appl. Environ. Microbiol. 42, 20–22.

Google Scholar

Schloss, P. D., Westcott, S. L., Ryabin, T., Hall, J. R., Hartmann, M., Hollister, E. B., et al. (2009). Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541. doi: 10.1128/AEM.01541-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Seedorf, H., Kittelmann, S., and Janssen, P. H. (2015). Few highly abundant operational taxonomic units dominate within rumen methanogenic archaeal species in New Zealand sheep and cattle. Appl. Environ. Microbiol. 81, 986–995. doi: 10.1128/AEM.03018-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, W., Moon, C. D., Leahy, S. C., Kang, D., Froula, J., Kittelmann, S., et al. (2014). Methane yield phenotypes linked to differential gene expression in the sheep rumen microbiome. Genome Res. 24, 1517–1525. doi: 10.1101/gr.168245.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Gylswyk, N. O. (1995). Succiniclasticum ruminis gen. nov., sp. nov., a ruminal bacterium converting succinate to propionate as the sole energy-yielding mechanism. Int. J. Syst. Bacteriol. 45, 297–300. doi: 10.1099/00207713-45-2-297

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinet, L., and Zhedanov, A. (2010). A “missing” family of classical orthogonal polynomials. Appl. Environ. Microbiol. 81, 4697–4710. doi: 10.1088/1751-8113/44/8/085201

CrossRef Full Text | Google Scholar

Wallace, R. J., Rooke, J. A., McKain, N., Duthie, C.-A., Hyslop, J. J., Ross, D. W., et al. (2015). The rumen microbial metagenome associated with high methane production in cattle. BMC Genomics 16:839. doi: 10.1186/s12864-015-2032-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Garrity, G. M., Tiedje, J. M., and Cole, J. R. (2007). Naive Bayesian classifier for rapid assignment of rRNA\nsequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267. doi: 10.1128/AEM.00062-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Wathes, D. C., Cheng, Z., Bourne, N., Taylor, V. J., Coffey, M. P., and Brotherstone, S. (2007). Differences between primiparous and multiparous dairy cows in the inter-relationships between metabolic traits, milk yield and body condition score in the periparturient period. Domest. Anim. Endocrinol. 33, 203–225. doi: 10.1016/j.domaniend.2006.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Weimer, P. J. (2015). Redundancy, resilience, and host specificity of the ruminal microbiota: implications for engineering improved ruminal fermentations. Front. Microbiol. 6:296. doi: 10.3389/fmicb.2015.00296

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, M. Y., Sun, H. Z., Wu, X. H., Liu, J. X., and Guan, L. L. (2020). Multi-omics reveals that the rumen microbiome and its metabolome together with the host metabolome contribute to individualized dairy cow performance. Microbiome 8:64. doi: 10.1186/s40168-020-00819-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Shi, H., Wang, Y., Li, S., Cao, Z., Ji, S., et al. (2017). Effect of dietary forage to concentrate ratios on dynamic profile changes and interactions of ruminal microbiota and metabolites in holstein heifers. Front. Microbiol. 8:2206. doi: 10.3389/fmicb.2017.02206

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Z., Kristensen, L., Difford, G. F., Poulsen, M., Noel, S. J., Abu Al-Soud, W., et al. (2018). Changes in rumen bacterial and archaeal communities over the transition period in primiparous Holstein dairy cows. J. Dairy Sci. 101, 9847–9862. doi: 10.3168/jds.2017-14366

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Z., Noel, S. J., Difford, G. F., Al-Soud, W. A., Brejnrod, A., Sørensen, S. J., et al. (2017). Community structure of the metabolically active rumen bacterial and archaeal communities of dairy cows over the transition period. PLoS One 12:e0187858. doi: 10.1371/journal.pone.0187858

PubMed Abstract | CrossRef Full Text | Google Scholar