Estimation of the genetic milk yield parameters of Holstein cattle under heat stress in South Korea

Objective The objective of this study was to investigate the genetic components of daily milk yield and to re-rank bulls in South Korea by estimated breeding value (EBV) under heat stress using the temperature-humidity index (THI). Methods This study was conducted using 125,312 monthly test-day records, collected from January 2000 to February 2017 for 19,889 Holstein cows from 647 farms in South Korea. Milk production data were collected from two agencies, the Dairy Cattle Genetic Improvement Center and the Korea Animal Improvement Association, and meteorological data were obtained from 41 regional weather stations using the Automated Surface Observing System (ASOS) installed throughout South Korea. A random regression model using the THI was applied to estimate genetic parameters of heat tolerance based on the test-day records. The model included herd-year-season, calving age, and days-in-milk as fixed effects, as well as heat tolerance as an additive genetic effect, permanent environmental effect, and direct additive and permanent environmental effect. Results Below the THI threshold (≤72; no heat stress), the variance in heat tolerance was zero. However, the heat tolerance variance began to increase as THI exceeded the threshold. The covariance between the genetic additive effect and the heat tolerance effect was −0.33. Heritability estimates of milk yield ranged from 0.111 to 0.176 (average: 0.128). Heritability decreased slightly as THI increased, and began to increase at a THI of 79. The predicted bull EBV ranking varied with THI. Conclusion We conclude that genetic evaluation using the THI function could be useful for selecting bulls for heat tolerance in South Korea.


INTRODUCTION
In South Korea, it has become necessary to evaluate dairy bull performance based on quantifiable factors exerting a beyond the Korean Peninsula, because the country has started to export dairy bull semen to countries with distinctly different environments. One such factor is climatic changes, which affect bull performance. Previous studies have suggested that the Earth's climate has warmed during the past century, and this change has continued at an unprecedented and continuous rate in recent years [1]. These recent climate changes are challenging for the dairy industry because animals in high-temperature environments voluntarily decrease their feed intake to reduce digestive heat production, which in turn causes a reduction in milk yield and a negative energy balance [2][3][4]. Cows with high milking abilities, in particular, are more susceptible to such heat-stress conditions [5]. These conditions drive dairy breeding schemes toward improvements in genetic predisposition for higher milking ability under stressful conditions. Nearly two decades ago, a study by Misztal [6] introduced a temperature-humidity index (THI)-based methodology, combining test-day records with weather data collected from public weather stations, for the genetic evaluation of animals with respect to heat tolerance. A clear advantage of using public weather station data is their reliability, because these stations do not rely on individual measurements of body temperature or respiration, which are difficult to generate at large scales. Ravagnolo and Misztal [7] also suggested the use of public weather data for national evaluations. The THI is a bioclimatic index now commonly used to determine the degree of heat stress in studies of animals including dairy cattle [5,8,9]. The THI is a unit measurement that represents a combination of air temperature and humidity [5,9,10]. Generally, the THI function considers three scenarios. First, the daily production (milk yield) of a cow begins to decrease after a THI threshold. Such production drops are typically nonlinear, but can be linearized by the THI function if necessary. Second, there should be genetic variability in susceptibility to heat stress. Dairy production by individual cows may drop at one particular threshold (scenario 1), at various different thresholds (scenario 2), or under both scenarios (scenario 3). There may also be differential individual rates of production decline as indexed by the THI. However, variability in the differential rate of decline (scenario 1) is preferred when using the THI function due to its better theoretical tractability; it is theoretically easier to express than other scenarios. Third, heat stress should comprise both genetic and non-genetic components, where two components are changed by the same curve (THI function) [7]. Therefore, the objective of this study was to investigate the effect of THI on the genetic variance components of daily milk yield and probable re-ranking of bulls in South Korea by estimated breeding value (EBV) under heat stress conditions.

Data
This study was conducted using 125,312 monthly test-day records, collected from January 2000 to February 2017 for 19,889 Holstein cows from 647 farms in South Korea. Milk production data were provided by two agencies: the Dairy Cattle Genetic Improvement Center and the Korea Animal Improvement Association. We excluded records with extreme milk production (<2 or >61 kg), age at calving (<17 or >31 months), days in milk (DIM, <5 or >500 days) and THI (<61 on test-day) from the data set. Milk records for each cow with at least three test-day records on test days when the THI was >72 in each lactation were retained for analysis. The minimum size of each contemporary group (herd-year-season, HYS) was five. Each cow was milked twice per day.
Meteorological data were obtained from 41 regional weather stations with the Automated Surface Observing System (ASOS) installed throughout Korea. The data included the maximum, minimum, and average daily temperature, and average daily humidity records. Weather variables used in this study included daily maximum temperature and daily average humidity. Test-day THI records indicate the heat stress response, according to Ravagnolo et al [7].
The THI was calculated for each day and weather station using the following expression [10]: Where T and RH are the maximum daily temperature (°C) and average relative humidity (%), respectively. Table 1 shows the distribution of data by calving age, calving season, and THI.

Statistical analysis
Genetic components of heat stress and other variables were estimated using the following mixed-model equation: Where y ijkl is the milk yield of the lth animal of the ith herdyear-season, recorded at the jth calving age and the kth DIM; HYS i is the effect of the ith herd-year-season (i = 1 to 2,483); age j is the effect of the jth calving age (j = 1 to 2); DIM k is the effect of the kth DIM (k = 1 to 50); a l and p l are the effects of the additive genetic merit and permanent environment of the cow l (l = 1 to 56,132); f(i) is the heat stress function for the herd-test-day (HTD); v l and q l are the additive and permanent environmental effects of heat tolerance on cow l, respectively; and e ijkl is the residual effect. The variances are: This model can be used to illustrate the repeatability of the model with random regression of the heat stress function, which determines the effect of heat tolerance by calculating the relative change in production per unit of heat stress. To simplify the calculation, residual variance was assumed to be homogeneous. Our preliminary study showed that milk yield started to decrease at a threshold THI of 72. However, milk protein yield and milk fat yield did not meet the threshold because they continued to decline. Therefore, we considered the THI threshold to be 72 [11]. The heat-stress function was defined as [7]: The variance components were estimated using the EM-REML algorithm in the REMLF90 software module from the BLUPF90 suite of programs [12].
The heritability value used to determine production at heat stress level f(i) was calculated as: The genetic correlation between the genetic and heat stress effects was calculated as: , Figure 1 shows the least-square solutions for DIM classes.

RESULTS AND DISCUSSION
These solutions tend to increase before the peak day class (DIM >90) and then decrease. The shapes of the least-squares curves for the DIM classes were similar to a normal lactation curve. In genetic evaluations, test-day records may cause overor underestimation of parameter estimates due to the different scales of DIM records. Therefore, the test-day model must incorporate the general shape of the lactation curve [13]. Cows exhibit different DIM effects on the same test days. Test-day evaluation results were closely correlated (r = 0.87-0.97) to the 305-day evaluation in a previous study [13]. Although a sophisticated test-day model, i.e., a random regression model of the lactation curve, can capture heat stress variation throughout the lactation period, this method has much greater difficulty in estimating the parameters [7]. Therefore, the simple test-day model is more appropriate than the sophisticated test-day model. The variance component estimates of milk yield are presented in Table 2. The THI values indicate that the variance in genetic heat tolerance was small. The covariance between the genetic additive and heat tolerance effects (r a,v ) was -0.33, indicating that the selecting animals for high milk yield could produce animals susceptible to heat stress. According to a previous study, the total body heat load of lactating cows increased according to the metabolic heat production associated with milk production [14], in turn increasing the ability to maintain homeothermy under heat stress conditions [15]. Other studies have reported correlation coefficients of -0.33 [7] and -0.38 [5], and other similar values [16,17]. Figure 2 shows estimates of additive genetic variance and heat tolerance variance with increasing THI. Under the THI threshold (≤72, no heat stress), the variance in heat tolerance was zero. However, the heat tolerance variance began to increase as THI exceeded the threshold. The total variance was smaller than the genetic variance at THI <86 due to the influence of negative covariance. These findings are similar to those of previous studies [5,7,17,18]. Figure 3 shows estimates of heritability for milk yield at different THI values. These estimates ranged from 0.111 to 0.176 (average: 0.128). Heritability decreased slightly as THI increased, and then began to increase at a THI of 79. However, some of these changes may not have been biological, instead being consequences of the random regression model of the THI function [7]. Our estimates were lower than those reported in previous studies that applied the same model (Heritability ranges have been reported as 0.10 to 0.25 by Aguilar et al [19] and 0.16 to 0.21 by Ravagnolo and Misztal [16]) [5,7,17], suggesting that the characteristics of the South Korean Holstein     Table 3 shows the inbreeding coefficient for the data used in this study. Generally, heritability decreased slightly due to increased inbreeding. The South Korean dairy population is small, at nearly 400,000, with a high average inbreeding coefficient of approximately 2.03% in 2012 [19]. The average inbreeding coefficient in this study was 1.5%, and 2,355 of the 56,441 head in pedigree were sires, a smaller population size and higher average inbreeding coefficient than used in previous studies. In the test-day model, HTD is a more suitable variable than herd-year-season. However, the Holstein cattle population in South Korea is small; therefore, we cannot use HTD as a variable because the contemporary group size was insufficient due to the small dairy population size. These differences would lead to larger residual variance than reported in other studies. Hence, our heritability value was lower than that reported elsewhere.
Heat tolerance is among the most important adaptive traits in cattle [5,20,21]. If bulls have similar EBVs, we should consider the effect of heat tolerance on EBV when temperatures are high. Figure 4 shows the genetic correlation between predicted breeding values at a THI ≤72 (no heat stress) and breeding    values at a given THI for milk yields. Below the heat stress threshold, the correlation was 1 and started to decrease as THI increased, such that the bull EBV ranking differed according to THI. Table 4 shows the top 10 bulls, having at least 50 daughters, with respect to the EBV for milk yield at different THI values. One bull had negative heat tolerance, and its EBV for milk yield exhibited a greater decrease as THI increased. Bull s1 had the highest EBV for milk yield at THI = 0 (no heat stress). However, bull s1 ranked fourth, and bull s4 ranked first, at THI = 80, because these bulls had a negative and positive heat tolerance EBV, respectively. When heat stress is a factor influencing daily milk yields, selection for high productivity would be useful only if these animals are able to maintain high productivity under heat stress [5,22]. Bull s1 is a good example of this type of condition, because its ranking decreased slightly under heat stress. The heat stress function used meteorological data collected from public weather stations, rather than directly from farms; however, genetic evaluation should not be affected by this substitution, because the use of HYS compensates for farm effects.

CONCLUSION
We conclude that genetic evaluation using THI could be applied to select bulls in South Korea for thermotolerance. Because daily milk yield and heat tolerance are negatively correlated, selection for high milk production could decrease thermotolerance. Therefore, sire rankings should be changed to incorporate the effects of high temperature and humidity. Our results will be helpful when South Korean dairy cattle are adapting to environmental conditions in other countries, or are selected in a changing climate.

CONFLICT OF INTEREST
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.