Impact of parity differences on residual feed intake estimation in Holstein cows

Graphical Abstract Summary: Residual feed intake (RFI) is typically obtained as the difference between dry matter intake (DMI) observations and predictions from regression on known energy sinks, and effects of parity, days in milk, and cohort. Understanding the impact of parity on RFI estimation is important because variation across parities might affect the estimation of RFI phenotypes and, hence, the prediction of RFI breeding values. To investigate this, we calculated RFI considering the energy sinks nested or not nested (constant) within parity using weekly records of 5,813 lactating Holstein cows collected from 2007 to 2022 in the United States. In addition, genetic correlations between weekly RFI for parities 1, 2, and 3 were estimated using bivariate repeatability animal models as the genetic breeding values for RFI in each parity. Our results showed that nesting energy sinks within parity when computing RFI improves model goodness of fit, but the impact on the RFI phenotypes and the estimated breeding values appears to be minimal.

Abstract: Residual feed intake (RFI) has been used as a measure of feed efficiency in farm animals. In lactating dairy cattle, RFI is typically obtained as the difference between dry matter intake observations and predictions from regression on known energy sinks, and effects of parity, days in milk, and cohort. The impact of parity (lactation number) on the estimation of RFI is not well understood, so the objectives of this study were to (1) evaluate alternative RFI models in which the energy sinks (metabolic body weight, body weight change, and secreted milk energy) were nested or not nested within parity, and (2) estimate variance components and genetic correlations for RFI across parities. Data consisted of 72,474 weekly RFI records of 5,813 lactating Holstein cows collected from 2007 to 2022 in 5 research stations across the United States. Estimates of heritability, repeatability, and genetic correlations between weekly RFI for parities 1, 2, and 3 were obtained using bivariate repeatability animal models. The nested RFI model showed better goodness of fit than the nonnested model, and some partial regression coefficients of dry matter intake on energy sinks were heterogeneous between parities. However, the Spearman's rank correlation between RFI values calculated from nested and nonnested models was equal to 0.99. Similarly, Spearman's rank correlation between the RFI breeding values from these 2 models was equal to 0.98. Heritability estimates for RFI were equal to 0.16 for parity 1, 0.19 for parity 2, and 0.22 for parity 3. Repeatability estimates for RFI across weeks within parities were high, ranging from 0.51 to 0.57. Spearman's rank correlations of sires' breeding values were 0.99 between parities 1 and 2, 0.91 between parities 1 and 3, and 0.92 between parities 2 and 3. We conclude that nesting energy sinks within parity when computing RFI improves model goodness of fit, but the impact on the estimated breading values appears to be minimal.
I dentifying cows that can convert feed to milk more efficiently than their contemporaries is challenging due to the biological complexity of feed efficiency. Residual feed intake (RFI) has been used widely as a measure of feed efficiency in cattle, and in lactating dairy cows it is usually obtained from regression of observed DMI on known energy sinks, namely metabolic body weight (mBW), body weight change (ΔBW), and secreted milk energy (MilkE). In addition, most RFI models account for parity, DIM, and cohort effects. The residual portion of this model represents the RFI phenotypes (Berry and Crowley, 2013;Tempelman et al., 2015;VandeHaar et al., 2016). Therefore, for a given population, more efficient cows eat less than predicted and have lower RFI values compared with the average of their contemporaries.
In general, RFI estimation models consider parity effects independent of the energy sinks (i.e., partial regression coefficients of DMI on energy sinks are constant across parities). However, these partial regression coefficients could vary across parities, particularly for primiparous cows that are still growing, and ignoring these differences might affect the estimation of RFI phenotypes, and hence, the prediction of RFI breeding values. The objectives of this study were to evaluate the impact of parity on RFI estimation in lactating Holstein cows by comparing RFI models considering energy sinks nested within parity or constant across parities, and to estimate genetic variance components and genetic correlations for RFI across different parities.
Data consisted of 72,474 weekly RFI records of 5,813 lactating Holstein cows (37,211 records of 3,584 cows in parity 1, 23,836 records of 2,324 cows in parity 2, and 11,427 records of 1,110 cows in parity 3) collected from 2007 to 2022 at 5 research stations: Iowa State University (Ames, IA), Michigan State University (East Lansing, MI), University of Florida (Gainesville, FL), University of Wisconsin-Madison (Arlington, WI), and USDA-Agricultural Research Service Animal Genomics and Improvement Laboratory (Beltsville, MD). All procedures were approved by the corresponding Animal Care and Use Committees. Animals were housed in freestall or tiestall facilities, and daily feed intakes were measured via roughage intake control system (Hokofarm Group), Calan Broadbent Feeding System (American Calan), GrowSafe System (Vytelle), or manual weigh-back of refusals. Milk weights were obtained daily, milk samples were obtained weekly for determination of milk composition, and BW were obtained on 3 consecutive days at the beginning, middle, and end of the experimental period. Milk energy (MilkE; Mcal) was calculated weekly using the following equation (NRC, 2001): Since the experiments lasted only a few weeks, weekly BW were estimated using a simple linear regression of measured BW on days for each experiment. The mBW was calculated as weekly average BW 0.75 , and ΔBW was calculated as the difference in predicted BW at the end and beginning of each week.

Impact of parity differences on residual feed intake estimation in Holstein cows
Weekly RFI phenotypes were calculated using 2 alternative linear mixed models: where DIM represents the effect of days in milk with 9 levels (grouped by 16 d), Lact represents the effect of parity (lactation number) with 3 levels (1, 2, and 3), MilkE is secreted milk energy with partial regression coefficient b 1 , mBW is metabolic body weight with partial regression coefficient b 2 , ΔBW is the change in body weight with partial regression coefficient b 3 , block represents the random effect of experiment-treatment, week represents the random effect of week of experiment, and e is the random residual of the model, representing RFI. Random effects were assumed to follow a multivariate normal distribution, with block N 0 2 Iσ e ( ) where I is the identity matrix, and covariances between these effects equal to zero. Therefore, the difference between these 2 alternative models consisted of nesting regression coefficients for DIM and the 3 energy sinks within parity in the nested model, versus treating them as constant across parities in the constant model. The mean and standard deviation of RFI were 0 ± 2.0 kg in both models, whereas the means and standard deviations of DMI, DIM, MilkE, mBW, and ΔBW were 24.8 ± 4.2 kg, 115.0 ± 37.6 d, 29.9 ± 5.8 Mcal, 126.0 ± 11.8 kg 0.75 , and 2.7 ± 3.0 kg, respectively. Goodness of fit was assessed using Bayesian information criterion and Akaike information criterion. We also evaluated the Spearman's rank correlation coefficient between RFI phenotypes from the 2 models, and compared partial regression coefficients for energy sinks across parities.
The impact of parity on RFI estimation was also evaluated in a genetic context. Weekly RFI phenotypes were calculated fitting a separate RFI model within each parity (1, 2, and 3) using the constant model described above, but with no parity effect as RFI within each parity was considered a unique trait. Estimates of heritability, repeatability, and genetic correlations between weekly RFI values for parities 1, 2, and 3 were performed using a bivariate repeatability animal model: where y is a vector of weekly RFI records, β is a vector of fixed effects of week, u is a vector of random additive genetic effects, pe is a vector of random permanent environmental effects, and e is the vector of random residual effects. Matrices X, Z, and W are incidence matrices relating y to β, u, and pe, respectively. Random effects were assumed to follow a multivariate normal distribution, u pe e where G 0 and P 0 are the additive genetic and permanent environmental (co)variance matrices, respectively; A is the matrix of additive relationships between animals based on 5 generations of pedigree information; R 0 is the residual (co)variance matrix, and I is an identity matrix with suitable dimensions. Variance component estimates were obtained using REML with the AIREMLF90 software (Aguilar et al., 2018). The Akaike information criterion and Bayesian information criterion values were 315,725 and 315,881, respectively, for the constant model, and 315,365 and 315,705, respectively, for the nested model, differing statistically by a chi-squared test (P-value = 2.2 −16 ). Thus, the nested model provided superior fit. However, Spearman's rank correlation coefficient of RFI phenotypes obtained from the constant model and nested model indicated minimal re-ranking ( Figure 1A, correlation = 0.99). In addition, Spearman's rank correlation between the RFI breeding values derived from these 2 models was equal to 0.98, showing also a minimal re-ranking for the estimated breeding values.
Estimated partial regression coefficients of DMI on the energy sinks ranged from 0.376 to 0.382 kg of DMI per Mcal of MilkE, 0.107 to 0.112 kg of DMI per kg 0.75 of mBW, and 0.086 to 0.124 kg of DMI per kg of ∆BW (Table 1). Based on the corresponding confidence intervals, some partial regression coefficients were heterogeneous across parities (nested model, Table 1). For MilkE, there was no difference in the partial regression coefficients across parities. However, for mBW, the partial regression coefficient for parity 2 were larger than that of parities 1 and 3, which were similar. Last, for ∆BW, the partial regression coefficient for parity 3 was significantly larger than that of parities 1 and 2, which were numerically similar. On average, partial regression coefficients were similar to those obtained in other feed efficiency studies in Holstein cows, many of which used overlapping data sets (Tempelman et al., 2015;Khanal et al., 2022), and also in a study using data from European research stations (Guinguina et al., 2020). Partial regression coefficients for DMI on ∆BW can be different across studies and difficult to compare due to differences in energy content of feed intake across populations/studies. As reported by Tempelman et al. (2015), variation in ∆BW regression coefficients seems to be greater than that of the other energy sinks when compared across models, data sets, and research stations, perhaps due to imprecision in ΔBW phenotypes measured over relatively short time periods. Negative partial regression coefficient for ∆BW was reported by Davis et al. (2014) in a study with lactating Holstein-Friesian cows.
The next step of our study was to evaluate potential genetic differences in RFI among parities. For this, we estimated genetic parameters and genetic correlations among parities 1, 2, and 3. Table  2 shows heritability, repeatability, and genetic correlation estimates for RFI in each of these 3 parities. Heritability estimates for RFI were equal to 0.16 ± 0.02 in parity 1, 0.19 ± 0.03 in parity 2, and 0.22 ± 0.06 in parity 3. Moreover, RFI phenotypes were repeatable across weeks in all parities, with repeatability estimates ranging from 0.51 to 0.57. Estimated genetic correlation for RFI were 0.98 ± 0.06 between parities 1 and 2, 0.79 ± 0.21 between parities 1 and 3, and 0.85 ± 0.18 between parities 2 and 3, indicating poorer precision for pairwise combinations involving parity 3. Complementing these results, Figure 1 (B, C, D) depict the relationship between sires' estimated breeding values for each parity, where the Spearman's rank correlation coefficient was 0.99 between parities 1 and 2, 0.91 between parities 1 and 3, and 0.92 between parities 2 and 3.
Given fixed experiment station resources for measurement of DMI, MilkE, mBW, ΔBW, and RFI phenotypes for national genetic evaluations of feed efficiency in dairy cattle, the structure of the data is limited, in terms of the number of cows in multiple lactations. Measurement of individual feed intake records is expensive and time consuming, and it is more valuable to measure phenotypes on a greater number of unique cows than to measure  repeated phenotypes on the same cows in different lactations. The numbers of cows with repeated phenotypes in parities 1 and 2, parities 1 and 3, and parities 2 and 3 were 713, 222, and 266, respectively. This data structure resulted in larger standard errors (0.18 and 0.21) for the estimated genetic correlations between parities 2 and 3 and parities 1 and 3, and it may have contributed to the uncertainty of the genetic correlation (0.79) between parities 1 and 3, which presented the largest SE. Despite the importance of feed efficiency in dairy cattle, research on the potential differences in RFI phenotypes and predicted breeding values across parities have been limited. Khanal et al. (2022) reported a difference in heritability estimates of RFI for primiparous (0.25 ± 0.08) and multiparous (0.17 ± 0.05) Holstein cows, as well as differences in the partial regression coefficients. Li et al. (2017) observed variation in partial regression coefficients and genetic variance components between different periods within the first lactation. In agreement with Lu et al. (2017), they noticed greater variation in early and late lactation, with only small differences between 50 and 200 DIM, the minimum and maximum used in the present study. Also, Martin et al. (2021b) observed that the averaged RFI over the lactation was highly correlated with RFI measured in mid lactation. Although our data are mainly from mid-lactation cows, future studies implementing flexible models, as random regression models, may be interesting to investigate the variation across time within lactation. Moreover, heritability estimates for RFI in lactating Holstein cows have been reported from 0.15 to 0.32 (Gonzalez-Recio et al., 2014;Tempelman et al., 2015;Khanal et al., 2022).
Overall, our findings suggest that the nested model is more elegant and provides greater goodness of fit, and it has the added advantage of providing a more detailed biological interpretation of differences in partial regression coefficients of DMI on known energy sinks within each parity. For these reasons, the nested model should be preferred for analysis of the effects of dietary treatments and other experimental interventions on RFI in lactating dairy cows (Martin et al., 2021a). However, considering regressions on the energy sinks as constant or nested within parity had a minimal impact on RFI phenotypes and resulted in negligible changes in cows' rankings by predicted breeding values. Strong and positive genetic correlations were observed between parities 1 and 2, and sires' predicted breeding values were relatively consistent across lactations, but our ability to draw conclusive inferences about genetic correlations between parities 1 and 3, and to a lesser extent parities 2 and 3, was hampered by the data structure in which repeated phenotypes of the same cows across multiple lactations were limited by costs and facility constraints.