Introduction
As climate change intensifies, livestock industries face unprecedented environmental threats. Given the antagonistic relationship between production efficiency and environmental resilience (Roy and Collier 2012), heat stress has emerged as the major climatic challenge facing the dairy industry. Elevated temperatures coupled with high humidity impose substantial physiological challenges on dairy cows, which are further exacerbated by the high metabolic heat production of lactating cows (Becker et al.
2020). Because milk synthesis increases metabolic heat production, lactating cows are especially sensitive to heat stress (Tao et al. 2020; Yan et al.
2021). High-producing cows are further challenged by greater internal heat load, limiting their ability to maintain core body temperature and increasing susceptibility to heat stress (Collier et al. 2019).
Moreover, heat stress not only reduces milk yield and alters composition, but it can also impact immune function (Lengi et al. 2022) and udder health (Bertocchi et al. 2014; Lambertz et al.
2014; Nasr and El-Tarabany 2017). Milk somatic cell count (SCC) is a widely used indicator of subclinical udder inflammation in dairy cattle (Dohoo and Meek 1982) and has shown to be the most impactful variable when detecting subclinical mastitis (Pakrashi et al. 2023).
In genetic studies SCC is typically converted to somatic cell score (SCS) via a logarithmic transformation to normalize its distribution (Norman et al. 2000; Moore et al. 2024).
Previous studies have found an association of SCS with the temperature–humidity index (THI) and determined THI thresholds at which an increase in somatic cell levels is observed due to heat stress (Brügemann et al. 2012; Hammami et al. 2013; Negri et al.
2021). This further indicated a potential role of high temperatures and humidity on the incidence of subclinical mastitis. Negri et al.
(2021) found a THI threshold for SCS lower than the thresholds for milk production traits and suggested SCS as an early indicator of heat stress for dairy cattle in tropical environments. However, under different environmental conditions, Hammami et al. (2013) found a higher THI threshold for SCS compared to milk production traits.
Brügemann et al. (2012) also reported differences in the relationship between somatic cell levels and THI across production systems. Furthermore, other studies have not observed an increase in somatic cell levels under heat stress (Smith et al.
2013; Lengi et al. 2022) or have found differences between housing systems (Lambertz et al. 2014), which suggests that the increase in somatic cell levels is likely related to a multiple factors, such as management, housing type, and not exclusively to heat load.
Given the detrimental effects of heat stress on dairy cattle and the potential relationship between heat stress and somatic cell levels, it becomes important to investigate this association to improve herd management strategies and mitigate economic losses in the dairy industry. Alongside the implementation of heat-abatement measures, genetic selection can help mitigate heat stress effects by identifying and selecting animals that are genetically more heat tolerant. Combining THI with large scale test-day records provides an opportunity to quantify the impact of heat stress on somatic cell levels and to generate crucial information for performing genetic evaluation and estimation of breeding values for heat tolerance.
Estimated breeding values for heat tolerance of milk production traits have been successfully implemented worldwide (Nguyen et al. 2016; Finocchiaro et al. 2022; Van Kaam et al.
2025). Genetic analyses of heat stress typically use reaction norm (RN) models to account for an individual’s performance on a continuous environmental descriptor. Many studies have modeled the relationship between the traits and the environment, using THI as an environment descriptor (Ravagnolo et al.
2000; Bernabucci et al. 2014).
In Canada, initial studies have assessed the impact of heat stress and performed genetic analyses for milk production (Campos et al. 2022; Rockett et al. 2023a) and fertility traits (Dodd et al.
2026). Re-ranking of bulls’ EBV was observed for these traits under thermal comfort and heat stress conditions, indicating a genotype by environment (GxE) interaction. Building on this work, the objective of the present study was to assess the relationship between SCS and THI in different regions across Canada and to estimate genetic parameters to evaluate the potential for genetic selection for heat tolerance of SCS.
Materials and methods
Approval from the Animal Care and Use Committee was not needed for this study since the data were obtained from existing databases made available by Lactanet (Guelph, ON, Canada; www.lactanet.ca).
Somatic cell score and weather data
A data set consisting of approximately 7.6 million test-day records of SCC recorded between 2018 and 2021 from 554,857 first-parity Holstein cows of 4,591 herds across Canada was provided by Lactanet Canada (Guelph, ON). The data set was divided into five Canadian regions which were: Ontario, Quebec, British Columbia, the Prairies (which included the provinces Saskatchewan, Manitoba, and Alberta), and the Atlantic Maritime (which included the provinces Newfoundland and Labrador, Nova Scotia, Prince Edward Island, and New Brunswick).
Using the location of each herd (ie, latitude and longitude), hourly ambient temperature (AT, °C) and relative humidity (RH, %) were retrieved from the NASA POWER database (https://power.larc.nasa.gov/) for the period studied, using the latitude and longitude coordinates of each herd to match the meteorological data to each location. The hourly THI was calculated using the NRC (1971) formula as follows:
To merge production data set with the weather data set, the hourly THI was averaged to obtain a daily THI, and the average of the THI on the three previous days (Hammami et al. 2013) was assigned to each test-day record. Cows were required to have a minimum of three test-day records collected within the interval between 5 and 305 DIM and were removed from the analysis if they had records from multiple herds.
Herds were required to have records throughout all years, and herd-year combination was required to have a minimum of five records.
For the genetic analysis, pedigree information was available for all cows with SCC records and included 1,182,593 animals. The SCC was transformed into SCS using the following formula: SCS = log2(SCC/100,000) + 3 (Ali and Shook 1980). The final number of records, animals, and herds used in the analyses, as well as the averages and standard deviations of SCS per region, are shown in Table 1.
Relationship between heat stress and SCS and THI threshold determination
A linear mixed model was used to adjust the phenotypes for known sources of environmental variation to investigate the relationship between THI and SCS. This model was applied to each region separately and performed using BLUPF90+ program with the AIREML option (Lourenco et al. 2022).
Additionally, a data set containing all regions was also analyzed. The univariate linear mixed model used can be described as follows:
where yijklm is the mth test-day SCS record of the lth cow; HYi is the fixed effect of the ith combination of herd and year of recording; Aj is the fixed effect of jth age at calving (8 classes: = < 24, 25, 26, 27, 28, 29, 30–31, > 30 mo); DIMk is the fixed effect of kth days in milk (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); cl ∼ N(0, Aσc2) is the random effect of the lth animal; pel ∼N(0, Iσpe2) is the random permanent environmental effect of the lth animal; and eijklm∼ N(0, Iσe2) is the residual term. Where A is the numerator relationship matrix, and I is an identity matrix. The residual term includes unknown sources of variation not accounted for in the predicted observation (y^ijklm), which includes the THI.
Therefore, the residual of model [1] is referred to as the adjusted phenotype used in the subsequent analyses of the relationship between SCS and THI.
To describe the relationship between SCS (ie, adjusted phenotype) and THI and to determine the threshold (ie, breakpoint) at which SCS increases, the adjusted phenotype was plotted against the THI. To ensure enough data at the upper extreme of THI, the three highest THI values in each region were grouped into a single THI value. Adjusted phenotypes were smoothed using the “rollmean” function from the “zoo” package in R (R Core team 2020), applying a rolling average of three neighboring THI values to reduce data noise.
The THI thresholds were estimated using the “segmented” package (Muggeo et al. 2003), with initial values informed by visual inspection of the adjusted SCS over the THI. The segmented model assumes a continuous piecewise linear relationship, allowing for a change in slope at the estimated breakpoints.
For the estimation of the THI thresholds, throughout the entire period studied, only records for above THI 30 were considered. This THI value was set as it represented the average minimum THI value observed during the summer months (June, July, August, and September), and low THI values would not contribute to heat stress inferences and assessment of the relationship between THI and SCS.
Additionally, this investigation also analyzed SCS in the early lactation divided in three different data sets: 5 to 30 DIM (SCS30), 5 to 60 DIM (SCS60), and 5 to 100 DIM (SCS100). The early postpartum period of dairy cows is associated with reduced immune competence shown by the highest incidence of diseases, most of which related to udder health (De Vliegher et al. 2004; Carvalho et al.
2019). Thus, the relationship between SCS and heat stress indicated by the THI may be intensified. In this case, the fixed effect of DIM in model [1] was separated into five classes of 6, 12, and 19 d for SCS30, SCS60, and SCS100, respectively.
Genetic analyses of SCS considering heat load
After examining the behavior of SCS across the THI values, genetic analysis was carried out on a data set that included cows from all herds. The trait used in the analysis was SCS100 which was considered the best trait to capture the impact of heat stress and ensured enough records per cow. A repeatability test-day RN model was used fitting a heat stress function (HSf).
The model can be written as:
where yijklm is the mth test-day record of SCS100 of the lth cow; HYi is the effect of the ith combination of herd and year of recording; Aj is the jth fixed effect of age at calving (eight classes: ≤ 24, 25, 26, 27, 28, 29, 30–31, > 30 mo); DIMk is the fixed effect of the kth days in milk (19 classes: 5–10, 11–15, 16–20, 21–25, 26–30, 31–35, 36–40, 41–45, 46–50, 51–55, 56–60, 61–65, 66–70, 71–75, 76–80, 81–85, 86–90, 91–95, 96–100 d); γ(THI) is the fixed linear regression on the THI fitted as covariate; cl and pl are the random additive and permanent environmental effects of the lth animal under thermal comfort (ie, intercept); αl and δl are the random additive and permanent environmental effects of the lth animal under heat stress (ie, slope), and HSf is given by max[0,(THIlm-τ)], where τ is the THI threshold, and THIlm is the THI assigned to the mth test-day of the lth cow; and eijklm is the residual term. The random terms were assumed to be normally distributed with mean zero and (co)variance structure specified as:
where, I is an identity matrix, A is the numerator relationship matrix, σc2 and σα2 are the variances of the general additive genetic effects and the additive genetic effects on the slopes, respectively, σp2 and σδ2 are the variances of the general permanent environment effects and the permanent environment effects on the slopes, σca is the covariance between the general additive genetic effects and additive genetic effects on the slopes, respectively, σpδ is the covariance between the general permanent environment effects and permanent environment effects on the slopes, and σe2 is the error variance. The heritability estimates after the THI threshold were obtained as follows:
Genetic correlation between additive genetic effect under thermal comfort and additive genetic effect under heat stress was computed as: corr(c,α)=σc,ασc2σα2.
Bivariate analyses were also performed pairing SCS100 with either milk, fat, or protein yields in the whole lactation (5–305 DIM) including their heat tolerance component. In the bivariate analyses, model [2] was used for both traits, but with the fixed effect of DIM on milk production traits categorized into 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. The random terms of the bivariate model were assumed to be normally distributed with mean zero and (co)variance structure specified as:
where C is a vector containing the additive genetic effects cl and αl, and P a vector containing the permanent environmental effects pl and δl; matrices G, P, and R for SCS100 and milk, fat, or protein yields (yield) can be described as:
For production traits, the THI thresholds were 65, 55, and 60, for milk, fat, and protein yields, respectively, based on previous studies conducted in Canada on the same Holstein population (Campos et al. 2022; Rockett et al. 2023b).
Additionally, genetic correlation between additive genetic effect under heat stress between traits was computed as: corr(αscs,αyield)=σαscs,αyieldσαscs2σαyield2.
Relative breeding value and G × E interaction
The EBVs for SCS100 under heat stress obtained from the univariate analysis were calculated for different THI levels as EBVl = c^lscs + HSf * α^lscs, where c^lscs and α^lscs are the estimated solutions for the additive effect under thermal comfort and additive effect of heat tolerance (ie, slope) of the lth animal, respectively. In Canada, SCS is expressed in relative breeding values (RBV) which is the EBV standardized to have a mean of 100 and standard deviation of 5 and reversed in sign. The RBV of all bulls that had ≥ 50 daughters with phenotypes (879 bulls) and the top 10% bulls based on their ranking under thermal comfort were compared under heat stress conditions (20 units above the THI threshold) using Spearman rank correlation to investigate the presence of GxE interaction.
Results and Discussion
Effect of heat stress on SCS across Canada
Different climate profiles are encountered in Canada, which indicates that cows are exposed to different heat stress conditions throughout the summer (Table 2). The average THI ranged from 58.9 to 66.0, with maximum THI reaching 78. Considering the daily average THI during the summer for all regions, about 40% of the days were above THI 60 which exceeds the THI thresholds found for milk production traits (Campos et al.
2022). In Ontario and Quebec, this corresponded to an average of 109 and 124 days per year, respectively, above this threshold. These two provinces account for 70% of the national dairy production (CDC 2026), which indicates that most of the cows in Canada are exposed to heat stress conditions during the summer.
The adjusted SCS throughout the whole lactation (5 to 305 DIM) showed a marginal increase associated with heat stress conditions in all regions, except the Prairies (Figure 1). In the Prairies, no clear relationship between SCS and THI was observed, suggesting no effect of heat stress. This region presents a distinct climate profile characterized by low RH and mild AT throughout the summer, which might have influenced the effect of heat stress observed for the whole lactation SCS.
Historical climate studies in Canada have shown no significant upward trend in temperatures within the Prairies when compared to the broader national context (Vincent et al. 2018; VanderZaag et al. 2023a).
When all regions were investigated together, an upward trend of SCS over the THI was also observed (Figure 2). While some individual regions showed a decrease in SCS at high THI, this was not observed when all regions were combined. This inconsistency indicates that the apparent regional declines likely resulted from the limited number of observations at the upper extreme of the THI distribution.
The relationship between SCS in early lactation (SCS30, SCS60, and SCS100) and the THI showed that an increase in THI was associated with an increase in SCS across all regions (Figure 1). This association between heat stress conditions and SCS was more pronounced during early lactation, where highest SCS was observed, compared to the whole lactation (Figure 1). In the Prairies and Quebec regions, where a clear association between SCS over the entire lactation period and heat stress conditions was not evident, an increase in SCS30 with increased THI could be observed.
Interestingly, Lambertz et al. (2014) found a higher increase in SCS at THI > 60 associated with late (201 to 305 DIM) compared to early (0 to 100 DIM), and mid (101 to 200 DIM) lactation. In the present study, the relationship between SCS in mid or late lactation and THI did not differ substantially from SCS in the whole lactation (results not shown).
During early lactation, cows experience a distinct physiological state and respond differently to heat stress compared to mid- and late-lactation cows (Tao et al. 2020). Additionally, nutritional deficits and metabolic conditions during the early postpartum period can impair the immune system (Sordillo 2016), which can result in various health disorders (Carvalho et al.
2019). The results of this study also suggest that heat stress can further amplify the physiological challenges during this period of lactation indicated by the higher SCS values observed.
The THI values at which SCS100 started to increase ranged from 52 to 60 across regions, with an increase in SCS per unit of THI ranging from 0.005 to 0.012 (Table 3). When all regions were considered together, a greater increase in SCS100 was observed from THI 54. Reported THI thresholds for increased SCS under heat stress vary across studies, reflecting differences in regional climates and cow populations.
In Europe, THI thresholds at 57 with an increase of 0.013 per unit of THI (Italy, Bertocchi et al. 2014) and 60 and 66 with marginal increases (Germany, Brügemann et al. 2012; Hammami et al.
2013) were found. In South Korea a THI threshold at 66.7 was identified with an increase of 0.003 per unit of THI (Lee et al. 2023).
Under a tropical environment, the threshold for SCS was identified at THI of 70, with an increase rate of 0.09 per unit of THI (Brazil, Negri et al. 2021).
The thresholds determined for SCS100 in this study were lower than THI thresholds found in other studies that look at complete lactation. However, most of the studies were conducted in warmer climates than Canada. Cold winters and hot summers impose a short-term challenge on cows’ thermoregulation, in contrast to places with a more uniform climate profile (ie, long-term challenge).
Acclimatization plays a role in how cows are affected by heat stress conditions (Collier et al. 2019), which may translate into the different THI thresholds observed. Moreover, management practices that results in lower SCS averages across the country may also influence the observations in this study.
In Canada, the selection for low SCS associated with increased mastitis resistance has been taking place since 2013 (CDN 2013). The studied Canadian Holstein population has an average SCS of 1.80. This is considerably lower than the average values of other populations, for example, mean SCS 2.22 (Hammami et al.
2013) and 3.30 (Negri et al. 2021), which may result in a greater sensitive to heat stress effects compared to these other populations. Nonetheless, the THI thresholds found in the present study compare with those identified for the same Canadian Holstein population for milk components, particularly fat yield, which were lower than for milk yield (Campos et al.
2022). This indicates that, similar to milk components, SCS100 shows earlier the effect of heat stress in Canadian Holsteins.
Multiple studies have investigated the effect of heat stress on milk composition and udder health (Moore et al. 2023; Lemal et al. 2025).
Although various factors at the farm level, such as management, housing type, and nutritional programs, can contribute to seasonal changes in somatic cells, the increase during the summer can be largely attributed to heat stress (Tao et al. 2018). Across Canada, with THI levels reaching 70 in early lactation, an overall increase of about 6% in SCS is expected.
Notably, this estimate represents only the lower bound of the expected increase in SCS, as the THI used in this study was based on estimated outdoor conditions, whereas THI conditions inside the barn are typically higher (Mylostyvyi et al. 2020; VanderZaag et al. 2023b).
While Canada is mostly a cold climate country, its temperature is increasing twice as fast as the global average (Bush and Lemmen 2019; ECCC 2025). Therefore, it is essential to continue investigating the effects of heat stress in Canadian dairy cattle.
Genetic parameters for somatic cell score under heat stress
The genetic parameters estimated for SCS100 under thermal comfort and heat stress conditions are present in Table 4. The total additive genetic variance showed an increase relative to the permanent environmental variance at high THI levels, which resulted in an increase in the heritability estimates. The genetic correlation between the general additive genetic effect and additive genetic effect under heat stress was 0.42 (±0.29).
This moderate positive genetic correlation between intercept and slope indicates that bulls with high SCS100 under thermal comfort are likely those with high SCS100 under heat stress. Thus, selecting bulls with low SCS under thermal comfort might also lead to daughters that are genetically more heat tolerant for SCS.
The heritability estimates slightly increased as the THI increased. Under thermal comfort, heritability was 0.08 (±0.01), which was slightly lower compared with previous estimates using different models in Canada, which ranged from 0.11 to 0.18 (Koeck et al. 2012; Lactanet 2021; Interbull 2023).
Under heat stress, the heritability increased as THI increased and reached 0.10 (±0.01). As heat stress intensifies, there is a higher proportion of the phenotypic variance in SCS100 attributed to the genetic component. An increase in heritability estimates for SCS as THI increased was also observed in different Holstein populations (Hammami et al.
2015; Santana et al. 2016; Bohlouli et al. 2019) and could be attributed to the increase of SCS leading to greater expression of genetic variability (Vinet et al.
2023).
The genetic correlations estimated between SCS100 and milk production traits are present in Table 5. Under thermal comfort, low genetic correlations were found between SCS100 and milk production traits, which is consistent with the correlations reported in the literature (Samoré et al. 2008; Jamrozik et al.
2010; Jamrozik and Schaeffer 2012). The genetic correlations between SCS and yields under heat stress were close to zero. This indicates that the genes involved in SCS under thermal neutral conditions are genetically distinct from those regulating heat tolerance for production.
With respect to the genetic component of SCS heat tolerance, the genetic correlations were −0.13 (±0.16), −0.55 (±0.15), −0.36 (±0.14), for milk, fat, and protein yields, respectively. The negative genetic correlation suggests that selecting more heat-tolerant cows based on milk components traits will indirectly select cows that are less affected by heat stress impacting their SCS in early lactation. Hence, breeding for heat tolerance based on milk production, in combination with efficient management practices, is a sustainable strategy to tackle heat stress, which potentially contributes to reducing subclinical mastitis.
Nevertheless, the relationship between SCS and milk production traits is complex, largely due to the mathematical and biological dependence of somatic cell concentration on milk volume. Alternative traits have been proposed to address this concern, such as the daily total SCS (Lehew and Dechow 2021). Furthermore, as the present study focused exclusively on first-parity cows, future investigations including later-parity cows will provide additional insights on the effect of heat stress on udder health.
Relative breeding value and G × E interaction
In Canada, genetic evaluations for functional traits such as SCS are expressed in RBV, which allows producers to easily interpret and compare values for multiple traits using the same scale and identify elite bulls (Lactanet 2021). Thus, for SCS100, higher RBV indicates bulls with daughters relatively more resistant to mastitis and lower SCC during the early lactation. The Spearman rank correlation between the RBV from the univariate analysis of SCS100 of all bulls with ≥ 50 daughters and for the top 10% bulls under thermal comfort and heat stress conditions (ie, 20 THI units above the threshold) was 0.99 and 0.94, respectively.
The average reliabilities of the intercept and slope RBVs were 0.65 and 0.33, respectively. These results suggest that, overall, the genetic component for SCS remained relatively consistent, indicating no substantial re-ranking of bulls under heat stress conditions. However, it is worth noting that minor re-ranking occurred within the top 10% bulls.
Moreover, the moderate positive genetic correlation (0.49) between the additive genetic component under thermal comfort and heat stress conditions (ie, intercept and slope, respectively) suggests that RBV for low SCS under thermal comfort conditions will continue to rank high under heat stress conditions.
Different traits have been defined for genetic evaluation of mastitis resistance based on SCC or SCS, and the genetic correlation between these traits and clinical mastitis in the Canadian Holstein have been found to range from 0.59 to 0.69 (Koeck et al. 2012, 2014). These traits based on SCC or SCS exhibit high genetic correlations (> 0.90), indicating they are genetically similar (Koeck et al.
2012; Narayana et al. 2018). Considering future research assessing the effect of heat stress on udder health and estimating the genetic component for heat tolerance, alternative traits should be explored.
For instance, the increasing adoption of automatic milking systems enables large-scale recording of milk electrical conductivity. Electrical conductivity is a useful indicator of intramammary infection and generally exhibits higher heritability than SCS (Wethal et al. 2020; Pedrosa et al.
2023). This trait could support the development of better phenotypes for assessing the impact of heat stress on udder health.
Conclusions
Somatic cell score is a widely used indicator of udder health and is included in genetic evaluations for mastitis resistance. In the context of rising temperatures, this study assessed the relationship between SCS and THI across Canadian regions. Whole-lactation SCS increased only marginally as THI rose, and the association was not consistent across all regions.
In early lactation (first 100 DIM), the SCS–THI relationship was stronger, suggesting that heat stress has a more pronounced effect on SCS during this critical postpartum period. The heritability estimates increased as THI increased, indicating an increase genetic variability of SCS under heat stress conditions. However, little re-ranking of bulls under heat stress conditions was observed.
Therefore, selecting bulls with high genetic merit for low SCS can be performed irrespective of the environmental conditions. Furthermore, the genetic correlations between heat tolerance components of SCS and milk production traits suggest that selection for heat tolerance based on production traits may indirectly favor low SCS under heat stress.