Heat stress resilience in Holstein cows: milk yield response to THI

Animals, housing and management

Data was collected on a commercial farm located in Eastern Croatia. As the study was based on routinely recorded farm data, ethical approval was not required. From a farm of 1,440 lactating cows, data from the most productive group of 169 Holstein cows were used in this study. This specific sub-group was selected in coordination with the farm management based on its high metabolic demand and lactation stage, representing the most sensitive physiological window for evaluating heat stress resilience. Furthermore, the group was housed in a dedicated barn section equipped with high-precision biometric monitoring equipment. Cows were in their second to sixth lactation. Parity distribution consisted of second (n=93; 55.0 %), third (n = 46; 27.2 %), fourth (n=19; 11.2 %), fifth (n=9, 5.3 %), and sixth (n=2; 1.2 %) lactations. The average daily milk yield in the analysed group was 38.72 kg, while the average days in milk (DIM) were 124.3 days (median 123.0 days). The study period spanned from March 1 to August 19, 2025, comprising 28,285 daily milk yield records. Daily milk yield records were exported by Delpro herd management software (DeLaval, Tumba, Sweden). Total daily production was calculated as the sum of the three individual milkings. Records with missing values were excluded from the database prior to statistical analysis. Cows were housed in a free-stall barn with cubicles, allowing unrestricted movement and ad libitum access to the feed bunk, automatic waterers, and resting areas. The barn was equipped with a cooling system comprising fans and misters; the ventilation system was programmed to activate when the air temperature exceeded 12 °C, while the misters along the feed bunk operated according to a predefined schedule. Milking was performed three times daily in a 64-stall rotary milking parlour at 06:00, 11:30, and 19:00 h to ensure consistent milking intervals. Individual cow identification was performed automatically using neck collars.

Microclimate monitoring and THI calculation

Air temperature and relative humidity were continuously monitored using data loggers (ExTech RHT20, Extech Instruments, USA) positioned at three farm locations: the feed bunk in the barn (Location 1), the outdoor passage exposed to sunlight (Location 2), and the area within the rotary milking parlour (Location 3). Recorded values were saved every hour, resulting in 24 daily records for each location. For statistical analysis, microclimatic data from Location 1 were used as a representative indicator of the continuous environmental exposure in the facility. The use of barn-level measurements for heat stress analyses is recommended in the literature as a more realistic indicator of the actual conditions that animals experience compared to data from external weather stations (Brügemann et al., 2012; Ouellet et al., 2019). Two forms of the THI were defined within the analysis, THI and THI_MAX, representing standard bioclimatic indicators of heat stress in dairy cows (Bernabucci et al., 2014; Moore et al., 2023). These two indicators enabled a comparison between indices based on average (THI) and extreme daily temperature and humidity values (THI_MAX), encompassing a wide range of heat load conditions (Ravagnolo and Misztal, 2000; Brügemann et al., 2012). The THI_MAX model relies on the assumption that daily maximum temperature combined with minimum relative humidity more accurately reflects the actual heat load animals are exposed to during the day (Ravagnolo and Misztal, 2000; Brügemann et al., 2012; Ouellet et al., 2019). Daily average THI values were calculated according to the National Research Council (NRC, 1971) formula, using daily average temperature (T, °C) and daily average relative humidity (RH, %):

THI = (1.8 × T + 32) – (0.55 – 0.0055 × RH) × (1.8 × T – 26) /1/

THI_MAX was calculated using the same mathematical formula as THI. First, hourly THI values were calculated from hourly temperature and relative humidity records, after which the maximum THI value for each day was selected to represent peak daily heat-stress conditions.

Statistical analyses

Statistical analysis was performed in the R programming environment (R Core Team, 2025) using two models (RRM and LHR) on 28,285 filtered observations. RRM analysis was performed using the lme4 (Bates et al., 2015) package, while the LHR analysis was performed using ggeffects package (Lüdecke, 2018). Results were visualized using ggplot2 (Wickham, 2016) package. Statistical significance level was set at p<0.05. The RRM was used to assess individual variability in cow responses to heat stress using the following model:

y_ijk = μ + β₁THI_ij + β₂DIM_ij + Parity_k +c_i + b_iTHI_ij + ε_ijk /2/

where y_ijk represents the daily milk yield of cow i on day j within parity class k, μ is the overall mean, THI_ij is the temperature-humidity index recorded for cow i on day j, DIM_ij is days in milk for cow i on day j, Parity_k is the fixed effect of parity class (ranging from the second to the sixth lactation), β₁ and β₂ are regression coefficients for THI and DIM, respectively, c_i is the cow-specific random intercept, b_i is the cow-specific random regression coefficient describing the individual response to THI, and ε_ijk is the residual error. THI and THI_MAX models were analysed separately.

The second method, LHR, was used to quantify the time-lagged response in milk yield following identified periods of acute heat stress, using the following model:

milk_dev_{ijk =} μ + day_rel_j​ + c_i​ + e_k ​+ ε_ijk​ /3/

where milk_dev_ijk represents the deviation of daily milk yield per cow i on relative day j within heat-stress event k, calculated relative to the reference milk yield on day -1, μ represents the overall mean, day_rel_j is the fixed effect of relative day, describing the temporal response pattern from three days before the event to five days after exposure, ranging from -3 to +5, cᵢ is the random cow effect, accounting for repeated observations and baseline differences among cows, eₖ is the random heat-event effect, accounting for variability among individual heat-stress episodes, εᵢⱼₖ is the residual error.

As previously mentioned, the event window was defined from day -3 to day +5 to capture both pre-event baseline stability and delayed post exposure production responses. Days -3 and -2 were included to evaluate whether milk yield was already changing before the identified thermal shock. Day -1 was used as the reference baseline because it immediately preceded the acute increase in THI or THI_MAX. Days +1 through +5 were selected to capture delayed biological responses in milk yield following heat exposure.

An event was defined as a day on which the THI or THI_MAX value exceeded a specified threshold while increasing by at least two units compared with the previous day. Each new event was identified only after the index first fell below the threshold and subsequently exceeded it again, ensuring that repeated days within the same heat episode were not treated as independent events. The lagged heat-stress response analysis comprised six thresholds: three based on THI thresholds (≥ 66, ≥ 70 and ≥ 74) and three based on THI_MAX thresholds (≥72, ≥76 and ≥80).

Microclimatic analyses

During the analysed period of nearly six months, significant seasonal variability in microclimatic conditions was observed. Daily ambient temperatures averaged 18.70 °C (S =6.52), while the average maximum daily temperature reached 23.91 °C (SD=7.17 °C; max=41.80 °C). Bioclimatic indices showed that average THI values were recorded at 63.40 (SD = 9.29; range 37.98-78.37), whereas average THI_MAX was 69.70 (SD=8.81; range 45.92-86.71). Higher values for THI_MAX compared to average THI directly support the assumption that utilizing maximum daily temperature and minimum relative humidity results in possibly more suitable bioclimatic indicator values (Ravagnolo and Misztal, 2000; Ouellet et al., 2019). Seasonal dynamics clearly illustrated a progression in heat load, with average THI values increasing systematically from 51.06 in March to 71.31 in August, while extreme THI_MAX peaks reached 86.71. Based on heat stress classifications, THI values between 68 and 74 indicate mild stress, while levels exceeding 75 indicate high stress which leads to a sharp decline in production performance (Armstrong, 1994; Polsky and von Keyserlingk, 2017; Herbut et al., 2018). Early studies established a THI of 72 as the critical threshold for milk yield decline (Armstrong, 1994) whereas subsequent research has identified a lower onset of 68 in high-yielding Holstein cows (Zimbelman et al., 2009; Collier et al., 2012; Bernabucci et al., 2014). THI≥66 was exceeded on 87 days (48.6 % of the study period), THI≥70 on 57 days (31.8 %), and THI≥74 on 12 days (6.7 %). THI_MAX was above the threshold of ≥72 on 87 days (48.6 %), THI_MAX≥76 on 49 days (27.4 %), and THI_MAX≥80 on 13 days (7.3 %).

Assessment of THI effect on daily milk yield

The first analytical approach, based on RRM, revealed substantial heterogeneity within the analysed high-production group of our study. The herd-level coefficient for THI was estimated to be numerically positive (β=0.180±0.017 kg; t=10.41), and a similar trend was observed for THI_MAX (β=0.152±0.015 kg; t=10.25). Although this positive association may initially appear counterintuitive, a more detailed interpretation confirms that it reflects a temporal confounding between lactation-related changes in milk yield and the seasonal progression of environmental indices (Bernabucci et al., 2014).

DIM, representing the current stage of lactation in days, provided a significant explanation for the observed milk yield dynamics. As advancing lactation is naturally associated with a gradual decline in daily production, DIM was included in the model as a fixed continuous covariate to account for the lactation curve, showing a significant negative influence (β=-0.054±0.002; t=-30.59). The observed result is further explained by the distribution of 28,285 records across all lactation phases. Phase 1 (0-100 DIM), indicating the absolute peak of lactation (average milk yield 40.5 kg), contributed with 11,010 records, while Phase 2 (101-200 DIM), representing a sustained high-production phase associated with stable milk yield (average 38.6 kg), was the most prevalent with 14,318 records. Conversely, Phase 3 (>200 DIM), where milk production naturally declines (average 33.7 kg), was the least represented with only 3,877 records. The average DIM in the analysed period was 124.3 days, confirming that the herd was predominantly in its peak or stable high-production stages during the observation period.

Furthermore, the environmental context of the study likely contributed to the positive statistical momentum. The study started in March (average THI=51.06), a period already within the thermoneutral zone and characterized by relatively stable spring conditions without extreme winter cold. Due to the absence of very low THI values (winter data), the model lacked a cold-weather counterbalance. Since conditions remained largely stable and thermoneutral from March through early June, lactation-related increases in milk yield numerically outweighed the heat-induced declines recorded later during the summer window (NRC, 1971). While the aggregate herd-level trend remained numerically positive, this momentum was challenged during the summer months of the study when acute thermal stress events became more frequent.

Beyond these average effects, the RRM successfully identified that 53.8 % of individual cows in the THI_MAX model actually suffered significant yield losses as indices increased, with individual slope deviations ranging from -0.447 to +0.628 kg (Figure 1). Notably, the strong negative correlation between intercepts and slopes (r=-0.93) confirmed that cows with higher production exhibit stronger negative responses to rising THI indices. Such findings imply that high-production cows, specifically those at the peak or in stable production phases, were the most sensitive to heat stress. This pattern aligns with the biological mechanism linking elevated milk production to increased metabolic heat production and reduced thermoregulatory capacity (Kadzere et al., 2002; Berman, 2005; Lambertz et al., 2014).

These findings highlight substantial individual variability in heat stress resilience. The discrepancy between the marginal R² (0.063) and conditional R² (0.304) in the THI_MAX model suggests that fixed herd-level effects explained only a limited proportion of the variance, whereas cow-specific random effects accounted for a substantial additional part of the variability in milk yield (over 30 %). This supports the presence of marked individual differences in both baseline milk production and sensitivity to increasing THI_MAX values.

To further investigate individual response patterns across different exposure levels, the following section applies a lagged heat-stress approach to quantify delayed yield declines following acute thermal spikes.

Figure 1. Distribution of cow-specific responses in daily milk yield to THI and THI_MAX

Temporal response dynamics

The second analytical approach, LHR, based on lagged heat-stress response modelling, was used to analyse the temporal dynamics of daily milk yield response following identified acute heat stress episodes. These episodes were defined by two criteria: the corresponding index (THI or THI_MAX) had to exceed a predefined threshold while simultaneously exhibiting an acute increase (spike) of at least 2 units compared to the previous day. The applied approach, based on the methodological framework of Ravagnolo and Misztal (2000) and Brügemann et al. (2012), was chosen to isolate the physiological impact of sudden heat shocks from general seasonal trends, ensuring that each identified event represented a distinct period of acute stress. Production deviations were analysed within a dynamic temporal window spanning from three days before the event (Day -3) to five days post-exposure (Day +5), using the day prior to the shock as the reference baseline for all comparisons. Results of this analysis are summarized in Table 1 and Figures 2 and 3, which visualize daily milk yield response dynamics across relative days.

Table 1. Temporal dynamics of daily milk yield changes (kg) following heat stress events at different THI and THI_MAX thresholds

Cow variance = cow random effect variance (kg²); Event variance = heat stress event random effect variance; Residual SD = residual standard deviation; Maximum decline (kg) = maximum statistically significant milk yield decline; *p<0,05.

The analysis indicated no immediate negative production response on the day of the thermal event (Day 0) across any of the analysed thresholds. At the peak exposure level (THI_MAX ≥ 80), a temporary increase of +0.53 kg was observed on Day 0, as shown in Table 1. However, this slight increase indicates that the biological response to acute heat stress is delayed rather than instantaneous, reflecting temporary metabolic adaptation before the actual onset of production decline. Values recorded on Day -3 and Day -2 provided a basis for assessing baseline production stability and model validity. In most threshold scenarios, pre-event deviations remained close to zero and were not statistically significant, such as 0.01 kg and -0.01 kg in the THI ≥ 66 model. In contrast, significant negative deviations were detected at the THI_MAX≥72 threshold on Day -3 (-0.53 kg) and Day -2 (-0.38kg). This response suggests a potential pre-exposure effect, as described by Ouellet et al. (2019), implying that animals exposed to lower stress thresholds may have already experienced unfavourable environmental conditions before the identified two-unit THI spike. Moreover, pre-event values at the most extreme threshold (THI_MAX≥80) were positive and non-significant (0.50 and 0.06 kg), indicating that these animals entered the acute stress period from a period of stable high productivity.

Results across average THI thresholds revealed a progressive increase in the negative impact of heat stress as threshold intensity increased. At the THI≥66 threshold (n=22), a negative trend was observed from Day +3 (-0.37 kg) and became statistically significant on Day +4 at -0.50 kg (p<0.05). Production losses intensified at the THI≥70 threshold (n=16), reaching -0.76 kg (p<0.05) on Day +4. Maximum production loss among THI models was observed at the THI≥74 threshold (n=6), where milk yield declined by -1.08 kg on Day +4. Across all THI levels, production decline consistently occurred four days after exposure, indicating that higher thermal loads resulted in more substantial and prolonged production losses. Further evidence was provided by the continued significant depression on Day +5 (-0.92 kg) in the THI≥74 model.

Figure 2. Milk yield response to heat stress events at different THI thresholds (THI = 66, 70 and 74)

Comparable, but more acute, response patterns were observed for THI_MAX thresholds. According to Ouellet et al. (2019), THI_MAX values more precisely capture short-term peak heat intensity than average THI values. Negative responses at the THI_MAX≥72 model developed progressively, reaching -0.59 kg on Day +4 and -0.66 kg on Day +5 (p<0.05). At higher thresholds, production losses became more pronounced, with maximum declines reaching -0.97 kg at THI_MAX≥76 and -1.07 kg at THI_MAX≥80 on Day +4. Unlike average THI thresholds, the THI_MAX models showed earlier indications of production decline, including a significant reduction already observed on Day +1 (-0.34 kg) at THI_MAX≥80. The persistence of significant reductions on Day +5, including -0.77 kg at THI_MAX≥76 and -0.60 kg at THI_MAX≥80, additionally indicates incomplete recovery following high-intensity acute heat stress. Overall, the findings indicate that more intense thermal exposure was associated with greater milk yield losses and a slower return to pre-event production levels.

Figure 3. Milk yield response to heat stress events at different THI thresholds (THI_MAX = 72, 76 and 80)

Across all models, the maximum statistically significant production decline consistently occurred four days after the thermal exposure (Day +4). The observed delay is slightly longer than the 1-to-3-day lag responses reported by Gan Li et al. (2021) in similar dynamic temporal analyses following peak thermal exposure. Physiologically, this 4-day lag effect appears to reflect the time required for heat-induced disruptions to fully impair milk synthesis, consistent with findings by Bernabucci et al. (2014). Although the thermal load is immediate, subsequent reductions in dry matter intake (DMI) and alterations in endocrine signalling, as discussed by Tao et al. (2018) and Morales- Piñeyrúa et al. (2022), require several days to fully manifest in mammary gland metabolism and milk production.

Variance components presented in Table 1, revealed a progressive increase in individual animal variability at higher heat-stress intensities, providing important insights into individual resilience. Cow variance increased from 1.094 kg² (SD=1.05 kg) at THI≥66 to 3.902 (SD=1.98kg) at THI≥74, and further to 5.151 kg² (SD=2.27 kg) at the THI_MAX ≥ 80 threshold. Conversely, event-level variance declined noticeably at higher thresholds, decreasing from 0.706 kg² (SD=0.8402) to as low as 0.021 kg² (SD=0.14 kg) in the THI_MAX≥80 model. Event-level variance reflects the proportion of variation associated with the environmental conditions of individual heat events, such as event duration or nocturnal cooling. The results indicate that under extreme thermal conditions, differences among individual heat events become less determinant than the biological response of the animal itself. The marked increase in cow-specific variability under intense heat stress indicates that differences in thermal tolerance become the dominant source of variation, highlighting the potential for selecting cows with enhanced phenotypic resilience to acute heat stress.

Despite the high resolution of the individual records and the application of an advanced dynamic framework, several limitations of this study should be acknowledged. The analysis was conducted on a single commercial farm over a six-month period, primarily due to the predefined research timeline and technical limitations of the research project. The focus on a specific group of 169 high-producing cows enabled the use of high-precision monitoring that was not available for the entire herd. Although data were collected from a single herd, minimizing the herd-year-season effect and improving environmental consistency, the findings may have limited applicability across different management systems. Replication across multiple herds and varied production environments is recommended to validate and extend these results. Nevertheless, the high-resolution individual-level data from this study, provide a valuable foundation for future multi-herd studies and for the development of more targeted strategies to improve heat-stress resilience in high-producing dairy cows.