Phenotypic analysis of heat stress in Holsteins using test-day production records and NASA POWER meteorological data

Weather station data and test-day production records can be combined to quantify the effects of heat stress on production traits in dairy cattle. However, meteorological data sets that are retrieved from ground-based weather stations can be limited by spatial and temporal data gaps. The National Aeronautics and Space Administration Prediction of Worldwide Energy Resources (NASA POWER) database provides meteorological data over regions where surface measurements are sparse or nonexistent. The first aim of this study was to determine whether NASA POWER data are a viable alternative resource of weather data for studying heat stress in Canadian Holsteins. The results showed that average, minima, and maxima ambient temperature and dewpoint temperature as well as 4 different types of temperature-humidity index (THI) values from NASA POWER were highly correlated to the corresponding values from weather stations (regression R 2 > 0.80). However, the NASA POWER values for the daily average, minima, and maxima wind speed and relative humidity were poorly correlated to the corresponding weather station values (regression R 2 = 0.10 to 0.49). The second aim of this study was to quantify the influence of heat stress on Canadian dairy cattle. This was achieved by determining the THI values at which milk, protein, and fat yield started to decline due to heat stress as well as the rates of decline in these traits after the respective thresholds, using segmented polynomial regression models. This was completed for both primiparous and multiparous cows from 5 regions in Canada (Ontario, Quebec, British Columbia, the Prai-ries, and the Atlantic Maritime). The results showed that all production traits were negatively affected by heat stress and that the patterns of responses for milk, fat, and protein yields to increasing THI differed from each other. We found 3 THI thresholds for milk yield, 1 for fat yield, and 2 for protein yield. All thresholds marked a change in rate of decrease in production yield per unit THI, except for the first milk yield threshold, which marked a greater rate of increase. The first thresholds for milk yield ranged between 47 and 50, the second thresholds ranged between 61 and 69, and the third thresholds ranged between 72 and 76 THI units. The single THI threshold for fat yield ranged between 48 and 55 THI units. Finally, the first and second thresholds ranged between 58 and 62 THI units and 72 and 73 THI units for protein yield, respectively.


INTRODUCTION
Homeotherms require effective thermoregulatory mechanisms to maintain a stable body temperature and subsequently preserve normal physiological form and function (Kadzere et al., 2002).However, most mammals can only adequately retain this state within a small range of environmental parameters.This range is defined by their tolerance limits (Miller and Stillman, 2012).Heat stress occurs when an animal is not able to cope well with environmental heat loads outside their tolerance limits (Das et al., 2016).In dairy cattle herds, this can have several negative health, economic, and welfare implications (Kadzere et al., 2002).
Furthermore, dairy cattle are extremely prone to heat stress because of their high metabolic heat production, which is associated with their rumen fermentation and high milk production (Collier and Gebremedhin, 2015).Improvements in animal nutrition and biotechnology as well as genetic selection have greatly increased the amount of milk produced per cow.Because increase in milk production is positively correlated with increase in metabolic heat production, these improvements in milk production capacity have most likely affected the thermoregulatory profile of dairy cattle and potentially reduced their ability to adequately cope with heat stress (Kadzere et al., 2002).Additionally, if anthropogenic activities continue at the current rate, it is likely that the global mean surface temperature will increase by another 0.5°C by 2040, reaching a total of 1.5°C above pre-industrial levels (Hoegh-Guldberg et al., 2018).If greenhouse gas emissions increase at a greater rate than current emissions, it is likely that global mean surface temperature could increase to at least 2°C above preindustrial temperatures.
It is also important to note that these changes in global mean surface temperature will not be uniform across the globe (Hoegh-Guldberg et al., 2018).Surface warming in North America and Europe will likely be amplified, resulting in significantly higher average temperatures than the global mean temperature.The mean ambient temperature may be as high as 2.0°C above pre-industrial temperatures for North America by 2050 if the current rate of greenhouse emissions continues (Romero-Lankao et al., 2014).Evidence also suggests that the number of warm days and nights have increased and that the length and frequency of warm spells over 20% of the globe, including Canada, have significantly increased from 1951 to 2003 (Alexander et al., 2006).Heat stress in dairy cattle should be investigated even in temperate regions such as Canada due to these changes in climate patterns and the potential reduction in thermal tolerance limits of dairy cattle.
To better understand the thermal tolerance limits of dairy cattle, numerous studies have combined weather station data and test-day production records to define the temperature-humidity index (THI) thresholds at which production and fertility traits begin to be negatively affected (Ravagnolo et al., 2000;Vitali et al., 2009;Zimbelman et al., 2009).However, a common challenge for these studies is the lack of consistent and reliable weather data, because weather stations can have numerous temporal data gaps and a low spatial resolution.Automatic weather stations (AWS) collect and transmit weather observations without the need of volunteers.Relative to volunteer-operated or manned stations, AWS improve the reliability of weather station observations (WMO, 2008).It is more difficult for manned weather stations to maintain routine record keeping, and they can also have significant site-to-site variability due to the use of human observers (Holder et al., 2006).This can lead to inaccurate or inconsistent data collection and data set gaps.White et al. (2008) found a high proportion of missing data reported by manned weather stations in the United States, due to the diminished ability of volunteer staff to report values on weekends and holidays.Additionally, it may not be possible for manned weather stations to collect weather data throughout the whole day, potentially skewing daily averages calculated from hourly observations (Colston et al., 2018).Wu et al. (2005) showed significant discrepancies in daily observations between 2 weather stations that were located 10 km apart.They speculated this may have been due to observation and sensor error, variation between times of data collection, and differing microclimates.Due to significant inconsistencies between different stations, caution should be taken when combining data sets, especially between AWS and manned stations.Environment and Climate Change Canada has a network of 585 fully automated stations, 465 volunteer climate stations, 304 aviation monitoring stations, and 381 other stations.The latter 2 groups included a mix of both AWS and manned stations in 2016 (Mekis et al., 2018).Overall, the large number of manned stations and the lack of consistency in the type of stations used in the Canadian weather station network could make it difficult to generate a large data set of reliable meteorological data.
The absence of weather stations in rural areas is another common challenge when studying heat stress in dairy cattle.Stations are typically located near large population centers and airports (Colston et al., 2018).Therefore, accurate weather data sets may not be available for remotely located herds.Overall, most weather station networks will likely have low spatial coverage over rural areas, consequently diminishing the number of herds that can be used in a study.It is also challenging to determine the maximum acceptable distance between a herd and a weather station that should be used.
A possible alternative weather data resource is the National Aeronautics and Space Administration Prediction of Worldwide Energy Resources (NASA POWER).NASA POWER provides solar and meteorological data sets for any location from 1981 to near-real time (Stackhouse, 2022).The overall aim of this study was to determine whether NASA POWER data represent a viable alternative resource of weather data for dairy cattle heat tolerance studies, as well as to quantify the effects of heat stress on Canadian dairy cattle.Therefore, the 2 specific objectives for this study were (1) to compare weather parameters and THI values collected from weather stations against NASA POWER estimates and (2) to estimate the THI thresholds and the subsequent production decays due to heat stress in primiparous and multiparous dairy cattle located in 5 regions in Canada.

MATERIALS AND METHODS
Because existing data sets were used in this research, approval by an Institutional Animal Care and Use Committee was not necessary.

Weather Data Resource Comparison
The Environment and Climate Change Canada weather station database (https: / / climate .weather.gc.ca/ ) provided weather data collected between 2009 and 2019 from 1,272 weather stations located across the Canadian provinces.This weather data included the hourly averages, minima, and maxima for relative humidity (RH, in %), dewpoint temperature (DP, in °C), wind speed (WS, in km/h), and ambient temperature (AT, in °C).Observations for WS were converted from kilometers per hour to meters per second.These hourly variables were then converted into daily averages, minima, and maxima.Stations that reported less than 85% of the number of hourly observations that could possibly be recorded within each day or year were removed from the analysis.Stations were also required to have at least 2 observations within each of the 4 time intervals for at least 85% of the days in the year, to account for diurnal changes in weather.These 4 daily time intervals were defined as 0000 to 0600 h, 0600 to 1200 h, 1200 to 1800 h, and 1800 to 2400 h.Therefore, each yearly weather data set needed to have a minimum of 310 d with at least 21 observations spread across the 4 different time periods throughout the whole day.Weather stations also needed to meet these requirements for at least 5 consecutive years between 2009 and 2019.
Because NASA POWER data are averaged over 0.5-degree latitude by 0.625-degree longitude geographical grids, the weather stations were subset into grid groups.The station closest to the centroid of each grid box was selected for the following analysis.In total, data from 575 weather stations located within 536 grid boxes were used in this comparison study.The locations of these weather stations are shown in Figure 1.The locations of the weather stations are spread across Canada, with a larger concentration of stations in the southern part of the provinces, where the majority of dairy herds in Canada are typically located.The daily averages, minima, and maxima for RH (in %), DP (in °C), WS (in m/s), and AT (in °C) from 2009 to 2019 were also downloaded from NASA POWER (https: / / power .larc.nasa.gov/ ) for each weather station location.
Four different types of THI values were calculated using these daily weather parameters.The 2 equations used to calculate these THI values were as follows: These THI equations are both commonly used in dairy cattle heat stress studies (Ravagnolo et al., 2000;Li et al., 2009).The 4 different types of THI values were the average THI using the first equation, the average THI using the second equation, the maximum THI using the first equation, and the maximum THI using the second equation.Each of these THI values were calculated twice for each weather station location, once using the weather station data and a second time using NASA POWER data.The average THI values using the first equation and the average THI values using the second equation were calculated using daily average AT and daily average RH or daily average DP.The maximum THI values using the first equation and the maximum THI values using the second equation were calculated using the daily maximum AT and daily minimum RH or daily minimum DP.The maximum THI values were calculated because Ravagnolo et al. (2000) found that different combinations of parameters had various correlations with production traits.They determined that the combination of maximum daily air temperature and minimum daily RH in the THI equation best quantified heat stress in dairy cattle.Linear regression was used to compare the daily NASA POWER parameters against the corresponding variables from the selected weather stations.The parameters compared included the 4 THI values and the average, minima, and maxima AT, RH, WS, and DP.

Phenotypic Analyses of Heat Stress
The initial data set for this analysis had 21.7 million test-day records of milk yield (kg), fat yield (g), and protein yield (g) measured between 2010 and 2019 from 1.2 million animals in approximately 8,500 herds located across Canada.These production data were provided by the Canadian Dairy Network (a member of Lactanet Canada, Guelph, Canada) and were grouped into 5 regional subsets.The regions were classified as Ontario, Quebec, British Columbia, the Prairies, and the Atlantic Maritime.The Prairies region included the provinces Saskatchewan, Manitoba, and Alberta.The Atlantic Maritime region included the provinces Newfoundland and Labrador, Nova Scotia, Prince Edward Island, and New Brunswick.Every herd was required to have known geographical coordinates and records from at least 5 different animals per year for at least 5 yr.
Each record was required to be from parity 1, 2, or 3, collected between 5 and 305 DIM, and having no missing trait measurements.Finally, animals were required to have at least 5 test-day records for each parity and were removed from the analysis if they had records from multiple herds.The addresses for each dairy herd were converted to latitude and longitude coordinates using Google Maps Geocoding (https: / / cloud .google.com/blog/ products/ maps -platform/ address -geocoding -in -google -maps -apis).The number of records, animals, and herds, as well as the averages and standard deviations for milk, fat, and protein yields, are listed in Table 1.
The production records were then adjusted for individual cow genetic effect and known environmental sources of variation.This was achieved by fitting a linear model to the data using the ASReml software (Gilmour et al., 2015).The equation for this linear model was as follows: where y ijklm is the mth test-day record, of either milk, fat, or protein yield, of the lth cow, μ is the overall mean, AD i is the fixed effect of the ith combination of age at calving and DIM, M j is the fixed effect of milking frequency (j = 1 to 3 classes), HY k is the fixed effect of the kth combination of the herd and year, c l is the random cow additive genetic effect, and e ijklm is the residual error term.For age at calving, 8 classes (<24, 25, 26, 27, 28, 29, 30, and >31 mo), 6 classes (<36, 37, 38, 39-40, 41-42, >43 mo), and 5 classes (<50, 51-52, 53-54, 55, >56 mo) were defined for the first, second, and third parities, respectively.For DIM 11 classes were defined (5-29, 30-60, 61-90, 91-120, 121-150, 151-180, 181-210, 211-240, 241-270, 271-300, and 301-305 d).The resultant residual term captures the variation in the production records that is not explained by the effects in the model but is caused by unknown fac- tors such as THI.In other words, the residuals are the phenotypic records that have been adjusted for known fixed and random effects and will henceforth be referred to as the adjusted phenotypes.
The daily minimums for RH (%) and daily maximum AT (°C) for each test-day record from 2010 to 2019 were retrieved from NASA POWER based on each herd's latitude and longitude.A THI value was calculated for the test-day and the 2 d before each test-day.Those 3 THI values were then averaged into a single value and assigned to a 2-d time lag variable labeled as THI 2d .This value was used to account for the possible delayed effect of heat stress caused by the thermal environment of the days before the test-day (West, 2003).The equation used to calculate THI values was the same as THI1: The adjusted phenotypes were averaged across herds within each THI unit starting at THI 30 for each region.The production averages were then plotted against THI to visually determine the average shape of the relationship between each production trait and an increasing environmental heat load.The threshold at which the production traits began to increase or decrease at a different rate were initially estimated by visually assessing each graph.In total, 3 breakpoints were found for milk yield, 1 for fat yield, and 2 for protein yield.These estimated breakpoints were then used to define several segmented polynomial regression models, which were fitted to each data set.As an example, a simple segmented polynomial model with 2 linear functions and 1 breakpoint could be described as follows (Toms and Lesperance, 2003): where y i is the average adjusted production yield for the ith THI 2d value, x i is the ith THI 2d value, B o is the intercept, x t is the THI 2d value that was selected as the estimated threshold, and B 1 and B 2 are the regression coefficients that describe the slope of the linear functions, with B 2 being the difference between the 2 slopes.The thresholds for each trait were adjusted systematically, and the resulting change in the coefficient of determination (R 2 ) value was monitored.The segmented polynomial model with the highest R 2 value was considered the best fit for that data set.

Weather Data Resource Comparison
Most studies on heat tolerance in dairy cattle require accurate and comprehensive AT and RH data sets collected over numerous years for multiple rural locations.Typically, such meteorological data are retrieved from a network of ground weather stations.Stations that are maintained properly can provide accurate measurements of ground conditions, but they can also have several limitations.These include sensor and human error, expensive routine maintenance, temporal and spatial data gaps, instrumental calibration drift, biased measurements due to the station location, and inconsistencies between stations due to different collection methods, recording, and instruments (Hubbard and Hollinger, 2005;Mahmood et al., 2006;Stackhouse, 2022).
Overall, this study found that several weather stations from the Canadian network had yearly and daily data gaps.For instance, only 575 stations out of the 1,272 Canadian weather stations with hourly data could be used in this analysis.Most of the stations were removed either due to a lack of data for any of the weather parameter throughout the year, because the observations collected on a specific day did not adequately represent the 4 time periods of the day, or because the station failed to report any weather parameter consistently for at least 5 years.This highlights how restrictive weather station data can be when studying heat stress on rural dairy farms.
A possible alternative to weather station data is NASA POWER, which provides long-term estimates of meteorological variables for any location in the world.The Modern Era Retrospective-Analysis for Research and Applications (MERRA-2) assimilation model provides meteorological data from 1981 to near-real time for the NASA POWER database (Stackhouse, 2022).The MERRA-2 model assimilates and optimizes observational data to estimate atmospheric variables.The observational data are collected from several sources such as remotely sensed information from satellites, ocean and land surface measurements, rawinsondes, aircraft and ship reports, and space-borne radar systems.The meteorological parameters reported by NASA POWER are obtained directly from MER-RA-2 or calculated using MERRA-2 products.The meteorological data are averaged over a geographical grid with a spatial resolution of 0.5° latitude by 0.625° longitude (Stackhouse, 2022).NASA POWER could provide meteorological data over regions where surface measurements are sparse, and could replace weather station data sets with substantial data gaps.However, the accuracy of NASA POWER meteorological estimates in Canada should first be assessed.The results from the weather data resource comparison analysis are shown in Table 2.
Overall, several weather parameter values from NASA POWER were highly correlated to the corresponding values collected from weather stations (regression R 2 > 0.85).These parameters included daily average, minima, and maxima AT and DP.However, the NASA POWER values for daily average, minima, and maxima WS and RH were poorly correlated to the corresponding weather station values (regression R 2 = 0.10 to 0.49).Other studies examining the accuracy of NASA POWER data have found similar results (White et al., 2008;Aboelkhair et al., 2019;Stackhouse, 2022).White et al. (2008) determined that NASA POWER temperature data showed good agreement with manned weather stations in the United States.They found that daily maximum AT and minimum AT between the 2 weather data resources were well correlated and reported a R 2 value of 0.88 for both variables.NASA POWER also conducted their own worldwide validation study by comparing surface station observations with the corresponding MERRA-2 estimates for every third year between 1981 to 2014.The R 2 values from comparing average AT, minimum AT, maximum AT, RH, DP, and average WS were 0.96, 0.93, 0.95, 0.61, 0.95, and 0.55, respectively (Stackhouse, 2022).
Similarly, Aboelkhair et al. ( 2019) found that NASA POWER estimates of RH did not have a high correlation with the weather station observations in Egypt (R 2 = 0.38).They found that weather stations located on the coast had a lower R 2 value than stations located more inland.White et al. (2008) also discussed adjusting for seasonal and regional variation.They speculated that proximity to water may bias NASA POWER temperature estimates.The estimates for grid cells that were located near an ocean or large lake had higher minimum AT values and lower maximum AT values.They also found that November to March temperatures had a greater error than other months.This may be a result of the greater variability in temperature during the winter months compared with the summer temperatures (White et al., 2008).Monteiro et al. (2018) also found that WS and RH estimation by NASA POWER requires further improvement.Nevertheless, the 4 types of THI values calculated from NASA POWER data were highly correlated to the corresponding weather station values (regression R 2 > 0.82).This was despite the low correlation found for the RH comparison.A possible explanation for this result is the relatively small weighting on RH in the THI equation that was used for this analysis.Therefore, changes in RH may not contribute to the final THI value to the same extent as would changes in AT.Some evidence suggests that this THI equation can adequately explain changes in production due to heat stress in less humid climates, but is less applicable to areas with greater humidity (Bohmanova et al., 2007).Finally, the NASA POWER estimates for both THI1 estimates and the average THI2 estimates had different scaling than the weather station estimates (regression coefficients somewhat different from 1, Table 2).This would make it difficult to combine weather data sets of these parameters from both weather sources into a single data set.The THI values for the phenotypic analyses were then calculated using only NASA POWER data.The replacement of weather station observations with NASA POWER data eliminated the restrictions imposed by weather station networks and potentially increased the number of herds that could be used in the following phenotypic analyses of heat stress.

Phenotypic Analyses of Heat Stress
The main factors that define a dairy cow's thermal environment are AT, humidity, WS, and solar radiation (Bianca, 1961).In summary, a high AT can increase the thermal load on an animal, as well as interfering with the dissipation of body heat, whereas high humidity interferes with evaporative heat loss by diminishing a cow's capacity for sweating and panting.Wind can also critically increase heat loss by increasing evaporative cooling through convection, and solar radiation can impose a substantial heat load on dairy cattle that are housed outside (Bianca, 1961).The THI is a heat stress index used to illustrate environmental heat load fluctuations.It is commonly studied alongside changes in production and fertility traits to measure the impact of heat stress on livestock.This index was originally developed to measure the combined effect of AT and RH at low WS on human comfort (Thom, 1959).Although THI has its limitations, it is advantageous compared with biological heat stress indicators such as body temperature and respiration rate, which can be expensive and difficult to measure in a large population.
It is evident from the results of these analyses that the productivity levels of Holsteins in all 3 parities from all 5 regions were negatively affected by increasing THI.The optimal thermal environment for maintaining homeostasis with the least amount of effort is rare, especially for domesticated animals confined to areas that are different from their natural habitats (Collier and Collier, 2012).Therefore, it is not surprising that dairy cattle across Canada, commonly raised in intensive housing systems, experience heat stress.The relationships between THI and the production traits for each region and parity are shown in Figures 2, 3, and  4. It is evident from these figures that each trait has a different pattern of response to increasing THI.Three THI thresholds were identified for milk production.For fat yield, 1 threshold was found, and 2 thresholds were identified for protein yield.For milk yield, the first threshold indicated a greater rate of increase in milk yield per unit THI on average, whereas the other 2 milk thresholds and the thresholds for fat and protein yield indicated a greater rate of decrease in the production traits per unit THI on average.
The reasons for and implications of the rate of increase in milk yield at low THI are unknown, but could be related to decreased water intake or to the effects of cold stress on milk yield.However, not enough evidence exists in the literature to support these speculations, and further research is required to effectively investigate this phenomenon.The differences in the patterns of production loss between the traits are also difficult to explain, but certain indirect and direct factors that contribute to the mechanisms of heat stress are discussed subsequently.For instance, a major factor that leads to a decline in milk production associated with heat stress is decreased feed intake.Decrease in feed intake causes subsequent changes in insulin action, glucose metabolism, adipose tissue mobilization, and skeletal muscle metabolism, which cause a negative energy balance.However, it is likely that decreased feed consumption is responsible for only about 35 to 50% of the decline in productivity associated with heat stress.Other major factors include changes in nutrient partitioning due to heat stress that occur independent of the amount of feed consumed (Sammad et al., 2020).There is less knowledge about the mechanisms involved in the decline in fat and protein yield due to heat stress (Becker et al., 2020).It is likely that a direct effect of heat stress, instead of reduction in feed intake, is more related to the observed decrease in protein yield.Carabaño et al. (2016) also found that the thresholds for milk components were lower than for total milk  yield.Ouellet et al. (2019) suggested that this higher sensitivity in Canada is related to the Canadian genetic selection program, which focuses on improving fat and protein yield because payment to dairy producers is based on milk component yields.This superior performance of fat and protein yield may confer a higher heat stress sensitivity.Second, it is important to provide adequate water supplies for dairy cattle due to their high milk production.A major physiological reaction to heat stress is increased water intake, to compensate for the water lost through milk, sweat, evaporation, and defecation (Becker et al., 2020).It is possible that providing adequate amounts of water during heat stress helps maintain cows' milk production for a lon-ger period (Ouellet et al., 2019).Overall, the direct and indirect effects of heat stress and the mechanisms responsible for the changes in productivity may differ between production traits, which could have resulted in the different thresholds and rate of changes per unit THI observed in this study.
The THI thresholds and average rates of change per unit THI for each trait are shown in Tables 3, 4, and 5.The first, second, and third thresholds for milk yield ranged between 47 to 50, 61 to 69, and 72 to 76 THI units, respectively.The regions with the lowest second and third milk thresholds were British Columbia and the Atlantic Maritime region.This may indicate that dairy cattle from these regions may be more sensitive to heat stress than cows from other regions.The regions with the highest milk thresholds were Quebec and the Prairies.The average effect on milk yield per THI unit above the second milk threshold ranged between −0.02 and −0.06 kg/d, whereas the average effect per THI unit above the third milk threshold ranged between −0.04 and −0.23 kg/d.All regions had similar declines in milk yield after the second threshold.However, British Columbia had a greater decline in milk yield after the third threshold relative to the other regions.
The THI thresholds for fat yield ranged between 48 and 55 THI units.All regions had a THI threshold for fat yield that was either 54 or 55 THI units, except for British Columbia.The average effect on fat yield per THI unit above the respective fat threshold ranged between −1.52 and −3.91 g/d.Heat stress caused a greater effect on fat yield per unit THI in dairy cattle from the Atlantic Maritime region.For protein yield, the first threshold ranged between 58 and 62 THI units, whereas the second threshold ranged between 72 and 73  3.Estimated rates of change (α) and SE for the average adjusted milk yields (kg) after 3 temperaturehumidity index thresholds (T1, T2, and T3) for all 3 parities in the 5 studied regions of Canada (QU = Quebec, ON = Ontario, BC = British Columbia, AM = Atlantic Maritime, and P = Prairies) THI units.All regions had similar thresholds for protein yield.The average effect on protein yield per THI unit above the first protein threshold ranged between −0.77 and −3.07 g/d, whereas it ranged between −0.96 and −7.17 g/d per THI unit above the second protein threshold.
The Atlantic Maritime region and British Columbia had greater rates of decline in protein yield per THI unit above the second threshold than the other regions, especially for second-and third-parity cows.In most regions, the estimated rates of decline per THI unit in these production traits were greater in parity 2 and 3 than in parity 1, which may indicate that multiparous cows are more sensitive to heat stress.Other studies have also concluded that later-parity cows may be more susceptible to heat stress (Aguilar et al., 2009;Bernabucci et al., 2014).This is likely due to their greater milk production potential than primiparous cows.Finally, the segmented polynomial models had a moderate to high R 2 value, indicating that the models adequately explained the variation in phenotypic data.The R 2 values ranged from 0.51 to 0.91, 0.92 to 0.99, and 0.90 to 0.98 for the milk yield, fat yield, and protein yield models, respectively.
The THI thresholds and subsequent rates of decline in production traits differed between regions.The reasons for this can only be speculated about in this study.However, regional variation in climate, housing and cooling systems, and nutritional management may have caused these differences.For instance, several cooling strategies, such as fans and sprinklers, can diminish the effects of heat stress in dairy cattle.Researchers have attempted to evaluate various cooling systems and have found that not all systems have the same efficiency (West, 2003).Furthermore, the thermal load on a dairy cow is affected by the type of housing system, stock density, ventilation, housing materials, and shading (Kadzere et al., 2002).Therefore, regional differences in the extents and types of cooling systems and housing structures could have resulted in different thermal strains at the same THI.Nutritional management may also cause variation in thermal strain at the same THI.For example, high water availability is essential for animals to cope well with heat stress.Dairy cattle in some regions may be more susceptible to heat stress if water resources are scarcer (West, 2003).Overall, it would be extremely beneficial to know more about which housing systems, cooling abatement strategies, and types of nutritional management are used for each herd, to explain the differences in THI breakpoints between regions.
Nevertheless, examining the climate of these 5 regions may help explain some of these differences.The average and maxima summer RH, AT, and THI from 2010 to 2019 for each region are shown in Figure 5. Dairy herds in Canada are found in 7 distinct ecozones, which are mapped in Figure 6.An ecozone is a classification used to define a region by abiotic and biotic ecological factors, including climate (Lands Directorate, 1986).The average winter and summer temperature ranges as well as the average annual rainfall for each ecozone are listed in Table 6.
British Columbia is composed of 2 ecozones, the Montane Cordillera and the Pacific Maritime.The Pacific Maritime has one of the warmest and wettest climates in Canada, whereas the Montane Cordillera ecozone is typically defined as having long cold winters and short warm summers (Lands Directorate, 1986).The region classified as the Prairies in the preceding analysis is also Rockett et al.: HEAT STRESS THRESHOLDS FOR CANADIAN HOLSTEINS Table 5.The estimated rates of change (α) and SE for the average adjusted protein yields (g) after 2 temperature-humidity index thresholds (T1, T2) for all 3 parities in the 5 studied regions of Canada (QU = Quebec, ON = Ontario, BC = British Columbia, AM = Atlantic Maritime, and P = Prairies) 1 composed of 2 different ecozones: the Boreal Plains and the Prairies.These ecozones can have a subhumid to semi-arid climate, due to the Rocky Mountains blocking moisture-bearing winds from West Canada.The Prairies region is also known for having high winds, with an average annual WS of 18 to 21 km/h (Lands Directorate, 1986).This increases evaporation and contributes to the dryness of the region.In contrast, the Atlantic Maritime ecozone is typically characterized as having a cool, moist, and moderate climate.Most dairy herds in Ontario and Quebec are found within 2 ecozones: the Mixedwood Plains and the Boreal Shield.The climate of the Boreal Shield ecozone can vary greatly but is usually described as having long cold winters and short warm summers.Herds in the Mixedwood Plains region experience cool winters and warm summers (Lands Directorate, 1986).Mixedwood Plains has the highest number of dairy herds in Canada, with 5,114 out of the 8,613 herds with known addresses being located within this ecozone.
The average number of days per year with a THI that negatively influences milk, fat, or protein yield between 2010 and 2019 are listed in Table 7. Overall, the average number of days per year with a detrimental THI varied depending on the trait as well as between regions.For milk yield, the average number of days per year with a THI between the second and third breakpoint ranged between 76 and 108 d.The average number of days per year above the third milk yield breakpoint ranged between 5 and 46 d.For fat yield, the average number of days per year with a THI above the fat breakpoint ranged from 180 to 211 d. moderate temperature.This is supported by the fact that weather within the Pacific Maritime ecozone does not vary a lot throughout the year.The effects of heat stress can be amplified if heat events are prolonged for multiple days (Hahn et al., 2009).Therefore, dairy cattle that are exposed to moderate THI values for a longer period may not be able to cope as well, resulting in lower THI thresholds.Furthermore, Brügemann et al. (2012) also found that thresholds for milk yield were lower in coastal regions.Similarly to British Columbia, these regions had the lowest percentage of production records at an extreme THI.The milk THI thresholds that Brügemann and colleagues found for herds in a coastal region, at pasture, and in an indoor environment were 67, 73, and 74, respectively.
The Atlantic Maritime had a higher summer average RH than other regions.As mentioned previously, dairy cattle in this region had a greater rate of decline in fat and protein yield compared with other regions.This could have been due to their diminished ability to lose heat through sweating and panting, especially at high THI values.On the contrary, the Prairies region typically has a drier climate with a low RH, high AT, and high WS.It is possible that dairy herds in this region are more likely to suffer from dehydration due to the low rainfall and high temperatures in the summer.Finally, Ontario and Quebec had very similar THI thresholds and production decays for all traits.This may be because most herds in these provinces are found within the same ecozone and experience a similar climate.Overall, the THI threshold at which an animal is no longer able to cope with heat stress may be determined by different environmental factors for each region, resulting in variations in THI breakpoints.
It is difficult to compare thresholds between studies because no single agreed method to calculate THI values exists.For instance, THI values can be calculated using different THI equations or by inputting hourly weather variables instead of daily averages.Furthermore, some other studies use average AT and RH to calculate THI values.Carabaño et al. (2016) found that average THI improved the statistical quality of their models compared with using the maximum values, whereas Ravagnolo et al. (2000) showed that different combinations of parameters had various correlations with production traits.They determined that the combination of maximum daily AT and minimum daily RH in THI was the best way to quantify heat stress (Ravagnolo et al., 2000).
Several limitations to quantifying heat stress in dairy cattle using the above methods also exist.For instance, this study does not fully account for the effects of prolonged heat stress or the lack of nighttime cooling.These factors could have influenced the results, as prolonged and intense heat events can amplify the negative effects of heat stress, and nighttime cooling is essential for dairy cattle to dissipate heat gained from the day (Holter et al., 1996).It is also difficult to interpret these results due to the lack of information available on the housing and cooling systems as well as the nutritional management of the individual herds.Herds that have implemented more extensive heat stress abatement strategies may skew THI thresholds and misrepresent the sensitivity of dairy cattle to heat stress.Finally, there may be a better indicator of heat stress other than the THI equation used in this study.Bohmanova et al. (2007) found that indices with a large weight on humidity were better predictors for cattle in regions with high humidity and indices with a small weight on humidity were better for cattle in low-humidity regions.Therefore, it may be beneficial to apply a different THI equation to regions that have higher average humidity, such as the Atlantic Maritime.regional variations in NASA POWER estimates should be further evaluated.This study also identified and examined several THI thresholds and rates of decline per unit THI for 3 production traits in primiparous and multiparous Holsteins from 5 regions in Canada.It was evident that all 3 production traits were negatively affected by heat stress in all regions of Canada.However, these traits did not have the same response to increasing THI, and the impact of heat stress varied between regions and parities.It was difficult to explain this regional variation, due to lack of information available on the management practices of each dairy herd.Finally, several days each year had THI above each trait's respective thresholds, indicating that dairy cattle could be experiencing heat stress during a large portion of the year.Although this study has limitations, it expands our current understanding of the effects of heat stress on Canadian Holsteins.

Figure 1 .
Figure 1.Locations of the weather stations (blue dots) in Canada used for the comparison between weather station data against National Aeronautics and Space Administration Prediction of Worldwide Energy Resources (NASA POWER) estimates.

Figure 2 .Figure 3 .
Figure 2. Relationships between temperature-humidity index (THI) and average adjusted milk yields (kg) as well as estimated THI thresholds (dotted lines) and fitted segmented polynomial functions (red lines).BC = British Columbia; AM = Atlantic Maritime region.

Figure 4 .
Figure 4. Relationships between temperature-humidity index (THI) and average adjusted protein yields (g) as well as estimated THI thresholds (dotted lines) and fitted segmented polynomial functions (red lines).BC = British Columbia; AM = Atlantic Maritime region.
Rockett et al.: HEAT STRESS THRESHOLDS FOR CANADIAN HOLSTEINSTable Figure 5. Average (solid lines) and maxima (dotted lines) ambient temperature, relative humidity, and temperature-humidity index (THI) values for each studied Canadian region from the months of June to August from 2010 to 2019.Quebec = red; Ontario = blue; British Columbia = pink; Atlantic Maritime = orange; Prairies = purple.

Table 1 .
Rockett et al.:HEAT STRESS THRESHOLDS FOR CANADIAN HOLSTEINS Numbers of records, cows, and herds, as well as means and SD for milk, fat, and protein yields for each parity by region of Canada studied (QU = Quebec, ON = Ontario, BC = British Columbia, AM = Atlantic Maritime, and P = Prairies)

Table 2 .
Rockett et al.: HEAT STRESS THRESHOLDS FOR CANADIAN HOLSTEINS Results from the linear regression analyses comparing National Aeronautics and Space Administration Prediction of Worldwide Energy Resources (NASA POWER) weather parameter estimates and temperature-humidity index (THI) values with the corresponding weather station values 1