Rumen Microbiome from Steers Differing in Feed Efficiency

The cattle rumen has a diverse microbial ecosystem that is essential for the host to digest plant material. Extremes in body weight (BW) gain in mice and humans have been associated with different intestinal microbial populations. The objective of this study was to characterize the microbiome of the cattle rumen among steers differing in feed efficiency. Two contemporary groups of steers (n=148 and n=197) were fed a ration (dry matter basis) of 57.35% dry-rolled corn, 30% wet distillers grain with solubles, 8% alfalfa hay, 4.25% supplement, and 0.4% urea for 63 days. Individual feed intake (FI) and BW gain were determined. Within contemporary group, the four steers within each Cartesian quadrant were sampled (n=16/group) from the bivariate distribution of average daily BW gain and average daily FI. Bacterial 16S rRNA gene amplicons were sequenced from the harvested bovine rumen fluid samples using next-generation sequencing technology. No significant changes in diversity or richness were indicated, and UniFrac principal coordinate analysis did not show any separation of microbial communities within the rumen. However, the abundances of relative microbial populations and operational taxonomic units did reveal significant differences with reference to feed efficiency groups. Bacteroidetes and Firmicutes were the dominant phyla in all ruminal groups, with significant population shifts in relevant ruminal taxa, including phyla Firmicutes and Lentisphaerae, as well as genera Succiniclasticum, Lactobacillus, Ruminococcus, and Prevotella. This study suggests the involvement of the rumen microbiome as a component influencing the efficiency of weight gain at the 16S level, which can be utilized to better understand variations in microbial ecology as well as host factors that will improve feed efficiency.


Introduction
The bovine gastrointestinal (GI) tract is a complex system that is responsible for animal nutrient uptake and overall health. Rumen function and the microbiology of culturable rumen microorganisms has been well-studied, but it has only been since the advent of next-generation sequencing technology in the last decade that research on total microbial diversity, microbial community function, and their effects on the host could be quantitatively examined at a higher resolution. These new technologies have enabled a thorough examination of the host microbiota, independent of culture-based methods.
The gastrointestinal tract is host to a diverse microbial ecosystem that can vary depending on both host genetic and environmental factors [1]. Studies have shown that even minor shifts in these populations can have a tremendous impact on livestock nutrition and productivity [2][3][4]. Changes in the composition and diversity of the ruminal microbiota have been linked with diet and age [5]. Extremes in body weight (BW) gain in mice and humans have been associated with different intestinal microbial populations [6][7]. However, the majority of research in cattle has focused on the microbial responses to changes in external parameters, such as diet influences, and management practices [8][9]. Little has been done to examine the effect of differing microbial populations on host phenotypes, specifically related to performance and feed efficiency. Furthermore, many of these studies lack adequate sample size; examining only a few subjects. With advances in sequencing technology, researchers can examine these relationships with greater depth and coverage than has previously been reported.
Feed costs remain the largest variable cost in beef production [10]. Optimizing feed efficiency in cattle has long been an effort devoted to host genetics, management, and diet. Yet, data supports possible microbial interactions within the host, influencing a multitude of processes pertaining to digestion and the host GI tract. We hypothesize that variable microbial populations have differential abilities to degrade feedstuffs in the rumen into nutrients available for absorption, as a possible route to impacting feed efficiency. This study aimed to characterize the microbiome, specifically the bacterial community, of the cattle rumen among steers differing in feed intake and growth, in order to assess the association of the ruminal microbial community profile with variation in bovine feed efficiency.

Experimental design and rumen sampling
Steers selected for this study came from a populations of cattle being developed to have a high percentage of the following breeds: Angus, Beefmaster, Brahman, Brangus, Braunvieh, Charolais, Chiangus, Gelbvieh, Hereford, Limousin, Maine Anjou, Red Angus, Salers, Santa Gertrudis, Shorthorn, Simmental, South Devon, and Tarentaise. Each year heifers and cows were artificially inseminated with semen from prominent industry bulls of their dominant breed. This program resulted in offspring ranging from 50% to 75% of the same breed as their sire with the exception of Angus and Hereford which ranged from 50% to 100% of the same breed as their sire. Individual feed intake was measured using an Insentec feeding system (Marknesse, The Netherlands). Steers were fed a ration (dry matter basis) of 57.35% dry-rolled corn, 30% wet distillers grain with solubles, 8% alfalfa hay, 4.25% supplement (containing 0.772 g/kg monensin), and 0.4% urea. Feed intake and BW gain were measured over a 63 day period [11]. Steers were selected from two contemporary groups. Group 1 (n = 148) were spring-born calves that were 371 ± 1 d of age and weighed 522 ± 4 kg at the start of the feed intake measurement. Group 2 (n = 197) were fall-born calves that were 343 ± 1 d of age and weighed 448 ± 4 kg at the start of the feed intake measurement. At the end of each feeding period, average daily BW gain (ADG) was regressed on average daily dry matter intake (ADFI). The four steers within each Cartesian quadrant were sampled (n = 16/group) from the bivariate distribution of ADG and ADFI. The result was a 2 X 2 factorial design consisting of high and low ADFI, and high and low ADG (Fig 1). At the end of the feeding period, steers were harvested, and rumen fluid was sampled and strained through 4 layers of cheesecloth. The 2 feeding studies yielded 32 animals for analysis. Due to the high-concentrate diet and whole-sample straining, rumen fluid samples were used for this study. Samples were individually stored in buffered peptone water (BPW, pH 7.0) + 15% glycerol stock for processing and kept at -70°C for long-term storage post-processing.
DNA extraction, amplification and sequencing DNA was extracted from rumen samples using a repeated bead beating plus column (RBB+C) method [12]. Briefly, 0.3g of sample was centrifuged for 5 min at 16,000× g to pellet solids including bacterial cells, then resuspended in 0.2 mL TE (Tris-EDTA, pH 8.0) buffer. Cell lysis was achieved by bead beating 0.15 g of the resuspended sample in ZR BashingBead Lysis Tubes (Zymo Research Corp, Santa Ana, CA, USA) using the TissueLyser II system (Qiagen, Hilden, Germany) for 3 min at 21Hz, in the presence of 4% (w/v) sodium dodecyl sulfate (SDS), 500 mM NaCl, and 50 mM EDTA. After mechanical and chemical lysis, ammonium acetate was used to precipitate and remove the impurities and the SDS, along with isopropanol precipitation for the recovery of the nucleic acids. RNA and proteins were removed or degraded using RNase and proteinase K, followed by the use of QIAamp columns from the Qiagen DNA Stool Mini Kit (Qiagen, Hilden, Germany). Genomic DNA concentration was determined using a Nanodrop 1000 spectrophotometer (ThermoScientific, Wilmington, DE, USA).
Amplicon library preparation was performed by PCR amplification of the V1-V3 region of 16S rRNA gene, using modified universal primers 27F (5'-Adapter / Index / AGAGT TTGATCCTGGCTCAG) and 519R (5' Adapter / Index / GTATTACCGCGGCTGCTG) including TruSeq adapter sequences and indices, as well as AccuPrime Taq high fidelity DNA Polymerase (Life Technologies, Carlsbad, CA). Amplification consisted of 20 cycles, with an annealing temperature of 58°C. Products were purified using AmPure bead purification (Agencourt, Beverly, MA) and all libraries were quantified by the PicoGreen dsDNA quantitation kit (Invitrogen, Carlsbad, CA) and by real-time PCR on the LightCycler 480 system (Roche, Mannheim, Germany). The PCR amplicon libraries were sequenced using the 2x300, v3 600-cycle kit and the Illumina MiSeq sequencing platform (Illumina, San Diego, CA). Feed Efficiency Sampling Model. A 2x2 factorial design consisting of high and low average daily feed intake (ADFI), and high and low average daily BW gain (ADG). Within a contemporary group, ADG was regressed on ADFI. Steers within each Cartesian quadrant were sampled, as represented by red dots (n = 16/group). Contemporary group 1(n = 197) ADFI was 10.33 ± 0.09 kg/d, and ADG was 1.54 ± 0.02 kg/d. Contemporary group 2 (n = 148) ADFI was 10.26 ± 0.12 kg/d, ADG was 1.90 ± 0.02 kg/d.

Sequence read processing and analysis
All sequences were processed using the QIIME-1.8.0 software package. Paired reads were joined using fastq-join [13] and filtered for quality (Q25) using the Galaxy server [14]. Sequences that contained read lengths shorter than 400 bp were removed and adapters/index sequences were trimmed. Chimeric sequences were checked using ChimeraSlayer [15]. All cleaned sequences were classified into taxa using the Greengenes 16S rRNA Gene Database [16]. Operational taxonomic units (OTUs) were calculated using the uclust program at 0.03 dissimilarity [17]. After calculating richness for each quadrant, singletons were removed from further diversity analyses. Based on rarefaction curves, the number of OTUs was normalized via random subsampling of 25,000 sequences from each rumen sample. A phylogenic tree was built with FastTree [18] to determine alpha and beta diversity metrics.

Statistical analysis
The mean abundances (n = 8) of data metrics and each taxon were compared among the feed efficiency groups using a model of contemporary group and Cartesian quadrant ((high ADG, high ADFI (ADG High -ADFI High ); high ADG, low ADFI (ADG High -ADFI Low ); low ADG, low ADFI (ADG Low -ADFI Low ); low ADG, high ADFI (ADG Low -ADFI High )) as fixed effects. The Benjamini-Hochberg method was used for multiple-testing correction [19]. Multiple-testing corrections were made for the number of phyla, the number of OTU groups, and other classified taxa groups. Significant difference was determined at P < 0.05. Linear contrasts were then applied to quadrant to separate whether microbial populations varied by low vs. high ADG, low vs. high ADFI, or their interaction (P < 0.01). Principal coordinates analysis (PCoA) was performed using weighted and unweighted UniFrac analyses [20].

Diversity of Ruminal Bacterial Communities
After quality control, chimera detection, and removal, the sampled ruminal contents of 32 steers, grouped into 4 feed efficiency phenotypes yielded a total of 5,565,909 cleaned reads with an average read length of 500 bp. The total cleaned reads represented individual samples ranging from 25,094 to 888,867 reads. OTUs were defined as a read sharing 97% nucleotide sequence identity. From the total sequences, 68,655 OTUs were detected with an average of 17,164±3619 OTUs per quadrant. From each Cartesian quadrant, data ranged from 13,304-21,971 OTUs, with the greatest number of OTUs stemming from the ADG High -ADFI High group and the fewest from the ADG Low -ADFI High group. The quadrants representing ADG High -ADFI Low and ADG Low -ADFI Low yielded 17,337 and 16,043 OTUs, respectively. The number of singletons accounted for approximately 38% of all OTUs. Coverage ranged from 98.31 to 99.56% as determined by Good's coverage estimator. Richness (Shannon) was slightly reduced, but not significant in the ADG High -ADFI High group at 6.32 compared to the other 3 groups (6.89-7.11).
In order to normalize individual samples for comparison among phenotype classes, each sample OTU table was rarefied to 25,000 reads, which was based upon optimizing both minimum sample read output and sample rarefaction curves. Each normalized sample was then analyzed and compared using the sample means within each quadrant. The normalized sequence reads were examined via bacterial diversity (Shannon), richness (Chao-1), and coverage (Good's coverage estimator). The average number of OTUs from the normalized samples did not differ (P > 0.05) among the 4 feed efficiency groups, which had an average of 1098±382 OTUs per group ( Table 1). The data representing richness similarly did not indicate any significance (P > 0.05) among the number of OTUs, which ranged from 1298-1601. Though some variability was observed, the diversity indices for each group were not significantly different. Each group was adequately covered, up to 98.84% in the ADG Low -ADFI High group and as low as 98.37% in the ADG High -ADFI High group.
The data was then reduced to an OTU-centric method, principal coordinates analysis (PCoA), to determine any separation into sample clusters. This was achieved by applying the phylogeny-based method, UniFrac, to the data. This method is a β-diversity measure that takes the phylogenic divergence between the OTUs into account aiding in identification of differences among microbial communities [21]. In both the weighted (quantitative) and unweighted (qualitative) UniFrac distances [22], there was no separation into clusters observed in the PCoA (Fig 2).
Analysis of the relative taxonomic abundance across all feed efficiency groups revealed differences among phyla and genera ( Table 2 and Fig 3). The data were analyzed to reflect the mean of the relative abundance (reads of a taxon/total reads in a sample) in each feed efficiency group. The phyla Firmicutes (P = 0.0364) and Lentisphaerae (P = 0.0339), and genera Dialister (P = 0.0062), Lactobacillus (P = 0.0419), Acidaminococcus (P = 0.0306), Anaerovibrio (P = 0.0291) Lysobacter (P = 0.0462), Janibacter (P = 0.0161), and Leucobacter (P = 0.0215) demonstrated significant differences in abundance between the feed efficiency groups. Additionally, several lower orders of classification were significant between the groups, which included the families Lachnospiraceae (P = 0.0197), Veillonellaceae (P = 0.0058), and Helicobacteraceae (P = 0.0160). Examination of the taxonomic profiles for the relative phylum-level and genuslevel abundances demonstrated changes between the efficiency groups. This was most evident within the Bacteroidetes, Firmicutes, and Prevotella identities (Fig 3). However, the high variability permitted limited determination of any visual trends within the profiles, specifically at lower abundances.

Effect of Gain and Intake
The association of the microbial populations with ADG and ADFI were analyzed similarly to the aforementioned feed efficiency data in order to determine whether microbial populations differed by low vs. high ADG, low vs. high ADFI, or their interaction. This was conducted to determine the association of microbial population factors contributing to feed efficiency. Tables 4 and 5 list the significant relative abundances of taxa and OTUs between ADG and ADFI groups (P<0.01). As expected, many of the same taxa and OTUs listed in the feed efficiency data were present, with the majority of the populations associating with ADG. Gain-associated taxa included the phylum Lentisphaerae (P = 0.00688), families Veillonellaceae (P = 0.00184) and Lachnospiraceae (P = 0.00214), and the genera Dialister (P = 0.00005) and Acidaminococcus (P = 0.00460). Although many gain-associated OTUs were identified, only three OTUs were intake-associated, which included order Clostridiales (P = 0.0097) and families Coriobacteriaceae (P = 0.0006) and Veillonellaceae (P = 0.0024).

Discussion
In this study, the V1-V3 rrs region was chosen to analyze the microbial diversity within the rumen of steers. A myriad of studies report microbial ecology data utilizing 16S variable regions, and due to read lengths available for the most common next generation sequencing platforms, a subset of variable segments are targeted, rather than the entire length of the gene. To this extent, data regarding each specific ecological community can be confounded and biased based upon the variable region(s) selected [24][25][26]. To that end, careful consideration was taken to best represent the microbial communities in the rumen and lower GI. The region spanning the V1-V3 regions was selected for the present study, based on analysis of ruminal and GI studies to maximize ability to compare across this and previous studies.
The study was able to recover 98 to 99% of all OTUs calculated at 0.03 dissimilarity, as determined by Good's coverage estimator. Rarified samples to 25,000 sequences/sample demonstrated adequate depth, as demonstrated in previous studies [27]. However, greater coverage may be warranted in certain cases, as noted by an in silico study in which the microbial diversity in the rumen was compared against all curated 16S rRNA gene sequences in the RDP database [28]. In this comparison, coverage up to 80,000 reads/sample was estimated to be minimal in order to adequately capture the OTUs at species level (0.03 phylogenetic distance). However, the percent coverage metric used for this conclusion was not the same metric used in this study or in that of Jami and Mizrahi (2012) [27], and thus the requisite number of reads cannot be directly compared.
Greater sequencing depth and number of animals per group in the present study, compared to an earlier attempt to correlate microbiome composition with feed efficiency phenotypes [29], allowed numerous taxonomic and OTU classification links to be established in our experiment. Although the present study had greater observed OTUs, more in keeping with other estimates of rumen microbial richness, diversity analyses of the OTUs identified across all samples and groups, revealed no differences in observed OTUs, richness (Chao1), or diversity (Shannon ; Table 1). However, the data cannot be accurately contrasted to the previous study, due to analysis of different 16S variable regions, as noted previously. Moreover, finer comparisons between studies are impeded by the inherent variability observed in ruminal microbial community datasets. Within-sample variation was minimal based on technical replicates from the same rumen fluid sample, at approximately 0.2% of taxa representing 1% of the reads, while animal-to-animal variation in the ruminal bacterial community within feed efficiency group was as large as 3.0% in the current study. It has previously been reported that variation between rumen metagenome profiles of different cattle is greater than of rumen metagenome profiles from repeated samples on the same animal [30]. We had hypothesized  that increases in sample size could overcome such variation within the microbial communities in the rumen, but this was not the case. In addition to the diversity analyses, the weighted and unweighted UniFrac PCoA examining the phylogenetic divergence between the OTUs did not separate or cluster, further supporting the similarities between the microbial communities within each group (Fig 2). It is possible that the diversity analyses may not be the correct level at which to discover microbiome variability impacts on feed efficiency phenotypes. The lack of differences observed between feed efficiency groups at the level of diversity analyses may simply indicate that the important variation in microbial communities lie at a finer resolution. For example, variation among or within specific taxa and OTUs provided by partial 16S sequencing may be more informative rather than the number and diversity of taxa and OTUs. Nevertheless, the lack of a noted change in diversity (Shannon) between the feed efficiency phenotypes was unexpected, given that obesity or extremes in BW have been associated with changes in the microbial ecology and diversity of intestinal microbial populations in humans and mice [6][7]. In contrast, other studies have shown that host specificity of the ruminal bacterial community may account for some of the similarities observed in microbial diversity, as well as contribute to the substantial animal-to-animal variation [31], which may play an important role in this study. There is limited research examining the effect of the rumen microbiome on feed efficiency, especially with regard to the phenotypes delineated in our model (Fig 1). As noted previously, there were no observable differences between the groups when evaluating OTU alpha and beta diversity (Table 1, Fig 2). However, the relative taxonomic abundance may still be affected, which was the case when examining the relative taxonomic profiles at the phylum and subphylum levels ( Table 2). Many of the changes were identified within the phylum Firmicutes, which included families Lachnospiraceae and Veillonellaceae, and genera Acidaminococcus, Dialister, and Anaerovibrio. Interestingly, increases in the abundance of Firmicutes, altering the Firmicutes-to-Bacteroidetes ratio, have been shown to affect energy harvesting and were correlated with increases of fat [23]. Our study found that increases in Firmicutes were associated with animals of greater ADG (Table 2), and many of the genera belonging to Firmicutes were correlated with high ADG ( Table 5). The great abundance of Firmicutes within the rumen suggests that these shifts may play a significant role in affecting feed efficiency.
The rumen content from the four feed efficiency groups revealed significant differences in the relative abundance of specific taxa (Table 2). Many members of the family Lachnospiraceae have cellulolytic activity and are closely associated with other cellulose-degrading bacteria [32,33]. Acidaminococcus are amino acid-fermenting bacteria [34] and are also related to the butyrate producing Butyrivibrio, which was also detected from changes in feed efficiency (Table 3). In addition, Anaerovibrio has been associated with succinate and propionate production, as well as lipid hydrolysis and metabolism in ruminant animals [35]. Increases in Dialister populations have been associated with hyposalivation [36], which may play a role in altering the buffering capacity of the rumen and fluid turnover. The phylum Lentisphaerae was also associated with changes in feed efficiency, and differences in Lentisphaerae populations have been observed in the rumen during subacute ruminal acidosis [37]. Additionally, lipopolysaccharide challenge has been shown to linearly decrease Bacteroidetes and, in part, Lentisphaerae populations possibly leading to a decrease in fermentative activity [38]. However, little data has been presented on their role of feed efficiency in the rumen, and their significance and ruminal residence remains to be determined. Differences in many OTU abundances were also detected, specifically in the genera Succiniclasticum, Lactobacillus, Ruminococcus, many Prevotella, and the family Veillonellaceae (Table 3). Succiniclasticum was detected at greatest abundance in the ADG Low -ADFI High group, and its relevance lies in being part of order Selenomonadales and specializing in fermenting succinate and converting it to propionate [39]. Veillonellaceae are also known to produce propionate as a major fermentation product [40]. The observation that Succiniclasticum were at higher levels in the ADG Low -ADFI High group, while Veillonellaceae were decreased in abundance within the same group, may indicate resource competition among the organisms. A higher abundance of propionate-producing bacteria may divert H 2 away from methanogenesis reducing methane emissions [8]. Ruminococci are known to be cellulolytic as well as active in acetate, formate, and hydrogen production [41]. Finally, changes in the abundances of Lactobacillus were not completely unexpected, as these organisms are high when concentrate diets are fed [42]. It must be noted however, that although the changes in taxa and OTUs and their putative functions in the rumen can be correlated with the observed differences in the phenotypes, it is not clear whether changes in the microbiome are contributing to differences in feed efficiency or host factors are driving changes in the microbiome.
Trends were observable within the relative taxonomic abundance of the efficiency groups, but did not follow previously observable taxonomic affects. For instance, Prevotella abundance has been shown to be greater in inefficient animals [43]. Yet, this was not adequately demonstrated in the feed efficiency microbial profiles. Prevotella was in abundance at 57.2, 45.4, 49.6, and 47.6% for the ADG High -ADFI High , ADG High -ADFI Low , ADG Low -ADFI Low , and ADG Low -ADFI High groups, respectively (Table 2). Our study did however parallel the accepted abundance of Prevotella seen in other studies regarding ruminants at 42% to 60% of the bacterial 16S rRNA gene sequences [44,45]. Prevotella has been reported as the most abundant ruminal genus, which is in agreement with our study [45].
The unassigned taxa accounted for approximately 9% of the reads. This data has yet to be resolved, but may be important in understanding the entire ruminal microbiome. Species associated with these unassigned taxa may play a yet to be established, important role in feed efficiency. A community-based microbial database would be the most promising option in determining an accurate microbiome. In such cases, not only would database hits increase, limiting unassigned taxa, but accuracy in reporting microbial changes and variability from altered parameters would enhance analysis as well. A current study demonstrated the importance of establishing such a system [30]. The authors developed a reference metagenome to compare rumen metagenomic profiles for individual cattle. When the reads from the study were aligned to a rumen metagenome reference, the rumen metagenome profiles were repeatable (P < 0.00001) within sample regardless of location of sampling rumen fluid. In addition, repeatability was approximately 9%, however because of the small number of cattle used in the study, there was a higher standard error. They were also able to detect variability between samples from different cattle, demonstrating the effectiveness and precision of utilizing a community-specific reference microbial database.
Changes in the microbial populations were also evident when examining ADG and ADFI individually or their interaction. Interestingly, the significant changes in microbial abundances were primarily affected by ADG, with many of the same organisms identified in the feed efficiency study. These data permit further evaluation of the microbial abundances, implicating not only the eventual microbial effect in that of feed efficiency, but also its components, relating information on how population shifts affect feed efficiency.
The majority of taxa and OTUs indentified as associating with changes in feed efficiency in this study were related to cellulolytic, fermentative, and metabolic activities in the cattle rumen. Alterations in their abundances therefore correlate to changes in observed feed efficiency phenotypes, and necessitate further study to examine, alter, and optimize microbial populations for improved feed efficiency.
Although others have demonstrated differences in feed efficiency using community-based PCR-denaturing gradient gel electrophoresis (DGGE) and targeted quantitative real-time PCR analysis [46], this was one of the first studies to directly examine the effect of rumen microbial communities and abundances on feed efficiency, ADFI, and ADG of cattle at the 16S level utilizing next-generation sequencing technologies. Additionally, the study had sufficient power and depth, in contrast to previous research. Although, differences in microbiomes among feed efficiency phenotype class could not be detected at the phylum or genus level by PCoA, diversity indices, or observed OTU count, relative taxonomic abundance and OTU classifications indicated many significant changes in microbial populations as a function of feed efficiency. These populations are also responsible for many of the transformations known to occur in the rumen. However, other analyses may better determine microbial effects on feed efficiency. A metagenomic approach could allow for greater resolution analyses, for example at the level of cellulase genes present and/or expressed [47]. Finally, analyzing the microbial communities throughout the lower GI (jejunum, cecum, and colon) as a function of feed efficiency should also aid in understanding the complete relationship between microbial communities and feed efficient cattle.