The possible role of climate change in late Quaternary animal extinctions is hotly debated, yet few studies have investigated its direct effects on animal physiology to assess whether past climate changes might have had significant impacts on now-extinct species. Here we test whether climate change could have imposed physiological stress on the Tasmanian devil (Sarcophilus harrisii) during the mid-Holocene, when the species went extinct on mainland Australia. Physiological values for the devil were quantified using mechanistic niche models of energy and water requirements for thermoregulation, and soil-moisture-based indices of plant stress from drought to indirectly represent food and water availability. The spatial pervasiveness, extremity, and frequency of physiological stresses were compared between a period of known climatic and presumed demographic stability (8000-6010 BP) and the extinction period (5000-3010 BP). We found no evidence of widespread negative effects of climate on physiological parameters for the devil on the mainland during its extinction window. This leaves cultural and demographic changes in the human population or competition from the dingo (Canis dingo) as the main contending hypotheses to explain mainland loss of the devil in the mid-Holocene.

Macroclimate

For paleoclimate assessment we downloaded bias-corrected data on monthly temperature (min and max), precipitation, and relative humidity using PaleoView, a free software that generates paleoclimate data at different temporal scales at a 2.5° × 2.5° resolution (Fordham et al. 2017). We selected an area cropped to Longitude 112.5°E- 155°E, Latitude 45°S to 10°S. We took ten-year averaged monthly climate variables centred on each ten-year period for 8000-2000 BP; i.e., the years 5010 BP (averaged from 5015 BP to 5006 BP) and 5000 BP (averaged from 5005 BP to 4996 BP) were taken. The PaleoView data were generated from the TRaCE21ka experiment (Liu et al. 2009; Otto-Bliesner et al. 2014) that uses the Community Climate System Model version 3 (Otto-Bliesner et al. 2006; Yeager et al. 2006), a global coupled atmosphere-ocean-sea ice-land general circulation model that includes a dynamic global vegetation module. The microclimate model that we used also requires cloud cover and wind speed data. This monthly data at yearly intervals—also generated by the TRaCE21ka experiment—was downloaded from the National Center for Atmospheric Research website (see Data Availability) and resampled from 3.75° × 3.75° to 2.5° × 2.5°, using the raster package (Hijmans 2019), then averaged over ten-year time periods to be consistent with the extracted PaleoView data. The cloud and wind data consist of 26 pressure levels representing slices of the atmosphere. We selected the lowest altitude level (~60 m) for wind speed, corrected to a reference height of 1.2 m by the equation:

v / v_o = (h / h_o)^α 

where v is the wind speed at height h (m/s), v_o is the wind speed at height h_o (m/s), and α is the wind shear exponent (0.15 to represent open grassland) (Campbell & Norman 1998). The mean value of mid-atmosphere cloud fraction (pressure levels 18-25 corresponding to 238- 6173 m altitude) was used for cloud cover and a diurnal pattern was imposed by multiplying by three to obtain daily maxima and by 0.5 such that cloud cover had an asymptotic relationship with rainfall at ~100% cloud cover. Diurnal variation in humidity was imposed by obtaining the daily minimum and maximum relative humidity from the mean daily humidity and air temperature via the WETAIR function in NicheMapR, on the assumption that the absolute relative humidity remained constant through the day.

Microclimate and animal models

Mechanistic niche models consist of two sub-models: a microclimate model and an animal model (Porter & Kearney 2009; Kearney et al. 2021a) (see Figure 1). We used the microclimate and endotherm model of the NicheMapR package (Kearney & Porter 2017; Kearney et al. 2021a) for the R programming environment (R Core Team 2018). The microclimate model downscales daily macroclimate data into hourly environmental conditions (air temperature, wind speed, relative humidity) at the height of the animal of interest as well as computing long and short-wave radiation fluxes and heat and water budgets for the substrate. The microclimatic variables are then used as the input by the animal (endotherm) model which solves the heat and mass balance equations for the animal given its functional physiological, morphological and behavioural traits (for example see (Kearney et al. 2021b)) .

We modified the micro_global function from NicheMapR to prepare the PaleoView/ TRaCE21ka data for input to the NicheMapR microclimate model. We interpolated (periodic spline) the monthly minimum and maximum air temperature, relative humidity and cloud cover for the decadal-averaged palaeoclimate data across 365 days and then ran the model for 730 days to allow sufficient spin-up time for the soil moisture calculations to reach steady state. For precipitation, we allocated monthly rainfall across the days of the month assuming the present-day pattern of rainy days per month (from New et al. (2002)), with the assumption that 50% of the rain fell on the first day of each month. This provided realistic annual cycles of soil moisture. Bulk density was set at a standard value for soil(1.3 g/cm³)(Campbell & Norman 1998) and hydraulic properties of a loam were assumed for soil moisture calculations. We simulated local air temperature, wind speed and relative humidity at a height of 30 cm – the approximate mid-point of a standing devil.

The environmentally imposed heat stress and associated energy and water requirements of the devil was simulated using the default version of the endoR function (Kearney et al. 2021a) which solves heat budgets for endotherms given their functional traits according to a specific sequence of morphological, behavioural and physiological thermoregulatory responses. The required parameters include the target core body temperature, basal metabolic rate, presence of fat, fur properties, and the size and shape of the species. These were obtained from the literature and the predictions compared against laboratory data on devil metabolic rates at different air temperatures (see Supplementary material for details). The resultant model is best considered as a general model of stress for an endothermic species of the devil’s size and shape because it does not include species-specific diet or behaviour besides denning behaviour. Values on the lower end of the weight range (6.5 kg) were chosen because stress on females, which are smaller than males, would likely have a greater negative effect on reproduction and therefore population growth and viability. The endoR model applies an ordered sequence of changes in behaviour (posture change) and physiology (piloerection, change flesh conductivity, allow core temperature to rise, pant, sweat) to maintain its specified core temperature given the minimum permissible metabolic heat production. Under cold conditions the model finds the metabolic rate required to maintain body temperature; under hot conditions the model finds the required water loss rate, contingent on the thermoregulatory options. The model thereby quantifies environmental stress in terms of energy and water requirements contingent on thermoregulatory responses. If the endotherm model cannot find a solution given the stated parameters, the model “fails”. This does not imply the conditions would have been certainly fatal, but rather that survival would be unlikely given the limited behaviours programmed in the model.

Microclimate inputs for our endotherm simulations made a distinction between foraging and denning environments. As devils are nocturnal and rest in burrows during the day (Rose et al. 2017), for night-time conditions (defined by the microclimate model predictions of solar zenith angle) the devils were assumed to experience conditions in the open at 30cm above the ground, and for day-time conditions we assumed that the den was a burrow at 50cm below ground with temperature-adjusted humidity and a wind speed of 0.1 m/s. Different sleeping conditions—a humid burrow (90% relative humidity and 0.1 m/s wind speed) and above ground at 30cm under the cover of thick vegetation (90% shade and half the wind speed at 30cm to simulate obstruction from the vegetation)—were also tested (see Supplementary material).

Stress indices

We formulated three stress metrics—plant, energy, and water stress—using the outputs from the microclimate model and animal model. Plant stress, as an indirect index of food availability for devils, was assumed to occur if the computed soil water potential at 10 cm below ground fell below -1500 J/kg in a given hour (a standard value for the permanent wilting point of plants (Campbell & Norman 1998)). Energy stress was defined as computed energy requirements being above the basal energy requirements multiplied by a conservative activity factor of 2x. Previous studies on mammals have used activity multipliers of 2.3 (Wang et al. 2018) and 2.25-4.5 (Mathewson et al. 2020). We chose a conservative value to compensate for our devil-specific behaviour dictating that devils remain active the entire night when they rest periodically during the night (Andersen et al. 2020). Water stress was assumed to occur when the panting multiplier exceeded one, indicating that the animal would need to pant (and thus lose water) to rid itself of excess heat. Hours when water stresses coincided with plant stress were also calculated, indicating a scenario of drought-like conditions with low food supply and high-water requirements. We calculated two composite stress metrics: physiological stress, the combined number of energy and water stress hours; and gross stress, physiological and plant stress combined. We quantified these six stress metrics in three ways resulting in 18 indices: the total number of stressful hours to occur in a year (total stress hours), the maximum number of consecutive hours when stress occurred in a year (cumulative stress hours), and the median number of stressful hours to occur in a year (frequent stress hours).

Analysis of stress indices

The stress indices were analysed to test for an increase in the interval preceding and including the devils’ mainland extinction. First, we obtained the ratio of summed stress in the extinction period (defined as 5000-3010 BP) compared to the thermal optimum (8000-6010 BP), expressed as a percentage (Figure 2) (e.g., 150% indicates an increase of stress by 50% during the extinction period). We then tested whether the changes in stress that occurred from 7000-2010 BP had a large enough effect size (paired difference in mean) to distinguish them from the baseline (8000-7010 BP) (Figure 3). We also present the bootstrap confidence interval (95%) of these effect sizes. The bootstrapping and effect sizes were calculate using the R package dabestr (Ho et al. 2019). Deserts and xeric shrublands were excluded from this analysis as it is unlikely that devils inhabited these areas (see Supplementary material for area defined as desert and analysis including desert, Fig. S12). Finally, we mapped the year that had the largest value for each stress index (Figure 4). If the largest value tied, the earliest year was chosen.

This is the R code to run all the analyses in the manuscript. However, you will need to download the freely available paleoclimate data (from PaleoView and the National Center for Atmospheric Research) and access to a supercomputer to run all the analyses but the final output dataset is provided and the figures in the manuscript can be recreated.