Using publicly available weather station data to investigate the effects of heat stress on milk production traits in Canadian Holstein cattle

The institutional research ethics board of the University of Guelph did not require ethics approval for this study. Weather station data are publicly available from Environment and Climate Change Canada (www.ec.gc.ca) and the production data used in this study was made available by Lactanet (Guelph, ON, Canada; www.lactanet.ca), which are recorded by producers and producer organizations following the Code of Practice for the care and handling of farm animals — Dairy Cattle — of the Canadian National Farm Animal Care Council (https://www.nfacc.ca/codes-of-practice/dairy-cattle/code).

A total of 2.6 million test-day (TD) records for milk, fat, and protein yields from 170,439 Holstein cows raised in 1342 herds in Ontario, and 3.8 million records from 226,863 Holstein cows raised in 1962 herds located in Quebec, were provided by Lactanet. In Quebec and Ontario, the average size of herds is 76 and 95 cows, respectively, with mainly two types of barns — tie stall and free stall (CDIC 2021a). The data collected spanned a 10 yr period (from 2010 to 2019) and included observations from the first three parities. Only herds with records in at least five years and only observations collected within the interval from 5 to 305 d in milk (DIM) were used for analyses. In each parity, cows were required to have a minimum of four records.

Meteorological data were obtained from Environment Climate Change Canada (ECCC) using the “weathercan” (LaZerte and Albers 2018) package available in R (R Core Team 2021). The data consisted of hourly measurements (i.e., 24 records per day) of ambient temperature (AT) and relative humidity (RH). Each day was divided into three periods (i.e., morning (5 am–11 pm), afternoon (12 pm–8 pm), and night (9 pm–4 am)), and weather stations were required to have a minimum of one reading per period in each day over the entire 11 yr period. A total of 30 and 37 weather stations in Ontario and Quebec were linked to more than one farm. Therefore, the nearest weather station was assigned to each farm based on their latitude and longitude using the “geosphere” (Hijmans et al. 2019) package of R (R Core Team 2021). Weather records were collected from a total of 33 and 43 weather stations located within a maximum distance of 20 km from each herd (average distance of 13.1 ± 5.2 km and 12.6 ± 4.6 km) in Ontario and Quebec, respectively.

Records of AT and RH were incorporated into the National Research Council (1971) formula to calculate the THI as follows:

Two alternatives of THI were used in this study: (1) THI_avg, where daily average AT and daily average RH were used in the formula; (2) THI_max, where maximum AT and minimum RH were used. These two alternate THI allowed fair comparisons with other studies that have assessed heat stress.

The THI formula currently used does not include wind speed and solar radiation, which may influence heat stress conditions experienced by cattle, especially if animals are kept outside (e.g., feedlot cattle; Mader et al 2006; Gaughan et al 2008). However, considering a predominantly housing system (i.e, tie or free-stall) and the lack of pasture access in Canada (Beaver et al 2020), the influence of wind speed and solar radiation on heat stress conditions may be limited. Moreover, wind speed measurements from weather stations may not represent well the farm conditions as it has a greater dependency on local topography compared with temperature and humidity (Dunn et al. 2014).

In general, THI values of 2 d prior to the TD record have been reported as having the greatest effect on reducing milk yield and feed intake (West et al. 2003; Spiers et al. 2004; Bohmanova et al. 2008). Thus, considering the delayed effect of heat load on production traits, the value assigned to each test day, for both THI definitions (i.e., THI_avg and THI_max), was the average of the THI on the day of milk recording and the THI of two previous days (i.e., up to 48 h before milk recording). The total number of TD records per THI value by lactation, after merging the production data with the meteorological data, is shown in Fig. 1.

To identify the THI threshold in each parity at which milk, fat, and protein yields start to decline, the following steps were carried out: (1) a linear mixed model was fitted to adjust the phenotypes for known sources of environmental variation; (2) the adjusted phenotypes were plotted against the THI by lactation (first, second, or third) to provide visual information on the relationship between production and THI; and (3) a segmented polynomial was used to describe the shape of the curve describing the relationship between production and THI values, and the thresholds (i.e., breakpoints) for the polynomials were defined along the THI gradient.

In the first step, the analysis was performed for each trait (milk, fat, and protein yields) and each lactation (first, second, or third) using single-trait animal models using ASReml software (Gilmour et al. 2015). The linear mixed model used can be described as follows:

where y_ijklm is the mth test-day record (milk yield in kg, fat and protein yields in g) of the lth cow; μ is the overall mean; HY_i is the effect of the ith combination of herd and year of recording; AD_j is the jth combination of age at calving (eight classes for first lactation: = <24, 25, 26, 27, 28, 29, 30–31, >30 mo; six classes for the second lactation: = <36, 37, 38, 39–40, 41–42, >43 mo; and five classes for the third lactation: = <50, 51–52, 53–54, 55–56, >56 mo) and DIM (11 classes = 5–29, 30–60, 61–90, 91–120, 121–150, 151–180, 181–210, 211–240, 241–270, 271–300, 301–305 d); F_k is the effect of kth frequency of milking (k = 1–3); a_l is the random additive genetic effect of the lth animal; and e_ijklm is the residual term. Frequency of milking was included in the statistical model to account for the fact that animals that are milked more frequently tend to produce more milk. The residual term includes unknown sources of variation in y_ijklm, such as THI, in which , where is the predicted phenotypic observation given by model [1]. Therefore, the residuals of model [1] are the phenotypic observations adjusted for the known fixed environmental effects (i.e., HY, AD, and F) and the random additive genetic effect of the animal.

Hereafter, the residual term of model [1] will be referred to as adjusted phenotype (y_adj), which was plotted against the THI to display the relationship between production and THI. Previous knowledge about THI thresholds for milk production (e.g., Bernabucci et al. 2014; Carabaño et al. 2016) and visual inspection of the plot provided information on the most likely number and location of the breakpoints, i.e., where the pattern of milk, fat, and protein yields would change considering the THI gradient. Visual inspection and previous knowledge of the data to be fit have been used to determine the number and location of breakpoints for segmented polynomials (Madár et al. 2003; Misztal 2006; Gorniak et al. 2014). A total of three, two, and one breakpoints were used for the segmented polynomials for milk, protein, and fat yield, respectively, to determine the THI thresholds for each trait and lactation.

Segmented polynomials describe a non-linear relationship between variables by connecting different functions through the breakpoints (Rice 1969). Thus, with a prior determination of the number and regions of breakpoints and the adjusted phenotype from model [1], segmented polynomials allowed to find THI thresholds where production started to decrease. Segmented polynomials with different linear functions were fitted using the PROC REG procedure of the SAS software (version 9.2; SAS 2011), and different breakpoint locations were tested. The general segmented polynomial model used in this study can be described as follows:in which elements of z are given by:

where y_adj is the vector of adjusted phenotypes obtained from model [1]; α is the intercept; 1 is the vector of ones; β₁ is the linear regression coefficient; x is the vector of THI values; γ₁,…, γ_n are linear regression coefficients, which give the differences in the slopes in relation to the preceding linear segments; z_1i,…,z_ni are covariates as defined above; and C₁,…, C_n is the breakpoint values. The segmented polynomial model used for milk, protein, and fat yield had three, two, and one breakpoints, respectively. The final location of the breakpoints (i.e., THI value) was determined based on the adjusted coefficient of determination (R²) of the model and the statistical significance (P value) of the segments describing the change in the adjusted phenotypes. The linear regression coefficient for each breakpoint (γ) gives the rate of change in production, from which production loss can be estimated.

A simplified estimation of economic loss (EL) associated with decreased fat and protein yield was calculated as follows:

where comp$ is the milk component price, which was calculated assuming the average prices paid to producers based on fat or protein yields (i.e., $8.28·kg⁻¹ for butterfat and $6.38·kg⁻¹ for protein yield) following the 2019–2020 Dairy Farmers of Ontario’s Annual Report (Dairy Farmers of Ontario 2020), HS_days is the average number of days above the THI threshold per year, γ is the average decline rate (converted to kg·unit⁻¹ of THI_max) estimated for the three parities and both provinces of either fat or protein yield, HS_THI is the average THI_max value above the THI threshold, and THI_threshold was set to THI_max = 58, based on the average of thresholds identified for fat and protein yields. A single THI threshold for fat and protein was assumed because the average THI threshold identified for each component across lactations did not differ considerably (<5 THI units). The EL for protein and fat yields was then multiplied by the number of milking dairy cows in Ontario and Quebec, and the resulted value summed up for the total EL.

The average production in Ontario and Quebec, for each trait in the three lactations, is shown in Table 1. Ontario and Quebec comprise the Central Canada region, characterized by a humid continental climate (Dfb; Ontario) and continental subarctic climate (Dfc; Quebec), according to the Köppen–Geiger system (Peel et al. 2007), with Ontario having warmer and longer summers compared with those in Quebec. These provinces have a similar level of production and are major contributors to the dairy cattle industry in Canada. Over the studied years, daily production averages considering both provinces were 32.07 kg (±8.25), 1312.76 g (±315.42), and 1095.87 g (±229.01) of milk, fat, and protein yield, respectively (Supplementary Table S1¹). These average productions were similar to the average production in other Holstein populations used to study the effect of heat stress worldwide (e.g., Bohmanova et al. 2007; Bernabucci et al. 2014; Carabaño et al. 2016).

Physiological mechanisms determined mainly by the endocrine system are responsible for maintaining the animal’s fitness in the environment (Roy and Collier 2012). Under environmental conditions that impair heat exchange, the animal’s body temperature increases, and physiological responses are triggered. Two important factors that impact the heat exchange between the environment and an animal are AT and RH (Hahn et al. 2009). Despite the predominance of cold weather in Canada, the summer months (from June to September) were characterized by elevated AT and high RH, with maximum AT of 37.6 and 36.1 °C in Ontario and Quebec, respectively (Table 2). In both provinces, cows were exposed to AT higher than 28.4 °C, which is considered as the upper limit at which they can maintain a stable body temperature (Dikmen and Hansen 2009). Moreover, the RH levels also exceeded the upper limit of 70%, which is needed for adequate cow comfort, as above this limit, evaporative cooling systems are hindered (Maust et al. 1972; Berman 2006). Therefore, the AT and RH levels reported in this study suggest that Canadian cows raised in Ontario and Quebec might experience heat stress during the summer months.

During the summer periods between 2010 and 2019, the average THI_max (THI_avg) was 70.0 (64.6) and 67.9 (61.9) in Ontario and Quebec, respectively. Assuming a THI_max threshold of 62 based on the average threshold for milk, fat, and protein yields found in the current study, cows most likely experienced heat stress (THI_max > 62) in, on average, 115 and 106 d per year in Ontario and Quebec, respectively (Fig. 2).

The use of meteorological information from the nearest weather station to each herd is a limitation for the accuracy of describing on-farm climatic conditions. The maximum distance allowed between the weather station and a farm varies among studies. For instance, the maximum distance was 5 km in Italy (Bernabucci et al. 2014), 60 km in Australia (Nguyen et al. 2016), and 60 km in Germany (Brügemann et al. 2012). Generally, the further away a weather station is from a farm, the less likely the weather station will accurately predict the weather conditions on-farm. However, it is important to consider that weather stations very close to the farm are not always available. Studies on the relationship between climatic conditions on-farm and measures from the closest weather station have found that THI calculated with weather station data underestimates the heat stress occurring within the barn (Schüller et al. 2013; Shock et al. 2016; Ouellet et al. 2019a). In this regard, studies showed that when using weather station data as a surrogate to on-farm data, THI_max describes better the within barn conditions compared with THI_avg (Ravagnolo et al. 2000; Ouellet et al. 2019a).

The observed relationship between the adjusted phenotype and THI values for the three lactations and three traits are shown in Figs. 3–5. Additionally, the number of observations available per THI_max (THI_avg) are presented in Supplementary material (Supplementary Tables S1–S4¹). A decline in production as THI values increased was observed after the first threshold for all traits except milk yield (an increase in milk production was observed from the first to the second threshold). The first THI_max (THI_avg) thresholds after which milk yield declined ranged from 67 (64) to 71 (66) in Ontario and 64 (61) to 69 (63) in Quebec across lactations. A second threshold describing a further severe decline in milk yield was observed when THI_max (THI_avg) values exceeded 76 (74) in both provinces across the three lactations. For fat yield, a single THI_max (THI_avg) of 56 (50) in Ontario and 58 (50) in Quebec was found, at which the production started to decrease. For protein yield, the first THI_max (THI_avg) threshold ranged from 57 (57) to 62 (60) in Ontario and from 59 (55) to 62 (58) in Quebec across lactations. A second THI_max (THI_avg) threshold, describing a greater decline in protein yield, was similar for the two provinces and ranged from 77 (71) to 80 (74) across lactations. Protein and fat yield had a similar response pattern to THI with one or two thresholds describing a negative slope, whereas for milk yield, the first breakpoint described a positive slope. Cows under heat stress respond to evaporative heat loss requirements with increased water intake (Hall et al. 2018; Collier et al. 2019), which may partially explain the observed increase in milk at relatively high levels of THI compared with the other two components. The reduction in milk fat and protein production has also been reported in other studies (Gorniak et al. 2014; Carabaño et al. 2016; Ouellet et al. 2019a).

In the present study, the THI_max (THI_avg) thresholds estimated for Ontario and Quebec, at which milk yield started to decline, was on average 68 (64) across the three lactations. Under comparable environmental conditions, THI_max (THI_avg) thresholds were identified at 75 (69) in Belgium and 73 (62) in Luxembourg (Hammami et al. 2013; Carabaño et al. 2016), at THI_avg = 65 in Germany (Pinto et al. 2020), and THI_avg = 68 in Poland (Herbut et al. 2015). Differences in heat stress thresholds across studies may be due to many factors, which can be environmental in nature, such as climate, measurement conditions (climate chambers vs. commercial conditions), or related to animals, such as parity, production level, and breed.

The average THI_max (THI_avg) thresholds across the three lactations for fat and protein yield were 57 (50) and 60 (58) in Ontario and Quebec, respectively. Therefore, lower THI thresholds were observed for milk components compared with milk yield, indicating greater sensitivity of milk components to heat stress. This finding agrees with other studies, which also found that milk components are more sensitive to heat stress than milk yield (e.g., Hammami et al. 2015; Carabaño et al. 2016; Ouellet et al. 2019b). Changes in milk composition during heat stress are associated with decreased milk protein synthesis and could partly explain this greater sensitivity (Hu et al. 2016). Furthermore, studies have shown that increased THI values are associated with increased somatic cell score (SCS) (Hammami et al. 2013; Lambertz et al. 2014; Negri et al. 2021), which need to be further evaluated in Canada. In general, the greater sensitivity of milk components to heat stress can be an additional challenge for the Canadian dairy industry, as farmers are currently paid based on fat and protein yields.

The THI thresholds identified for fat and protein were higher in Quebec than in Ontario. However, the opposite was observed on the THI thresholds for milk yield. Possible explanations for this fact are related to factors such as types of barns, stocking density, and level of productivity, which differ considerably between the two provinces. For instance, in Quebec, more than 90% of the barns are tie-stall barns, whereas in Ontario, this type of barn represents less than 70% (CDIC 2021a, 2021b). Improving cow’s comfort to ensure an acceptable level of welfare is paramount, especially for producers using a tie-stall system, as there is a positive association between welfare and higher productivity and profitability in this farm (Villettaz Robichaud et al. 2019). This might indicate that a higher number of farms in Quebec are equipped with some type of cooling system when compared with Ontario. Further studies should consider the factors altering the micro-climate within the barn in each province. Moving forward, it is recommendable that barn climatic conditions should be recorded and stored in dairy cattle breeding databases for future studies.

Assuming the THI threshold of 58 for fat and protein yields, the average number of days per year exceeding this threshold in both provinces was 156 d, with an average of 10 THI units above the threshold. As mentioned before, according to THI thresholds estimated in this study, milk components are more sensitive than milk yield to heat stress. Thus, these components may be affected by high temperature and humidity in a greater number of days in both provinces. The THI_max and THI_avg thresholds and their respective average productive losses associated with heat stress for the two provinces and three traits are shown in Table 3 for Ontario, and in Table 4 for Quebec. The rate of decline is expressed in kg·d⁻¹ per unit of THI above the threshold for milk yield and in g·d⁻¹ per unit of THI above the threshold for fat and protein yield. There was no considerable difference between the estimated losses from THI_avg and THI_max, except for the case of extreme THI values (e.g., THI_max (THI_avg) = 80 (71)) for protein yield. Collier and Zimbelman (2007) and West et al. (2003) also reported a linear decline in milk yield above the THI threshold.

In Ontario, the average impact of a unit of THI_max or THI_avg above the first threshold at which a decrease in milk yield was observed ranged from −0.05 to −0.08 kg·d⁻¹. Severe heat stress was picked up by a second threshold, which resulted in milk yield losses ranging from −0.16 to −0.31 kg·d⁻¹. In Quebec, the average decline rate ranged from −0.01 to −0.05 kg·d⁻¹ per unit of THI_max or THI_avg, with a greater decline observed after the second THI threshold ranging from −0.21 to −0.48 kg·d⁻¹. With comparable average daily milk production, milk losses due to heat stress were estimated to be up to −0.14 kg·d⁻¹ per unit of THI_max or THI_avg in Belgium (Carabaño et al. 2016). In Italy, milk losses ranged from −0.91 to −1.27 kg·d⁻¹ per unit of THI_max for the first, second, and third lactation (Bernabucci et al. 2014). In Germany, based on different production systems, the decline of milk yield ranged from −0.16 to −0.47 kg·d⁻¹ per unit of THI_max (Brügemann et al. 2012).

The rate of decline found for protein yield ranged from −1.81 to −2.57 g·d⁻¹ per unit of THI_max or THI_avg above the threshold in Ontario, and from −1.07 to −1.98 g·d⁻¹ in Quebec. Maximum declines above the second THI_max (THI_avg) threshold of −10.83 and −9.08 g·d⁻¹ were observed in Ontario and Quebec, respectively. Average estimated fat yield losses ranged from −2.09 to −2.40 g·d⁻¹ per unit of THI_max or THI_avg in Ontario, and from −1.79 to −2.28 g·d⁻¹ in Quebec. The results of the current study were similar to losses estimated in Europe. Average protein losses were reported ranging from −2.6 to −8.4 g·d⁻¹ per THI_max unit, with losses up to −12.1 g·d⁻¹ per THI_max unit when the top 1% of animals were considered (Carabaño et al. 2016). In the same study, fat yield losses of −2.5 to −5.1 g·d⁻¹ per unit of THI_max or THI_avg were observed at extreme THI values (i.e., THI_max(THI_avg) = 72 (82)). By contrast, Hammami et al. (2013) estimated fat and protein yield losses of −20 and −13 g·d⁻¹ per unit of THI_avg, respectively. In Italy, protein and fat yield losses for the first three lactations were on average −50 and −30 g·d⁻¹ per unit of THI_max, respectively (Bernabucci et al. 2014). In another study using the Canadian Holstein cow population, losses of −20 g·d⁻¹ were observed for fat yield when cows were exposed to one to two consecutive days of heat stress (Ouellet et al. 2019b). It seems that the rate of decline is quite sensitive to a number of factors, such as different THI scaling and thresholds (Bohmanova et al. 2007), the period of time that cows are exposed to high THI (Ouellet et al. 2019b), and different production systems (Brügemann et al. 2012).

Parity differences for heat tolerance were observed, with first-parity cows having lower average rates of decline for the three traits when compared with second- and third-parity cows. This result is in line with other studies of heat stress that evaluated animals across the first three lactations (Aguilar et al. 2009; Bernabucci et al. 2014; Gantner et al. 2017), showing that first-parity cows exhibited lower sensitivity to heat stress. This lower sensitivity can be explained by the different metabolic states of the animals, which results in lower milk yield of first-parity cows compared with multiparous cows (Wathes et al. 2007).

Assuming the THI threshold of 58 for fat and protein yields, cows would be under heat stress on an average 156 d per year, with an average of 10 THI units above the threshold. Considering the average decline per unit of THI exceeding the threshold for each component across lactations and provinces, the overall reduction in fat and protein yields was estimated to be 2.98 and 3.79 kg·cow⁻¹ yr⁻¹, respectively. Assuming the average prices paid to producers (i.e., $8.28·kg⁻¹ for butterfat and $6.38·kg⁻¹ for protein yield; Dairy Farmers of Ontario 2020), the impact of heat stress on the productivity of cows results in a potential EL of $50.40 per cow per year. For an average annual production of a Holstein cow (i.e., 426.5 kg of fat and 354 kg of protein yield, Agriculture and Agri-Food Canada 2021), and the price paid to producers, the EL represents about 1% of the total income per cow per year. Considering the total number of cows (i.e., 222,129 in Ontario and 247,166 in Quebec), an EL of $34.5 million per year might occur due to heat stress in Ontario and Quebec.

In the United States, St-Pierre et al. (2003) estimated averaged annual ELs due to heat stress ranging from $897 million to $1.5 billion. Studies have also estimated annual losses associated with milk yield reduction if dry cows were exposed to hot environmental conditions without cooling to be $87 per cow (Ferreira et al. 2016), and losses of daughters born from dams exposed to heat stress during late gestation to be $39 per daughter (Laporta et al. 2020). The differences in EL estimates are mostly related to the size of the respective dairy industry, the method used to estimate the impact, and the different environmental conditions across countries (e.g., number of days over the thermal comfort). In the current study, only the decrease in fat and protein yields was considered to estimate the ELs, which results in a partial estimation of the economic impact of heat stress. Future studies on ELs due to heat stress in the dairy industry should account for factors, such as disease incidence (Nardone et al. 2010), reproductive performance (Pinedo and de Vries 2017), cow mortality (Bishop-Williams et al. 2015), and other important economic traits, such as SCS (Negri et al. 2021).