Influence of periparturient and postpartum diets on rumen methanogen communities in three breeds of primiparous dairy cows

Background Enteric methane from rumen methanogens is responsible for 25.9 % of total methane emissions in the United States. Rumen methanogens also contribute to decreased animal feed efficiency. For methane mitigation strategies to be successful, it is important to establish which factors influence the rumen methanogen community and rumen volatile fatty acids (VFA). In the present study, we used next-generation sequencing to determine if dairy breed and/or days in milk (DIM) (high-fiber periparturient versus high-starch postpartum diets) affect the rumen environment and methanogen community of primiparous Holstein, Jersey, and Holstein-Jersey crossbreeds. Results When the 16S rRNA gene sequences were processed and assigned to operational taxonomic units (OTU), a core methanogen community was identified, consisting of Methanobrevibacter (Mbr.) smithii, Mbr. thaueri, Mbr. ruminantium, and Mbr. millerae. The 16S rRNA gene sequence reads clustered at 3 DIM, but not by breed. At 3 DIM, the mean % abundance of Mbr. thaueri was lower in Jerseys (26.9 %) and higher in Holsteins (30.7 %) and Holstein-Jersey crossbreeds (30.3 %) (P < 0.001). The molar concentrations of total VFA were higher at 3 DIM than at 93, 183, and 273 DIM, whereas the molar proportions of propionate were increased at 3 and 93 DIM, relative to 183 and 273 DIM. Rumen methanogen densities, distributions of the Mbr. species, and VFA molar proportions did not differ by breed. Conclusions The data from the present study suggest that a core methanogen community is present among dairy breeds, through out a lactation. Furthermore, the methanogen communities were more influenced by DIM and the breed by DIM interactions than breed differences.


Background
In the United States, enteric methane emissions from ruminants are the second largest anthropogenic source of methane, contributing to 25.9 % of all methane emissions and global warming [1]. Methane production caused by rumen archaea (i.e., methanogens) leads to a 2-12 % net loss of the dairy cow's gross energy intake [2]. This loss contributes to a significant economic loss for farmers as it increases the quantity of feed needed to meet milk production demands.
The rumen is an anaerobic environment that houses a microbiome consisting of bacteria, protozoa, fungi, phages, and archaea. The bacteria, protozoa, and fungi (e.g., yeast) ferment feedstuff consumed by the host, and produce VFA. Acetate, butyrate, and propionate, the predominant VFA in the rumen, are the main energy sources for the host animal. Fermentation byproducts such as carbon dioxide, formate, hydrogen gas, methanol, and methylamines are used by methanogenic archaea for methane production. The majority of the methane is eructated and exhaled out by the ruminant into the environment. For methane mitigation strategies to be successful, it is important to identify factors that may influence the rumen environment and thus, affect the methanogen density and diversity.
Holstein and Jersey dairy cattle are the two most common dairy breeds used in the United States. Holstein cows are recognized for their high milk production, whereas Jersey cows are recognized for their increased fertility and higher milk components. Additionally, there is global interest in Holstein-Jersey crossbreeds to compensate for the decreased fertility in Holsteins and milk production in Jersey cows. It has been demonstrated that first generation Holstein-Jersey crosses have dry matter intakes, milk yields and solids in between those measured in Holstein and Jersey cows, respectively [6].
The transition period from a diet high in neutral detergent fiber (NDF) to a diet high in starch is a challenge for lactating dairy cattle. Prior to parturition, the NDF content in the diet is elevated, while after parturition the energy content is increased with higher starch and fat levels. Kumar et al. [7] showed no difference in archaeal Shannon diversity or taxa when cows were transitioned from a high-fiber pre-partum diet to a low-fiber postpartum diet, however, no studies described the rumen methanogen community across a lactation period. When quantifying VFA, another study observed that concentrations of total VFA, acetate, and propionate were decreased during the transition period in comparison to 100 days in milk (DIM) [8].
Previous research focused on rumen bacteria in preand post-partum dairy cattle, but rumen methanogens have not been identified or quantified under these conditions. Furthermore, the rumen methanogens of Holstein and Jersey dairy cattle with different parities and DIM were identified with limited data generated from pooled PCR samples using clone libraries. The present study focused on Holstein-Jersey crossbreeds, used nextgeneration sequencing (NGS), and animals of the same age, DIM, and parity. Given previous investigations into the rumen methanogen community in relation to breed and what is known about transitioning dairy cattle from one diet to another, we hypothesized that the rumen methanogen diversity and rumen VFA proportions in primiparous dairy cattle are affected by both breed and DIM, while methanogen densities do not vary. The objectives of the present study were to (1) measure the rumen VFA, (2) use NGS techniques to identify rumen methanogens, (3) distribute the archaeal 16S rRNA gene sequence reads into operational taxonomic units (OTU), (4) quantify the rumen methanogens, and (5) correlate VFA with specific rumen methanogen taxa from each breed during early (3 DIM), peak (93 DIM), mid-(183 DIM), and latelactation (273 DIM).

Results
The 16S rRNA gene sequence data set is accessible through NCBI's Sequence Read Archive, under the study accession number [SRP058775].

Rumen volatile fatty acids
Breed and breed x DIM differences in total VFA concentrations or in individual VFA molar proportions were not observed. Total VFA concentrations were highest at 3 DIM (P < 0.01). Propionate proportions were lowest at 273 DIM and highest at 3 and 93 DIM (P < 0.05). Relative to 3 DIM, acetate proportions were higher at 183 and 273 DIM (P <0.001). Isobutyrate and lactate proportions were highest at 273 DIM (P < 0.01). Isovalerate proportions did not differ by DIM (Table 1).

OTU-based analyses
The 16S rRNA gene sequence reads clustered into 403 (3 DIM), 383 (93 DIM), 590 (183 DIM), and 547 OTUs (273 DIM). Breed and breed by DIM did not affect rumen methanogen diversity measures (Table 3). Good's coverage, Shannon diversity and Inverse Simpson indices were affected by DIM. The Inverse Simpson indices were highest at 3 and 183 DIM (P < 0.01), while Good's coverage and the Shannon Diversity indices were highest at 3 DIM (P < 0.05). The most and least OTUs shared between all animals were at 93 and 273 DIM, respectively (

Discussion
The present study is the first to investigate the rumen methanogen community across a lactation period in three dairy cattle breeds. The purpose of this experiment was to provide more knowledge about the rumen methanogen community and rumen parameters at 3, 93, 183, and 273 DIM in Holstein, Jersey, and Holstein-Jersey crossbreeds. This study identified the core methanogen    community with NGS technologies, quantified rumen VFA and methanogen densities, and correlated rumen methanogen species to one another and to VFA.

Volatile fatty acids
VFA are the main energy source provided to lactating dairy cattle and are the by products of carbohydrate fermentation by rumen bacteria, protozoa, and fungi. Generally, propionate is a precursor to glucose and is increased when animals are provided a high-starch diet or provided the ionophore, monensin. Relative to 183 and 273 DIM, proportions of propionate were increased at 3 and 93 DIM, suggesting a greater demand for glucose by the cow during early lactation. Although the animals were provided 0.06 % monensin pre-partum versus 0.02 % post-partum, it is not possible to correlate the increase in propionate at 3 and 93 DIM with this additive. An effect from monensin would be more plausible at 3 DIM, when the cows were transitioning from a prepartum to a post-partum diet, but this would not explain why propionate was also increased at 93 DIM. Furthermore, the increase in total VFA concentrations observed at 3 DIM suggests an increase in carbohydrate fermentation at the start of lactation. Johnson et al. [9] stated that the fermentation of fiber is favored, providing insight into why VFA concentrations were elevated at 3 DIM versus 93, 183, and 273 DIM. In contrast, Danielsson et al. [10] found that total VFA concentrations did not vary in cannulated mid-lactation dairy cattle consuming 500:500 and 900:100 g/kg dry matter forage to concentration diets.
The present study is the first to compare the methanogen densities in three breeds of dairy cattle and by DIM. In agreement with our hypothesis, the methanogen densities did not vary by breed or DIM. Previously reported methanogen densities (i.e., log 10 mcrA gene copies) from Between H-J c 12 13 11 9 Between H-X 11 13 12 8 Between J-X bulls on high-fiber (9.02) and starch diets (9.07) [11] were higher than what was observed in our study [11]. However, differences between the methanogen densities were not observed between the two diet groups [11]. In a study by Zhou et al. [12], the use of an exogenous fibrolytic feed enzyme additive did not affect the methanogen densities, yet, affected the methanogen community and methane production of lactating Holstein cows. Therefore, it appears that methanogen densities are not markedly affected by these specific diet alterations.
Previous work in dairy [4,7,10] and beef cattle [11,13] also showed the genus Mbr to be the most predominant genus. As methanogens belonging to the genus Mbr use the rumen fermentation byproducts, such as hydrogen and carbon dioxide as substrates for methanogenesis, it is thought that the high levels of these byproducts in the rumen enable these methanogens to thrive over other species that rely on scarce substrates such as methylamines, methanol, or acetate [14,15].
Although Mbr is the most abundant archaeal genus in ruminants, there are several species that are distributed into two different phylogenetic clades (i.e., SGMT and RO). In the present study, the SGMT clade was the most dominant branch by breed and DIM. Both Mbr. smithii and Mbr. thaueri made up the majority of the SGMT clade, while Mbr. ruminantium made up the majority of the RO clade. Previous research, using the same archaeal forward primer (Met86F), suggested a difference between SGMT-RO clade distributions between Holstein and Jersey cows [4]. However, the study revealed several limitations in the interpretation of the results. Animals were not blocked by parity, DIM, or age, the PCR products were pooled by breed and a clone library was constructed for each breed, and a limited number of clones were sequenced. It is conceivable that these variables, but not primer bias, may have contributed to the observed breed differences. Another study showed a prevalence of the RO clade in both corn-fed Hereford crossbreed and potatofed Hereford feedlot cattle in Canada [13]. Finally, the present study showed a strong negative correlation between the two clades suggesting that ruminants possess either a high abundance of SGMT or of RO and that dairy breed and DIM do not impact these proportions.
Because the four methanogen species Mbr. smithii, Mbr. thaueri, Mbr. ruminantium, and Mbr. millerae were identified in each breed and at each DIM time point investigated, our data showed the presence of a core methanogen community. Jeyanathan et al. [16] identified a common methanogen community between Holstein-Jersey crossbreeds, sheep, and red deer. Finding a core rumen methanogen community will enable further investigations into targeting specific species that are key contributors to methane production. Future work could isolate these species and determine which species produce the most methane.
Three out of the four methanogen species in the present study were previously identified in both Holstein and Jersey cows, while Mbr. thaueri was not. Recently, Mbr. thaueri was identified with the same primer pair used in the present study, at a high abundance in wild impalas from South Africa [17]. Omission of Mbr. thaueri from previous studies could be due to lack of sequencing depth or diet of the animal. The higher abundance of Mbr. thaueri in Jersey cows at 93 DIM (i.e., peak lactation) may be a result of a higher dry matter intake (DMI), but future studies are needed to draw a clear link between DMI, milk yield, and the rumen methanogen species identified. It is conceivable that Mbr. thaueri and the other three methanogen species persisted because the rumen environment and the substrates created by bacteria, protozoa, and fungi enabled these methanogens to thrive.
Throughout the lactation period and by breed, the methanol-utilizing genus, Msp. was identified in low abundances. Its mean % abundance was highest around peak lactation (93 DIM) and lowest at 3, 183, and 273 DIM. Similarly, Kumar et al. [7] compared the methanogen diversity between Holsteins at four weeks before calving and 1-5 days after calving and observed no differences in the genera Mbr or Msp. At 1-5 DIM, the methanol-utilizing genus Msp was more abundant (4.5 %) in primiparous Holsteins than those from the present study (<1 %). Like in the present study, animals were stomach tubed 2-3 h post-feeding and received a diet before calving with the same NDF content (44 %, [7]). Msp is typically more abundant in animals consuming feeds with elevated pectin levels [14]. Although not analyzed, it is possible that the diet in the present study had a lower quantity of pectin.
Although both OTU and 16S rRNA gene sequence classifications identified a core methanogen community, certain methanogen species were more abundant at different DIM time points. Relative to 93, 183, and 273 DIM the proportions of Mbr. millerae and Mbr. woesei were highest at 3 DIM. While previous research has not focused on rumen methanogen communities pre-and post-partum, one study suggested that the highly cellulolytic anaerobic fungi were more prevalent in pre-partum dairy cows, while rumen protozoa were less prevalent [18]. This suggests that as the dairy cows transition to a diet typically fed post-partum, there is a shift in the rumen microbiome and likely in the methanogen community.
The OTU coverage in the present study was almost 100 %, indicating a sufficient sampling effort. OTU distribution did not vary by breed or DIM and the majority of the sequences were distributed into four main OTUs. The methanogen diversity in the present study was not influenced by breed, but by DIM. King et al. [4] reported a higher Shannon diversity index and number of OTUs in lactating Holstein cows than in Jersey cows. However, the limited number of sequences from the cloned libraries, pooled samples by breed, and parity may have influenced the results. According to Kumar et al. [7], multiparous Holstein cows exhibit a higher Shannon diversity index than primiparous.
The Shannon diversity indices (i.e., species evenness and abundance) from the present study were highest at 3 DIM, while the 16S rRNA gene sequences reads clustered during this time as well. At 93, 183, and 273 DIM the sequence reads were mixed in one cluster. Previous research also demonstrated that the Shannon diversity indices of methanogens increased in bulls on a high-fiber diet (0.95) and decreased with a high-starch diet (0.79) [11]. These data suggested that a high-fiber diet leads to a more diverse methanogen community when compared to a high-starch diet. In another study, the Shannon diversity of Holsteins at 4 weeks before calving and 1-5 days after calving did not differ, but was most likely because there was not enough time between sampling [7]. Belanche et al. [19] suggested that the increased amount of cellulose and other heteropolysaccharides in a diet high in fiber leads to a more diverse microbial community. Therefore, a more diverse bacterial, fungal, or protozoal community may provide different substrates that enable the presence of a more diverse methanogen community.

Conclusions
The data presented here are the first to characterize the rumen methanogen communities in three dairy cattle breeds across a lactation period. NGS produced over 1 million sequence reads and demonstrated that diversity was different at 3 DIM. Notably, a core methanogen community persisted and consisted of four species, Mbr. smithii, Mbr. thaueri, Mbr. ruminantium, and Mbr. millerae. These methanogens may play a significant role in methanogenesis and in the utilization of substrates from bacterial, protozoal, and fungal fermentation. However future work is required to better delineate these relationships. The SGMT-RO clades did not vary by breed or DIM, instead, the SGMT clade was dominant in all three breeds. Although our results show that breed does not affect the rumen methanogen taxa per se, more studies are needed to clarify if this finding is consistent in other geographic locations and in dairy cattle consuming varying diets. was excluded from all data analyses because of post-partum health concerns. The Institutional Animal Care and Use Committee at the UVM approved all animal sample collection methods under protocol # 13-031. All animals calved within a 2-month period. At 3, 93, 183, and 273 DIM, whole rumen digesta samples (50 mL) were collected 2-3 h post-feeding at 0900 h via stomach intubation from each animal. To collect rumen samples, a flexible milk hose (2.54 cm diameter) was passed through a speculum to the esophagus and to the rumen. The hose was marked at 200 cm to indicate the approximate location of the rumen. Once the tube contacted the fiber mat and rumen sounds were heard, a 600 cc livestock drench gun (Labelvage, France) collected the digesta. Whole digesta samples were immediately frozen at −20°C to minimize microbial activity.

Diet
Prior to calving, all animals consumed a pre-partum total mixed ration (TMR) diet. Within 24 h post-partum, each cow was transitioned to a diet that was higher in starch and lower in NDF in comparison to the diet fed pre-partum (Table 5). Through out the study, a 70:30 forage to concentrate TMR was fed. Prior to calving the forages included corn silage (51.   [20] demonstrated that monensin does not alter the rumen methanogen diversity or density in lactating dairy cows. TMR samples were collected weekly for three consecutive days and composited at the end of each week. Because cows calved within two months of one another, the mean and the standard error of the diet provided before calving are reported. Cumberland Valley Analytical Services (Hagerstown, MD) analyzed the feed samples and provided nutrient composition data.

Volatile fatty acid analysis
The whole rumen digesta samples were spun in a Beckman J2-21 centrifuge at 10,000 x g for 20 min at 4°C and the resulting supernatant was filtered through a 25 mm hardened ashless filter (Whatman Inc. Clifton, NJ) to remove debris. Subsequently, the samples were diluted 1:1 with 0.06 M oxalic acid containing 50 μM trimethyl acetic acid (internal standard) and analyzed for VFA by gas chromatography (Varian 3800 GC, Walnut Creek, CA) coupled with a flame ionization detector and a customized packed column (2 m x 2 mm ID glass) with 4 % carbowax and 80/ 120 Carbopac B-DA (Sulpeco, Bellefonte, PA) with nitrogen as carrier gas (15 mL/min flow rate). The other gases were purified air at 300 mL/min and hydrogen makeup gas at 30 mL/min. The column was operated at 175°C; the total run time was 25 min. Both the injector and detector temperature were kept at 200°C. The injection volume was 1 μL. The identification of VFA was based on retention times against known standards using software Star Chromatography v5 (Varian) and quantified for their concentrations using respective VFA standards. Results are expressed in mM VFA.

Microbial DNA extraction
Across 4 time points (3,93,183, and 273 DIM), 87 individual whole rumen digesta samples were collected. The previously frozen samples were thawed overnight at 4°C. Each sample was vortexed for 30s to homogenize the sample and break up the solid particles that settled to the bottom of the conical tube. The microbial DNA was extracted using the repeated-bead beating plus column (RBB + C) method [21] and followed previously described procedures [17].

Real-time PCR amplification
The log 10 of the copy number of the mcrA gene per mL of rumen digesta was determined by real-time PCR, each sample was amplified in triplicate. Each real-time PCR included: 12.5 μL of SYBR® Green Mix, 6.5 μL of double distilled water, 2.5 μL of the methanogen-specific primer pair, mcrA-F and mcrA-R [22], and 1 μL of either diluted template DNA (10 ng/μL), positive (mixture of microbial DNA extract from King et al. [4]) or negative controls (doubledistilled water). The mcrA gene was amplified in a Bio-Rad C1000 Touch Thermal Cycler (Bio-Rad, Hercules, CA) under previously published conditions [22]. The five mcrA gene standards were created by serial dilutions of 10 ng/μL of purified PCR product. The range of concentrations for the standard curve was from 0.001-10 ng/μL. An acceptable standard curve had an R 2 value greater or equal to 0.997. Using the BioRad CFX Manager (v.3.0) a model equation of the standard curve (y = mx + b) was used, where "x" was the log of the starting quantity, "b" the yintercept of the quantification value (Cq), and "m" the slope of the line for log starting quantity versus Cq values. Individual densities were calculated using previously established methods [23].

PCR Amplification of the 16S rRNA gene
The archaeal-specific primer pair, Met86F [24] and Met471R [17] were used to amplify the V1-V3 hypervariable regions of the 16S rRNA gene via PCR, following previously published procedures [17]. Purified archaeal amplicons (25 μL) were sent to Molecular Research DNA Laboratories (Shallowater, TX) and sequenced with the Illumina MiSeq version 3 NGS platform.

Bioinformatics workflow used to analyze MiSeq sequences
The program MOTHUR, version 1.33.3, was used to perform bioinformatics analyses in-house [25]. Each time period was analyzed separately because of the lack of computing power and memory. Prior to using MOTHUR, a Perl script was used to trim the sequences to 350 bp at the reverse primer and each sequence was quality checked. All sequences with a Phred quality score of 25 or above were kept for further analyses. The command, trim.seqs removed barcodes and created a file that identified which sample belonged to which sequences.
To determine the number of unique sequences in the data set, the command unique.seqs was used. A Needleman-Wunsch pairwise alignment and an aligned reference file of known rumen methanogen 16S rRNA gene sequences were used with the command align.seqs to align the unique sequences. Chimeric sequences were identified with UChime [26] and removed with MOTHUR.
In order to detect any bacterial 16S rRNA gene sequences, the sequences were classified with the 16S rRNA reference Ribosomal Database Project (RDP) files provided by MOTHUR. These files contained known bacteria and methanogen sequences from kingdom to genus taxonomic levels. The files were modified to contain species-level names. Any bacteria sequences (<0.01 %) were removed. The online RDP Classifier was used at a 95 % confidence threshold to further quality check the sequences. Sequences with an unknown root were removed. Taxonomy and fasta files containing 765 known archaeal species names and 16S rRNA gene sequences were used to classify the sequence reads into taxonomic species. The command cluster.split, was used with a 2 % cutoff to cluster the 16S rNA gene sequences into OTUs. Once OTUs were formed, they were classified with the command, classify.otu. The command summary.single calculated OTU-based alpha diversities, Shannon Diversity Index, Chao I richness estimator, Inverse Simpson index, and Good's coverage.
The subsample parameter in MOTHUR was used to analyze the same number of sequences per individual sample. Shared OTUs within and between breeds were counted with get.sharedseqs. To correlate OTUs (e.g., OTU 1 to OTU 2 at 3 DIM), the otu.association command used a Pearson correlation. The command, merge.file combined the fasta files of unique sequences. A subsample of 100,000 sequences was taken and the distances between sequences with a 2 % cutoff were calculated with dist.seqs. Once a phylip distance file was created, the clearcut command created a phylogenetic tree file. The output from clearcut was used with Fast Unifrac to perform a PCoA [27]. The PCoA determined if sequences from each sample cluster were based on breed or DIM.

Statistical analyses
All data were analyzed with a repeated measures ANOVA model in SAS 9.4 (SAS Inst. Inc., Cary, NC) with PROC MIXED. The model included breed, DIM, and breed by DIM interactions as fixed effects and used the Kenward-Roger method to determine the degrees of freedom. A Pearson correlation, to determine the relationship between VFA and methanogen taxa, was performed with the PROC CORR. The online data visualization tool, Plotly, was used to generate a heatmap of the correlation values (r). Statistical significance was declared at P < 0.05 and trends were declared at 0.05 ≤ P ≤ 0.10.

Ethics approval and consent to participate
The present study was performed in accordance with the Institutional Animal Care and Use Committee at the University of Vermont under protocol # 13-031.