Trends in the risk of heat stress to Canadian dairy cattle in a changing climate

Abstract Canada's climate is warming faster than the global average, but the warming is unevenly distributed. This study analyzes historical and future climate change in dairy-producing regions across Canada to better understand how Canada's dairy cows are affected. Historical changes (i.e., 1960–2019) were assessed using temperature and humidity data from 29 weather stations across the country. The temperature–humidity index (THI) was used as an indicator of dairy cattle at risk of heat stress, and three THI metrics evaluated the frequency, severity, and duration of potential heat stress. Future scenarios were investigated using five global climate models to project daily THI under three Shared Socioeconomic Pathways (SSPs). Projections were grouped into three time periods (2020–2049, 2040–2069, and 2060–2089). Historical climate trends show an increase in temperature, humidity, and THI exceedance in most west coast and eastern Canada locations, affecting 84% of the national dairy herd. Future scenarios project that 90% of the national herd will experience a large increase in the frequency, severity, and duration of THI exceedance under all but the most optimistic SSP. These findings highlight the need for Canadian dairy farmers to consider heat-stress adaptation strategies.


Introduction
Dairy cattle are sensitive to heat stress, which occurs when an animal's thermal balance cannot be maintained (Hansen 2020).Heat stress can negatively affect milk production (Wheelock et al. 2010) and milk composition (Campos et al. 2022a), reproduction (Schüller et al. 2014), and health of dairy cows (West 2003), resulting in lost revenue (St-Pierre et al. 2003).To study heat stress in dairy cows, the temperaturehumidity index (THI), an environmental indicator of heat stress, is commonly used.Multiple THI formulas have been developed over the past half-century, all of which use air temperature (T a ) in combination with a measure of humidity (relative humidity = RH, wet-bulb temperature, or dew-point temperature).Although the parameter weighting differs in each equation, there is overall agreement that THI increases as temperature and humidity increase (Wang et al. 2018).The most widely used THI formula uses T a and RH and was published by NRC (1971).Milk production decreases when THI increases above a certain threshold, although multiple thresholds have been identified.For example, Zimbelman et al. (2009) identified a THI threshold of 68.For the Canadian provinces of Ontario and Quebec, Campos et al. (2022a) identified the lower thresholds of 64, 50, and 58 for a decline in milk, fat, and protein yield, respectively, based on estimated outdoor THI.Considering all Canadian provinces, Rockett et al. (2023) identified thresholds as low as 47-50 for milk yield, 48-55 for fat yield, and 58-62 for protein yield.In addition to intensity of heat stress, the duration of heat stress (i.e., 1-8 consecutive days of THI > 65) also has a negative impact on milk production, decreasing fat, and protein content (Ouellet et al. 2019).
Dairy production in Canada takes place in varied climates, including cooler humid climates in coastal regions, warmer humid conditions in Southern Ontario, and arid climate in the prairies.The nearly 10 000 Canadian dairy farms (∼930 000 lactating cows) are distributed across every province with particularly dense concentrations in British Columbia (Fraser Valley), Central Alberta, and the corridor from Southern Ontario to Eastern Quebec along the St-Lawrence River (Fig. 1, VanderZaag et al. 2023a).Genetically, ∼93% of Canadian dairy cows are Holstein and most farms keep lactating cows in barns year-round (Bewley et al. 2017;AAFC 2021;CDIC 2021).Studies have shown that in-barn THI conditions tend to be higher than outdoor conditions (Schüller et al. 2013;Shock et al. 2016;Mylostyvyi et al. 2020;VanderZaag et al. 2023b).If available, THI conditions experienced in-barn by the cows should therefore be considered when quantifying the impact of heat stress.Global climate models (GCMs) have projected increased temperature in Canada in the 21st century (Li et al. 2018;Sobie et al. 2021) and temperatures are increasing at about double the rate of the global increase (Bush and Lemmen 2019;IPCC 2021).Considering all parts of the country, there was an overall increase in the number of days with average temperature >25 • C from 1948 to 2016 (Vincent et al. 2018).However, given the regional variation among the areas where dairy farming is concentrated, broad conclusions for Canada have limited use for informing dairy producers about changes in their specific regions.Agriculture-specific climate change analyses have primarily focused on crop production (e.g., Rötter et al. 2018;Qian et al. 2019).Of those that focused on livestock, Key and Sneeringer (2014) used modelling to project that US milk production will decrease by 0.60 to 1.35% by 2030 relative to the projected increasing trend (i.e., increase at a slower rate) as a direct result of climate change (assuming no mitigation or countermeasures taken).Canadian studies of climate change specific to dairy cattle are limited to one study in Quebec (Ouellet et al. 2020).They noted increased frequency and duration of high THI in Eastern and Southwestern Quebec in two future periods.Their results showed that compared to the reference period , the number of days when THI max > 68 will increase by an average of 18 and 39 days in the 2020-2049 horizon, and by 35 and 63 days for the 2050-2079 horizon under a moderate greenhouse gas (GHG) emissions pathway, for the Eastern and Southwestern regions, respectively.
The objective of this study is to quantify the historical and projected trends in the THI heat-stress indicator in key dairyproducing regions across Canada.Historical data from 1960 to 2019 and projections until 2089 are used for analysis.The future period considered is within the long-term management time-horizon of the sector considering dairy barn lifespans and the timespan for selective breeding.By focusing the analysis on dairy-producing regions, this analysis provides insight that is highly relevant to the Canadian dairy industry.

Historical Environment and Climate Change Canada weather dataset
Although historical weather observations have been made at hundreds of locations across Canada, relatively few stations collect long-term RH a data, which constrained the choice of stations available to perform our analyses.We selected 29 Environment and Climate Change Canada (ECCC) weather stations located in relevant dairy agricultural areas and with a long historical record of hourly relative humidity (RH a ), air temperature (T a ), and air pressure (P) data.Since RH a and T a are counter-cyclical, outdoor specific humidity (SH a ; g water kg −1 air) was also calculated and used as an indicator of humidity.This parameter (SH a ), the mixing ratio of water vapour and air, was calculated based on observed RH a , T a , and P from ECCC weather stations using the equations in Appendix A. The location of each of the stations is shown in Fig. 1 and additional details, including station ID and geographic coordinates, are provided in Table A1.When hourly data (T a and/or RH a ) were missing, matching entries from a nearby weather station (within 50 km) were used.This only affected 433 out of 627 132 daily records, with the mostimpacted stations being Pilot Mound (249 records), Swift Current (41 records), and Mont-Joli (37 records), while 17 stations did not have any missing data.THI a was calculated hourly and used to calculate daily mean THI a .If more than 7 h of data were missing in a day, then that day was considered incomplete.The only exception to this rule being the Pilot Mound station (station ID 5022125) at which, from 1960 to 1986, observations were made four times daily at 00 h, 06 h, 12 h, and 18 h.Overall, the historical weather dataset had 627 132 daily records over the period from 1 January 1960 to 31 December 2019.For assessment of historical changes, the means between 30-year periods of 1960-1989 and 1990-2019 were compared.Thirty-year periods were used because this is considered long enough to reflect long-term changes while being sufficient to eliminate year-to-year variation.This approach aligns with the World Meteorological Organization guidelines for calculation of climate normals (WMO 2017).

Modelled weather dataset
Climate scenario data (daily T a and RH a ) used in this study were obtained from phase 3b of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP3b; https://www.isimip.org/protocol/3/).The ISIMIP3b climate scenarios were derived from CMIP6 (Coupled Model Intercomparison Project Phase 6) output on 0.5 • grids using the bias adjustment and statistical downscaling method ISIMIP3BASD (Lange 2019), based on the observational dataset W5E5 that is version 2.0 of WFDE5 over land merged with ERA5 over the ocean (Cucchi et al. 2020).
The five primary GCMs included in ISIMIP3b were used in this study.They are as follows: (i) GFDL-ESM4 (Dunne et al. 2020), (ii) IPSL-CM6A-LR (Bonnet et al. 2021), (iii) MPI-ESM1-2-HR (Gutjahr et al. 2019), (iv) MRI-ESM2-0 (Yukimoto et al. 2019), and (v) UKESM1-0-LL (Sellar et al. 2019).ISIMIP3b currently includes climate scenarios under three Shared Socioe-conomic Pathways (SSPs; Riahi et al. 2017).Due to uncertainties associated with climate models, the results from all five models were combined as an ensemble mean and standard deviation (Qian et al. 2021).Each of the five models was run once under three SSPs, which represent "pathways" that the world could follow in the future: SSP1 (which has an increased radiative forcing of 2.6 W•m −2 ), SSP3 (7.0 W•m −2 ), and SSP5 (8.5 W•m −2 ), henceforth referred to simply as SSP1 or "low end," SSP3 or "medium-high end," and SSP5 or "high end."SSP1 refers to a world of sustainability-focused growth and equality, taking a green road with low consumption and population growth.SSP3, which takes a pathway called a rocky road with high challenges for mitigation and adaptation, a resurgence in nationalism, concerns about competitiveness and security, and regional conflicts push countries to increasingly focus on domestic or, at most, regional issues.In contrast, SSP5 depicts a world of fossil-fueled development with high challenges to mitigation and a focus on high economic growth and consumption.

Calculation of THI indicators and statistical analyses
Climate projections of daily mean T a and RH a were used to calculate daily THI a at each location using: where T a is outdoor air temperature ( • C) and RH a is outdoor relative humidity (%) measured at 2 m.The equation was developed by the National Research Council: a guide to environmental research on animals (NRC 1971).T a and RH a were either obtained from weather station data or outputs from the ISMIP3b climate scenarios.
In addition, daily THI conditions expected inside a naturally ventilated dairy barn (THI b ) were estimated using an equation developed by VanderZaag et al. (2023b): where THI b is the daily average in-barn THI and THI a is the daily average outdoor THI.The equation had an R 2 of 0.95 over a range of THI a from 55 to 78 (VanderZaag et al. 2023b).
Changes in the prevalence and extent of heat stress were assessed in three ways: (i) By calculating the number of days per year that exceeded a THI threshold of 65 or that fall within specific THI ranges (i.e., 60 < THI ≤ 65, 65 < THI ≤ 70, THI > 70, similar to Ouellet et al. 2020); (ii) by comparing the THI load degree days (cumulative THI exceedance over 65), which combines the magnitude and duration of elevated THI exposure; and (iii) by determining changes in the duration of heat-stress periods (HSP max ).The calculation of these indicators is described below.
The THI ranges were selected strategically to describe sub-threshold THI (60 < THI ≤ 65), as well as moderate (65 < THI ≤ 70), and high (THI > 70) threshold exceedance.As an indicator of cumulative exposure to elevated THI, THI load (in degree days) was determined based on Key and Sneeringer (2014).THI load was calculated as the sum of exceedance over the mean daily THI a threshold of 65.The cumulative sum of these values was used to determine annual THI load in degree days.Although THI itself is unitless, the THI load was assigned the units of degree days as an intuitive term for this accumulation over a threshold metric, estimated in a similar manner as growing degree days for crops.Maximum heatstress period duration (HSP max ) was calculated by counting the number of consecutive days where mean daily THI > 65, and then determining the maximum value (i.e., the longest heat-stress period) in each year at each location.
Certain calculations were restricted to the "warm period" of the year, which was defined as the 6-month period inclusive of 15 April-14 October.Likewise, for brevity, some analyses were focused primarily on the "medium-high end" SSP (SSP3--7.0W•m −2 ).Calculations, tabulation, and statistical analyses were done in JMP 16.2 (SAS Institute, Inc., Cary, NC, USA).Statistical differences were tested using Student's t tests at α = 0.05 and error bars represent the 95% confidence interval (95% C.I.) of the means.

Bias correction for individual stations
Historical weather data were used to further correct possible bias in the climate projections for each of the 29 station locations in Canada, while maintaining GCM-projected changes in future periods.Each of the 5 GCMs for each of the 29 locations were used to generate a historical dataset , which was then used to calculate mean THI a and THI a indicator values (>65 THI days, THI Load, and HSP max ) for this 30-year period.Indicators derived from each GCM were compared to the baseline observed ECCC mean values from each of the 29 locations for the same 30-year period.A single offset, derived from comparing the 30-year observed and modelled historical climate, was then added or subtracted to correct the difference between modelled and measured indicators for each site/GCM combination (Qian et al. 2022).The same offset was then applied to each indicator projected for future 30-year periods (2020-2049, 2040-2069, and 2060-2089).

Historical climate trends
Climate data (T a , RH a , THI a , and SH a ) were compared between two 30-year periods (1960-1989 and 1990-2019) to target the decades of the 1970s and 2000s, respectively.Focusing only on the "warm period" (15 April-14 October), we observed that mean T a has increased significantly in Canadian coastal regions (East and West) as well as in Eastern Canada (Table 1).On the other hand, "warm period" mean T a has not increased significantly in Central (MB) and Western Canada (SK, AB) (0.152 < p < 0.996).In other words, T a has increased significantly in Canada's more humid climatic zones.This trend was also observed when considering mean annual T a (as opposed to "warm period" only) between the two reference periods, except in this case Calgary, AB, Swift Current, SK, and Winnipeg, MB saw significant increases in T a , whereas Gander and Deer Lake (NL) did not see any significant changes (data not shown).
Many of the prairie locations saw significant increases in RH a between the two historical periods.In ecozones with already very high RH a (i.e., Atlantic Canada), changes in RH a were not observed between the two time periods, because RH is a relative measure that depends not only on how wet the air is, but also on T a .The temperature-independent measure of humidity, SH a , increased significantly at all but three locations, indicating that the air became wetter as T a increased (which maintained RH a over time), as expected.Only two locations in Central Alberta (Red Deer, Edmonton) had no significant changes in any parameter (T a , RH a , THI a , and SH a ).Significant changes in warm period THI a were observed in all locations except the prairies (AB, SK, and MB), which coincides with locations where there were significant changes in T a .This is consistent with the fact that although RH a has an influence on THI a , most THI formulas are more heavily weighted toward temperature.
The historical data show an increase in the number of days with THI a > 65 over time at several locations (Table 2).The change is greatest in the most humid locations (e.g., Halifax, NS), where the mean number of days when THI > 65 has increased from 22.9 day year −1 to 36.0 day year −1 between the 1970s and the 2000s (net increase of 13.1 days or 57%).In both periods, Windsor (ON) has the highest mean annual number of days of THI a > 65 (>75 days), and has increased significantly, but the increase is less sharp (net increase of 10.6 days or 13%).Calgary, Edmonton and Red Deer (AB) were the coolest locations, not exceeding 10 mean annual THI a > 65 days in either of the periods.Similar to T a , by comparing the mean annual number of days of THI a > 65 between the two periods (1960-1989 and 1990-2019), we see that means are significantly different in southern coastal regions (East and West) as well as Southern Ontario and Quebec, but not in the prairies (AB, SK, MB) (Table 2).At one prairie location--Estevan, SK--a significant decrease in heat-stress indicators was observed.Generally, significant differences in the THI a > 65 indicator were also observed in the THI a load indicator, except in the case of Sherbrooke and Mont-Joli, QC, where THI a load did not change significantly between the periods.
Fewer locations had significant changes in HSP max than the previously discussed indicators.The locations where significant mean differences were observed are in humid coastal regions (BC, NS) and warmer humid in-land locations (ON).The locations that saw significant increases were Vancouver, BC, Windsor, ON, Greenwood, NS, and Halifax, NS.This is consistent with the warming trends observed at these locations.The increase in Windsor, ON, was particularly striking as HSP max rose from 24.4 days year −1 (1960-1989) to 35.4 days year −1 (1990-2019), indicating the average year has over 1 month of consecutive days with daily mean THI a > 65.This finding is pertinent as there is a high concentration of dairy cows in Southern Ontario (Fig. 1); however, the more inland station in London, ON, did not show a significant HSP max increase (Table 2).

Climate change projections
We looked at projected data using three partially overlapping 30-year periods (2020-2049, 2040-2069, and 2060-2089) to compare the 2030s, 2050s, and 2070s, respectively.Under SSP1, there were no significant changes in the mean number of THI a > 65 days between the 2030s and the 2070s for any of the locations (data not shown).However, focusing on a "medium-high end" pathway (SSP3--7.0W•m −2 ), significant differences in the number of THI a > 65 days between the 2030s and the 2070s were observed at all locations (Table 3).
Table 1.Comparison of historical "warm weather" (15 April-14 October) climate data between two 30-year periods (1960-1989 and 1990-2019), representing the 1970s and 2000s, respectively.The projected changes in THI a > 65 days also coincided with a shift into higher THI a categories as can be seen when data are grouped by THI ranges.Under SSP5, the mean differences are even more pronounced (data not shown).
In a relatively cool climate-stable region such as Central Alberta (e.g., Edmonton), under SSP3, it is expected that there will be 6.8 days year −1 with THI a > 70 in the 2030s and 31.0 days year −1 in the 2070s.Under the same assumptions, in a more rapidly changing climate (e.g., Halifax) the number of >70 THI a days are predicted to be 19.4 days year −1 in the 2030s, rising to 54.0 days year −1 by the 2070s.
Under SSP3, mean THI a load projections are expected to increase significantly from the 2030s to the 2070s in some regions, and most regions under SSP5 (Fig. 2).No significant changes were observed between periods in any locations under SSP1.Although increases are expected in many locations, the magnitude of THI a load increases vary between the 2030s and 2070s.This variability under SSP3 is represented graph-ically in Fig. 3 and identifies Ontario, Canada's second most populous dairy region, as a hotspot for THI a load increases.In addition to Ontario, substantial increases in mean THI a load between the 2030s and 2070s are also projected for sites in Manitoba and Saskatchewan.London (ON) and Windsor (ON) are expected to have an increase of 373 and 439 THI degree days between these two periods (2030s and 2070s), respectively, compared to an increase of <200 THI degree days in parts of Alberta (Calgary) and Newfoundland (Deer Lake and Gander).
Another important factor when considering the effect of increased THI is the expected duration of sustained elevated THI a periods.No significant increase in HSP max were observed at any of the locations under SSP1, but under SSP3 and SSP5, significant increases in HSP max are observed in the 2030s, 2050s, and 2070s.Focusing on SSP3, predicted trends in mean HSP max show similarly sloped linear increases over time in most locations (Fig. 4).The only exception being cooler regions that currently have the lowest THI a load and >65 THI a Table 2. Comparison of outdoor temperature-humidity index (THI a ) load (degree days year −1 ), >65 THI a (days year −1 ), and maximum heat-stress period (HSP max )(days) between two 30-year periods (1960-1989 and 1990-2019)  days year −1 (Deer Lake, NL, Gander, NL, Calgary, AB, Red Deer, AB, and Edmonton, AB).In all other sites, HSP max is predicted to increase by 5.6-8.7 days per decade from the 2030s to the 2070s under SSP3 and by 8.6-11.4days per decade under SSP5.Trends in the mean HSP max (days year −1 ) are presented in Fig. 4.

Consideration of In-Barn conditions
Compounding the effect of the increasing THI a trends is the observation that in barn THI (THI b ) tends to be higher than THI a (VanderZaag et al. 2023b).Using equation 2, we can expect THI b to be ∼5%-10% greater than THI a .The effect of this increase on THI load is illustrated on a selection of locations in Table 4.For example, in Windsor, ON, an additional 164-171 degree days (depending on period and SSP) are expected in-barn compared to outdoors.In the case of a cooler location (Red Deer, AB), the in-barn degree days are higher by 57.1-122 degree days (depending on period and SSP).In this case, the magnitude of the difference is more heavily dependent on the future pathway.Under SSP1, the THI load is greater inbarn compared to outdoors by 57.1 degree days (2020-2049) and 73.2 (2060-2089), but under SSP5 the difference is similar for 2020-2049 (62.8 degree days) but climbs rapidly to 122.1 degree days by 2060-2089.
The increases outlined above assume constant barn management over time, which could be considered a worst-case scenario, as new innovations and practices are likely to be developed and implemented that could reduce THI b relative to the baseline THI a .The same analyses were done using equations for mechanically ventilated barns (not shown), which also resulted in more THI load in-barn than outdoors.

Historical changes in climate
Previous analyses of historical temperature data have shown that Canada's annual mean temperature has increased by 1.7 • C from 1948 to 2012 (Vincent et al. 2015).However, Table 3. Mean outdoor temperature-humidity index (THI a ) > 65 days (days year −1 ), as well as breakdown of >60 THI a days (%) into three intensity categories (60 < THI a ≤ 65; 65 < THI a ≤ 70; THI a > 70), comparing 2020-2049 (2030s)  more historical warming has occurred in winter than in summer (Zhang et al. 2011;Vincent et al. 2018).While winter warming has many impacts on dairy farming--such as growing season length, overwintering of pests--it does not increase the prevalence of heat stress.Spatially, previous analyses found the highest increase in the mean annual temperature in northern Canada (primarily non-agricultural), yet summertime warming was greatest in southern Canada (Zhang et al. 2011).
Past analyses of global surface humidity data from 1973 to 2003 have found an "overall moistening of the globe" as indicated by SH a (Hadley Centre and Climate Research Unit humidity dataset; Willet et al. 2008;Willett and Sherwood 2012).Our analyses are consistent with this finding and show that the historical increases in SH a occurred in many provinces.This increasing trend in humidity means that increases in THI can be driven both by rising temperature (in some regions) and by rising humidity (in all provinces).Rising humidity also reduces the efficacy of evaporative cooling, a potential adaptation strategy (Fournel et al. 2017).
Our analysis of historical THI showed an increase in warmseason THI from the 1970s to 2000s in most locations in BC and Eastern Canada, but not in the prairies.Stable temperature with increased humidity in the prairies is consistent with a widespread change in land-use, where crop acreage increased and summerfallow decreased (Gameda et al. 2007).This spatial distribution of warming is consistent with changes in temperature reported by Vincent et al. (2018), who found little historical summertime warming (and some cooling) in parts of the prairies (notably SK and central AB), as indicated by the number of cooling degree days and the number of hot days.They and others (Shah et al. 2022) have also noted increasing warm nights in the ON-QC corridor, which agrees with our findings related to increased THI in Southern and Eastern Ontario.The increase in warm nights reduces the opportunity for dairy cows to have Fig. 2. Mean (95% C.I.) THI a degree days per year by 30-year period: 2020-2049 (2030s), 2040-2069 (2050s), and 2060-2089 (2070s).Bars represent the annual mean of five global climate models (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL) for each combination of period and socioeconomic pathway (SSP1--2.6W•m −2 , SSP3--7.0 W•m −2 , and SSP5--8.5 W•m −2 ).Error bars are constructed using a 95% confidence interval for each 30-year period.Under SSP1, no significant differences between periods were detected (Student's t test at α = 0.05).Under SSP3, significant differences were detected for all locations except Red Deer (AB), Calgary (AB), Deer Lake (NL), Edmonton (AB), Gander (NL), Quebec (QC), and Mont-Joli (QC), which were marginally significant (0.05 < p < 0.10).Under SSP5 significant differences were detected for all locations except for Deer Lake (NL), which was marginally significant (0.05 < p < 0.10).
reprieve from heat stress.This trend suggests that understanding the role of night-time cooling (e.g., THI min ) on dairy cattle productivity may be important as this metric may change more than THI max in the future.

Projected future changes in climate
We observed some variation in projected THI among the 5 GCMs, but the main source of variation in future projections is the SSP.This is logical since each SSP represents changes in GHG emissions, and thus the projected climate change.All three SSPs resulted in increasing frequency and extent of THI exceedance, but future change was not statistically significant at many locations under SSP1--the most optimistic pathway.The more moderate pathway (SSP3) resulted in substantial changes in future conditions.For example, in the 2070s, Abbotsford, Winnipeg, Sherbrooke, and Halifax are projected to experience a greater number of days with THI >65 than the warmest location in the country had in the 1970s (Windsor).Changes under SSP5 were even more striking.Although dairy farmers can take action to reduce their GHG emissions, no individual sector can control the global SSP that will unfold in the future.Thus, it seems prudent to prepare for a substantial increase in the frequency and extent of THI exceedance.
The high proportion of days that fall in the 60 < THI a ≤ 65 category highlights the importance of research to determine the appropriate THI threshold at which milk production is negatively affected in Canadian dairy farms (Rockett et al. 2023).It also highlights the need to consider regions where THI a is consistently stable just below heat-stress thresholds.In the future, farms that have not previously considered heat stress as part of barn design and management should put more emphasis on this aspect.Fig. 3. Map showing the mean annual THI a load (degree days) for observed historic periods  as well as projected future periods (2020-2049, 2040-2069, and 2060-2089).THI load for future periods was generated from an ensemble of five global climate models (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL) assuming a "mediumhigh end" future pathway (SSP3--7.0W•m −2 ).
According to the projections in our study, the frequency of THI exceedance will rise significantly from the 2030s to 2070s in all regions, including all four locations in Quebec.This is consistent with those observed by Ouellet et al. (2020) for two regions in Quebec.Their SIM6 model (based on 8.5 W•m −2 scenario) predicted 176 days of THI max > 65 for the period of 2050-2079 for Southwestern Quebec compared to our estimate of 122 days of THI avg > 65 in St-Hubert for the period 2060-2089 under SSP5 (8.5 W•m −2 ).The differences between these values are reasonable since the daily THI max will exceed 65 more often than THI avg .Likewise, for the duration of heat-stress periods, Ouellet et al. (2020) found an increase in the proportion of days in the 7-8 days consecutive heatstress category (their highest category).Our analysis showed the same trend, namely HSP max , increased from 16 in the 1970s to 28 in the 2030s, and 57 in the 2070s under SSP3 for St-Hubert, QC (31 and 71, respectively, under SSP5).The potential impact of both increased THI, and duration of THI periods have been shown to negatively impact the yield of energy-corrected milk, as well as the fat and protein content (Ouellet et al. 2019).

Impact on Canada's dairy cattle
Approximately 33% of Canada's dairy herd is located in Ontario (VanderZaag et al. 2023a), with the majority in the southern parts of the province.This region has the highest frequency and duration of heat-stress indicators in both historical and future analyses.Significant increases in THI indicators are projected for this region.For example, under SSP5, southwestern Ontario is projected to experience a THI load approaching levels reported for the southern United States (Key and Sneeringer 2014), with sustained heat-stress periods of more than 1 month (e.g., HSP max of 48.6 for London, ON in the 2050s under SSP3).This projected warming may pose new challenges to the provincial dairy sector's ability to increase production year-over-year to match population growth, and to maintain profit margins and production throughout the warm season.On the other hand, while Ontario clearly faces the greatest heat stress, the current heat-stress challenges of Southern Ontario could be used as an early adaptation strategy incubator for regions that are likely to track similar climate patterns in the future.For example, Eastern Ontario will approach the heat levels currently experienced in Southwestern Ontario by the 2030s and exceed these levels by the 2070s under SSP3.To give another example, parts of Atlantic Canada are expected to reach the heat-stress levels currently experienced in Eastern Ontario by the 2030s and those currently experienced in Southwestern Ontario by the 2070s.
Alberta and Newfoundland are the two provinces with the least THI exceedance, but account for only a modest 8.5% and 0.6% of the national dairy herd, respectively (CDIC 2023a).Although historical data show no increase in THI at stations in Alberta, future projections suggest that a significant increase in both provinces may occur under SSP3 and SSP5.As a result, Red Deer, AB in the 2070s (under SSP3 or SSP5) may be like Ottawa, ON was in the 2000s.While the resulting frequency and extent of THI exceedance are still low compared to other locations, this may necessitate changes in the way barn ventilation and fans are designed and/or managed in these areas to give more consideration for cooling the barn in the warm season in addition to retaining heat in the winter.For example, a previous study included a naturally ventilated dairy barn in central Alberta, which did not have a single fan for air circulation nor for ventilation--a design that would likely need to change in the future (VanderZaag et al. 2023b).
The rest of the dairy-producing regions fall between these extremes.Quebec, with ∼37% of the national dairy herd (CDIC 2023a) has had modest historical increase in THI exceedance and projected to have significant increases in all THI indicators in most future scenarios.Current and pro-jected increases in THI will be highest in the highly productive agricultural region in southwestern Quebec, represented by St-Hubert.On the west coast, British Columbia has approximately 8% of the national herd and has had a significant historical increase in THI, which is projected to continue at a slower pace, resulting in a moderate increase under most future scenarios.The east coast maritime provinces (New Brunswick, Nova Scotia, PEI) account for around 5% of the national herd and have a similar historical trend as BC, but with more warming projected for the future.Saskatchewan and Manitoba combine for around 7% of the national herd.These provinces have not had a historical increase in THI, but this is projected to change in the future with significant increases in all THI-related indicators.
Taken all together, 84% of the national dairy herd has experienced significant historical increases in THI, and over 90% of the national herd is projected to experience a high level of increase in the future.The overlap between Canada's most intensive dairy regions and regions expected to experience Table 4. Comparison of temperature-humidity index (THI) load (degree days) between outdoor (THI a ) and in-barn (THI b ) in future periods (2020-2049, 2040-2069, and 2060-2089) under various future pathways (SSP1, SSP3, and SSP5) for select locations.
Outdoor, THI a load (degree days) In-Barn, THI b load (degree days) While the projected changes are substantial, the conditions are not extreme by global standards.For example, Ouellet et al. (2021) observed that THI was typically >70 from May to October in Florida.Flamenbaum and Lavon (2020) noted that the percentage of the year with THI > 68 was nearly 50% in Florida, and >40% in Israel.Areas of the Southern United States see sustained periods of THI > 70 (Bohmanova et al. 2007) and annual THI loads >1700 (Key and Sneeringer 2014).In Italy, Bertocchi et al. (2014) reported the average monthly THI max was above 79.6 for 3 months during 2003.In Mexico, Rodriguez-Venegas et al. ( 2022) reported THI > 80 for 31.2days year −1 .On the other hand, thermotolerance is affected by climate and by heat exposure i.e., a cow in Florida will tolerate higher heat load than a cow raised in Canada (Bohmanova et al. 2005).Moreover, older high producing cows will be impacted by lower THI (Kamal et al. 1970).

The need for adaptation
Canadian dairy farms strive to maintain consistently high milk production year-round, which is another reason farmers may want to avoid production declines caused by heat stress (Campos et al. 2022a).Canadian dairy cows have a strong track record of continuous improvements increasing the genetic potential for high production, and this has resulted in long-term increases in per-cow milk production and a decrease in the overall number of dairy cows in Canada (Jelenski et al. 2015;VanderZaag et al. 2023a; CDIC 2023b).
However, increasing milk production is also associated with increased sensitivity to heat stress (Bohmanova et al. 2005) and a negative genetic correlation between heat tolerance and production was noted previously (Ravagnolo and Misztal 2000;Campos et al. 2022b).This implies continuous selection for production may further decrease tolerance to heat stress.In addition, the Canadian dairy herd is composed primarily of Holstein cows, which are more sensitive to climaterelated stress than other breeds (Sharma et al. 1983).Therefore, there may be a long-term benefit to including heat-stress tolerance as a breeding selection criterion, alongside production and other traits (Carabano 2016;Garner et al. 2016;Caranano et al. 2019).To that end, the use of multi-omics tools to identifying selection criteria and markers for heat tolerance would be a valuable area of future research in Canada.
Strategies related to management and barn design are also important.Some dairy farmers in Canada have already begun adapting by replacing old barn designs with new ones that have greater ventilation, lower stocking density, and high-speed fans over the stalls.An area for further work is to decrease the gap between in-barn THI and outdoor THI, potentially through practices like increasing ventilation at night when outdoor THI tends to be lower than indoor, especially when there is less wind to drive ventilation naturally (VanderZaag et al. 2023b).It is notable that Shah et al. (2022) found a decrease in the wind speed in Ontario consistent with a broader trend of atmospheric stilling.This may mean that hybrid barn designs combining mechanical and natural ventilation may become more appropriate in Ontario and some other Canadian locations.It will be important to adapt activation thresholds for ventilation systems to better represent thermotolerance of Canadian Holstein cows.While the sector may be able to maintain production using adaptation strategies, these additions will also increase the cost of production (Brouk et al. 2001;Key and Sneeringer 2014).
In the broader context of dairy production, the THI and warm-season heat stress is only one aspect of weather and climate change.Changes in water availability, water demand, crop production, pest overwintering, and extreme weather events are all critical to understand and develop adaptation strategies.For example, Qian et al. (2022) point to an expansion of range favourable to soybean, but at the same time a projected precipitation deficit in semi-arid regions, such as the prairies, is likely to make crop production in these regions water-constrained without expansion of irrigation.Likewise, THI is only one heat-stress indicator and it does not include the role of wind (important for barn ventilation and heat dissipation) and solar exposure (important for pastured cows).There are many areas in need of further research to better understand the intersection of climate change projections and dairy production to inform proactive management and policies.

Conclusions
Historical climate trends show an increase in temperature, humidity (relative, and specific), and increased THI exceedance in most locations in BC and Eastern Canada, affecting about 84% of the current national dairy herd.Models of future conditions indicate that 90% of the national herd is projected to experience large increases in the frequency, severity, and duration of THI exceedance.The largest source of uncertainty is the future emission scenario--which farmers cannot control--and significant warming is projected in all except the most optimistic future scenarios.In a moderate future pathway, the frequency of THI exceedance in many locations across the country in the 2070s is projected to be greater than what the warmest location in the country was in the 1970s.Thus, adaptation will be important for most dairy farms.These findings highlight the need for Canadian dairy farmers to consider appropriate strategies to adapt to elevated THI.Given the long-time horizon for climate change, both individual (e.g., barn design, management) and collective (e.g., breeding) strategies will be beneficial.

Fig. 1 .
Fig. 1.Map of the 2016 dairy cow population in Census Agricultural Regions, and the 29 weather station locations (circles) used for analysis of historical and future climate.The map was generated in ArcGIS using the 2016 Census agricultural regions boundary files from Statistics Canada (2019), the Lambert Conformal Conic projection, and the GCS North American 1983 coordinate system.Dairy cow population is based on the 2016 Census of Agriculture by Statistics Canada.
representing the 1970s and 2000s.Mean > 65 THI a (days year −1 ) M e a n T H I a Load (degree days year −1 ) M e a n H S P max (days) Note: Calculated from historical climate data measured at ECCC climate stations.Significant differences were analysed using student's t test at α = 0.05.
That being said, most locations are still in the middle THI range of 60-70 most of the time and have annual THI loads below 500 degree days.Only locations in Southern Ontario are projected to exceed annual THI loads of 1000 degree days in the future.

Table A1 .
Table of Environment and Climate Change Canada weather stations used in the study.The period of record was between 1960 and 2020.* Climate ID provided is the most recent one, as many ID numbers changed between 1960 and 2020. Note: