Investigation into the effect of divergent feed efficiency phenotype on the bovine rumen microbiota across diet and breed

The relationship between rumen microbiota and host feed efficiency phenotype, for genetically divergent beef cattle breeds is unclear. This is further exacerbated when different growth stages, chemically diverse diets and production systems are considered. Residual feed intake (RFI), a measure of feed efficiency, was calculated for individually fed Charolais (CH) and Holstein–Friesian (HF) steers during each of four 70-day (excluding adaptation) successive dietary phases: namely, high-concentrate, grass silage, fresh zero-grazed grass and high-concentrate again. Rumen fluid from the ten highest- (HRFI) and ten lowest-ranking (LRFI) animals for RFI, within breed, during each dietary phase was collected using a trans-oesophageal sampler and subjected to 16S rRNA amplicon sequencing and metabolic profiling. The datasets were analysed to identify microbial and rumen fermentation markers associated with RFI status. Age, dietary phase and breed were included in the statistical model. Within breed, for each dietary phase, mid-test metabolic weight and average daily gain did not differ (P > 0.05) between HRFI and LRFI steers; however, for the initial high-concentrate, grass silage, fresh grass herbage and final high-concentrate dietary phases, HRFI HF steers consumed 19, 23, 18 and 27% more (P < 0.001) than their LRFI counterparts. Corresponding percentages for CH HRFI compared to CH LRFI steers were 18, 23, 13 and 22%. Ten OTUs were associated with RFI (q < 0.05) independent of the other factors investigated. Of these Methanomassiliicoccaceae, Mogibacteriaceae and the genus p-75-a5 of Erysipelotrichaceae and were negatively associated (q < 0.05) with RFI. The results gave evidence that microbial species could potentially be an indicator of RFI in ruminants rather than broader microbiome metrics; however, further research is required to elucidate this association.

The agricultural sector contributes to about 10-12% of global anthropogenic emissions with enteric methane accounting for 32-40% of this 1 . As a result, improving the efficiency of cattle production will be key to adhering reduction in individual national greenhouse gas emissions under the United Nations Framework Convention on Climate Change Paris Agreement 2 . Beef production can use up to five times more biomass for producing 1 kg of animal protein than dairy 3 contributing to both an environmental and economic challenges. Feed provision is classed as the single largest variable cost associated with beef production; reducing feed costs through identification of more feed efficient cattle should increase profitability and environmental and economic sustainability of beef farming enterprises 4 . Improvement of feed efficiency (FE) in cattle is not only critical to improve environmental but also for economic sustainability 5,6 . Residual feed intake (RFI) is an established measure of FE in beef cattle, defined as the difference between an animal's actual compared to its predicted feed intake, usually based on weight and growth 7 . Animals with a low RFI (LRFI) phenotype, eat less than expected and are deemed to be Scientific RepoRtS | (2020) 10:15317 | https://doi.org/10.1038/s41598-020-71458-0 www.nature.com/scientificreports/ efficient, typically consuming up to 20% less feed than their inefficient high RFI (HRFI) contemporaries, while supporting the same bodyweight and tissue growth [8][9][10] . The rumen of cattle is a complex fermentation vat inhabited by a multitude of microbial species acting synergistically to convert plant materials into nutrients, primarily short chain fatty acids (SCFAs) and microbial protein 11 . This process enables ruminants to convert human indigestible plant polysaccharides into high quality meat and milk products suitable for human consumption. The predominate SCFAs, acetate, propionate and butyrate, provide up to 70% of the host's energy requirements 12,13 . Determination of the relationship between the rumen microbiome and host FE has the potential to facilitate selection of livestock with enhanced nutrient utilization or enable the manipulation of the rumen microbiome to enhance its energy harvesting capacity [14][15][16] . Methanogenic archaea utilize by-products of rumen fermentation to produce methane which is an energetically wasteful process for the host animal 17 and a potent greenhouse gas 18 . Diversion of these by-products away from methanogensis and towards alternative biochemical pathways, more beneficial to the host ruminant has been hypothesized as one mechanism to influence host FE 11,14 . Therefore, enhancing feed efficiency has the dual benefit of reducing both feed costs and anthropogenic methane production 19,20 . Currently, however, there is a lack of understanding of the ruminal metabolic and physiological mechanisms underlying variation for RFI in cattle and how this might be impacted upon by genetic makeup of the animal as well as the chemical composition of the diet fed.
The main impediment to genetic progress and adoption of selection for feed efficiency is the logistics and expense of measuring individual animal intake necessary to accurately identify feed efficient animals 21 . Kelly et al. 8 and Coyle et al. 22 have demonstrated that RFI phenotype is moderately repeatable for cattle offered the same diet over different phases of their lifetime; however, repeatability is reduced when contrasting diets are fed successively (i.e. forage versus concentrate-based diets), as per commercial farming practice 23 . These results may have implications for beef genetic evaluations and breeding programmes. Shabat et al. 14 reported that rumen microbiome genes and species could predict the variation in a dairy cow's RFI phenotype (91% accuracy). Thus, if a microbial marker for FE could be deduced there is potential to utilize it for genomic selection of feed efficient cattle. However, the study by Shabat et al. 14 was conducted across a single timepoint, using Holstein Friesian dairy cattle and on a concentrate diet, leaving the resilience of microbial markers across breed, dietary phase and physiological age of the ruminant yet to be investigated.
The current study investigated the effect of RFI phenotype of both beef and dairy bred steers across common dietary phases of Irish pastoral-based beef production systems on the rumen microbiome and its metabolites. This aims to elucidate if particular microbial taxa influence host FE phenotype independent of stage of production, diet and host breed, leading to the potential for rumen microbiome manipulation or selection for FE based on microbial markers.

Results
Animal performance. Descriptive  Sequence analysis. Amplicon sequencing generated 84,462,720 total reads giving an average of 272,460 ± 69,596 reads per sample. This reduced to 217,817 ± 55,519 when sequences were merged and quality filtered. The average number of counts per sample that were assigned to an OTU (post filtering) was 175,304 ± 74,272. Positive control samples were subjected to the same bioinformatics pipeline as rumen 16S rRNA libraries; generating an average of 203,612 ± 19,378 sequence reads, 141,873 ± 15,724 merged sequences and 132,567 ± 14,397 final counts per sample. The positive controls were positively correlated with the theoretical composition of the microbial community standard DNA (P < 0.05) deeming the laboratory and informatic methodology satisfactory. Negative controls generated 215.3 ± 107 final counts per sample. Table 1. Dry matter intake (DMI) (kg/day), residual feed intake (RFI), average daily gain (ADG) (kg) and Mid-test metabolic weight (MBW) (kg) for CH and HF steers ranked low (LRFI) and high (HRFI) RFI offered; high-concentrate, grass silage and zero-grazed grass and a second high-concentrate diet. P values are derived using a Wilcoxon rank sum test to assess the differences between treatments. www.nature.com/scientificreports/ Rumen fermentation profile. Rumen SCFA analysis showed that LRFI CH steers had lower concentrations of propionate and total SCFA in their rumen fluid in comparison to HRFI CH steers (P < 0.05) when offered a high-concentrate diet during their growth phase ( Table 2). In HF steers offered the same diet, LRFI steers had a higher concentration of acetate, increased total SCFA concentration in comparison to HRFI steers ( Table 2). No differences in any metabolite concentrations existed between HRFI and LRFI steers offered a grass silage diet for either the CH or HF breed (P > 0.05) ( Table 2). LRFI CH steers had reduced butyrate concentration (P < 0.01) HRFI CH steers offered a zero-grazed grass diet, while there was an increase in acetate, propionate, isobutyrate, valerate and total SCFA concentrations in LRFI steers in comparison to HRFI HF steers offered a zero-grazed grass diet (P < 0.05) ( Table 2). There was an increase in isovalerate concentrations in LRFI CH steers in comparison to HRFI CH (P < 0.1) and no differences in SCFA concentrations existed between LRFI and HRFI HF steers fed a high-concentrate diet during the finishing phase (P > 0.05, Table 2). Overall no consistent rumen fermentation profile was observed between divergent RFI phenotypes across the phases investigated in the study ( Table 2). The association between SCFA concentration and RFI was investigated, highlighting dietary phase, age and breed as potential confounding factors. No association between RFI and rumen fermentation profile was observed (q > 0.05).
Ordination analysis of microbial community structure. Principal component ordination analysis (PCoA) analyses ( Fig. 1) were used to visualise the similarities and diversity of bacterial and archaeal populations of HRFI and LRFI steers across dietary phase and breed. The plot showed a lack of clustering based on RFI phenotype. Samples clustered tightly based on diet offered and dietary phase. Pairwise comparisons using Bray Curtis and Unifrac dissimilarity metrics between RFI phenotype within each diet and breed and comparisons between breed with the same RFI phenotype within dietary phase are summarised in Table 3 and visualised in Figs. 1 and 2. In summary, bacterial community structure did not differ between RFI phenotype in steers in any experimental phase irrespective of breed (P > 0.05) ( Table 3).
Alpha diversity. Similar to the beta diversity results, there was no difference in alpha diversity observed in the rumen microbiome of steers with divergent RFI phenotypes in any phase or for either breed (Table 4). Maas-lin2 was also used to investigate if alpha diversity index showed an association with RFI phenotype. However, no alpha diversity metrics were significantly associated with RFI regardless of phase, breed or age (q < 0.05).

Effect of RFI phenotype on the archaeal and bacterial populations. No differentially abundant
microbial taxa were identified between RFI phenotypes. In total there were 25 OTUs associated with RFI (q < 0.05) with 10 of these associated with RFI independent of any other factor investigated (i.e. Phase, Breed, Age). These are summarised in Table 5. Table 2. Short chain fatty acid concentrations (mmol/L) (A) CH and (B) HF steers ranked low (LRFI) and high (HRFI) RFI offered common dietary phases of Irish pastoral-based beef production systems; highconcentrate, grass silage and zero-grazed grass and a second high-concentrate diet during finishing phase. P values are derived using a Wilcoxon rank sum test to assess the differences between treatments.

Discussion
Feed efficiency is an environmentally and economically important driver sustainable beef production 24 . The rumen microbial community mediates energy available to the host ruminant through its fermentative activity, which suggests it plays a fundamental role in feed efficiency 21 . Previous studies have identified that host FE phenotype influences the rumen microbiome and rumen fermentation 14,25,26 . However, a key challenge is to identify whether these relationships are consistent across animal breed, physiological age and dietary phase. This study therefore examined the effect of RFI phenotype on ruminal fermentation profile and bacteria and archaeal populations in CH and HF steers, across four dietary phases. Ruminal SCFA parameters are indicative of active bacterial fermentation and host rumen epithelial absorption 27 . It is presumed that due to the central role of ruminal digestion in the supply of nutrients for post absorptive metabolic processes, that differences in rumen fermentation profile would be observed between HRFI and LRFI phenotypes 14,28 . Despite this, there was no consistent relationship observed across breed and dietary phase between individual rumen metabolites and RFI rank for steers. Furthermore, fermentation profiles of the steers of the same breed and the same feed efficiency phenotype, offered the same diet during different production phases (i.e. high-concentrate 1 and high-concentrate 2) also showed no consistency in fermentation profile. This result is in agreement with previous studies which highlight the lack of consistency of the rumen fermentation profiles of cattle divergent for feed efficiency, as reviewed by Kenny, et al. 21 . Results do however, indicate that high concentrate (C1), grass silage (GS) and zero grazed grass (ZGG) and a second high-concentrate diet (C2). This is based on similarity of OTU composition of each sample calculated using Bray-Curtis similarity index and plotted using Principal component ordination analysis.  14 hypothesised that a LRFI phenotype supports an efficient rumen microbiome which has lower diversity and richness than its HRFI counterparts, which is also confirmed by the work of Guan et al. 31 . The efficient rumen microbiome is less diverse and produces a smaller range of "relevant" output metabolites suitable for energy and carbon channelling to the animal i.e. meeting the ruminant's energetic needs, while lowering methane emissions to the atmosphere 14 . There is some ambiguity, however, in the literature regarding this hypothesis, as it is contradicted by the studies of Paz et al. 25 , Li 16 , Myer et al. 32 and indeed our own previous findings 33 which reported no differences in microbial diversity between HRFI and LRFI animals. In agreement with much of the literature, no decrease in diversity or richness was observed for any steers deemed as feed efficient relative to their inefficient contemporaries of the same breed offered the same diet in the current study. These results indicate that a lower diversity microbiome may not necessarily equate to a more feed efficient ruminant and that microbial diversity is influenced by the chemical composition of the prevailing diet and breed as previously reported in the literature 30 .
Although there have been some reports of differences in overall community structure in relation to RFI 14,31 , research from this study indicates that it is more likely that individual microbial species or strains may impact the efficiency of the animal more than an overall ruminal community shift. The association analysis provided some evidence of the influence of RFI phenotype on rumen microbial populations, across all dietary phases and breed. Rumen methanogen populations have been associated with RFI in ruminants 19,34,35 , with efficient livestock identified to produce lower methane yields compared to their inefficient contemporaries 19,20 . Methanogen populations in the rumen are responsible for methane production; an energetically wasteful process for the host ruminant, therefore it has been hypothesized that increased abundance of methanogenic archaea reduce the feed efficiency status of the host ruminant 36 . In contrast, in the current study a negative relationship was observed between an archaeal members of the Methanomassiliicoccaceae family exhibited a negative relationship with the RFI trait at genus and OTU level. The Methanomassiliicoccaceae family are obligate H 2 -dependent methylotrophs, utilising methyl groups from methanol and methylamines (mono-, di-, and tri-methylamine) and methyl thiols for the production of methyl coenzyme 37 and provide additional NH 4 + to rumen bacteria. Li 16 found that Methanomassiliicoccales tended to be more abundant in low compared with high RFI cattle. In different experimental ruminant models, McGovern et al. 38 and Danielsson et al. 39 Table 3. Effect of residual feed intake (RFI) phenotype and breed on beta diversity for Charolais (CH) and Holstein Friesian (HF) steers ranked low (LRFI) and high (HRFI) RFI offered common dietary phases of Irish pastoral-based beef production systems; high-concentrate, grass silage and zero-grazed grass during growth phases and a second high-concentrate diet. Distance based permutation multivariate analysis of variance (PERMANOVA) was performed to test the null hypothesis that there were no differences in the microbial community structure across treatment at a significance level of P = 0.05 based on 999 permutations.  42,43 . Methanomassiliicoccaceae, similar to the RO clade, thrive in a lower H 2 environment. Therefore, it is more likely that the composition of the archaeal community rather than total methanogen abundance may have greater significance with respect to feed efficiency 38 . It has also been hypothesized that the Methanomassiliicoccales order may have a protective effect in the rumen of LRFI steers, eliminating the potential negative effects of high methylamine concentrations 44,45 and also providing additional NH 4 + for nitrogen cycling 46 . The association between bacterial genera and RFI less clear with very few of the bacterial species shown to be associated with the trait. A Mogibacteriaceae-affiliated unnamed OTU which was found to be negatively correlated with RFI, Mogibacterium are asaccharolytic, Gram positive bacteria 47 and are hypothesised to be involved with ammonia assimilation 48 . Mogibacterium utilises ammonia and produces phenylacetate, a precursor for phenylalanine 47 . Previously a Mogibacteriaceae-affiliated genus was reported to be associated with feed efficiency in beef cattle with multiple genetic breeds 32 . Despite the lack of information pertaining to this bacterial genera and its contribution to increased FE, there is evidence of increased nitrogen digestion in LRFI cattle 19,49 . Abundance of members in this family were negatively correlated with body mass index (BMI) in humans 50,51 . Cattle with LRFI phenotype have been observed to have leaner carcasses 52 , suggesting the higher abundance of Mogibacteriaceae in L-RFI individuals may correspond to a lean phenotype. It was also observed that an OTU Figure 2. Principal component ordination analysis (PCoA) plot indicating similarity bacterial and archaeal community of Charolais (CH) and Holstein Friesian (HF) divergent for residual feed intake (RFI) steers offered; high concentrate (C1), grass silage (GS) and zero grazed grass (ZGG) and a second high concentrate diet (C2). This is based on similarity of OTU composition of each sample calculated using Weighted Unifrac distance metric and plotted using principal component ordination analysis.  53 also investigating the impact of breed on the microbiome and its association with feed efficiency, it was reported that abundance of p-75-a5 was elevated in LRFI Charolais steers in comparison to their HRFI contemporaries, however, this elevation was not observed in the other breeds investigated.

Conclusion
Research to date suggests there is potential to enhance nutrient utilization from feed and improve FE by altering the rumen microbiome. However, there is a lack of information surrounding the effects of RFI phenotype on the rumen microbiome, particularly during varying stages of production when cattle, are offered contrasting diet types. The results of rumen metabolite profiling were inconsistent and showed little to no consistency with RFI phenotype. This study provides evidence that the composition of the methanogenic community present in the rumen may prove to be an indicator of host RFI; however, further comprehensive investigation focusing on archaeal communities at multiple time points will be required to fully elucidate relationship. These results emphasise that targeted metabolite profiling or molecular surveying of the phylogenetic diversity within the rumen may not comprehensively identify microbial or fermentation factors which influence RFI phenotype. This indicates that global high depth metagenomic shotgun sequencing approach is most likely required to elucidate the intricacies of both the diversity and functionality of the rumen microbiome in relation to FE.

Materials and methods
All procedures involving animals in this study were approved by the Teagasc Animal Ethics Committee and conducted under an experimental license approved by the Irish Health Products Regulatory Authority, project authorisation number AE19132/P029. All methods and experimental protocols in this study were performed in accordance with relevant guidelines and regulations.
Animal model. This experiment was conducted as part of a larger study designed to examine the within-animal repeatability of intake, growth and FE between the growing and finishing stages of the lifespan of Charolais (CH) and Holstein-Friesian (HF) beef steers offered either the same diet or diets contrasting in energy density and chemical composition 10,22,23 . Briefly, 167 cattle comprising of 90 CH and 77 HF were used in this study. Following a dietary adaptation period, individual dry matter intake (DMI) and growth were measured over four 70-day feeding phases (Fig. 3). Mean BW (SD) and age (SD) at the start of feeding phase 1 (high-concentrate 1) were 394 kg (37.5) and 283 days (18.3), and 294 kg (41.8) and 307 days (7.7), for CH and HF, respectively. Corresponding BW at the start of phase 2 (grass silage) were 485 kg (37.6) and 519 kg (38.3), phase 3 (zero-grazed fresh grass) were 516 kg (37.3) and 440 kg (41.8) and phase 4 (high-concentrate 2) were 676 kg (50.0) and 611 kg (49.1). During the two high-concentrate feeding phases 1 ('yearling') and 4 ('two-year old') steers were individually offered the same concentrate diet (860 g/kg rolled barley, 60 g/kg soya bean meal, 60 g/kg molasses, and 20 g/kg minerals and vitamins) ad libitum plus a restricted allowance of grass silage daily. Phase 2 consisted of offering grass silage (precision-chop harvested from a primary growth sward, which comprised mainly of perennial ryegrass) to appetite. During phase 3 steers were individually offered zero-grazed grass (DM 196 g/kg) ad libitum. The grass herbage was harvested (without chopping) twice daily from Lolium perenne dominant swards using a 'zero-grazer' (Model AB70 Zero Grazer, Dromone, Oldcastle, Co. Meath). At the midpoint of each experimental phase after a dietary adaptation perios a single rumen fluid sample was collected from each animal via stomach intubation (Flora Rumen Scoop, Profs-Products, Guelph, Canada) approximately 2 to 4 h post-feeding 20 , for metabolite profiling and microbial analysis. Rumen fluid samples were immediately snap frozen in liquid nitrogen and subsequently stored at -80 °C for molecular analyses.

RFI calculation.
At the end of each dietary phase, average daily live weight gain (ADG) of each steer was computed as the coefficient of the linear regression of BW (kg) on time (days) using the GLM procedure of SAS 9.3 (SAS Inst. INC., Cary, NC). Mid-test metabolic weight (MBW) was determined as BW 0.75 35.5 days prior to the end of the test, which was estimated from the intercept and slope of the regression line after fitting a linear regression through all BW 0.75 observations. RFI was calculated for each animal as the difference between actual DMI and expected DMI. Expected DMI was computed within breed for each animal using a multiple regres- www.nature.com/scientificreports/ sion model, regressing DMI on MBW and ADG. Animals were ranked according to RFI coefficient and the 10 highest HRFI and 10 lowest LRFI ranking animals for RFI, within breed, during each dietary phase were used for this study. For logistical reasons seven samples could not be obtained during the sampling process; leaving a sample size of n = 153 Phase one: 10 HRFI CH, 10 LRFI CH, 10 HRFI HF, 10 LRFI HF Phase two: 9 HRFI CH, 9 LRFI CH, 9 HRFI HF, 9 LRFI HF Phase three: 10 HRFI CH, 9 LRFI CH, 9 HRFI HF, 10 LRFI HF Phase four: 10 HRFI CH, 10 LRFI CH, 10 HRFI HF, 9 LRFI HF.
Rumen fermentation profiling. Short chain fatty acid concentrations in rumen fluid samples were measured using a gas chromatograph (model 3800 Varian gas chromatograph). No results were obtained for eight rumen fluid samples leaving a total of n = 145 (Phase one: 10 HRFI CH, 10 LRFI CH, 10  DNA extraction. Approximately 20 g of frozen rumen liquid sample from each animal was considered as representative. Each sample was homogenised to a fine frozen powder under liquid nitrogen using a pestle and mortar and stored at − 80 °C. Approximately 250 mg of homogenized frozen powder was used for DNA extraction. DNA was extracted using the repeated bead beating and column purification method 54 . DNA extractions were also performed on molecular grade water as a negative control. DNA quality was assessed on an agarose gel. DNA extract yield and purity were assessed with two consecutive readings on the Nanodrop 1000 spectrophotometer. 16S rRNA amplicon library preparation. Amplicon libraries (n = 155) targeting the V4 region of the 16S rRNA gene in bacteria and archaea were generated by PCR amplification. PCR reactions were performed for amplicon libraries with 20 ng of rumen microbial DNA and 515F forward and 806R reverse primers 55  Amplicons were pooled in equal concentrations and gel purified to remove unwanted products using the Qiagen Gel Extraction Kit (Qiagen, Manchester, UK). An extra purification step using the QIAquick purification kit (Qiagen, Manchester, UK) was performed to remove residual agarose. The pooled purified libraries were measured for purity and quantity on the Nanodrop 1000 and further quantified using the KAPA SYBR FAST universal kit with Illumina Primer Premix (Roche Diagnositics, West Sussex, UK). The library pool was then diluted and denatured as recommended by Illumina MiSeq library preparation guide. The sequencing was conducted using 500 cycle MiSeq reagent kits (Illumina, San Diego, CA, USA). Sequence analysis. Two samples (HRFI HF, LRFI CH) from phase 3 were removed from the analysis due to low sequence output. The 16S rRNA gene forward and reverse reads were imported into Qiime2 56 . The DADA2 pipeline 57 was used for detecting and correcting Illumina amplicon sequences, removal of primers and chimeric reads, and assembly into sequence variants (SV)/operational taxonomic units (OTUs) 58 . Taxonomy was assigned using a naïve Bayes classifier trained on the RefSeq database 59 . Sequence files associated with each sample have been submitted to the NCBI Sequence Read Archive (Accession no. PRJNA483745).
Rumen samples were rarefied to a sampling depth of 49,075 sequences. Alpha-diversity metrics investigated included Faith's Phylogenetic Diversity (PD), Evenness, Observed Operational Taxonomy Units (OTUs) and Shannon's diversity index which were calculated using qiime2-q2-diversity. Beta-diversity metrics was calculated using qiime2. Distance metrics calculated included weighted and unweighted UniFrac50 and Bray-Curtis dissimilarity index. Data was visualized using principal coordinate analysis (PCoA) plots and boxplots which were generated using ggplot2 60 within R version 3.5.2.
A Kruskal-Wallis test with Dunn's post-hoc test was used to identity individual taxa differentially represented across treatment groups and assess variations in alpha diversity between treatment groups. Distance based permutation multivariate analysis of variance (PERMANOVA) 61 was performed to test the null hypothesis that there were no differences in the microbial community structure across treatment at a significance level of P = 0.05 based on 999 permutations. Resulting P values were false discovery rate (FDR) corrected 0.05 as the cut off. The corrected P values are presented as q-values.
MaAsLin2 is a statistical framework that performs boosted, additive general linear models between clinical data and bacterial abundance to find associations between clinical data (categorical or continuous) and microbial taxa 62 . MaAsLin2 was used to find associations between bacterial taxa and RFI phenotype. Dietary phase, age (days) and breed were included in the statistical model in order to elucidate microbes which were associated with RFI independent of any other these factors. A FDR P value of < 0.10 was used as significance cut-of. MaAs-Lin2 was also performed to show the potential association between SCFA and alpha diversity metrics with RFI.