Impact of the rumen microbiome on milk fatty acid composition of Holstein cattle

Background Fatty acids (FA) in bovine milk derive through body mobilization, de novo synthesis or from the feed via the blood stream. To be able to digest feedstuff, the cow depends on its rumen microbiome. The relative abundance of the microbes has been shown to differ between cows. To date, there is little information on the impact of the microbiome on the formation of specific milk FA. Therefore, in this study, our aim was to investigate the impact of the rumen bacterial microbiome on milk FA composition. Furthermore, we evaluated the predictive value of the rumen microbiome and the host genetics on the composition of individual FA in milk. Results Our results show that the proportion of variance explained by the rumen bacteria composition (termed microbiability or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{B}^{2}$$\end{document}hB2) was generally smaller than that of the genetic component (heritability), and that rumen bacteria influenced most C15:0, C17:0, C18:2 n-6, C18:3 n-3 and CLA cis-9, trans-11 with estimated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{B}^{2}$$\end{document}hB2 ranging from 0.26 to 0.42. For C6:0, C8:0, C10:0, C12:0, C16:0, C16:1 cis-9 and C18:1 cis-9, the variance explained by the rumen bacteria component was close to 0. In general, both the rumen microbiome and the host genetics had little value for predicting FA phenotype. Compared to genetic information only, adding rumen bacteria information resulted in a significant improvement of the predictive value for C15:0 from 0.22 to 0.38 (P = 9.50e−07) and C18:3 n-3 from 0 to 0.29 (P = 8.81e−18). Conclusions The rumen microbiome has a pronounced influence on the content of odd chain FA and polyunsaturated C18 FA, and to a lesser extent, on the content of the short- and medium-chain FA in the milk of Holstein cattle. The accuracy of prediction of FA phenotypes in milk based on information from either the animal’s genotypes or rumen bacteria composition was very low. Electronic supplementary material The online version of this article (10.1186/s12711-019-0464-8) contains supplementary material, which is available to authorized users.


Background
The bovine rumen provides an anaerobic environment for the acquisition and promotion of a complex microbial community that consists of bacteria, archaea, fungi and protozoa, with bacteria being the most abundant [1]. The relative microbiome composition may differ between cows [2]. The rumen microbiome composition has been shown to have an important impact on the production of methane. Difford et al. [2] reported that variation in the relative rumen microbiome composition could explain up to 13% of the variation in total methane emission from dairy cattle and that the microbial community in the cows' genetic make-up is associated with methane emission. The function of the metabolic activity of these microbial symbionts is to digest complex fibrous feed substrates into volatile fatty acids (VFA) and microbial proteins that are used directly by the host for maintenance, growth and lactation [3]. Furthermore, these fermentation products have an effect on milk composition [4,5], thus the rumen microbial content has the potential to influence variation in host phenotypes [6].
Bovine milk contains many different fatty acids (FA), from essential FA such as linoleic (C18:2) and α-linolenic (18:3) FA to human health promoting FA such as conjugated linoleic acid (CLA) (C18:2) [7]. The short and medium chain FA (C4:0 to C14:0) are synthesized in the mammary gland. The long chain FA derive mainly from the feed but can be modified in the rumen and the mammary gland [8]. C16:0 can be derived both from feed and de novo synthesis in the mammary gland. To date, studies have focused on the genetic factors that affect the synthesis of FA in the milk and have shown that they are the main factors that influence the synthesis of saturated FA, whereas feed is more important for the synthesis of unsaturated FA [9][10][11]. Within the context of the digestion of feed by the microbial community in the rumen, it is interesting to investigate how much of the variation in milk FA can be explained by variation in the rumen microbiome.
Thus, our study focused on the influence of bacterial composition in the rumen microbiome on the synthesis of specific milk FA in Holstein cows by evaluating the value of the rumen microbiome in addition to the host genetics in predicting differences in the composition of individual FA phenotypes in milk.

Animals
Morning milk samples were collected from 339 Holstein cows from three different herds. Each herd had one to three visits within a short time period, during which rumen samples were collected from individual cows immediately after milking (one sample per cow). The cows were kept indoors, fed a total mixed ration (TMR) ad libitum, and milked individually using automated milking systems. They were divided into two groups according to parity: 120 cows for parity 1 and 209 cows for parity 2, and days in milk (DIM) ranged from 3 to 877. For all herds, a standardized TMR recipe was provided, which consisted primarily of rolled barley, corn silage, grass clover silage, rapeseed meal, soybean meal, and up to 3 kg of concentrate supplement given during milking. The procedures for collecting biological samples from the animals were approved by the Danish Veterinary and Food Administration under the Ministry of Environment and Food of Denmark and performed in accordance with the National Guidelines for Animal Experimentation and the guidelines of the Danish Animal Experimental Ethics Committee.

Rumen sampling, DNA extraction and analysis of the bacterial community
Rumen sampling, DNA extraction, library generation and analysis of the bacterial community were performed as previously described in detail by Difford et al. [2]. Briefly, approximately 40 g of rumen content (both liquid and particle matter) were sampled using a flora rumen scoop and a representative subsample was immediately frozen at − 80 °C for DNA extraction [14]. The profile of the bacterial community was assessed by sequencing the V1-V3 region of the 16S rRNA gene. DNA extraction, sequence library construction and paired-end sequencing were performed by a commercial company (GATC Biotech, Constance, Germany). After quality control and processing of sequence reads, they were clustered into operational taxonomic units (OTU) based on a similarity threshold of 97% using the LotuS pipeline [15] via the UPARSE algorithm [16]. Low quality read ends were trimmed, and reads were quality filtered based on the following criteria: (1) a minimum average phred score of 27 and (2) a maximum accumulated error (E max ) of 0.75, which represents the average number of errors per read, i.e. on average 0.75 wrongly sequenced basepairs per read. In order to focus on representative samples that had a sufficient sequencing depth to be able to describe less abundant taxa as well, we removed samples with less than 50,000 reads after quality control [17]. Finally, OTU that contained less than 10 sequences were filtered from the OTU table. Taxonomy was assigned to each OTU using the RDP classifier with a confidence level of 0.8 using greengenes (gg_13_8_otus) as the reference database (http://www.metag enomi cs.wiki/tools /16s/ qiime /insta ll/green genes ).

Genotyping
Genomic DNA was extracted from ear tissue of the 339 Holstein cows and genotyped with the Bovine single nucleotide polymorphism (SNP) 50 beadchip (http:// www.illum ina.com/Docum ents/produ cts/datas heets / datas heet_bovin e_snp5O .pdf ) using an Illumina ® Infinium II Multisample assay device. The iScan and the Beadstudio version 3.1 software were used to scan the SNP chips. SNPs were selected based on the following quality parameters: a minimum call rate of 80% for individuals and a minimum call rate of 95% for loci. SNPs with a minor allele frequency (MAF) lower than 1% were excluded. The quality of the SNPs was assessed using the Illumina GenCall data analysis software. Individuals with average GenCall scores lower than 0.65 were excluded as described in [18]. SNP positions were based on the Bos taurus genome assembly (ARS-UCD1.2) [19]. Finally, 39,121 SNPs remained for analysis.

Cleaning of the animal data
We retained only the cows that had matched rumen samples, milk samples and genotypes (n = 310), and removed from these, those that were in lactation for more than 400 days, which resulted in 292 animals available for the analysis.

Statistical analysis Cluster analysis
In order to detect (biological) patterns in the bacterial OTU, a cluster analysis was performed based on the count data of the individual OTU. In total, 8515 bacterial OTU with 97% sequence similarity were identified. OTU with a total sequence count less than 10 were removed from the data, which resulted in 3055 bacterial OTU available for this analysis. The sequence count data of the bacterial genomes were centered (by subtracting the column means from their corresponding columns) and scaled (by dividing the (centered) columns by their standard deviations), resulting in the B matrix of 292 (number of cows) by 3055 (number of OUT). The distance matrix for cows based on B was calculated using the Dist function of the R package amap with the Pearson distance measure (https ://cran.r-proje ct.org/web/packa ges/amap/ index .html). Hierarchical clustering based on the distance matrix was calculated using the hclust function of the R package stats (R version 3.2.3).

Variance estimation
To estimate the variance explained by the host's genetics and rumen bacteria, we used the REML approach in DMU [20]: where y ijkl is the FA phenotype of individual l , µ is the fixed mean effect, herd i is a fixed effect ( i = 1, 2, 3 ), parity j is a fixed effect ( j = 1, 2 ), b 1 and b 2 are the regression coefficient for DIM k DIM k is a covariate of days in milk (d3-d398) and is modelled according to Wilmink [21], g i is the random additive genetic effect of the animal, m i is the random effect associated with bacterial composition, and e ijkl is the random residual effect. The random effects are assumed to be independent normally distributed values described as follows: g ∼ N 0, Gσ 2 g , m ∼ N 0, Mσ 2 m and e ∼ N 0, Iσ 2 e . The genomic relationship matrix (GRM) G was computed using all the SNPs based on the first method described by VanRaden [22]. The metagenomics relationship matrix ( M ) was computed based on the bacteria counts as previously described by Difford et al. [2]: M = BB ′ /c , where B is the centered and scaled bacteria count matrix, and c is the number of bacterial OTU.
The proportion of the total variance explained by the genetic effect of the animal ( h 2 g ) was defined as: where σ 2 g is the genetic variation, and σ 2 e is the residual variation.
The proportion of the total variance explained by the effect associated with bacteria count ( h 2 B ), which is also called microbiability [2], was defined as: where σ 2 m is the bacterial variation, and σ 2 e is the residual variation.

Cross-validation study to compare models
To assess the influence of the animal's genome and of the metagenome on milk and FA traits, we performed a cross-validation study by using a two-step procedure. In the first step, phenotypic records on milk FA, PP and FP were adjusted for relevant factors (herd, parity, DIM) using the following linear model: where y ijkl is the phenotype of individual l , and the factors are as for Model M1.
In the second step, we built a model in several steps for each adjusted phenotype (i.e. residuals obtained from Model M3) using the following linear mixed models: where ỹ ij is the adjusted phenotype for each milk or FA trait, g j is the random additive genetic effect of the animal, m i is the random effect associated with bacteria count, m i × g j is the Hadamard product between m i and g i , and e i ( e ij ) are the random residual effects. The random effects are re-estimated and are assumed to be independent normally distributed values described as follows: g ∼ N 0, Gσ 2 g , m ∼ N 0, Mσ 2 m and e ∼ N 0, Iσ 2 e . The analysis was performed using a genomic feature best linear unbiased prediction (GFBLUP) for each model separately, as implemented in the R library qgg (http://psoer ensen .githu b.io/qgg/) [23,24]. The features in Models M5, M6 and M7 are the bacteria count data and the genetic data (39,121 SNPs). The ability of each of the models (M4-M7) to predict phenotypes was assessed using a cross-validation procedure. The validation set was generated by randomly sampling 50 out of 292 animals and the remaining 242 animals were used as the training set. In the validation procedure, we estimated the model parameters based on the observations of the cows in the training data and predicted the phenotypes of the cows in the validation data. Then, we calculated the correlation between the predicted phenotype and the observed phenotype. This was repeated 10 times and the predictive ability of a model was defined as the average correlation of these 10 validations.
The predictive ability of each model was compared using a Welch's t test (i.e. unequal variance t test) as implemented in R (R version 3.2.3). We tested whether the average correlation of these 10 validations was significantly higher than the average correlation obtained from another model. We performed the following comparisons: M4 versus M6 (bacteria are important); M5 versus M6 (genetics is important); M6 versus M7 (interaction between bacteria and genetics is important). A Bonferroni corrected P value lower than 0.05 was considered significant (n = 57 comparisons).

Discussion
In this study, we evaluated the influence of the bacterial composition of the rumen microbiome on the biosynthesis of specific FA in bovine milk. FA composition in milk is the result of a complex process in which FA originate from different sources. Short chain and medium chain FA are synthesized de novo in the mammary gland, while long chain FA derive mainly from feed [7]. In general, the heritability of de novo synthesized FA is high, whereas that of long chain FA is low to moderate [8]. Poulsen et al. [12] reported that the feeding regime has a clear influence on the composition of long chain FA, but also that they are modified considerably in the rumen and by desaturation in the mammary gland. Even if cows have the same feeding regime, individual feed intake per cow varies because its energy requirement differs depending on parity and lactation stage. Furthermore, cows sort out their feed even when fed a TMR. Thus, one can argue that the microbial composition in the rumen could be a proxy for feed type, rather than influencing milk fat composition. However, Sasson et al. [25] showed that abundance of some bacteria in the rumen is heritable [25]. Thus, variation in the rumen microbiome may influence the biosynthesis of milk FA and the modifications of feed-derived FA.

Core microbiome
Henderson et al. [1] have suggested that a 'core microbiome' of dominant bacteria exists in the rumen at the genus level or higher and consists of Prevotella, Butyrivibrio, Ruminococcus as well as unclassified Lachnospiraceae, Ruminococcaceae, Bacteroidales and Clostridiales. Our results are in line with the results of Henderson et al. [1], i.e. they show that the clusters of OTU with the highest counts were mainly assigned to Prevotella and the unclassified group Bacteroidales, whereas the clusters with the lowest count data also contained the (unclassified) groups Clostridiales, Lachnospiraceae, Ruminococcaceae, Ruminococcus, Treponema, and RF39.

Heritability, microbiability and sample size
In our study, estimated heritabilities for fat % and protein % were low (0.19, and 0.18, respectively) compared to previous reports on fat %, i.e. 0.24 in [9] and 0.47 in [10] and on protein %, i.e. 0.47 in [26] and 0.53 in [27]. The discrepancy between our current results and those in the literature could be due to the use of different models to estimate the variances or simply to the relative small sample size used here. Ideally, the variance parameters in this study would be estimated jointly, combining models M1 and M2, but the sample size in our study limits the possibility of estimating parameters from the full model in which both G and M are fitted simultaneously. Nevertheless, we believe that our results show that some of the milk FA are affected by the rumen microbiome, and that our study could serve as a basis for future studies on the influence of the microbiome on milk FA composition.
Our results show that the short and medium chain FA have a stronger genetic component, compared to the long chain C18 FA, which is in line with the expectation that long chain FA derive from feed. However, our results do not confirm the trend that saturated FA have a higher heritability than unsaturated FA as previously reported [9,10]. Comparison of our results to an earlier independent study on milk FA composition in Holstein cattle [9] indicates that the estimated heritabilities are in line with those in [9] except for C14:0 and C18:3 n-3, for which we detected no genetic component. Interestingly, both C14:1 cis-9 and C16:1 cis-9 had high heritabilities, which reflect a strong genetic influence related to the desaturase activity in the mammary gland. The conversion of C14:0 and C16:0 into their mono-unsaturated counterparts is catalyzed by stearoyl-CoA desaturase 1 (SCD1). The SCD1 gene is located on B. taurus chromosome (BTA) 26 and is associated with milk FA composition [28].
We detected no influence of the rumen bacteria composition on the C6-C12 and C16 FA group, which suggests that they have little or no effect on de novo synthesis in particular, but also on C16:0 derived from feed. This is further reflected by the absence of influence on C14:1 cis-9, and C16:1 cis-9. In general, the influence of bacteria in the rumen on the formation of milk FA is more pronounced for the odd-chain FA (C15:0, C17:0) and polyunsaturated C18 FA (C18:2 n-6, C18:3 n-3 and CLA cis-9, trans-11), as expected based on the origin of the FA in milk [29]. These FA originate partly from hydrogenation processes in the rumen by specific microorganisms that affect feed-derived C18 FA and partly from odd-chain FA that are mainly synthesized microbially in the rumen [30].
For C17:0 and many of the C18 FA, the variation in rumen bacterial content explained more of the host phenotypic variation than the host genome variation, albeit marginally, and with a large standard error. Although variation in these FA is often associated with differences in the feed, in our study, feed variation was minimized since animals were fed a standardized TMR. Furthermore, these findings confirm the existence of variation in the Table 3 Mean predictive ability a (PA) and between brackets root mean square error of the prediction (RMSE) for models M4, M5, M6 and M7 and P value for the differences in PA between M4 versus M6, M5 versus M6, and M6 versus M7 *Significant after Bonferroni correction for multiple testing at the P < 0.05 level (n = 57 comparisons) a Average correlation between predicted and real FA phenotype in the validation set over 10 replicates of the validation  [31]. That study also showed that the abundance of 49 genera of the pig gut microbiome was heritable [31].

Prediction of fatty acids in the milk
Recent findings in dairy cattle suggest that the host genetics has some influence on the rumen microbiome [25,32,33]. There is considerable interest in understanding the interactions between the host and the rumen microbiome for directing desired changes to the host phenotype. If associations between heritable rumen bacteria and phenotypes of the host exist, there is scope for research in selective breeding. Selective breeding takes generations to induce genetic changes in livestock populations, whereas interventions on the rumen microbiome can act rapidly within a few generations. Rumen transfaunation using the cud from a healthy donor cow to treat a recipient cow with a digestive disorder has long been an effective method of implementing a desired change to the host cow phenotype [34]. Recent work with repeated inoculations in beef heifers on a poor strawbased diet with rumen contents from bison has shown that the rumen bacteria and the host metabolic activity are altered, which result in increased protein and nitrogen retention [35].
It is relevant to investigate whether the bacterial and genetic information can contribute to predicting the milk fat composition of cows. In general, we found that the predictive value for FA was low, which could be due to the relatively small sample size used in our study. Genotype information is of higher value for the prediction of specific FA in the milk than information on the rumen bacterial composition. We found that genotype information was relevant for the prediction of the medium chain FA C6 to C12, C14:1 cis-9 and C16:0, C16:1 cis-9. This is in line with the fact that these FA are fully or partly de novo synthesized or reflect desaturase activity in the mammary gland.
Interestingly, information on the rumen bacterial community was better at predicting the C15:0 and C18:3 n-3 content in milk (r = 0.29-0.34) than genotype information. The predictive value of rumen bacteria for C15:0 and C18:3 n-3 was in line with that of gut bacteria in pigs for feed intake, daily gain and feed conversion (r = 0.33-0.41) [30] and that of methane emission in dairy cows using full rumen metagenomics sequences (r = 0.466) [35]. The estimated h 2 g and h 2 B for C15:0 were similar, which is as expected since C15:0 is known to be both synthesized de novo and derived from blood [36,37]. However, information on the microbiome is preferable for predicting the C15:0 content in milk. Regarding the prediction of C18:3 n-3 content in milk, C18:3 n-3 derives from the feed and thus depends highly on the feeding regime of the cow [37] and on the microbiome, which plays an important role in degrading the feed in the rumen and thereby regulates the C18:3 n-3 content in milk through rumen biohydrogenation. Our results also show that, when h 2 g or h 2 B are equal to 0 (Table 2), the mean predictive value is negative (Table 3), which indicates that the correlation between the observed and predicted phenotypes are biased downwards.

Conclusions
Our findings show that variation in rumen microbiome composition has a pronounced influence on the content of odd chain FA and the polyunsaturated C18 FA, and to a lesser extent, on the content of the short and medium chain FA in milk. The prediction of individual FA content in milk based on information of either the animals' genotypes or the rumen bacteria was low. The results can be explained from a biological point of view, e.g. the h 2 g for the saturated FA was generally higher than h 2 B and the predictive ability of the model fitting the GRM was better for the saturated FA than the model fitting the microbial relationship matrix. Nevertheless, standard errors of the heritability estimates were large and, in some cases, the predictive values tended to be biased downwards, which indicated that it is necessary to increase the sample size in a future study.

Additional file
Additional file 1. Annotation of the most abundant bacteria OTU.