INTRODUCTION
Heat stress (HS), caused by high temperature and humidity, is a harmful issue in dairy farms. Under HS conditions, cows show decreased feed intake and milk yield, but increased somatic cell counts (SCC) [1-4], thereby deteriorating milk quality. Hammami et al. [5] reported a reduction in milk composition (5.7 kg of milk fat/year, and 4.2 kg of protein/year) in a European dairy farm that had undergone severe HS. Similar to global climate change, an increase in air temperature has been observed in South Korea. It has been reported that the increase in annual average temperature has been 0.5°C since the decade between 2010 and 2019 [6]; therefore, the decrease in dairy productivity is an emerging issue in South Korea. The temperature-humidity index (THI) is a value estimated using temperature (dry or wet) and humidity (direct or relative) to measure the degree of HS. It was first applied to monitor human health [7]; however, recently, livestock nutritionists have also been widely used THI to monitor the relationship between critical THI and the emergence of HS in animals. For example, in dairy cattle, a decrease in milk yield and feed intake [2,8] was observed when the THI for dairy cattle was over 72. Recently, high-producing cows were found to be more sensitive to HS; therefore, a decrease in milk yield was observed at a THI of 68 in the US [9]. Bohmanova et al. [10] reported that the critical THI at which HS was observed in dairy cattle varied according to regional characteristics.
Recently, big data has been used to elucidate the association between HS and dairy production, along with technological advances in data accumulation, computer hardware, and sensors to monitor environmental parameters. For example, Hagiya et al. [11] developed 2-phase linear model using 17,245,709 test-day records from 2,018,406 cows and demonstrated that Japanese Holstein cows exhibited drastic milk depression at a THI of 70.4. Similar to a previous study [11], both dairy milking and climate data were provided from public data centers in South Korea; however, to the best of our knowledge, little research has been conducted to demonstrate the association between high THI and milk performance in Korean dairy cattle systems. Therefore, the aim of this study was to investigate how HS adversely affects milk production and traits in South Korea using big data (daily milking, air temperature, humidity, and region of farms). The threshold values of THI for milk performance were also estimated using segmented linear regression analysis.
MATERIALS AND METHODS
Test-day records of milk yield, fat- and protein-corrected milk (FPCM), milk composition (fat, protein, and urea nitrogen [MUN]), and SCC were collected from the Dairy Cattle Improvement Center (NongHyup Agribusiness Group, Korea) in Korea. Among the milk composition data, SCC was log-transformed into somatic cell scores (SCS) using the following equation [12]:
and FPCM was estimated as suggested by [13]:
Data were collected from 215,276 Holstein cows (122,087 primiparous cows and 93,189 multiparous [2–4 parities]) at 1–305 days in milk (DIM) for three years (2017–2020). The total number of herds included in this study was 2,419. In total, 1,498,232 records were included in this study. The descriptive statistics used in this study are presented in Table 1. Weather records (July 2017 to April 2020) from 600 weather stations near the respective herds were acquired from the Korea Meteorological Administration (KMA) website. To collect weather records for each herd, the weather records of the automatic weather checking station located closest to the respective herd were used, and the daily THI was calculated using the following equation [14]:
where “T” is the daily average temperature (°C), and “RH” is the relative humidity (%).
Test-day records for milk traits were linked to the daily average THI calculated from weather records. The effects of HS on milk traits were estimated using the following statistical model:
where Yijklmn is an observation of test-day records for milk traits or SCS; Ri is the fixed effect of the region (five subclasses); Yj is the fixed effect of year at calving (4 subclasses); Mk is the fixed effect of month (12 months); Pl is the fixed effect of parities (primiparous and multiparous); DIMm denotes the days in milk m (305 subclasses); THIn is the index HS as expressed by THI (40 subclasses); and eijklmn represents the vector of random residual effects. The distribution of data for the fixed effects used in the model is shown in Fig. 1.
Considering that HS had significant results linearly, the break point (BP) of THI was evaluated using segmented regression analysis on R [15] Segmented package [16]. The least-squares mean (LSM) of milk traits in different THI, which was calculated from Equation (1), was used as the dependent variable. The slope of the segmented linear regression was assumed as follows:
where y*i is the LSM of milk traits in different THI; a is an intercept; b1 is a regression coefficient on THI Xi when THI Xi is lower than the BP; b2 is a regression coefficient on THI Xi when THI Xi exceeds the BP; ei is the random residual term, and BP is the break point, defined as the appropriate threshold value of THI. Linear regression was applied to the heat stress effect.
RESULTS AND DISCUSSION
The average milk yield and milk composition are presented in Table 1. The average daily milk yield, milk fat, and protein concentrations were 33.9 kg/d, 3.9%, and 3.2%, respectively. These results were in agreement with the annual domestic dairy production statistics for 2021 (milk yield, 34.1 kg/d; fat, 3.9%; and protein, 3.2%) [17]. Milk traits (milk yield, FPCM, milk fat, milk protein, MUN, and SCC) were expressed as LSM values adjusted for region, age, month, parity, DIM, and THI (Fig. 2). The BP of THI for milk yield was 71.9 and 70.6 for primiparous and multiparous cows, respectively (Fig. 2A). Similar to milk yield, FPCM, milk fat, and milk protein yields also had specific BP (Figs. 2B, 2C, and 2D). For all milk yield parameters, the BP in primiparous cows was higher than that in multiparous cows. After exceeding the BP, milk yield gradually decreased with increasing THI, regardless of parity. A decrease in each milk component was also observed when the THI exceeded the BP. The results of BP for milk and component yield suggested that dairy cows in South Korea started to exhibit thermal HS at around THI 70 (Figs. 2A, 2B, 2C, and 2D). Different THI thresholds have been suggested in previous studies. Traditionally, it is believed that the threshold value between thermal comfort and mild HS is 72 [18], while De Rensis et al. [19] reported that cows exhibited signs of mild HS when the experimental environment was designed to have a THI over 68. Hammami et al. [5] reported that a reduction in milk yield was observed after THI of 62. In Japan, the THI-BP for milk yield was 70.4 [11], similar to the results of this study. The difference in critical THI for dairy production might be explained by several environmental factors (cooling system, climatic region, and degree of genetic improvement) influencing the relationship between THI and production [20-23]. Higher THI thresholds for milk performance traits were found in primiparous cows than in multiparous cows (71.9 and 70.6 in milk yield; 71.2 and 70.5 in milk FPCM yield; 71.0 and 69.5 in milk fat yield; 71.5 and 69.9 in milk protein yield, respectively), indicating that multiparous cows were more susceptible to HS. This finding might be explained by the fact that primiparous cows generate far less metabolic heat because they have a greater surface area compared with internal body mass and lower milk production. Because of metabolic heat production during milk synthesis, a severe response to HS is more easily observed in high-yielding cows than in normal cows [24].
The MUN concentration increased drastically after BP of THI (74.6 in primiparas, and 74.3 in multiparous; Fig. 2E). Under HS, cows experience low ruminal pH and inefficient microbial crude protein production, thereby increasing ruminal ammonia concentrations [25]. Cowley et al. [26] suggested that the reduction in milk protein from heat-stressed cows is the result of downregulation of mammary protein synthetic activity, but also increased MUN due to catabolized muscle tissue, a process that is intensified under HS conditions. Similar to a previous study, a reduction in milk protein yield was observed at approximately a 70 of THI (Fig. 2D). In periods of HS, Koch et al. [27] observed stable transcription rates of enzymes influencing the urea cycle, and the group explained the increased MUN concentrations with lower kidney perfusion in periods of HS. Moreover, the observed high MUN levels indicated an energy deficiency due to reduced feed intake during HS.
The BP of THI on SCS was 66.7 in primiparous cows and 53.6 in multiparous cows (Fig. 2F). When the threshold was exceeded, the SCS sharply increased in primiparous cows; however, a similar trend was not observed in multiparous cows. In Japan, a dramatic increase in SCS was also observed when the THI exceeded 68.5 [11], and the THI calculation method was the same as that used in this study [14]. Pragna et al. [24] explained that HS might impair the dairy immune system, which is eventually linked to udder infection. Based on the data obtained in this study, the increase in udder infection was not known; however, it is possible to speculate that HS over BP can influence udder health and increase SCS, especially in primiparous cows. In conclusion, this study represents the first comprehensive research on the effect of THI on Holstein cattle raised in South Korea. The negative effect of HS on dairy cattle performance, which has been reported in previous studies abroad, was also confirmed in this study by analyzing big data on dairy performance and regional climates. Therefore, more detailed feeding management and construction of a farm cooling system are required to prevent HS from dairy cows when the THI is over the specific BP.