Authors: M. Schären,*† J. Frahm,* S. Kersten,* U. Meyer,* J. Hummel,‡ G. Breves,§ and S. Dänicke,*
*Institute of Animal Nutrition, Friedrich-Loeffler-Institute (FLI), Federal Research Institute for Animal Health, Bundesallee 50, 38116 Brunswick, Germany
†Clinic for Ruminant and Swine, Veterinary Faculty, University of Leipzig, An den Tierkliniken 11, 04103 Leipzig, Germany
‡Ruminant Nutrition, Department of Animal Sciences, Faculty of Agricultural Sciences, Georg-August University Göttingen, Kellnerweg 6, 37077 Göttingen, Germany
§Department of Physiology, University of Veterinary Medicine Hanover, Bischofsholer Damm 15, 30173 Hannover, Germany
State of publication: under revision at the Journal of Dairy Science (December 2017) Contribution of authors:
Head of organization and execution: SD, UM, GB, MS
Trial and project design: SD, UM, SK, JF, MS
Trial implementation and sample collection: MS
Sample analysis: MS
Data analysis and interpretation: MS, JF, SK, UM, JH
Writing of manuscript: MS
Revision of manuscript: JF, SK, UM, JH, GB, SD
Interrelations between the rumen microbiota and production, behavioral, rumen fermentation, metabolic and immunological attributes of dairy cows
M. Schären,*† J. Frahm,* S. Kersten,* U. Meyer,* 1 J. Hummel,‡ G. Breves,§ and S. Dänicke,*
*Institute of Animal Nutrition, Friedrich-Loeffler-Institute (FLI), Federal Research Institute for Animal Health, Bundesallee 50, 38116 Brunswick, Germany
†Clinic for Ruminants and Swine, Faculty of Veterinary Medicine, University of Leipzig, An den Tierkliniken 11, 04103 Leipzig, Germany
‡University of Goettingen, Department of Animal Sciences, Faculty of Agricultural Sciences, 37077 Göttingen, Germany
§Department of Physiology, University of Veterinary Medicine Hanover, Bischofsholer Damm 15, 30173 Hannover, Germany
1Corresponding author: ulrich.meyer@fli.de
Submitted August 2017 at the Journal of Dairy Science (currently under revision)
Abstract
Different studies have shown a strong correlation between the rumen microbiome and a range of production traits (e.g. feed efficiency, milk yield and components) in dairy cows. Underlying dynamics concerning cause and effect are, however, still widely unknown and warrant further investigation. The aim of the current study was to describe possible functional interrelations and pathways using a large set of variables describing the production, the metabolic and immunological state as well as the rumen microbiome and fermentation characteristics of dairy cows in early lactation (n = 36, day 56 ± 3d in milk). It was further hypothesized that the feed intake associated behavior may influence the ruminal fermentation pattern and a set of variables describing these individual animal attributes was included. Principal component analysis (PCA) as well as Spearman´s Rank correlations were conducted including a total of 265 variables. The attained plots describe several well-known associations between metabolic, immunological and production traits. Main drivers of variance within the dataset included milk production and efficiency, and rumen fermentation and microbiome diversity attributes, whereas behavioral, metabolic and immunological variables did not exhibit any strong interrelations with the other variables. The previously well documented strong correlation of production traits with distinct prokaryote groups was confirmed. This mainly included a negative correlation of OTU ascribed to the Prevotella genus with milk and fat yield and feed
efficiency. A central role of the animals´ feed intake behavior in this context could however not be affirmed. Furthermore, different methodological and interpretability aspects concerning the microbiome analysis by 16S rRNA gene sequencing, such as the discrepancy between taxonomic classification and functional communality, as well as the comparability with other studies, are discussed. It is concluded that, to further investigate the driving force that causes the difference between efficient and inefficient animals, studies including more sophisticated methods to describe phenotypical traits of the host (e.g. rumen physiology, metabolic and genetic aspects) as well as the rumen microbiome (e.g. Metagenome, transcriptome, -proteome and Metabolome analysis) are needed.
Key words: rumen microbiota, feed efficiency, behavior, immunology
Introduction
The cow, being a ruminant, lives in a symbiotic relationship with her rumen microbiota. By ingesting feed, she delivers new substrates to the microorganisms, which in return produce valuable nutrients through fermentation and form a nutrient source themselves (Mizrahi, 2013).
Early studies showed that the rumen microbial composition is very much determined by the feed composition, as well as feed intake pattern of the host (Bryant and Burkey, 1953; Warner, 1962; Mackie et al., 1978; Leedle et al., 1982). It was further described that especially the rumen protozoal population exhibits a host individuality (Kofoid and MacLennan, 1933; Eadie, 1962). With the advent of DNA fingerprinting and sequencing techniques this finding was extended to the rumen prokaryotes (Li et al., 2009; Kong et al., 2010; Welkie et al., 2010; Jami and Mizrahi, 2012). The current understanding of rumen microbial dynamics is that the rumen microbiota consists of a core and a variable microbiota (Wu et al., 2012; Creevey et al., 2014;
Henderson et al., 2015). The core microbiota is found across a wide geographical range and consists of different taxa that increase or decrease in their abundance according to the diet fed (Henderson et al., 2015). It therefore constitutes a key element in the survival strategy of ruminants, by allowing fast and appropriate adaptation to new diets (redundancy and resilience: Weimer (2015), Solden et al. (2016b), Dieho et al. (2017a), Schären et al. (2017b)).
It is thought that the variable or individual microbiota is a result of inter-animal variation in behavioral and genetic attributes, as well as environmental influences (Mizrahi, 2013; McCann et al., 2014a; Henderson et al., 2015; Weimer, 2015; Malmuthuge and Guan, 2017). Different studies have shown interrelations between production variables, such as feed efficiency (Guan et al., 2008; Zhou et al., 2009; Hernandez-Sanabria et al., 2010; Zhou et al., 2010; Carberry et al., 2012; Hernandez-Sanabria et al., 2012; Rius et al., 2012; McCann et al., 2014b; Myer et al., 2015; Jewell et al., 2015; Shabat et al., 2016; Li and Guan, 2017) and milk production
and composition (Jami et al., 2014; Lima et al., 2015), and the rumen microbiota. However, the underlying dynamics concerning cause and effect still need to be elucidated (Weimer, 2015; Malmuthuge and Guan, 2017).
Therefore, the aim of the current study was to investigate the associations between the rumen microbiota and a large set of variables describing the production, as well as the metabolic and immunological state of dairy cows in early lactation, plus behavioral attributes, attempting to describe possible functional interrelations and pathways.
Material and Methods
Experimental work was conducted at the experimental station of the Institute of Animal Nutrition (Friedrich-Loeffler-Institute) in Brunswick, Germany. The experiment was carried out in accordance with the German Animal Welfare Act approved by the LAVES (Lower Saxony State Office for Consumer Protection and Food Safety, Germany).
Experimental Design
The data was collected in a trial investigating the influence of monensin and a blend of essential oils on production, rumen fermentation and metabolic variables (Drong et al., 2016a), immunology (Drong et al., 2016b), and rumen microbiome (Schären et al., 2017a) of transition dairy cows. Sixty German Holstein cows were divided into a low (n=15) and a high condition group (n=45) at the beginning of the dry period according to their body condition score (BCS).
The animals in the high condition group were further divided into a control group and two treatment groups (either receiving a blend of essential oils or a monensin controlled release bolus, n=15). The animals of the low condition group were then fed a normal transition ration (80 % roughage (50 % maize silage, 50 % grass silage) and 20 % concentrate based on DM content) during the dry period, and after calving a TMR with an initial concentrate feed proportion of 30 %, which was increased stepwise to 50 % of the daily ration within 2 weeks (details in Drong et al. (2016a)). The animals in the high condition group were exposed to a ketogenic ration by an oversupply with energy during the dry period (concentrate feed proportion of 60 %) and a subsequently decelerated increase in concentrate feed proportion p.p. (from 30 % to 50 % in 3 instead of 2 weeks, animal model described in Schulz et al.
(2014)). Production data, blood, liver and rumen fermentation samples were collected troughout the trial at different points in time and showed that monensin increased the energy availability in the animal by increasing the ruminal propionate production, whereas the blend of essential oils failed to elicit any positive effect (Drong et al., 2016a; b). At day 56 postpartum (p.p.) oral rumen liquid samples were collected from 48 animals (n=12) to investigate the underlying microbial alterations. The PCR-single-strand-conformation-polymorphism (SSCP) and 16S rRNA gene amplicon sequencing analysis revealed alterations in the rumen
microbiota due to the monensin treatment, and corresponding with the results on animal level no effect of the blend of essential oils was observed (Schären et al., 2017a). At that stage of lactation (day 56 p.p.) no significant difference between the control and the essential oil groups on production, metabolic and rumen fermentation level was observed (except higher serum protein concentrations in the high condition control group), whereas the monensin treated animals exhibited lower serum BHB as well as higher rumen propionate concentrations (Schären et al., 2017a). Because the present study aimed to investigate the interrelations between phenotypic traits of dairy cows and their rumen microbiota at a normal physiological state, it was decided to exclude the animals of the monensin group in the current analysis. The animals in the two control and the essential oils group (n = 36) were considered as being representative for dairy cows at this stage of lactation with an average milk production of 30.7
± 6.0 kg/d (fat %: 4.3 ± 0.7, protein %: 2.8 ± 0.2, means ± SD), DMI of 17.7 ± 2.9 kg/d, a BCS of 2.9 ± 0.4, and all measured metabolic variables in physiological ranges (e.g. serum BHB of 0.86 ± 0.36 mmol/L and fatty acids of 0.48 ± 0.29 mmol/L) (Schären et al., 2017a).
Sample Collection and Analysis
All variables presented in the manuscript were assessed in samples collected at day 56 p.p., a detailed description of sample collection and analysis has been published in Drong et al.
(2016a; 2016b) and Schären et al. (2017a). The production data (milk yield, milk components and body weight (BW)) were summarized for the period day 56 ± 3 d p.p. and have previously been published in Schären et al. (2017a).
Production and behavior variables. The cows were milked twice per day (0530 h and 1530 h) and the BW was recorded thereafter. The TMR was fed ad libidum (offered fresh daily at approximately 1100 h) and individual intakes were continuously recorded using electronic balance troughs (Insentec, B.V., Marknesse, The Netherlands). Additionally, a small amount of concentrate (2.18 ± 0.73 kg DM/cow/d, depending on the individual daily TMR intake to adjust to a concentrate feed proportion of 50 %) was fed using automated feeding stations (Insentec, B.V., Markenesse, The Netherlands). The cows were able to access the first half of their daily ration at midnight, the second half at 0800 h. The concentrate was fed in 100 g portions (a new portion was fed whenever the cow remained in the station, until ration limit was reached). The number and sizes of the ingested portions were recorded for TMR [number of meals per day (accessions when no feed was ingested excluded) = TMR Freq, average TMR intake per meal = TMR Size, variation in meal size = SD of TMR Size = TMR SD, average duration of a meal (time eating) = TMR TE, number of times cow accessed a weighing trough without eating (zero counts) = TMR ZC, average time spent at weighing trough when no feed was ingested (time zero) = TMR TZ], concentrate intake [CI, total CI from feeding station per day = CI Tot, number of meals per day at feeding station (accessions when no concentrate
was ingested excluded) = CI Freq, average CI per meal = CI Size, variation in meal size = SD of CI Size = CI SD, average duration of meals = CI TE, how often the cow was standing in the feeding station without concentrate allowance/not eating = CI ZC, average time spent in feeding station when no feed was ingested = CI TZ] and water intake [WI, times cow accessed water trough to drink (accessions when no water was ingested excluded) = WI Freq, average WI per drinking = WI Size, variation in amount of water ingested = SD of TMR Size = WI SD, average duration of drinking = WI TD, number of times cow accessed a water trough without drinking = WI ZC, average time spent at water trough when no water was ingested = WI TZ].
The BCS was assessed according to a 5-point scale (Edmonson et al., 1989). Milk samples were collected twice per week (weighed means of Monday evening & Tuesday morning and Thursday evening & Friday morning, stored at 4 °C until analysis). Milk component analysis included fat, protein, lactose and urea concentrations using an infrared milk analyzer (Milkoscan FT 6000, Foss Electric A/S, Hillerød, Denmark) and weighted daily means were calculated.
From the different production variables three feed efficiency estimators were calculated. For the first variable, the daily DMI was divided by the daily milk yield (FE1 = inverse of gross feed efficiency = DMI/milk yield). For the second variable, the daily protein intake was divided by the daily milk protein production (FE2 = protein intake/protein yield). The two variables exhibited a mean of 0.56 and 2.74 with an SD of 0.09 and 0.42, respectively. Also, the residual energy intake (REI) was calculated for each animal according to Hurley et al. (2016).
Therefore, a piecewise regression construction was applied including different possible energy sinks to estimate the estimated net energy intake (NEIntakeEstimated), resulting in the following equation (r2 = 0.614):
NEIntakeEstimated = -85.24+(1.45*MilkKG)+(38.40*MilkProtein%)+(0.16*BW)+(-12.23*BCS) Thereafter, the actual net energy intake (NEIntakeActual) was calculated from energy content of the diet multiplied by the individual DMI. The individual REI was then calculated by subtraction of the two latter variables: REI = NEIntakeEstimated - NEIntakeActual. The REI exhibited a mean of -3.06 with a SD of 13.04.
Blood variables. Blood samples were collected from a Vena jugularis externa in a serum separating, an EDTA, and a heparin containing blood tube.
A complete blood count of each sample was performed within 2 h after sampling using an automated hematology analyzer [Celltac alpha MEK-6450, Nihon Kohden Corporation, Tokyo, Japan; including red blood cell count (RBC) and associated variables (MCV, MCH, MCHC), hemoglobin (HGB), hematocrit (HCT), white blood cell population (lymphocyte (LY), monocyte
(MO), eosinophile (EO) and granulocyte (GR) count, and the respective proportions (LY%, MO% and GR%)), Schären et al. (2016a)].
The serum was separated and stored at -80 °C before chemical analysis using an automatic clinical chemistry analyzer (Eurolyser CCA180, Eurolyser Diagnostica GmbH, Salzburg, Austria; including serum glucose, beta-hydroxy-butyrate (BHB), fatty acids (FA), urea, albumin, total protein, cholesterol, aspartate transaminase (AST), γ-glutamyltransferase (γ-GT), total bilirubin, glutamate dehydrogenase (GLDH) and triglyceride, Schären et al. (2016a)).
Kynurenine (Kyr) and tryptophan (Trp) concentrations were determined via HPLC (Shimadzu, Kyoto, Japan; as described in detail in Drong et al. (2017)) and their ratio (Kyr:Trp) calculated.
Peripheral blood mononuclear cells (PBMC) were isolated from whole blood samples (heparinized blood) by gradient centrifugation and stored at -80 °C until analysis (Renner et al., 2011; Drong et al., 2016b). Cell metabolic activity and Concavalin A (ConA)-stimulated proliferation of PBMC were evaluated using the Alamar Blue (AB) assay (described in detail in Drong et al. (2016b)). The proliferation of PBMC (stimulation index ex vivo) was then calculated by the ratio between fluorescence of ConA-stimulated PBMC (ABstim) and non-stimulated PBMC (ABunstim): SI = (Fluorescence of ConA - Fluorescence of ABstim)/(Fluorescence of ABunstim) (Drong et al., 2016b). T-cell phenotyping was performed on whole blood samples by monoclonal antibodies staining (for CD4 and CD8, or isotype controls) and subsequent flow cytometry analysis (FACSCantoII, BD Bioscience, San Jose, CA, USA).
The lymphoid populations were then identified according to their side- and forward-scattering properties and an estimated number of T-cells of each phenotype as well as their ratio (CD4+/CD8+) were calculated using the percentages obtained by the flow cytometer (Drong et al., 2016b). The capacity of polymorphonuclear leukocytes (PMN) and leukocytes to release reactive oxygen species (ROS) was measured using an assay based on the oxidation of the nonfluorescent dihydrorhodamine 123 (DHR) to the fluorescent metabolite rhodamine 123 (R123) by quantifying the mean fluorescence intensity (MFI,flow cytometry (FACSCantoII)), resulting in the basal amount of PMN and lymphocytes that oxidase DHR (R123%unstim and R123+ Lym%) and the amount of radicals that are produced per cell on average (R123 MFIunstim and R123+ Lym MFI). To elicit and quantify the maximum oxidative burst capacity the PMN were additionally stimulated in parallel with phorbol-12-mystristate-13-acetate (PMA), resulting in the population of PMN performing an oxidative burst (R123%stim) and the mean capacity per cell for the oxidative burst (R123 MFIstim) (Drong et al., 2016b).
Rumen fermentation variables. Rumen liquid samples (ca. 750 mL) were collected in the morning after milking (prior to feeding) orally using an oral rumen tube and a hand vacuum pump. Immediately after collection, pH was measured using a glass electrode (model: pH 525;
WTW, Weilheim, Germany). Until further processing (approximately 1-2 h after sample
collection), samples for ammonia (NH3-N), volatile fatty acids (VFA, total VFA concentration = VFA, acetate (C2), propionate (C3), butyrate (C4), isovalerate (iC5), valerate (C5), and their respective proportions (C2%, C3%, C4%, iC5%, C5%), and samples for LPS concentration were cooled to 4 °C. Volatile fatty acids were determined according to Koch et al. (2006) using a gas chromatograph (Gaschromatograph 5890 II, Hewlett Packard®, Böblingen, Germany) and NH3-N was determined using steam distillation according to the Kjeldahl method (DIN38406-E5-2, Anonymous (1998)). Rumen lipopolysaccharide (LPS) concentrations were measured spectrophotometrically after centrifugation, heating and storage (-20 °C), using the Limulus amebocyte lysate (LAL) assay (Kinetic-QCLTM, Lonza, Walkersville, MD, USA;
following the manufacturer´s instructions) as described in Schären et al. (2016b) and Gozho et al. (2005)). Due to technical issues at the day of sample collection and storage, 11 rumen liquid samples for fermentation variable analysis were lost. The dataset therefore includes fermentation variables of n = 25 instead of the total 36 animals.
Rumen microbiome analysis. Rumen liquid samples for protozoal density assessment were fixed 1:1 with methylgreen-formalin solution (stored at 4 °C) directly after sample collection and protozoa were counted using a Fuchs-Rosenthal chamber under an optical microscope and differentiated into entodiniomorpha and holotrichia (Ogimoto and Imai, 1981). The rumen liquid samples for microbial analysis were stored at -20 °C and DNA extraction was performed as described in Schären et al. (2017a). Prokaryotes were separated from feed and debris trough several centrifugation steps (1 x 5 min at 600 g, 4 x 20 min at 27’000 g (4 °C) with resuspension between centrifugation steps in 0.9 % NaCl) and liquid shock frozen as little droplets. DNA extraction involved mechanical lysis by a bead beating method, incubation steps with lysozyme, sodiumdodecylsulfate, proteinase K and cetyltrimethylammoniumbromide, a protein purification step with (phenol)-chloroform-isoamylalcohol, and washing steps using the peqGold Tissue-Kit (peq lab, Erlangen, Germany). Samples were stored at 4 °C until further analysis. For sequencing, gDNA samples were sent to Microsynth AG (Balgach, Switzerland).
For 16S sequencing library preparation the primers A519F (S-D-Arch-0519-a-S-15) and 802R (S-D-Bact-0785-b-A-18) were chosen (Klindworth et al., 2013) and amplified using the HiFi HotStart PCR Kit (Kapa Biosystems, Wilmington, MA, USA). Illumina Nextera Libraries were prepared according to the manufacturer’s instruction (Illumina, San Diego, USA) and sequencing was performed on the Illumina MiSeq Sequencing System using the Illumina MiSeq reagent Kit v2 (2 x 250bp). De-multiplexing (using the Illumina MiSeq v2.5.1.3. reporter and cutadapt v1.8.1 software package (Martin, 2011)), read stitching (using FLASH v1.2.11 (Magoč and Salzberg, 2011)), de novo Chimera detection, identification and removal (using the Uchime v4.2 (Edgar et al., 2011) and Usearch v8.1.1861 (Edgar, 2010)), operational taxonomic units (OTU) clustering (based on 97 % sequence similarity, using Uclust (Edgar, 2010)), and taxonomic assignment (using QIIME v1.9.1 (Caporaso et al., 2010) and the SILVA
rRNA database v111 (Quast et al., 2013)) was performed. Singeltons were removed from the dataset to reduce bias introduced by sequencing errors. For downstream analysis, only OTU with a relative abundance of at least 0.1% were considered to guarantee a solid differentiation between artifact and true organism. The amplicon sequencing resulted in 12’206 ± 3’424 reads (after filtering, mean ± SD) per sample. In total, a number of 177 different species-level OTUs were identified (167 ± 6 (mean ± SD) different OTUs per sample). A list of the detected OTU, their corresponding taxonomic classification and proportional abundance is given in Table 1.
Most OTUs could be taxonomically classified to the family level, while their genus or species level affiliation were “uncultured bacterium or archaeon” in many cases.
Statistical Analysis
All statistical analyses were performed in STATISTICA 12.0 (StatSoft, Inc. 2014, Tulsa, OK, USA). If variables were recorded more than once a week, means were calculated per cow and week (day 56 ± 3 p.p.) prior to statistical evaluation. To obtain a normal distribution, rumen LPS concentrations were logarithmically transformed prior to statistical analysis. Correlation coefficients between different variables were estimated using Spearman´s Rank correlation and results were considered significant at P < 0.05. Principal component analysis (PCA) was performed including all measured variables. Thereafter, for some variables (which clustered due to their close biological interrelation) only one representative variable was chosen, to obtain a clearly arranged graph. This was the case for following variables (representative variable chosen, followed by deleted variables in brackets): DMI (Starch Intake, Concentrate Intake, Protein Intake, NEIntakeActual), HCT (HGB, MCV, MCH, MCHC), and LY, MO, GR (LY%, MO%, GR%, respectively).
Results and Discussion
All data points were summarized in one “master-table” and a Spearman´s Rank correlation with all 265 variables was conducted, resulting in a table with 70,225 correlations (Appendix 10, available for download in supplements). Because the interpretation and the deduction of biological dynamics from such a table is difficult, one main PCA plot (“master-plot”) was created, to gain a better overview of the data and the different interrelations (Figure 1, interactive document available for download in supplements). The “master-plot” is interpreted as follows: each dot represents one variable (grouped by color). The closer two variables are, the higher they are positively associated. The farther away two variables are, the more negatively they are associated. The more variables are located towards the middle of the plot, the less they contribute to the variation within the data and the smaller is their association with
All data points were summarized in one “master-table” and a Spearman´s Rank correlation with all 265 variables was conducted, resulting in a table with 70,225 correlations (Appendix 10, available for download in supplements). Because the interpretation and the deduction of biological dynamics from such a table is difficult, one main PCA plot (“master-plot”) was created, to gain a better overview of the data and the different interrelations (Figure 1, interactive document available for download in supplements). The “master-plot” is interpreted as follows: each dot represents one variable (grouped by color). The closer two variables are, the higher they are positively associated. The farther away two variables are, the more negatively they are associated. The more variables are located towards the middle of the plot, the less they contribute to the variation within the data and the smaller is their association with