A comprehensive meta-analysis of genetic parameters for resilience and productivity indicator traits in Holstein cattle

Selection for resilience indicator ( RIND ) traits in Holstein cattle is becoming an important breeding objective as the worldwide population is expected to be exposed to increased environmental stressors due to both climate change and changing industry standards. However, genetic correlations between RIND and productivity indicator ( PIND ) traits, which are already being selected for and have the most economic value, are often unfavorable. As a result, it is necessary to fully understand these genetic relationships when incorporating novel traits into selection indices, so that informed decisions can be made to fully optimize selection for both groups of traits. In the past 2 decades, there have been many estimates of RINDs published in the literature, albeit in small populations. To provide valuable pooled summary estimates, a random-effects meta-analysis was conducted for heritability and ge-netic correlation estimates for PIND and RIND traits in worldwide Holstein cattle. In total, 926 heritability estimates for 9 PIND and 27 RIND traits, along with 362 estimates of genetic correlation (PIND x RIND traits) were collected. Resilience indicator traits were grouped into the following sub-groups: Metabolic Diseases, Hoof Health, Udder Health, Fertility, Heat Tolerance, and Other. Pooled estimates of heritability for PIND traits ranged from 0.201 ± 0.05 (energy corrected milk) to 0.377 ± 0.06 (protein content) while pooled estimates of heritability for RIND traits ranged from 0.032 ± 0.02 (incidence of lameness, incidence of milk fever) to 0.497 ± 0.05 (measures of body weight). Pooled estimates of genetic correlations ranged from −0.360 ± 0.25 (protein content vs. milk acetone concentration) to 0.535 ± 0.72 (measures of fat-to-protein ratio vs. milk acetone concentration). Additionally, out of 243 potential genetic correlations between PIND and RIND traits that could have been reported, only 40 had enough published estimates to implement the meta-analysis model. Our results confirmed that the interactions between PIND and RIND traits are complex, and all relationships should be evaluated when incorporating novel traits into selection indices. This study provides a valuable reference for breeders looking to incorporate RIND traits for Holstein cattle into selection indices.


INTRODUCTION
Genetic selection for increased production performance in dairy cattle has successfully improved multiple traits associated with production, such as lactation milk yield (MILK) and composition (Miglior et al., 2017).Additionally, the more recent introduction of genomics as a component of genetic selection has resulted in faster genetic progress for several breeds (Brito et al., 2021;Cole and VanRaden, 2018).When considering dairy cattle, Holstein cattle have become the most productive breed for milk production in intensive production systems (Brito et al., 2021).The high genetic gain for milk, fat, and protein yield that Holstein cattle have experienced is a major advantage for the dairy industry, as these traits have driven a large portion of the economic incentive (Cole and VanRaden, 2018).Subsequently, Holstein cattle have become the primary pure breed utilized by dairy producers around the world, especially in developed countries (Agriculture and Agri-Food Canada, 2021;CDCB, 2021).Moreover, the use of frozen semen and artificial insemination since the middle of the 20th century has allowed the dispersal of superior genetic material on an international scale (Brotherstone and Goddard, 2005).With the widespread adoption of the Holstein breed for dairy production and the spread of genetic material on a global scale, worldwide Holstein cattle have come to share a high degree of genetic relatedness (Norman and Powell, 1999;Sørensen et al., 2005;Hammami et al., 2007;Zenger et al., 2007;Danchin-Burge et al., 2011;García-Ruiz et al., 2015).
While selection for increased production has favorably affected multiple traits considered essential for the dairy industry, there have been unexpected side effects for other unconsidered traits (e.g., Rauw et al., 1998;Rauw et al., 2017;Brito et al., 2021).This is illustrated by the past selection that has been primarily focused on milk-related traits, which has been accompanied by a reduction in fertility rates (Berry et al., 2016;Fleming et al., 2019;Brito et al., 2021), higher incidence of metabolic diseases (Pryce et al., 2016), increased environmental sensitivity (Becker et al., 2020), and a relatively short productive life (De Vries and Marcondes, 2020;Souza et al., 2023).Further investigation has revealed that all these traits compete with production traits for energy allocation, and the primary selection for increased production has allowed a disproportionate allocation of energy resources toward production traits, leaving inadequate energy for unexpected stressors (Rauw et al., 1998).Therefore, modern populations of dairy cattle have become less resilient, especially in suboptimal environments with increased stressor exposure.In this context, resilience can be defined as the ability of an animal to cope with exposure to environmental stressors and efficiently return to pre-stressor homeostasis (Colditz and Hine, 2016).
Animal agriculture will face a myriad of changes in the near future as new social pressures and challenging climatic conditions force the refinement of breeding goals.These challenges are highlighting the need to select more resilient animals capable of maintaining adequate production levels in suboptimal conditions, as future production environments will expose animals to a considerable variation in environmental stressors (Rojas-Downing et al., 2017).To select for these important traits, we need to better understand their genetic background and relationship with traits already under selection.Currently, multiple publications have investigated resilience indicator (RIND) traits in Holstein cattle and published estimates of genetic parameters, but with a reduced sample size (Amin et al., 2000;Luo et al., 2021) or only a representation of particular Holstein subpopulations (Irano et al., 2014;Köck et al., 2019).In animal breeding, similar challenges have been overcome by utilizing a meta-analysis and have provided pooled genetic parameter estimates for traits within the Holstein breed (Zanton and Heinrichs, 2005;Hu et al., 2017;Wang et al., 2023a) and for indicator traits of resilience in other ruminant species (Mucha et al., 2022).However, no meta-analysis studies have focused on indicators of resilience and productivity in Holstein cattle.A meta-analysis is a statistical analysis of data from independent studies investigating the same question that produces a single, pooled (or weighted) estimate of effect, which can be used when designing or refining Holstein breeding programs (Gopalakrishnan and Ganeshkumar, 2013;Mikolajewicz and Komarova, 2019;Borenstein et al., 2021).Therefore, the primary objectives of this study were to: 1) verify if there is heterogeneity between studies and evaluate the factors that impact the differences in heritability (h 2 ) and genetic correlation estimates of productivity indicator (PIND) and RIND traits; 2) provide summaries of pooled h 2 and genetic correlation estimates based on a comprehensive meta-analysis for various PIND and RIND traits in worldwide Holstein cattle; and 3) present genetic parameters for some novel RIND traits with few reported estimates per trait.

Ethics Statement
No Animal Care Committee approval was necessary for the purposes of this study, as all information required was obtained from published studies.

Scope and Evaluated Traits
A literature search was conducted using Google Scholar (https: / / scholar .google.com)and PubMed (https: / / pubmed .ncbi.nlm.nih.gov) to identify published estimates of h 2 and genetic correlation for Holstein cattle.Two groups of traits were considered in this study: 1) PINDs and 2) RINDs (including welfare).Search terms used to conduct the literature search included ("heritability(ies)/genetic correlation(s)/genetic parameters"), (<name of the trait of interest > ), and ("Holstein/Holstein-Friesian/dairy cattle").The results generated from these search terms were then individually evaluated based on the following criteria: i) published between 2000 to 2022; ii) English language; iii) purebred Holstein (or Holstein-Friesian) cattle; and iv) provided single/average estimates (not a range) of h 2 / genetic correlation.Additionally, for estimates of the same parameter and trait derived from the same population but reported in multiple studies, only the most recent one was included in the analyses.To avoid bias, traits were only included in the meta-analysis if there were a minimum of 3 parameter estimates from at least 2 independent studies, a method recently utilized in a similar meta-analysis by Mucha et al. (2022).
After filtering the results obtained, a total of 986 h 2 estimates (270 PINDs and 716 RINDs) along with 362 genetic correlation estimates were retained for further analyses.These estimates came from 209 independent studies on 6 continents.Table 1 contains a full list of the traits included in the meta-analysis and their abbreviations.A larger number of traits were collected within the RIND group, so further sub-groups were created, including Metabolic Diseases, Hoof Health, Udder Health, Fertility, Heat Tolerance, and Other.A full listing of references along with genetic parameter estimates, standard error (SE), and number of records is available in Supplementary Table S1 (h 2 estimates) and Supplementary Table S2 (genetic correlation estimates) (https: / / doi .org/ 10 .6084/m9 .figshare.24239056.v2).Additionally, some RIND traits were evaluated but could not be incorporated within the framework of this analysis due to the small number of reported estimates per trait.Therefore, their genetic parameter estimates are presented in Table 5 for completeness.

Data Recorded and Quality Control
The final data set included h 2 and genetic correlation estimates along with their corresponding SE, if reported.Additionally, the following information was also recorded for each trait when available: publication year, name of publication venue, country of study, parity, number of records, type of quantitative trait (i.e., continuous or categorical), number of traits included in the statistical model (i.e., single-trait or multi-trait), variance component estimation method (i.e., Restricted Maximum Likelihood methods -REML or Bayesian methods), and the use of genomic data (Yes or No).Some studies did not report an associated SE for the genetic parameter estimates.In this case, there are 2 options to consider: either remove these estimates from the analyses or estimate SEs for these estimates.After consideration, the second option was chosen to retain the information provided by these estimates.Subsequently, approximate SEs were derived using the combined variance method (Sutton, 2000), which has been previously implemented in similar meta-analysis studies investigating genetic parameters for production animals (e.g., Akanno et al., 2013;Rojas de Oliveira et al., 2018;Hossein-Zadeh, 2021;Ladeira et al., 2021).The approximation equation is as follows: where SE ij is the predicted SE for the ith trait in the jth study that does not have a reported SE, s ik is the reported SE for the parameter estimate for the ith trait in the kth study with a reported SE, n ik is the sample size used to predict the reported parameter estimate for the ith trait in the kth study that has a reported SE, and n' ij is the sample size used to predict the published parameter estimate for the ith trait in the jth study that does not have a reported SE.
To overcome the absence of normal distribution, genetic correlation estimates were transformed to an approximate normal scale using Fisher's Z transformation (Steel and Torrie, 1960;Koots et al., 1994) as follows: , where r gij is the reported genetic correlation estimate for the ith trait in the jth study.Associated SEs were also transformed as follows: where n is the number of records used to calculate the SE.Following analyses, transformed parameter estimates were returned to the original scale using the following equation: Where r gij * is the back-transformed genetic correlation for the ith trait in the jth study, and Z ij is the Fisher's Z transformation.

Estimation of Pooled (Weighted) Means
A 3-level meta-analysis model was fitted for h 2 and genetic correlation estimates using the metafor package (Viechtbauer, 2010) implemented in the R software (R Core Team, 2021).In most cases, fitting a randomeffects meta-analysis model with only 2 levels would be prudent.However, this type of model assumes independence and normality of data, and there were multiple estimates included in the analyses originating from the same study and same population but being differentiated by factors such as parity, type of model, and day in lactation which were reported as separate traits.Knowing this violation of independence, we expanded the model to 3 levels to account for the dependence between the estimates.The 3-level model used for each individual trait (Harrer et al., 2021) is described as follows: where θij is the ith reported parameter estimate in the jth study (j is the unit of cluster), μ is the pooled population parameter mean, ζ ij is the within-cluster component of the deviation from the mean, assumed as ( ) where σ W 2 is within-cluster variance, ζ j is the between-cluster component of the deviation from the mean, assumed as ζ σ j B N , , ∼ 0 2 ( ) where σ B 2 is the between-cluster variance, and ϵ ij is the sampling error, assumed to be ϵ ij ~ N(0,v ij ).
The 95% lower and upper limits for the estimated parameters were computed as follows: . , and where LL μ is the 95% lower limit, UL μ is the 95% lower limit, μ is the pooled population parameter mean, and SE μ is the estimated SE of the pooled population parameter mean.

Heterogeneity Test
The I 2 index was used to quantify the proportion of total variation within-I W 2 ( ) and between-I B 2 ( ) cluster using the dmetar package (Harrer et al., 2019) implemented in the R software (R Core Team, 2021).The I 2 statistic was calculated as (Cheung, 2014): 2 and v are the estimated within-cluster heterogeneity, estimated between-cluster heterogeneity, and "typical" within-study variance, respectively.The I 2 estimates of quantified heterogeneity are reported on a scale of 0 to 100%.

Moderator Analyses
In an effort to explain and adjust for observed heterogeneity, moderator analyses were employed within the 3-level random-effects analyses for each trait.The objective of a moderator analysis is to find and account for variables (i.e., moderators) which influence the strength and/or direction of the relation between the independent and dependent variables (Harrer et al., 2021).In the current study the following moderators were investigated: country of study, parity, type of quantitative trait, number of traits included in the statistical model, variance component estimation method, and the use or not of genomic data.It is worth noting that the number of studies in each subgroup was small due to the reduced number of published estimates per moderator for many of the RIND and PIND traits included in the study, which can affect the estimate of the between-study heterogeneity across subgroups.

Pooled Estimates of Heritability
Pooled estimates of h 2 for PIND traits ranged from 0.201 ± 0.05 (lactation energy corrected milk) to 0.377 ± 0.06 (lactation protein content; PROT%), which may be considered moderately heritable (Table 2).Similar ranges of pooled estimates of dairy production traits have been reported in meta-analyses for other dairy species (Medrado et al., 2021;Mucha et al., 2022).For instance, the pooled estimate of h 2 for MILK in this study (0.225 ± 0.02) was similar to those found in Mucha et al. (2022) for dairy goats (0.27 ± 0.02) and dairy sheep (0.24 ± 0.02).Additionally, while models and methods to estimate traits such as DMI, ECM, and RFI were different across many of the studies, the heterogeneity observed for all of these traits was minimal (<0.00%;Table 2).
Pooled estimates of h 2 for RIND traits were generally low.This is consistent with most traits of health, fertility, and longevity due to their complex nature and the large impact of environmental effects on these traits.However, there were some traits: measures of electrical conductivity (EC; 0.381 ± 0.18), measures of stature (0.454 ± 0.02), and measures of body weight (0.497 ± 0.05), which were highly heritable (Table 3).Measures of EC may be highly heritable due to the reduced environmental noise during the phenotypic measurement of the trait, which is facilitated by the use of automatic milking systems and repeated records throughout the day and lactation (Norberg, 2005;Piwczyński et al., 2021;Pedrosa et al., 2023).The other 2 highly heritable traits, measures of stature and measures of body weight, are measurements of body conformation and are similar in magnitude to estimates of h 2 for other conformation traits (Němcová et al., 2011;Banos and Coffey, 2012;Khansefid et al., 2021).Excluding these 3 estimates, all other pooled estimates of h 2 for RIND traits (included in the meta-analysis) ranged from 0.032 ± 0.02 (incidence of milk fever/hypocalcemia) to 0.159 ± 0.04 (measures of acetone concentration; ACE; Table 3).
Additional observations can be elicited when examining individual subgroups of RIND traits.For traits included in the subgroup of Metabolic Diseases, some traits have already been incorporated in some Holstein selection indices worldwide, such as direct selection for incidence of clinical milk fever/hypocalcemia, incidence of clinical ketosis (KET), and incidence of clinical displaced abomasum (CRV, 2020;VanRaden et al., 2021;NAV, 2022).Additionally, there are examples of pooled estimates of h 2 for indicator traits, such as measures of milk β-hydroxybutyrate concentration (BHBA; 0.137 ± 0.02), and ACE (0.159 ± 0.04), both of which are indicator traits of KET (0.036 ± 0.02).The moderate heritability observed for these traits provides an alternate and potentially more effective route for genetic selection against KET which is not as heritable but is genetically correlated (van der Drift et al., 2012).As the dairy industry continues to focus increased resources toward selection for health traits, incorporating both types of traits (direct measures and indicators) into selection indices can help to reduce the incidence of metabolic diseases and increase overall resilience.
The next subgroup, Hoof Health, also contains important indicators of resilience that, if left unchecked, can cause negative welfare outcomes along with incurred extra costs and decreased profits for producers (Malchiodi et al., 2020).Additionally, the prevalence of hoof health problems in dairy cattle may motivate breeders to focus their attention on these traits (Chapinal et al., 2013;Malchiodi et al., 2017).While the pooled h2 estimates for these traits were low (no estimates of pooled h 2 exceeded 0.075; Table 3), direct genetic selection for improved hoof health is a viable avenue for long-term genetic improvement of these traits (Chapinal et al., 2013;Häggman and Juga, 2013).Thus, it may be beneficial to consider selecting for the most prevalent hoof health disorders, such as incidence of digital dermatitis, incidence of sole ulcer, or incidence of sole hemorrhage (Häggman and Juga, 2013;Malchiodi et al., 2017;Köck et al., 2019).
For Udder Health traits, an obvious target for genetic selection is incidence of clinical mastitis (MAST), which has severe indirect impacts on milk production and its derivatives (Heringstad et al., 2000;Martin et al., 2018).Similar to other disease traits, estimates of h 2 for direct selection on MAST in the past have mostly been lowly heritable (Miglior et al., 2015;Govignon-Gion et al., 2016;Lopes et al., 2020) with the pooled estimate being 0.044 ± 0.01 (Table 3).This may be due to the large impact of the environment and the complex nature of infection patterns with differences in infecting microorganisms, the severity and the duration of infection (Heringstad et al., 2000;Sender et al., 2013).Alternatively, indicator traits for MAST, such as measures of somatic cell score (0.153 ± 0.02) and EC (0.381 ± 0.18), represent traits that are comparatively more heritable, are specific measures of an immune response, and have the ability to identify subclinical cases of mastitis (Table 3).They represent valuable avenues that, when added to a selection index, can improve genetic gain when selecting for resistance to MAST.Specifically, EC is a highly heritable indicator trait that may merit further research as it can be easily collected through automatic milking systems and is highly genetically correlated to MAST (Norberg et al., 2006).However, it should also be considered that while the pooled h 2 estimate of EC was high, the SE was also relatively.Considering this, EC needs further investigation to fully understand its potential as an indicator trait of MAST.
The Fertility subgroup also represents a group of traits affecting welfare and profitability interests for the dairy industry, with consequences on rebreeding that might result in animal culling (Pritchard et al., 2013;Parker Gaddis et al., 2014).Due to the significant interactions between genetic and environmental factors for traits in this subgroup, h 2 estimates for these traits have been relatively low (no pooled estimates larger than 0.074) and continued research to find indicators that help to explain more of the genetic variation is needed (Walsh et al., 2011; Table 3).It should be recommended that, along with selection for improved resistance to these diseases, dairy producers should also work to improve the environments that their cattle are interacting with to realize the most improvement when selecting to increase cattle fertility.
Only one trait was included in the Heat Tolerance subgroup, which was measures of rectal temperature which had a pooled estimate of h 2 = 0.067 ± 0.04.This trait can be used as a direct indicator of heat stress and may provide valuable information when attempting to select animals that are less affected by elevated temperature exposure (Dikmen et al., 2012).While this was the only trait that could be included in the meta-analysis, there are multiple other novel indicators of heat tolerance that have been examined, especially based on reaction norm models for performance and reproduction traits (e.g., Ravagnolo and Misztal, 2000;   Shi et al., 2021).However, there have been relatively few studies that provide estimates of h 2 for other indicators of heat stress response which will be discussed later (Obioma et al., 2019;Luo et al., 2021).Traits included in the Longevity subgroup are a representation of a cow's ability to both stay healthy and productive (i.e., delay culling from the herd) over a measured period of time (Zavadilová and Štípková, 2012).All pooled estimates of heritability for traits in this subgroup were lowly to moderately heritable (Table 3).These traits have been of interest due to their economic significance (i.e., impact on profitability), MFEV = incidence of milk fever/hypocalcemia; KET = incidence of ketosis; BHBA = milk β-hydroxybutyrate concentration; ACE = milk acetone concentration; DA = incidence of displaced abomasum; DD = incidence of digital dermatitis; ID = incidence of interdigital dermatitis; HHE = incidence of heel horn erosion; IH = incidence of interdigital hyperplasia; WLD = incidence of white line disease; IP = incidence of interdigital phlegmon; SH = incidence of sole hemorrhage; SU = incidence of sole ulcer; LAME = incidence of lameness; MAST = incidence of mastitis; SCS = measures of somatic cell score; EC = measures of electrical conductivity; RP = incidence of retained placenta; MET = incidence of metritis; CO = incidence of cystic ovaries; ENDO = incidence of endometritis; CR = measures of conception rate; NRR = measures of nonreturn rate; RT = measures of rectal temperature; HLIFE = measures of herd life; PLIFE = measures of productive life; STAY = measures of stayability; PTUB = incidence of paratuberculosis; BRD = incidence of bovine respiratory disease; PERS = measures of persistency of lactation milk yield; BCS = measures of body condition score; STA = measures of stature; WT = measures of body weight.
but their impact on improved health can also indirectly improve overall cow resilience (Kašná et al., 2022).The final subgroup, Other, contains a variety of traits less related to the previously discussed trait groups but may be valuable for breeders when looking to holistically improve resilience in Holstein cattle.First, paratuberculosis (0.064 ± 0.02) and bovine respiratory disease (0.105 ± 0.12) are infectious diseases that have significant economic impacts on dairy production (Kirkpatrick and Lett, 2018;Quick et al., 2020).Estimates of pooled h 2 for both traits were low (Table 3).However, selection for resistance to either of these 2 diseases could benefit producers concerned about the impact of either disease and/or interested in increasing the longterm resistance to these diseases in their herds.Next, measures of lactation milk persistency has been a trait given focus as genetic improvement for this trait can result in both increased profitability and reduced stress experienced by dairy cows during lactation (Gengler, 1996).There were many studies providing estimates of h 2 for this trait (n = 75) and the pooled estimate of h 2 was low (0.120 ± 0.02).Additional discussion on the varied calculations used to define this trait will be mentioned further when addressing the heterogeneity of the estimates.Finally, measures of body condition score, measures of stature, and measures of body weight are all helpful measures that quantify body size.Knowing about these traits could be important in the context of heat tolerance, as there is evidence that larger and higher-producing animals tend to generate more metabolic heat (Kadzere et al., 2002).It may be useful to consider these traits, and maybe other indicators of metabolic heat production, when selecting for heat tolerance in dairy cattle.Body size traits are also related to reproductive performance and productive efficiency (Berry et al., 2003).

Weighted Estimates of Genetic Correlation
A heatmap of the numbers of estimates for all potential genetic correlations between PIND and RIND traits is presented in Figure 2.Many of the potential genetic correlations between traits in these 2 groups did not have a sufficient number of available estimates to conduct a meta-analysis (Figure 2).Immediately, this observation illustrates the need for more research that presents these correlations, as selection based on only single traits may inadvertently introduce unfavorable selection pressure on unconsidered traits (e.g., Brito et al., 2021).While researchers and producers may mainly be concerned with genetic interactions between PIND and RIND traits, it should also be considered that knowing the interactions between different RIND traits could also be of unique value (e.g., genetic cor-relations between KET and BHBA/ACE; Klein et al., 2020).This may allow breeders to focus their selection pressure on individual RINDs while incorporating progress for other RIND traits at the same time to produce animals that have better overall resilience.As a final point, the observed variation among favorable and unfavorable genetic correlations seen in this study elucidates the need to consider all potential interactions between selected traits if researchers want to make informed decisions when constructing selection indices (Table 4).
Pooled genetic correlation estimates for comparisons of PIND with RIND traits ranged from −0.361 ± 0.25 (PROT%_ACE) to 0.535 ± 0.72 (FPR_ACE), and some interesting trends were present when all pooled estimates of genetic correlation for individual RIND traits were compared (Table 4).First, PERS had a moderate, positive correlation with both MILK and FAT which signifies that favorable progress can be more easily made when one decides to include these traits in the same selection index.Second, WT was moderately and positively correlated with both FAT and DMI, and producers can obtain favorable genetic progress if looking to increase FAT or DMI through selection for WT.However, WT was lowly and negatively correlated with MILK and lactation energy corrected milk, which suggests that selection for WT will have a negligible effect on these traits.Next, KET, ACE, and BHBA were moderately and positively correlated with FAT% and FPR while being moderately and negatively correlated to PROT% (Table 4).This helps to illustrate the need to consider all potential genetic interactions as differing breeding goals (e.g., higher FAT% vs. higher PROT%) could result in different considerations when trying to realize favorable genetic progress for traits related to Metabolic Diseases and different PIND traits at the same time.Additionally, BHBA and ACE were lowly and negatively correlated to MILK which suggests negligible genetic correlations (Table 4).The pooled correlation between FPR and ACE should also be mentioned as FPR is classified as a PIND within the confines of this study.However, in many cases FPR is considered as an indicator of metabolic disease and the high correlations between it and other indicators show a reasoning behind its usage do to the ability to explain sizeable portions of the same genetic variation (Ranaraja et al., 2018;Klein et al., 2020).Finally, RT was positively and moderately correlated with MILK, FAT, and PROT.These correlations are unfavorable for animal breeders looking to increase production output of MILK, FAT, and PROT as their increase could results in increased body temperatures and subsequently increased incidence of heat stress, which negatively effects the production of the same traits (Dikmen et Maskal et al.: GENETICS OF RESILIENCE AND PRODUCTIVITY al., 2012;Luo et al., 2021).Furthermore, this elucidates the need to carefully consider selection decisions when trying to incorporate resilience traits (especially those related to heat tolerance) that require negative selection pressure for favorable genetic progress that would result in unfavorable progress in economically important traits.This finding suggests that using this direct indicator of heat tolerance may provide valuable information for creating heat tolerance indexes along with more studied indirect indicators.

Heterogeneity of Heritability Estimates
Within-and between-cluster heterogeneity for PIND traits were generally low, with the only values larger than 25% being I B 2 for FAT% (42.83%),PROT% (33.44%), and FPR (60.62%;Table 2).The increased heterogeneity observed for FAT% and PROT% may be due to a variety of factors, such as variation in the number of estimates of test-day records used to calculate lactation estimates of these traits (monthly vs. daily; Abdullahpor et al., 2013, Bastin et al., 2013), differences in breeding goals driving selection influenced by both location and time (Miglior et al., 2005), and the more recent implementation of selection for these traits in comparison to other yield traits (Miglior et al., 2017).The other pooled estimate with a high I B 2 was FPR, which may be influenced by variation among estimates of fat and protein yields among worldwide herds.However, further investigation of this analysis revealed that this pooled estimate appeared to be largely influenced by a single estimate originating from Jamrozik and Schaeffer (2012) (Figure 1).Within this study, they also observed a difference between models that expressed h 2 estimates for FPR on a daily basis (<0.4), which could not be incorporated in this metaanalysis as it was only reported as a range of estimates, and the h 2 estimate of FPR as an average of daily estimates, which was included in this meta-analysis (Jamrozik and Schaeffer, 2012).Within this study, it was hypothesized that these differences might be due to the inability of the daily estimate model to account for larger influences of short-term environmental effects (Jamrozik and Schaeffer, 2012).An additional analysis was done in which this single estimate from the original analysis (pooled estimate = 0.252 ± 0.05) was removed, and a new estimate of pooled h 2 from this resulting analysis was 0.198 ± 0.03 (95% CI: 0.131-0.264;I W 2 = 0.00%; I B 2 = 17.62%).Although the difference between both pooled estimates is not large, the estimate considering all the literature estimates (0.252 ± 0.05) is recommended to be used.

Overall, I W
2 and I B 2 for pooled estimates of RIND traits were generally low, with the only observed I B 2 greater than 25% being ACE (44.63%) and DA (32.61%;Table 3).As with the increased I B 2 observed for some PIND traits, this observation suggests variation in estimates between studies included in the meta-analysis.Additionally, it shows that this heterogeneity needs to be accounted for when conducting a meta-analysis, which in the case of this study has been done by utilizing a 3-level random-effects model to more conservatively estimate results (Sutton, 2000;Safari et al., 2005;Akanno et al., 2013).
There are multiple definitions of PERS as researchers have continually attempted to best define this trait (Atashi and Shahrbabak, 2006;Biassus et al., 2010;Otwinowska-Mindur et al., 2016).In this meta-analysis, all these definitions were considered a single trait in an attempt to best provide a representative estimate of PERS.The lack of heterogeneity observed among estimates (I W 2 = 0.00%; I B 2 = 18.97%) suggests that these methods provided a robust pooled estimate (Table 3).

Heterogeneity of Genetic Correlation Estimates
When examining the heterogeneity of genetic correlation estimates, there were many moderate and high values of I W 2 and I B 2 .With evidence for the need for a random-effects model being presented for meta-analyses of h 2 estimates in this study, these finding illustrate the need to implement a random-effects model for metaanalyses of genetic correlation estimates so that heterogeneity between estimates both with and between journals can properly be accounted for (Harrer et al., 2021).It is worth highlighting that the number of estimates was reduced, which can also contribute to less precise pooled estimates and greater heterogeneity.Additionally, accompanying SEs for many genetic correlations could be considered quite large.Knowing this, it will be beneficial for future researchers to calculate these estimates in larger populations to obtain more precise estimates of these genetic correlations.Ideally, having large data sets with records for all traits of interest measured in the same animals would be very beneficial for obtaining more complete and accurate estimates of genetic correlations among all trait pairs.

Moderator Analyses
When testing for moderators in h 2 meta-analyses, the only significant moderator effect that was detected was the type of statistical model for incidence of displaced abomasum.Results of this h 2 meta-analysis including this moderator are presented in the Supplementary Ta-  .v2).A larger number of traits had a significant moderator effect for genetic correlation meta-analyses with country being a significant moderator for 14 analyses and parity being a significant moderator for 7 analyses.The differences in estimates from countries might be due to differences in trait definitions across countries, recording procedures, and data collection protocols, which tend to be more standardized within each country.Results of genetic correlation meta-analyses, including these moderators are also presented in the Supplementary Table S4 (https: / / doi .org/ 10 .6084/m9 .figshare.24239056.v2).These differences in findings between pooled estimates of h 2 and genetic correlation are sensible when observing the larger number of pooled estimates with increased heterogeneity for genetic correlations which was previously discussed in this study.Additionally, the differences that are often seen between countries for breeding goals, the smaller sample size, genetic resources, trait definitions, and data collection protocols might explain the high heterogeneity observed for pooled estimates of genetic correlation.After more estimates of genetic correlations become available in the future, the analyses performed here could be revisited and the traits better separated based on their definitions and similarities.

Review of Other Novel Resilience Indicators
As phenotyping technologies continue to improve (Norton et al., 2019), novel traits are derived and greater emphasis is put on selection for traits related to resilience (Brito et al., 2020;Aquilani et al., 2022;Friggens et al., 2022).Multiple RIND could not be incorporated in the present meta-analysis due to a lack of estimates presented in Tables 5 and 6.Some interesting traits have been developed due to the abundance of data being obtained from automatic milking systems (Table 5).By investigating fluctuations in milk yield from calculated lactation curves, researchers have been able to develop traits based on these deviations (e.g., Poppe et al., 2020;Wang et al., 2022;Chen et al., 2023).A recent study by Wang et al. (2023b) quantified phenotypic variability in milk fluctuation periods related to diseases and explored milk fluctuation traits as indicators of disease resilience.All milk fluctuation traits evaluated were heritable with heritability estimates ranging from 0.01 to 0.10, and moderate to high genetic correlations with milk yield (0.34 to 0.64), milk loss throughout the lactation (0.22 to 0.97), and resilience (0.39 to 0.95), confirming the high impact of diseases on milk yield variability.Additionally, similar approaches have been recently implemented for data on daily step count to investigate fluctuations as well (Table 5).Finally, there are some metabolic diseases and indicators of hoof health that have been defined (e.g., Pryce et al., 2016;Malchiodi et al., 2017Malchiodi et al., , 2020)), however, there are only a few measures of genetic parameters for these traits available nowadays (Table 5).More emphasis will need to be allocated to these traits of their genetic correlations with other similar traits may be utilized if breeders want to emphasize resistance to these traits.Most RIND traits have been recorded in cows.However, there have been more recent studies evaluating RIND measured in calves, including survival from 3 to 60 d (h 2 = 0.088 ± 0.005), survival from 61 to 365 d (h 2 = 0.166 ± 0.006), survival from 366 d to the first calving (h 2 = 0.102 ± 0 0.006), calf diarrhea (h 2 = 0.048 ± 0.003), calf pneumonia (h 2 = 0.063 ± 0.004), and calf serum total protein (h 2 = 0.170 ± 0.019) (Zhang et al., 2022).
Climate change has received increasing attention in recent years, particularly in regard to its effects on  animal agriculture.Subsequently, researchers have begun to emphasize this in their research and have begun developing methods to select for improved tolerance to elevated temperatures (e.g., Nguyen et al., 2017).When discussing this subject, 2 main categories of indicators can be used to investigate heat tolerance, direct and indirect traits.Direct indicators are measurements of traits that are directly impacted by elevated temperatures, such as rectal temperature (Table 3) and respira-tion rate (Table 6) and can be used to provide accurate indicators of heat tolerance.In Holstein cattle, there are still relatively few heat tolerance indicators, while concern for heat stress continues to increase.By conducting more studies to obtain these estimates, we can start to better implement selection schemes to combat heat tolerance.Additionally, indirect indicators of heat tolerance have been proposed.These indicators are developed by investigating the effect of thermal heat index on common production traits that are already commonly collected (Table 7).This has the benefit of a lower implementation cost while still providing valuable information to breeders looking to genetically improve the heat tolerance of their dairy cattle.Both types of indicators of heat tolerance have advantages and disadvantages for making genetic progress when trying to improve heat tolerance, however, both provide valuable information that can be used to improve heat tolerance in Holstein cattle and merit further research to fully develop methods to implement them in breeding programs.

CONCLUSIONS
Dairy cattle breeding programs have successfully increased productivity over the past decades, especially in the Holstein breed.While the focus of past breeding programs has been primarily on traits related to PIND, it is also important to consider traits that are indicators of resilience.This is especially important as traits related to PIND often have antagonistic interactions with RIND traits.Pooled estimates of heritability for PIND traits ranged from 0.201 ± 0.05 (lactation energy corrected milk) to 0.377 ± 0.06 (PROT%) while pooled estimates of heritability for RIND traits ranged from 0.032 ± 0.02 (incidence of lameness, incidence of milk fever) to 0.497 ± 0.05 (measures of body weight).Pooled estimates of genetic correlations ranged from −0.360 ± 0.25 (PROT% vs. ACE) to 0.535 ± 0.72 (measures of fat-to-protein ratio vs. ACE).Knowing this information, it is necessary to pool the limited number of estimates of RIND traits and their genetic correlations with traditional PIND traits which are available in cur- rent literature so that researchers working to integrate new resilience traits into breeding programs can have more accurate estimates of these genetic parameters which otherwise are subject to many hurdles for accurate estimation with the current limited amount of data.Our results confirmed that the interactions between PIND and RIND traits are complex, and all relationships should be evaluated when incorporating novel traits into selection indices.

Figure 1 .
Figure 1.Forest plot for meta-analysis of fat-to-protein-ratio in Holstein cattle populations.
milk yield; FAT = lactation fat yield; PROT = lactation protein yield; FAT% = lactation fat content; PROT% = lactation protein content; DMI = measures of dry matter intake; FPR = measures of fat-to-protein ratio; ECM = lactation energy corrected milk; KET = incidence of ketosis; BHBA = milk β-hydroxybutyrate concentration; ACE = milk acetone concentration; DA = incidence of displaced abomasum; LAME = incidence of lameness; MAST = incidence of mastitis; SCS = measures of somatic cell score; RP = incidence of retained placenta; MET = incidence of metritis; PTUB = incidence of paratuberculosis; PERS = measures of persistency of lactation milk yield; STA = measures of stature; WT = measures of body weight; RT = measures of rectal temperature.

Figure 2 .
Figure 2. Heatmap of the number of estimates of genetic correlations between productivity and resilience indicator traits in Holstein cattle.
Maskal et al.: GENETICS OF RESILIENCE AND PRODUCTIVITY

Table 1 .
Maskal et al.:GENETICS OF RESILIENCE AND PRODUCTIVITY Description of the traits included in the meta-analysis and their abbreviations

Table 2 .
Pooled estimates of heritability for productivity indicator traits in Holstein cattle PROT% = lactation protein content; DMI = measures of dry matter intake; RFI = measures of residual feed intake; ECM = lactation energy corrected milk; FPR = measures of fat-to-protein ratio.

Table 3 .
Maskal et al.: GENETICS OF RESILIENCE AND PRODUCTIVITY Pooled estimates of heritability for resilience indicator traits in Holstein cattle

Table 4 .
Maskal et al.: GENETICS OF RESILIENCE AND PRODUCTIVITY Pooled estimates of genetic correlations for comparisons between productivity and resilience indicator traits in Holstein cattle

Table 6 .
Maskal et al.: GENETICS OF RESILIENCE AND PRODUCTIVITY Novel direct indicator traits of heat tolerance in Holstein cattle

Table 5 .
Heritability estimated for novel indicator traits of resilience in Holstein cattle