Overview of the modelling approach

The primary data required to build chronological analyses are derived from dates of fossilised remains (or dated human artefacts expressing the timing of first human occurrence)37,38. Because of sporadic fossilisation processes, the youngest dated fossil record in a time series is unlikely in the extreme to reflect the true timing of a species’ demise, just as much as the oldest dated fossil is unlikely to indicate the true timing of a species’ arrival (i.e., the Signor-Lipps effect)13. A whole suite of statistical models can be applied to correct for this and estimate an unbiased timing of arrival or extinction (see ref. 14), but the lack of spatial resolution from biased sampling arising from a restricted spatial distribution of fossil sites constrains the inferences of extinction timing either to specific locations or to broad-scale (e.g., continental), spatially averaged phenomena39,40,41,42.

To create continuous maps of unbiased timings of megafauna extirpation and human arrival given the inherent rarity of fossil sites, and to investigate the geographic patterns of the timing of regional extinctions, we designed and implemented a new statistical approach to infer spatial patterns of the regional timing of megafauna extirpations and initial human appearance across south-eastern Australia. Our approach explains these dynamics as a function of the spatial patterns of palaeoclimate change simulated over the last 120,000 years. We successfully validated the approach against hundreds of scenarios of simulated datasets (Supplementary Fig. 2), and the error map presents a confidence interval of the duration of the window of coexistence/non-coexistence of <3000 years for most of the study area; this indicates robust extrapolation performance of the method based on few data (Supplementary Fig. 3). We restricted our study to south-eastern Australia (Fig. 1a) owing to the scarce and patchy distribution of high-quality data available for other parts of the continent (see ‘Megafauna and human datasets’ below).

Climate variables

We used a climate reconstruction based on the three-dimensional Earth-system LOVECLIM model of intermediate complexity19,39. The model includes representations of the atmosphere, ocean and sea ice, land surface (including a vegetation model that simulates the dynamics of two main terrestrial plant functional types [trees and grasses], as well as deserts), ice sheets, and the carbon cycle. We simulated 1000-year average climates over the past 120,000 years using LOVECLIM. The original spatial resolution of LOVECLIM is 5.625 × 5.625°, but we downscaled the output resolution to 1 × 1° using bilinear interpolation because it retains the integrity and limitations of the original model output, where orography is highly smoothed relative to the real world40. For each grid cell and each temporal snapshot, we extracted mean annual temperature, mean annual precipitation, mean annual available freshwater (calculated as the difference between the mean annual evapotranspiration and mean annual precipitation), mean annual desert fraction, and mean annual net primary production. We then calculated the anomaly of each of these variables by substracting their value every 1000 years for the last 120,000 years from their respective average values calculated over the interval 50–30 ka (i.e., time window of the main megafauna extinction events in Australia7) to capture the temporal variation in climate during the period of megafauna extinctions. These five climate variables characterise both natural resources for megafauna species41 and constraints to human appearance and movement23,39.

Megafauna and human datasets

From the 2138 records of extinct fauna species in the FosSahul database42, we quality-rated each fossil date following the A* to C scale (from ‘high quality’ to ‘unreliable’) based on objective criteria43, including reliability in sample pretreatment and measurement (see Supplementary Table 1). We only used A* and A quality-rated dates (Fig. 1) 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 6349 archaeological records (ref. 44, Supplementary Table 1). 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 ka from Lake Mungo, because it is highly contested45.

Estimating the timing of extirpation and initial arrival

We applied a maximum-likelihood method to correct for the Signor-Lipps effect first developed by Solow46 that we adapted for spatial inference of both megafauna extirpation and human appearance patterns. The original (non-spatial) version of the model assumes that the true ages (e.g., estimated using radiocarbon techniques) of specimens within a time series, T 1 ,…, T n , are independent and uniformly distributed over the interval (T ext , γ), which correspond to the extirpation (T ext ) and settlement time (γ) for the taxon under investigation, respectively. In contrast to the original version of this model that assumes a constant dating error across samples, we assumed that the radiometric errors associated with each estimated age were approximately normally distributed47. Thus, \(\widehat T_k\) estimates the kth age assuming the true age T k follows:

$$\widehat T_k\sim g\left( {T_k,\sigma _{\left( k \right)}^2} \right)$$ (1)

with g being the Gaussian probability density function describing the radiometric error σ. The estimated timing of megafauna extinction \(\widehat T_{{\mathrm{ext}}}\) and the estimated timing of initial human arrival \(\widehat \gamma\) are then calculated by numerically maximising the log-likelihood ℒ (ϑ) over ϑ:

$${\cal{L}}\left( \vartheta \right) = \sum_{k = 1}^{n} {\log h\left( {\widehat T_{k}} \right)}$$ (2)

with ϑ representing either \(\widehat T_{{\mathrm{ext}}}\) or \(\widehat \gamma\) and h the probability density of \(\widehat T_k\)

$$h\left( {\widehat T_k} \right) = \int_\gamma ^{T_{{\mathrm{ext}}}} {g\left( {t,\,\sigma _k^2} \right)\frac{{dt}}{{T_{{\mathrm{ext}}} - \gamma }}}$$ (3)

To adapt this approach to questions of spatial inference, let us define W as the spatial landscape, and T ϑ(x) as either the date of megafauna extirpation or initial human arrival at a given grid cell x in W. We gridded W at a spatial resolution of 1 × 1° (i.e., the same resolution as the climate data for consistency), and then used Eq. (3) to estimate T ϑ(x) in each grid cell of W. However, to estimate T ϑ(x) in each grid cell, even for those cells without data, we used the entire dataset of W and attributed a weight for each age on T ϑ(x) as a function of the distance of the age relative to x. Thus, the closer the age to x, the more weight this age will have on T ϑ(x) . Calculating the distance of each age relative to a given grid cell x in W is similar to calculating a local density of age around x.

First, we assumed that (\(\widehat T_1\),…, \(\widehat T_n\)) are individual point estimates of the taxon’s presence (i.e., megafauna or human) based on independently discovered dated specimens found at location sites (x 1 ,…, x N ). According to Eq. (1), the estimated timing of human appearance or megafauna extirpation \(\widehat T_\vartheta \left( x \right)\) in a grid cell x is calculated by numerically maximising ℒ x (T ϑ ) over ϑ across W:

$${\cal{L}}_{x}\left( {T_{\vartheta} } \right) = \sum_{k = 1}^{n} {\log h_{\gamma} \left( {\widehat T_{k}} \right)^{w\left( {x - x_{k}} \right)}}$$ (4)

where w(x−x k ) is a weighting factor, such that Σw(x−x k ) = 1. The standard procedure in estimating local density is to select a weighting factor proportional to:

$$w\left( z \right) = \frac{1}{b} \times g_0\left( {\frac{z}{b}} \right)$$ (5)

where g 0 = the density of the standardised Gaussian distribution and b is an optimised bandwidth (see details hereafter) so that the larger the bandwidth, the more dated neighbour specimens from x over W are considered in the approximation of the distribution48,49.

Including a spatial component such as estimating local density generates a spatial bias in \(\widehat T_\vartheta \left( x \right)\)48 that depends on the size of the bandwidth b: wider b take specimens farther away from x into account in w(z), adding some bias into \(\widehat T_\vartheta \left( x \right)\), whereas narrower b take fewer specimens into account in w(z), thus increasing the variance of \(\widehat T_\vartheta \left( x \right)\). Here, the idea is to find an optimal size for b to obtain the lowest values of both bias and variance into \(\widehat T_\vartheta \left( x \right)\). We corrected for this bias by using a simulation-based spatial bias-correction procedure (Supplementary Fig. 1). This procedure estimates the bias generated by the model at each grid cell (along with its variance across space) given an optimised bandwidth size. The first step (model inference, see Supplementary Fig. 1 step 1) calculates a ‘preliminary’ spatial estimate from the model for each grid cell for a predefined bandwidth size (arbitrarily set to one tenth the maximum pairwise distance between the data at each location).

In the following steps we assumed this ‘preliminary’ estimate to be a true date of extinction or arrival (T ϑ(p) ) for each grid cell. Based on these T ϑ(p) , we generated n simulated time series (n = 100, Supplementary Fig. 1 step 2) following the same spatial pattern and characteristics (i.e., number of ages, laboratory dating error; see details in ref. 18) as the dated record. We then inferred for each of the n simulation time series a spatial estimate of timing of extinction per grid cell using the model (Supplementary Fig. 1 step 3), and we calculated the average estimate for each grid cell across the n simulation time series (Supplementary Fig. 1 step 4). By comparing the average estimate with T ϑ(p) , we calculated a mean bias (i.e., \(B_{({{\widehat{T}}_{\vartheta}} ({x}))}\)) in each grid cell x and the associated variance (i.e., \({{\sigma}^{2}}_{ {{\widehat{T}}_{\vartheta}} ({x})}\)) across space that are used to approximate the integrated mean-squared error across the landscape W (E):

$$E \simeq \sum_{x} \left( {{\sigma ^{2}}}_{{\widehat{T}}_\vartheta (x)} + {B_{(\widehat T_\vartheta(x))}}^{2} \right)$$ (6)

We then repeated steps 1 to 4 (Supplementary Fig. 1) using different bandwidth sizes (assuming that the size of the bandwidth is the same for each grid cell) until we obtained the lowest E (Supplementary Fig. 1 steps 5 and 6). In the final step, we subtracted this bias spatially from T ϑ(p) in each grid cell to provide a final and spatially unbiased estimate of extinction timing (Supplementary Fig. 1 step 7) given the optimal bandwidth.

Bearing of the spatial gradients

We calculated spatial gradients of (i) the timing of megafauna extirpation, (ii) timing of initial human arrival, and (iii) the five climate-driven variable generated by LOVECLIM from a 3 × 3 grid-cell neighbourhood using the average maximum technique modified to accommodate different cell widths at different latitudes. Following ref. 50, we calculated the average ‘north-south’ and ‘west-east’ gradients for the focal cell, excluding any missing values (usually along coastlines) using weightings of 1 and 2 for cells diagonal and adjacent, respectively, to the focal cell51. We calculated spatial gradients as the vector sum of the ‘north-south’ and ‘west-east’ gradients, with the associated vector angle giving the bearing of the gradient50.

Explanatory factor analyses

We constructed and compared generalised least-squares models to determine which predictor among the climate variables (see details in ‘Climate variables’ above) and initial human appearance best described spatial variation in the timing and bearing of megafauna extirpations. Generalised least-squares models have the advantage of accounting for spatial autocorrelation in predictor variables (i.e., model residuals are correlated and each observation does not contribute a full degree of freedom)52. Because of the large number of possible combinations of predictor variables (including their interactions), we used a phasic approach. We first computed and compared the Akaike’s information criterion weights (corrected for small sample size: wAIC c )53 of 29 and 60 generalised least-squares models assuming a Gaussian spatial autocorrelation structure, in areas with no coexistence between humans and megafauna (29 models) and areas with coexistence (60 models), respectively—these correspond to all combinations of predictor variables excluding their interactions. We took climate values from LOVECLIM19,39 for the time period in each grid cell corresponding to the estimated timing of megafauna extirpation and its confidence interval.