DOI: 10.5061/dryad.2v6wwpzpx
title: Too hot for the devil? Did climate change cause the mid-Holocene extinction of the Tasmanian devil (Sarcophilus harrisii) from mainland Australia?

description: 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.
description: <b>Macroclimate</b> 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: <i>v / v<sub>o</sub> = (h /
 h<sub>o</sub>)<sup>α</sup> </i>
 where v is the wind speed at height h (m/s),
 <i>v<sub>o</sub></i> is the wind speed at height
 <i>h<sub>o</sub></i> (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.
 <b>Microclimate and animal models</b>
 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.<b> </b>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).
 <b>Stress indices</b> 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). <b>Analysis of stress
 indices</b> 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.
description: 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.