Megafauna extinction and first human occurrence

We used six established frequentist methods (see full model descriptions in Supplementary Methods, Supplementary Table 1 and ref. 17) to infer the timing of extinction for 28 identified megafauna species (plus 8 indeterminate species identified only to genus) grouped into 16 genera (other genera do not present enough high-quality-rated fossil ages to compute the extinction window; Table 1) and first human occurrence. We applied these models to 659 megafauna fossil records and 438 archaeological records extracted from the new FosSahul database (see author information) by first aggregating species into genera and then computing the extinction window for each genus. Each method returns an extinction window (temporal confidence interval) for each taxon; thus, we calculated a window of cross-model agreement through time (that is, for every year from 120 to 0 kyr, we calculated how many models predicted extinction for a given taxon) under the assumption that higher cross-model agreement decreases uncertainties in extinction-window estimates43. From the 2,138 records of extinct fauna species in the database, each fossil was quality rated following an A* to C scale (from ‘high quality’ to ‘unreliable’) based on objective criteria20, including reliability in sample pretreatment and measurement (see details below). We used only A* and A quality-rated data and further excluded all data for which ages were obtained from materials in depositional context below or above the fossil(s) of interest. We used the same methodology to extract data from the 6,349 archaeological records. There are few human skeletal remains in the database; thus, human presence is largely restricted to artefacts related to human activities (for example, hearth charcoal, shell middens, stone tools and rock art). We excluded an age estimate of 62±6 kyr from Lake Mungo, because it is highly contested44,45. Including this age would push back our first human appearance estimate to about 62 kyr (see Supplementary Fig. 3), but does not affect our conclusion that megafauna extinctions are not correlated with climate variation. Its inclusion would result in a longer period of human–megafauna coexistence (18.5 kyr), thus extending the timeframe by about 5 kyr for humans to extirpate megafauna species continent wide. Similarly, the youngest Genyornis ages should be treated with caution, because 14C ages on carbonate materials are sensitive to exchange with younger carbon; such ages are potentially ‘minimum’ ages and the true ages could be older than the estimated peak of extinction, but this complication does not affect our conclusions. We calibrated all 14C dates to calendar years before present using the Southern Hemisphere Calibration curve (ShCal13) and Marine Calibration curve (Marine13 for marine shells) from the OxCal radiocarbon calibration tool Version 4.2 (see https://c14.arch.ox.ac.uk).

Quality rating of fossil ages

The quality-rating scheme is based on a two-stage set of objective criteria (see full method description and justification in ref. 20). An initial assessment is made of the reliability of an age, based on the dating procedure used, followed by an evaluation of the confidence in the stratigraphic relationship (referred to as ‘association’) of the dated material to the target megafaunal taxon or archaeological event.

The ages obtained using five commonly used geochronology techniques in Sahul were rated for their quality20, with reliable ages coded as m* or m.

Radiocarbon dating. Reliable ages can be obtained by applying ultrafiltration, ninhydrin or XAD-2 protocols to individual amino acids or well-preserved bone and dentin collagen, and strong oxidation reagents (acid chlorate or dichromate) to charcoal. Reliable ages can also be obtained from cellulose isolated from wood, seeds, macrofossils and gut contents, and from the carbonate fractions of shells and corals with insignificant recrystallization, as assessed by X-ray diffraction.

Amino acid racemization dating. Reliable ages can be obtained from eggshells and otoliths of the target species, provided they have acted as chemically closed systems (i.e., no exchange in amino acids following burial of remains). Multiple analyses should be reproducible with low uncertainties and age calibration requires the availability of independent (and reliable) age estimates.

Uranium-series dating. Reliable ages can be obtained from materials that behave as chemically closed systems or as open systems when combined with modelling of uranium-migration processes. Reliable ages can be obtained from fossil teeth using coupled uranium-series/electron-spin resonance dating.

Electron-spin resonance dating. Reliable ages for tooth enamel can be obtained when combined with uranium-series modelling or if the internal dose rate due to uranium in dentine and enamel is <10% of the total dose rate.

Luminescence dating. Reliable ages can be obtained using the optically stimulated luminescence signal from individual grains and multi-grain aliquots of quartz or feldspar if the sediments were well bleached by sunlight before deposition. Reliable single-grain ages can also be obtained for partially bleached sediments using appropriate models.

The second step of the quality-rating scheme attributes a final A*/A rating to each m*/m age ranked in the first step (that is, to ages that satisfy the high-quality criteria for sample pretreatment and dating), based on the strength of association of the dated material to the target species of megafauna or the archaeological remains. Direct dating of megafaunal fossils can give A* or A ages, provided they satisfy the above criteria, as association is assured. However, indirect ages require careful checking for the strength of association. Such ages can be reliable if samples are recovered from contexts with high stratigraphic integrity and where fossils or other dated materials show no evidence of reworking (for example, burial sediments containing the articulated skeleton of a target species).

Climate variables

The EPICA ‘Dome C’ record provides an index of glacial/interglacial temperature fluctuations in Antarctica and has been used previously to explore potential drivers of megafauna extinctions in Sahul7. The EPICA initiative has provided the longest ice-core climate record yet, by drilling through 3,270-m-thick ice at a site known as ‘Dome C’ in central East Antarctica22. As the time resolution of the core increases towards the present fluctuations on short time scales are more visible in recent parts of the record and can lead to misinterpretation of climate variability if the sampling bias is not corrected, we therefore resampled at an interval of 3 kyrs (to ensure at least 5 temperature values were available to calculate mean deviation per window) over 1,000 iterations and we calculated the mean temperature deviation for each interval and iteration.

The climate of a large part of eastern and northern Sahul is sensitive to changes in ENSO activity46; thus, metrics describing variation in the strength of ENSO events (such as ENSOp) are a good proxy for shifts in precipitation and temperature, and for changes in aridity. We extracted modelled ENSOp variation directly from Tudhope et al.23, who estimated ENSOp variability from application of the Zebiak–Cane coupled ocean–atmosphere model forced only by changing orbital parameters. ENSOp is an extracted sea surface temperature signal from the Nino 3 region, for which the higher power inference means more important effects of El Niño events (that is, increase in rainfall across the east-central and eastern Pacific and drier-than-normal conditions over northern Australia, Indonesia and the Philippines).

We also included temperature and precipitation proxies based on a hindcast of the HadCM3 global circulation model focussed specifically on Sahul. This model consists of linked atmospheric, ocean and sea ice models at a spatial resolution of 2.5° latitude × 3.75° longitude resampled at a 1 × 1° resolution47. The temporal resolution of the raw data is 1 kyr slices back to 22 kyr, 2 kyr from 22 to 80 kyr and 4 kyr to 120 kyr; thus, we ran our analyses using a 4-kyr running mean (Fig. 1e, f and Supplementary Fig. 1). We tested whether both the EPICA record and ENSOp data are representative of Sahul’s palaeoclimate by calculating the spatial correlation of HadCM3 model outputs (temperature and precipitation) over the last 120 kyr for each gridcell to both EPICA and ENSOp. The strong correlation between these variables justifies the use of both EPICA and ENSOp as appropriate measures of climate variation in Sahul (Supplementary Figs 5 and 6). We also used the HadCM3 palaeoclimate simulation data to calculate vectors of climate velocity (km per year) for mean annual temperature (°C) and mean annual precipitation (mm per day) following established methods34. Climate velocity describes the rate (and direction) that an organism would need to migrate to maintain an isocline of a temperature and precipitation variable34. This proxy is often considered synonymous with the rate of climate displacement for a species, such that the higher the velocity, the less probable a species will keep pace with its shifting climatic niche and thus survive. Climate velocity is considered to be more biologically relevant than the use of climate anomalies, because it accounts for regional changes in climate and the ability of topographic heterogeneity to buffer biota against these changes34. We divided the rate of climate change through time (°C per year) by the spatial gradient in climate at that location (°C km−1).

Information-theoretic ER

We fitted two regression models to the estimated number of extinct genera as a function of climate variables (described by our six proxies; Fig. 1d–h) at the same date. The first model includes a slope parameter that assumes a linear relationship between the number of genera going extinct and the climate conditions at a time step of 1,000 years. The second model is a simpler, mean-field approach, with only an intercept that assumes no climate–extinction relationship (the ‘null’ model). The ER provides a way of explicitly evaluating the parameter bias-corrected likelihood of the null hypothesis (that is, no correlation between the number of extinct genera and climate variation) against the alternative (i.e., correlation between climate variation and the number of extinct genera). The ER is calculated as the sample size-adjusted AIC c weight of the slope model divided by AIC c weight of the competing model (in our case, the intercept-only model). Higher ERs (>3, see the ER interpretation scale in ref. 27) provide stronger support for the alternative hypothesis (model with slope) versus the null hypothesis, meaning that more genera went extinct as climate variation increased; for example, a strongly correlated climate variable might yield an ER>150. We also accounted for a possible temporal lag between climate conditions and megafauna extinctions, and for uncertainties in climate variables. We calculated the temporal lag by regressing against climate from 0 to 20,000 years (at a 2 kyr time step=11 temporal lag scenarios) for the period earlier than 35 kyr ago, which has the maximum number of extinct genera (Fig. 1a). For each of these scenarios, we generated 1,000 new sets of climate values resulting in 1,000 random resamples of the climate data within their confidence intervals (11 scenarios of ‘temporal lag’ × 1,000 random resamples of each climate-lag data set within their confidence intervals=11,000 sets of climate values for each climate proxy). We refitted both linear models to the estimated number of genera going extinct as a function of these new sets of climate values and recalculated their ERs.

Data availability: Models and data used in analyses are available on request to F.S. (frederik.saltre@adelaide.edu.au) from the Global Ecology Laboratory, School of Biological Sciences at the University of Adelaide. Sahul’s megafauna data set is also available online in the AEKOS data repository (DOI:10.4227/05/564E6209C4FE8).