# Climate-human interaction associated with southeast Australian megafauna extinction patterns

Frédérik Saltré, Joel J. Chadoeuf, Katharina Peters, Matthew C Mcdowell, Tobias Friedrich, Axel Timmermann, Sean Ulm, Corey Bradshaw

► To cite this version:

Frédérik Saltré, Joel J. Chadoeuf, Katharina Peters, Matthew C Mcdowell, Tobias Friedrich, et al. Climate-human interaction associated with southeast Australian megafauna extinction patterns. Nature Communications, 2019, 10 (1), 10.1038/s41467-019-13277-0. hal-02877577

HAL Id: hal-02877577

https://hal.inrae.fr/hal-02877577

Submitted on 22 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Climate-human interaction associated with southeast Australian megafauna extinction patterns

Frédérik Saltré ID 1*, Joël Chadoeuf², Katharina J. Peters ID <fnref n="1" />, Matthew C. McDowell ID ³, Tobias Friedrich ID ⁴, Axel Timmermann ID ⁵,⁶, Sean Ulm ID ⁷ & Corey J.A. Bradshaw ID <fnref n="1" />

The mechanisms leading to megafauna (>44 kg) extinctions in Late Pleistocene (126,000–12,000 years ago) Australia are highly contested because standard chronological analyses rely on scarce data of varying quality and ignore spatial complexity. Relevant archaeological and palaeontological records are most often also biased by differential preservation resulting in under-representated older events. Chronological analyses have attributed megafaunal extinctions to climate change, humans, or a combination of the two, but rarely consider spatial variation in extinction patterns, initial human appearance trajectories, and palaeoclimate change together. Here we develop a statistical approach to infer spatio-temporal trajectories of megafauna extirpations (local extinctions) and initial human appearance in south-eastern Australia. We identify a combined climate-human effect on regional extirpation patterns suggesting that small, mobile Aboriginal populations potentially needed access to drinkable water to survive arid ecosystems, but were simultaneously constrained by climate-dependent net landscape primary productivity. Thus, the co-drivers of megafauna extirpations were themselves constrained by the spatial distribution of climate-dependent water sources.

<footnote n="1">¹Global Ecology, College of Science and Engineering and ARC Centre of Excellence for Australian Biodiversity and Heritage, Flinders University, GPO Box 2100Adelaide, SA 5001, Australia. ²UR 1052, French National Institute for Agricultural Research (INRA), Montfavet, France. ³Dynamics of Eco-Evolutionary Pattern and ARC Centre of Excellence for Australian Biodiversity and Heritage, University of Tasmania, Tasmania 7001, Australia. ⁴University of Hawai'i at Mānoa, Honolulu, HI, USA. ⁵Center for Climate Physics, Institute for Basic Science, Busan 46241, Korea. ⁶Pusan National University, Busan 46241, Korea. ⁷ARC Centre of Excellence for Australian Biodiversity and Heritage, College of Arts, Society and Education, James Cook University, PO Box 6811 Cairns, QLD 4870, Australia. *email: frederik.saltre@flinders.edu.au</footnote>
The Late Pleistocene megafauna (usually defined as animals >44 kg) extinction wave was a rapid, large-scale ecological phenomenon of unprecedented magnitude in the Quaternary¹. However, the causes and consequences of these extinctions are still unresolved². Many conceptual models have been developed to explain the causes of megafauna extinctions, which can be classified into the following three themes: (i) climatic changes increasingly restricted suitable habitats for megafauna species³; (ii) the spread of humans—a new and efficient top predator—from Africa into other continents negatively affected megafauna either by exploiting species' naivety, or profoundly modifying their ecosystems⁴,⁵; or (iii) a possible combination of human exploitation of populations already compromised by climate-driven environmental changes (or vice versa)⁶.

Most scientific contributions to the relative weight of these three themes rely exclusively on chronological analyses and ignore spatial patterns entirely—the general approach is to compare the timing of both megafauna extinctions and the initial arrival of humans (associated with the age of the last and first palaeontological and archaeological records, respectively) to the reconstruction of contemporaneous, relative climate time series¹,³,⁷⁻¹⁰. Such analyses are challenging because robust data are rare, and inferences are often biased because older events are under-represented due to differential preservation of evidence over time¹¹, and both the oldest and youngest ages of a palaeo time series are inaccurate reflections of an event¹². This phenomenon is exacerbated in Australia where the average density of dated archaeological sites is <1 per 4000 km² and only 1 per 10000 km² in lesser-studied regions¹³. Statistical methods that attempt to overcome these chronological biases (reviewed in ref. ¹⁴) are either not spatially explicit and disregard spatial variation in extinction patterns, initial human appearance trajectories, and palaeoclimate change⁷,¹⁵, or they generate new spatial biases via arbitrary geographic binning⁸⁻¹⁰ or spatially continuous estimation¹⁶ that do not account for the uncertainty arising from sampling and taphonomic biases¹⁷, inherent dating errors¹⁸, or spatial biases generated when interpolating a linear chronology from unevenly spaced age estimates.

We overcome these methodological limitations by developing and applying a statistical approach to the available data describing megafauna extirpations and patterns of human appearance in south-eastern Australia over the last 120,000 years to test three explicit hypotheses: (1) megafauna extirpations followed changing climate conditions, (2) megafauna extirpations followed the arrival of Aboriginal populations in the region, or (3) megafauna extirpations followed a combination of both climate change and human arrival. To test these hypotheses, we generate maps of continuous regional timings of megafauna extirpations and initial human appearance, and apply generalised least-squares models to explore their relative roles in explaining both the regional timing and directions of the spatial gradients (bearings) of megafauna extirpations across south-eastern Australia. We also account for possible temporal lags, because extinctions could have been driven by earlier climatic variation⁷ (see Methods). We show that (i) >80% of south-eastern Australia had a period of human-megafauna coexistence lasting from 1000 to >15,000 years, and (ii) the pattern of megafauna extirpation in these areas is best explained by an additive effect of the patterns of human spread and freshwater availability across the region. These findings suggest that only by adding spatial complexity into standard chronological analyses can we explicitly identify both humans and climate change as the most likely drivers of these extinctions.

## Results

**Regional patterns of megafauna extirpation and human spread.** We first identified the regional locations where and

periods during which humans had putatively coexisted with megafauna species within south-eastern Australia by comparing their respective regional chronologies (i.e., extirpation and initial arrival) and their spatial bearings. To estimate an accurate regional chronology for these patterns, we applied a spatio-temporal statistical method (see details in Methods and Supplementary Table 1) to high-quality (i.e., suitable dates; see quality-rating approach in Methods) megafauna fossil records and to the archaeology of initial human appearance of Sahul (Fig. 1a, Supplementary Table 1 and Supplementary Data 1). Our method returns for each $1 \times 1°$ terrestrial grid cell in south-eastern Australia a mean estimate of both a timing of megafauna extirpation and initial human arrival along with their respective standard deviations (see details in Methods and Supplementary Note 1).

The regional pattern of megafauna extirpations started around 48 ka (1 ka = 1000 years ago) in the southeast coastal plains (location 7; Fig. 1b) and the southern coast of mainland Australia (location 9). Megafauna extirpations mostly radiated (i) 0–60° north from location 7 (Fig. 1c) to reach locations 5 and 8 between ~46 and 42 ka, via location 9 (Fig. 1a), (ii) 110 to 180° south from locations 7 and 9 to reach location 10 (Tasmania) between ~46 ka (Fig. 1b), and continued both (iii) 210 to 240° south from location 1 and (iv) 330 to 0 north from location 6 (Kangaroo Island) to reach locations 2 (Lake Eyre) and 4 by ~42 and 40 ka. In contrast, the pattern of the peopling of the region started >50 ka from location 5 (Fig. 1d), then mostly radiated northward (310° to the northwest; Fig. 1e) to location 2 by ~44 ka via location 4 (Port Augusta), and location 1 by ~41 ka (Fig. 1d) via location 3 (Darling Downs, ~44 ka). Humans reached location 10 (Tasmania) by ~42 ka from the northwest via location 7 (~44 ka) and later (~42 ka) from the northeast via the southern coast of mainland Australia at location 9.

**Drivers of coexistence/non-coexistence areas.** The regional chronology analysis indicated that 81.1% (ranging from 69 to 84%, depending on the criteria used; see Supplementary Note 2 for details) of the $1 \times 1°$ grid cells (Fig. 2a) had a period of human-megafauna coexistence lasting between 1000 (east of location 3, and between locations 6 and 8; Fig. 2b) and 8000 years (between locations 2 and 5; Fig. 2a). In these areas of temporal overlap, we built 60 generalised least-squares models to explain variation in the timing of megafauna extirpations ($Ext_t$) (Table 1) as a function of all possible combinations of mean annual temperature anomaly, mean annual precipitation anomaly, mean annual water availability anomaly, mean annual net primary production anomaly, mean annual percentage of desert anomaly, and the timing of initial human arrival (see model description in Methods). We obtained the corresponding climate data from a transient LOVECLIM¹⁹ Earth-system model experiment, which uses time-varying orbital forcings, ice-sheet volume, and greenhouse gases as a forcing covering the past 784,000 years. We calculated the anomalies relative to the period from 50 to 30 ka (see details in Methods).

None of the generalised least-squares models we tested explained any variation in the timing alone of megafauna extirpation in areas of human-megafauna coexistence (i.e., ~0% of variance explained). However, the top-ranked model explained variation in the spatial pattern of timing of extirpation (bearing = $Ext_b$; Table 1), such as the bearing of the timing of human arrival, the bearing of mean annual precipitation, and water availability (>81.2% of variance explained; Table 1). Both the bearing of the timing of human arrival and the bearing of water availability had the strongest effects on the change in estimated likelihood for this model (Fig. 3a). This indicates that these
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Spatial pattern of timing of megafauna extirpation and first human arrival. **a** Regional estimates in south-eastern Australia of the **b** timing and **c** bearing of megafauna extirpation and **d** timing and **e** bearing of initial human arrival, from high quality-rated ages of megafauna fossils (blue circles) and archaeological specimens (red stars) calculated at a spatial resolution of 1 × 1°. Grey arrows represent the spatial gradient of megafauna extirpation and human appearance from the oldest to the youngest point estimates. The angle-histogram polar plots in **c**, **e** show the distribution of bearings (calculated in degrees) of the spatial gradient of timings of **c** megafauna extirpation and **e** initial human arrival. These angle-histogram polar plots are normalised so that the height of each bar represents the number of grid cells in the bin (specific bearing) divided by the total number of observations multiplied by the width of the bin. The area of each bar is the number of grid cells, and the sum of the bar areas = 1. The spatial component is not included in these plots, i.e., the centre of these angle-histogram polar plots do not relate to any particular geographic location on the map. Each region is named as follows: Sandy Hollow Creek (1), Lake Eyre (2), Eastern Great Dividing Ranges (3), Port Augusta (4), Willandra Lakes (5), Kangaroo Island (6), southeast coastal plains (7), Sydney Basin (8), southern coast of mainland Australia (9), and Tasmania (10).</figcaption></figure>

variables are the main contributors to the goodness of fit (i.e., how well this model fits the data of timing of megafauna extirpations), even accounting for up to 5000 years of possible temporal lags in climate variation.

Regions such as locations 7, 9 and 10 (19% of the grid cells ranging from 16 to 31%; Fig. 2a, and Supplementary Figs. 4 and 5) indicated no period of human-megafauna coexistence because the first humans arrived there between 5000 (e.g., locations 5 and
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Spatial comparison of estimated timing of megafauna extirpation and initial human arrival at a spatial resolution of 1 × 1°. **a** Areas of coexistence (brown grid cells) and non-coexistence (blue grid cells) between humans and megafauna. For each grid cell, we indicated whether the timing of human arrival preceded (i.e., coexistence = humans arrived before megafauna were extirpated) or followed (i.e., no coexistence = humans arrived after megafauna were extirpated) megafauna extirpation. We accounted for confidence intervals around estimated timings so that 'coexistence' indicates when the lower confidence limit of the timing of initial human arrival is older than the lower confidence interval limit of megafauna extirpation. We discuss the sensitivity of the size of these areas to the approach we used to take into account these confidence intervals in Supplementary Note 1. **b** Duration of the window of coexistence/non-coexistence between humans and megafauna. For each grid cell, we subtracted from the estimated timing of extirpation the timing of initial human arrival. Negative values (from light yellow to dark green) indicate how long humans coexisted with megafauna (i.e., the darker the green, the longer they coexisted), whereas positive values (from light yellow to dark red) indicate how long megafauna had been gone prior to the arrival of the first humans (i.e., the darker the red, the longer megafauna had been absent).</figcaption></figure>

<table><caption>**Table 1** Performance values for the top-ranked generalised least-squares models.</caption><thead><tr><th>Model</th><th>Area</th><th>Lag(ka)</th><th>k</th><th>DE(%)</th><th>AIC<sub>c</sub></th><th>wAIC<sub>c</sub></th><th>LogLik</th></tr></thead><tbody><tr><td>Ext<sub>t</sub> ~ P<sub>a</sub> + EminP<sub>a</sub> + DF<sub>a</sub> + (P<sub>a</sub> × EminP<sub>a</sub>) + (DF<sub>a</sub> × EminP<sub>a</sub>)</td><td>n/coexist</td><td>0</td><td>8</td><td>99.4</td><td>634.4</td><td>0.6</td><td>-306.7</td></tr><tr><td>Ext<sub>t</sub> ~ P<sub>a</sub> + EminP<sub>a</sub> + DF<sub>a</sub></td><td>n/coexist</td><td>1</td><td>6</td><td>99.9</td><td>645.8</td><td>0.4</td><td>-315.5</td></tr><tr><td>Ext<sub>t</sub> ~ P<sub>a</sub> + EminP<sub>a</sub> + DF<sub>a</sub></td><td>n/coexist</td><td>2</td><td>6</td><td>48.8</td><td>643.7</td><td>0.3</td><td>-314.5</td></tr><tr><td>Ext<sub>t</sub> ~ P<sub>a</sub> + EminP<sub>a</sub> + DF<sub>a</sub></td><td>n/coexist</td><td>3</td><td>6</td><td>10.1</td><td>644.8</td><td>0</td><td>-315.1</td></tr><tr><td>Ext<sub>t</sub> ~ none</td><td>n/coexist</td><td>4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ext<sub>t</sub> ~ none</td><td>n/coexist</td><td>5</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ext<sub>t</sub> ~ none</td><td>coexist</td><td>0,1,2,3,4,5</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ext<sub>b</sub> ~ NPP<sub>b</sub></td><td>n/coexist</td><td>0</td><td>4</td><td>98.8</td><td>444.9</td><td>0.3</td><td>-217.8</td></tr><tr><td>Ext<sub>b</sub> ~ T<sub>b</sub> + P<sub>b</sub></td><td>n/coexist</td><td>1</td><td>6</td><td>74.7</td><td>1875.9</td><td>0.6</td><td>-931.7</td></tr><tr><td>Ext<sub>b</sub> ~ T<sub>b</sub> + NPP<sub>b</sub> + (T<sub>b</sub> × P<sub>b</sub>)</td><td>n/coexist</td><td>2</td><td>7</td><td>73.1</td><td>1865.5</td><td>0.9</td><td>-925.4</td></tr><tr><td>Ext<sub>b</sub> ~ T<sub>b</sub> + NPP<sub>b</sub> + (T<sub>b</sub> × P<sub>b</sub>)</td><td>n/coexist</td><td>3</td><td>7</td><td>70.4</td><td>1866.9</td><td>0.8</td><td>-926.1</td></tr><tr><td>Ext<sub>b</sub> ~ T<sub>b</sub> + NPP<sub>b</sub></td><td>n/coexist</td><td>4</td><td>5</td><td>66.6</td><td>1873.1</td><td>0.6</td><td>-931.3</td></tr><tr><td>Ext<sub>b</sub> ~ T<sub>b</sub> + NPP<sub>b</sub></td><td>n/coexist</td><td>5</td><td>5</td><td>68.9</td><td>1876.9</td><td>0.9</td><td>-933.3</td></tr><tr><td>Ext<sub>b</sub> ~ H<sub>b</sub> + P<sub>b</sub> + EminP<sub>b</sub></td><td>coexist</td><td>0</td><td>6</td><td>81.2</td><td>1886.9</td><td>0.4</td><td>-937.2</td></tr><tr><td>Ext<sub>b</sub> ~ H<sub>b</sub> + NPP<sub>b</sub></td><td>coexist</td><td>1</td><td>5</td><td>99.9</td><td>422.7</td><td>0.8</td><td>-205.4</td></tr><tr><td>Ext<sub>b</sub> ~ H<sub>b</sub> + NPP<sub>b</sub></td><td>coexist</td><td>2</td><td>5</td><td>99.9</td><td>422.9</td><td>0.7</td><td>-205.5</td></tr><tr><td>Ext<sub>b</sub> ~ H<sub>b</sub> + NPP<sub>b</sub></td><td>coexist</td><td>3</td><td>5</td><td>99.9</td><td>423.1</td><td>0.7</td><td>-205.6</td></tr><tr><td>Ext<sub>b</sub> ~ H<sub>b</sub> + NPP<sub>b</sub></td><td>coexist</td><td>4</td><td>5</td><td>99.9</td><td>424.2</td><td>0.7</td><td>-206.2</td></tr><tr><td>Ext<sub>b</sub> ~ H<sub>b</sub> + NPP<sub>b</sub></td><td>coexist</td><td>5</td><td>5</td><td>99.9</td><td>424.6</td><td>0.7</td><td>-206.4</td></tr></tbody></table>

These top-ranked generalised least-squares models contain climate (T = mean annual temperature, P = mean annual precipitation, EminP = mean annual freshwater availability, NPP = mean annual net primary production and DF = fraction of desert within the grid cell), and human predictors (H) to describe (i) the timing of megafauna extirpation (Ext<sub>t</sub>) in human-megafauna non-coexistence areas (n/coexist) and (ii) in human-megafauna coexistence areas (coexist), (iii) the bearing of timing of megafauna extirpation (Ext<sub>b</sub>) in human-megafauna non-coexistence areas (n/coexist) and (iv) in areas with coexistence (coexist). For each of these four scenarios we included five temporal lags (Lag) between the climate from 0 to 5 ka (at a 1 ka-year time step, with ka = 1000 years) for the period earlier than the estimated timing of megafauna extirpation in each grid cell. Predictor variables subscripted a (P<sub>a</sub>, EminP<sub>a</sub>, DF<sub>a</sub>) indicate that we used the mean annual anomaly relative to the period 50-30 ka for these variables, whereas predictor variables subscripted b (T<sub>b</sub>, P<sub>b</sub>, NPP<sub>b</sub>, EminP<sub>b</sub>, H<sub>b</sub>) indicate that we used the directional bearing. Shown are the number of parameters (k), metric of the model's structural goodness of fit (%DE), minimised negative log-likelihood (LogLik), weight scaled to a sum of 1 (wAIC<sub>c</sub>) and the Akaike's information criterion corrected for small sample sizes (AIC<sub>c</sub>) for the highest-ranking model. See details of the full list of generalised least-squares models describing all the combinations among variable predictors in Supplementary Data 2
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Most likely drivers of megafauna extirpation timings and patterns. Relative importance of predictor variables for the top-ranked generalised least-squares model assuming a Gaussian spatial autocorrelation structure best describing a the spatial gradient (bearings) of megafauna extirpation (Ext<sub>b</sub>) in human-megafauna coexistence areas, b the timing of megafauna extirpation (Ext<sub>t</sub>) in human-megafauna non-coexistence areas and c the bearings of megafauna extirpation (Ext<sub>b</sub>) in human-megafauna non-coexistence areas. Human-megafauna coexistence and non-coexistence areas are described in Fig. 2a and in the top panel of each barplot in Fig. 3 along with red arrows indicating the bearing of megafauna extinction gradients. The relative importance of each predictor variable is calculated as the change in the full model likelihood when one of its predictor variables is removed. Climate predictors were from LOVECLIM<sup>19,39</sup> for the time period in each grid cell corresponding to the estimated timing of megafauna extirpation and its confidence interval. Predictor variables of the model describing Ext<sub>t</sub> are mean annual temperature anomaly (T), mean annual precipitation anomaly (P), mean annual freshwater availability anomaly (EminP), mean annual net primary production anomaly (NPP), and the fraction of desert anomaly within the grid cell (DF). Climate anomalies are calculated relative to the 50–30 ka mean time period (see Methods). Predictor variables subscripted b (T<sub>b</sub>, EminP<sub>b</sub>) indicate that we used the directional bearing of these climate variables, including the directional vectors for the timing of initial human arrival (H<sub>b</sub>), to build the model describing the spatial pattern bearing of megafauna extirpation (Ext<sub>b</sub>). For clarity, we did not present the results describing Ext<sub>t</sub> for the areas with human and megafauna coexistence because we did not have any relevant model (i.e., percentage of variance explained by the models ~ 0%) in those areas. For each variable, error bars represent the standard deviation of the relative importance of predictor variables for the top-ranked generalised least-squares model accounting for the temporal lag by regressing extirpation against climate from 0 to 5000 years (at a 1000-year time step = 5 temporal-lag scenarios, i.e., 6 values of relative importance in total) per grid cell. Source data are provided as a Source Data file.</figcaption></figure>

9) and >9000 years (location 10; Fig. 2b) after megafauna species went locally extinct. A generalised least-squares model built on mean annual precipitation anomaly (P), mean annual water availability anomaly (EminP), the percentage of desert fraction (DF), and their respective interactions (i.e., P × EminP and EminP × DF) as predictive variables provided the highest support to describe the timing of extirpation in these areas (>99.4% of the variance explained and the lowest AIC<sub>c</sub>; Table 1). These three predictive variables provide the highest contribution in this model (Fig. 3b), even accounting for up to 5000 years of possible temporal lags in climate variation. The coefficients of the highest-ranked generalised least-squares model (Table 1) indicate that the oldest extinction events are associated with both lowest mean

annual precipitation anomalies, mean annual desert anomalies (i.e., both with negative coefficients; Table 2), and mean annual freshwater availability anomalies (i.e., positive EminP means higher evapotranspiration relative to precipitation influx, which means less available freshwater; Table 2). However, the bearing of net primary productivity (NPP<sub>b</sub>) provides the highest support to the spatial pattern of timing of extirpation (bearing = Ext<sub>b</sub>; Table 1 and Fig. 3c) in these areas, while the bearing of mean annual temperature became an important explanatory variable when accounting for the temporal lag in climate. The negative NPP<sub>b</sub> coefficient (Table 2) means that the regional pattern of megafauna extinction (Ext<sub>b</sub>) and net primary production (NPP<sub>b</sub>) shifted in an opposite direction.
<table><caption>Table 2 Summary statistics for the top-ranked generalised least-squares models.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Area</th><th rowspan="2">Variable</th><th rowspan="2">Coefficient value</th><th rowspan="2">std error</th><th colspan="5">Correlation</th></tr><tr><th>P<sub>a</sub></th><th>EminP<sub>a</sub></th><th>Df<sub>a</sub></th><th>P<sub>a</sub> × EminP<sub>a</sub></th><th>Df<sub>a</sub> × EminP<sub>a</sub></th><th>NPP<sub>b</sub></th><th>H<sub>b</sub></th><th>P<sub>b</sub></th><th>EminP<sub>b</sub></th></tr></thead><tbody><tr><td>Ext<sub>t</sub> - P<sub>a</sub> + EminP<sub>a</sub> + DF<sub>a</sub> + (P<sub>a</sub> × EminP<sub>a</sub>) + (DF<sub>a</sub> × EminP<sub>a</sub>)</td><td>n/coexist</td><td>P<sub>a</sub></td><td>-4098.61</td><td>555.92</td><td>1</td><td>0.01</td><td>0.63</td><td>-0.66</td><td>-0.28</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td>EminP<sub>a</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>DF<sub>a</sub></td><td>378.48</td><td>784.24</td><td>1</td><td></td><td>0.29</td><td>0.42</td><td>-0.18</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td>P<sub>a</sub> × EminP<sub>a</sub></td><td>-513.70</td><td>1044.20</td><td></td><td>1</td><td></td><td>-0.39</td><td>-0.69</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td>DF<sub>a</sub> × EminP<sub>a</sub></td><td>3592.36</td><td>1340.64</td><td></td><td></td><td>1</td><td></td><td>-0.29</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td>NPP<sub>b</sub></td><td>-10519.01</td><td>3238.69</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ext<sub>t</sub> - NPP<sub>b</sub></td><td>n/coexist</td><td>NPP<sub>b</sub></td><td>-0.38</td><td>0.09</td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Ext<sub>t</sub> - H<sub>b</sub> + P<sub>a</sub> + EminP<sub>b</sub></td><td>coexist</td><td>H<sub>b</sub></td><td>0.34</td><td>0.09</td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.04</td><td>0.33</td><td>-</td></tr><tr><td></td><td></td><td>P<sub>b</sub></td><td>-0.18</td><td>0.09</td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>-0.18</td><td>1</td></tr><tr><td></td><td></td><td>EminP<sub>b</sub></td><td>-0.13</td><td>0.07</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Coefficients estimated for each predictor variable of the top-ranked generalised least-squares model, standard error, and correlation between predictor variables to describe the timing of megafauna extirpation (Ext<sub>t</sub>) and its bearing (Ext<sub>b</sub>) in human-megafauna non-coexistence (n/coexist) and coexistence (coexist) areas. Climate variables are mean annual temperature (T), mean annual precipitation (P), mean annual freshwater availability (EminP), mean annual net primary production (NPP), and the fraction of desert within the grid cell (DF), and the human predictor is the estimated timing of initial human arrival (H). Predictor variables subscripted a (P<sub>a</sub>, EminP<sub>a</sub>, Df<sub>a</sub>) indicate that we used the mean annual anomaly relative to the period 50–30 ka for these variables, whereas predictor variables subscripted b (P<sub>b</sub>, NPP<sub>b</sub>, EminP<sub>b</sub>, H<sub>b</sub>) indicate that we used the directional bearing. For clarity, we did not present the results of the model best describing Ext<sub>t</sub> in the areas with human and megafauna coexistence because we did not have any model with a percentage of variance explained >0% in those areas.

## Discussion

Our results support an additive effect of available freshwater and human arrival on megafauna extirpations in south-eastern Australia, at least in terms of the complex regional pattern (i.e., bearing) of change (Table 1 and Fig. 3a). The combination of climate change and human pressure has been suggested in global and continental analyses as the most plausible set of conditions explaining megafauna extinctions<sup>8,9,15,20</sup>. Here, the loss of suitable habitat driven by rapid climate shifts (i.e., a shift in temperature and precipitation) would have been important precursor conditions exacerbating the effects of humans on megafauna. Similar results are also supported by more regional and local palaeoecological analyses in South America concluding that megafaunal extinctions did not occur until human appearance and warming coincided<sup>21</sup>, with humans potentially limiting species recovery after climate-induced stress<sup>22</sup>.

This result in such a continuous, spatially explicit context cannot be explored by using standard chronological approaches (i.e., without integrating a spatial component—the bearing—that connects all grid cells together instead of considering them as independent entities), because none of the 60 generalised least-squares models we tested explained any variation in the timing alone of megafauna extirpation in areas of human-megafauna coexistence (Table 1). Such an interaction between the regional patterns of a climate-driven variable (available freshwater) and those of human appearance suggests that the first humans in this region migrated across Australia by exploiting palaeolakes and other sources of drinkable water connecting the drier regions in between<sup>23</sup>. Given that approximatively three-quarters of Australia is semi-arid or arid, it is plausible that megafauna species were attracted to the same freshwater resources as were humans<sup>24</sup>, thus increasing the chances of encountering humans after the latter had arrived. We contend that these water points became hotspots of interspecific competition for water resources and human hunting that affected megafauna persistence depending on the particular configuration of freshwater resources in the local landscape.

Despite this effect of human appearance on the regional pattern of megafauna extirpations (Table 1), we found no effect of the duration of human-megafauna coexistence per region on the specific timing of megafauna extirpation. Thus, our results support the idea that human arrival in south-eastern Australia was an additional stressor to climate perturbations, rather than humans being a ‘super predator’<sup>4</sup>. Indeed, human populations were most likely small at the time of arrival due to low growth rates that are consistent with high mobility<sup>25</sup> compared to the period when human population growth increased following the Last Glacial Maximum (~23–19 ka) and the mid-Holocene climatic optimum (~7–5 ka)<sup>26,27</sup>.

An alternative hypothesis is that human populations were kept small by the landscape’s low carrying capacity typically experienced by species expanding into new environments<sup>28,29</sup>, as well as the rapid loss of a main food source following the demise of megafauna species<sup>30</sup> in the regions with the shortest windows of temporal coexistence. However, we found that low carrying capacity was mostly climate-driven due to lower mean annual precipitation and freshwater availability (see Supplementary Figure 6); this, along with high mean annual temperature, led to both lower net primary production and harsher desert conditions than during the Last Glacial Maximum (see Supplementary Fig. 7a’–e’). We argue that as for many other species<sup>31</sup>, human population density likely followed a slow growth pattern, quickly reaching a density plateau after establishing in the new environment due to competition with other species occupying the same landscape.
This alternative scenario of extinction is even more relevant in areas where climate was the only plausible driver of megafauna extinctions—in areas where there was an absence of temporal human-megafauna coexistence such as in Tasmania (Fig. 2a, blue areas, and ref. 32) because mean annual precipitation, mean freshwater availability, and mean annual desert fraction best explained the timing of megafauna extinctions there (Table 1 and Fig. 3b). In addition to local drought conditions, the pattern of primary production best explained the pattern of megafauna extirpations (Table 1, Fig. 3c and Supplementary Fig. 6d) during the prevailing cold conditions (i.e., up to 4°C lower than the conditions dominating during the Last Glacial Maximum; Supplementary Fig. 7) and low atmospheric CO₂ concentrations that would have limited plant growth³³ at that time.

We acknowledge that our approach has some potential limitations, such as disregarding the possible feedbacks of human-induced fire³⁴ and the complexity of interspecific interactions at the community level³⁵. Despite some available data describing fire regimes in Australia³⁶, our approach cannot distinguish human-from climate-driven fire patterns. Further, the low spatial coverage of high-quality megafauna fossil data in Australia still makes it impossible to generate outputs at the species level. However, by adding a spatial component (i.e., the bearing) we demonstrated a combined effect of climate and humans on the regional patterns of megafauna extirpations. In addition, we found that a shift in net primary production was the main driver of the regional pattern of megafauna extirpations in areas where humans were absent.

Our results therefore (i) highlight the limitations of standard (i.e., assuming no explicit spatial connectivity between sites) chronological analyses that are inadequate for describing the multidimensional complexity of megafauna extirpations in Australia. Our work also revealed that (ii) climatic effects on ecosystem resources are more complex than a simple interpretation that rainfall instead of temperature was the main climatic variable driving extirpations in Australia³⁷,³⁸. Here, freshwater availability essential to plant and animal life in an arid continent is the result of intricate physical processes between rainfall, evapotranspiration and geology. Finally, our work (iii) challenges the notion built from simple chronological analyses that anatomically modern humans were the principal drivers of extinctions¹,⁶,¹⁰ in Australia, and that climate variation was at best a secondary contributor⁹. Instead, we have revealed at least for south-eastern Australia a more complex scenario where climate change could have limited the amount of available resources for species, but that human appearance was likely another important and necessary contributor to explain megafauna extirpations in many parts of the landscape.

## Methods

**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)³⁷,³⁸. 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)¹³. 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 phenomena³⁹⁻⁴².

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 complexity¹⁹,³⁹. 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 world⁴⁰. 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 subtracting 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 Australia⁷) to capture the temporal variation in climate during the period of megafauna extinctions. These five climate variables characterise both natural resources for megafauna species⁴¹ and constraints to human appearance and movement²³,³⁹.

**Megafauna and human datasets.** From the 2138 records of extinct fauna species in the FosSahul database⁴², we quality-rated each fossil date following the A* to C scale (from 'high quality' to 'unreliable') based on objective criteria⁴³, 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. ⁴⁴, 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 contested⁴⁵.

**Estimating the timing of extirpation and initial arrival.** We applied a maximum-likelihood method to correct for the Signor-Lipps effect first developed by Solow⁴⁶ 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₁, ..., Tₙ, are independent and uniformly distributed over the interval (Tₐₕₓ, y), which correspond to the extirpation (Tₐₕₓ) and settlement time (y) 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 distributed⁴⁷. Thus, ̂Tₖ estimates the kth age assuming the true age Tₖ follows:

$$ \hat{T}_{k} \sim g(T_{k}, \sigma_{k}^{2}) \qquad (1) $$

with g being the Gaussian probability density function describing the radiometric error σ. The estimated timing of megafauna extinction ̂Tₐₕₓ and the estimated timing of initial human arrival ̂y are then calculated by numerically maximising the log-likelihood L(θ) over θ:

$$ \mathcal{L}(\theta) = \sum_{k=1}^{n} \log h(\hat{T}_{k}) \qquad (2) $$

with θ representing either ̂Tₐₕₓ or ̂y and h the probability density of ̂Tₖ

$$ h(\hat{T}_k) = \int_y^{T_{\text{ext}}} g(t, \sigma_k^2) \frac{dt}{T_{\text{ext}} - y} \qquad (3) $$

To adapt this approach to questions of spatial inference, let us define W as the spatial landscape, and Tₐₕₓ 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ₐₕₓ in each grid cell of W. However, to estimate Tₐₕₓ 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ₐₕₓ 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ₐₕₓ. 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 $(\hat{T}_1, ..., \hat{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 $\hat{T}_\theta(x)$ in a grid cell $x$ is calculated by numerically maximising $\mathcal{L}_x(T_\theta)$ over $\theta$ across $W$:

$$ \mathcal{L}_x(T_\theta) = \sum_{k=1}^{n} \log h_y(\hat{T}_k)^{w(x-x_k)} \quad (4) $$

where $w(x-x_k)$ is a weighting factor, such that $\Sigma w(x-x_k) = 1$. The standard procedure in estimating local density is to select a weighting factor proportional to:

$$ w(z) = \frac{1}{b} \times g_0\left(\frac{z}{b}\right) \quad (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 distribution$^{48,49}$.

Including a spatial component such as estimating local density generates a spatial bias in $\hat{T}_\theta(x)^{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 $\hat{T}_\theta(x)$, whereas narrower $b$ take fewer specimens into account in $w(z)$, thus increasing the variance of $\hat{T}_\theta(x)$. Here, the idea is to find an optimal size for $b$ to obtain the lowest values of both bias and variance into $\hat{T}_\theta(x)$. 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_{\theta(p)}$) for each grid cell. Based on these $T_{\theta(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_{\theta(p)}$, we calculated a mean bias (i.e., $B(\hat{T}_\theta(x))$) in each grid cell $x$ and the associated variance (i.e., $\sigma^2(\hat{T}_\theta(x))$) across space that are used to approximate the integrated mean-squared error across the landscape $W$ (E):

$$ E \approx \sum_x \left( \sigma^2 \hat{T}_\theta(x) + B \hat{T}_\theta(x) \right)^2 \quad (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_{\theta(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 \times 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 cell$^{51}$. 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 gradient$^{50}$.

**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: $w_{AIC_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

LOVECLIM$^{19,39}$ for the time period in each grid cell corresponding to the estimated timing of megafauna extirpation and its confidence interval.

We selected the top-ranked model as the model with the highest $w_{AIC_c}$ and tested this model for the effect of all possible interactions among predictor variables. We computed $w_{AIC_c}$ for these newly constructed models and ranked them against the top-ranked model without interactions to select the final top-ranked model. From this final model, we evaluated the relative importance of each predictor variable by calculating the change in likelihood (i.e., a metric for the model's goodness of fit) when this predictor is left out of the full model. This procedure is equivalent to a sensitivity analysis of the final model to each of its predictors to quantify and rank the relative importance of each predictor in the full model. We also accounted for a possible temporal lag between climate conditions and megafauna extirpation. Here, we calculated the temporal lag by regressing extirpation against climate from 1000 to 5000 years (at a 1000-year time step = 5 temporal-lag scenarios) for the period earlier than the estimated timing of megafauna extirpation in each grid cell (Fig. 1b). For each new time-lag scenario, we reconstructed and compared a new set of generalised least-squares models to determine which predictor among these best described spatial variation in the timing and bearing of megafauna extirpations.

## Data availability

The data supporting the Fig. 1 of this study are available within the paper as Supplementary Table 1. All data generated and analysed during this study as well as the raw data values underlying all reported averages in graphs and charts (i.e., Fig. 3 and Supplementary Figs. 4 and 5) are available at github.com/FredSaltre/SEOZ_megafauna_extirpation. The authors declare that the source data underlying Figs. 1a, 2a–d, 6d, h and 7c and Supplementary Figs. 1a and 5d are provided as a Source Data file.

## Code availability

The code is freely available on github.com/FredSaltre/SEOZ_megafauna_extirpation.

Received: 13 February 2019; Accepted: 25 October 2019;

Published online: 22 November 2019

## References

1. Prescott, G. W., Williams, D. R., Balmford, A., Green, R. E. & Manica, A. Quantitative global analysis of the role of climate and people in explaining late Quaternary megafaunal extinctions. Proc. Natl Acad. Sci. USA 109, 4527–4531 (2012).

2. Johnson, C. N. Ecological consequences of Late Quaternary extinctions of megafauna. Proc. R. Soc. Lond. B. 276, 2509–2519 (2009).

3. Wroe, S. et al. Climate change frames debate over the extinction of megafauna in Sahul (Pleistocene Australia-New Guinea). Proc. Natl Acad. Sci. USA 110, 8777–8781 (2013).

4. Martin, P. S. & Klein, R. G. Quaternary extinctions: a prehistoric revolution (University of Arizona Press, 1984).

5. Flannery, T. F. Pleistocene faunal loss: implications of the aftershock for Australia's past and future. Archaeol. Ocean. 25, 45–55 (1990).

6. Koch, P. L. & Barnosky, A. D. Late quaternary extinctions: state of the debate. Annu. Rev. Ecol. Evol. S 37, 215–250 (2006).

7. Saltré, F. et al. Climate change not to blame for late Quaternary megafauna extinctions in Australia. Nat. Commun. 7, 10511 (2016).

8. Araujo, B. B. A., Oliveira-Santos, L. G. R., Lima-Ribeiro, M. S., Diniz-Filho, J. A. F. & Fernandez, F. A. S. Bigger kill than chill: the uneven roles of humans and climate on late Quaternary megafaunal extinctions. Quatern. Int. 431, 216–222 (2017).

9. Bartlett, L. J. et al. Robustness despite uncertainty: regional climate data reveal the dominant role of humans in explaining global extinctions of Late Quaternary megafauna. Ecography 39, 152–161 (2016).

10. Sandom, C., Faurby, S., Sandel, B. & Svenning, J.-C. Global late Quaternary megafauna extinctions linked to humans, not climate change. Proc. R. Soc. Lond. B. 281, 20133254 (2014).

11. Raup, D. M. Taxonomic diversity during the phanerozoic. Science 177, 1065–1071 (1972).

12. Signor, P. W. & Lipps, J. H. in Geological implications of large asteroids and comets on the Earth (Silver, L. T. & Schultz, P. H. eds) (Geological Society of America, 1982).

13. Williams, A. N. et al. A continental narrative: human settlement patterns and Australian climate change over the last 35,000 years. Quat. Sci. Rev. 123, 91–112 (2015).

14. Wang, S. C. & Marshall, C. R. Estimating times of extinction in the fossil record. Biol. Lett. 12, 20150989 (2016).

15. Cooper, A. et al. Abrupt warming events drove Late Pleistocene Holartic megafaunal turnover. Science 349, 602–606 (2015).
16. Emery-Wetherell, M. M., McHorse, B. K. & Byrd Davis, E. Spatially explicit analysis sheds new light on the Pleistocene megafaunal extinction in North America. *Paleobiology* **43**, 1–14 (2017).

17. Benton, M. J. *Vertebrate Palaeontology* (Wiley, 2004).

18. Saltré, F. et al. Uncertainties in specimen dates constrain the choice of statistical method to infer extinction time. *Quat. Sci. Rev.* **112**, 128–137 (2015).

19. Goosse, H. et al. Description of the Earth system model of intermediate complexity LOVECLIM version 1.2. *Geosci. Model Dev.* **3**, 603–633 (2010).

20. Řičánková, V. P., Horsák, M., Hais, M., Robovský, J. & Chytrý, M. Environmental correlates of the Late Quaternary regional extinctions of large and small Palaearctic mammals. *Ecography* **41**, 516–527 (2018).

21. Metcalf, J. L. et al. Synergistic roles of climate warming and human occupation in Patagonian megafaunal extinctions during the Last Deglaciation. *Sci. Adv.* **2**, e1501682 (2016).

22. Raczka, M. F., Bush, M. B. & De Oliveira, P. E. The collapse of megafaunal populations in southeastern Brazil. *Quat. Res* **89**, 103–118 (2018).

23. Bird, M. I., O'Grady, D. & Ulm, S. Humans, water, and the colonization of Australia. *Proc. Natl Acad. Sci. USA* **113**, 11477–11482 (2016).

24. Prideaux, G. J. et al. Extinction implications of a chenopod browse diet for a giant Pleistocene kangaroo. *Proc. Natl Acad. Sci. USA* **106**, 11646–11650 (2009).

25. Tobler, R. et al. Aboriginal mitogenomes reveal 50,000 years of regionalism in Australia. *Nature* **544**, 180–184 (2017).

26. Williams, A. N., Ulm, S., Turney, C. S. M., Rohde, D. & White, G. Holocene demographic changes and the emergence of complex societies in prehistoric Australia. *PloS One* **10**, e0128661 (2015).

27. Reeves, J. M. et al. Palaeoenvironmental change in tropical Australasia over the last 30,000 years – a synthesis by the OZ-INTIMATE group. *Quat. Sci. Rev.* **74**, 97–114 (2013).

28. Tallavaara, M., Eronen, J. T. & Luoto, M. Productivity, biodiversity, and pathogens influence the global hunter-gatherer population density. *Proc. Natl Acad. Sci. USA* **115**, 1232–1237 (2018).

29. O'Connell, J. F. & Allen, J. The restaurant at the end of the universe: Modelling the colonisation of Sahul. *Austral Archaeol.* **74**, 5–31 (2012).

30. Goldberg, A., Mychajliw, A. M. & Hadly, E. A. Post-invasion demography of prehistoric humans in South America. *Nature* **532**, 232–235 (2016).

31. Arim, M., Abades, S. R., Neill, P. E., Lima, M. & Marquet, P. A. Spread dynamics of invasive species. *Proc. Natl Acad. Sci. USA* **103**, 374–378 (2006).

32. Lima-Ribeiro, M. S. & Diniz-Filho, J. A. F. Obstinate overkill in Tasmania? The closest gaps do not probabilistically support human involvement in megafaunal extinctions. *Earth Sci. Rev.* **135**, 59–64 (2014).

33. Johnson, C. N. et al. Geographic variation in the ecological effects of extinction of Australia's Pleistocene megafauna. *Ecography* **39**, 109–116 (2016).

34. Johnson, C. N. Fire, people and ecosystem change in Pleistocene Australia. *Aust. J. Bot.* **64**, 643–651 (2016).

35. Bliege Bird, R., Bird, D. W., Codding, B. F., Parker, C. H. & Jones, J. H. The "fire stick farming" hypothesis: Australian Aboriginal foraging strategies, biodiversity, and anthropogenic fire mosaics. *Proc. Natl Acad. Sci. USA* **105**, 14796–14801 (2008).

36. Mooney, S. D. et al. Late Quaternary fire regimes of Australasia. *Quat. Sci. Rev.* **30**, 28–46 (2011).

37. Kershaw, P., Moss, P. & Van Der Kaars, S. Causes and consequences of long-term climatic variability on the Australian continent. *Freshw. Biol.* **48**, 1274–1283 (2003).

38. Rule, S. et al. The aftermath of megafaunal extinction: ecosystem transformation in pleistocene Australia. *Science* **335**, 1483–1486 (2012).

39. Timmermann, A. & Friedrich, T. Late Pleistocene climate drivers of early human migration. *Nature* **538**, 92–95 (2016).

40. Fordham, D. A. et al. PaleoView: a tool for generating continuous climate projections spanning the last 21 000 years at regional and global scales. *Ecography* **40**, 1348–1358 (2017).

41. Johnson, C. N. *Australia's mammal extinctions: a 50,000 year history* (Cambridge University Press, 2006).

42. Peters, K. J. et al. FosSahul 2.0, an updated database for the Late Quaternary fossil records of Sahul. *Sci. Data* (2019), https://doi.org/10.1038/s41597-019-0267-3.

43. Rodriguez-Rey, M. et al. Criteria for assessing the quality of Middle Pleistocene to Holocene vertebrate fossil ages. *Quat. Geochronol.* **30**, 69–79 (2015).

44. Williams, A. & Ulm, S. AustArch: a database of 14C and non-14C ages from archaeological sites in Australia—composition, compilation and review. *Internet Archaeology* **36** (2014).

45. Gillespie, R. & Roberts, R. G. On the reliability of age estimates for human remains at Lake Mungo. *J. Hum. Evol.* **38**, 727–732 (2000).

46. Solow, A. R., Roberts, D. L. & Robbitt, K. M. On the Pleistocene extinctions of Alaskan mammoths and horses. *Proc. Natl Acad. Sci. USA* **103**, 7351–7353 (2006).

47. Walker, M. *Quaternary Dating Methods* (Wiley, 2005).

48. Härdle, W., Werwatz, A., Müller, M. & Sperlich, S. *In Nonparametric and Semiparametric Models* (Springer, Berlin Heidelberg, 2004).

49. Hu, F. & Zidek, J. The weighted likelihood. *Can. J. Stat.* **30**, 347–371 (2002).

50. Ordonez, A., Williams, J. W. & Svenning, J.-C. Mapping climatic mechanisms likely to favour the emergence of novel communities. *Nat. Clim. Chang.* **6**, 1104–1109 (2016).

51. Burrows, M. T. et al. The pace of shifting climate in Marine and Terrestrial Ecosystems. *Science* **334**, 652–655 (2011).

52. Fox, G. A., Negrete-Yankelevich, S. & Sosa, V. J. *Ecological Statistics: Contemporary Theory and Application* (Oxford University Press, 2015).

53. Burnham, K. P. & Anderson, D. R. *Model selection and multimodel inference: A practical information-theoretic approach* 2nd edn. (Springer-Verlag, 2010).

## Acknowledgements

We are grateful to Alejandro Ordonez Gloria (Aarhus University) for providing the resources to calculate the bearing of the spatial gradients. This research was supported by Australian Research Council (ARC) Centre of Excellence for Australian Biodiversity and Heritage (CE170100015) to S.U. and C.J.A.B. A.T. was supported by the Institute for Basic Science (IBS) under IBS-R028-D1.

## Author contributions

All authors contributed to drafting of the paper; F.S. designed the study, developed the model, analysed the data and led the writing; J.C. helped develop the models and statistical analyses; K.J.P. and M.M. compiled, vetted and cleaned the fossil dataset; T.F and A.T. provided the climate data from LOVECLIM and helped interpret the climate simulations; S.U. vetted the archaeological records and helped interpret the results; C.J.A.B. helped design the study and interpreted the results; all authors gave final approval for publication.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-019-13277-0.

Correspondence and requests for materials should be addressed to F.Sé.

Peer review information *Nature Communications* thanks Rainer Grun and Matheus Lima-Ribeiro for their contribution to the peer review of this work.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

Open Access This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made. The images or other third party
material in this article are included in the article's Creative Commons license, unless
indicated otherwise in a credit line to the material. If material is not included in the
article's Creative Commons license and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this license, visit http://creativecommons.org/
licenses/by/4.0/.

© The Author(s) 2019