Genetic evaluation of heat tolerance in Holstein cows in Japan

Abstract We used test‐day records and daily records from provincial weather stations in Japan to evaluate heat tolerance (HT) in Holstein cows according to a random regression test‐day model. Data were a total of 1,641,952 test‐day records for heritability estimates and 17,245,694 test‐day records for genetic evaluation of HT by using milk yield and somatic cell score (SCS) in Holstein cows that had calved for the first time in 2000 through 2015. Temperature–humidity index (THI) values were estimated by using average daily temperature and average daily relative humidity records from 60 provincial Japanese weather stations. The model contained herd–test‐day, with lactation curves on days in milk within month–age group as a fixed effect. General additive genetic effect and HT of additive genetic effect were included as random effects. The threshold value of THI was set to 60. For milk yield, estimated mean heritabilities were lower during heat stress (THI = 78; 0.20 and 0.28) than when below the heat stress threshold (THI ≤ 60; 0.26 and 0.31). For SCS, heritability estimates (range 0.08–0.10) were similar under all heat stress conditions. Genetic trends of HT indicated that EBVs of HT are changing in an undesirable direction.

Weather-wise, Japan is quite humid, especially in summer, and Holstein cows in southern areas experience severe HS during this season (Hagiya et al., 2017). Several studies have reported various effects of HS in Holstein cows in Japan (Atagi et al., 2018Hagiya et al., 2019;Nagamine & Sasaki, 2008). According to a simulation study of global warming that estimated temperature according to various RCP (Representative Concentration Pathway) scenarios for CO 2 conditions, average annual temperatures will increase by 0.3°C to 4.8°C in the 21st century (Ministry of the Environment, 2018). Therefore, in Japan, genetic improvement of heat tolerance (HT) is gaining interest in the dairy industry. Ravagnolo and Misztal (2000) proposed a model to estimate additive genetic variances caused simultaneously by general production and HT through random regression on THI and containing repeatability of animal effects on test-day milk yields. Bohmanova, Misztal, Tsuruta, Norman, and Lawlor (2008) applied their genetic evaluation methodology for HT to 56 million test-day records from first-parity Holstein cows in the United States. Aguilar, Misztal, and Tsuruta (2009) extended the model to evaluate HT by using milk, fat, and protein yields in a random regression test-day model on days in milk (DIM). Bernabucci et al. (2014) suggested that HT selection should be included in selection objectives, owing to existing unfavorable relationships between HT and milk production. Carabaño, Bachagha, Romón, and Díaz (2014) evaluated HT by using complex models containing higher than one-order Legendre polynomial regression on temperature hot and dry conditions in Spain. Nguyen et al. (2017) reported implementation of a breeding value for HT of dairy cattle in Australia and found that genetic trends declined over time. Atagi et al. (2018) reported genetic parameters of HT by using a random regression model on DIM and THI for milk, fat, protein, and somatic cell score (SCS) in southern Japan.
The aim of our current study was to estimate the genetic parameters of cows for HT, as indicated by changes in milk yield and SCS according to THI, by using a random regression model that applied information from daily weather records collected at provincial weather stations throughout Japan. In addition, we compared the estimated breeding values (EBVs) of HT between estimates based on milk and on SCS.

| Data
Data contained test-day records of milk and SCC in Holstein cows that had calved for the first time in 2000 through 2015. Records were collected through the Dairy Herd Improvement Program and were provided by the Livestock Improvement Association of Japan (Tokyo, Japan). These data included records of test-day milk yield and SCC for first-lactation cows from all over Japan. SCCs were log-transformed into SCSs by using the formula SCS = log 2 (SCC/100,000) + 3 (Ali & Shook, 1980). Weather records from 60 provincial weather stations were obtained from the website of MeteoCrop DB (National Institute for Agro-Environmental Science, NARO, 2017). THI values were estimated by using average daily temperature and average daily relative humidity, according to Hagiya et al. (2019). First, THI was estimated by using the following formula (NRC, 1971): where t is temperature in degrees Celsius and rh is relative humidity as a percentage.
Test-day records were linked to data from each of the nearest weather stations of the 14 on Hokkaido, Japan's northernmost island, and from stations in the other 46 prefectures. Hagiya et al. (2019) reported that the lag in response to a heat effect was about 3 days for milk yield and 8 to 10 days for SCC under Japanese weather conditions. Therefore, we set the length of the lag in response on test day to represent heat effects as 3 and 8 days for milk yield and SCS, respectively. To estimate genetic parameters, 10 subsets of data were chosen randomly according to herd; each subset included more than 140,000 records and 40,000 animals (Table 1). Data for genetic evaluation contained 17,245,694 test-day records from 2,018,402 cows. The pedigree file was provided by the Holstein Cattle Association of Japan (Tokyo, Japan) and comprised data on 2,987,536 animals, including cows with test-day records and whose lineage could be traced back five generations.

| Model
Genetic and environmental (co)variances for milk yield and SCS were estimated by using a random regression animal model. The general additive genetic and general permanent environmental (PE) effects were modeled by using the Legendre polynomial

Implications
Heat stress in Holstein cattle is a growing concern in Japan because of its negative effects on dairy cows and associated economic losses in the dairy industry. To improve cow' heat tolerance, we investigated the heritabilities of heat tolerance based on milk yield and cow health. Correlations between genetic performances of heat tolerance estimated by using milk yield and those based on a cow health indicator were close to zero. Therefore, to improve heat tolerance in dairy cattle, we recommend selecting animals that are concurrently superior according to heat tolerance based on milk yield and as well as cow health. function, with third-order functions for milk yield and secondorder functions for SCS. We chose fixed effects by referring to the national genetic evaluation for production traits in Japan (Hagiya, 2019), as follows:

TA B L E 1 Number of records, mean, and standard deviation (SD) in each data subset
where y ijklt = observation of test-day milk or SCS; HTD i = fixed effect of herd × test-day i; AM jl = fixed regression coefficients l for lactation curves of age × month j (15 age groups × 12 calendar months = 180 subclasses); u kl = the general random additive genetic effect of animal k corresponding polynomial coefficients (covariates) l for lactation curve; uh k = the random additive genetic effect of HT of animal k corresponding linear regression coefficient (covariate) for HS in herd i; pe kl = the general random PE effect of animal k corresponding to polynomial coefficients (covariates) l for lactation curve; peh k = the random PE effect of HT of animal k corresponding linear regression coefficient (covariate) for HS in herd i; e ijklt = vectors of random residual effects; 05t , or the fourth-order Legendre polynomial and Wilmink's function (Wilmink, 1987)  where G = (m + 1) × (m + 1) additive genetic (co)variance matrix for u; g = additive genetic covariance vector between u and uh; 2 uh = additive genetic variance for uh; ⊗ = the Kronecker product; A = the additive genetic relationship for animals; P = (m + 2) × (m + 2) PE (co)variance matrix for pe and peh; I = the identity matrix for cows; and R = the residual variance.
We tested two models for random regression on DIM. Model 1 contained m = 0 (i.e., additive genetic variance and PE variance were constant throughout the lactation period), and Model 2 contained m = 3 for milk yield and m = 2 for SCS (i.e., random regression was adopted for variances changing by DIM). In the preliminary study, we tried third-order Legendre polynomial functions for genetic and PE effects on SCS, but most subsets did not yield appropriate convergence.
Genetic correlations between HT on THI i and milk yield or SCS at DIM t were obtained by Genetic parameters for the random regression models were estimated by using the AIREMLF90 program (Misztal et al., 2002), and EBVs from Model 1 and from Model 2 were estimated by using the MMEF program (T. Osawa, personal communication, 19 March 2020) with the preconditioned conjugate gradients method, which is used to estimate the national genetic evaluation in Japan. We compared the correlations between EBVs of THI from Model 1 with those from Model 2.

| Summary statistics
Means of milk yield and SCS are given in Table 1. Daily milk yield ranged from 26.5 to 27.1 kg and daily SCS from 2.31 to 2.41 (Table 1).
THI ranged from 4 to 84 (Figure 1). The average THI was 50.6 and the mode was 62.

| Changes in heritabilities on THI
Genetic parameter estimates of test-day milk yield were obtained from all subsets for Model 1 and from 8 subsets in Model 2 (  When the values of THI exceeded a threshold (THI = 60), heritability estimates and genetic variances from Model 1 decreased with increasing THI for test-day milk yield but did not differ for SCS ( Figure 4). For both traits, PE variances increased with increasing THI. Heritability estimates for milk yield were nearly constant at 10 DIM ( Figure 5) but tended to decrease with increasing THI (DIM = 200). For SCS, heritability estimates and genetic variances did not vary with increasing DIM (Figure 6).

| Genetic correlations
Estimated genetic correlations between general and HT additive genetic effects for milk yield were -0.44 according to Model 1 and ranged from -0.13 to -0.46 with Model 2 (Figure 7). In Model 2, trends for genetic correlations between general and HT additive genetic effects on DIM for milk yield decreased with increasing DIM. In early lactation, estimated genetic correlations for milk yield tended to be weaker for Model 2 than Model 1, but estimates from both models were similar during mid-lactation. Estimated genetic correlations for SCS were close to zero throughout lactation in both models.

| Correlations between EBVs for heat durability
Estimated EBV showed a high correlation (0.77) between models within each milk yield and SCS class (Table 4). Correlation between models for EBVs for milk yield and SCS were negative but close to zero (range, -0.01 to -0.18).

| Genetic trends of heat durability
Genetic trends of HT estimated by using milk yield decreased over time in both models ( Figure 8). Those estimated by using SCS increased with increasing year. According to estimates from both

| D ISCUSS I ON
Daily milk yield (range, 26.9 to 27.2 kg) and SCS (range, 2.32 to 2.34) in the current study were similar to those recently reported (daily milk yield, 26.9 kg; SCS, 2.3 to 2.5) for Holstein cows in Japan (Hagiya et al., 2017Yamazaki et al., 2016). The shapes of the histograms (Figure 1) were almost the same as those in previous work (Hagiya et al., 2019). Approximately 30% of total days were THI > 60: that is, Holstein cows in Japan are under HS for 120 days each year.
We were unable to estimate genetic parameters for some subsets (Tables 2 and 3). When initial values for a subset were appropriate, estimation involved fewer rounds of calculation than when initial values were suboptimal. That is, when initial values were poor, estimation required numerous rounds of calculation and excessive computational time. Moreover, in some cases, we were unable to find any adequate initial values at all and therefore could not obtain appropriate estimates.
For both models, estimated heritability for milk yield was higher in the absence of HS (THI ≤ 60) than under HS conditions Our results agree with theirs for milk yield and are slightly higher for SCS, but the trends for estimates without and with HS were similar ( Figure 4). In addition, our heritability estimates were constant with increasing THI during early lactation and decreased with increasing THI later during lactation ( Figure 5). Although these changes were not statistically significant, they suggest that HS might have an increased effect during late lactation. Heritability estimates for SCS did not differ throughout lactation. Estimated adversely affects not only milk yield and SCS but also milk fat, milk protein, female fertility, and other traits (e.g., Atagi et al., 2018;Hagiya et al., 2017). Therefore, genetic improvement of HT is very important, especially in southern Japan.
The correlation between EBV of HT was 0.77 for each trait and in both models, thus yielding similar EBVs (Table 4). However, correlations between EBVs of HT based on milk yield and on SCS were both close to zero. This finding suggests that HT EBVs estimated by using milk yield and those estimated by using SCS reflected different traits. EBV of HT with milk yield indicates HT effects on production performance, and that with SCS indicates HT effects on cow health. Nguyen et al. (2017) showed an example of EBV HT that used weighted EBVs of HT based on milk, fat, and protein. For calculating the total EBV of HT, we recommend weighting EBVs of HT based on SCS as a health indicator and EBVs of HT as production indicators based on yield traits.

| CON CLUS ION
We investigated genetic parameters in the responses of milk yield and SCS to HS in cows by using random regression models that incorporated daily weather records from provincial weather stations throughout Japan. In addition, we compared the EBVs of HT between estimates based on milk with those for SCS. Mean estimated heritabilities were lower during HS than when below the HS threshold for milk yield but were similar regardless of THI for SCS.
Moderate negative genetic correlations between general additive genetic effects and those for HT were estimated for milk yield, but SCS showed no such genetic correlations. Correlations between EBVs of HT based on milk yield and those based on SCS were negative but close to zero.