Global risk of heat stress to cattle from climate change

On 27 February 2020, a search of the Web of Science Core Collection for terms relating to different livestock species and THI yielded 1138 results, of which 298 were included after systematic screening (see table S1 for details of the databases included, the search terms and inclusion/exclusion criteria); 29 additional studies were included from citation chaining. While the original intention was to examine heat stress in different livestock species, because most of the articles were on cattle (255 of 327), and these were also most comparable in their reporting, the scope of the study was reduced to focus only on this species. Of these 255 articles, 26 reported the results of 105 underlying studies (henceforth referred to as 'records') in sufficient detail for their data to be used in a meta-analysis determining the threshold for the onset of heat stress-related impacts, and 27 articles included 59 records that quantified the impact of heat stress on milk yield. Because there was some overlap between these two groups, in total 40 articles were included (figures S2, S3 and supplementary data). THI equations and their sources, breed, study location, management system, impact category (fertility, production or mortality), and the magnitude of impact (litres or kilograms of milk per cow per day per THI) were extracted.

Because the underlying studies used different methods to identify THI thresholds, to identify and compare onset of heat stress in a standardized way, we used the WebPlotDigitizer [51–53] to obtain the individual data points from the 105 records that plotted changes in mortality, fertility or production against THI. We then tested for a threshold for onset of heat stress impacts using the same method for each record, using piecewise regression to select the number of breakpoints and the Bayesian Information Criterion in the 'segmented' [54] package in R. As a small number of studies also included thresholds for onset of cold stress, a model with two breakpoints was required to fit these data. Of the 105 records, 5 (5%) showed no relationship of outcomes with THI, and 10 (9.5%) showed no breakpoint, instead having strictly increasing impacts from the lowest reported THI. A sensitivity analysis checking whether including the ten monotonically increasing records altered the overall thresholds showed them to have little effect (see supplementary materials; confidence intervals are all bootstrapped); therefore, we only used the 90 records with an inflection point for our analysis. Of these 90 records, 85 used local weather station information and 5 used on-farm measurements. We discuss the limitations and strengths of these two sources of weather data in the supplementary materials.

For consistency, when comparing THI thresholds across studies, we converted all THI specifications to Kelly 1971 specification [25] (equation (1)) as follows: using the THI equation from each study and fixing RH at 50%, we calculated the temperature from the THI threshold identified during our breakpoint analysis. This calculated temperature and 50% RH were then inputted into the Kelly [25] equation to get the normalized THI threshold. We used this conversion method because the alternative THI equations in our sample (e.g. [55]) are comparatively flat in humidity space (that is, humidity is weighted less than in Kelly [25]). Converting at 50% RH minimizes conversion bias.

$$\begin{align*}{\text{THI}} & = \left( {1.8 \times T + 32} \right) - \left[ \left( {0.55-0.0055{ } \times {\text{ RH}}\% } \right)\right.\notag\\ &\quad \left.\times \left( {1.8{ } \times { }T - 26} \right) \right]\end{align*}$$

equation (1) [25] (T in degrees Fahrenheit, RH in percent)

We used daily mean THI in our calculations, which has been shown to be a good indicator of heat stress, and was most commonly reported in the included studies [34, 42] (see supplementary materials). While daily maximum THI predicts the most severe heat stress at that moment, daily mean THI indirectly incorporates night-time cooling, or whether the cattle experience relief from heat stress at night [44]. For those studies that used daily maximum THI, these were corrected to daily mean THI using the baseline climate (1985–2014) for their georeferenced location from HadISD v3.1.0.2019f (HadISD is 'a sub-daily, station-based, quality-controlled dataset designed to study past extremes of temperature, pressure and humidity and allow comparisons to future projections' [56], see supplementary materials).

To test whether the onset of heat stress increases from cool to hot locations, a potential indicator of adaptation or tolerance, we compared thresholds against their current (1985–2014) climatological mean annual daily THI obtained from ERA-5 (ERA-5 is the fifth generation European Centre for Medium-Range Weather Forecasts reanalysis of global climate and weather, and provides hourly estimates of atmospheric, land and oceanic climate variables covering the Earth on a 30 km grid) [57]. For records where cattle are distributed across large regions (e.g. New Zealand or the state of Florida), we used the spatially weighted mean of the annual daily THI based on the observed cattle distribution in the region [58]. Regression models included THI thresholds as the response variable and THI of the current climate as the predictor variable, with the country of the study and impact type (production, fertility and mortality) as fixed effects. Despite finding a positive and weakly statistically significant relationship between THI thresholds and mean THI of the current climate at a site (p = 0.1; table S6 and figure S16), which provides some evidence of adaptation in hotter climates, and matches findings for other agricultural systems [59], particularly crops [60, 61], for this analysis, we elected not to account for spatially explicit thresholds when determining the future risk of heat stress and for quantifying impacts of exposure to heat. This was because: (i) the apparent effect of current climate on heat stress threshold may overestimate adaptation in hotter regions due to the limited number of studies, (ii) for projected future risk, there is no evidence of the time scale over which this adaptation may take place, and (iii) confounding with differences in breeds/types/management systems. This is a limitation that future work should attempt to address.

We also compared the THI thresholds for cattle kept indoors with those kept outside and found there was no statistically significant difference (indoor median = 69.0 (s.d. 7.7), outdoor median = 67.2 (s.d. 8.8), Moods median test = 1.65, p = 0.19).

Gridded reanalysis hourly temperature and dewpoint temperature were taken from ERA-5 [57]. Daily THI mean and maximum values were calculated for each 0.25° grid cell using equation (1). RH was calculated using the equation from Stull [62] and Bolton [63]. Heat stress duration was calculated as the mean number of days in a year above the THI threshold (THI = 68.8). Severity of heat stress was calculated as the mean amount by which daily THI exceeded the THI threshold for all days above the THI threshold in a year. A description of the results of a sensitivity analysis testing our heat stress projections for cattle with high heat tolerance is presented in the supplementary materials.

Climate model historical simulations (1985–2014) and future projections under the Shared Socioeconomic Pathways (SSPs, representing five scenarios where global society, demographics and economics change to different extents over the next century) were taken from the coupled model intercomparison project phase 6 (CMIP-6) archive [64, 65]. CMIP-6 consists of a large set of general circulation models (dynamically coupled to ocean and land models) which simulate climate under harmonized conventions including model initialization, setup and forcings over a wide variety of control simulations and future economic pathways. Future climate projections rely on illustrative scenarios that combine each SSP with one of four possible future greenhouse gas concentration trajectories known as representative concentration pathways (RCPs).

Daily mean near surface temperatures (tas) and specific humidity (huss) were used to calculate daily mean THI using equation (1).

Future heat risk was calculated for each 0.25° grid cell as the number of THI units over threshold for each day, summed across all days in a decade (sum of max(0, THI − THI threshold) for all days in a decade), and multiplied by the number of cattle in the grid cell in that decade, which returns units of cattle-THI. Percentage increases in risk were calculated as the percentage increase in heat risk in 2095–2100 compared to 1985–2014 in cattle-THI.

For details of the climate models used, the steps taken to calculate climate and heat stress projections for end of century for each SSP-RCP scenario, country-level aggregation, and geographic cattle density weightings, please refer to the supplementary materials.

Here we focused only on SSP1-RCP2.6 ('Sustainability' scenario) [66] and SSP3-RCP7.0 ('Regional Rivalry' scenario) [67]. We selected these two scenarios as they represent among largest projected divergence in future livestock consumption globally, and also in the expected challenges to climate change mitigation and adaptation. For each SSP, we used the so-called marker model scenario for land use change projections (LUH2 [68], see Popp et al [69] and Stehfest et al [70] for details about model assumptions; 'reforming economies' are not included).

Projected expansion or contraction of future cattle numbers was calculated at the regional level, based on 2015 numbers [8] and applying the projected percent change in production from the integrated assessment models (IAMs) used for SSP analyses for each decade until 2090 [71]. Cattle were added or removed from individual countries based on those countries' share of change in agricultural land (cropland + pasture + rangeland) from the same baseline IAM. We did not account for density variation within a country, or move existing cattle within a country. For further detail please see the supplementary materials.

The linearized reduction in milk yield per cow per day per unit THI over threshold was extracted from studies reporting changes in milk yield in kg or l per unit THI per day, yielding 59 published slope coefficients. Because there is evidence of differing heat tolerance between Bos taurus (taurine) and Bos indicus (indicine) cattle, both anecdotal and shown experimentally [72], to account for the possible effect of type on changes in milk yield, we used genotype data to determine the dominant type of dairy cow in each country [73] (figure S9).

The mean daily increase in THI units over the overall median threshold was calculated for SSP5-RCP8.5 and SSP1-RCP2.6 in 2045–2055, with the number and location of cattle in each country held fixed at 2010 levels. Dairy cows are assumed to be evenly spatially distributed across all cattle occurrences within a country (no dairy-only dataset is available). A sensitivity test was also applied where dairy cows were assumed to be evenly distributed across the coldest third of grid cells with cattle in each country (figure S15). To determine the reduction in milk yield in litres per cow per day, the mean daily increase in THI units over threshold at the country level was multiplied by the slope coefficients from the literature for either taurine or indicine cattle, depending on the dominant type in each country. The percent reduction in milk yield at a country-level was then calculated from the baseline daily milk yield production data from the Food and Agriculture Organization of the United Nations (FAO) statistical database (FAOSTAT) [8]. Absolute regional production values were aggregated based on country-level total milk production and the country-level percent reduction.

We separated uncertainty in projected impacts due to: (i) uncertainty in published statistical estimates of the effect of heat stress on milk yield, and, (ii) uncertainty in climate projections. To quantify (i) we used the median climate model projection for a given SSP-RCP scenario, and projected milk yield losses for the bootstrapped 95% confidence interval around the median published slope estimate. To quantify (ii) we used the median slope coefficient and reported the range of impacts for each region across the 11 global climate models.

Peer-reviewed studies quantifying heat stress impacts on cattle fertility, production and mortality are mostly concentrated in temperate climates with high cattle density, particularly Europe and North America, with cattle in India, South America and Africa relatively understudied (figure 1(a)). Nevertheless, the studies in our sample cover a wide range of climate conditions (figure 1(b)), spanning 20 years (1999–2020) and 22 countries.

Zoom In Zoom Out Reset image size

Combining all studies globally, we find a median threshold for onset of negative impacts at THI of 68.8 (95% C.I. 67.3, 70.7; figure 1(c)). We did not find any statistically significant difference between median THI thresholds for mortality (65.8; 95% C.I. 63.4, 68.1), fertility (68.4; 95% C.I. 61.8, 75.0), and production (mostly milk yield, 70.5; 95% C.I. 68.1, 72.9) (Mood's median test = 3.3, p = 0.20; figure 1(d)). We find that 83% of the THI threshold values from published studies are associated with the onset of negative impacts on cattle at lower temperature-humidity conditions than the frequently cited thresholds used in warning systems for cattle, such as the Livestock Weather Safety Index [74] (figure 1(c)).

Combining the median heat stress threshold (THI > 68.8) and the current global distribution of cattle, we estimate that 77% of cows are already exposed to climate conditions likely to cause heat stress for at least 30 days each year, with 20% of cattle (those in the tropics) exposed to heat stress conditions year-round (figure 2(a)). Holding global cattle distribution constant, by 2100, under a high emissions scenario (SSP5-RCP8.5, that is >4°C global warming above 1850-1900 levels for the set of climate models used in this study), these percentages are projected to increase to 90% of cows exposed for at least 30 days each year, and 34% experiencing year-round heat stress.

Zoom In Zoom Out Reset image size

The projected increase in heat stress duration is greatest in sub-tropical and temperate regions (figures 2(b) and (c)). By 2100, under a high emissions scenario, year-round heat stress expands into southern Brazil, southern Africa, northern India, northern Australia, and central America (figure 2(c)). Multiple months of heat stress are projected across Europe and the United States, with the emergence of over 180 days of heat stress across China and southern Japan, Mexico and the southern United States, southern Europe, and southern Australia (figures 2(c), S6). Under a low emissions scenario (SSP1-RCP2.6, that is ∼2°C global warming above 1850-1900 levels for the set of climate models used in this study), the expansion of 180 days of heat stress is limited to sub-tropical regions (figure 2(b)). Uncertainty in projections due to differences across multiple climate models indicates a large range in the potential spatial extent of regions exposed to two-month heat stress across Eurasia and North America under a high emissions scenario (figure S5).

The severity of heat stress is also projected to increase with climate change. Heat stress severity is currently greatest in the tropics, as well as northern India and Pakistan (figure 2(d)). Tropical regions that already endure year-round heat stress for cattle are projected to experience substantially more severe heat stress under climate change (figures 2(e) and (f)). Although absolute temperature increases are projected to be greater toward the poles, increases in the severity of heat stress for cattle are projected to be greatest in the tropics where THI conditions already exceed the heat stress threshold (figures 2(g) and (h)).

Reductions in greenhouse gas emissions will reduce the projected severity and duration of heat stress in the future. Holding current cattle densities and spatial distributions constant, for a high emissions scenario (SSP5-RCP8.5) in 2070–2100, the median cow is projected to experience an additional two months (63 days) of heat stress at a mean severity of 9.3 (±3.9 s.d.) units over the heat stress threshold. In contrast, a low emissions pathway (SSP1-RCP2.6) is projected to reduce this exposure to an increase of only 24 days with a mean severity of 6.8 (±3.1 s.d.). Excluding regions already enduring year-round heat stress, a low emissions scenario avoids on average 48 days (mean ± 34 s.d.) of exposure for the median cow.

Multiple regions have both high severity of heat stress under current climate, as well as high cattle densities, including: India, Brazil, the Sahel, southern United States, and China (figure 3(a)). However, although increases in temperature in these regions increase risk of heat-related impacts, changes in cattle exposure to heat due to future land-use decisions and associated changes in livestock numbers are also potentially important drivers of risk. For instance, in some regions, under a future scenario of 'Regional Rivalry' (SSP3-RCP7.0), with high population growth in developing regions, high greenhouse gas emissions, and uneven cooperation for addressing environmental concerns [71, 75], livestock production is projected to increase substantially. Assuming cattle remain a large proportion of all livestock, we project cattle numbers will nearly double in Asia, more than double in Latin America, and increase over four-fold in Africa for SSP3-RCP7.0 (figures S11 and S12).

Zoom In Zoom Out Reset image size

Yet, under this scenario of high emissions and low environmental protection, this expansion in cattle production is projected to occur in many of the same regions and countries where heat stress is also projected to increase the most, thereby exacerbating risk (figure 3(b)). If increases in cattle production follow patterns of overall agricultural expansion, our projections indicate that development pathways that increase cattle numbers in India and expand cattle production into the Congo Basin and the Peruvian and Brazilian Amazon will expose over 537 million more cattle to heat stress conditions by 2090 compared to 2010 cattle numbers (figure 3(b)), with regions of Peru at high risk as early as 2030 (figure S13). Globally, under the SSP3-RCP7.0 scenario, the number of cattle exposed to a 1 unit THI increase over threshold is projected to increase by 1.078 billion compared to 2010 cattle numbers (figure 3(b)).

This combination of rising temperatures and projected changes in the number and location of cattle drives risk across both developed and developing country regions. Under the 'Regional Rivalry' scenario, the highest increase in risk is projected for Africa and Latin America, followed by Asia and the Organisation for Economic Cooperation and Development (OECD) countries (figures 4(a)–(d)). By 2090–2100, these increases correspond to a 200% increase of heat risk to cattle in Asia and the OECD, a 300% in Latin America, and a 1000% increase in Africa compared to current conditions. However, the relative contribution to heat risk that is from warming versus changes in the number and location of cattle differ strongly among these regions. In Asia and Latin America, the projected increase in risk is due almost equally to the increasing heat hazard from warming and increased exposure from increasing cattle numbers in hot regions. In contrast, in Africa, the overwhelming majority of the projected risk increase is driven by increasing cattle numbers and projected expansion of cattle farming into more heat-stressed regions (figures 4(e)–(h)).

Zoom In Zoom Out Reset image size

Compared to a high emissions, low environmental protection scenario (SSP3-RCP7.0), a low emissions and higher environmental protection scenario (SSP1-RCP2.6) is projected to reduce heat risk for cattle by at least 50% in Asia, 63% in South America, and 84% in Africa (figure 4).

Global milk production is projected to be adversely affected by increased heat stress (figure 5). Compared with 2015 levels and holding cattle numbers fixed, global milk yields are projected to decline by around 11 million tonnes under a high emissions scenario by 2050, but are reduced by only around 6.5 million tonnes under a low emissions scenario (figure 5(b)). Fractional losses (>50%) are projected to be highest in Africa (due to low current milk yields), while absolute losses of milk yield are largest in Africa and Asia (figures 5(a) and (b)). However, climate uncertainty in these projections is high, with projected losses in Africa, Latin America and Asia potentially ten times greater than the median response (figure 5(b); table S7).

Zoom In Zoom Out Reset image size

A sensitivity test for climate impacts on milk production considered that dairy farming is sometimes more concentrated in cooler regions of a country. Placing dairy cows in the coldest third of each country's climate resulted in OECD countries seeing the greatest reduction in projected impacts, with projected losses in milk production due to future increases in heat stress decreased by almost two-thirds (although these countries were already projected to be least impacted). Projected impacts were decreased by a third in Latin America, but Asia, the Middle East and Africa showed smaller benefits for decreasing impacts (figure S15; table S7).