Advertisement

Heterogeneity in genetic and nongenetic variation and energy sink relationships for residual feed intake across research stations and countries

Open ArchivePublished:January 09, 2015DOI:https://doi.org/10.3168/jds.2014.8510

      Abstract

      Our long-term objective is to develop breeding strategies for improving feed efficiency in dairy cattle. In this study, phenotypic data were pooled across multiple research stations to facilitate investigation of the genetic and nongenetic components of feed efficiency in Holstein cattle. Specifically, the heritability of residual feed intake (RFI) was estimated and heterogeneous relationships between RFI and traits relating to energy utilization were characterized across research stations. Milk, fat, protein, and lactose production converted to megacalories (milk energy; MilkE), dry matter intakes (DMI), and body weights (BW) were collected on 6,824 lactations from 4,893 Holstein cows from research stations in Scotland, the Netherlands, and the United States. Weekly DMI, recorded between 50 to 200 d in milk, was fitted as a linear function of MilkE, BW0.75, and change in BW (ΔBW), along with parity, a fifth-order polynomial on days in milk (DIM), and the interaction between this polynomial and parity in a first-stage model. The residuals from this analysis were considered to be a phenotypic measure of RFI. Estimated partial regression coefficients of DMI on MilkE and on BW0.75 ranged from 0.29 to 0.47 kg/Mcal for MilkE across research stations, whereas estimated partial regression coefficients on BW0.75 ranged from 0.06 to 0.16kg/kg0.75. Estimated partial regression coefficients on ΔBW ranged from 0.06 to 0.39 across stations. Heritabilities for country-specific RFI were based on fitting second-stage random regression models and ranged from 0.06 to 0.24 depending on DIM. The overall heritability estimate across all research stations and all DIM was 0.15 ± 0.02, whereas an alternative analysis based on combining the first- and second-stage model as 1 model led to an overall heritability estimate of 0.18 ± 0.02. Hence future genomic selection programs on feed efficiency appear to be promising; nevertheless, care should be taken to allow for potentially heterogeneous variance components and partial relationships between DMI and other energy sink traits across environments when determining RFI.

      Key words

      Introduction

      Lifetime feed efficiency of dairy cows has increased over the past 60 yr, as cows have produced more milk and thus decreased the portion of feed needed for maintenance. Elite Holsteins currently produce at more than 4 times their maintenance requirement, so the advantage of higher production at the population level to improve biological efficiency in dairy cattle genetics due to the dilution of maintenance costs has been mostly exploited (
      • VandeHaar M.J.
      • St-Pierre N.
      Major advances in nutrition: Relevance to the sustainability of the dairy industry.
      ). Further increases in feed efficiency must focus on selecting cows directly for their ability to convert feed to milk production. Feed intake is expensive to collect relative to other economically important traits (e.g., milk yield or reproduction variables); hence, extensive collaboration among researchers is required (
      • Banos G.
      • Coffey M.P.
      • Veerkamp R.F.
      • Berry D.P.
      • Wall E.
      Merging and characterising phenotypic data on conventional and rare traits from dairy cattle experimental resources in three countries.
      ). Estimates of genetic parameters (i.e., heritabilities) require large enough sample sizes to facilitate reliable inferences (
      • Klein T.W.
      Heritability and genetic correlation: Statistical power, population comparisons, and sample size.
      ), and probably much larger than the few hundred cows that have characterized many previous studies on feed efficiency. Finally, a robust scope of inference necessitates the contribution of data from a broad range of management systems and environments. Understanding the relationships among feed intake and energy outputs in cows across environments is critical to the possible development of successful selection strategies for feed efficiency.
      Residual feed intake (RFI) has been proposed as a measure of dairy feed efficiency as an alternative to ratio-based (i.e., input/output) measures, partly because of its ability to take into account changes in energy dynamics over the course of a lactation due to BW changes, for example (
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      ). The determination of RFI is typically a 2-step modeling process. The first stage model is an energy sink model (
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      ) whereby feed intake, whether recorded as DMI (kg) or energy intake (Mcal), is typically specified as partial linear regression functions of key energy sinks such as milk energy (MilkE), metabolic BW (MBW) defined as BW0.75, BW changes (ΔBW), and cohort effects. The resulting estimated residuals from this first stage model are deemed to be RFI phenotypes. These RFI phenotypes are then used as the response variables in a second stage quantitative genetics model to estimate heritabilities and, subsequently, breeding values for feed efficiency.
      Heritability estimates for RFI in lactating dairy cattle range widely, from 0.0 to 0.38 in Holsteins, as reported in a review by
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      ; these estimates are typically smaller than heritability estimates for RFI in growing heifers (
      • Williams Y.J.
      • Pryce J.E.
      • Grainger C.
      • Wales W.J.
      • Linden N.
      • Porker M.
      • Hayes B.J.
      Variation in residual feed intake in Holstein-Friesian dairy heifers in southern Australia.
      ;
      • Pryce J.E.
      • Arias J.
      • Bowman P.J.
      • Davis S.R.
      • Macdonald K.A.
      • Waghorn G.C.
      • Wales W.J.
      • Williams Y.J.
      • Spelman R.J.
      • Hayes B.J.
      Accuracy of genomic predictions of residual feed intake and 250-day body weight in growing heifers using 625,000 single nucleotide polymorphism markers.
      ;
      • Gonzalez-Recio O.
      • Pryce J.E.
      • Haile-Mariam M.
      • Hayes B.J.
      Incorporating heifer feed efficiency in the Australian selection index using genomic selection.
      ), as further noted by
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      . Each lactating RFI study cited in that review was based on limited cow numbers, ranging from 204 (
      • Veerkamp R.F.
      • Emmans G.C.
      • Cromie A.R.
      • Simm G.
      Variance components for residual feed intake in dairy cows.
      ) to 970 (
      • Vallimont J.E.
      • Dechow C.D.
      • Daubert J.M.
      • Dekleva M.W.
      • Blum J.W.
      • Barlieb C.M.
      • Liu W.
      • Varga G.A.
      • Heinrichs A.J.
      • Baumrucker C.R.
      Heritability of gross feed efficiency and associations with yield, intake, residual intake, body weight, and body condition score in 11 commercial Pennsylvania tie stalls.
      ). Furthermore, some studies were limited by infrequent (i.e., monthly) recording of DMI (
      • Vallimont J.E.
      • Dechow C.D.
      • Daubert J.M.
      • Dekleva M.W.
      • Blum J.W.
      • Barlieb C.M.
      • Liu W.
      • Varga G.A.
      • Heinrichs A.J.
      • Baumrucker C.R.
      Heritability of gross feed efficiency and associations with yield, intake, residual intake, body weight, and body condition score in 11 commercial Pennsylvania tie stalls.
      ). A recent study (
      • Berry D.P.
      • Coffey M.P.
      • Pryce J.E.
      • de Haas Y.
      • Lovendahl P.
      • Krattenmacher N.
      • Crowley J.J.
      • Wang Z.
      • Spurlock D.
      • Weigel K.
      • Macdonald K.
      • Veerkamp R.F.
      International genetic evaluations for feed intake in dairy cattle through the collation of data from multiple sources.
      ) involving nearly 7,000 lactating cows determined an average estimated heritability of 0.34 for DMI; however, heritability estimates for DMI tend to be larger than that for RFI (
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      ;
      • Gonzalez-Recio O.
      • Pryce J.E.
      • Haile-Mariam M.
      • Hayes B.J.
      Incorporating heifer feed efficiency in the Australian selection index using genomic selection.
      ), likely because DMI is strongly genetically correlated with MilkE and BW. Finally, in contrast to some recent work involving DMI (
      • Manzanilla Pech C.I.V.
      • Veerkamp R.F.
      • Calus M.P.L.
      • Zom R.
      • van Knegsel A.
      • Pryce J.E.
      • de Haas Y.
      Genetic parameters across lactation for feed intake, fat- and protein-corrected milk, and liveweight in first-parity Holstein cattle.
      ), most quantitative genetic analyses of dairy feed efficiency have not modeled potential DIM-specific genetic effects nor the gradually changing temporal effects that can be inferred upon using random regression model analyses, now considered to be the global standard for genetic evaluation of milk production traits (
      • Schaeffer L.R.
      Application of random regression models in animal breeding.
      ).
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      further cautioned that in some cases RFI may not represent true feed efficiency but can be due, in part, to biases in estimated regression coefficients in the first-stage model resulting from measurement or prediction errors in the covariates. This can consist of several different scenarios, including feed intake not being recorded daily, and different recording frequencies or different types of recording equipment for various traits. Furthermore, the partial regression relationships between DMI and MilkE, MBW, and ΔBW can be highly heterogeneous across different trials (
      • Davis S.R.
      • Macdonald K.A.
      • Waghorn G.C.
      • Spelman R.J.
      Residual feed intake of lactating Holstein-Friesian cows predicted from high-density genotypes and phenotyping of growing heifers.
      ). Finally, heterogeneous variances might exist across regions or herds such that, without accounting for it, the estimated genetic merit of RFI could be inflated in high-variance regions but be unduly muted in low-variance herds (
      • Hill W.G.
      On selection among groups with heterogeneous variance.
      ).
      Our long-term goal is to develop breeding strategies for dairy feed efficiency. As a step toward that goal, we pooled data from several research stations across 3 countries to address several key research questions pertinent to a quantitative genetic analysis of RFI. Because of its implications for deriving RFI phenotypes, our first objective was to characterize any potential heterogeneity in the partial regression coefficients of DMI on various energy sinks (i.e., MilkE, MBW, ΔBW) across research stations or studies. Our second objective was to determine how genetic and permanent environmental variability in RFI might vary across stage of lactation, separately by country, and jointly together.

      Materials and Methods

      Data

      Milk yields (MY), DMI, BW, fat (FAT%), protein (PROT%), and lactose (LACT%) components, were collected from Holstein cows on research stations within 3 countries: the United States (US), the Netherlands (NL), and the United Kingdom (UK). Data from US were derived from 6 research stations: Iowa State University (ISU; Ames), Michigan State University (MSU; East Lansing), the University of Florida (UF; Gainesville), the University of Wisconsin-Madison (UW), the United States Dairy Forages Research Center (FRC; Madison, WI), and the USDA Animal Genomics and Improvement Laboratory (AGIL; Beltsville, MD). Some of the data provided by ISU has been analyzed previously and described in more detail by
      • Spurlock D.M.
      • Dekkers J.C.
      • Fernando R.
      • Koltes D.A.
      • Wolc A.
      Genetic parameters for energy balance, feed efficiency, and related traits in Holstein cattle.
      and
      • Yao C.
      • Spurlock D.M.
      • Armentano L.E.
      • Page C.D.
      • VandeHaar M.J.
      • Bickhart D.M.
      • Weigel K.A.
      Random Forests approach for identifying additive and epistatic single nucleotide polymorphisms associated with residual feed intake in dairy cattle.
      ), whereas some of the UW data has been analyzed previously as well (
      • Ferraretto L.F.
      • Shaver R.D.
      • Espineira M.
      • Gencoglu H.
      • Bertics S.J.
      Influence of a reduced-starch diet with or without exogenous amylase on lactation performance by dairy cows.
      ;
      • Ferraretto L.F.
      • Shaver R.D.
      • Bertics S.J.
      Effect of dietary supplementation with live-cell yeast at two dosages on lactation performance, ruminal fermentation and total-tract nutrient digestibility in dairy cows.
      ;
      • He M.
      • Perfield K.L.
      • Green H.B.
      • Armentano L.E.
      Effect of dietary fat blend enriched in oleic or linoleic acid and monensin supplementation on dairy cattle performance, milk fatty acid profiles, and milk fat depression.
      ). Data from AGIL were previously analyzed by
      • Connor E.E.
      • Hutchison J.L.
      • Norman H.D.
      • Olson K.M.
      • Van Tassell C.P.
      • Leith J.M.
      • Baldwin R.L.
      Use of residual feed intake in Holsteins during early lactation shows potential to improve feed efficiency through genetic selection.
      .
      Data from NL were derived from 4 primary experimental herds and studies: (1) an experimental herd ‘t Gen (TGEN) in Lelystad previously described by
      • Veerkamp R.F.
      • Oldenbroek J.K.
      • Van Der Gaast H.J.
      • van der Werf J.H.J.
      Genetic correlation between days until start of luteal activity and milk yield, energy balance, and live weights.
      and
      • Banos G.
      • Coffey M.P.
      • Veerkamp R.F.
      • Berry D.P.
      • Wall E.
      Merging and characterising phenotypic data on conventional and rare traits from dairy cattle experimental resources in three countries.
      , (2) the Nij Bosma Zathe (NBZ) herd located near Leeuwarden and also previously described by
      • Banos G.
      • Coffey M.P.
      • Veerkamp R.F.
      • Berry D.P.
      • Wall E.
      Merging and characterising phenotypic data on conventional and rare traits from dairy cattle experimental resources in three countries.
      , (3) a third study (ZOM) based on the work by
      • Zom R.L.G.
      • André G.
      • van Vuuren A.M.
      Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
      , and (4) a compilation of studies labeled NLN based on data collected from various nutritional experiments. Data on all variables (i.e., DMI, MY, BW, FAT%, PROT%, and LACT%) from each of these 4 NL studies were summarized on a weekly basis before further analysis.
      Data from the UK were originally derived from the Langhill (LAN) farm near Edinburgh from 1992 to 2001 and from the Scottish Agricultural College (SAC) Dairy Research Centre based at Crichton Royal Farm near Dumfries with data collection from 2003 to 2011. These data are described in more detail by
      • Banos G.
      • Coffey M.P.
      • Veerkamp R.F.
      • Berry D.P.
      • Wall E.
      Merging and characterising phenotypic data on conventional and rare traits from dairy cattle experimental resources in three countries.
      and various citations therein. These data, and the NBZ herd from NL, were the basis for selection experiments, each involving 2 different selection lines for fat plus protein yield.
      The experimental designs from which the data originated varied from single ration studies completely dedicated to investigating genetic aspects of feed efficiency to completely randomized designs (i.e., cows never assigned to more than 1 ration), simple crossover designs, and multiple treatment Latin square designs. In total, 185 different rations were identified across the 12 different research stations.

      Phenotypes

      Only data between 50 and 200 DIM were used for this analysis, as ΔBW was conjectured to be relatively more stable within this window than outside it such that errors in determining RFI due to large ΔBW would be minimal. Data from a particular lactation were used only if at least 28 daily records on each of DMI and MY were available within that lactation window. Because the NL data were already provided on a weekly basis, this edit was modified to 4 weekly records on each of DMI and MY between 50 and 200 DIM. There were periods of time where SAC cows were recorded to be on pasture; these data were deleted before invoking the aforementioned edits on that particular data source. Only cows ≥75% Holstein were included in the current study.
      Recording frequencies for all other traits (e.g., BW, FAT%, PROT%, LACT%) varied widely among stations and even by different experiments within stations. Frequency of BW recording ranged from daily (SAC, FL, MSU, UW) to weekly (ISU, MSU, FRC, LAN, TGEN, NBZ, ZOM, NLN) or biweekly (AGIL, UW), with some hybrid protocols such as BW being recorded on 2 consecutive days every 1 or 2 wk. Milk components were either recorded daily (FL), weekly (LAN, SAC, AGIL, ISU, MSU, UW, FRC, TGEN, NBZ, ZOM, NLN), or biweekly (UW). A summary of the number of cows and lactations by station, after the subsequently described edits, are provided in Table 1. Some stations (UF, FRC, TGEN, and NBZ) had only 1 lactation per cow whereas other stations (ZOM, LAN, and SAC) averaged nearly 2 or more. The LACT% data were not available from AGIL, LAN, or SAC.
      Table 1Frequency of number of cows and lactations by station
      CountryStation/study
      UF=University of Florida, ISU=Iowa State University, MSU=Michigan State University, FRC=USDA Forage Research Center, UW=University of Wisconsin-Madison, AGIL=USDA Animal Genomics Improvement Laboratory, TGEN=‘t Gen, NBZ=Nij Bosma Zaathe, ZOM=data from Zom et al. (2012), NLN=miscellaneous compilation of studies, LAN=Langhill, SAC=Scottish Agricultural College.
      No. of

      cows
      No. of

      lactations
      No. of

      weekly records
      Average number

      of weekly records

      per lactation
      Years of data

      collection
      United StatesUF2052201,6637.62009–2013
      ISU6927298,30611.42008–2013
      MSU1631921,8079.42011–2013
      FRC1151171,28010.92009–2011
      UW4124635,37211.62007–2013
      AGIL2884483,6178.12007–2011
      NetherlandsTGEN67167111,76017.51991–1998
      NBZ1001005255.32003–2004
      ZOM7831,23210,4018.41991–2001
      NLN4125315,79510.82003–2012
      United KingdomLAN5901,21724,24219.91990–2001
      SAC46290415,78617.42003–2011
      TOTAL4,8936,82490,55413.3
      1 UF = University of Florida, ISU = Iowa State University, MSU = Michigan State University, FRC = USDA Forage Research Center, UW = University of Wisconsin-Madison, AGIL = USDA Animal Genomics Improvement Laboratory, TGEN = ‘t Gen, NBZ = Nij Bosma Zaathe, ZOM = data from
      • Zom R.L.G.
      • André G.
      • van Vuuren A.M.
      Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
      , NLN = miscellaneous compilation of studies, LAN = Langhill, SAC = Scottish Agricultural College.

      Statistical Analyses

      A separate random regression model using SAS PROC HPMIXED (SAS Institute Inc., Cary, NC) was fitted to each of the 6 variables (DMI, MY, FAT%, PROT%, LACT%, and BW) for each of the 12 different stations or study compilations (i.e., ISU, UW, FRC, MSU, AGIL, UF, LAN, SAC, TGEN, NBZ, ZOM, and NLN) unless it was unavailable, as previously noted. The statistical model used in all cases was
      y=parity×k=05DIMk+rationexpt+test week+animal×k=0nkDIMk+e,
      [1]


      where y is any of the 6 variables; parity×k=05DIMk are fixed effects of parity (primiparous vs. multiparous) and specific fifth-order Legendre polynomial regressions of y on DIM, whereas ration(expt) are the fixed effects of ration within experiment; test week is the random effect of the midweek calendar date, whereas animal×k=0nkDIMk is the random effects Legendre polynomial of order nk ≤ 3 of y on DIM for the animal effect. Finally, e is a normally distributed residual. Parity was not fitted for TGEN data, as all TGEN cows were primiparous. The variance of each individual animal×k=0nkDIMk term was an nk × nk unstructured covariance matrix with independence assumed across animals. For 18 of the 69 (12 stations × 6 variables − 3 stations missing LACT%) analyses, it was required to specify nk < 3 in Equation [1], as convergence was not otherwise attained; for the other 51 analyses, nk = 3. Upon all such analyses, records of each of the 6 key variables were designated to be residual outliers if the absolute values of the corresponding externally studentized residuals were greater than 5, with 1 exception: if the absolute studentized residual exceeded 5 for DMI and the absolute studentized residual for MY exceeded 3 with the same sign and vice versa (i.e., exceeding 5 for MY and exceeding 3 for DMI, both with the same sign), then both DMI and MY records were kept as provided. This exception was intended to keep biological outliers in the data so as to best reflect energy balance between DMI and its energy sinks (i.e., as based on the other 5 variables), whereas the critical intent of data editing was to correct residual outliers (<0.5% of the data) that were most likely to be recording errors. Outlier date effects for each trait were also targeted if the corresponding date estimate exceeded 3.5 date effect standard deviations (based on the square root of the estimated variance component for date). Any records on the corresponding trait for an outlier date were corrected relative to the estimated effect of the most previous nonoutlier date to preserve any potential temporal trends typically associated with date effects.
      For all traits, observations deemed to be residual outliers were set equal to their respective BLUP based on the aforementioned random regression analyses, as were any missing daily records for each of the 6 traits between the first and last recorded DIM on all traits within any particular cow parity. Weekly averages for each trait were then created from these daily records and BLUP. To be retained for further analysis, weekly trait averages were required to include at least 1 actual daily DMI and MY record.
      MilkE was determined from weekly MY, FAT%, PROT%, and LACT% using MilkE = (0.0929 × FAT% + 0.0563 × PROT% + 0.0395 × LACT%) × MY (National Research Dairy Council, 2001). When LACT% was missing, as was true for AGIL, LAN, and SAC, a value of 4.85 was substituted as conventionally recommended in such cases (National Research Dairy Council, 2001). Weekly MBW was computed as the weekly average BW0.75. The weekly ΔBW was determined based on a simple linear regression of actual or projected BW on days from the beginning to the end of a week. The ration assignment for any particular week was based on the ration that was fed at the midpoint of the week.
      Determination of RFI typically requires a 2-step modeling process. The first model is an energy sink model (
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      ) based on weekly DMI, which was conducted separately for each of the 12 different research stations or studies. This was done in the recognition that these energy sink relationships could vary substantially across stations or studies; furthermore, variability in DMI and, hence, RFI might be rather heterogeneous as well across studies, whether because of management differences or differences in recording protocols. The first-stage model was
      DMI=parity×k=05DIMk+b1MilkE+b2MBW+b3ΔBW+ rationexpt+test week+ϵ,
      [2]


      where parity×k=05DIMk are the fixed effects of parity (primiparous vs. multiparous) and specific fifth-order Legendre polynomial regressions of DMI on DIM, b1 is the partial regression coefficient of DMI on MilkE, b2 is the partial regression coefficient of DMI on MBW, b3 is the partial regression coefficient of DMI on ΔBW, ration(expt) is the random effect of experiment-specific rations, whereas test week is the random effect of test week with ε being the residual. The estimated residuals ϵˆ from the corresponding analysis were determined to be the weekly RFI for use in a subsequent second stage quantitative genetic analysis.
      Pedigrees went at least 2 generations on each cow. However, because the data as a whole spanned 2 decades, even within some stations (e.g., LAN and SAC), the number of generations of pedigree information on many animals were generally went beyond 2 generations. Genetic groups for animals with unknown parents were created based on 5-yr birthdate increments within country separately for males and females (
      • Quaas R.L.
      Additive genetic model with groups and relationships.
      ). The second-stage (variance components) model was
      RFI=μ+cow×k=0naDIMk+cow×parity×k=0npDIMk+ cow+e,
      [3]


      where RFI is fitted as a function of overall mean μ with cow×k=0naDIMk being the order na Legendre polynomial on DMI for random additive genetic effects (with correlation based on the numerator relationship matrix with genetic groups), cow×parity×k=0npDIMk being the order np Legendre polynomial on DMI for cow-parity-specific lactations (based on the correlation matrix being the identity matrix) to model within random lactation permanent environmental effects, cow being the random intercept of cow (based on the correlation matrix being the identity matrix), and e being the residual with variance σe2. This last term (cow) was required to separate between-lactation permanent environmental effects from within-lactation permanent environmental effects when multiple lactations per cow occurred. The Legendre polynomial orders were constrained up to na = np = 3 until the next polynomial was no longer statistically significant based on simple Wald tests on the variance components, or until REML convergence was no longer attainable.
      Separate analyses for quantitative genetic variance components were run for each country to assess potential heterogeneity in genetic variability across countries; furthermore, a joint analysis involving all of the data was also run. In all cases, heterogeneous residual variances across the 12 different research stations were fitted. Overall estimates of heritabilities (i.e., across all DIM) were based on reduced model subsets (random intercepts only for genetic and both sets of permanent environmental effects) of the aforementioned random regression models.
      Finally, we also considered a 1-step joint model analysis of the data whereby the 2 previously described first- and second-stage models (Equations [2] and [3]) were combined together as one, thereby excluding ε in Equation [2]. In other words, DMI was modeled jointly as a function of all the terms in Equations [2] and [3], with the exception that only random intercept terms (i.e., na = np = 0) were fitted to animal-specific effects to infer an overall heritability estimate in each case. We also considered additional terms in this joint model:
      DMI=Equation2RHS+Equation3RHS+rationexpt×MilkEMBWΔBW+station×interceptMilkEMBWΔBW+e,
      [4]


      where Equation [2]RHS and Equation [3]RHS pertain to terms on the right hand side of Equations [2] and [3], respectively. Furthermore, ration(expt) × (MilkE MBW ΔBW) are the effects of ration and its interaction effects with MilkE, MBW, and ΔBW. These interaction effects, along with the random intercept term ration(expt) from Equation [2], were modeled as correlated random effects to more reliably infer the potential (co)variability in the energy sink partial regression coefficients across rations while jointly accounting for animal sources of variability (i.e., genetic and permanent environment). This joint model additionally specified station × (intercept MilkE MBW ΔBW) as the random effects of research stations and their interactions with MilkE, MBW, and ΔBW to allow inference on station-specific partial regression coefficients just as we implicitly conducted with the station-specific first-stage energy sink model analyses in Equation [2]. We also specified station-specific residual variances. These 1-step Equation [4] analyses were conducted separately for each country and jointly. Variance components were estimated by REML using the ASREML package (
      • Gilmour A.R.
      • Gogel B.J.
      • Cullis B.R.
      • Thompson R.
      ). The DIM-specific genetic variances σg2, within-lactation permanent variances σp2, and between-lactation permanent environmental variances σb2 were estimated using DIM-based covariance functions (
      • van der Werf J.H.J.
      • Goddard M.E.
      • Meyer K.
      The use of covariance functions and random regressions for genetic evaluation of milk production based on test day records.
      ) on the corresponding variance-covariance matrices of the random effects as listed in Equation [3], such that DIM-specific heritabilities were defined as σg2/σg2+σp2+σb2+σe2.

      Results

      Summary of Key Variables

      In total, 90,544 weekly records were derived from 6,824 lactations on 4,893 cows. The breakdown of these numbers by research station is provided in Table 1, with visual summaries of raw averages of weekly DMI, MilkE, MBW, and ΔBW provided in Figure 1. In general, DMI were lowest for LAN and NBZ, perhaps because these data were older than most of the other stations (Table 1) and because 1 diet was intended to be a low-energy diet in both herds (
      • Banos G.
      • Coffey M.P.
      • Veerkamp R.F.
      • Berry D.P.
      • Wall E.
      Merging and characterising phenotypic data on conventional and rare traits from dairy cattle experimental resources in three countries.
      ). Both DMI and MilkE were greater in SAC compared with LAN, likely because SAC cows were descendants of LAN cows. Cows from US research stations tended to have higher DMI and MilkE relative to UK research stations or older NL data (i.e., not NLN). However, substantial variability existed between stations for MBW, both within NL and the US compared with the 2 stations from the UK.
      Figure thumbnail gr1
      Figure 1Raw weekly means of (a) DMI, (b) milk energy, (c) metabolic BW, and (d) BW changes as a function of midweek DIM for USDA Animal Genomics Improvement Laboratory (AGIL; ∆), University of Florida (UF; ◊), Iowa State University (ISU; ♦), Michigan State University (MSU; □), USDA Forage Research Center (FRC; ■), University of Wisconsin-Madison (UW; ▲), ‘t Gen (TGN; +), Nij Bosma Zaathe (NBZ; ⊕),
      • Zom R.L.G.
      • André G.
      • van Vuuren A.M.
      Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
      ; ZOM; *), miscellaneous compilation of studies (NLN; ×), Langhill (LAN; ○), and Scottish Agricultural College (SAC; ●). Color version available online.

      Estimated Partial Regression Relationships from First Stage Energy Sink Model

      The first-stage or energy sink model, used to derive the RFI phenotypes, was fitted separately to each of the 12 research stations to potentially accommodate any heterogeneities in the partial regression relationships between DMI and each of the 3 key energy sinks, namely MBW, MilkE, and ΔBW. The corresponding estimates are provided in Table 2. The estimated partial regression coefficients on MilkE were somewhat heterogeneous between the research stations or studies, varying from 0.29 (NLN) to 0.47 (NBZ) kg/Mcal, and all were statistically significant (P < 0.05) from 0. Furthermore, some of these standard errors were small enough to suggest that some of the partial regression coefficients were different from each other using a simple z-test. For example, to compare any 2 stations j and j’ for b1 (partial regression on MilkE), one could readily determine the z-test statistic as
      z=bˆ1jbˆ1j'/sebˆ1j2+sebˆ2j2


      and declare a significant difference (P < 0.05) if z > 1.96. Nevertheless, these standard errors are badly understated, as discussed herein, and hence formal separation of these coefficients based on P-values is not provided in Table 2.
      Table 2Estimated partial regression coefficients (±SE) of DMI on each of several energy sinks (MilkE = milk energy, MBW = metabolic BW, ΔBW = change in BW) based on first-stage energy sink model
      CountryStation/study
      UF=University of Florida, ISU=Iowa State University, MSU=Michigan State University, FRC=USDA Forage Research Center, UW=University of Wisconsin-Madison, AGIL=USDA Animal Genomics Improvement Laboratory, TGEN=‘t Gen, NBZ=Nij Bosma Zaathe, ZOM=data from Zom et al. (2012), NLN=miscellaneous compilation of studies, LAN=Langhill, SAC=Scottish Agricultural College.
      MilkEMBWΔBWIntercept
      Expressed relative to a DIM=125, MilkE=25Mcal, MBW=121.3kg0.75, and ΔBW=0. N/A refers to estimates with unusually large SE or predictions beyond allowable parameter space (i.e., negative intercepts) because DIM range for that station did not include 125 d.
      United StatesUF0.35 ± 0.0160.16 ± 0.0070.13 ± 0.03N/A
      ISU0.36 ± 0.0050.10 ± 0.0030.30 ± 0.0322.39 ± 0.15
      MSU0.40 ± 0.0090.10 ± 0.0040.36 ± 0.0422.26 ± 0.21
      FRC0.34 ± 0.0130.12 ± 0.0070.18 ± 0.0622.74 ± 0.25
      UW0.44 ± 0.0080.11 ± 0.0040.27 ± 0.0423.39 ± 0.24
      AGIL0.37 ± 0.0080.08 ± 0.0040.39 ± 0.0822.43 ± 0.16
      NetherlandsTGEN0.40 ± 0.0050.10 ± 0.0020.08 ± 0.0123.29 ± 0.11
      NBZ0.47 ± 0.0100.07 ± 0.0100.27 ± 0.06N/A
      ZOM0.34 ± 0.0050.06 ± 0.0020.07 ± 0.0120.34 ± 0.12
      NLN0.29 ± 0.0050.07 ± 0.0030.18 ± 0.0320.40 ± 0.42
      United KingdomLAN0.35 ± 0.0030.06 ± 0.0010.35 ± 0.0217.83 ± 0.06
      SAC0.36 ± 0.0030.07 ± 0.0010.06 ± 0.0119.18 ± 0.14
      1 UF = University of Florida, ISU = Iowa State University, MSU = Michigan State University, FRC = USDA Forage Research Center, UW = University of Wisconsin-Madison, AGIL = USDA Animal Genomics Improvement Laboratory, TGEN = ‘t Gen, NBZ = Nij Bosma Zaathe, ZOM = data from
      • Zom R.L.G.
      • André G.
      • van Vuuren A.M.
      Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
      , NLN = miscellaneous compilation of studies, LAN = Langhill, SAC = Scottish Agricultural College.
      2 Expressed relative to a DIM = 125, MilkE = 25 Mcal, MBW = 121.3 kg0.75, and ΔBW = 0. N/A refers to estimates with unusually large SE or predictions beyond allowable parameter space (i.e., negative intercepts) because DIM range for that station did not include 125 d.
      We also observed universally statistically significant (P < 0.05) partial regression coefficients for DMI on MBW, ranging from 0.06 to 0.16 kg/kg0.75. The partial regression coefficients of DMI on ΔBW were also all deemed to be statistically significant (P < 0.05) for all research stations. Furthermore, all such partial regression coefficients were positive, intuitively suggesting that as ΔBW increases, DMI increases and vice versa; nevertheless, there was substantial range (0.07 to 0.36 kg of DMI/kg of ΔBW) in these estimates. These coefficients were particularly heterogeneous within the 2 UK stations (0.06 to 0.35). One could use the same z-test to suggest that some of these partial regression coefficients appeared to be different from each other. Again, this heterogeneity could be a reflection of differences in energy or statistical considerations, as addressed herein.
      Finally, station-specific intercepts are also reported in Table 2 relative to a baseline of DIM = 125, MilkE = 25 Mcal, MBW = 121.3 (i.e., BW = 600 kg), and ΔBW = 0 (i.e., these intercepts represent the overall estimated DMI for those particular values of DIM, MilkE, MBW, and ΔBW). Remarkable consistency was noted in these estimates (23.08 to 24.18 kg) among all US research stations, whereas these intercepts were estimated to be generally 2 to 5 kg less for some of the older data derived from the UK and NL. Note that using the NRC prediction equation (
      NRC
      Nutrient Requirements of Dairy Cattle. Natl.
      ), DMI would be predicted to be 23.85 kg for these same specifications of DIM, MilkE, and MBW; this prediction is thereby within 1.0 kg reported for each of the US research stations and a single Dutch station TGEN, but again substantially larger than the other estimated intercepts reported in Table 2.

      Quantitative Genetic Analysis of RFI from Two-Stage Analyses

      Estimated residuals from the first-stage energy sink model (Equation [2]) were used for the derivation of RFI for the second-stage quantitative genetic analyses (Equation [3]). The Legendre polynomial order fitted on the random regression terms on DIM for each of the additive and within-lactation permanent environmental effects were, respectively, up to na = 1 and np = 2 on UK, na = 0 (i.e., just intercept) and np = 1 on US, na = 2 and np = 3 on NL data, and na = 2 and np = 3 on the overall data. Based on these fits, DIM-specific genetic variances, within-lactation permanent environmental variances, and heritabilities were estimated. In addition, a between-lactation permanent environmental variance component σb2 was also estimated because some cows had multiple lactations. This term was important for the joint analysis σˆb2=0.33±0.05kg2 and for each country-specific analysis, but appeared to be substantially larger in the US σˆb2=0.51±0.17kg2 relative to NL σˆb2=0.22±0.06kg2 or UK σˆb2=0.32±0.07 kg2.
      Estimated genetic variances σg2 by DIM for each of the 3 countries and overall are plotted in Figure 2 as a function of DIM. Because no random polynomial regression term was significant for genetic variation in US data, σg2 was estimated to be constant across all DIM as σˆg2=0.56±0.15. Genetic variances from UK data were estimated to be lower than that for US data at earlier DIM, but then slightly higher at later DIM (>175 d). Genetic variability for NL was generally estimated to be substantially lower than that for US or UK, yet was rather stable as a function of DIM.
      Figure thumbnail gr2
      Figure 2Genetic variance for residual feed intake (RFI) as a function of DIM based on country-specific analyses [United Kingdom (○), United States (◊), and the Netherlands (∆)] and joint analyses of all data (□).
      Because of other estimated heterogeneous variance components, estimated heritabilities were not necessarily directly proportional to estimated genetic variances (Figure 3). For US and NL cows, heritabilities averaged around 10% throughout lactation, although the NL profile was more stable. However, the estimated heritabilities for UK data increased from just near 15% at 50 DIM to just less than 25% at later DIM. Higher heritabilities in the UK were generally due to estimated residual variances being substantially smaller in UK herds relative to the US and NL. However, lower heritabilities in the US were also partly due to large between-lactation permanent environmental variance components σb2, as reported earlier, as well as rather large within-lactation permanent environmental variances, particularly in the latter part (>150 d) of the DIM range (Figure 4). Note, however, the data in this range were relatively sparse for US (16% of all weekly records) compared with NL (24%) or UK (30%). Overall estimates of heritabilities ± SE across all DIM based on simpler random intercept versions (i.e., no random polynomial effects on DIM for additive genetics nor for within-lactation permanent environment) of these random regression models were determined to be 0.15 ± 0.03 (US), 0.09 ± 0.02 (NL), and 0.19 ± 0.04 (UK) with an overall heritability estimate of 0.15 ± 0.02 based upon a joint analysis across all countries.
      Figure thumbnail gr3
      Figure 3Heritability estimates for residual feed intake (RFI) as a function of DIM based on country-specific analyses [United Kingdom (○), United States (◊), and the Netherlands (∆)] and joint analyses of all data (□).
      Figure thumbnail gr4
      Figure 4Estimates of within-lactation permanent environmental variance for residual feed intake as a function of DIM based on country-specific analyses [United Kingdom (○), United States (◊), and the Netherlands (∆)] and joint analyses of all data (□).
      In spite of the relatively low heritabilities, within-lactation repeatabilities were remarkably high across all stages of lactation (Figure 5), ranging from as low as 0.55 for midlactation in NL cattle to as high as 0.90 in UK cattle with repeatability based on the joint analysis hovering consistently just under 0.80. However, as should be also noted from Figure 5, the repeatability of RFI across lactations (determined as the proportion of all variability that was due to genetics and between-lactation cow variability) was substantially lower (10 to 35%). This intuitively implies that RFI is more highly repeatable within lactations than across lactations.
      Figure thumbnail gr5
      Figure 5Estimated within-lactation (solid line) and between lactation (dotted line) repeatabilities for residual feed intake (RFI) as a function of DIM country-specific analyses [United Kingdom (○), United States (◊), and the Netherlands (∆)] and joint analyses of all data (□).

      One-Step Model Analysis

      As an alternative to the 2-stage approach, we also conducted a single-step model analysis based on DMI as a function of all effects in the first- (Equation [2]) and second-stage (Equation [3]) models and additional terms (Equation [4]) as described previously, focusing only on random intercept effects of additive genetics and between- and within-lactation permanent environment effects (i.e., no random regression terms). Estimated heritabilities based on this model (using the same definition as provided for the previous second stage analysis) were determined to be 0.18 ± 0.05 for US, 0.18 ± 0.04 for NL, 0.28 ± 0.05 for UK, and 0.18 ± 0.02 based on a joint analysis of all of the data. These estimates are somewhat greater than those previously reported for the more classical RFI analysis. The differences may be partly due to the joint accounting of all effects together in 1 model rather than in 2 separate models as well as modeling the additional interaction effects of ration with the partial regression coefficients on MilkE, MBW, and ΔBW.
      The station-specific estimated partial regression coefficients and intercepts based on the joint 1-stage analysis are provided in Table 3. This table is similar to Table 2 except that these estimated partial regression coefficients are adjusted for more effects (i.e., second-stage model effects) than were reported previously for the first stage model from Table 2. The standard errors in Table 3 are appreciably larger than those in Table 2, likely because of the proper accounting for all sources of (co)variation (i.e., the standard errors in Table 2 ignore the correlation between repeated measurements on cows). Taking that into account, the heterogeneity in these coefficients and intercepts between some research stations do not appear to be as statistically significant in Table 3 as they do in Table 2. The point estimates of the partial regression coefficients generally agree for the 2 approaches, with some stark exceptions. For example, the partial regression coefficients of DMI on MilkE are nearly halved for the 2 UK herds in Table 3 compared with Table 2.
      Table 3Estimated partial regression coefficients (±standard errors) of DMI on each of several energy sinks (MilkE = milk energy, MBW = metabolic BW, ΔBW = change in BW, BCS) based on 1-stage model analysis of all data
      CountryStation/study
      UF=University of Florida, ISU=Iowa State University, MSU=Michigan State University, FRC=USDA Forage Research Center, UW=University of Wisconsin-Madison, AGIL=USDA Animal Genomics Improvement Laboratory, TGEN=‘t Gen, NBZ=Nij Bosma Zaathe, ZOM=data from Zom et al. (2012), NLN=miscellaneous compilation of studies, LAN=Langhill, SAC=Scottish Agricultural College.
      MilkEMBWΔBWIntercept
      Expressed relative to a DIM=125, MilkE=25Mcal, MBW=121.3kg0.75, and ΔBW=0.
      United StatesUF0.34
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.033
      0.13
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.01
      0.14
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      23.64
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.39
      ISU0.34
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      0.11
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      0.21
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.06
      23.82
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.70
      MSU0.36
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      0.14
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.01
      0.25
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      23.82
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.38
      FRC0.37
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.03
      0.12
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      0.10
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.06
      23.74
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.39
      UW0.38
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      0.10
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.01
      0.08
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      24.18
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.29
      AGIL0.30
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.08
      0.11
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      0.06
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.14
      23.08
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 1.36
      NetherlandsTGEN0.26
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.08
      0.16
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      0.32
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.12
      24.35
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 1.36
      NBZ0.46
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.08
      0.02
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.04
      0.29
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.11
      22.27
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 1.17
      ZOM0.28
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      0.07
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.01
      0.08
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      21.80
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.26
      NLN0.18
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      0.10
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.01
      0.14
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.02
      20.92
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.27
      United KingdomLAN0.20
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.06
      0.07
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.03
      0.27
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.08
      17.83
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.97
      SAC0.18
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.06
      0.05
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.03
      0.04
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.08
      19.11
      Estimates not sharing the same letter within a column are concluded to be different (P<0.05) from each other.
       ± 0.97
      a–g Estimates not sharing the same letter within a column are concluded to be different (P < 0.05) from each other.
      1 UF = University of Florida, ISU = Iowa State University, MSU = Michigan State University, FRC = USDA Forage Research Center, UW = University of Wisconsin-Madison, AGIL = USDA Animal Genomics Improvement Laboratory, TGEN = ‘t Gen, NBZ = Nij Bosma Zaathe, ZOM = data from
      • Zom R.L.G.
      • André G.
      • van Vuuren A.M.
      Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
      , NLN = miscellaneous compilation of studies, LAN = Langhill, SAC = Scottish Agricultural College.
      2 Expressed relative to a DIM = 125, MilkE = 25 Mcal, MBW = 121.3 kg0.75, and ΔBW = 0.
      One might argue that ration and associated interaction effects should not be treated as random; however, the small number of records for some rations and partial confounding of ration effects with other model terms makes this a far more efficient modeling strategy than treating such effects as fixed. The estimated variance components for ration (i.e., intercept), ration × MBW, ration × MilkE, and ration × ΔBW were all statistically significant at 1.69 ± 0.21 kg2, 0.0012 ± 0.00022 (kg/kg0.75)2, 0.0060 ± 0.00095 (kg/Mcal)2, and 0.013 ± 0.0026, respectively. The estimated standard deviations of these effects can then be determined as the square roots of these variance component estimates; for example, 1.69=1.3 kg for intercept, 0.0012=0.035 kg/kg0.75 for the partial regression on MBW, 0.0060=0.077 kg/Mcal for partial regression on MilkE, and 0.013=0.11 for the partial regression on ΔBW. Thus, assuming that ration-specific partial regression coefficients on MilkE are normally distributed, at a mean of around 0.40 kg/Mcal (i.e., based on roughly averaging the second columns of either Tables 2 and 3), one should then roughly anticipate a range of ±2 SD about the mean between various rations [i.e., 0.40 ± 2 (0.077) = (0.25, 0.55) kg/Mcal for the partial regression of DMI on MilkE].
      A reinforcing illustration of this heterogeneity of the intercept and these partial regression coefficients across rations is obtained by expressing the random effect solutions relative to the corresponding research station baseline solutions previously provided in Table 3. These solutions are illustrated in Figure 6. These solutions partly indicate that in spite of the heterogeneity in these partial regression coefficients and intercepts between rations, the greatest level of heterogeneity still appears to be among research stations.
      Figure thumbnail gr6
      Figure 6Estimated ration-specific partial regression coefficients of DMI on (a) metabolic BW, (b) milk energy, and (c) BW changes for USDA Animal Genomics Improvement Laboratory (AGIL; ∆), University of Florida (UF; ◊), Iowa State University (ISU; ), Michigan State University (MSU; □), USDA Forage Research Center (FRC; ■), University of Wisconsin-Madison (UW; ▲), ‘t Gen (TGN; +), Nij Bosma Zaathe (NBZ; ⊕),
      • Zom R.L.G.
      • André G.
      • van Vuuren A.M.
      Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
      ; ZOM; *), miscellaneous compilation of studies (NLN; ×), Langhill (LAN; ○), and Scottish Agricultural College (SAC; ●). Color version available online.
      We further investigated the implications of the 2 different strategies for modeling this heterogeneity in station-specific partial regression coefficients and intercepts on prediction of genetic merit for feed efficiency by estimating the correlation in EBV for animals having phenotypes between the classical (2-stage) RFI approach and the 1-stage approach. Correlation between the 2 sets of EBV based on the joint analysis was 0.80; however, for each of the country-specific analyses, this correlation was as high as 0.96 for US data and as low as 0.80 for the UK data, whereas the correlation was 0.88 for the NL data.

      Discussion

      This work represents a comprehensive quantitative genetic analysis of dairy feed efficiency based on a large international data set. Although much attention has been directed toward the wide range in heritability estimates across different studies, there has been little discussion as to how heterogeneous the partial regression coefficients of DMI on the various key energy sinks (e.g., MilkE, MBW, and ΔBW) may be across various management systems. Finally, little consideration has been taken as to how genetic and nongenetic sources of variability for RFI may differ across lactation, even though some effort to identify potential windows within a lactation upon which to focus selection for RFI has occurred (
      • Connor E.E.
      • Hutchison J.L.
      • Olson K.M.
      • Norman H.D.
      Opportunities for improving milk production efficiency in dairy cattle.
      ).

      Potential Heterogeneity in the Partial Regression Coefficients

      Regarding our first objective, we concluded that the partial regression coefficients of DMI on various energy sinks are somewhat heterogeneous across research facilities or studies. Smaller partial regression coefficients (b1) of DMI on MilkE for some research stations might imply that these cows may be more efficient at converting DMI to MilkE. Nevertheless, partial regression coefficients for DMI on MilkE for most other stations or studies were in reasonable agreement with those presented by the National Research Council (NRC, 2001). On average, the partial regression coefficients for MilkE from the first-stage model (Table 2) are less than those reported on page 4 of
      NRC
      Nutrient Requirements of Dairy Cattle. Natl.
      ), which reports a partial regression coefficient for DMI (kg) on FCM of 0.372 (which roughly translates to 0.5 kg DMI/Mcal of MilkE). However, it should also be noted that the NRC equation is based on forcing a 0 intercept. At first glance, this heterogeneity suggests that cows at some stations are more efficient at converting DMI to MilkE; however, other possibilities include differences in diet energy densities, given that DMI was recorded in kilograms and not megacalories of energy, or potential systematic measurement errors.
      We also note that partial regression coefficients on the other energy sinks (MBW and ΔBW) were also heterogeneous; although, on average, estimated partial regression coefficients of DMI on MBW in Table 2 were generally consistent with that reported by
      NRC
      Nutrient Requirements of Dairy Cattle. Natl.
      ), which suggested a linear relationship between DMI and MBW based on a partial regression coefficient of 0.0968 kg/kg0.75 at midlactation. Our results imply that great care may be needed to create the RFI variable for subsequent quantitative genetic analysis. This recognition was made by
      • Davis S.R.
      • Macdonald K.A.
      • Waghorn G.C.
      • Spelman R.J.
      Residual feed intake of lactating Holstein-Friesian cows predicted from high-density genotypes and phenotyping of growing heifers.
      , who derived RFI separately for each of 4 consecutive trials of New Zealand Holstein Friesian cattle, noting substantial heterogeneity in these partial regression coefficients across trials even though cows were all presumably fed the same diet across these trials. Furthermore, some of their significant partial regression coefficients were somewhat nonintuitive (i.e., they noted significant negative relationships between DMI and ΔBW). As some researchers are now starting to further explore genetic variability in these partial relationships in RFI studies involving beef cattle (
      • Savietto D.
      • Berry D.P.
      • Friggens N.C.
      Towards an improved estimation of the biological components of residual feed intake in growing cattle.
      ) and poultry (
      • Aggrey S.E.
      • Rekaya R.
      Dissection of Koch's residual feed intake: Implications for selection.
      ), it will be even more important to accurately characterize such heterogeneities in future studies.
      Some partial coefficients might be biased if some of the key energy sink covariates (e.g., BW) are measured with some error (
      • Robinson D.L.
      Accounting for bias in regression coefficients with example from feed efficiency.
      ). Caution might also be warranted as to the potential existence of nonlinear associations (
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      ) between DMI and the various energy sinks, although our data did not suggest strong deviations from partial linear relationships (results not shown). For instance, estimated station-specific intercepts were made with reference to particular baseline on DIM (125 d), BW (600 kg), MilkE (25 Mcal), and ΔBW (0), although some of those specifications might be extrapolations for some of the older cows in the data set. Furthermore, the NRC prediction equation appeared to be in better agreement with analyses based on some of the newer data. Other closely aligned efforts are currently taking place, such as the gDMI (global Dry Matter Initiative) project that is focusing on DMI rather than RFI as a key response variable (
      • Berry D.P.
      • Coffey M.P.
      • Pryce J.E.
      • de Haas Y.
      • Lovendahl P.
      • Krattenmacher N.
      • Crowley J.J.
      • Wang Z.
      • Spurlock D.
      • Weigel K.
      • Macdonald K.
      • Veerkamp R.F.
      International genetic evaluations for feed intake in dairy cattle through the collation of data from multiple sources.
      ). Some discussion on the arguments in favor and against the use of either trait in selection programs is provided by
      • Berry D.
      • Pryce J.
      Feed efficiency in growing and mature animals.
      , whereas the quantitative genetic modeling equivalence between the use of DMI versus RFI within a multiple trait analysis framework is discussed by

      Lu, Y. F., M. J. Vandehaar, K. A. Weigel, L. E. Armentano, D. M. Spurlock, C. R. Staples, E. E. Connor, and R. J. Tempelman. 2014. An alternative approach to modeling genetic merit of feed efficiency in dairy cattle. Proc. 10th World Congr. Genet. Appl. Livest. Prod., Vancouver, Canada. Am. Soc. Anim. Sci., Champaign, IL.

      . Our results on heterogeneity of partial regression coefficients would indirectly suggest that there should be substantial heterogeneity in genetic or residual correlations between DMI and various energy sink traits in multiple trait analysis extensions of that work.

      Genetic and Permanent Environmental Variability in RFI

      Genetic variances for RFI were somewhat heterogeneous across the 3 countries, estimated as being substantially higher in the US relative to the other countries between 50 and 125 DIM but then being lower after 150 DIM. Genetic variances did not show a consistent trend across DIM except for a general tendency to increase at later DIM in the other countries. Our inferences on genetic variance components did not track well with estimated heritabilities for RFI across lactation in part due to a high degree of heterogeneity in residual variances across research stations as well as large across-country differences in between-lactation permanent environmental variances and in within-lactation permanent environmental variances across DIM. Heritabilities in NL and US tracked between 0.05 and 0.15 and were relatively stable across DIM, whereas jointly estimated and UK heritabilities appeared to be somewhat larger, especially later in lactation. These heritability estimates are in general agreement with some previous work but somewhat larger than the average of 0.04 for lactating RFI based on an average of 11 studies as reported by
      • Berry D.P.
      • Crowley J.J.
      Genetics of feed efficiency in dairy and beef cattle.
      .
      A relatively large difference was noted between somewhat low heritability estimates yet rather high repeatability estimates, suggesting that large components of nonadditive genetic or common environmental variability create consistency in RFI values within cows throughout lactation beyond that due to additive genetic inheritance. The shape of the repeatability curves provided in Figure 4 were primarily driven by the shape of the within-lactation permanent environmental variance, tending to be larger at the more extreme DIM with a nadir at intermediate DIM. However, the results should be treated with some caution given that random regression models can lead to unusual behavior in estimated genetic and permanent environmental variances at the extremes if the data are sparse in those areas (
      • Misztal I.
      Properties of random regression models using linear splines.
      ). Recall that the US data had proportionally far less data (16% of all weekly records) between 150 and 200 DIM compared with UK (30%) or NL (24%) data, such that inferences on variance components for the US in that DIM range should be treated with caution.
      • Connor E.E.
      • Hutchison J.L.
      • Norman H.D.
      • Olson K.M.
      • Van Tassell C.P.
      • Leith J.M.
      • Baldwin R.L.
      Use of residual feed intake in Holsteins during early lactation shows potential to improve feed efficiency through genetic selection.
      recently reported a moderately high repeatability of 0.47.
      We also considered a 1-stage model analysis to further complement our 2-stage model analysis. We noticed that heritabilities increased slightly based on the 1-step approach, perhaps reflecting the fact that a 1-step model accounts for all sources of variation somewhat more efficiently, or because it included ration interaction effects with MilkE, MBW, and ΔBW. Differences in heritability estimates may occur because some effects or factors in the first-stage (energy sink) model are not orthogonal (i.e., balanced) with respect to factors in the second-stage (RFI) model. In our specific application, this would imply that cow effects (genetics and permanent environment) are not orthogonal to other systematic effects in the first-stage model. This is not too surprising given that cow is an important source of variability for some of the energy sink terms (e.g., MilkE and MBW) in the first-stage model. As another example, a 2-step approach might suboptimally adjust for test week effects, as genetic relationships between animals having records in different test weeks would not be accounted for in the first-stage model. This lack of orthogonality may be partly the reason why some estimated partial regression coefficients were estimated to be dramatically different between the 2 approaches. Furthermore, it appears that rations are important sources of variability that further influence the partial relationships between DMI and the various energy sinks beyond that due to station. These issues will be addressed in future research based on alternative modeling approaches for feed efficiency in dairy cattle.

      Conclusions

      Based on an international study involving 4,893 lactating Holstein cows, we concluded that genetic variation in RFI exists and that it is moderately heritable at about 15 to 18% overall depending upon the modeling approach, thereby suggesting that genomic selection strategies on RFI will be fruitful. Furthermore, heritability in RFI varies across lactation and is moderately different from one region to the next. Moreover, the partial regression relationships used to derive RFI are heterogeneous across countries and across stage of lactation. The heterogeneity in variation of these relationships could lead to bias in future genetic or genomic inferences for RFI, if ignored, and should be considered in future attempts to select dairy cattle for efficiency. Other potential sources of heterogeneity (e.g., interaction of partial regression coefficients of DMI on energy sinks with rations) should be explored routinely as well.

      Acknowledgments

      The authors acknowledge funding from USDA-NIFA Grant Number 2011-68004-30340 . CRV (Arnhem, the Netherlands) and the WhyDry-project (Wageningen University, the Netherlands) are acknowledged for providing part of the Dutch data. UK data collection was funded by Scottish Government .

      References

        • Aggrey S.E.
        • Rekaya R.
        Dissection of Koch's residual feed intake: Implications for selection.
        Poult. Sci. 2013; 92: 2600-2605
        • Banos G.
        • Coffey M.P.
        • Veerkamp R.F.
        • Berry D.P.
        • Wall E.
        Merging and characterising phenotypic data on conventional and rare traits from dairy cattle experimental resources in three countries.
        Animal. 2012; 6: 1040-1048
        • Berry D.
        • Pryce J.
        Feed efficiency in growing and mature animals.
        in: Proc. 10th World Congr. Genet. Appl. Livest. Prod., Vancouver, Canada. Am. Soc. Anim. Sci, Champaign, IL2014
        • Berry D.P.
        • Coffey M.P.
        • Pryce J.E.
        • de Haas Y.
        • Lovendahl P.
        • Krattenmacher N.
        • Crowley J.J.
        • Wang Z.
        • Spurlock D.
        • Weigel K.
        • Macdonald K.
        • Veerkamp R.F.
        International genetic evaluations for feed intake in dairy cattle through the collation of data from multiple sources.
        J. Dairy Sci. 2014; 97: 3894-3905
        • Berry D.P.
        • Crowley J.J.
        Genetics of feed efficiency in dairy and beef cattle.
        J. Anim. Sci. 2013; 91: 1594-1613
        • Connor E.E.
        • Hutchison J.L.
        • Norman H.D.
        • Olson K.M.
        • Van Tassell C.P.
        • Leith J.M.
        • Baldwin R.L.
        Use of residual feed intake in Holsteins during early lactation shows potential to improve feed efficiency through genetic selection.
        J. Anim. Sci. 2013; 91: 3978-3988
        • Connor E.E.
        • Hutchison J.L.
        • Olson K.M.
        • Norman H.D.
        Opportunities for improving milk production efficiency in dairy cattle.
        J. Anim. Sci. 2012; 90: 1687-1694
        • Davis S.R.
        • Macdonald K.A.
        • Waghorn G.C.
        • Spelman R.J.
        Residual feed intake of lactating Holstein-Friesian cows predicted from high-density genotypes and phenotyping of growing heifers.
        J. Dairy Sci. 2014; 97: 1436-1445
        • Ferraretto L.F.
        • Shaver R.D.
        • Bertics S.J.
        Effect of dietary supplementation with live-cell yeast at two dosages on lactation performance, ruminal fermentation and total-tract nutrient digestibility in dairy cows.
        J. Dairy Sci. 2012; 95: 4017-4028
        • Ferraretto L.F.
        • Shaver R.D.
        • Espineira M.
        • Gencoglu H.
        • Bertics S.J.
        Influence of a reduced-starch diet with or without exogenous amylase on lactation performance by dairy cows.
        J. Dairy Sci. 2011; 94: 1490-1499
        • Gilmour A.R.
        • Gogel B.J.
        • Cullis B.R.
        • Thompson R.
        ASREML User's Guide Release 3.0. VSN International, Hemel Hempstead, U.K2009
        • Gonzalez-Recio O.
        • Pryce J.E.
        • Haile-Mariam M.
        • Hayes B.J.
        Incorporating heifer feed efficiency in the Australian selection index using genomic selection.
        J. Dairy Sci. 2014; 97: 3883-3893
        • He M.
        • Perfield K.L.
        • Green H.B.
        • Armentano L.E.
        Effect of dietary fat blend enriched in oleic or linoleic acid and monensin supplementation on dairy cattle performance, milk fatty acid profiles, and milk fat depression.
        J. Dairy Sci. 2012; 95: 1447-1461
        • Hill W.G.
        On selection among groups with heterogeneous variance.
        Anim. Sci. 1984; 39: 473-477
        • Klein T.W.
        Heritability and genetic correlation: Statistical power, population comparisons, and sample size.
        Behav. Genet. 1974; 4: 171-189
      1. Lu, Y. F., M. J. Vandehaar, K. A. Weigel, L. E. Armentano, D. M. Spurlock, C. R. Staples, E. E. Connor, and R. J. Tempelman. 2014. An alternative approach to modeling genetic merit of feed efficiency in dairy cattle. Proc. 10th World Congr. Genet. Appl. Livest. Prod., Vancouver, Canada. Am. Soc. Anim. Sci., Champaign, IL.

        • Manzanilla Pech C.I.V.
        • Veerkamp R.F.
        • Calus M.P.L.
        • Zom R.
        • van Knegsel A.
        • Pryce J.E.
        • de Haas Y.
        Genetic parameters across lactation for feed intake, fat- and protein-corrected milk, and liveweight in first-parity Holstein cattle.
        J. Dairy Sci. 2014; 97: 5851-5862
        • Misztal I.
        Properties of random regression models using linear splines.
        J. Anim. Breed. Genet. 2006; 123: 74-80
        • NRC
        Nutrient Requirements of Dairy Cattle. Natl.
        Acad. Press, Washington, DC2001
        • Pryce J.E.
        • Arias J.
        • Bowman P.J.
        • Davis S.R.
        • Macdonald K.A.
        • Waghorn G.C.
        • Wales W.J.
        • Williams Y.J.
        • Spelman R.J.
        • Hayes B.J.
        Accuracy of genomic predictions of residual feed intake and 250-day body weight in growing heifers using 625,000 single nucleotide polymorphism markers.
        J. Dairy Sci. 2012; 95: 2108-2119
        • Quaas R.L.
        Additive genetic model with groups and relationships.
        J. Dairy Sci. 1988; 71: 1338-1345
        • Robinson D.L.
        Accounting for bias in regression coefficients with example from feed efficiency.
        Livest. Prod. Sci. 2005; 95: 155-161
        • Savietto D.
        • Berry D.P.
        • Friggens N.C.
        Towards an improved estimation of the biological components of residual feed intake in growing cattle.
        J. Anim. Sci. 2014; 92: 467-476
        • Schaeffer L.R.
        Application of random regression models in animal breeding.
        Livest. Prod. Sci. 2004; 86: 35-45
        • Spurlock D.M.
        • Dekkers J.C.
        • Fernando R.
        • Koltes D.A.
        • Wolc A.
        Genetic parameters for energy balance, feed efficiency, and related traits in Holstein cattle.
        J. Dairy Sci. 2012; 95: 5393-5402
        • Vallimont J.E.
        • Dechow C.D.
        • Daubert J.M.
        • Dekleva M.W.
        • Blum J.W.
        • Barlieb C.M.
        • Liu W.
        • Varga G.A.
        • Heinrichs A.J.
        • Baumrucker C.R.
        Heritability of gross feed efficiency and associations with yield, intake, residual intake, body weight, and body condition score in 11 commercial Pennsylvania tie stalls.
        J. Dairy Sci. 2011; 94: 2108-2113
        • van der Werf J.H.J.
        • Goddard M.E.
        • Meyer K.
        The use of covariance functions and random regressions for genetic evaluation of milk production based on test day records.
        J. Dairy Sci. 1998; 81: 3300-3308
        • VandeHaar M.J.
        • St-Pierre N.
        Major advances in nutrition: Relevance to the sustainability of the dairy industry.
        J. Dairy Sci. 2006; 89: 1280-1291
        • Veerkamp R.F.
        • Emmans G.C.
        • Cromie A.R.
        • Simm G.
        Variance components for residual feed intake in dairy cows.
        Livest. Prod. Sci. 1995; 41: 111-120
        • Veerkamp R.F.
        • Oldenbroek J.K.
        • Van Der Gaast H.J.
        • van der Werf J.H.J.
        Genetic correlation between days until start of luteal activity and milk yield, energy balance, and live weights.
        J. Dairy Sci. 2000; 83: 577-583
        • Williams Y.J.
        • Pryce J.E.
        • Grainger C.
        • Wales W.J.
        • Linden N.
        • Porker M.
        • Hayes B.J.
        Variation in residual feed intake in Holstein-Friesian dairy heifers in southern Australia.
        J. Dairy Sci. 2011; 94: 4715-4725
        • Yao C.
        • Spurlock D.M.
        • Armentano L.E.
        • Page C.D.
        • VandeHaar M.J.
        • Bickhart D.M.
        • Weigel K.A.
        Random Forests approach for identifying additive and epistatic single nucleotide polymorphisms associated with residual feed intake in dairy cattle.
        J. Dairy Sci. 2013; 96: 6716-6729
        • Zom R.L.G.
        • André G.
        • van Vuuren A.M.
        Development of a model for the prediction of feed intake by dairy cows: 1. Prediction of feed intake.
        Livest. Sci. 2012; 143: 43-57