In this study we address prior limitations by using regionally resolved global climate reconstructions and an improved database of megafaunal last appearance dates. Crucially, we focus not only on assessing the relative importance of the two extinction drivers at a global level, but also on identifying geographic areas where event chronologies are poorly understood. This approach is informed by three linked questions: 1) What are the absolute and relative explanatory powers of human colonisation and climatic changes as predictors of megafaunal extinction patterns? 2) How sensitive are these results to uncertainties in human arrival dates and last appearances of megafaunal genera? 3) Where and when do human and climatic factors accurately predict extinction patterns, and in which regions do they fail to do so?

It is possible to test the sensitivity of the results by repeating the analysis over a large number of scenarios. We can therefore explicitly take the uncertainty into account in the analysis, rather than trying to reconstruct unfeasibly precise chronologies. This approach was recently used to investigate the role of colonising dingoes Canis lupus dingo in Australian extinctions (Prowse et al. 2014 ), successfully demonstrating that findings were robust to uncertainties in the underlying data (Roberts 2014 ). Similarly, we used this technique to quantify the relative role of climatic and anthropogenic megafaunal extinctions at a global level (Prescott et al. 2012 ); however, the previous analysis was limited by coarse temporal resolution, restricted geographic coverage and a lack of region specific climate proxies (McGlone 2012 ).

Reconstructing a global chronology of events is challenging. Extinction and human arrival dates are difficult to estimate due to inherent uncertainties in dating techniques and the scarcity of megafaunal and human records (Barnosky et al. 2004 , Prescott et al. 2012 , Stuart 2014 ). Challenging taphonomic conditions or limited sampling efforts across much of the world mean many megafauna are poorly documented and so subject to large Signor–Lipps effects (Field et al. 2013 ), whilst there has been extensive debate on the reliability of last appearance dates even for well‐represented species, e.g. Coleodonta antiquitatis (Lister and Stuart 2013 ). Similar problems also occur when estimating human arrival dates. To avoid limiting the scope of our study to well researched regions, we use an analytical approach that explicitly accounts for these uncertainties in the dating of extinctions and human arrival.

To be considered an adequate explanation, any driver must account for the spatial and geographic patterns of extinction observed. To evaluate the importance of different drivers we therefore need reliable chronologies of megafaunal extinction, human colonisation and climatic change. This has been attempted for small groups of species and for single geographic regions, with conclusions from these studies failing to support a universal explanation of the extinctions (Brook and Bowman 2004 , Burney and Flannery 2005 , Koch and Barnosky 2006 ). Given that the extinctions are observed across most of the globe, we feel a larger scale analysis is more likely to yield conclusions that are robust to uncertainties and idiosyncrasies in the palaeoecological record. Such an approach can also highlight the regions where further study would be most constructive.

Two broad drivers of extinction have been proposed and extensively debated: Late Pleistocene and Holocene climatic change, or the arrival of anatomically modern humans (Grayson and Meltzer 2003 , Fiedel and Haynes 2004 , Burney and Flannery 2005 , Koch and Barnosky 2006 , Wroe and Field 2006 , Wroe et al. 2006 , Ugan and Byers 2007 , Pushkina and Raia 2008 , Nogués‐Bravo et al. 2010 , Haynes 2013 ). Despite five decades of research and debate (Leakey 1966 , 1967 , Martin 1966 , 1967 ), the relative importance of these drivers across the globe remains contentious (Boulanger and Lyman 2014 , Flores 2014 , Lima‐Ribeiro and Diniz‐Filho 2014 , Yule et al. 2014 ).

Our world has lost most of the large terrestrial animals present 100 kyr ago (Barnosky et al. 2004 ). Although their extinctions occurred over a remarkably short period of geological time (Martin and Wright 1967 , Martin and Klein 1984 , MacPhee 1999 ), they were asynchronous across the globe (Barnosky et al. 2004 ). There has been little ecological replacement of these megafauna, resulting in vacant ecological niches and physiological anachronisms in surviving animals (Lindstedt et al. 1991 ) and plants (Guimarães et al. 2008 , Johnson 2009 ). Parallel extinctions are not seen in small animals, plants and the marine realm (Koch and Barnosky 2006 ), indicating a high degree of selectivity, further narrowed by common life history traits and ecology amongst extinct species (Johnson 2002 ).

Finally, to identify consistency of predictive ability between different geographic regions, we performed a leave‐one‐out cross validation (LOOCV). We excluded each region in turn from the dataset, fitted the best performing model to the remaining regions, and used this model to predict extinction chronologies for the excluded region. We estimated goodness of fit for the ‘left out’ region as the correlation between its observed extinction proportions for each time interval and those predicted from the model fitted for other regions.

For each combination of human arrival and extinction scenarios, we compared models using qAICc (Aikake information criterion corrected for small sample size and based on quasilikelihood). Models combining climate and human arrival were the most informative in all cases, so we quantified the relative explanatory power of these two classes of predictors based on Nagelkerke's R 2 (a measure of explained variation accounting for the models’ non‐Gaussian error structure).

For each combination of extinction and arrival scenarios, we fitted a set of 12 GNMs with six combinations of predictor variables (Table 2 ). These included as predictors: no predictors; focal climate only; focal and lagged climate; human arrival only; focal climate plus human arrival; focal and lagged climate plus human arrival. All models were fitted with either a global intercept (assuming a single background rate of extinction for all regions), or an intercept of the form d /logArea (assuming the background absolute extinction rate in each region to be a non‐linear function of its area). No interaction effects between predictors were accommodated for in the models.

To account for uncertainty in last appearance dates we generated 1000 datasets (‘extinction scenarios’) by randomly sampling dates from the ranges obtained from the literature. For each extinction scenario we then calculated the proportion of genera going extinct in each region in each 4 kyr time interval, and fitted a generalised non‐linear model (GNM) to these data using the ‘gnm’ package in R (Turner and Firth), with a quasibinomial error structure to account for overdispersion. As predictors in the models, we used four climatic variables and human arrival. Climate variables were those described above during the focal 4 kyr time interval (‘focal climate’) and the previous time interval (‘lagged climate’), allowing for a lag effect of climatic conditions. The effect of human arrival was modelled as a function of landmass size using a Ricker function, of the form:where bgives the maximum effect of human arrival, whereas the exponential governs the speed at which this maximum is reached, and how quickly the effect dissipates. This allows the impact of colonisation to range from quickly reaching its maximum and then rapidly decaying, to rising gradually and then also dissipating slowly – Fig. 3 shows a representative scenario of this. We would expect the effect of human arrival to be faster and stronger in small regions (e.g. islands) compared to a slower and weaker effect in larger continents. However, we did not constrain the coefficients, allowing the model to fit all possible effects, including human arrival decreasing extinction rates.

We used a climate reconstruction based on the HadCM3 circulation model driven by changes to orbital configuration, atmospheric greenhouse gas concentrations, and ice‐sheet extent and sea level, reconstructed from a variety of palaeo‐archives (Singarayer and Valdes, Eriksson et al.). Climates over the past 80 kyr were reconstructed at a 1° × 1° scale in 2 kyr snapshots. For each grid cell, we extracted annual mean temperature for the snapshot, averaged for the 4 k time interval, and temperature change. Temperature change was quantified as the sum of the change during the two 2 kyr temperature spans between the three snapshots (0, 2 k, 4 kyr) in each 4 kyr interval; summarised as an equation of the form:whererepresents one of the three estimates of annual mean temperature in a given interval. This measure of change ensures that time intervals with both increases and decreases in temperature over the 4 kyr would have a large temperature change value, as well as those showing consistent warming or cooling. We then obtained region‐specific estimates of mean annual temperature and temperature change by taking the mean of each across all cells within a region. The temporal resolution of these 2 k snapshots, and the 4 k time intervals our analysis operated at, meant that we did not use simulations which included Heinrich or other millennial‐scale events. Further, simulations from the HadCM3 model including these events are only available up until the last glacial maximum (Singarayer and Valdes) and so would not be consistent across our 80 kyr analysis period.

We consulted published literature for evidence of human arrival in our 14 regions (Turney et al. 2001 , O'Connell and Allen 2004 , Cupper and Duncan 2006 , Bulbeck 2007 , Goebel et al. 2008 , Armitage et al. 2011 , Benazzi et al. 2011 , Higham et al. 2011 , Gillespie et al. 2012 , Kaifu and Fujita 2012 , Bueno et al. 2013 , Cooke et al. 2013 , Dewar et al. 2013 , Latorre et al. 2013 ) and constructed eight representative global colonisation scenarios to capture the range of plausible arrival dates (Table 1 ). Our climate model necessitated 4 kyr time intervals throughout the analysis (see below), so our arrival scenarios use dates rounded to their nearest 4 kyr interval. The breadth of these time intervals captured most of the variation in plausible arrival dates, with our scenarios differing mainly in arrival dates for Sahul and the Americas (Table 1 ).

Published dates of remains vary in reliability, leading some previous studies to employ rigorous selection criteria (Roberts et al. 2001 , Lister and Stuart 2008 ). However, considering the rarity of finds for many genera, such criteria can limit sample sizes and make inclusion of some regions very difficult. By explicitly accounting for uncertainties in dates, we were able to relax our selection criteria and include any published date established directly from remains. Where possible, we included the measurement and calibration uncertainties of the quoted date, either directly from the calibrated calendar date published, or by including any measurement uncertainty in our own calibrations (which were undertaken when dates were only published as uncalibrated). Where no directly dated finds were available, we used previously published broad extinction range estimates based on alternative indirect analyses (Supplementary material Appendix 1, Table A1). All genera in the database were present before the start of our 80 kyr analysis period (Prescott et al. 2012 ).

We compiled a database of last appearance dates of megafaunal genera through a comprehensive literature review. We searched for all published records of dated remains or extinction estimates for terrestrial animal genera potentially present in the past 80 kyr and with a maximum mass > = 40 kg (Supplementary material Appendix 1, Table A1). We used genera to avoid complications arising from ambiguity in identifying remains to the species level. We mapped last appearance dates and the presence of extant, sufficiently massive genera onto 14 different geographic regions (Fig. 1 ). These regions broadly followed terrestrial biogeographic provinces (Udvardy 1975 ) but were divided into higher resolutions where possible, or used country boundaries where limitations in tracing the precise location of dated fossils necessitated. Some islands were excluded due to the unreliability of reconstructing climate conditions for very small landmasses.

Leave one out cross validation of model performance across the 14 geographic zones, showing goodness of fit to a region of the model generated when that region is left out of model parameterisation (box plot) compared to the median of when it is included (red dot). The spread of the box plot represents variation across 1000 extinction scenarios generated to account for the uncertainty in the extinction dates.

A representative fit of human arrival scenario 1 (global early arrival), comparing models that only include climate (green line), only human arrival (blue line), and combining both effects (red line). Time of arrival of humans in different geographic zones is marked by a vertical dashed black line and yellow shading of the period after arrival (note that Africa is completely shaded, as anatomically modern humans were present before 80 kyr ago).

There was considerable variation between regions in how well models’ predictions matched the observed pattern (Fig. 5 ). In Europe, Tasmania, and to lesser extents Japan, Canada and Alaska, model predictions closely matched observed patterns (Fig. 5 ) across all human arrival and extinction scenarios (Fig. 6 ). In the regions with the shortest, most recent, and most severe extinctions (New Zealand, Madagascar and parts of the Americas) predictions were accurately timed but underestimated observed losses (Fig. 5 ). This observation varies little across extinction and arrival scenarios (Fig. 6 ), even in South America, where uncertainty in last appearance dates is highest (Supplementary material Appendix 1, Table A1). In regions with few extinctions, most notably central Asia and Indo‐Malaya, the models performed badly and overestimated levels of extinction (Fig. 5 ) regardless of arrival or extinction scenario (Fig. 6 ). Only mainland Australia and New Guinea, which saw the earliest extinction events in the analysis, show peaks of extinction at times not predicted by the model (Fig. 5 ). These are also the only regions where model performance shows a higher degree of sensitivity to human arrival and extinction scenarios (Fig. 6 ). Overall, the LOOCV analysis showed that model performance in any individual region is largely unaffected by the exclusion of that region when fitting the model (Fig. 6 ).

The relative explanatory power of human arrival and climate variables compared across eight different human colonisation scenarios. Variation in extinction probability solely attributed to human arrival in red, solely to climate in yellow, and explained by both human arrival and climate in orange. Error bars represent standard deviations across 1000 extinction scenarios generated to account for the uncertainty in the extinction dates.

Overall, the top‐performing models consistently explained a high proportion of the global variation in extinction patterns, (Nagelkerke R 2 ∼ 75%, Fig. 4 ). The majority of the models’ explanatory power was uniquely ascribable to human arrival, accounting for approximately 60% of the Nagelkerke R 2 . A smaller proportion (approximately 25%) was uniquely ascribable to climatic predictors with the remaining proportion (approximately 15%) ascribable to either climate or human predictors (Fig. 4 ). These findings were again highly consistent across both arrival and extinction scenarios.

In all models, the effect of human arrival on extinction rates was consistent with our expectations. On islands (∼ 62 400–∼ 786 000 km 2 ), our models predicted large increases in extinction rates peaking after approximately 8 kyr and decaying to very low levels after 16 kyr. However, on continental landmasses our models predicted prolonged periods of elevated extinction rates peaking at 10 k–12 kyr after arrival and persisting beyond 30 kyr. Figure 3 shows some representative fits.

All models with more than negligible support from qAICc included the effect of human arrival as well as both focal and lagged climate, irrespective of the extinction or arrival scenario considered (Fig. 2 ). The best supported models (qAICc weights ∼ 0.7) had a global intercept (regions had the same baseline extinction rate); the second best supported models (qAICc weights ∼ 0.3) included a region specific intercept affected by the area of a region, which behaved mostly as expected, with either higher probabilities of extinction in smaller regions, or with no consistent distinguishable effect across extinction scenarios (Supplementary material Appendix 2, Table A3.1–3.8). Thus, the roles of both limatic and anthropogenic drivers of extinction are strongly supported, even when the uncertainties in the data are accounted for. A summary of the variation of model parameter values across scenarios in these two models is presented in Supplementary material Appendix 2.

Discussion

The top‐performing model identifies a combination of human colonisation, focal and lagged climate as the most important predictors of extinctions. Importantly, our models are able to explain the data well (R2 70%), capturing most of the worldwide variation in timing and extent of extinctions despite the uncertainties of the palaeo‐archaeological record. While our previous study (Prescott et al. 2012) had slightly higher explanatory power, this was due to a coarser temporal and geographic resolution (10 kyr time intervals and six geographic regions).

The majority of our models’ explanatory power is uniquely attributable to human colonisation, with a large minority uniquely attributable to climate. The considerably higher explanatory power of human colonisation supports theories that favour global expansion of anatomically modern humans as the principal driver of extinctions, in agreement with previous global analyses (Koch and Barnosky 2006, Prescott et al. 2012, Sandom et al. 2014). However our analysis provides evidence that climate was an additional and important contributor to the extinction event.

Our analysis is the first of its kind to investigate the megafaunal extinction event using high resolution global climate reconstructions. In our view this is crucial to understanding climatic effects, given many of the changes of the Late Pleistocene. This is demonstrated by those regions where extinctions are well explained by our models through temporally separated, sequential impacts of climate changes and human colonisation (Europe, Tasmania, Japan, Canada and Alaska, Fig. 5). A higher temporal resolution also allowed us to identify a lagged effect of climate, providing evidence that the full impact of climatic changes on megafauna could take several thousand years to be realised. This finding can be of considerable importance to understanding the effects of ongoing anthropogenic climate change on extant species.

Whilst our models are generally very good at predicting the timings of extinction episodes, the strength of such episodes is underestimated for a few regions in two situations. Firstly, the model does not predict the complete megafaunal extinctions that occurred on Madagascar and New Zealand (Fig. 5). It is uncontroversial that human colonisation was the critical factor in most of these island extinctions (Worthy and Holdaway 2002, Crowley 2010, Allentoft et al. 2014). There are two likely explanations for this underestimation: our human impact curve might be inappropriate, either for the smallest landmasses or most recent extinctions; or ecological naivety, the inability to adapt to introduced novel predators after evolving in isolation and having lost defensive adaptions, could be important in determining extinction intensity as is seen on islands in the modern day (Courchamp et al. 2003). The models’ strong performance on other islands (Japan, Tasmania) supports the naivety explanation, as these islands have been far less evolutionarily isolated compared to Madagascar or New Zealand. However human arrival occurred much earlier on Japan and Tasmania, and it is possible that subsequent changes in human hunting behaviour led to the greater severity of the more recent island extinctions.

The model also under‐predicts extinction rates in the North American and South American regions, although again correctly captures the timing based both on climate and human arrival. The higher extinction rates observed when human and climatic stressors coincide are evidence that there may have been synergistic effects between the two, a process that is not accounted for by our model. This idea is well supported by ecological theory (Boulanger and Lyman 2014), and has previously been suggested as an important factor in accounting for megafaunal extinctions (Burney and Flannery 2005, Barnosky and Lindsey 2010, Lorenzen et al. 2011). A variety of mechanisms for such synergies have been proposed, including climate mediating human colonisation (Eriksson et al. 2012), and climatic stress reducing populations to sizes or ranges that are more vulnerable to overexploitation by humans (Lima‐Ribeiro et al. 2013).

For a few regions, model performance depends on arrival and extinction scenarios. In Australia, the models predict a human‐driven extinction peak earlier than is observed for the early arrival scenarios, but performs better for later ones. Climatic factors add little predictive power in this region. This could be due to a number of factors specific to Australia. Firstly, rainfall rather than temperature may be the most important climatic variable driving extinctions in Australia (Kershaw et al. 2003, Pack et al. 2003, Hesse et al. 2004, McGlone 2012, Rule et al. 2012). It is also possible that changes in early human culture, such as changing usage of fire, might have caused secondary extinction peaks (Webb 2008). This latter explanation could also apply to New Guinea. Notably, these are the earliest extinction peaks observed in the analysis and the high level of uncertainty in the dataset in these regions and may contribute to the poorer model performance, an explanation which is supported by the higher sensitivity of model performance in these regions to the different scenarios. These higher dating uncertainties have led some Australian studies to employ rigorous date selection criteria for remains (Roberts et al. 2001) and for some Australian dates to come under intense scrutiny (Gillespie et al. 2006, Brook et al. 2007). Given the early predicted extinction peak in some of our arrival scenarios for Australia, we repeated the analysis using more conservative last appearance dates for Australia from Wroe et al. (2013) (Supplementary material Appendix 3, Table A4). Examination of the predicted extinctions (Supplementary material Appendix 3, Fig. A1) and LOOCV analysis (Supplementary material Appendix 3, Fig. A2) from this repeated analysis shows the same qualitative results as our original dataset, the only notable difference being in the smoother observed extinction peak for Australia (which remains well predicted only for later arrival scenarios). The sensitivity of our findings to differences in inclusion criteria for remains is therefore very small.

Only across mainland Asia and Africa does the model perform poorly by predicting higher than observed extinction rates, with low levels of observed extinction in Indo‐Malaya, Africa and central Asia. Uniquely, ‘colonisation’ never occurred in Africa, with some recent studies rejecting the hypothesis of humans as drivers of extinction in Africa (Faith 2014). However, there is evidence for hominin impacts occurring in Africa before the start of our analysis, coinciding with the evolution of earlier hominins (Werdelin and Lewis 2013) and their technologies. The lower levels of extinction observed in regions with histories of earlier hominin populations (Wells and Stock 2007) may suggest this exposure reduced the impact of final colonisation by anatomically modern humans, an idea supported by other global analyses (Sandom et al. 2014). However, this explanation of lower Asian extinction rates is speculative, and our models’ performances in these regions could instead reflect the uncertainties of megafauna population extents and last appearances across Asia. We suggest that concentrated archaeological study showing where and when megafauna lived in these regions should be a priority for this field; currently, it is across temperate and tropical Asia that our understanding of these extinctions seems to break down.

Our results are highly robust to uncertainty in the palaeological record. Plausible variation in extinction and arrival scenarios has little impact on our conclusions on the absolute and relative explanatory power of humans or climate. On a regional basis, only parts of Sahul showed moderate sensitivity across different scenarios, despite large uncertainties in the datasets of other regions, e.g. South America. No specific regions were shown to be unusual in our LOOCV (Fig. 6), demonstrating a degree of consistency in the model's behaviour across all regions. Overall, the high degree of robustness across our results means that the overall conclusion of this analysis is unlikely to change qualitatively with improvements in data. Future alterations to specific human arrival or megafaunal extinction chronologies, both due to new archaeological finds or new analyses, will have little effect on our results and the nature of our conclusions, although they may affect the detailed narrative for some regions.