Emphasis on resilience in dairy cattle breeding: Possibilities and consequences

This study aimed to investigate dairy cattle breeding goals with more emphasis on resilience. We simulated the consequences of increasing weight on resilience indicators and an assumed true resilience trait (TR). Two environments with different breeding goals were simulated to represent the variability of production systems across Europe. Ten different scenarios were stochastically simulated in a so-called pseudogenomic simulation approach. We showed that many modern dairy cattle breeding goals most likely have negative genetic gain for TR and promising resilience indicators such as the log-transformed, daily deviation from the lactation curve (LnVAR). In addition, there were many ways of improving TR by increasing the breeding goal weight of different resilience indicators. The results showed that adding breeding goal weight to resilience indicators, such as body condition score and LnVAR, could reverse the negative trend observed for resilience indicators. Loss in the aggregate genotype calculated with only current breeding goal traits was 12 to 76%. This loss was mainly due to a reduction in genetic gain in milk production. We observed higher genetic gain in beef production, fertility, and udder health when breeding for more resilience, but from an economical point of view, this was not high enough to compensate for the reduction in genetic gain in milk production. The highest genetic gain in TR was obtained when adding the highest breeding goal weight to LnVAR or TR, both with 0.29 genetic standard deviation units. The indicators we used, body condition score and LnVAR, can be measured on a large scale today with relatively cheap methods, which is crucial if we want to improve these traits through breeding. Economic values for resilience have to be estimated to find the most optimal breeding goal for a more resilient dairy cow in the future.


INTRODUCTION
To date, no dairy cattle breeding goals (BG) have been developed or implemented that include resilience as a breeding goal trait.Scheffer et al. (2018) defined resilience as a risk of system failure ranging from organs to entire organisms.This definition will require new dynamic indicators of resilience, such as through wearable electronics that monitor the status of the organism more frequently than before.It has been shown that animals differ in the ability to cope with environmental disturbances (Friggens et al., 2016;Revilla et al., 2019).A resilient dairy cow can avoid being affected by or quickly recover from environmental disturbances and return to its normal status (Berghof et al., 2019;Poppe et al., 2020).These disturbances can be caused by changes in quantity or composition of feed intake, pathogens, or heat waves.If we want to breed for a more resilient dairy cow, we need ways of measuring resilience.To date, it is difficult to define and measure resilience in a suitable way for breeding.Because disturbances can be of different nature, resilience is a composite trait (Berghof et al., 2019).Several resilience indicators in dairy cattle have been proposed.Body condition score reflects the cow's energy balance and is related to disease resistance and tolerance (Calus and Veerkamp, 2003;König and May, 2019).Log-transformed, daily deviation from the lactation curve (LnVAR) has been proposed as a valuable resilience indicator (Elgersma et al., 2018;Poppe et al., 2020).In addition, the ability of an animal to consume DM has been proposed as a resilience indicator (van Dixhoorn et al., 2018).Resilience indicators, such as LnVAR and BCS, have been shown to have an unfavorable genetic correlation with milk production (Kadarmideen, 2004;Poppe et al., 2020), whereas LnVAR has favorable genetic correlations to health, fertility, and longevity (Poppe et al., 2020).Furthermore, Elgersma et al. (2018) suggested that a resilient dairy cow needs less attention from the farmer because it had, for example, better health and fertility.Hence, it is expected that a BG aiming to improve resilience will emphasize functional traits and, together with resilience indicators, will improve true resilience.
Holstein cattle dominate the global dairy industry, and it is also the most common dairy cattle breed in Europe (FAO, 2021).Holstein cattle are used in different production systems across Europe, from more intensive indoor production systems to more extensive systems utilizing the grazing and implementing seasonal calving.Furthermore, BG for Holstein cattle are different across Europe due to the different importance of traits.Genotype by environment interaction (GxE) in breeding programs has been an issue in animal breeding because it often reduces the response to selection when multiple environments are part of the breeding program.When breeding for more resilient cows, GxE has been proposed to be valuable, because it is a source of genetic variation for adaptation to the environment (Mulder, 2016).Slagboom et al. (2019) showed that genomic selection improves the possibility of applying multiple breeding schemes in various environments.Hence, breeding collaborations across Europe could be more important when breeding for more resilient cows.Likewise, the possibilities for this are increased with the implementation of genomic selection.
This study aimed to investigate dairy cattle BG with more emphasis on resilience.We simulated the consequences of increasing weight on resilience indicators and our definition of true resilience within 2 environments.We hypothesized that we can breed for more resilient animals through resilience indicators, but this will have a cost in terms of achieved genetic gain for traits in current BG.

MATERIALS AND METHODS
In this study, genetic progress of resilience in dairy cattle was studied in different scenarios that included GxE interactions and differed in breeding goal weights for resilience indicator traits and true resilience.Two environments with different BG were simulated to represent the variability of production systems across Europe.The first environment reflected an intensive indoor production system, referred to as the West Atlantic environment (WA).The second environment reflected a grass-based extensive production system, referred to as the central mountainous environment (CM).

Breeding Goals
An overview of traits included in this study can be found in Table 1.All traits were defined so their preferred direction has a positive economic value.These traits were included because they are of economic importance or were expected to be of importance to a BG that has the aim of increasing resilience in dairy cattle.In addition, a trait was included that represented true resilience in dairy cows (TR).Our basic scenario, which mimicked current WA and CM BGs, included milk production (MPR) and beef production (BP) as production traits and DMI, fertility (FERT), and udder health (UH) as functional traits.All 3 functional traits have been suggested as important traits for resilience (Elgersma et al., 2018;van Dixhoorn et al., 2018).From now on, MPR, BP, DMI, FERT, and UH are referred to as current breeding goal traits.We are aware that modern dairy cattle BG often include more traits, such as the Nordic total merit index (NTM; Sørensen et al., 2018).The main reason for not including all traits in the NTM was that the simulations would require too much computer power and take too long time to run.A BG was set up for each subpopulation and each environment; in the WA environment, the BG (WA-BG) was more focused on improving milk production, and in the CM environment the BG (CM-BG) was more focused on improving functional traits (Table 2).Breeding goal weights were initially based on economic values calculated for usage in the NTM for Holstein cows (Sørensen et al., 2018).Because the BG used in this study include a smaller number of functional traits than the NTM, the BG weights were adjusted so that the correlation between milk production and the total BG matched current BG in Denmark, Finland, and Sweden.By setting the correlations between MPR and the overall breeding goal at a predefined level we ensured that the balance between production and functional traits were at a realistic level.The correlation between each trait and the total BG was calculated was calculated using equation [1]: where r i,j = the genetic correlation between trait i and BG j (j = WA-BG, CM-BG); b i = a vector containing a 1 for trait i and 0 for the other traits; b j = a vector containing the BG weights for BG j; and G = a matrix containing the genetic correlations between BG traits (Table 3).
The correlation between milk production and WA-BG matched that of the yield index and the NTM for conventional dairy production (0.69; Nordic Cattle Genetic Evaluation, 2020).The correlation between milk production and CM-BG matched that of the yield index and the NTM for organic dairy production (0.40; Sørensen et al., 2018).No weight was given to BCS, LnVAR, and TR in either environment in the basic scenario (Table 2), but BG weights on these traits increased gradually in the alternative scenarios.In the scenario named BCS1, we used the BG weights in the basic scenario but changed the BG weight for BCS to Trait: milk production (MPR), beef production (BP), fertility (FERT), udder health (UH), log-transformed variance of deviations from lactation curve (LnVAR), true resilience (TR).All traits were defined so their preferred direction has a positive economic value (less DMI = positive economic value, less LnVAR = positive economic value).References: Van Arendonk et al., 1991;Dematawewa and Berger, 1998;Pryce et al., 2001;Berry et al., 2002;Søndergaard et al., 2002;Calus and Veerkamp, 2003;Kadarmideen, 2004;Buch and Norberg, 2008;Weller and Ezra, 2008;Vallimont et al., 2010;Pabiou et al., 2011;de Haas et al., 2015;Lu et al., 2015;Li et al., 2016Li et al., , 2020;;Hardie et al., 2017;Team Avlsvaerdivurdering, 2019;Poppe et al., 2020;Interbull, 2021. 3 Genetic SD are in trait units: ECM production per milking day (MPR), growth of bull and heifer calves in g/d (BP), DMI in kg/d, pregnancy rate in cows in % (FERT), incidence of clinical mastitis per lactation (UH), BCS on a scale of 1-5, log-transformed variance of deviation from the lactation curve (LnVAR), and the genetic SD of TR was set to 1.00.
50 and kept the weights on LnVAR and TR at zero.In scenario BCS2, the BG weight for BCS was increased to 100, and in BCS3 it was further increased to 150.In a similar way we increased the BG weight of LnVAR to 50, 100, and 150 in the scenarios named LnVAR (1-3) and we increased the BG weight of TR in scenarios TR (1-3).

Genetic Parameters
A genetic and a phenotypic correlation matrix between traits within one environment were constructed from literature (Table 3) to be used in both WA and CM.Residual correlations were calculated based on genetic and phenotypic correlations and variances.Missing correlations were based on expert opinions gathered within the GenTORE project of the European Union's Horizon 2020 Research and Innovation Program.In some cases, missing phenotypic correlations were assumed to be the same as genetic correlations.
The genetic correlation between a trait that was measured in both environments was used as a measure for GxE interaction (Table 3; r g WA-CM).Genotype by environment interactions were defined via correlations between the 2 environments and were based on Interbull GxE estimates between Denmark, Finland, and Sweden, and Ireland for the traits MPR, FERT, and UH (Interbull, 2021).We assumed the r g for DMI to be 0.6 in our simulated environments (Berry et al., 2014;de Haas et al., 2015).The r g for BP measured in the 2 environments was set to 0.80 (Pabiou et al., 2011); GxE estimates for BCS, LnVAR, and the TR were not available.The r g for BCS, LnVAR, and TR was set to 0.70 based on expert discussions and a mean of r g for other important traits for resilience, DMI, UH, and FERT.
Genetic correlations between a trait that is measured in one environment and another trait that is measured in the other environment were calculated using equation [2] (Buch et al., 2009): where r i1,j2 = r g between trait i in environment 1 and trait j in environment 2; r i1,i2 = r g between trait i in environment 1and environment 2 (r g WA-CM); r i,j = r g between trait i and trait j within an environment; r i1,j2 = r g between trait j in environment 1and environment 2 (r g WA-CM).
In some scenarios the genetic or residual correlation matrix was not positive-definite and thus it was bent using the method by Hayes and Hill (1981).The correlations shown in Table 3 are not bent.

True Resilience
True resilience was defined and developed based on expert opinions within the GenTORE innovation program.Although there is currently no consensus on a satisfactory measurement of true resilience, it was defined in the simulation as a measurable trait, and its meaning was described by the genetic correlation of this trait to each of the other traits.Arguments for each correlation with the true resilience were as follows: • Milk production had an unfavorable and moderate r g with true resilience.The r g was lower than the r g for MPR and LnVAR because LnVAR is based on milk data and therefore more directly linked to MPR than true resilience.• Beef production had a low and favorable r g with true resilience due to higher body reserves leading to higher resilience.• Due to the absence of any prior information, the correlation between DMI and TR was estimated based on the genetic covariance matrix by using the softImpute package in R statistical software (Hastie and Mazumder, 2015).The correlation that was found by this method was moderate and favorable (i.e., higher DMI -higher TR).• Fertility and body condition score can be interpreted as indirect, long-term indicators of resilience.A highly resilient cow will have enough body reserves to allocate to fertility.In situations where feed is restricted due to environmental disturbances or illness for example, a higher body condition score means a cow has reserves to use to respond to these environmental disturbances.The correlations of FERT and BCS to TR were therefore favorable and moderate.• Udder health can be interpreted as a direct, shortterm indicator of resilience.Sickness occurs more easily when resilience is low.This genetic correlation was therefore higher than the correlation with FERT and BCS.• The LnVAR was believed to be the best indicator for true resilience out of the traits included in this study; therefore, the r g was favorable and high.• Phenotypic correlations between TR and the other traits were assumed to be the same as the genetic correlations.

Simulation
The 10 different scenarios were simulated in a socalled pseudogenomic simulation in stochastic simulation program ADAM (Pedersen et al., 2009).This means that for every true breeding value (TBV) of a trait in the simulation, a correlated trait was added that represented the direct genomic value (DGV) for that trait (Dekkers, 2007;Buch et al., 2012).The genetic correlation between the DGV and the TBV represented the accuracy of prediction using DGV for that specific trait.Thus, no QTL or markers were simulated, and we set a fixed accuracy of genomic predictions in each year of the simulation.This accuracy was calculated using the deterministic method by Goddard (2009) and adjusted for multiple information sources by Buch et al. (2012) (Table 3).A historical effective population size of 750 was used to calculate DGV accuracies, along with a genome length of 30 Morgans and 38,000 informative markers.The reference population size was assumed to be 180,000 for MPR based on the EuroGenomics reference population (Koivula et al., 2021).For other traits we were more conservative with 85,000 animals in the reference population except for DMI, which was set to 10,000 based on available feed intake data that is assumed to be available in the near future.The DGV accuracies ranged from 0.47 (DMI) to 0.54 (MPR, BP).
The breeding scheme was the same in all scenarios (Figure 1).Each environment consisted of 10,000 cows in the breeding nucleus in the age of 1 to 5. Every year, 2,000 bull calves and 2,000 heifer calves were selected within each environment for genotyping.This (truncation) selection step was based on parent average EBV and included only one-year-old animals.Out of the genotyped heifers, the 200 best were selected as donors for multiple ovulation and embryo transfer.Each donor was mated to 2 different bulls and each mating resulted in 3 offspring giving a total of 6 offspring per donor.Out of the 2,000 genotyped bull calves, the 100 best were selected for breeding to all cows in the breeding nucleus and to all donors.Bulls were selected in 2 selection steps; one for each environment and corresponding breeding goal.Selected bulls could originate from both environments.Every year, 15% of all animals were randomly culled, as well as unselected males.Every scenario was simulated for 30 years and replicated 30 times.Breeding value estimation was performed in DMU (Madsen and Jensen, 2013).All cows had phenotypic observations for the traits in their own environment, as each trait was simulated separately for every environment.Most traits were realized at the age of 3 after cows had been selected for breeding, except for BP; animals had phenotypic observations for this trait at the age of 1.

Statistical Analysis
All statistical analysis were performed in R statistical software (https: / / r -project .org).Annual genetic gain was calculated in genetic standard deviation (σ A ) units for every trait in every environment and scenario by using genetic levels in the last 10 yr of the simulation.Genetic gain in the aggregate genotype (H) was calculated for the current breeding goal traits (Table 1; i.e., all traits excluding BCS, LnVAR, and TR).Genetic gain in σ A units was multiplied with the BG weight in Table 2 for each current BG trait and this was then all summed together to calculate the genetic gain in H.
The genetic gain in H in the alternative scenario compared with the basic scenario was calculated by using equation [3]: An ANOVA test was performed to find significant differences in genetic gain per trait and environment between all scenarios.A significant result of the ANO-VA test was followed by Tukey's honest significant difference test to find all pairwise differences, by use of the Agricolae package in R (de Mendiburu, 2020).In addition, significant differences in inbreeding per generation between scenarios were calculated in the same way.We allowed bulls to be selected across WA and CM environment, and a 95% confidence interval was calculated for the average proportion of sires per year selected by the other environment.

Genetic Gain in Resilience
In both environments, genetic gain in BCS, LnVAR, and TR significantly improved in all alternative scenarios [BCS (1-3), LnVAR (1-3), and TR (1-3)] compared with the basic scenario.Genetic gain in BCS, LnVAR, and TR was negative in the basic scenario in environment WA (Figure 2).The highest genetic gain in BCS, LnVAR, and TR was observed when adding the most BG weight on the respective trait.Positive genetic gain was seen in TR and BCS in scenarios BCS (1-3), Ln-VAR (1-3), and TR (1-3).The highest genetic gain in TR was observed in scenario TR3 and LnVAR3 both with 0.29 σ A units (Figure 2).The highest genetic gain in BCS, LnVAR, and TR was observed when adding the most BG weight on the respective trait.Genetic gain in BCS, LnVAR, and TR was larger in all scenarios in the CM environment compared with genetic gain in these traits in the same scenario in the WA environment (Figure 3).In the CM environment, the basic scenario achieved a positive genetic gain in TR, but a negative genetic gain in BCS and LnVAR.In addition, all scenarios with a BG weight on BCS, LnVAR, or TR resulted in a genetic gain in BCS, LnVAR, and TR.

Genetic Gain in the Current Breeding Goal Traits
Genetic gain for the current BG traits and for H were calculated for the WA environment (Table 4) and for the CM environment (Table 5).The largest gain for H was observed in the basic scenario in both environments.A significantly lower genetic gain in H was seen in all alternative scenarios [BCS (1-3), LnVAR (1-3), and TR (1-3)] compared with the basic scenario.The differences between scenarios for genetic gain in H were large, especially for the scenarios with the most  2).In the other scenarios, basic scenario breeding goal weights and breeding goal weights (50/100/150) for BCS = BCS (1-3) or LnVAR = LnVAR (1-3) or TR = TR (1-3).The SD were in the range 0.01-0.02.emphasis on resilience and resilience indicators (BCS3, LnVAR3, TR3).However, this decrease was less profound in the CM environment compared with the WA environment.Some clear differences in current breeding goal traits were also found between scenarios.The largest genetic gain in MPR was seen in the basic scenario for both environments.Compared with the basic scenario, a significantly larger genetic gain was seen for BP, FERT, and UH in scenarios BCS (1-3), LnVAR (1-3), or TR (1-3).Genetic gain in DMI was negative in all scenarios in the WA environment (Table 4), but genetic gain in DMI became slightly positive in the LnVAR3 scenario in the CM environment (Table 5).The CM environment had more genetic gain in the functional traits and beef production and less genetic gain in milk production compared with the WA environment.In general, shifting to alternative scenarios yielded more gain in functional traits compared with the basic scenario.2).In the other scenarios, basic scenario breeding goal weights and breeding goal weights (50/100/150) for BCS = BCS (1-3) or LnVAR = LnVAR (1-3) or TR = TR (1-3).The SD were in the range 0.01-0.02.Genetic gain in aggregate genotype was calculated for the current breeding goal traits (milk production, beef production, DMI, fertility, udder health).

Bulls Selected Across Environments
The strongest correlation between the BG in the WA and the CM environment was observed in the basic scenario (0.81) and the weakest in LnVAR3 (0.70).At equilibrium (year 21-30), the average proportion of sires per year selected from the other environment was the highest in the basic scenario, 0.12 (0.03, 0.21) selected by WA environment and 0.02 (0.00, 0.05) selected by CM environment.Conversely, the lowest proportion of bulls selected from the other environment was found in the LnVAR3 scenario, 0.04 (0.00, 0.08) selected by the WA environment and 0.01 (0.00, 0.03) selected by the CM environment.
There were no significant differences in the rate of inbreeding between scenarios, because the same breeding scheme was used.The rate of inbreeding per generation in the WA environment ranged from 0.09 to 0.10%, and in the CM environment 0.12 to 0.15%.

DISCUSSION
We explored genetic gain in resilience in 2 different environments by adding different BG weights to resilience and the resilience indicators BCS and LnVAR.The results showed that adding weight to resilience or resilience indicators in the BG would adversely affect genetic gain in H and MPR.We saw a significantly higher genetic gain in true resilience and resilience indicators, and due to favorable genetic correlations, higher genetic gain in BP, FERT, and UH.From an economical point of view the increased gain in functional traits was however not high enough to compensate the reduction in genetic gain in MPR.

Breeding for Better Resilience
In this study, dairy cattle BG with more emphasis on resilience were studied.The results showed that increasing the economic weight on true resilience in the breeding goal improved true resilience the most, as was to be expected (scenario TR3, Figures 2 and  3).The BG with additional weight on the resilience indicators show what practically can be implemented when breeding for more resilience.The results showed that BCS and LnVAR can be used to improve TR, where scenario LnVAR3 yielded the highest gain in TR (Figures 2 and 3).The reason for this is the high and favorable correlation between LnVAR and TR.The LnVAR is a trait that is easily measured from daily milk yield data and is thus readily available from most modern milking equipment.As was predicted by the genetic correlations of this traits by Poppe et al. (2020), LnVAR is well suitable as a potential selection criterion when aiming to improve resilience.
Another important result of the current study is the effect of breeding for improved resilience on genetic gain in H.This was calculated with the economic weights of the current breeding goal traits and thus the improvement in resilience did not have an economic value.In scenario LnVAR3 the loss in H was high; 61% in the CM environment and 76% in the WA environment.Such a loss in H is unlikely to be accepted by any breeding company for implementation in practical 1 SD were in the range 0.01-0.02.
2 Genetic gain in aggregate genotype was calculated for the current breeding goal traits (milk production, beef production, DMI, fertility, udder health).
3 Basic = scenario with breeding goal weights for milk production, beef production, DMI, fertility, udder health (Table 2).In the other scenarios, basic scenario breeding goal weights and breeding goal weights (50/100/150) for BCS = BCS (1-3) or log-transformed variance of daily deviation from the lactation curve (LnVAR) = LnVAR (1-3) or our defined true resilience (TR) = TR (1-3). 4 where ∆G H change , = change in genetic gain in H in percentages compared with the basic scenario, and ∆G H basic , = genetic gain in H in the basic scenario.
breeding programs, and therefore these scenarios are mostly useful for scientific purposes.In scenario Ln-VAR1, the loss in H was only 13 to 14% compared with the basic scenario.This scenario could be a starting point for improving resilience without losing too much in H. Furthermore, an alternative approach to avoid or even achieve some genetic improvement in resilience could be to adjust the weights of the current BG traits using a desired gains approach (Kariuki et al., 2017).

Indicator Traits for Resilience
In addition to BCS and LnVAR, some of the current breeding goal traits have also been proposed as resilience indicators.For example, FERT and UH are resilience indicators that we measure and often breed for today in dairy cattle breeding.However, we observed a negative genetic gain in TR in the basic scenario even if FERT and UH were included based on their economic values.Thus, to improve resilience through these indicator traits, the BG weight of FERT and UH needs to be increased to a higher weight than what is calculated based on present economic models.Dry matter intake has also been proposed as a resilience indicator (van Dixhoorn et al., 2018).However, current BG aim for lower DMI and higher MPR.When breeding for resilience, the opposite is desirable.
Longevity has also been proposed as a resilience indicator (Schuster et al., 1992).However, other authors have argued that it is better to improve functional traits such as health and fertility than longevity (Schuster et al., 1992;Schmidtmann et al., 2021).Cows with greater production capacity and better health and fertility are better at avoiding involuntary culling and are therefore more valuable for farmers.We did not include all important functional traits due to limited computer power.However, by setting the correlations between MPR and the overall breeding goal at a predefined level, we ensured that the balance between production and functional traits were at a realistic level.In this way, UH can be seen as a proxy for multiple health traits.Thus, there will still be genetic improvement of longevity through the direct selection of FERT, health, and resilience.

Definition of TR
We developed the TR trait with experts in resilience research within the GenTORE innovation program.True resilience was defined as a measurable trait, and its meaning was here described through the genetic correlation of this trait to each of the other traits.The aim was to reflect the overall genetic gain of resilience and include this as a trait in the breeding goal.We used a heritability of 0.10 for TR, similar to Berghof et al. (2019).One of the most used definitions of resilience is that a cow can avoid being affected by or quickly recover from environmental disturbances and return to its normal status (Berghof et al., 2019;Poppe et al., 2020).The results of this study show that the definition of TR that we have used will increase the ability of a cow to avoid or quickly recover from environmental disturbances.For example, genetic gain in TR, BCS, and FERT changed from negative in the basic scenario to positive in scenario TR1 in the WA environment (Table 4 and Figure 2).Furthermore, genetic gain in UH changed from 0.03 to 0.12, whereas the genetic level in LnVAR remained stable.Thus, adding breeding goal weight on TR will change the breeding goal in a direction that is in line with the definition for resilience in previous literature.

Difference in Environments
The 2 environments, WA and CM, were chosen to represent 2 mainstream production environments in Europe today.Our results show that the WA environment had a negative genetic trend for BCS, LnVAR, and TR when ignoring resilience in the breeding goal.This is due to the high economic value for milk production and the unfavorable genetic correlations between MPR and resilience traits.The CM environment had a negative genetic trend for BCS and LnVAR but a slightly positive genetic trend for TR.The BG in this environment was more focused toward improving functional traits and less focused on milk compared with the BG in the CM environment.Due to favorable genetic correlations between functional traits and resilience, the CM environment showed a better genetic trend for the resilience indicators.Thus, choosing a more functional BG compared with a more dairy-focused BG will positively influence resilience, even without considering resilience as a BG trait.

Resilience and Efficiency
Another trait that is important to look at in dairy cattle breeding is feed efficiency.The traditional way of defining feed efficiency is by calculating the ratio of kilograms of milk produced per kilogram of DMI consumed.The genetic trend for DMI was negative (i.e., higher feed intake) in all scenarios in this study, except in scenario LnVAR3 in the CM environment.At the same time, we observed less genetic gain or even negative genetic gain in MPR when breeding for more resilience.Hence, if feed efficiency would be calculated as the ratio between MPR and DMI, one could argue that we are breeding for less feed efficient dairy cows in the scenarios where DMI and MPR both show a decrease in genetic gain.However, we see more genetic gain in BCS even in the alternative scenarios with no weight on BCS.The antagonism between MPR and BCS can be better understood in the light of the resource allocation theory.More resources can be allocated to body reserves if there is less demand of resources for MPR.When fewer resources are needed for MPR, more resources can be used to start a new gestation or cope with disturbances.Therefore, we argue that efficiency should be seen over a lifetime and that a resilient animal is likely to have a longer productive lifespan than a less resilient animal.Mechanistic models linking the main aptitudes related to resilience and efficiency in the dairy cow are promising tools to better understand the selection response expected on the different components leading to higher lifetime lactation efficiency (Puillet et al., 2021).Additionally, we observed more genetic gain in BP when breeding for more resilient animals.Indeed, BP was measured for growing animals whereas all other traits were measured at later physiological stages.This could explain why genetic correlations between BP and other traits are only low to moderate.The genetic correlation was also weaker between DMI and BP than DMI and BCS.More research is necessary to better understand the relationship between efficiency and resilience.

Collaboration Between Breeding Schemes
We allowed bulls to be selected across WA and CM environment.The correlation between the BG in the WA and the CM environment decreased when giving BCS, LnVAR, and TR weight in the BG.In general, we saw minor usage of bulls between environments, the correlation between the 2 environments was below the split-point correlation (Mulder and Bijma, 2006).Cao et al. (2020) found the split-point correlation to be 0.85 to 0.875 in a genomic selection breeding scheme, and our observed correlations were weaker (0.70-0.81), which explains why we only saw minor usage of bulls between environments.
We used the same breeding scheme for both environments and only minor differences in the rate of inbreeding per generation was observed.The rate of inbreeding per generation (0.09-0.15%) was low and far below recommended 1% rate of inbreeding (FAO, 2013).We believe that it is reasonable to have 2 breeding schemes with the size used in our study, because Holstein is also the most common dairy cattle breed in Europe.

Implementation
The results of our study can be implemented by breeding organizations.We believe that it is unlikely that breeding organizations or farmers would adopt a breeding goal with a negative genetic gain in MPR, which would implicate a large loss in H.However, there were scenarios where we achieved genetic gain in MPR and resilience traits simultaneously.These scenarios are more suitable to implement in practice.However, for a proper conclusion as to which breeding goal is the most optimal in terms of genetic gain in H, economic values for resilience have to be estimated.
New resilience indicators are likely to become available in the future with the current development of big data (Berghof et al., 2019).In addition, promising resilience indicators, such as LnVAR, could be combined and improved with more daily data, such as recovery from clinical mastitis (Adriaens et al., 2021).The breakeven price for recording indicators depends on the value of the traits we want to improve (Hansen Axelsson et al., 2015).Both BCS and LnVAR are traits that can be measured on many farms today.The LnVAR can be measured through modern milking equipment such as automatic milking systems (Poppe et al., 2020), and BCS can be measured through camera equipment or traditional conformation scoring (Fischer et al., 2015;Nir et al., 2018).Hence, it should be possible to obtain a relatively large amount of phenotypes for BCS and LnVAR without high costs.New methods for large-scale measuring of DMI are also becoming available, such as camera equipment (Lassen et al., 2018;Bezen et al., 2020).It is more challenging to measure DMI in a grazing-based system such as the CM environment in our study, but it may be possible by using drones (Herlin et al., 2021).In general, we agree with Berghof et al. (2019) that the current development of big data is very positive in determining new resilience indicators.Thus, the future looks promising in terms of improving resilience in dairy cattle breeding.

CONCLUSIONS
We showed that many modern dairy cattle BG most likely have negative genetic gain for our definition of TR and promising resilience indicators such as the logtransformed, daily deviation from the lactation curve (LnVAR).Further, there are many ways of improving resilience in dairy cattle breeding by putting more weight on different resilience indicators.The indicators we used, BCS and LnVAR, can be measured on a large scale today with relatively cheap methods, which is crucial if we intend to improve these traits through breeding.Adding weight to resilience indicators could reverse the negative trend for resilience indicators and our definition of true resilience, but with a cost in genetic gain in milk production.Economic values for resilience have to be estimated to find the most optimal breeding goal for a more resilient dairy cow in the future. 2 Figure 1.Schematic overview of the simulated breeding scheme.MOET = multiple ovulation and embryo transfer.
Figure 3. Genetic gain per year in σ A units for BCS Log-transform variance of daily deviation from the lactation curve (LnVAR), and true resilience (TR) in each scenario in the central mountainous environment.Basic = scenario with breeding goal weights for milk production, beef production, DMI, fertility, udder health (Table2).In the other scenarios, basic scenario breeding goal weights and breeding goal weights (50/100/150) for BCS = BCS (1-3) or LnVAR = LnVAR (1-3) or TR = TR (1-3).The SD were in the range 0.01-0.02.

Table 1 .
Bengtsson et al.: EMPHASIS ON RESILIENCE IN DAIRY CATTLE BREEDING Traits included in this study

Table 2 .
Bengtsson et al.: EMPHASIS ON RESILIENCE IN DAIRY CATTLE BREEDING Breeding goal 1 weights and correlations with the whole breeding goal for all traits in the basic scenario 2

Table 3 .
Genetic Bengtsson et al.: EMPHASIS ON RESILIENCE IN DAIRY CATTLE BREEDING

Table 4 .
Genetic gain 1 in σ A units and change in genetic gain (%) in the aggregate genotype 2 (H) compared with the basic scenario ∆G H change , 2

Table 5 .
Bengtsson et al.: EMPHASIS ON RESILIENCE IN DAIRY CATTLE BREEDING Genetic gain 1 in σ A units and change in genetic gain (%) in the aggregate genotype 2 (H) compared with the basic scenario ∆G H change , Bengtsson et al.: EMPHASIS ON RESILIENCE IN DAIRY CATTLE BREEDING