Genomic evaluation of feed efficiency in US Holstein heifers

There is growing interest in improving feed efficiency traits in dairy cattle. The objectives of this study were to estimate the genetic parameters of residual feed intake (RFI) and its component traits [dry matter intake (DMI), metabolic body weight (MBW), and average daily gain (ADG)] in Holstein heifers, and to develop a system for genomic evaluation for RFI in Holstein dairy calves. The RFI data were collected from 6,563 growing Holstein heifers (initial body weight = 261 ± 52 kg; initial age = 266 ± 42 d) for 70 d, across 182 trials conducted between 2014 and 2022 at the STgenetics Ohio Heifer Center (South Charleston, OH) as part of the EcoFeed program, which aims to improve feed efficiency by genetic selection. The RFI was estimated as the difference between a heifer’s actual feed intake and expected feed intake, which was determined by regression of DMI against midpoint MBW, age, and ADG across each trial. A total of 61,283 SNPs were used in genomic analyses. Animals with phenotypes and genotypes were used as training population, and 4 groups of prediction population, each with 2,000 animals, were selected from a pool of Holstein animals with genotypes, based on their relationship with the training population. All traits were analyzed using univariate animal model in DMU version 6 software. Pedigree information and genomic information were used to specify genetic relationships to estimate the variance components and genomic estimated breeding values (GEBV), respectively. Breeding values of the prediction population were estimated by using the 2-step approach: deriving the prediction equation of GEBV from the training population for estimation of GEBV of prediction population with only genotypes. Reliability of breeding values was obtained by approximation based on partitioning a function of the accuracy of training population GEBV and magnitudes of genomic relationships between individuals in the training and prediction population. Heifers had DMI (mean ± SD) of 8.11 ± 1.59 kg over the trial period, with growth


INTRODUCTION
Emerging technologies and innovative research methodologies continue to provide dairy producers with new opportunities to improve production efficiency and, therefore, profitability of their operations.One opportunity of immense interest today involves identification of animals that eat less feed while maintaining similar production to their herdmates.Such an animal has a superior ability to efficiently convert feed energy into usable products such as meat or milk.Given that feed accounts for over 50% of total farm costs (USDA ERS, 2022), the efficient use of feed by dairy cattle is critical for the economic success of dairy operations.Feed conversion efficiency also directly affects the environmental sustainability of farms, as animals that eat less while maintaining production have reduced manure and methane emissions per unit of product, thereby providing farmers with a solution to increase both profitability and sustainability (Nkrumah et al., 2004;Hafla et al., 2012;Basarab et al., 2013).Accordingly, significant research efforts have been dedicated to improving the feed conversion efficiency of dairy cattle (Islam et al., 2020;Li et al., 2020;Khanal et al., 2022).
Improvement in feed conversion efficiency has been limited by the relatively high cost and difficulty of accurately quantifying individual feed intake in cattle.
One method that has gained traction in recent years is genomic selection for feed conversion efficiency, which would be cumulative across generations and permanent (Pryce et al., 2015;Connor et al., 2019;Li et al., 2020).As with other hard-to-measure traits, genomic selection offers a unique opportunity to generate reliable breeding expectations for a population based on data collected from a representative subset.Limitations to this technique include the challenges of collecting accurate phenotypes from a large reference population, to provide high-accuracy predictions on the larger population (Pryce et al., 2015).However, given the wide adaptation of genomic testing, this technique provides the most realistic opportunity to enable widespread selection for feed conversion efficiency across the dairy industry.
Traditionally, genetic selection for feed efficiency has focused on favorable selection for ratio traits such as feed-to-gain and feed conversion ratio.Favorable selection for ratio-based traits can lead to greater mature body size and subsequently greater maintenance energy requirements due to their strong association with growth traits (Herd and Bishop, 2000;Arthur et al., 2001).Therefore, such selection practices can result in animals with greater maintenance feed intake requirements and subsequently higher methane production, given the strong ties between feed intake and methane production (Shibata and Terada, 2010).Instead, residual feed intake (RFI) is a more suitable trait for use in selection programs for improved feed efficiency, as it is moderately heritable (VandeHaar et al., 2016;Li et al., 2020) and accounts for the variation in individual animals' feed efficiency, independent of growth and production.
Measuring RFI in growing heifers, as opposed to dairy cows, is attractive because cows undergo the problem of negative energy balance and tissue mobilization during lactation (Williams et al., 2011;Pryce et al., 2014).Accordingly, ideal models for RFI in cows should include adjustments for body reserves and information on milk components, milk yields (Pryce et al., 2014), and status of pregnancy, which could have potential effects on energy partitioning (Khanal et al., 2022).Obtaining records of all these variables is difficult and expensive for cows, as opposed to heifers.Given that raising dairy heifers costs about 25% of total farm's production cost and feed consists of 50% of the total cost of raising heifers (Gabler et al., 2000;Akins, 2016), RFI of heifers is an economically relevant and logistically sensible trait to investigate.Some studies have reported that selection for RFI as heifers can carry on as mature cows (Nieuwhof et al., 1992;Macdonald et al., 2014).However, Connor et al. (2019) reported the phenotypic correlation of 0.37 between RFI of growing heifers and their later cow RFI in the United States.Efficient or low-RFI dairy heifers consumed 13 to 20% less feed (Williams et al., 2011;Connor et al., 2019) compared with their inefficient or high-RFI counterparts, with no effect on growth and performance.Further, selection for favorable RFI has been shown to decrease methane emissions in beef (Hegarty et al., 2007;Basarab et al., 2013) and dairy (Manzanilla-Pech et al., 2022) cattle.
The objectives of this study were to explore the genetic background of feed efficiency traits in Holstein heifers, estimate the genetic parameters of RFI and its component traits of Holstein heifers, and develop a system for genomic evaluation for RFI in Holstein dairy calves.

Experimental Animals and Design
No Animal Care Committee approval was necessary for the purposes of this study because all the information required was obtained from pre-existing databases.Data were collected from 6,563 growing Holstein heifers (initial BW = 261 ± 52 kg; initial age = 266 ± 42 d) across 182 trials conducted between 2014 and 2022 at the STgenetics Ohio Heifer Center (South Charleston, OH).This facility is a continuous-flow research dairy in which calves are managed according to usual farm practices from birth through their productive life.Upon reaching 6 to 8 mo of age, heifers were moved to the feed conversion testing facility and placed in warm-up pens, where they began adaptation to a corn silage-based test ration (NE G = 0.95 Mcal/kg DM; CP = 12.0%DM).During adaptation, heifers were evaluated for health and body size to establish groups of 40 to 64 animals that exhibited good health and adequate body size to freely consume feed from specialized feed bunks in the testing pens.Once a group was established, they were moved to a testing pen equipped with 8 electronic feed bunks (GrowSafe Systems Ltd., Airdrie, AB, Canada).Once heifers consumed the test ration for a minimum of 21 d, feed intake and performance were measured daily for a minimum of 70 d.Throughout the adaptation and testing period, heifers were offered ad libitum access to feed and clean drinking water.During each trial, BW was recorded using a chute weighing system (Tru-Test Inc., Mineral Wells, TX), either biweekly (every 2 wk) or twice during the first and last week of the testing period.For each animal, DMI was recorded daily during the testing period.

Computation of Traits
Individual animal feed intake was computed using a subroutine of the GrowSafe 4000E software (Process feed intakes) as described by Parsons et al. (2020).Dur- ing this study approximately 20% of the feed intake data was lost due to system failure.When data were lost, daily intake values were estimated by linear regression of DMI on day of trial using the Standard Least Squares procedure of JMP 16 (SAS Institute Inc., Cary, NC).Animal records were removed if the initial age of a heifer was less than 180 d or more than 400 d, if the feed intake was less than 2.27 kg for 10 or more days during the testing period.Additionally, animal records for metabolic BW (MBW) and ADG were removed if fewer than 3 usable body weights were recorded for an animal or if fewer than 40 d between the initial and final body weight measurements were available for an animal.Outliers were determined as values beyond 1.5 × interquartile range for DMI and BW.A linear regression of serial BW on day of trial was calculated for each animal using the Standard Least Squares procedure of JMP to determine the midpoint MBW (BW 0.75 , MBW = 59.5 ± 7.7 kg) and ADG (ADG = 1.09 ± 0.26 kg).Midpoint DMI was determined as the mean of all daily DMI over the measurement period.

Genotype Information
Animals from research herds and commercial herds were genotyped with different versions of commercial or STgenetics customized chips (GP4 and VM2; proprietary information).The contribution of GP4 (~28,000 SNPs) and VM2 (~61,000 SNPs) were 15% and 85%, respectively.Animals from different commercial herds that provided permission to STgenetics to use their data for research and development purposes were included in the analyses, and these animals form the part of prediction population (see Genomic Prediction).Quality control procedures were applied by removing sex chromosomes and within chip by removing the SNPs that had call rates of less than 95%, minor allele frequencies less than 0.05, and all animals with SNP call rates lower than 95%.The final number of SNP markers after quality control was 61,283.Genotypes were imputed in 3 steps.First, genotypes were pre-phased, and the missing markers were imputed using Eagle v2.4.1 (Loh et al., 2016). Second, shapeit v2.r904 (Delaneau et al., 2012) was used to convert the output from Eagle (HapMap) to VCF.Finally, minimac 4 1.0.2(Das et al., 2016) was used to perform the imputation.The imputed VCF files were then converted to plink files.The combination of software allowed for improved performance and accuracy of imputation, reaching 99.9%.

Calculation of Residual Feed Intake
Residual feed intake was computed within each trial as the difference between actual and expected DMI from the linear regression of mean DMI on MBW and ADG, as described by Koch et al. (1963): RFI = DMI − (µ + age + ADG + MBW), [1] where µ was overall mean effect for each trial, MBW was the midpoint metabolic body weight, ADG was the average daily gain, and age was the midpoint age of heifer fitted as covariate to account difference in DMI and MBW due to age of heifers.Here, each trial referred to the combination of date, barn, and pen.We estimated the average for the top 10% and bottom 10% ranked heifers to see the difference between high-and low-RFI heifers.

Estimation of Variance Components
Variance components of RFI component traits (DMI, MBW, and ADG) were estimated using DMU version 6 (Madsen and Jensen, 2013) according to the following model: where y ijk was the component traits of RFI for each animal, µ was overall mean effect, Trial i was the fixed effect of ith trial (182 levels), b 1 AGE jk was the linear regression on age of animal at midpoint, and animal k was the random additive genetic effect, which was assumed normally distributed with animal ~, , N a 0 2

Aσ ( )
where A was the pedigree relationship matrix built on pedigree traced back 10 generations and σ a 2 was the estimated additive genetic variance.e ijk was the vector of random residuals with e N e ~, , 0 2 Iσ ( ) where σ e 2 was the residual variance and I was the identity matrix.In the pedigree, the base population was considered unrelated, and unknown parent groups were not used.
The RFI was only modeled as the function of genetic effects (a) in addition to second residual (e) and overall mean (µ) to estimate the variance components because age and trial were already accounted for during estimation of RFI.

Genomic Prediction
The predictions of genomic EBV (GEBV) for the training population were obtained by solving the mixed model equations.For this, we fixed the variance components to the values obtained using pedigree relationship matrix (A) and used a genomic relationship matrix (G) to estimate the genomic breeding values.
Suppose that M represents the n g × n m matrix of n m genotypes for the n g genotyped animals, coded as 0, 1, and 2 for number of copies of reference or first allele for each genotype.The corresponding relationship matrix G was created according to VanRaden (2008).Here, where Π was a matrix with all the elements in column j containing the frequency of the second allele for SNP marker j for the base population.Once the GEBV of training population (animals with genotype and phenotype) were estimated, the GEBV of prediction population (animals with only genotype but not phenotype) were estimated as explained by VanRaden (2008).In brief, ˆˆ, where û V was the vector of breeding values of prediction population, û T was the vector of GEBV of the training population, G VT was the relationship between training and prediction population, and G TT −1 was the inverse of the genomic relationship of the training population.Reliability was calculated as follows: where C 22 was part of the inverse of the left-hand side matrix of the mixed-model equation (Taylor, 2014).
For validation, we selected 4 groups of prediction populations from commercial partner herds, each with 2,000 animals, in such a way that group 1 contains animals that share the same sires as those of the training population, group 2 contains animals that only share the same grandsire but not sire as those of the training population, group 3 contains animals that only share the same great-grandsire but not sire and grandsire as those of the training population, and group 4 contains animals that do not share the same sire, grandsire, and great-grandsire as those of the training population.

Estimating Standard Errors of Variance Components
The DMU software provided an average information matrix for the estimates of additive and residual effects.The asymptotic variance-covariance matrix of (co)variance estimates were based on the inverse of the average information matrix.As explained by Meyer and Houle (2013), 1,00,000 multivariate normal random vectors were drawn from a multivariate normal distribution with mean vector based on estimated variance components and variances based on the asymptotic variance-covariance matrix.Then, standard errors (SE) were determined from these random draws as standard deviation (SD) of the random draws.

Correlation of RFI Genomic PTA With Economically Relevant Traits
To ascertain the associations between the RFI and genomic PTA (gPTA) for the traits in the US national genetic evaluation for Holsteins (CDCB, 2023), approximated genetic correlations were calculated using the method of Calo et al. (1973), as described in Parker Gaddis et al. (2014): where ˆ, r g1 2 is the approximate genetic correlation between traits 1 and 2; REL 1i and REL 2i are reliabilities of (g)PTA for traits 1 and 2, respectively, for animal i; and r 1,2 is the Pearson (product-moment) correlation between gPTA for traits 1 and 2.

RESULTS AND DISCUSSION
The summaries of RFI, DMI, age, MBW, and ADG of all heifers are shown in Table 1.A histogram of phenotypic RFI for all calves (Figure 1) showed that RFI was normally distributed.The trait was centered at zero, with animals with lower RFI representing the most efficient animals because their feed intake was lower than expected based on their performance and body size.The observed phenotypic SD of RFI was 0.86 kg of DM per day, which was about 10% of average DMI.This variability was greater than that of Williams et al. (2011), who reported a phenotypic SD of 0.43 kg of DM per day.The greater variability in our data set compared with that of Williams et al. (2011) was due to the older calves in our study, which would increase the feed intake and hence the variability.The variability of RFI in our data set was lower than that reported for Holstein cows in the United States (Cavani et al., 2022; phenotypic SD of RFI of 1.99 kg of DM per day, about 7% of average DMI).The variation of RFI in growing heifers might be less compared with mature cows because more feed is consumed by cows than by heifers.Furthermore, RFI in cows is influenced by various factors such as pregnancy status, BCS, change in BCS, and more, which have not been accounted for in the model described by Cavani et al. (2022).Summaries of RFI and its component traits for the top and bottom 10% of heifers are shown in Table 2.The average RFI for the most efficient and least efficient heifers were −1.68 kg of DM per day and 1.54 kg of DM per day, respectively.This suggested a significant degree of variation in RFI, which could be used for genetic selec-  2), demonstrating the independence of RFI from growth performance and body size.We found a difference (P < 0.001) in DMI between the most and least efficient heifers, with efficient heifers consuming 33% less than their inefficient counterparts (6.70 vs. 10.04 kg/d, respectively).This was in concordance with Green et al. (2013) and Williams et al. (2011), who reported that the most efficient heifers ate less compared with the least efficient heifers.Similar results have been reported in beef cattle (Durunna et al., 2011;Elolimy et al., 2018).
The heritability estimate of RFI was 0.24 ± 0.02 (Table 3), which was comparable with the range of results (0.27 ± 0.12 to 0.33 ± 0.05) from other studies of Holstein heifers (Williams et al., 2011;Gonzalez-Recio et al., 2014).In beef heifers, Freetly et al. (2020) and Novo et al. (2021) reported RFI heritabilities of 0.25 ± 0.11 and 0.20 ± 0.08, respectively.The moderate RFI heritability measured in our study indicated that significant genetic variation existed within the commercial Holstein population and that genetic selection to improve feed conversion efficiency is possible.
The heritability estimate for DMI was 0.27 ± 0.02.This was within the range of results reported in the literature for Holstein heifers, ranging from 0.17 to 0.45 (Williams et al., 2011;Lin et al., 2013).The heritability of ADG was 0.19 ± 0.02.Our estimate was close to the estimate for Holstein heifers as reported by Williams et al. (2011) and lower than reported heritability estimates from beef heifers of 0.26 to 0.53 (Rolfe et al., 2011;Freetly et al., 2020;Esfandyari and Jensen, 2021).
The difference in heritability from beef heifers could be attributed to differences in the rations provided to dairy versus beef heifers (Freetly et al., 2020).Generally, beef heifers are fed a higher-energy and lower-fiber diet to achieve higher ADG compared with dairy heifers.The increased fiber content of the dairy heifer ration limits consumption as a result of increased gut fill compared with the lower-fiber and higher-concentrate ration fed to beef heifers.In addition to ration differences, the observed difference in heritability estimates could be due to differences in sample size, population characteristics (breeds, e.g., Holstein vs. Senepol; genetic diversity), and differences in trait measurement protocols (e.g., taking initial BW vs. midpoint BW).
Measuring RFI in growing heifers is complex and expensive, limiting the possibilities of implementing an industry-wide strategy to improve feed conversion efficiency based on phenotypic selection.Genomic selection is an appealing alternative.Genomic selection involves creating a reference population of individuals with phenotypes for the target trait (in this case, RFI),   genotyping these individuals for a panel of genomewide SNP, and then using this information to derive a prediction equation of the effects of SNP markers on the trait.This equation can then be used to predict breeding values for RFI in any genotyped animal (Meuwissen et al., 2001).Overall distribution of RFI gPTA and their reliabilities for training and different groups of prediction population are presented in Figure 2. The SD and range of gPTA of the training population were higher compared with that of different groups of the prediction population (0.27 kg vs. 0.23 kg and −0.9 to 0.75 kg vs. −0.82 to 0.72 kg; Figure 2A).A higher range of gPTA is preferred because it helps in segregation of high-and low-efficiency animals for selection purposes.Average reliability of the training population was 13% greater than the average reliability of the prediction population.The majority of the reliability fell below 50% for the prediction population (Figure 2B).Average reliability decreased from group 1 to group 4 as the relationship distance from the training population increased.The low reliability in group 4 was due to the poor connection with the training population, because these animals have smaller off-diagonal values in the G matrix with training population.Lower reliability of animals that were distant by age from the reference population was also reported by Li et al. (2020) in US Holstein bulls.Although genomic selection could be used to estimate gPTA of difficult-to-measure traits such as RFI, the size of the reference population and its relationship with the predicted population are the major determinants of prediction reliability.Table 4 shows the approximated genetic correlations of RFI and several key traits in the US national genetic evaluation provided by the Council of Dairy Cattle Breeding (CDCB, 2023): the lifetime net merit index ($NM), cow RFI (RFIc), productive life, milk yield, daughter pregnancy rate (DPR), SCS, fat yield, and protein yield.Overall, the approximated genetic correlations of RFI with the production traits were close to zero.Net merit, milk, fat yield, and protein yield had low positive correlations (0.06-0.14) with RFI, whereas DPR and SCS had low negative correlations (−0.05) with RFI.The low SE for those estimates is a function of the large number of samples used in this study; however, the biological significance of any of those estimates remained negligible.Weak genetic correlations of RFI with $NM, milk yield, fat yield, and protein yield in our study could be because our RFI was estimated from heifers, and we did not fit milk energy in the model to estimate RFI, whereas $NM included the cow RFI.The correlations between traits and RFI were highly affected by the energy sink traits that were fitted into the model for the estimation of RFI.In Australia, Gonzalez-Recio et al. ( 2014) reported weak genetic correlations of heifer RFI with milk yield, protein yield, and fat yield, albeit with high SE.We found nonsignificant genetic correlations of RFI with productive life, $NM, SCS, and DPR.The Pearson correlation of the breeding values for RFI heifers with RFIc breeding values was 0.17.Connor et al. (2019) reported a phenotypic correlation of 0.37 between RFI of growing heifers and their later cow RFI in the United States.However, the approximated genetic correlation of heifer RFI with RFIc was 0.43.It should be noted that approximated genetic correlations might be overestimated due to comparably lower reliabilities of heifer and cow RFI for animals used in Calo's formula (Calo et al., 1973).The moderate correlation of RFI for heifer and cow suggested that an index might be required to account for the relationship between RFI of cow and heifers, to select animals for lifetime production efficiencies.Genetic correlations between RFI and traits in the national genetic evaluation suggest that feed-efficient cows could be selected without affecting other major production traits.
Research has indicated that selection for low-RFI animals is a good strategy for reducing methane emissions in cattle, as significant relationships between methane production and RFI have been found, primarily as a function of the reduced DMI associated with low-RFI animals (Nkrumah et al., 2006;Hegarty et al., 2007;Johnson et al., 2019).With RFI having a low genetic correlation with other economically relevant traits, selection programs could use RFI in combination with other economically relevant traits to create highly productive progeny with increased feed conversion efficiency.More specifically, selection for heifers with genomic PTA of 1 SD (0.23 kg/d) below the mean could result in a feed cost savings of $20.48 per heifer based on average heifer consumption of 8 kg per day (Table 1) across heifer rearing period of 580 d, which was assumed to be 23 mo (age at first calving) minus the period from birth to the time (4 mo) when a heifer was fed growing heifer ration [$140 per US ton (907.20 kg) of DM].In addition to reducing overall feed cost, this would have direct effects on farm environmental sustainability, as animals that eat less while maintaining production have reduced manure and methane emissions per unit of product, thereby providing farmers with a solution to increase both profitability and sustainability (Nkrumah et al., 2004;Hafla et al., 2012;Basarab et al., 2013).

CONCLUSIONS
We explored the genetic basis of RFI and estimated the genetic parameters for RFI and its component traits for Holstein heifers.Our study suggests the existence of sufficient genetic variation for RFI among growing dairy animals to be used for selection.Genomic predictions of RFI for Holstein heifers provides new selection tools for feed efficiency improvement to commercial US dairy producers nationally.Future research should be directed toward integrating heifer RFI and cow RFI for selection of animals based on their overall lifetime production efficiency.

Figure 1 .
Figure 1.Distribution of phenotypic residual feed intake of US Holstein heifers.

Figure 2 .
Figure 2. Distribution of (A) genomic PTA (gPTA) and (B) their reliabilities of residual feed intake (RFI) for Holstein heifers on different groups of populations.The training group is the population with genotypes and phenotypes; group 1 contains animals that share the same sires as those of the training population; group 2 contains animals that only share the same grandsire but not sire as those of the training population; group 3 contains animals that only share the same great-grandsire but not sire and grandsire as those of the training population; and group 4 contains animals that do not share the same sire, grandsire, and great-grandsire as those of training population.The upper and lower edges of the box represent the third quartile and first quartile, respectively; the midline represents the median; the upper and lower whiskers represent values outside the middle 50%; and dots represent the potential outliers.
Khanal et al.: GENOMIC EVALUATION OF RESIDUAL FEED INTAKE tion among dairy heifers.As expected, no differences were observed in age, MBW, or growth rate among heifers with divergent RFI (Table

Table 1 .
Khanal et al.:GENOMIC EVALUATION OF RESIDUAL FEED INTAKE Summary statistics of residual feed intake (RFI) and its component traits at midpoint

Table 2 .
Characteristics [mean (SD)] of 10% of calves with lowest (most efficient) and highest (less efficient) residual feed intake (RFI) rankings

Table 4 .
Approximated genetic correlations of residual feed intake with US Council of Dairy Cattle Breeding genomic PTA for the traits in the national genetic evaluation(CDCB, 2023)