Comparative analysis of thermal indices for modeling cold and heat stress in US dairy systems

Quantifying the impact of thermal stress on milk yields is essential to effectively manage present and future risks in dairy systems. Despite the existence of numerous heat indices designed to communicate stress thresholds, little information is available regarding the accuracy of different indices in estimating milk yield losses from both cold and heat stress at large spatio-temporal scales. To address this gap, we comparatively analyzed the performance of existing thermal indices in capturing US milk yield response to both cold and heat stress at the national scale. We selected four commonly used thermal indices: the Temperature and Humidity Index (THI), Black Globe Humidity Index (BGHI), Adjusted Temperature and Humidity Index (THIadj), and Comprehensive Climate Index (CCI). Using a statistical panel regression model with observational and reanalysis weather data from 1981–2020, we systematically compared the patterns of yield sensitivities and statistical performance of the four indices. We found that the US state-level milk yield variability was better explained by the THIadj and CCI, which combine the effects of temperature, humidity, wind, and solar radiation. Our analysis also reveals a continuous and nonlinear responses of milk yields to a range of cold to heat stress across all four indices. This implies that solely relying on fixed thresholds of these indices to model milk yield changes may be insufficient to capture cumulative thermal stress. Cold extremes reduced milk yields comparably to those impacted by heat extremes on the national scale. Additionally, we found large spatial variability in milk yield sensitivities, implying further limitations to the use of fixed thresholds across locations. Moreover, we found decreased yield sensitivity to thermal stress in the most recent two decades, suggesting adaptive changes in management to reduce weather-related risks.


INTRODUCTION
Understanding the effects of thermal stress on milk yields is crucial to properly manage risks in dairy systems and ensure food security (Ji et al., 2020;Herbut et al., 2018;Wang et al., 2018).Thermal stress occurs when changes in weather push the thermal conditions of dairy cows beyond their thermoneutral zone, leading to not only physiological adjustments for temperature regulation (Rusche, 2020;Collier et al., 2019;Gaughan et al., 2009) but also to reduced milk production (Souza et al., 2023;Gisbert-Queral et al., 2021;Bohmanova et al., 2007;West, 2003;Kadzere et al., 2002) and exacerbated health and reproduction problems (Becker et al., 2020;Collier et al., 2019;Kadzere et al., 2002).While cows can acclimatize to thermal stress to a certain extent, the thermoneutral zones of the cows vary based on animal (e.g., breed, age, hair coat, level of production, and body condition), meteorological (e.g., temperature and humidity), and environmental management (e.g., barns, fans, and wind breakers) factors (Becker et al., 2020;Collier et al., 2019;Kadzere et al., 2002).Decades of selective breeding for high milk yield productivity has rendered dairy cows highly susceptible to heat stress, primarily due to increased metabolic heat production (Becker et al., 2020;Collier et al., 2019;Zimbelman et al., 2009;Kadzere et al., 2002).This challenge is significant for the US, the world's second-largest producer of milk and milk products (FAO, 2022), where economic losses due to heat stress are estimated to be approximately one billion dollars annually, despite the implementation of heat abatement strategies (Hutchins et al., 2023;Gunn et al., 2019;Ferreira et al., 2016;Key and Sneeringer, 2014;St-Pierre et al., 2003).Given the predictions of rising temperatures across the US (US-GCRP, 2023), an improved understanding of thermal impacts on dairy systems across the US is vital for enhancing climate resilience via more effective adaptation strategies.
The development of thermal indices for dairy cows has been critical for identifying stress responses, yet the diversity of indices, different versions, and different proposed thresholds for the indices can lead to confusion about what approach to utilize.The first thermal index in dairy systems was the Temperature Humidity Index (THI), which was adopted from the human discomfort index (Thom, 1959).Since the development of the THI, more comprehensive thermal indices have been proposed to incorporate consideration of additional weather variables, including wind speed and radiation flux.These improvements involve modifying the THI and refining apparent ambient temperature (i.e., dry-bulb temperature).The Black Globe Temperature Index (BGHI; Buffington et al., 1981), adjusted THI (THIadj; Mader et al., 2006), and the Comprehensive Climate Index (CCI; Mader et al., 2010) are major other thermal indices that incorporate not only temperature and humidity but also wind speed and radiation flux.The BGHI and THIadj were developed by substituting the dry-bulb temperature of the THI (Yousef, 1985) with black-globe temperatures (which are sensitive to solar radiation) and adding wind speed and solar flux into the THI (NRC, 1971), respectively.On the other hand, the CCI was designed to adjust dry-bulb temperature with humidity, wind speed, and radiation flux data.
The evolution of developing thermal indices has been continuing, with the THI remaining the one of the most widely used indices (Ji et al., 2020;Herbut et al., 2018;Wang et al., 2018).Nevertheless, the presence of 7 distinct versions of the THI complicates the establishment of consistent thermal assessments in the US.Each version, authored by different individuals or groups, offers varying heat stress thresholds, spanning from 68 to 83 (Wang et al., 2018;Bohmanova et al., 2007).Moreover, when utilizing a specific version of the THI (NRC, 1971), different studies have provided varying THI threshold values: 68 (Zimbelman et al., 2009), 72 (Bohmanova et al., 2007;Ravagnolo and Misztal, 2000), and 74 (Bohmanova et al., 2007).This is because the majority of these studies focus on animal, on-farm, or local scales and only at a few select points in time.The application of a certain index or stress threshold requires careful consideration of its experimental conditions including study periods (e.g., summer vs. all seasons), climate regions (e.g., temperate vs. tropics), physiological capacity of cattle (e.g., breed, level of production), and environmental management practices (e.g., barns vs. pasture) (Wang et al., 2018).Hence, a comprehensive evaluation of existing indices in estimat-ing the effects of thermal stress on milk yields is crucial to reduce inconsistencies in their applications.
A large knowledge gap also remains about the accuracy of different indices in capturing the effects of weather on milk yield at large spatio-temporal scales.The effectiveness of thermal indices in explaining milk yield changes at large spatio-temporal scales is influenced by the range of environmental management practices employed by farmers; regions with advanced housing and management would typically exhibit reduced thermal stress sensitivity (Becker and Stone, 2020;Hill and Wall, 2015).Yet, the reliability of applying animal-based thermal indices at state or national dairy sectors is questionable primarily due to their derivations from conditions specific to individual animals and farms.A study using the THI in a regression model accounting for local and time-varying environmental management practices (Gisbert-Queral et al., 2021) found US milk yield sensitivities to heat stress to be substantially lower than those reported in animalbased studies (Mauger et al., 2015).Furthermore, Gisbert-Queral et al., 2021 highlighted the similarly detrimental effects of both cold and heat stress on milk yields, underscoring the adverse effects of cold stress driven by high energy demands for body heat production (Rusche, 2020;Broucek et al., 1991;Young, 1983;MacDonald and Bell, 1958).Therefore, it is important to understand how different thermal indices estimate milk yield changes to both cold and heat stress at large spaito-temporal scales and what thermal index most accurately captures this response.We hypothesize that thermal indices considering a variety of weather variables would more accurately predict changes in milk yield due to both cold and heat stress.
This study aims to comprehensively evaluate how existing thermal indices estimate the response of US milk yields to both cold and heat stress from 1981 to 2020.Four major thermal indices were included in our analysis -THI, BGHI, THIadj, and CCI.By employing a statistical panel regression model with 2 high quality and bias-corrected weather data sets across 36 states, we quantify milk yield sensitivities to a range of cold to heat stress at different spatio-temporal scales.Specifically, the objectives of our study are (1) to assess the magnitude and patterns of milk yield sensitivity to both cold and heat stress at national, regional, and multidecadal levels; and (2) to quantify the statistical performance of the 4 indices, evaluating and identifying the most adequate index for large spatio-temporal scale assessments.

DATA
Our analysis examines daily milk yield across 36 states in the US from 1981 to 2020.We utilized publicly available milk yield data from the USDA National Agricultural Statistics Service (NASS; (USDA, 2022), which compiles monthly average milk production per cow at state level.Each state had different number of yield survey data, so our analysis was based on an unbalanced data set (Figure 1A).The daily milk yield was computed by dividing the monthly average milk yield by the number of days in each month to ensure comparability of milk productivity across time.The average daily milk yield over the available number of years per state are depicted in Figure 1B.Countylevel dairy cow population data were acquired from the USDA-NASS for the recent 5 census years (USDA, 2022(USDA, ): 1997(USDA, , 2002(USDA, , 2007(USDA, , 2012(USDA, , and 2017.Over time, national average milk yields (across millions of dairy cows) substantially increased throughout the US, from ~15.5 kg/head/day in 1981 to ~29.5 kg/head/day in 2020 (Figure 1C).
Our analysis uses daily weather data from 2 sources: PRISM (Oregon State Univ., 2022) and AgERA5 (Copernicus Climate Change Service, 2019) (Table 1).We chose these 2 data sets because they are highquality, bias-corrected surface level data available for 1981-2020.The PRISM data set is driven by a wide range of in situ observations and continuous quality control measures and exhibits fine spatial (4km) and temporal (daily) resolutions (Oregon State Univ., 2022).The AgERA5 data set is built upon the hourly ECMWF reanalysis of ERA5 data interpolated by the operational high-resolution atmospheric model at a 0.1° spatial resolution (Boogaard et al., 2022).The PRISM data set has been widely used in agricultural studies (Choi et al., 2023;Gisbert-Queral et al., 2021;Butler et al., 2018;Schlenker and Roberts, 2009) and the Ag- ERA5 data set is specifically designed for agriculture and agro-ecological studies (Boogaard et al., 2022).We spatially averaged each of gridded PRISM and AgERA5 data into the county-level across all grid boxes within county boundaries of the 2017 US Census Bureau Map.
We employed a range of surface weather variables that can affect thermal environments of dairy farms.We obtained daily minimum (T min ) and maximum (T max ) air temperatures, minimum (VPD min ) and maximum (VPD max ) vapor pressure deficits, dew-point temperature (T dp ), and precipitation (PPT) data from the PRISM data set.Daily wind speed (V) and downward solar radiation (SSRD) data were from the AgERA5 data set, as these 2 variables are not provided by the PRISM data set (Table 1).To capture the surface level wind speed, we adjusted the AgERA5′s wind speed data at 10m to 2m based on the conversion equation of the Food and Agricultural Organization (FAO) (Allen et al., 1998).We also derived net solar radiation after calculating both net shortwave and longwave radiative fluxes along with elevation data from the 2010 Global Multi-resolution Terrain Elevation Data (Danielson and Gesch, 2011).Summary statistics of daily countylevel weather data to calculate the 4 thermal indices are shown in Table 2.

FOUR MAJOR THERMAL INDICES
Temperature Humidity Index (THI) The calculation of THI in this study was based on equation [1] (NRC, 1971): where T max is the maximum temperature and RH min is minimum relative humidity calculated based on the RH equation of Table 1 by using maximum vapor pressure deficit (VPD max ) from the PRISM dataset.We used T max and RH min for the THI calculation as it was found to be more effective in capturing heat stress compared to other combinations that employ minimum or average temperature and relative humidity (Ravagnolo et al., 2000).This THI version has been widely used in previous studies (Souza et al., 2023;Gisbert-Queral et al., 2021;Fodor et al., 2018;Hill and Wall, 2015;Dunn et al., 2014;Brügemann et al., 2012;Zimbelman et al., 2009;Bohmanova et al., 2007;St-Pierre et al., 2003;Ravagnolo et al., 2000).However, apart from its optimal application to outdoor cattle (Zimbelman et al., Choi et al.: Comparative analysis of thermal… where SVP is saturated vapor pressure and defined as SVP e  Reference: (Allen et al., 1998).2009; Bohmanova et al., 2007), no information about the development of equation [1] is available (Ji et al., 2020;Wang et al., 2018;Dikmen and Hansen, 2009) (Table 3).Furthermore, the THI does not take into account the effects of wind speed and radiation (Wang et al., 2018) (Table 3).This potentially becomes a limitation of using the THI for dark-colored cattle in certain regions where wind speed or solar radiation can be important in thermal exchange between the cows and their surroundings (Ji et al., 2020;Wang et al., 2018;Silva et al., 2007;Mader et al., 2006).For example, higher wind speeds can reduce heat stress via the convective cooling of cows when their body temperature is higher than the ambient temperature (Becker and Stone, 2020;Ji et al., 2020;Mader et al., 2006;Kadzere et al., 2002).Conversely, higher solar radiation can increase the heat load of cows, particularly for black-colored cattle in tropical regions, open lots, or poorly insulated barns (Ji et al., 2020;Wang et al., 2018;Herbut and Angrecka, 2013;Silva et al., 2007;Mader et al., 2006).Among seven different versions of the THI, our study refers to the value of 72 as the THI stress threshold, which many studies have reported as the heat stress cutoff for reducing milk yields (Wang et al., 2018).

Black Globel Humidity Index (BGHI)
The calculation of BGHI in this study was based on equation [2] (Gaughan et al., 2002;Buffington et al., 1981): ., In equation [2], T dp is average dew-point temperature from the PRISM dataset and T bg is the black globe temperature calculated by using average air temperature, T avg , from the PRISM dataset and net solar radiation, R n , from the AgERA5 dataset.The unit of R n was converted from W m 2 to cal cm −2 min −1 .Given that the BGHI was developed by replacing the drybulb temperature of the THI (Yousef, 1985) with a field-measured black globe temperature, T bg needs to be measured by a black globe temperature sensor to accurately capture the combined effects of air temperature, convective cooling of wind, and solar radiation.However, it is not feasible to obtain direct T bg measurements across the US.We therefore referred to Gaughan et al., (2002) to calculate T bg as stated in equation [2].This derived T bg does not account for the effect of wind speed (Ji et al., 2020;Wang et al., 2018).It also cannot take negative air temperatures mathematically, thus we set all temperatures below 1 °C as a constant value of 1 °C.Nevertheless, we included the BGHI in our analysis because it is the only index that accounts for the impact of thermal radiation when assessing changes in milk yields due to heat stress (Table 3).Our study refers to the value of 70 as the BGHI stress threshold (Buffington et al., 1981) (Table 3).

Adjusted Temperature Humidity Index (THIadj)
The calculation of THIadj in this study was based on equation [3] (Mader et al., 2006): where THI is from equation [1] by using the PRISM dataset, and V and SSRD are mean wind speed at 2m and downward solar radiation from the AgERA5 dataset, respectively.The THIadj advances the THI by considering the effects of V and SSRD on panting scores (Table 3).While increasing V has a chilling effect on the THI, increasing SSRD has the heating effect on it (Wang et al., 2018;Mader et al., 2006).
Its development is based on data collected during the summertime in Nebraska between 2 and 5 pm, over a span of two years (Table 3).Our study refers to the value of 74 to 79 as the THIadj stress threshold (Mader et al., 2006) (Table 3).

Comprehensive Climate Index (CCI)
The calculation of CCI in this study was based on equation [4]: where  .
In equation [4], T a , RH adj , V adj , and RAD adj indicate the average air temperature, adjusted relative humidity, adjusted wind speed, and adjusted solar radiation, respectively.We used the data on T avg and RH from the PRISM dataset, and V and SSRD from the Ag-ERA5 dataset.Unlike other selected indices, the CCI enables measuring both cold and heat stress for an air temperature range of -30°C to +45°C (Table 3).It was built from previous thermal indices: the THIadj (Mader et al., 2006) and Heat Load Index (Gaughan et al., 2008) for hot conditions and Windchill Index (Ames and Insley, 1975) for cold conditions.The final continuous forms of combined hot and cold conditions were developed based on 15 years of the weather station data in Nebraska (Mader et al., 2010).However, it mainly adjusts air temperatures affected by the individual effect of RH, V, and SSRD.It does not consider their interactive effects such as wind and water vapor pressure in determining evaporative cooling (Kadzere et al., 2002).Our study refers to the value of 0 to -10 as the CCI cold stress and the value of 25 to 30 as the CCI heat stress threshold (Mader et al., 2010) (Table 3).

Cow population weighted thermal indices
To capture the potential nonlinear effects of thermal stress on milk yields, we computed the exposure time of cows to various levels of each index in different states.Each index was categorized into bins with a range of 4.9 to 5 units, using the 8th and 92nd percentiles as the boundaries for the edge bins.Since the indices were computed at daily county-level, they were aggregated into monthly state-level in each bin to serve as a cumulative measure of exposure (Schlenker and Roberts, 2009).To give greater consideration to counties with a large number of cows over those with a small or nonexistent cow population in each state, we used the average cow population data of each county, averaged over the USDA's five census years, as the weight in equation [6].For each thermal stress index, the weighted exposure to index bin j in state i in month m of year t is defined as X j,i,m,t , calculated as: where E j,n,i,m,t is the number of days on which cows in county n of state i in month m of year t exposed to index bin j and D m,t is the number of days in month m of year t.w n,j is the weight of county n in state i, calculated using equation [6] (Gisbert-Queral et al., 2021).
where C n,i is the average number of cows for county n in state i over the recent five agricultural census years and C k,i are the sum of the average number of cows in county k for state i.
In our analysis, we referred to the smallest index bin as the cold extreme and the largest index bin as the heat extreme.

Statistical Panel Regression Model
To quantify the marginal response of milk yields to a range of thermal indices (i.e., cold to heat stress), we applied the statistical panel regression model from Gisbert-Queral et al., 2021 (equation [7]).Panel regression models have been widely used to assess the relationship between weather and agricultural yields (Choi et al., 2023;Zhu et al., 2022;Lobell et al., 2020;Butler et al., 2018;Schlenker and Roberts, 2009).These models utilize fixed effects to control for unobserved, time-invariant characteristics of each panel unit (here, US states) as well as time-varying factors (such as technological advancements or average seasonality), thus allowing a precise identification of the impacts of weather variability on yields.Given that the long-term and seasonal variations in milk yields can be influenced by factors such as technical improvements, calving cycles, and forage availability in addition to thermal stress, our model estimates the impact of thermal stress on milk yields given those non-climatic local and time varying effects.For each thermal index, we fit the following panel regression model: where ) is the natural logarithm of average milk yield expressed in kilograms per cow per day in state i in month m of year t, where average monthly milk yield is divided by the number of days per month of each year.X j,i,m,t indicates the weighted exposure to an index bin j in state i in month m of year t, calculated from equation [5], with the corresponding coefficient vector β j .We excluded one bin out of the index bins as the 'optimal milk yield bin' to avoid perfect multi-collinearity of the explanatory variables.By multiplying β j by 100, the coefficient vector represents the percentage change of the average milk yield per cow for an additional day at a given index bin j, relative to the optimal condition.γ i,t are state-year fixed effects that control for factors changing yearly per state.Specifically, the state-year fixed effects include technical improvements (e.g., breeding, milking systems, housing, sprinklers, etc.) (Becker and Stone, 2020;Hill and Wall, 2015) and other year-to-year variations (e.g., feed price and policy) (Hahn, 2022) that may influence milk productivity.
Similarly, δ i,m are state-month fixed effects that control for the average seasonal cycle in each state.These statemonth fixed effects include calving cycles, forage availability, forage quality, (National Academies of Sciences and Medicine, 2021; Ritz et al., 2021;Ferreira et al., 2020) and other features (e.g., disease) (Fleischer et al., 2001) that could possibly affect seasonal variations in milk productivity.Note that forage availability, forage quality, and disease can be affected by thermal stress, thus possibly changing milk yield (Becker et al., 2020;Kadzere et al., 2002).The impacts of the state-year fixed effects can be inferred from yearly trends of daily milk yields in Figure 1C, while the impacts of the state-month fixed effects can be inferred from monthly cycles of detrended daily milk yield in Figure 2. ε i,m,t are an error term.Our approach does not require any data related to the fixed effects, and these dummy variables capture changes in technology and management for which we have limited data across space and time.It also does not differentiate the influence of parity (e.g., cows with only one birth vs. multiple births) on milk yields (Hutchins et al., 2023;Campos et al., 2022;Becker et al., 2020), as milk yield data in our study are state-level, averaged across all herds.Additionally, our model does not consider time-lagged effects of animal responses or consecutive days of exposure to thermal stress (Souza et al., 2023;West et al., 2003).The pattern of yield sensitivities remains robust to the number of bins and different percentile values to the edge bins.Confidence

Model performance evaluation metrics
To determine the effectiveness of the 4 major indices in quantifying milk yield variability to thermal stress at large spatio-temporal scales, we used 3 goodness of fit metrics: adjusted R 2 , within adjusted R 2 (Berge, 2018) and Bayesian Information Criterion (BIC) (Schwarz, 1978).The adjusted R 2 and within adjusted R 2 are used to estimate how much of the variance in the response variable (i.e., milk yields) is explained by predictors.While the adjusted R 2 examines the explanatory power of all predictors (e.g., index bins and fixed effect variables in equation [7]), the within adjusted R 2 evaluates the explanatory power of predictors that are detrended by fixed effect variables (Berge, 2018).In other words, the within adjusted R 2 computes the variance of the milk yields explained exclusively by thermal index predictors.In both R 2 metrics, the values span from 0 (i.e., worst fit) to 1 (i.e., best fit).Additionally, the BIC represents a penalized-likelihood information criterion that panelizes more complex models when comparing different models (Schwarz, 1978).Consequently, lower BIC values suggest a better model performance.

Evaluation of Thermal Stress Indices in Mod-
eling Milk Yield Changes.Our statistical models reveal that both cold and heat stress have significant and deleterious effects (P < 0.05) on US milk yields across all 4 indices (Figure 3).The damaging effect of heat stress on milk yields has been extensively discussed in prior studies (Souza et al., 2023;Collier et al., 2019;Gunn et al., 2019;Zimbelman et al., 2009;Bohmanova et al., 2007;West et al., 2003).Our findings demonstrate that US milk yields were similarly compromised by both cold and heat stress.Specifically, cold extremes (i.e., the bin containing the smallest thermal index value) reduced milk yields by approximately 6%, which was comparable to the impact of heat extremes (i.e., the bin containing the largest thermal index value).While the CCI is the only index that was explicitly formulated to cover the range of effects from cold to heat stress (Table 3), our analysis reveals that the other 3 indices also detected negative effects on milk yields from cold extremes.This result may reflect adverse effects of extreme cold temperatures on milk yield, primarily due to the diversion of energy from production to maintenance functions (Broucek et al., 1991;Young, 1983;MacDonald and Bell, 1958).By including year-round weather conditions into the 4 indices, we found that all 4 indices show a greater number of cold days (i.e., bins to the left of optimal thermal condition) than heat days (i.e., bins to the right of optimal thermal condition) across the US.
Despite similar yield sensitivities to cold and heat extremes, milk yield losses emerged more rapidly from increasing heat stress than from cold stress (Figure 3).We observed a larger gradient of milk yield losses as the thermal index increased toward the heat extreme.In contrast, the gradient of yield losses remains gradual and consistent as the thermal index decreased toward the cold extreme.Across all 4 indices, the steepest gradient in yield losses occurred in the maximum heat extreme condition compared with the adjacent thermal bin.For instance, the THI model exhibited a notable change (about -4%) in milk yield losses between the midpoint of bins 78 and 87 (Figure 3A).This pronounced sensitivity gradient in increasing heat stress might be associated with the vulnerability of highproducing lactating dairy cows to heat stress, possibly due to elevated metabolic heat output (Becker et al., 2020;Zimbelman et al., 2009;Kadzere et al., 2002) and decreased feed intake (Becker et al., 2020;Collier et al., 2019).
The implied optimal temperature conditions vary by thermal index (Figure 3).For each index, we identified the optimal thermal index bin, which is associated with the highest milk yield on a national average in the historical record (Methods).For the THI and CCI, optimal yield conditions are bins with a midpoint of around 19°C (approximately 67°F).In contrast, the BGHI and THIadj suggest optimal yield conditions at around 15.5°C (approximately 60°F).Previous studies have reported the thermoneutral zone at a THI of 57 (Collier et al., 2019) or within a temperature range of 5 to 25°C (Herbut et al., 2018;Kadzere et al., 2002).Our estimated optimal conditions may differ from these animal-level thermoneutral zones, since our approach incorporates the mediating effects of environmental management practices (e.g., housing, heating, shading, etc.) on the relationship between weather and milk yields (Methods).Furthermore, the optimal yield conditions in our study represent national and climatological average conditions rather than animal-specific conditions.
Overall, we observed continuous and nonlinear responses of milk yields to a range of cold to heat extremes for all 4 indices (Figure 3).We found that this trend diverges from the typical pattern of crop yield losses emerging abruptly above or below certain threshold levels (Schlenker and Roberts, 2009).Thresholds in ecological systems are generally defined as points where sudden changes occur in response to external or internal pressures or disturbances (Groffman et al., 2006).However, in our analysis, the milk yield sensitivities showed a gradual decrease from optimal to both cold and heat stress conditions, with no obvious threshold effect, although we noted some acceleration in losses for the most heat extremes.These continuously decreasing patterns of milk yield responses can consistently be found in animal-level studies using on-farm or nearby weather data for heat stress (Campos et al., 2022;Brügemann et al., 2012;Collier et al., 2012;Zimbelman et al., 2009) and cold stress (Fu et al., 2022;MacDonald and Bell, 1958).Yet, many previous studies have relied on fixed thresholds to define heat stress for dairy cows, such as 70 (Gunn et al., 2019;Dunn et al., 2014;St-Pierre et al., 2003) or 72 (Hutchins et al., 2023;Collier et al., 2019;Bouraoui et al., 2002).Our finding suggests that the use of such fixed stress thresholds may not effectively capture the complex dynamics of thermal loads on milk yield responses, as compared with the use of a continuous yield response model.
Additionally, we found limited applicability of the 4 indices to estimating national milk yield losses from thermal stress which is beyond their experimental conditions (Table 2 and Table 3).Notably, the BGHI is the sole index to estimate milk yield responses to thermal stress but is unable to take temperatures below 1°C (Equation [2]).This probably contributed to its lower explanatory power (BIC and adjusted within R 2 ) among the 4 indices (Table 4).The THI also exhibited poor performance in modeling milk yield sensitivities, possibly because it only considers the effect of temperature and humidity as thermal stress (Table 3).Conversely, the THIadj and CCI, which take account into additional factors such as wind speed and solar radiation, showed improved performance.Nevertheless, they were not designed to predict milk yield response, rather they were designed based on panting scores and/ or dry matter intake (Table 3).Despite positive correlations between physiological responses (e.g., panting scores) and milk yields, these can be transient indicators of heat stress before milk yield reductions (Thulani et al., 2019;Chang-Fung-Martel et al., 2017).Furthermore, none of the indices were developed using diverse sets of weather variables (e.g., vapor pressure deficit, precipitation) across various climate regimens.This potentially suggests the need to designing a complementary index, calibrated specifically to milk yield, that captures variation across a wide range of environmental conditions.
Our analysis also shows that a great amount of milk yield variability was influenced by non-climatic factors such as technological advancements and calving cycles.We accounted for these non-climatic factors as statespecific year (e.g., technological advance) and month (e.g., calving cycles, feed availability, and forage quality) fixed effects in our statistical model (Methods).All models including these fixed effects exhibited 0.99 of adjusted R 2 (Table 4).However, the ability of each index to solely explain the variance of milk yields through thermal conditions (within adjusted R 2 ) was reduced to around 0.03 (Table 4), indicating that the non-climatic variables maximized their goodness of fit.Yet, this does not necessarily imply that cows are insensitive to thermal stress.Rather, it is more likely attributable to the large influence of the non-climatic factors coupled with our conservative approach in accounting for them through year-to-year and month-to-month variations in each state.For example, the state-specific year fixed effects in our model captures the near-doubling of milk yields since 1981 related to technological progress, including improved genetics and environmental management practices (Becker and Stone, 2020;Hill and Wall, 2015) (Figure 1C).Similarly, the state-specific month effects control for average seasonal variability in milk yields related to calving cycles and forage availability and quality (National Academies of Sciences and Medicine, 2021;Ritz et al., 2021;Ferreira et al., 2020)  ure 2).Note, however, that changes in environmental management practices can also alter the sensitivity of milk yields to weather in general, a topic we investigate by examining spatial temporal differences in yield sensitivities in the following sections.Further, changes in forage quantity and quality can also be a mechanism for milk yield sensitivity to thermal stress, in addition to direct organismal response.This mechanism should be captured in our thermal stress coefficients, whereas the state-specific month effects will capture any other average patterns of seasonality.
Quantification of Spatial Differences in Yield Sensitivities.We found evidence of large variations in yield sensitivities to thermal stress between cool and warm states for all 4 indices (Figure 4).Consistent with a previous study (Gisbert-Queral et al., 2021), warm states tended to be more sensitive to cold stress and less sensitive to heat stress than cool states.As the thermal index decreased toward the cold extreme, the difference in milk yield sensitivities between the 2 regions increased.In contrast, the divergence in the milk yield sensitivities decreased as the thermal index increased toward the heat extreme.The large divergence in increasing cold stress may arise from different biological and management factors in the 2 regions.For example, cows can adjust or acclimate to cold climates by thickening their hairs or increasing fat layers (Rusche, 2020;Collier et al., 2019;Young, 1983).Access to barns or dry bedding can also help to lessen the adverse effects of cold stress on milk yields in the cool states (Fu et al., 2022;Rusche, 2020).On the other hand, the narrower divergence in increasing heat stress between the 2 regions might be attributed to the selective breeding of high milk producing cows across the US or higher percentage of operations with heat abatement systems in warm states than cool states (USDA, 2022).Our analysis shows that the geographic variation in milk yield response to thermal stress is robust across thermal indices.
The spatial variability in milk yield sensitivities in our results suggests that relying on fixed thresholds may lead to under-or over-estimation of milk yield losses from thermal stress at different locations (Figure 4).Regional-scale assessments of thermal stress for US dairy systems have been less explored than animalor local-scale assessments.While a few selected regional studies have estimated milk yield losses to heat stress, they relied on fixed index value as heat threshold across the US (e.g., THI >70) (Hutchins et al., 2023;Gunn et al., 2019;Mauger et al., 2015;Key and Sneeringer, 2014).Failing to account for region-specific threshold variations at large spatial scales could complicate the reliability and accuracy of using existing thermal indices for milk yield assessment.This has been highlighted by studies that found substantial variations in milk yield losses based on the choice of a THI threshold in different locations (Fodor et al., 2018;Bohmanova et al., 2007).Therefore, our finding of the yield variability between the cool and warm states provides an insight into a thorough evaluation of the indices at different spatial scales (e.g., region) for decision-makers to minimize spatial bias in yield assessments.
Our models of the 2 regions better explained milk yield variations in thermal stress compared with the national-level model (Table 4 and Table 5).Except for the BGHI, the 3 indices exhibited more than 20% increased performance (adjusted within R 2 ) for the warm states while the cool states remained similar or had lower performance levels of the national level.The warm state models may have higher explanatory power because they are more akin to the conditions in which the indices were developed.As the national sample was divided into the 2 regions, we found large confidence intervals, particularly in extreme cold and heat conditions.Nevertheless, exposure to the cold and heat extremes in the 2 regions continued to be significantly detrimental to the milk yields.The THIadj and CCI appeared to have different optimal yield conditions in both regions compared with the national-level optimal condition.However, these were within the confidence interval range of the national-level yield sensitivities, thereby not showing significant differences.
Quantification of Temporal Differences in Yield Sensitivities.Distinct heterogeneity in yield sensitivities to thermal stress was observed between 1981-2000 (early period) and 2001-2020 (recent period), with the impacts of heat extremes remaining consistently detrimental in both periods (Figure 5).Overall, we found that the yield in the recent period appeared to be relatively less sensitive to thermal stress than the early period, consistent with earlier work based on the THI (Gisbert-Queral et al., 2021).Particularly, the yield losses to cold extreme during the recent period was large enough to be substantially diverged from the early period (Figure 5A, B, and D).Nevertheless, milk yield sensitivity remained nearly unchanged under heat extreme conditions regardless of the time periods across all 4 indices.For example, in the CCI model (Figure 5D), the difference in yield loss between the 2 periods was approximately 2.8% in the cold extreme condition (i.e., a midpoint of bin -23) while its difference in the heat extreme condition (i.e., a midpoint of bin 39) was about 0.8%.Our finding suggests the tem-poral dependence of milk yield sensitivities to thermal stress.
The overall pattern of decreased sensitivity to thermal stress in the recent period may be linked to adaptive management (Figure 5).Prior studies have suggested that sustaining increasing yield trends over time (Fig- ure 1C) could be attributed to evolving management practices of farmers via improved feed intake or/and efficiency (Becker et al., 2020;Thulani et al., 2019;Kadzere et al., 2002), the use of housing (USDA, 2016;Hill and Wall, 2015), microclimate modification (e.g.,  sprinkler and shade) (Gunn et al., 2019;USDA, 2016;West, 2003;D. E. Buffington et al., 1983), refined nutritional management (e.g., water, protein contents and nutrient quality) (Becker et al., 2020;Thulani et al., 2019;West, 2003), etc. Changes in feed sourcing could be another potential driver for smaller yield sensitivities in recent decades (Hahn, 2022).Future research could examine what contributed to reduced milk yield sensitivity over time to support adaptation efforts and avoid climate damages.
The statistical performance of the models for the 2 periods was also improved compared with the nationallevel model (Table 5).The increase in adjusted within R 2 was more apparent in the models of the recent period, particularly for the THI and BGHI.Similar to the models of the 2 regions, both cold and heat extremes significantly and negatively affected milk yields in both periods.Significant changes in the optimal yield condition in the recent period were observed only for the THI (i.e., a median value of 72).While the development of the BGHI and THIadj were based on 2 to 3 years of summer records, the utilization of 15-year weather data in constructing the CCI seemed more conducive to better milk yield predictions.However, we did not find statistical consistency in the superiority of the CCI (Table 4 and Table 5).This could be partly because the 15-year weather data was predominantly focused on Wisconsin (Table 3) and their range of different weather variables was not comprehensive enough to cover conditions relevant to other US states (Table 2).A robust understanding of thermal indices for evaluating milk yield changes to thermal stress across the US is needed to reliably project thermal impacts and inform effective adaptation strategies.

CONCLUSION
Our study fills the gap in understanding the effectiveness of 4 major thermal indices -THI, BGHI, THIadj, and CCI -in assessing milk yield change due to cold and heat stress across large spatio-temporal scales.We found that the THIadj and CCI better explained for changes in milk yield to a range of cold to heat stress at national, regional, and decadal scales than the THI and BGHI.For all 4 indices, our findings of continuous and nonlinear responses of milk yield to thermal stress suggest potential limitations of using fixed thresholds for estimating yield losses at large spatio-temporal scales.Our study also reveals large spatial variations in milk yield sensitivity between cool and warm states, emphasizing the importance of understanding the entire system and regionally-specific impact assessments.Additionally, we found suggestive evidence of adaptive management practices, as milk yield sensitivities to thermal stress decreased in the recent period (2000-2020) compared with the earlier period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000).We note that our results are not purely representative of the organismal response to thermal stress.By utilizing the panel regression model with state-level yield and year-round weather data, our study quantifies climatologically and nationally representative yield sensitivities to thermal stress after accounting for the management practices that farmers have used.Furthermore, our study provides important insights into the applicability of animal-based thermal indices for large-scale milk yield assessment and need to advance data on calving cycles, feed intake, and management to refine evaluations of thermal stress impacts for future risk projections.

Figure 1 .
Figure 1.Spatial and temporal patterns of average daily milk yields (kg/head) at state-level: (A) the number of available years of yield data that was used in our analysis; (B) daily average milk yields over the available years between 1981 and 2020; and (C) temporal patterns of average daily milk yields from 1981 to 2020 at individual states (gray line) and the national average (black line) weighted by the average cow population of each state over the 5 census years.States in white on the spatial maps (A and B) indicate no available data.
Choi et al.:  Comparative analysis of thermal… intervals are estimated by using heteroskedasticity robust standard errors.

Figure 2 .
Figure 2. Change in milk yield at monthly state-level.The change in milk yield (kg/head/day) was estimated by removing a quadratic yearly trend in average daily milk yields, then calculating the monthly averages of these detrended yields from 1981 to 2020 for each state (gray line).The lines were smoothed using a polynomial approximation.We highlighted 2 states with high milk yields (~30 kg/head/day) and large cow populations from the cool region in blue and warm region in red.

Figure 3 .
Figure 3. National average milk yield sensitivities to thermal conditions ranging from cold to heat extremes by using (A) THI, (B) BGHI, (C) THIadj, and (D) CCI.The filled black circles indicate our model-estimated milk yield sensitivities (i.e., coefficients) to the thermal conditions of each index.These sensitivities are interpreted as the percentage change of average milk yield per cow for an additional day at a given index bin, relative to the optimal yield bin (Methods).All yield sensitivities appear negative since they are relative to optimal conditions (denoted by the dashed zero line).The yield sensitivities are connected in solid lines to illustrate the approximated smooth response of milk yields.The literature-based thresholds of each index are denoted in red for heat stress and in blue for cold stress, based on their respective references (Table 3).The gray shading indicates 95% confidence intervals of our modeled yield response based on robust standard errors.The green histograms at the bottom of each panel show the number of days exposed to each thermal index bin.
Choi et al.: Comparative analysis of thermal…

Figure 4 .
Figure 4. Milk yield sensitivities to thermal conditions between cool and warm states by using (A) THI, (B) BGHI, (C) THIadj, and (D) CCI.Colored lines indicate differential yield sensitivities in blue for cool states and in red for warm states.Shaded areas represent 95% confidence intervals.

Table 1 .
Summary of weather variables used to calculate thermal indices avg Average air temperature which was calculated from averaging the sum of T min and T max T dp Average dew-point temperature °C PRISM VPD Average vapor pressure deficit which was calculated from averaging the sum of minimum (VPD min ) and maximum (VPD max ) vapor pressure deficits.

Table 2 .
Summary statistics of daily county-level weather data across 36 states for 1981-2020

Table 4 .
Choi et al.:Comparative analysis of thermal… Summary statistical scores of thermal indices

Table 5 .
Adjusted within R 2 for spatio-temporal analysis