Genetic evaluation of heat tolerance in Holsteins using test-day production records and NASA POWER weather data

Heat stress is a prominent issue in livestock production, even for intensively housed dairy herds in Canada. Production records and meteorological data can be combined to assess heat tolerance in dairy cattle. The overall aim of this study was to evaluate the possibility of genetic evaluation for heat tolerance in Canadian dairy cattle. The 2 specific objectives were (1) to estimate the genetic parameters for milk, fat, and protein yield for Holsteins while accounting for high environmental heat loads, and (2) to determine if a genotype-by-environment interaction causes reranking of top-ranked sires between environments with low and high heat loads. A repeatability test-day model with a heat stress function was used to evaluate the genetic merit for milk, fat, and protein yield under heat stress and at thermal comfort for first parity in 5 regions in Canada. The heat stress function for each trait was defined using a specific temperature-humidity index (THI) threshold. The purpose of this function was to quantify the level of heat stress that was experienced by the dairy cattle. The estimated genetic correlation between the general additive genetic effect and the additive effect on the slope of the change in the trait phenotype for milk, fat, and protein yield ranged from −0.16 to −0.30, −0.20 to −0.44, and −0.28 to −0.42, respectively. These negative correlations imply that there is an antagonistic relationship between sensitivity to heat stress and level of production. The heritabilities for milk, fat, and protein yield at 15 units above the THI threshold ranged from 0.15 to 0.27, 0.11 to 0.15, and 0.11 to 0.15, respectively. Finally, the rank correlations between the breeding values from a repeatability model with no heat stress effect and the breeding values accounting for heat stress for the 100 top-ranked bulls indicated possible interaction between milk production traits and THI, resulting in substantial reranking of the top-ranked sires in Canada, especially for milk yield.


INTRODUCTION
An environment with high heat loads can negatively affect the reproduction, production, and other physiological functions of dairy cattle.Heat stress can also result in substantial economic losses for the dairy industry (St-Pierre et al., 2003).In addition, it is likely that a 0.4°C increase in the global mean surface temperature will occur over the next 2 decades.Moreover, there is evidence that heat waves will become more intense, last for longer durations, and occur at higher frequencies, and that near-surface humidity and evaporation over land will increase (Hoegh-Guldberg et al., 2018).This could result in an increase in the occurrence and severity of heat stress in dairy cattle, amplifying the negative effects associated with high heat loads.
Management strategies, such as the implementation of sprinklers and fans as well as modifying feed nutrition, are effective tools at diminishing heat stress (Baumgard and Rhoads, 2009;Nickerson, 2014).However, genetic selection provides certain advantages over management interventions.For instance, an animal with a high tolerance to heat stress may have inherited a greater capacity to grow, thrive, or maintain fitness despite being exposed to an unfavorable environmental factor and, therefore, would rely less on expensive and labor-intensive management interventions.The additive genetic component of heat tolerance could be passed down to succeeding generations resulting in positive and cumulative genetic gains in thermotolerance (Simms, 2000).Different species, breeds, and individuals have diverse abilities for coping with high heat loads.Their thermotolerance can vary due to individual or herd factors, ability to acclimatize, or a genetic adaptation (Collier and Collier, 2012).
In terms of production in dairy cattle, a heat tolerance cow can be defined as a cow that experiences a slower decline in production traits in response to increasing heat loads relative to a conspecific that is more sensitive to heat stress (Nguyen et al., 2018).Most heat tolerance studies in dairy cattle use the temperaturehumidity index (THI) to define the heat load in an animal's environment.Reaction norm (RN) models can be used to model changes in an individual's performance across a varying environmental parameter, such as THI to quantify their phenotypic plasticity (Hammami, 2009).A genotype-by-environment interaction (G×E) occurs when 2 or more genotypes respond differently to environmental variation, such as changes in THI (Hayes et al., 2016).A plastic animal would be more affected by changes in the environment, whereas a more robust animal would maintain production across environment conditions, in other words, the RN model would estimate stronger slopes across the varying environments for plastic animals (Hayes et al., 2016).This modeling technique makes it possible to evaluate the genetic merit of individuals under heat stress.
A significant G×E can cause the reranking of genotypes or heterogeneity of variances across environments.The reranking of genotypes is problematic because individuals with superior performance may have severely reduced productivity in a different environment (Hammami et al., 2015;Wakchaure et al., 2016).A significant G×E can complicate a breeding program by decreasing the accuracy of breeding values, reducing selection responses, and hindering the economic gain for a trait.Even though there is not strong evidence for the reranking of genotypes within countries due to heat stress (Carabaño et al., 2014;Hayes et al., 2016), the possibility of a significant G×E due to heat stress across all the dairy cattle production regions in Canada was never investigated.
Several previous studies have used weather station data to estimate the genetic parameters for heat tolerance in dairy cattle.However, ground-based weather station data sets can be severely limited by numerous factors, such as inconsistent data collection, low spatial resolution, and daily and yearly data set gaps.An alternative weather data resource is NASA POWER (https: / / power .larc.nasa.gov/), which provides meteorological data over regions where surface measurements are sparse or nonexistent.Rockett et al. (2022) showed that THI values based on weather data retrieved from NASA POWER were highly correlated with weather station based THI values for dairy herd locations in Canada.The data from the NASA POWER database are advantageous because they are available for any dairy herd location and requires less editing, especially with respect to missing or partial recording (Rockett et al., 2022).Therefore, the overall aim of this study was to evaluate the possibility of genetic evaluation for improved heat tolerance in the Canadian Holstein population using NASA POWER weather data to allow for broader coverage of herds and regions of dairy production in Canada, which are limited when using weather station data.The 2 specific objectives of this study were to (1) estimate the genetic parameters for 3 production traits in Holsteins for 5 regions in Canada while accounting for heat stress, and (2) determine if there is significant reranking of EBV for the top-ranked bulls at high and low THI values.

MATERIALS AND METHODS
Because existing data sets were used in this research, approval by an Institutional Animal Care and Use Committee was not necessary.The initial data set for this investigation had 21.7 million test-day records of milk yield (g/d), fat yield (g/d), and protein yield (g/d) measured between 2010 and 2019 from 1.2 million cows in approximately 8,500 herds located across Canada.The production data were provided by the Canadian Dairy Network (CDN, a member of Lactanet Canada) and was grouped into 5 regional subsets.The regions were classified as Ontario, Quebec, British Columbia, the Prairies, and the Atlantic Maritime.The Prairies region included the provinces Saskatchewan, Manitoba, and Alberta, and the Atlantic Maritime region included the provinces Newfoundland and Labrador, Nova Scotia, Prince Edward Island, and New Brunswick.These 2 regions were defined based on climatic similarities between provinces.The data set was split into different regions for this investigation because Canada does not have a uniform climate.Every herd was required to have known geographical coordinates and records from at least 5 different animals per year for at least 5 yr, while each record was required to be collected between 5 and 305 DIM.Only records from first lactation animals were included in this analysis.Animals were required to have at least a total of 5 test-day records in the parity one with at least one trait measurement recorded at a THI below the estimated heat stress threshold and at least 2 records above the same threshold.The number of records and cows removed based on these restrictions were between 222,100 and 810,370, and 26,283 and 105,796, respectively.The number of records, animals, and herds used in the analyses, as well as the averages and standard errors for milk, fat, and protein yields for each region are listed in Table 1.
The addresses for each dairy herd were converted to latitude and longitude coordinates using Google Maps Geocoding (https: / / cloud .google.com/blog/ products/ maps -platform/ address -geocoding -in -google -maps -apis).The daily minimum relative humidity (RH, %) and maximum ambient temperature (AT, °C) for each test-day record from 2010 to 2019 were retrieved from NASA POWER based on each herd's geographical coordinates.Weather data from NASA POWER are freely available and can be downloaded through their data access viewer or by including a POWER API URL in a scripting application.A THI value was calculated for the test-day and the 2 d prior for each production record.The 3 THI values were then averaged into a single value and assigned to a variable called THI 2d .This value was used to account for the possible delayed effect of heat stress caused by the thermal environment of the days before the test-day (Ravagnolo and Misztal, 2000).The equation that was used to calculate the THI values was THI = [1.8× (AT) + 32] -[0.55 -0.0055 To estimate the THI thresholds for dairy cattle across Canada, all regional production data sets were combined into a single data set.This data set was then adjusted using a linear model for known sources of environmental and genetic variation using the same method as Rockett et al. (2022).The model included fixed effects of classes of age at calving and DIM, milking frequency and herd-year, and random cow additive genetic effect.The adjusted phenotypes were averaged across herds within each THI unit starting at THI 30 and plotted against THI.The thresholds at which the production traits began to increase or decrease at a different rate were initially estimated by visual assessment.These thresholds were then used to define the segmented polynomial regression models that were fitted to the data set for each production trait.The THI values above and below the initially selected thresholds were also used to define the best fitting segmented polynomial knots based on maximization of the coefficient of determination (R 2 ) value.
A repeatability test-day RN model was then used to estimate the genetic parameters for the 3 production traits records from first lactation animals while accounting for heat stress for each individual region using ASReml software (Gilmour et al., 2015).The equation for the repeatability test-day model was where y ijklm is the mth test-day record, which was either the milk, fat, or protein yield, for the lth cow, µ is the overall mean, AD i is the fixed effect of the ith combination of age at calving and DIM classes, M j is the fixed effect of milking frequency (j = 1 to 3 classes), HY k is the random effect of the kth combination of the herd and year, β is the coefficient of the fixed regression on T m , and T m is the THI 2d value for the mth test-day record, a l is the general random additive genetic effect of the lth cow, p l is the general random permanent environment effect of the lth cow, α l is the random additive genetic effect on the slope of the change in the phenotype of the lth cow, π l is the random permanent environment effect on the slope of the change in the phenotype of the lth cow, f(m) is the heat stress covariate for the mth test-day record (as defined below), and e ijklm is the random error term.For age at calving, there were 8 classes (<24, 25, 26, 27, 28, 29, 30, >31 mo) defined for the first parity.There were 11 classes defined for DIM (5-29, 30-60, 61-90, 91-120, 121-150, 151-180, 181-210, 211-240, 241-270, 271-300, 301-305 d).The heat stress covariate was defined as , where f(m) is the heat stress covariate for the mth testday record, THI 2d is the average THI for the test-day and the 2 previous days, and THI t is the THI threshold value for the specific production trait.The random terms for the above repeatability model were assumed to have a mean of zero and follow a normal distribution.
The variance-covariance structure for this model was , where I is an identity matrix, A is the numerator relationship matrix, σ a 2 and σ ∝ 2 are the variances of the general additive genetic effects and the additive genetic effects on the slopes, σ p 2 and σ π 2 are the variances of the general permanent environment effects and the permanent environment effects on the slopes, σ a∝ is the covariance between the general additive genetic effects and additive genetic effects on the slopes, σ π p is the covariance between the general permanent environment effects and permanent environment effects on the slopes, σ h 2 is the herd-year effect variance, and σ e 2 is the error variance.A large random sample of approximately 75% of the herds was taken from each regional data set to estimate the variance components using ASReml software (Gilmour et al., 2015).The number of animals in each subset are shown in Table 4.The entire data sets were used to estimate the breeding values accounting for heat stress using the above repeatability model.Heritability (h 2 ) estimates for the production traits at f(m) points of THI 2d above the threshold were calculated using .
The genetic correlations between the general additive genetics (a) and the additive effects on the slope ∝ ( ) of the change in the phenotype were calculated using Breeding values without the heat stress component will be referred to as the traditional breeding values.These were calculated for all 3 production traits using another repeatability model without the effect of THI using AS-Reml software (Gilmour et al., 2015).The equation for this repeatability test-day model was All the terms in this model have already been defined above.The EBVs from this repeatability model without the heat stress function will be referred to as EBV_TM with TM referring to the traditional model.
The EBV from the repeatability model with the heat stress function for animals under heat stress or in other words experiencing a THI above their threshold will be referred to as EBV_HS (= a + f(m) × α), with HS meaning heat stress.The EBV from the repeatability model with the heat stress function for animals within their thermal comfort zone or in other words experiencing a THI below their threshold will be referred to as EBV_CZ (= a) with CZ referring to a cow's comfort or thermoneutral zone.Spearman correlations between the EBV_TM and EBV_CZ, as well as EBV_TM and EBV_HS for each trait were then calculated twice, once using all bulls with a reliability of 70% or greater and once using the top 100 ranked bulls with 70% or greater reliability.Reliability was calculated based on the prediction error variance from the traditional model for each trait, using the method proposed by Van Vleck (1993):

RESULTS AND DISCUSSION
Livestock breeding programs have become increasingly more interested in traits such as robustness, environmental sensitivity, resilience, and environmental flexibility.Although there are several different definitions for these terms in the literature, robustness and environmental sensitivity are commonly defined using an RN model.A more robust or less environmental sensitive genotype in terms of heat stress is one with a close to zero slope in response to an increasing heat load or a higher threshold that marks the late onset of heat stress relative to other individuals (Ravagnolo and Misztal, 2000;Hayes et al., 2016).Due to this functionality, a repeatability model fitting a linear RN was used to investigate the environmental sensitivity of Canadian dairy cattle to heat stress in this study.
An RN can be defined by 2 parameters: the thermal or THI threshold, and the subsequent decline in the trait of interest.The THI values were used to quantify the heat load within a dairy cow's environment.For simplicity, a universal threshold was assumed for all animals.The THI thresholds at which milk, fat, and protein yield in Canadian dairy cattle started to change at a different rate per unit THI are presented in Table 2.The relationships between THI and each production trait are also shown in Figure 1.There were 3 thresholds identified for milk yield, one threshold identified for fat yield, and 2 thresholds identified for protein yield.The thresholds that were chosen to mark the onset of heat stress and consequently define the heat stress functions for milk, fat, and protein yield were 66, 55, and 53 THI, respectively.These thresholds were selected because they marked the lowest THI at which the average rate of change in the adjusted production traits per unit of THI was negative.The third and second threshold for milk and protein, respectively, marked the THI at which production traits started to decline at a faster rate per unit of THI, on average.However, this model assumed that each production trait started to decline linearly after the selected THI threshold.This was a simple, but adequate, approximation of the actual pattern of change in these production traits as THI increases.However, there was some loss of information while assuming this static universal threshold and linear decay.The average THI for spring, summer, and fall between 2010 and 2019 for each region are shown in Table 3.This shows that the average THI during this time period is above at least one of the traits THI thresholds in at least one region of Canada for all 3 seasons examined.
The genetic parameters were examined at 15 units above the respective thresholds for each trait to ad- *The thresholds with an asterisk were used to define the heat stress function for each trait.THI1, THI2, and THI3 are the first, second, and third temperature-humidity index thresholds, respectively.equately study the full effect of heat stress in the summer months in Canada.At 15 units of THI above the threshold, Holstein cows are under a level of HS where production traits decrease substantially with a corresponding THI value for milk solids (fat and protein) that would happen in about 100 d in the provinces with most dairy cattle in Canada (Campos et al., 2022).The THI values, 81, 70, and 68, were 15 units above the respective thresholds for milk, fat, and protein yield, respectively.The estimates for the (co)variance components at comfort zone and 15 units above the THI thresholds for milk, fat, and protein yield are shown in Table 4.The permanent environment variances of the slope were generally higher than the additive genetic variances of the slope for all traits, but with some exceptions.This may indicate that the sensitivity to heat stress is to a great extent acquired instead of inherited.
Descriptive statistics for the EBV of bulls are presented in Table 5 for milk, fat, and protein yield.The mean EBV for bulls with daughters under heat stress (EBV_HS) evaluated at 15 THI units above the threshold was lower than the mean EBV for bulls with daughters within their thermal comfort zone (EBV_CZ) and the mean EBV from the traditional repeatability model (EBV_TM) for all traits.The EBV_CZ describes a cow's genetic potential for production within their thermal comfort zone.The thermal comfort zone refers to THI values below the respective THI threshold, whereas the EBV_TM describes their genetic potential for production without accounting for heat stress.The mean EBV_HS for each trait was negative for all regions, except for milk yield in Quebec and Prairies.This indicates that heat stress has a negative impact on the genetic potential for production in Canadian Holsteins.Additionally, all regions had similar EBV_HS, EBV_TM, and EBV_CZ means and ranges with some exceptions.Milk yield Quebec 1 σ a 2 = general additive genetic variance; σ p 2 = general permanent environment variance; σ ∝ 2 = additive genetic variance of the slope; σ π 2 = permanent environment variance of the slope; σ a ∝ = additive genetic covariance; σ π p = permanent environment covariance all evaluated at 15 units above the temperature-humidity index threshold (i.e., at 15 × α and 15 × π), where α and π are the slope of the additive genetic effect and permanent environmental effect above the threshold.Table 5.The mean and range of the EBV for bulls from the traditional repeatability model without the heat stress effect (TM) and from the repeatability model with the heat stress function at a temperature-humidity index (THI) value within the thermal comfort zone (CZ) and at 15 THI units above the heat stress threshold (HS) for milk, fat, and protein yields (g/d) Heritability estimates at 15 units above the THI threshold are shown in Table 6.These heritabilities for milk, fat, and protein yield ranged from 0.15 to 0.27, 0.11 to 0.15, and 0.11 to 0.15, respectively.Heritabilities were similar across regions and traits.These low to moderate heritability estimates were similar to the results reported by other studies using a similar repeatability model with a heat stress function.For instance, Ravagnolo and Misztal (2000) found that the heritability for milk and protein yield ranged between 0.16 and 0.21, and Bernabucci et al. (2014) found that heritability for milk, fat, and protein yield ranged between 0.07 and 0.23.Additionally, Aguilar et al. (2009) reported heritabilities for the same 3 traits ranging from 0.10 to 0.24 using a test-day random regression model with linear splines.The heritabilities at the THI threshold from the model that included THI effect and the heritabilities estimated using the variances from the model without THI effect are shown in Table 6.These heritabilities were similar for all 3 production traits with no visible trends.

Rockett et al.: GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINS
The genetic correlations between the general additive genetic effects and the additive genetic effects for the slope are shown in Table 7.The estimated genetic correlations for milk, fat, and protein yield ranged from −0.16 to −0.30, −0.20 to −0.44, and −0.11 to −0.15, respectively.Therefore, there is a low to moderate antagonistic relationship between the decline in production due to heat stress and the general level of production.This means that sustained selection for production traits without consideration for heat stress will likely continue to decrease the ability of dairy cattle to cope with heat stress.Ravagnolo and Misztal (2000) also found a moderate antagonistic genetic correlation between general production and less sensitivity to heat stress ranging from −0.23 to −0.35 for milk, fat, and protein yield.Similarly, Bernabucci et al. (2014) found that genetic correlations between general and less sensitivity to heat stress ranged from −0.24 to −0.56 for 4 production traits.This negative relationship is most likely due to the highly active metabolism that is associated with high production.A highly active metabolism produces a substantial amount of body heat, which can increase an animal's core body temperature in a hot environment if their heat dissipation processes are not able to compensate for the internal and external heat gained (Collier and Collier, 2012).
Table 7 also presents the genetic correlations between the heat stress estimated breeding values (EBV_HS) and general additive genetic effect (a), and between EBV_HS and additive genetic effect on slope (α) both evaluated at 15 units above the THI threshold (i.e., EBV_HS = a +15 × α) for milk, fat, and protein yield.In general, EBV_HS was more strongly correlated with general additive genetic effects (r ~0.79) than with additive genetic effects on the slope (r ~0.19) for both fat and protein yields, but the difference was less pronounced for milk yield (r ~0.65 and 0.50, respectively), with the Atlantic Maritime and Prairies regions showing a higher correlation between EBV_HS and additive genetic effects on the slope compared with the general additive genetic effects (r = 0.71 vs. 0.52 and r = 0.65 vs. 0.59, respectively), which indicates a possible stronger G×E interaction for milk yield.
Genotype-by-environmental interaction can be grouped into 3 categories.This includes no interaction, noncrossover interaction, and crossover interaction (Baye et al., 2011).There is no interaction when one genotype constantly outperforms another genotype by the same quantity across different environments.A noncrossover interaction is when one genotype outperforms another genotype across a varying environment, but the difference in performance between genotypes is not consistent, but does not result in the reranking of genotypes.A crossover interaction occurs when the difference in performance between genotypes across different environments is not consistent causing a reranking of genotypes (Baye et al., 2011).A crossover interaction or in other words, significant reranking reduces the selection response of a trait when animals are reared in a different environment from the selection environment.
A scaling effect on variances can also occur due to G×E and can negatively affect genetic gain if not properly accounted for in genetic evaluations.If a scaling effect is not accounted for, animals evaluated in an environment with lower genetic variance will be selected less often than animals evaluated in an environment with a higher genetic variance (Kolmodin, 2003).Spearman rank correlations between the EBV_TM and EBV_CZ, as well as EBV_HS at 15 units above the THI threshold for milk, fat, and protein yield are presented in Table 8.The rank correlations in Table 8 were calculated from all the bulls with a reliability of 70% or greater based on the traditional repeatability model with no heat stress effect.Generally, a genetic correlation between 2 trait measurements in different environments that is less than 0.80 is an indicator of a G×E that can cause substantial reranking of genotypes (Hayes et al., 2016).The results of the current study showed very low reranking (high rank correlation between EBV) for fat and protein yields and modest reranking for milk yield.For milk yield the rank cor-relation between EBV_TM and EBV_HS was 0.85, whereas it was 0.96 and 0.98 for fat and protein yield, respectively.However, because producers tend to only use a small number of top-ranked sires for breeding purposes, this result has less implications.
Therefore, the rank correlations between the EBV for the 100 top-ranked bulls may be more of an accurate representation of the impact of G×E on production traits in Canada, which are shown in Table 8.This number of bulls was chosen because CDN publishes the ranking of 100 bulls in their top bull list.For milk yield, an average rank correlation of 0.53 across all regions was observed between EBV_TM and EBV_HS, whereas for fat yield and protein the average rank correlations were 0.80 and 0.86, respectively.In conclusion, there was strong evidence for G×E effect on milk yield under heat stress, causing substantial reranking of the top-ranked sires in Canada.There was evidence of smaller G×E effect on fat and protein yields resulting in less reranking of top-ranked sires.Bernabucci et al. (2014) also found that the inclusion of a heat stress function in their model caused the reranking of sires in Italian Holsteins, whereas Bohmanova et al. (2007) reported no changes in bull rankings between cooler and hotter regions in the United States.Therefore, there is conflicting evidence for G×E effect on production traits within countries.
The relationships between the EBV for milk, fat, and protein yield and THI for the 10 top-ranked bulls from the repeatability model fitting heat stress effect are depicted in Figure 2. Generally, the bulls did not have the same response to a high THI because the slope of the EBV differed from each other after the THI thresholds.Bulls with a slope closer to zero would have daughters less sensitive to heat stress, whereas daughters of bulls with a steeper slope would be more sensitive to heat stress.Overall, this variability in heat stress response shows that some sires that have superior performance in a thermoneutral environment may not perform as well at a high THI, whereas other top-ranked sires can maintain adequate performance despite experiencing heat stress.Therefore, bulls that are used for breeding purposes do not have the same daughters' heat stress response in terms of production traits.
Overall, this study also had some limitations.For instance, the heat stress function does not account for other specific factors that affect the incidence and severity of heat stress, such as heat stress management, barn stocking density, nutrition, or previous exposure to heat stress.As a result, the estimated genetic parameters for heat tolerance could have been affected.For example, the genetic parameters may have been underestimated because heat stress would be less evident in herds with extensive management and cooling devices (Freitas et al., 2006).
Furthermore, this model does not account for the effects of prolonged periods of heat stress.Heat waves can have a debilitating effect on dairy cattle due to the lack of nighttime cooling and high heat intensity (Igono et al., 1992).The THI equation used in these analyses may not optimally quantify heat stress in every region as well.For instance, it may underestimate heat stress in humid environments, such as Atlantic Maritime and Ontario.This is due to the relatively small weighting on RH compared with the weighting on AT (Bohmanova et al., 2007).Furthermore, this THI equation was developed to estimate human comfort during heat stress and the correlations with physiological parameters in dairy cattle have not been extensively researched.Additionally, test-day records only capture a fraction of the effect of heat stress on productivity and may not reflect a dairy cow's long-term response to heat stress.These record types cannot account for the cumulative effect of heat stress between test-days.Furthermore, there needs to be enough test-day records under heat stress to achieve high reliability of sire evaluations (Bohmanova et al., 2005).
Last, there is evidence that dairy cattle can have different thresholds for the onset of heat stress.Therefore, it may be problematic to assume a universal threshold.Sánchez et al. (2009) found that individual variation occurs in the onset of heat stress.However, they also reported that the production decline and individual thresholds were highly positively correlated.Therefore, breeding for a smaller slope after the threshold would simultaneously increase the threshold for the onset of The number of bulls that were used for each comparison are also listed.heat stress.Still, they concluded that a model with both selection criteria is still preferable.This is problematic because the estimation of individual thresholds is more computationally demanding (Sánchez et al., 2009).Furthermore, the model used in this study assumed that production traits decline linearly after the THI threshold, which may not reflect the actual pattern of production decay.A different modeling approach, such as a random regression model with a higher order polynomial, would have allowed for greater flexibility, because it does not require defining a static threshold level, and changes in genetic parameters could have been studied across the entire trajectory of a different covariate (Brügemann et al., 2011).However, these models are harder to implement due to their greater model complexity and poorer biological interpretation of resultant parameters (Ravagnolo and Misztal, 2000;Carabaño et al., 2014).
Several studies have investigated the relationship between heat tolerance and other traits, such as reproduction and risk of mortality in dairy cattle (Ravagnolo and Misztal, 2002;Vitali et al., 2009;Biffani et al., 2016).There is a negative genetic relationship between heat tolerance and reproduction (Biffani et al., 2016;Ravagnolo and Misztal, 2002).Vitali et al. (2009) also found that the summer had the highest frequency of death in dairy cattle in Italy.Therefore, the complete impact of heat stress on the profitability and sustainability of dairy farms is not fully captured by solely analyzing production traits such as in this study.Furthermore, maintaining homeostasis is a dynamic process that requires changes in some functions to stabilize others (Rauw and Gomez-Raya, 2015).It is important to consider the effect of reducing plasticity in production traits on other physiological processes.For instance, sires with improved heat tolerance transmit higher fertility, productive life, and stature and strength, but also lower production potential (Bohmanova et al., 2005;Aguilar et al., 2010).It may also be important to consider environmental flexibility.Environmental flexibility refers to the ability of animal to respond adaptively to environmental challenges by modifying their physiology and behavior.Phenotypic plasticity on the whole organism level is essential for environmental flexibility and adaptation.Low flexibility results in low survival, impaired fertility, high levels of culling, and loss of economic resources (Misztal, 2016).
Although maintaining production levels are a priority in the dairy industry, future studies should consider if a more holistic approach to studying heat tolerance in dairy cattle is possible to avoid altering the resilience of other traits.Future research could also focus on alternatives to THI as an environmental heat stress indicator, estimate the genetic parameters for heat tolerance in multiple parities, implementing a more sophisticated prediction models, such as a random Legendre polynomial regression model, or accounting for prolonged heat events.
This study did not estimate the genetic parameters for heat tolerance using both weather station and NASA POWER data to compare the 2 approaches, as THI values were shown to be highly correlated between the 2 weather data sources (Rockett et al., 2022).However, this could be done in a subset of herds that have both sources of weather data in a future study.

CONCLUSIONS
A novel genetic evaluation of heat tolerance in dairy cattle across Canada was completed with NASA POWER weather data instead of weather station data.The results of this study showed that heat tolerance has a low to moderate heritability and that variation in thermotolerance exists in the Canadian Holstein population.Therefore, genetic selection for heat tolerance is possible for Canadian dairy cattle, but the genetic gain may be slow due to the low heritability.This study also showed that breeding programs for Canadian dairy cattle should consider incorporating heat tolerance into future breeding goals because there is a moderate negative correlation between production traits and the ability of dairy cows to cope with heat stress.Overall, this study provided insight on the possibility of improving thermotolerance in Canadian Holsteins at national level using NASA POWER weather data, avoiding the spatial and temporal data gaps present in most weather station data sets.

Figure 1 .
Figure 1.Relationship between the temperature-humidity index (THI) and the average adjusted milk yield (g/d), average adjusted fat yield (g/d), and average adjusted protein yield (g/d) in Canada, as well as the estimated THI thresholds (dashed lines).
2 BC = British Columbia; AM = Atlantic Maritime.3 Number of cows used in these analyses.

2Figure 2 .
Figure 2. Relationship between the temperature-humidity index (THI) and the EBV for milk (first column), fat (second column), and protein yield (third column), for the top 10 bulls from parity 1. BC = British Columbia; AM = Atlantic Maritime.
Rockett et al.: GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINS

Table 1 .
Rockett et al.: GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINSThe number of records, cows, and herds used for the analysis of milk yield (g/d), fat yield (g/d), and protein yield (g/d), as well as the means and SD for each region Rockett et al.: GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINSwhere Rel i , PEV i , and F i are the reliability, prediction error variance, and inbreeding coefficient of the ith animal, respectively, and σ a 2 is the additive genetic variance of the trait.Last, the genetic correlations between EBV_HS and general additive genetic effects, and between EBV_HS and additive effects on the slopes were calculated for EBV_HS evaluated at 15 units above the THI threshold (EBV_HS 15 = a + 15 × α) for each trait:

Table 2 .
Rockett et al.: GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINSThe estimated rates of change (α) and SE for the average adjusted milk (g/d), fat (g/d), and protein (g/d) yields above the respective temperature-humidity index (THI) thresholds 1

Table 3 .
The average temperature-humidity index (THI) for spring, summer, and fall from 2010 to 2019 for each region as well as the number of herd locations that were used in the calculations 1 BC = British Columbia; AM = Atlantic Maritime.*The averages with an asterisk are above at least one of the selected THI thresholds for milk, fat, or protein yield.

Table 4 .
(Co)variance component estimates for the general additive genetic and permanent environmental effects, and for the slope of the additive genetic and permanent environmental effects evaluated at 15 units above the respective temperature-humidity index threshold for milk, fat, and protein yields (g/d) 1

Table 6 .
Rockett et al.:GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINS Heritability estimates for milk, fat, and protein yield with SE when temperature-humidity index (THI) was not included in the model (traditional model = TM), at the THI threshold for the model that included THI (heat stress model, under comfort zone = CZ), and at 15 units above THI threshold (heat stress model, under heat stress = HS)

Table 7 .
Rockett et al.: GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINSThe genetic correlations between the general additive genetic effects and the additive genetic effects on the slopes of change in trait phenotypes (r a, α ), heat stress EBV [EBV_HS; EBV_HS = a + f(m)α] 1 and general additive genetic effect (r EBV_HS, a ), and between heat stress EBV_HS and additive genetic effect on the slopes (r EBV_HS , α ) evaluated at f(m) = 15 units above the respective temperature-humidity index threshold for milk, fat, and protein yield ± SE ) is the number of units above the threshold and α is the slope of the additive genetic effect above the threshold.

Table 8 .
Rockett et al.:GENETIC EVALUATION OF HEAT TOLERANCE IN HOLSTEINS Rank correlations between the EBV for all the bulls with a reliability of 70% or greater and for the top 100 bulls with a reliability of 70% or greater from the traditional repeatability model (TM) and the EBV for the same bulls within the thermal comfort zone (CZ) as well as under heat stress (HS) at 15 units above the temperature-humidity index threshold for milk, fat, and protein yield 1 Region 2