Relative demographic susceptibility does not explain the extinction chronology of Sahul’s megafauna

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Appendix 6
Appendix 7
Appendix 8
Appendix 9
Data availability
References
Article and author information
Metrics

Abstract

The causes of Sahul’s megafauna extinctions remain uncertain, although several interacting factors were likely responsible. To examine the relative support for hypotheses regarding plausible ecological mechanisms underlying these extinctions, we constructed the first stochastic, age-structured models for 13 extinct megafauna species from five functional/taxonomic groups, as well as 8 extant species within these groups for comparison. Perturbing specific demographic rates individually, we tested which species were more demographically susceptible to extinction, and then compared these relative sensitivities to the fossil-derived extinction chronology. Our models show that the macropodiformes were the least demographically susceptible to extinction, followed by carnivores, monotremes, vombatiform herbivores, and large birds. Five of the eight extant species were as or more susceptible than the extinct species. There was no clear relationship between extinction susceptibility and the extinction chronology for any perturbation scenario, while body mass and generation length explained much of the variation in relative risk. Our results reveal that the actual mechanisms leading to the observed extinction chronology were unlikely related to variation in demographic susceptibility per se, but were possibly driven instead by finer-scale variation in climate change and/or human prey choice and relative hunting success.

Introduction

The myriad mechanisms driving species extinctions (Sodhi et al., 2009) are often synergistic (Brook et al., 2008), spatially variable (Saltré et al., 2019), phylogenetically clumped (Fritz and Purvis, 2010), correlated with population size (O'Grady et al., 2004), and dependent on biotic interactions (Strona and Bradshaw, 2018). This complexity means that even in contemporary settings involving closely monitored species, identifying the ecological mechanisms underlying the causes of a particular extinction can be difficult (Caughley and Gunn, 1996; Fagan and Holmes, 2006). This challenge is considerably greater for palaeo-extinctions because of the restricted ecological knowledge of extinct species. In the case of prehistoric extinctions, we can only infer the environmental conditions likely operating at the estimated time of the events from rare and sparsely distributed proxies (Johnson et al., 2016).

The rapid and widespread disappearance of megafauna in the late Quaternary on most continents is one of the best-studied extinction events of the past. Many plausible causes of megafaunal extinction have been proposed (Koch and Barnosky, 2006). The main drivers of these extinctions appear to differ depending on taxon, region, and time period (Lorenzen et al., 2011; Metcalf et al., 2016; Wan and Zhang, 2017), but there is growing consensus that multiple drivers were involved, including the interactions between climatic shifts and novel human pressure as dominant mechanisms (Bartlett et al., 2016; Johnson et al., 2016; Saltré et al., 2019; Sandom et al., 2014; Villavicencio et al., 2016). This consensus mostly relies on approaches examining extinction chronologies relative to indices of temporal and spatial environmental variation. While such correlative approaches can suggest potential causes of extinction, they cannot by themselves provide strong inference on the plausible ecological processes involved. Instead, approaches that construct mechanistic models of environmental and other processes that drive extinctions could reveal the relative susceptibility of species over the course of a large extinction event (Alroy, 2001; Timmermann, 2020). Revealing such mechanisms not only provides evidence-based explanations of how and why past extinction events unfolded, it also assists contemporary and future ecological analyses by describing extinction processes over longer timeframes than historical records permit. This can further contextualise baselines for conservation and management targets, and assist in predicting the magnitude and sequences of future extinctions (Pardi and Smith, 2012; Willis and Birks, 2006; Wingard et al., 2017) as the current crisis unfolds (Bradshaw et al., 2021).

Existing mechanistic models applied to megafauna systems differ in their complexity, ranging from predator-prey models (Frank et al., 2015; Nogués-Bravo et al., 2008), to fully age-structured stochastic models (Prowse et al., 2013), or stochastic predator-prey-competition functions (Prowse et al., 2014). If sufficiently comprehensive, such models can be useful tools to compare the likelihood of the processes driving extinctions, even in the deeper past. Although measuring the demographic rates of long-extinct species is impossible, robust estimates can be approximated from modern analogues and allometric relationships derived from extant species (Prowse et al., 2014; Prowse et al., 2013). This makes it possible to construct stochastic demographic models of both extinct and related extant species, and compare their relative susceptibility to perturbations by mimicking particular environmental processes in silico.

Despite these methodological advances, unravelling the causes of the disappearance of megafauna from Sahul (the combined landmass of Australia and New Guinea joined during periods of low sea level) is still a major challenge. This is because of the event’s antiquity (Johnson et al., 2016) and the sparse palaeo-ecological information on megafauna extinctions compared to other parts of the world (Johnson et al., 2016; Saltré et al., 2019). However, based on the expectation that if high demographic susceptibility is an important feature of a species’ actual extinction dynamics, the most susceptible species should have gone extinct before more resilient species did.

Stochastic demographic models can quantify the relative support for the hypothesized mechanisms potentially involved in the megafauna disappearances in Sahul (summarized in Figure 1), which could assist in explaining the observed extinction chronology. These mechanisms include: (i) There is a life history pattern in which the slowest-reproducing species succumbed first to changing conditions (e.g., novel and efficient human hunting and/or climate change) (Cardillo et al., 2005; Purvis et al., 2000; Sodhi et al., 2009). This mechanism assumes that human hunting or climate change, even if non-selective, would differentially remove species that were more demographically sensitive to increased mortality (González-Suárez et al., 2013; Johnson, 2002). (ii) The most susceptible species were those whose life histories conferred the highest sensitivity to human hunting, such as species most sensitive to the loss of juveniles (Brook and Johnson, 2006). (iii) Bottom-up processes drove the extinctions, as supported by differences in the extinction timing between carnivores and their herbivore prey. Under this mechanism, prey-specialist carnivores should be more susceptible than their prey (i.e. because they depend on declining prey populations), whereas more generalist carnivores that could switch food sources would be less susceptible than their main prey (Chamberlain et al., 2005; Ripple and Van Valkenburgh, 2010). (iv) Species susceptible to temporal variation in climate succumbed before those most able to adapt to changing conditions. Under this mechanism, we expect the largest species — that is, those possessing traits associated with diet/habitat generalism (Monaco et al., 2020), physiological resilience to fluctuating food availability (Herfindal et al., 2006; Morris and Doak, 2004), high endurance, and rapid, efficient dispersal away from stressful conditions (Johnson, 2006) — would persist the longest in the face of catastrophic environmental change, independent of the intensity of human predation or changing climates. However, there are many other ways that climate change can alter demography (reviewed in Selwood et al., 2015). For example, heat stress can disadvantage juveniles by reducing parent capacity to raise neonates, or by affecting relatively heat-intolerant juveniles directly; species with more-sensitive juveniles are therefore expected to go extinct before more heat-tolerant species (Selwood et al., 2015). An increasing frequency of climate-induced catastrophes can also drive relatively smaller populations toward extinction faster, meaning large-bodied species with smaller populations are potentially more susceptible (Cardillo et al., 2005; Liow et al., 2008; Tomiya, 2013). (v) If none of the aforementioned mechanisms explains the extinction event’s chronology, non-demographic mechanisms such as differential selection of or ease of access by human hunters could have played more important roles.

Figure 1

Download asset Open asset

Description of five dominant mechanisms by which megafauna could have been driven to extinction and the associated seven perturbation scenarios examined.

To examine the relative support of these hypotheses, we developed the first stochastic, age-structured demographic models ever constructed for 13 extinct megafauna species in Sahul broadly categorized into five functional/taxonomic groups: (i) 4 vombatiform herbivores, (ii) 5 macropodiform herbivores, (iii) 1 large, flightless bird, (iv) 2 marsupial carnivores, and (v) 1 monotreme invertivore. We also built demographic models for 8 of some of the largest, local, extant species, including representatives from each of the functional/taxonomic groups described above for comparison. Our null hypothesis is that these extant species should demonstrate higher resilience to perturbations than the extinct species, given that they persisted through the main extinction event to the present. Subjecting each species’ model stochastically to different types of demographic perturbations, we tested seven scenarios (described in more details in Materials and methods) regarding the processes that could lead to extinction (see also Figure 1): an allometric relationship between the time of extinction and species’ body mass and/or generation length (Scenario ‘LH’), reducing juvenile survival (Scenario ‘↓S_j’) (Brook and Johnson, 2006; Munn and Dawson, 2001), reducing fertility (Scenario ‘↓F’) (Gittleman and Thompson, 1988; Miller et al., 2016; Oftedal, 1985), reducing survival across all ages (Scenario ‘↓S’), individual offtake from the population via hunting (Scenario ‘↓ind’), increasing environmental variability generating extreme climate events causing catastrophic mortality (Scenario ‘↑cat_f’), and increasing the magnitude of environment-driven, catastrophic mortality events (Scenario ‘↑cat_M’) (Reed et al., 2003).

We hypothesize that one, or several, of these types of perturbations would provide a better match than the others between relative demographic susceptibility and the continental-scale chronology of extinctions as inferred from the fossil record. Identifying which, if any, of the scenarios best matches the chronology would therefore indicate higher relative support for those mechanisms being the most likely involved in driving the observed extinctions. We first compared the expectation of larger (and therefore, slower life-history; Scenario LH, Figure 1) species more likely to go extinct than smaller species when faced with novel mortality sources (Brook and Bowman, 2005), followed by the outcomes of all other scenarios (Figure 1) to test if sensitivity to specific demographic changes supported other mechanisms.

Results

There was no indication that relatively heavier (Figure 2a) or slower life history (longer-generation; Figure 2b) species went extinct before lighter, faster life-history species (Scenario LH), even considering that the two mid- and small-sized carnivores Thylacinus and Sarcophilus went extinct on the mainland late in the Holocene at approximately the same time (~3200 years before present; Figure 2; White et al., 2018).

Figure 2

Download asset Open asset

Relationship between estimated date of species extinction (across the entire continent) and (a) body mass (M, kg) or (b) generation length (G, years) (Scenario LH).

Extinction-timing windows are estimated based on the agreement among six different models that correct for the Signor-Lipps effect (described in Materials and methods) in chronologies of quality-rated (Rodríguez-Rey et al., 2015) fossil dates for the studied taxa described in Peters et al., 2019. Here, we have depicted *Sarcophilus* as ‘extant’, even though it went extinct on the mainland >3000 years ago. Also shown are the approximate major climate periods and transitions: Marine Isotope Stage 3 (MIS 3), MIS 2 (including the Last Glacial Maximum, Antarctic Cold Reversal, and Younger Dryas), and the Holocene (including the period of sea level flooding when Tasmania separated from the mainland, and the relatively warm, wet, and climatically stable Holocene optimum).

To distil overall extinction risk per species in each scenario, we progressively increased the relevant perturbation and calculated the proportion of 10,000 stochastic model runs where the final population size fell below a quasi-extinction threshold (E_q) of 50 females (Frankham et al., 2014). The area under the resulting quasi-extinction-probability curves provides a scenario-specific representation of extinction risk across the entire range of the specific perturbation. The quasi-extinction curves for each species differed markedly in each perturbation scenario (Appendix 6—figure 1), although there were some similarities among scenarios. For example, in all scenarios except for fertility reduction (Scenario ↓F) and offtake (Scenario ↓ind), the smallest, extant carnivore Dasyurus was the least susceptible, whereas Genyornis was one of the most susceptible in four of the six scenarios (Appendix 6—figure 1).

Across all species, log₁₀ body mass explained some of the variance in the total area under the quasi-extinction curve (Figure 3) for the individual-removal (↓ind: evidence ratio [ER]=41.22, R² = 0.35; Figure 3d) and catastrophe-magnitude (↑cat_M: ER = 19.92, R² = 0.31; Figure 3f) scenarios, but less variance for the scenarios with reductions in juvenile (↓S_j: ER = 2.04, R² = 0.14; Figure 3a) and all-ages survival (↓S: ER = 5.05, R² = 0.21; Scenario iv; Figure 3c). There was little to no evidence for a relationship in the fertility-reduction (↓F) and catastrophe-frequency (↑cat_f) scenarios (ER ≲ 1; Figure 3b,e). The relationships between area under the quasi-extinction curve and log₁₀ generation length (G) were generally stronger (Figure 4). The strongest relationships here were for all-ages survival-reduction (↓S) and magnitude-of-catastrophe (↑cat_M) scenarios (ER >490, R² ≥0.49; Figure 4c,f), followed by weaker relationships (ER <11, R² ≤0.26) for Scenarios ↓S_j (Figure 4a), ↓ind (Figure 4d), and ↑cat_f (Figure 4e), and no evidence for a relationship in Scenario ↓F (ER <1; Figure 4b).

Figure 3

Download asset Open asset

Area under the quasi-extinction curve (from Appendix 6—figure 1) — extinction risk: ∫Pr(E_q)d(x) — as a function of body mass (M, kg) for (a) (Scenario ↓S_j) decreasing juvenile survival, (b) (Scenario ↓F) decreasing fertility, (c) (Scenario ↓S) decreasing survival across all age classes, (d) (Scenario ↓ind) increasing number of individuals removed year⁻¹, (e) (Scenario ↑cat_f) increasing frequency of catastrophic die-offs per generation, and f: (Scenario ↑cat_M) increasing magnitude of catastrophic die-offs.

Shown are the information-theoretic evidence ratios (ER) and variation explained (R²) for the lines of best fit (grey dashed) in each scenario. Here, we have depicted *Sarcophilus* as ‘extant’, even though it went extinct on the mainland >3000 years ago.

Figure 4

Download asset Open asset

Area under the quasi-extinction curve (from Appendix 6—figure 1) — extinction risk: ∫Pr(E_q)d(x) — as a function of generation length (G, years) for (a) (Scenario ↓S_j) decreasing juvenile survival, (b) (Scenario ↓F) decreasing fertility, (c) (Scenario ↓S) decreasing survival across all age classes, (d) (Scenario ↓ind) increasing number of individuals removed year⁻¹, (e) (Scenario ↑cat_f) increasing frequency of catastrophic die-offs per generation, and f: (Scenario ↑cat_M) increasing magnitude of catastrophic die-offs.

Shown are the information-theoretic evidence ratios (ER) and variation explained (R²) for the lines of best fit (grey dashed) in each scenario. Here, we have depicted *Sarcophilus* as ‘extant’, even though it went extinct on the mainland >3000 years ago.

Allometric scaling of extinction risk was also apparent within most taxonomic/functional groups. For example, Diprotodon had the highest extinction risk among the extinct vombatiform herbivores in every scenario except fertility reduction (Scenario ↓F; Figures 3b and 4b). Most species were relatively immune even to large reductions in fertility, except for Sarcophilus, Dasyurus, Vombatus, Alectura, Dromaius, and Genyornis (Figures 3b and 4b). In the offtake sub-scenario where we removed only juvenile individuals (↓ind_j), the results were qualitatively similar to Scenario ↓S_j where we progressively decreased juvenile survival (Appendix 7—figure 1), although the relative susceptibility for most species decreased from Scenario ↓S_j to ↓ind_j. However, susceptibility increased for the extinct carnivores Thylacinus and Thylacoleo, and remained approximately the same for Dasyurus and Notamacropus (Appendix 7—figure 1).

In the fertility-reduction sub-scenario (↓F_e) where we progressively increased the mean number of eggs removed per female per year in the bird species to emulate egg harvesting by humans, there was a progressively increasing susceptibility with body mass (Figure 5). Genyornis was clearly the most susceptible to extinction from this mechanism compared to the other two bird species (Figure 4).

Figure 5

Download asset Open asset

Inset: Increasing extinction risk for birds — quasi-extinction probability: Pr(E_q) — as a function of increasing the mean number of eggs harvested per female per year (Scenario ↓F_e).

The main graph shows the area under the quasi-extinction curve — extinction risk: ∫Pr(E_q)d(x) — as a function body mass (M, kg).

When the extinction dates are viewed relative to the extinction risk (quasi-extinction integrals) calculated for each scenario, there is no indication that the most susceptible species went extinct earlier in any perturbation scenario (Figure 6). Taking the sum of the quasi-extinction integrals across scenarios indicated that five of the eight extant species examined (Sarcophilus [extinct on mainland; extant in Tasmania], Dromaius, Alectura, Vombatus, Tachyglossus) had extinction risks that were equivalent or higher than most of the extinct species (Figure 7a). Taking the median rank of the quasi-extinction integral across scenarios generally indicated the lowest susceptibility in the macropodiformes (although the small, extant carnivore Dasyurus was consistently the least susceptible for all scenarios except fertility reduction; Figure 7a), followed by the carnivores (except Dasyurus), monotremes, vombatiform herbivores, and finally, large birds (Figure 7b). The carnivores had susceptibility ranks spread across most of the spectrum (Figure 7b).

Figure 6

Download asset Open asset

Relationship between estimated date of species extinction (across entire mainland) and area under the quasi-extinction curve (from Appendix 7—figure 1) — extinction risk: ∫Pr(E_q)d(x) — for (a) (Scenario ↓S_j) decreasing juvenile survival, (b) (Scenario ↓F) decreasing fertility, (c) (Scenario ↓S) decreasing survival across all age classes, (d) (Scenario ↓ind) increasing number of individuals removed year⁻¹, (e) (Scenario ↑cat_f) increasing frequency of catastrophic die-offs per generation, and f: (Scenario **↑cat_M**) increasing magnitude of catastrophic die-offs.

Extinction-timing windows are estimated based on the agreement among six different models that correct for the Signor-Lipps effect (described in Materials and methods) in chronologies of quality-rated (Rodríguez-Rey et al., 2015) fossil dates for the studied taxa described in Peters et al., 2019.

Figure 7

Download asset Open asset

Extinction susceptibility of the 21 modelled species.

(a) Sum of the areas under the quasi-extinction curve for each of the six scenarios considered — total extinction risk: Σ∫Pr(E_q)d(x) — for each of the 21 modelled species (^†extinct; ^♀extant; scenario abbreviations: ↓S_j = reducing juvenile survival; ↓F = reducing fertility; ↓S = reducing survival; ind = reducing number of individuals; ↑cat_f = increasing frequency of catastrophe; ↑cat_M = increasing magnitude of catastrophe); (b) median susceptibility rank across the six scenarios considered (where higher ranks = higher susceptibility to extinction) for each species (red = extinct; black = extant; outline-only bars = macropodiformes; solid bars = vombatiforms; angled crosshatching = birds; vertical crosshatching = carnivores; brick crosshatching = monotremes); (c) median susceptibility rank as a function of log₁₀ generation length (G, kg) — there was a weak correlation including all species (solid grey line 1), but a strong relationship removing *Sarcophilus* and *Dromaius* (dashed gray line 2) (information-theoretic evidence ratio [ER] and variance explained [R²] shown for each); (d) estimated date of species extinction (across entire continent) as a function of median susceptibility rank; taxonomic/functional groupings are indicated by coloured symbols and convex hulls (macropodids: brown; monotremes: grey; vombatiforms: green; birds: blue; carnivores: red). Extinction-timing windows are estimated based on the agreement among six different models that correct for the Signor-Lipps effect (described in Materials and methods) in chronologies of quality-rated (Rodríguez-Rey et al., 2015) fossil dates for the studied taxa described in Peters et al., 2019.

Expressed as a function of log₁₀ generation length, there was evidence for a moderate relationship across all species (ER = 12.07, R² = 0.27; Figure 7c), but removing the outlier species Sarcophilus and Dromaius resulted in a much stronger relationship (ER = 3.61 × 10⁵, R² = 0.76; Figure 7c). Susceptibility also tended to increase with body mass within a group, except for carnivores (Sarcophilus being the anomaly) and monotremes (Figure 7c). There was no indication of a pattern when extinction date is plotted against median susceptibility rank (Figure 7d).

Discussion

The megafauna species of Sahul demonstrate demographic susceptibility to extinction largely following expectations derived from threat risk in modern species — species with slower life histories have higher demographic risk to extinction on average (Cardillo et al., 2005; Liow et al., 2008; Purvis et al., 2000; Tomiya, 2013). Indeed, our models show a convincing relationship between a taxon’s overall relative demographic susceptibility to extinction and its generation length. While demography clearly must have played a role in the extinction of Sahul’s megafauna given that mainly the largest species succumbed (Lyons et al., 2016), our stochastic models revealed no clear relationship between relative demographic susceptibility and the order in which the taxa we considered went extinct (even after considering alternative approaches to calculate the window of extinction — Appendix 8—figure 1–3). In particular, we did not find that the most demographically susceptible species went extinct before more resilient species. We can only conclude that the actual extinction chronology must have instead been an emergent property of many interacting demographic rates, temporal and spatial variation in population abundance, particular environmental contexts, community interactions, and likely other traits (Sodhi et al., 2009). As different perturbations compromise different aspects of a species’ life history, its relative susceptibility to extinction compared to other species in its community varies in often unpredictable ways.

We can therefore reject the hypothesis that the continent-wide extinction chronology is explained by species’ relative demographic susceptibility and non-selective hunting by humans. If demographic susceptibility coupled with non-selective hunting were the primary causes of these extinctions, we would expect a relationship between vulnerability and extinction date whereby the more vulnerable species went extinct earlier — but we found no such relationship (Figure 7d). However, the comparison of susceptibility under increasing intensities of egg harvesting revealed the highest demographic risk from this type of activity for the extinct Genyornis compared to the extant Dromaius and Alectura birds, supporting the notion that egg harvesting (and not hunting of adults) might have been at least partially responsible for the demise of Genyornis (Miller et al., 2016). We can also reject the hypothesis that the largest, and therefore the most physiologically buffered and mobile species, were the most resilient given the lack of relationship with the inferred extinction chronology (Figures 2 and 7d). Neither did the species with the highest sensitivities to reductions arising from the various perturbation scenarios succumb earlier (Figures 5 and 6d).

This lack of relationship to the chronology, combined with the result that many of the extant species had some of the highest extinction susceptibilities, suggest that no obvious demographic properties can explain the taxon-specific timing within the Sahul extinction event of the Late Pleistocene. This opens the possibility that the chronology instead reflects either a random set of circumstances — that species succumbed to circumstantial combinations of stressors depending on the local perturbations experienced by particular populations (Saltré et al., 2019) — or even that the chronology is still insufficiently resolved. However, this conclusion does not accord well with the notion that Late Pleistocene megafauna extinctions were non-random and occurred at a much higher pace than background extinction rates (Johnson, 2006; Koch and Barnosky, 2006).

Another possibility is that in the case of human hunting, preferences for selecting or avoiding particular species, such as targeting larger species for more efficient returns (Broughton et al., 2011), could have overridden or interacted with intrinsic demographic susceptibility. Indeed, Lyons et al., 2016 concluded that either large-bodied mammals were selectively targeted by humans during the Late Quaternary, or that these species were relatively more vulnerable to human hunting than smaller-bodied species, or both. In addition, specific behavioural adaptations could potentially have rendered demographically high-risk species in fact less vulnerable to human hunting, such as the behavior of Vombatus to dig and defend burrows that were difficult to access by humans (Garvey et al., 2016) compared to larger burrowing or non-burrowing vombatiformes. In the case of the macropodiformes, interspecific variation in the type of locomotion — a trait not captured by our demographic models — could have contributed more to their relative susceptibility to human hunters. For example, the ability to hop at high velocities as in Osphranter rufus could have given it an escape advantage over the relatively slower sthenurine macropodiformes that likely employed more bipedal striding than hopping (Janis et al., 2020). Similar hypotheses regarding risk-persistence mismatches in multispecies simulations have been proposed for the Late Pleistocene megafauna extinction event in North America (Alroy, 2001).

Although marsupials are widely included in studies estimating the types of mammalian demographic relationships like those we used here (de Magalhães et al., 2007; Healy et al., 2014; McCarthy et al., 2008), more explicit consideration of their reproductive differences compared to placentals might further improve the resolution of future models. In particular, marsupials are born at the extreme altricial stage and complete most of their development ex utero through lactation (Tyndale-Biscoe, 2005). This might change the way the cost of reproduction is borne, because the unusually long period of marsupial lactation can reduce the cost of raising offspring per unit time (Cork and Dove, 1989; Weisbecker and Goswami, 2010). There is also an increasing number of small marsupial and placental mammals known to have gone extinct during the Pleistocene in Sahul (e.g. Cramb et al., 2009; Cramb et al., 2018). The inclusion of these smaller species in future demographic susceptibility/chronology of extinction analyses will likely be insightful, but unfortunately estimates of extinction dates — which require at least 8–10 reliability dated specimens (Bradshaw et al., 2012; Saltré et al., 2015) — are not yet available for these animals.

Our models, although age-structured, stochastic, and incorporating compensatory density feedback, are still simplified expressions of a species’ particular ecological and environmental contexts. As stated, our models are aspatial, yet we know that spatial processes are correlated with local extinctions across the landscape (Saltré et al., 2019). For example, large proboscideans like mammoths managed to persist well into the Holocene on island refugia despite having a high intrinsic extinction risk (Nogués-Bravo et al., 2008). It is therefore plausible that more localized measures of extinction risk, timing, and particular climate and habitat contexts (see Appendix 9 for an examination of demographic susceptibility relative to hindcasted climate trends) could reveal subtler demographic processes at work (Chase et al., 2020). However, Sahul’s fossil record is still generally too sparse at a regional level to test this properly (Peters et al., 2019; Rodríguez-Rey et al., 2016), nor do we have data indicating how spatial variation might have altered local expressions of demographic rates in long-extinct species.

Our models also ignore biotic dependencies such as predator-prey, plant-herbivore, and competition relationships that could have modified relative susceptibility in different ways depending on the community in question (Brook and Bowman, 2004; Choquenot and Bowman, 1998). Trophic community networks constructed for south-eastern Sahul show that bottom-up processes most strongly affect lower trophic levels, with their influence diminishing at higher trophic levels, although extinct carnivores were more vulnerable to coextinction than extant carnivores (Llewelyn et al., 2020).

The particulars of the Genyornis extinction are also still debatable given the possibility that the egg shells used to date the species (Miller et al., 2005) are potentially confounded with an extinct Progura megapode (Grellet-Tinner et al., 2016; Shute et al., 2017). However, removing ‘Genyornis’ from the extinction chronology makes no difference to our overall conclusions, but it is problematic for comparing relative extinction risk between the extinct Genyornis and the extant Dromaius and Alectura. In fact, by including an extant megapode (Alectura) in our model simulations, we determined that this much smaller (2.2 kg) and faster-reproducing species had a much lower extinction susceptibility than both Genyornis and Dromaius.

That we found no clear patterns among the extinct megafauna of Sahul to explain their relative extinction chronology supports the notion that, at least for mammals, risk can be high across all body masses depending on a species’ particular ecology (Davidson et al., 2009), even if relative extinction risk appears to follow allometric expectations (Brook and Bowman, 2005; Cardillo et al., 2005; Liow et al., 2008; Purvis et al., 2000; Tomiya, 2013) as we demonstrated clearly here (Figures 3, 5 and 7). By definition, the megafauna were generally large (>45 kg) species, yet neither their body mass or their correspondingly higher relative demographic susceptibility explains the extinction chronology in Sahul. Our approach also provides a template for assessing relative demographic susceptibility to extinction for other Sahul species that we did not consider here, and for those in other continents, that could reveal previously underappreciated dynamics and drivers. Nonetheless, more spatially and community-dependent models are still needed to provide a more complete picture of the dynamics of Late Pleistocene megafauna extinctions.

Materials and methods

Choice of species

Request a detailed protocol

Our first step was to choose enough extinct and extant species from the Sahul fossil record (Peters et al., 2019; Rodríguez-Rey et al., 2016) to represent a diversity of clades that were particularly affected during the main extinction event (estimated between 60 and 40 ka, where one ka = 1000 years ago) (Saltré et al., 2016). We also aimed to include at least one extant species within each functional/taxonomic group to compare extant with extinct species’ susceptibility. We settled on a total of 21 species (13 extinct; 8 extant) from five different functional/taxonomic groups: (i) 5 vombatiform herbivores, (ii) 7 macropodiform herbivores, (iii) 3 large birds, (iv) 4 carnivores, and (v) 2 monotreme invertivores. For a full list and justification of species chosen as well as the distribution of mean body masses, refer to Appendix 1.

Estimating demographic rates

Request a detailed protocol

To build age-structured population models for extinct taxa, we relied on different allometric, phylogenetic, and measured relationships to predict the plausible range of component demographic rates. For most extant marsupials, we relied mainly on the marsupial life-history database published in Fisher et al., 2001, but updated some values for some species with more recent sources (see below). A detailed description of how we estimated the necessary demographic rates and other ecological data to build the stochastic models is provided in Appendix 2, and a full table of all demographic values is provided in Appendix 2—table 1. We also provide a correlation matrix among demographic values across species (Appendix 2—table 2).

Age-structured (Leslie) population models

Request a detailed protocol

From the estimated demographic rates for each species, we constructed a pre-breeding, ω+1 (i) × ω+1 (j) element (representing ages from 0 to ω years old), Leslie transition matrix (M) for females only (males are demographically irrelevant assuming equal sex ratios). Fertilities (m_x) occupied the first row of the matrix, survival probabilities (S_x) occupied the sub-diagonal, and we set the final diagonal transition probability (M_i,j) to S_ω for all species except Vombatus, Thylacinus, and Sarcophilus for which we instead set the value to 0 to limit unrealistically high proportions of old individuals in the population, and the evidence for catastrophic mortality at ω for the latter two species (dasyurids) (Cockburn, 1997; Holz and Little, 1995; Oakwood et al., 2001). Multiplying M by a population vector n estimates total population size at each forecasted time step (Caswell, 2001). Here, we used n₀ = ADMw, where w = the right eigenvector of M (stable stage distribution), and A = the surface area of the study zone applied in the stochastic extinction scenarios — we arbitrarily chose A = 250,000 km² (500 km × 500 km; approximately 10% larger than the state of Victoria) so that the species with the lowest n₀ would have a population of at least several thousand individuals at the start of the simulations (see Appendix 3). We also included a compensatory density-feedback function in all simulations to avoid exponentially increasing populations (see Appendix 4).

Stochastic extinction scenarios

Request a detailed protocol

With the base M including density feedback tailored for each species, we perturbed various elements of their life histories to examine the relative support for hypotheses regarding plausible extinction drivers and pathways (see Figure 1). We first tested the relationship between extinction date and speed of life history as a baseline without any perturbation (Scenario LH), and then we generated six additional scenarios with perturbations. The second scenario (↓S_j) decreased juvenile (x = 0 to α−1) survival (plus a sub-scenario [↓ind_j] where we progressively removed individual juveniles from the population as we did for all individuals in Scenario ↓ind — see below). This scenario aims to emulate either food shortages of sufficient magnitude to make growing juveniles with higher relative energy and water demand (Munn and Dawson, 2001) succumb to environmental change more than adults, or from targeted hunting of juveniles by humans (Brook and Johnson, 2006). The third scenario (↓F) progressively reduces fertility to emulate food shortages lowering energetically demanding reproduction/lactation (Gittleman and Thompson, 1988; Oftedal, 1985). We also considered a sub-scenario (Scenario ↓F_e, see details below) for the category of large birds where we progressively increased the number of eggs harvested by humans (Miller et al., 2016). In Scenario ↓S, we progressively reduced survival across all age classes to examine the influence of an age-independent environmental stressor. Scenario ↓ind progressively removed individuals from the n population vector emulating offtake where animals are directly removed from the population to simulate human hunting (with age-relative offtake following the stable stage distribution of the target species). In Scenario ↑cat_f, we emulated how environmental variability would compromise populations via an increased relative (i.e. per generation) frequency of catastrophic die-offs by progressively increasing the number of catastrophic ~50% mortality events occurring per generation. Finally, Scenario ↑cat_M progressively increased the magnitude of the catastrophic mortality events to examine species’ responses to rising severity of catastrophes (Reed et al., 2003).

For Scenario ↓F_e, we estimated the egg-production component for Genyornis by calculating the proportion of total fecundity contributed by individual egg production in Dromaius (nest success of 0.406 × hatching probability of 0.419 = 0.17), and then multiplying this proportion by the total fertility estimated for Genyornis from Equation 11 — this produced an estimated per-individual annual egg production of 7.74 eggs for Genyornis (or, 7.74/2 = 3.87 eggs resulting in daughters). For Scenario ↑cat_M, we randomly allocated a multiplier of the expected frequency per generation (uniform sampling) derived from the species-specific range of multipliers identified in Scenario ↑cat_f (i.e. from one to the value where the species has an extinction probability = 1). In this way, we both standardized the relative risk among species and avoided cases where catastrophe frequency was insufficient to elicit any iterations without at least one extinction.

We ran 10,000 stochastic iterations of each model starting with allometrically predicted stable population size (see Appendix 3) divided into age classes according to the stable stage distribution. We projected all runs to $40 ⌊ G ⌉$ for each species (removing the first $⌈ G ⌋$ values as burn-in). In each scenario, we progressively increased the relevant perturbation and calculated the proportion of 10,000 stochastic model runs where the final population size fell below a quasi-extinction (E_q) of 50 female individuals (100 total individuals total assuming 1:1 sex ratios). This threshold is based on the updated minimum size below which a population cannot avoid inbreeding depression (Frankham et al., 2014). After calculating the per-increment probability of E_q in each of the seven scenarios, we calculated the total area under the quasi-extinction curve (integral) for each species as a scenario-specific representation of extinction risk across the entire range of the specific perturbation — this provides a single, relative value per species for comparison. Finally, we ranked the integrals among species in each scenario (lower ranks = higher resilience), and took the median rank as an index of resilience to extinction incorporating all scenario sensitivities into one value for each species.

Extinction dates

Request a detailed protocol

We compared the relative susceptibilities among all extinction scenarios, as well as the combined extinction-susceptibility rank of each species, to estimates of continental extinction times for the genera we examined. We took all estimates of continental extinction dates from the Signor-Lipps corrected values provided in Saltré et al., 2016; however, more recent continent-wide disappearance dates for Thylacinus and Sarcophilus are provided in White et al., 2018. The Signor-Lipps-correction of extinction timing is explained in detail in Saltré et al., 2016, but we briefly summarize it here. The mean date of extinction and its confidence intervals are derived from six frequentist methods to infer the timing of extinction — Strauss and Sadler, 1989, McCarthy, 1998, Marshall, 1997, Solow et al., 2006, McInerny et al., 2006, and the Gaussian-Resampled, Inverse-Weighted McInerny (GRIWM) model (Bradshaw et al., 2012) — from chronologies of dated fossils described in the FosSahul 2.0 database (Peters et al., 2019). All dates in a taxon’s chronology are first assessed for reliability, with dates of less than ‘A’ rejected for the calculation of the extinction window (Rodríguez-Rey et al., 2015). Each method calculates an extinction window for each species (taxon), and from these a window of cross-model agreement is calculated for every year (assuming that higher cross-model agreement indicates a greater likelihood of the true extinction date occurring during those times). We also examined the sensitivity of our overall results to uncertainty in extinction estimates by deriving a jack-knifed version of the GRIWM model (Bradshaw et al., 2012; Saltré et al., 2015) (Appendix 8).

We hypothesize that one, or several, of these types of perturbation scenarios would lead to a better match between the continental-scale chronology of extinctions as inferred from the fossil record compared to the simpler expectation of larger species with slower life-histories being more likely to go extinct than smaller species with faster life histories when faced with novel mortality sources (Scenario LH) (Brook and Bowman, 2005).

Data and code availability

Request a detailed protocol

All data and are R code needed to reproduce the analyses are available for download at github.com/cjabradshaw/MegafaunaSusceptibility.

Appendix 1

Choice of species and body mass distribution

Given data availability, we settled on a total of 21 species (13 extinct; 8 extant) from five different functional/taxonomic groups: (i) 5 vombatiform herbivores: Diprotodon optatum (2786 kg; extinct) (Wroe et al., 2004), Palorchestes azael (1000 kg; extinct) (Richards et al., 2019), Zygomaturus trilobus (500 kg; extinct) (Johnson, 2006), Phascolonus gigas (200 kg; extinct) (Johnson, 2006), and Vombatus ursinus (common wombat; 25 kg; extant) (McIlroy, 1996; Saran et al., 2011); (ii) 7 macropodiform herbivores: Procoptodon goliah (250 kg; extinct) (Johnson and Prideaux, 2004), Sthenurus stirlingi (150 kg; extinct) (Johnson and Prideaux, 2004), Protemnodon anak (130 kg; extinct) (Johnson, 2006), Simosthenurus occidentalis (120 kg; extinct) (Johnson, 2006), Metasthenurus newtonae (55 kg; extinct) (Johnson, 2006), Osphranter rufus (red kangaroo; 25 kg; extant) (McIlroy, 2008), and Notamacropus rufogriseus (red-necked wallaby; 14 kg; extant) (Strahan, 1991); (iii) 3 omnivorous (but primarily plant-eating) large birds: Genyornis newtoni (200 kg; extinct) (Johnson, 2006), Dromaius novaehollandiae (emu; 55 kg; extant) (Sales, 2007), and Alectura lathami (brush turkey; 2.2 kg; extant) (Jones et al., 1995); (iv) 4 carnivores: Thylacoleo carnifex (marsupial ‘lion’; 110 kg; extinct) (Johnson, 2006), Thylacinus cynocephalus (marsupial ‘tiger’; 20 kg; extinct) (Jones and Stoddart, 1998), Sarcophilus harrisii (devil; 6.1 kg; extinct in mainland Australia, but extant in Tasmania — see below) (Guiler, 1978), and Dasyurus maculatus (spotted-tail quoll; 2.0 kg; extant) (Belcher et al., 2008); and (v) 2 monotreme invertivores: Megalibgwilia ramsayi (11 kg; extinct) (Johnson, 2006), and Tachyglossus aculeatus (short-beaked echidna; 4.0 kg; extant) (Nicol and Andersen, 2007).

For each species, we identified the body mass of mature females. However, the sex of an extinct individual from its fossilized remains in many species is difficult to determine, especially when sample sizes are small (Alonso-Llamazares and Pablos, 2019). As such, we might have inadvertently assigned a female mass based on an estimated male mass, given evidence that there is a male bias in many fossil collections (Allentoft et al., 2010; Gower et al., 2019). For this reason, we attempted to cover a broad range of body masses among species to maximize the relative difference between them for comparison.

The two genera Thylacinus and Sarcophilus require special consideration in both the design of the analysis and the interpretation of the results. While Sarcophilus is extant, it is restricted to the island state of Tasmania that has been separated from mainland Australia since approximately 8–10 ka. However, the species went extinct on the mainland 3179 (±24) years ago, whether considering the species complex Sarcophilus harrisii and S. laniarius together or separately (White et al., 2018). Like Sarcophilus, Thylacinus went extinct on the mainland just over 3000 years ago (White et al., 2018) and persisted in Tasmania. However, unlike Sarcophilus, Thylacinus also went extinct in Tasmania in the 1930s. In our analyses we treat Sarcophilus as ‘extant’, with the proviso that it should be considered extinct on the mainland. Although we could also have treated Thylacinus as ‘extant’ in the sense that it persisted into historical times in Tasmania, we treat this species as extinct in our analyses.

The distribution of the body masses (M) across the 21 species (range: 1.68–2786 kg) was approximately log-Normal (Shapiro-Wilk Normality test on log₁₀M: W = 0.9804; p=0.9305; Appendix 1—figure 1). Median dates of continental (i.e., total species) extinction ranged from 51,470 (±3167 standard deviation) for Metasthenurus, to 3179 (±24) for Sarcophilus (mainland only; currently extant in Tasmania) (Saltré et al., 2016; White et al., 2018).

Appendix 1—figure 1

Download asset Open asset

Histogram of log₁₀ adult body masses (in kg) for the 21 species examined.

The distribution is approximately log-Normal (Shapiro-Wilk Normality test on log₁₀*M: W* = 0.9804; p=0.9305).

Appendix 2

Estimating demographic rates and population data as input parameters for the stochastic models

For each species, we first calculated the maximum rate of instantaneous population growth (r_m) using the following equation:

r_{m} = 10^{0.6914 - 0.2622 {log}_{10} M}

for mammals, where M = mass (g) (Hennemann, 1983). For birds, we used an optimization of the objective function based on age at first breeding (α, estimated as shown below), adult survival (s_ad, estimated as shown below), and the allometric constant for birds (Dillingham et al., 2016) a_rT = 1.107:

r_{m} = \log_{e} λ_{1} = \log_{e} (min | λ_{1} - e^{\frac{a_{r T}}{α + \frac{s_{ad}}{λ_{1} - s_{ad}}}} |)

We then calculated theoretical equilibrium population densities (D, km⁻²) based on the following:

D = 10^{4.196 - 0.740 \log_{10} M} / 2

for mammalian herbivores (M = body mass in g), where dividing by two predicts for females only (i.e., assumed 1:1 sex ratio) (Damuth, 1981). For large, flightless birds (i.e., Genyornis and Dromaius) (Latham et al., 2020), we applied:

D = 10^{3.65 - 0.82 \log_{10} M} / 2

and:

D = 10^{1.63 - 0.23 \log_{10} M} / 2

for omnivorous birds (i.e., Alectura) (Juanes, 1986) where M is in g. For mammalian carnivores, we applied:

D = e^{1.930 - 1.026 \log M} / 2

(M in kg), which we derived from Stephens et al., 2019. There were no specific invertivore or taxonomically specific equations to estimate D; however, we determined that the equation for the fitted 97.5 percentile in mammalian carnivores:

D = 10^{1.91 - 1.02 \log_{10} M} / 2

Stephens et al., 2019 provided is a reasonable D for female Tachyglossus = 9.9 km⁻². This is comparable to echidna densities measured for Kangaroo Island (4.4 females km⁻²) (Rismiller and McKelvey, 2000) and Tasmania (8.4 females km⁻²) (Nicol and Andersen, 2007). We therefore also used Equation 7 to predict D for Megalibgwilia. For a detailed description of the distribution and trends of equilibrium densities, population sizes, and biomasses across the modelled species, see Appendix 3.

We estimated the maximum age (ω) of each species according to:

ω = 10^{0.89 + 0.13 \log_{10} M}

for non-volant birds and mammals (M in g) (Healy et al., 2014), or

ω = {7.02 M}^{0.174}

for Alectura (M in g) (de Magalhães et al., 2007), the latter of which produces a ω that closely matches the maximum longevity of 25 years estimated for the similar-sized megapode Leipoa ocellata (malleefowl) (Bode and Brennan, 2011).

For other species, we made species- or group-specific adjustments to the estimates of ω: for Vombatus (McIlroy, 2008), we set ω = 26; the disparity between this and the ω derived from the allometric prediction (Equation 7 gives ω = 29) means we adjusted ω′=26/29ω for all vombatiform herbivores. Similarly for the macropodiformes, we scaled the predicted ω according to the degree of overprediction of the parameter from the equation for Osphranter rufus (for the latter, the equation predicted ω = 30, but in reality it is closer to 13) (Jonzén et al., 2010), meaning we adjusted ω′=13/30ω for all macropodiformes except Notamacropus; for Notamacropus rufogriseus, we took the mean maximum age (Grzimek, 1990) of 16. For Dromaius (Atlas of Living Australia, 2020),we set ω = 17; for Thylacoleo, we set ω = 17 to match female lion Panthera leo longevity (Packer et al., 1998); for Thylacinus (Corbett, 1995; Prowse et al., 2014), we set ω = 10, for Sarcophilus (Bradshaw and Brook, 2005) ω = 5, and Dasyurus ω = 4 given the catastrophic mortality at maximum lifespan characteristic of dasyurids (Cockburn, 1997; Holz and Little, 1995; Oakwood et al., 2001). In the case of Dasyurus (both D. maculatus and D. hallucatus), most females die in their third year, although some can persist into the fourth year (Cremona et al., 2017; Glen, 2008; Glen and Dickman, 2013; Moro et al., 2019) and maximum longevity in captivity can be up to 5 years (Way, 1988). For Tachyglossus, ω is extremely high compared to similar-sized marsupials or placentals: up to 50 years in captivity, and possibly 45 years in the wild (Nicol and Andersen, 2007); we set the latter value of ω = 45 to be conservative. For Megalibgwilia, we assumed that the underestimation for Tachyglossus according to Equation 8 would also apply, so here we set Megalibgwilia ω = (45/23)26 = 51 years.

We estimated fecundity (F; mean number of female neonates produced per year and per breeding female) for mammals (Allainé et al., 1987) as:

F = e^{2.719 - 0.211 \log M} / 2

dividing by two for daughters only (M in g). Although we used well-establish allometric relationships to derive our input parameters, most of these relationships are based on placental species. It is accepted that average life-history traits differ between similar-sized marsupials and placentals (Fisher et al., 2001), and we therefore estimated a correction factor for F and age at first breeding (α) (see Equation 12) for both the vombatiform and macropodiform herbivore groups separately (see Appendix 5 figures 1–3 for approach) using demographic data describing marsupial life histories (Fisher et al., 2001).

For Vombatus (McIlroy, 1996), we set F = 0.25 given a litter size = 1, an inter-birth interval of 2 years, and an assumed sex ratio of 1:1, from which we derived F for the extinct vombatiform herbivores (Appendix 5—figure 2). For Notamacropus rufogriseus, we set the average annual number of offspring = 1 multiplied by a 2.8% twinning rate (Catt, 1977), 1.3 based on an interbirth interval of 286 days, and an assumed 1:1 sex ratio. For Genyornis, we applied the following equation:

F = e^{2.35 + 0.17 \log M} / 2

(M in g) (Allainé et al., 1987). For Dromaius (Sales, 2007), we used the average of 6.7 eggs clutch⁻¹ and 3.4 clutches year⁻¹, an assumed sex ratio of 1:1, and nest success (0.406) and hatching probabilities (0.419) for ostriches (Kennou Sebei et al., 2009). For Alectura, we used the annual mean of 16.6 eggs pair⁻¹ for Leipoa ocellata (Bode and Brennan, 2011), a hatching success of 0.866 for Alectura lathami (Jones, 1988), and an assumed 1:1 sex ratio. For Thylacoleo, we applied the average litter size of 1 for large vombatiforms, and a 1:1 sex ratio; for Thylacinus and Sarcophilus, we applied the values of 3.42 progeny litter⁻¹ and the proportion (Lachish et al., 2009; Prowse et al., 2014) of adults reproducing year⁻¹ (0.91), and a 1:1 sex ratio. The allometric prediction (Equation 10) nearly matched the product of mean litter size (4.9) (Glen, 2008) and proportion females breeding (0.643) (Glen and Dickman, 2013) for Dasyurus, so we used the former. For Tachyglossus (Nicol and Morrow, 2012), we used the production of 1 egg breeding event (divided by two for daughters) multiplied by the probability of breeding = 0.55. For Megalibgwilia, we assumed that the overestimation for Tachyglossus according to Equation 10 would also apply, so here we set Megalibgwilia F = (0.275/0.659)0.532 = 0.222.

To estimate the age at first breeding (α), we used the following relationship for mammals (de Magalhães et al., 2007):

α = e^{- 1.34 + 2.14 \log M}

and

α = 0.214 M^{0.303}

for birds (de Magalhães et al., 2007) (M in g). For the macropodiforms, Equation 12 appeared to overestimate α by ~ 20% (see Appendix 5—figure 2), so we adjusted the extinct macropodiforms accordingly. For the vombatiform herbivores, Equation 12 performed more according to expectation. For example, the 1800–2000 kg female white rhinoceros (Ceratotherium simum) has α = 6–7 years (Wilson and Mittermeier, 2001), which is similar to the 7 years predicted by equation 12 for the 2786 kg Diprotodon optatum (but < α = 10–12 years for the > 6000 kg female African savanna elephant Loxodonta africana) (Asier, 2016). We also made species-specific adjustments to α for the extant species (or recently extinct in the case of Thylacinus). In the case of Vombatus (Roger et al., 2011), we set α = 2, α = 2 for Osphranter rufus (Jonzén et al., 2010), α = 1 for Notamacropus rufogriseus (Catt, 1977), α = 3 for Dromaius (Patodkar et al., 2009), α = 2 for Alectura based on data from Leipoa ocellata (Bode and Brennan, 2011; Frith, 1959), and α = 1 for Thylacinus (Lachish et al., 2009; Prowse et al., 2014), Sarcophilus (Bradshaw and Brook, 2005), and Dasyurus (Glen and Dickman, 2013). For Thylacinus and Sarcophilus, although α = 1, only a small proportion of females breed at this age (see below), so for most females α is in fact 2. But for Dasyurus, some females can become sexually mature and breed at < 12 months old (Moro et al., 2019), so we incorporated a modest capacity to reproduce in the year prior to age 1 (40% of total fertility). For Tachyglosssus, we set α = 3 based on evidence that echidnas take 3–5 years to reach adult mass (Nicol and Andersen, 2007), and only adults are observed to breed (Rismiller and McKelvey, 2000); as such, we set m = 0.5F in year 3, 0.75F in year 4, and m = F thereafter. As we did for ω and F, we estimated the bias between the allometrically predicted and measured α (equation 12) for Tachyglossus, and applied this to Megalibgwilia; however, rounding to the nearest year also means α = 3 for Meglibgwilia. The log₁₀ of the resulting α among species predicted their respective log₁₀ r_m well (R² = 0.73) (Appendix 5—figure 4), with the fitted parameters similar to theoretical expectation (Hone et al., 2010), and thus, supporting our estimates of r_m as realistic.

To estimate age-specific fertilities (m_x) from F and α, we fit a logistic power function of the general form:

m_{x} = \frac{a}{1 + {(\frac{x}{b})}^{c}}

where x = age in years, and a, b, c are constants estimated for each species, to a vector composed of (α−1) values at 0F, $⌊ α / 2 ⌉$ values at 0.75F, and for the remaining ages up to ω, the full value of F. This produced a continuous increase in m_x up to maximum rather than a less-realistic stepped series. For Sarcophilus, we instead used the parameters from an existing devil model (Bradshaw and Brook, 2005) to populate the m_x vector.

To estimate realistic survival schedules, we first used the allometric prediction of adult survival (s_ad) as:

S_{ad} = e^{- e^{- 0.5 - 0.25 \log M}}

for mammals, and:

S_{ad} = e^{- e^{- 1.78 - 0.21 \log M}}

for birds, where M = body mass (g) (McCarthy et al., 2008). For Tachyglossus, we used the mean S_ad = 0.96 based on upper and lower estimates of mortality for tagged individuals over 15 years in Tasmania (Nicol and Morrow, 2012), and applying the allometric-bias correction for Megalibgwilia as for ω, F, and α as described above. We then applied the Siler hazard model (Gurven and Kaplan, 2007) to estimate the age- (x-) specific proportion of surviving individuals (l_x); this combines survival schedules for immature, mature, and senescent individuals within the population:

l_{x} = e^{(\frac{- a_{1}}{b_{1}}) (1 - e^{- b_{1} x})} e^{- a_{2} x} e^{(\frac{a_{3}}{b_{3}}) (1 - e^{b_{3} x})}

where a₁ = initial immature mortality, b₁ = rate of mortality decline in immatures, a₂ = the age-independent mortality due to environmental variation, a₃ = initial adult mortality, and b₃ = the rate of mortality increase (senescence). From l_x, age-specific survival can be estimated as:

S_{x} = 1 - \frac{(l_{x} - l_{x + 1})}{l_{x}}

We estimated the component parameters starting with 1 – S_ad for a₁ and a₂, adjusting the other parameters in turn to produce a dominant eigen value (λ₁) from the transition matrix containing S_x (see Materials and methods) such that log_e λ₁ ≈ r_m. However, in many cases, marsupial and monotreme life histories were incapable of reproducing predicted r_m (see Appendix 2—table 1 below), although we attempted to maximize log_e λ₁ wherever possible. This appears to be biologically justified given the slower life histories of vombatiforms and monotremes in particular compared to macropodiforms. We also generally favoured a more pronounced senescence component of l_x in the longer-lived species given evidence for survival senescence in long-lived mammals (Turbill and Ruf, 2010). For Sarcophilus, we instead used the parameters from an existing devil model (Bradshaw and Brook, 2005) to populate the S_x vector.

Appendix 2—table 1

Predicted demographic values for each species (equations provided in Materials and methods and this appendix).

M = mass, r_m = maximum rate of instaneous exponential population growth predicted allometrically, $r_{m}^{'}$ = realised r_m predicted from the constructed matrix (see text), ω = longevity, F = fertility (daughters per breeding female per year), α = age at first reproduction (primiparity), S_ad = yearly adult survival, G = generation length. ^†extinct; ^♀extant. See Appendix 2—table 2 for rank correlations among demographic values across species.

Species	M (kg)	r_m	$r_{m}^{'}$	D (km⁻²)	ω (yrs)	F (n_♀yr⁻¹♀⁻¹)	α (yrs)	S_ad (yr⁻¹)	G (yrs)
vombatiform herbivores
Diprotodon^†	2786	0.100	0.061	0.134	48	0.1311	7	0.985	18.1
Palorchestes^†	1000	0.131	0.077	0.285	42	0.1705	6	0.981	15.1
Zygomaturus^†	500	0.157	0.095	0.476	39	0.2038	5	0.977	13.2
Phascolonus^†	200	0.200	0.121	0.938	34	0.2586	4	0.972	10.7
Vombatus^♀	25	0.345	0.119	4.370	26	0.2500	2	0.953	10.0
macropodiform herbivores
Proctoptodon^†	250	0.189	0.188	0.795	17	0.524	3	0.973	8.3
Sthenurus^†	150	0.216	0.215	1.161	17	0.617	3	0.970	8.1
Protemnodon^†	130	0.224	0.224	1.290	16	0.646	3	0.969	7.8
Simosthenurus^†	120	0.229	0.226	1.369	16	0.663	3	0.968	7.8
Metasthenurus^†	55	0.281	0.280	2.438	14	0.858	2	0.961	6.0
Osphranter^♀	25	0.345	0.343	4.370	13	0.750	2	0.953	5.5
Notamacropus^♀	14	0.402	0.351	6.712	16	0.668	1	0.993	6.3
large birds
Genyornis^†	200	0.041	0.041	0.101	38	0.658	9	0.987	20.0
Dromaius^♀	55	0.100	0.100	0.290	17	1.938	3	0.983	5.9
Alectura^♀	2.2	0.176	0.175	3.633	27	7.188	2	0.967	6.8
carnivores
Thylacoleo^†	110	0.234	0.201	0.028	14	0.500	4	0.967	9.1
Thylacinus^†	20	0.366	0.368	0.159	10	1.556	1	0.950	5.2
Sarcophilus^†,♀	6.1	0.500	0.094	0.539	5	1.205	1	0.820	3.1
Dasyurus^♀	2.0	0.701	0.644	2.023	4	1.582	1	0.910	2.3
monotremes
Megalibgwilia^†	11.0	0.307	0.107	3.522	51	0.222	3	0.977	16.4
Tachyglossus^♀	4.0	0.400	0.112	9.883	45	0.275	3	0.950	14.1

Sarcophilus harrisii is extinct in mainland Australia, but extant in the island state of Tasmania.

Thylacinus could also be treated like Sarcophilus in that Thylacinus survived in Tasmania until historical times (1930s).
In the case of the vombatiform and macropodiform herbivores, ω shown in the table is in fact the downscaled ω′ calculated for each group (see below). Likewise, both allometric predictions of F and α are corrected for these groups (see Supplementary Information Appendix 5).

Appendix 2—table 2

Kendall’s rank correlation (τ) among demographic values given in Appendix 2—table 1 across species.

	M (kg)	r_m	rm'	D (km⁻²)	ω (yrs)	F (n_♀yr⁻¹♀⁻¹)	α (yrs)	S_ad (yr⁻¹)
r_m	−0.684
r_m'	−0.350	0.507
D (km⁻²)	−0.476	0.499	0.358
ω (yrs)	0.390	−0.511	−0.605	−0.068
F (n_♀yr⁻¹♀⁻¹)	−0.484	0.278	0.486	0.148	−0.567
α (yrs)	0.671	−0.685	−0.583	−0.459	0.556	−0.562
S_ad (yr⁻¹)	0.552	−0.681	−0.394	−0.328	0.484	−0.327	0.572
G (yrs)	0.471	−0.446	−0.616	−0.201	0.762	−0.711	0.730	0.487

Appendix 3

Allometric predictions of equilibrium population density, total population size, and biomass

The allometric predictions of stable population size (N_stable) for each species in the 500 km × 500 km study area showed the largest populations for some of the smallest, extant species (e.g., N_stable > 600,000,000 for Vombatus, Osphranter, Notamacropus, Alectura, Dasyurus, Tachyglossus) (Appendix 3—figure 1a). There was a clear separation in the allometric predictions of N_stable among the species in each group (Appendix 3—figure 1b). When expressed as total biomass across the study area, the four carnivores had approximately equal biomasses (~10⁶ kg), as did the macropodiformes and monotremes (Appendix 3—figure 1c). For the large birds and herbivore vombatiformes, biomass increased with body mass (Appendix 3—figure 1c).

Appendix 3—figure 1

Download asset Open asset

Mass-predicted abundance and biomass for the 21 modelled species.

(a) Stable population sizes (N_stable) for each modelled species predicted from allometric estimates of population density for a 500 km × 500 km (250,000 km²) landscape (outline-only bars = macropodiformes; solid bars = vombatiforms; angled crosshatching = birds; vertical crosshatching = carnivores; brick crosshatching = monotremes); (b) N_stable plotted against body mass (M, in kg), showing the allometric scaling separating the vombatiform herbivores (green)/macropodiformes (brown), flightless birds (blue), carnivores (red), and monotremes (grey); (c) predicted landscape biomass (N_stable × M) for each species (outline-only bars = macropodiformes; solid bars = vombatiforms; angled crosshatching = birds; vertical crosshatching = carnivores; brick crosshatching = monotremes).

Here, we have depicted *Sarcophilus* as ‘extant’, even though it went extinct on the mainland >3000 years ago.

Appendix 4

Compensatory density feedback

To avoid an exponentially increasing population without limit generated by a transition matrix optimized to produce values as close to r_m as possible, we applied a theoretical compensatory density-feedback function. This procedure ensures that the long-term population dynamics were approximately stable by creating a second logistic function of the same form as m_x to calculate a modifier (S_mod) of the S_x vector according to total population size (Σn):

S_{mod} = \frac{a}{1 + {(\frac{\sum n}{b})}^{c}}

We adjusted the a, b, and c constants for each species in turn so that a stochastic projection of the population remained stable on average for 40 generations ( $40 ⌊ G ⌉$ ), where:

G = \frac{\log ({(v^{T} M)}_{1})}{λ_{1}}

and (v^TM)₁ = the dominant eigen value of the reproductive matrix R derived from M, and v = the left eigenvector (Caswell, 2001) of M. Although arbitrary, we chose a $40 ⌊ G ⌉$ projection time as a convention of population viability analysis to standardize across different life histories (Brook et al., 2006; Traill et al., 2010).

The projections were stochastic in that we β-resampled the S_x vector assuming a 5% standard deviation of each S_x and Gaussian-resampled the m_x vector at each yearly time step to $40 ⌊ G ⌉$ . We also added a catastrophic die-off function to account for the probability of catastrophic mortality events (C) scaling to generation length among vertebrates (Reed et al., 2003):

C = \frac{p_{C}}{G}

where p_C = probability of catastrophe (Reed et al., 2003) (set at 0.14). Once invoked at probability C, we applied a β-resampled proportion centred on 0.5 to the β-resampled S_x vector to induce a ~ 50% mortality event for that year (Bradshaw et al., 2013), as we assumed that a catastrophic event is defined as “… any 1 yr peak-to-trough decline in estimated numbers of 50% or greater” (Reed et al., 2003). Finally, for each species we rejected the first $⌊ G ⌉$ years of the projection as a burn-in to allow the initial (deterministic) stable stage distribution to stabilize to the stochastic expression of stability under compensatory density feedback.

Appendix 5

Deriving marsupial correction factors for fecundity (F) and age at first breeding (α)

Given that allometric predictions of various life-history traits in mammals are based primarily on data from extant placentals, we investigated the degree of potential bias in our estimates of longevity (ω) fecundity, (F) and age at first breeding (α) based on a comparison of theoretical and observed data for extant marsupials. We discuss the bias correction ω for in the main text for the vombatiform and macropodiform herbivores (see Appendix 2), so here we report how we derived group-specific corrections for F and α.

Fertility (F) correction

We first collected adult female mass, inter-birth interval (I_b), and age at first breeding data for twenty-three extant species in the database compiled by Fisher et al., 2001. We included all species for which data were listed in the genera: Macropus (Osphranter and Notamacropus), Dorcopsis, Lagorchestes, Petrogale, Thylogale, and Wallabia. We excluded the genus Dendrolagus (tree kangaroos) because they represent a distinct clade and differ ecologically from most other macropodids. We also excluded the tammar wallaby (Macropus eugenii) because it is a strongly seasonal-breeding species that potentially strong leverage on estimating the allometric slope.

To correct F, we first examined the relationship between I_b and body mass (M) for these species (Appendix 5—figure 1):

We therefore concluded that there was sufficient evidence for an allometric relationship between the two variables for this group, which we used to project the degree to which F was overestimated by the allometric relationship (Equation 10) used to estimate F for the extinct macropodidiform herbivores. Using the intercept and slope estimated in the relationship shown in Appendix 5—figure 1, we predicted an inter-birth interval of 384 days for Procoptodon, 363 days for Sthenurus, 357 days for Protemnodon, 354 days for Simosthenurus, and 322 days for Metasthenurus. These changed the allometrically predicted F by −4.9%, +0.5%, +2.2%, +3.1%, and +13.3%, respectively (corrected F shown in Appendix 2—table 1).

Appendix 5—figure 1

Download asset Open asset

Relationship between the logarithm of adult female body mass (M, kg) and inter-birth interval (I_b, in days) for 23 extant macropodid herbivores (Fisher et al., 2001).

The estimated parameters of the linear fit (y ~ *α + βx*) are: α = 159.3 ± 31.5 days (± SE) and β = 93.6 ± 36.4, explaining 24.2% of the variation ( $R_{adj}^{2}$ ), with the information-theoretic evidence ratio (ER) of the slope *versus* intercept-only model = 11.0.

For the vombatiforms, there are only four extant phascolarctid (koala Phascolarctos cinereus) and vombatiform herbivores (common or bare-nosed wombat Vombatus ursinus, northern hairy-nosed wombat Lasiorhinus krefftii, and southern hairy-nosed wombat L. latifrons). There were not enough species to estimate an allometric relationship that might predict the expected I_b for extinct vombatiform herbivores, so instead we assumed that the extinct vombatiform herbivores we considered would scale allometrically relative to Vombatus ursinus, which has a measured I_b of 730 days (Fisher et al., 2001). For this, we assumed the same slope as measured for the extant macropodiforms (Appendix 5—figure 1), and an intercept that aligned Vombatus with the relationship (Appendix 5—figure 2):

Appendix 5—figure 2

Download asset Open asset

Relationship between the logarithm of adult female body mass (M, kg) and inter-birth interval (I_b, in days) for four large, extant phoscolarctid and vombatiform herbivores (black circles): koala *Phascolarctos cinereus*; common wombat *Vombatus ursinus*; northern hairy-nosed wombat *Lasiorhinus krefftii*; southern hairy-nosed wombat *L. latifrons*.

Shown is the assumed relationship between I_b and log₁₀M setting the slope to that estimated for the extant macropodiforms (β = 93.6; Appendix 5—figure 1) and an intercept that aligned with *Vombatus* (α = 599 days) to estimate the inter-birth interval for the four extinct vombatiform herbivores (red circles) considered in this analysis.

Age at first breeding (α) correction

Next, we estimated the bias in the predicted age at first breeding (α) for the macropodiforms. A similar correction for the vombatiform herbivores was not warranted given that the allometric predictions were close to expectation for placentals of similar mass (see main text). We plotted α predicted from the allometric Equation 12 against those observed from the marsupial database (Fisher et al., 2001) for the extant macropodiform species as described above for F (Appendix 5—figure 3):

Appendix 5—figure 3

Download asset Open asset

Relationship between the predicted age (years) at first breeding (α_pred) and observed α (α_obs) for 23 extant macropodid herbivores (Fisher et al., 2001).

The allometric prediction over-estimated α by and average of ~20%. Also shown is the 1:1 line (dashed).

Calculating the average disparity between the predicted and observed α across species, the allometric prediction over-estimated α by 20% for the macropodiform herbivores. We therefore applied this correction factor to the estimated α for the extinct macropodiform species (corrected values in Appendix 2—table 1). The corrected α predicted maximum rate of population growth (r_m) approximately following theoretical expectation (see Appendix 5—figure 4).

Appendix 5—figure 4

Download asset Open asset

Negative relationship between the logarithm of the maximum rate of intrinsic population growth (log₁₀r_m) and the logarithm of the age at primiparity (log₁₀α) for the 21 species examined.

The estimated parameters of the linear fit (y ~ *α + βx*) including all species (black lines) are: α = -0.388 ± 0.075 (± SE) and β = -0.931 ± 0.146, and explaining 66.4% of the variation ( $R_{adj}^{2}$ ), with the information-theoretic evidence ratio (ER) of the slope *versus* intercept-only model = 4.052×10⁴. This relationship is similar to the theoretical expectation for the intercept = -0.15 and slope = -1.0 for mammals (Hone et al., 2010). Birds (*Alectura*, *Dromaius*, *Genyornis*; ◆) potentially fall outside this relationship, so just considering the remaining mammals (●), the parameters for the linear fit (red lines) become: α = -0.388 ± 0.083 (± SE), β = -0.875 ± 0.170, $R_{adj}^{2}$ = 60.1%, and ER = 1.567×10³.

Appendix 6

Quasi-extinction curves for each species in each of the six perturbation scenarios

Appendix 6—figure 1

Download asset Open asset

Increasing probabilities of quasi-extinction — quasi-extinction probability: Pr(E_q) — as a function of (a) decreasing juvenile survival (Scenario ↓S_j), (b) decreasing fertility (Scenario ↓F), (c) decreasing survival across all age classes (Scenario ↓S), (d) increasing number of individuals removed year⁻¹ (Scenario ↓ind), (e) increasing frequency of catastrophic die-offs per generation (Scenario ↑cat_f), and (f) increasing magnitude of catastrophic die-offs (Scenario ↑cat_M).

Species notation: DP = *Diprotodon optatum*, PA = *Palorchestes azael*, ZT = *Zygomaturus trilobus*, PH = *Phascolonus gigas*, VU = *Vombatus ursinus* (vombatiform herbivores); PG = *Procoptodon goliah*, SS = *Sthenurus stirlingi*, PT = *Protemnodon anak*, SO = *Simosthenurus occidentalis*, MN = *Metasthenurus newtonae*, OR = *Osphranter rufus*, NR = *Notamacropus rufogriseus* (macropodiformes); GN = *Genyornis newtoni*, DN = *Dromaius novaehollandiae*, AL = *Alectura lathami* (large birds); TC = *Thylacoleo carnifex*, TH = *Thylacinus cynocephalus*, SH = *Sarcophilus harrisii*, DM = *Dasyurus maculatus* (carnivores); TA = *Tachyglossus aculeatus*, MR = *Megalibgwilia ramsayi* (monotreme invertivores). Here, we have depicted SH as ‘extant’, even though it went extinct on the mainland >3000 years ago.

Appendix 7

Comparing increasing mortality in juveniles and increasing offtake of juvenile individuals

Appendix 7—figure 1

Download asset Open asset

Appendix 8

Sensitivity analysis of extinction dates using a jack-knife approach based on the Gaussian-Resampled, Inverse-Weighted McInerny (GRIWM) Signor-Lipps correction method

To test the sensitivity of our conclusion to uncertainty in the dates of final extinction for each taxon taken from Saltré et al., 2016 and White et al., 2018, we adapted the Gaussian-resampled, inverse-weighted McInerny (GRIWM) approach (Bradshaw et al., 2012; Saltré et al., 2015) to derive new estimates and uncertainties of extinction time. Here, we jack-knifed the estimates 10,000 times for each taxon, calculating the median and 95% lower and upper confidence bounds. We applied this new approach to the following taxa: Diprotodon, Genyornis, Megalibgwilia, Metasthenurus, Palorchestes, Phascolonus, Procoptodon, Protemnodon, Sismosthenurus, Sthenurus, Thylacoleo, and Zygomaturus (but we kept the original dates for Sarcophilus and Thylacinus).

GRIWM does not depend on distributional assumptions given that it is a resampling approach (Saltré et al., 2015), even if it often provides wider uncertainty windows (especially for older dates) compared to other approaches. The main effect of this approach was to push back the extinction date for Metasthenurus (from a median of 50 ka to 89 ka), and push forward some of the dates within the other macropodiformes and vombatiforms (Appendix 8 figures 1–3). However, the net effect when we related these new estimates to any of the different risk metrics — mass (Appendix 8—figure 1a), generation length (Appendix 8—figure 1b), vital rate-specific extinction risk (Appendix 8—figure 2), or median susceptibility rank (Appendix 8—figure 3) — was to maintain our general conclusion of no relationship with the chronology.

Appendix 8—figure 1

Download asset Open asset

Relationship between estimated date of species extinction (across the entire continent) based on a jack-knifed GRIWM approach (Bradshaw et al., 2012; Saltré et al., 2015) and (a) body mass (kg) or (b) generation length (years) (Scenario LH).

Extinction-timing windows are estimated based on the agreement among six different models that correct for the Signor-Lipps effect (described in Materials and methods) in chronologies of quality-rated (Rodríguez-Rey et al., 2015) fossil dates for the studied taxa described in Peters et al., 2019. Here, we have depicted *Sarcophilus* as ‘extant’, even though it went extinct on the mainland >3000 years ago. Also shown are the approximate major climate periods and transitions: Marine Isotope Stage 5 (MIS 5), MIS 4, MIS 3, MIS 2 (including the Last Glacial Maximum, Antarctic Cold Reversal, and Younger Dryas), and the Holocene (including the period of sea level flooding when Tasmania separated from the mainland, and the relatively warm, wet, and climatically stable Holocene optimum).

Appendix 8—figure 2

Download asset Open asset

Relationship between estimated date of species extinction (across the entire mainland) based on a jack-knifed GRIWM approach (Bradshaw et al., 2012; Saltré et al., 2015) and area under the quasi-extinction curve (from Fig. S7) — extinction risk: ∫Pr(E_q)d(x) — for (a) (Scenario ↓S_j) decreasing juvenile survival, (b) (Scenario ↓F) decreasing fertility, (c) (Scenario ↓S) decreasing survival across all age classes, (d) (Scenario ↓ind) increasing number of individuals removed year⁻¹, (e) (Scenario ↑cat_f) increasing frequency of catastrophic die-offs per generation, and (f) (Scenario ↑cat_M) increasing magnitude of catastrophic die-offs.

Extinction-timing windows are estimated based on the agreement among six different models that correct for the Signor-Lipps effect (described in Materials and methods) in chronologies of quality-rated (Rodríguez-Rey et al., 2015) fossil dates for the studied taxa described in Peters et al., 2019.

Appendix 8—figure 3

Download asset Open asset

(a) Sum of the areas under the quasi-extinction curve for each of the six scenarios considered — total extinction risk: Σ∫Pr(E_q)d(x) — for each of the 21 modelled species (^†extinct; ^♀extant; scenario abbreviations: ↓S_j = reducing juvenile survival; ↓F = reducing fertility; ↓S = reducing survival; ind = reducing number of individuals; ↑cat_f = increasing frequency of catastrophe; ↑cat_M = increasing magnitude of catastrophe); (b) median susceptibility rank across the six scenarios considered (where higher ranks = higher susceptibility to extinction) for each species (red = extinct; black = extant; outline-only bars = macropodiformes; solid bars = vombatiforms; angled crosshatching = birds; vertical crosshatching = carnivores; brick crosshatching = monotremes); (c) median susceptibility rank as a function of log₁₀ generation length (G, kg) — there was a weak correlation including all species (solid grey line 1), but a strong relationship after removing *Sarcophilus* and *Dromaius* (dashed grey line 2) (information-theoretic evidence ratio [ER] and variance explained [R²] shown for each); (d) estimated date of species extinction (across entire continent) as a function of median susceptibility rank; taxonomic/functional groupings are indicated by colored symbols and convex hulls (macropodids: brown; monotremes: grey; vombatiforms: green; birds: blue; carnivores: red).

Extinction-timing windows are on a jack-knifed GRIWM approach (Bradshaw et al., 2012; Saltré et al., 2015) that corrects for the Signor-Lipps effect in chronologies of quality-rated (Rodríguez-Rey et al., 2015) fossil dates for the studied taxa described in Peters et al., 2019.

Appendix 9

Comparing demographic susceptibility to climate variation

Failing to observe any relationship between overall demographic susceptibility and the extinction chronology for Sahul, we might alternatively expect species’ demographic susceptibility to align with increasing environmental stress expressed as hotter, drier climates ( Saltré et al., 2019). We therefore hypothesized that the most extreme (hottest/driest) climates of the past would eventually drive the most-resilient species to extinction, which would manifest as a negative relationship between demographic susceptibility and warming/drying conditions (i.e., only when conditions became bad enough did the least-susceptible species succumb).

To test this hypothesis, we compiled climate indices hindcasted for the estimated extinction windows for the species we considered here. To this end, we acquired four hindcasted, continentally averaged climate variables from the intermediate-complexity, three-dimensional, Earth-system model known as LOVECLIM (Goosse et al., 2010; Timmermann and Friedrich, 2016). LOVECLIM hindcasts various climatic conditions by incorporating representations of the atmosphere, ocean and sea ice, land surface (including a vegetation submodel), ice sheets, and the carbon cycle. These variables were — mean annual temperature (°C), mean annual precipitation (mm), net primary production (kg C m⁻² year⁻¹), and fraction of the landscape designated as ‘desert’ — all expressed as anomalies of their respective average values calculated relative to 120 ka (i.e., a time when all species we considered were extant). We downscaled the original spatial resolution of LOVECLIM (5.625°×5.625°) to an output resolution of 1°×1° using bilinear interpolation because it retains the integrity and limitations of the original model output data.

We then calculated the information-theoretic evidence ratios (ER) for all relationships between the mean value of the climate variable across the entirely of Sahul and the sum of the extinction integrals across scenarios as the Akaike’s information criterion (AIC) of the slope model: y = α + βx divided by the AIC of the intercept-only model: y = α (i.e., ER_mean = AIC_slope/AIC_intercept). To incorporate full uncertainty in the climate variables (y), we developed a randomization test where we uniformly resampled the y values between y_min and y_max, estimating the residual sum of squares of the resampled values at each iteration compared to a randomized order of these residuals. We then calculated the probability (p_u) of producing a randomly generated relationship between the climate variable and the sum of the extinction integrals as the number of iterations when the randomized order produced a residual sum of squares ≤ the residual sum of squares of the resampled (ordered) climate variables divided by the total number of iterations (10,000).

When we plotted the four climate variables against the sum of the quasi-extinction integrals across scenarios, there was evidence for a weak, negative relationship between mean annual precipitation anomaly and relative extinction susceptibility (ER_mean = 9.6; Appendix 9—figure 1b), and a weak, positive relationship with desert-fraction anomaly (ER_mean = 4.2; Fig. Appendix 9—figure 1d). There was no evidence for a relationship between the means for temperature (Appendix 9—figure 1a) and net primary production (Appendix 9—figure 1c) anomalies and quasi-extinction integrals. The relationship with desert fraction supports the hypothesis that a drier (more desert-like) environment might have been related to extinction susceptibility. However, when we took full uncertainty of the climate variables into account in a randomization least-squares regression, none of the relationships could not be differentiated from a random process (p_u >0.17; Appendix 9—figure 1).

Appendix 9—figure 1

Download asset Open asset

Sum of the areas under the quasi-extinction curve over the six scenarios considered — total extinction risk: Σ∫Pr(E_q)d(x) — for each of the 13 extinct (mainland only) modelled species relative to (a) mean annual temperature anomaly (°C): information-theoretic evidence ratio of the slope model relative to the intercept-only model (ER_mean) = 0.75 for the mean climate values; probability of a non-random slope relationship incorporating full uncertainty in the climate variable p_u = 0.603.

(b) mean annual precipitation anomaly (mm): ER_mean = 9.6; p_u = 0.461. (c) net primary production anomaly (kg C m⁻² year⁻¹): ER_mean <0.01; p_u = 0.411. (d) desert fraction anomaly: ER_mean = 4.2; p_u = 0.425. The dashed red line in panels a, b, and d indicate evidence for a slope model versus the intercept-only model for these variables (ER_mean >2). Error bars indicate ±1 standard deviation.

The lack of evidence for a relationship between extinction susceptibility and warming/drying conditions for the mean climate conditions across the continent contrasts with recent evidence that water availability potentially exacerbated mortality from novel human hunting (Saltré et al., 2019). However, there is too much uncertainty in the climate hindcasts to test this hypothesis definitively. Another weakness of this approach is that we were obliged to take continental-scale averages of average climate conditions, which obviously ignores spatial complexity previously established as an important element in explaining the chronology and directionality of megafauna extinctions, at least in south-eastern Sahul (Saltré et al., 2019). Thus, this enticing, but still unsupported hypothesis that warming and drying conditions were related to intrinsic extinction susceptibility, will require more precise estimates of extinction timing and hindcasted climate conditions, and perhaps greater sample sizes across more species.

Data availability

All data and are R code needed to reproduce the analyses are available for download at http://github.com/cjabradshaw/MegafaunaSusceptibility.

References

(1987) The relationship between fecundity and adult body weight in homeotherms
Oecologia 73:478–480.

https://doi.org/10.1007/BF00385268
- PubMed
- Google Scholar
(2010) Highly skewed sex ratios and biased fossil deposition of moa: ancient DNA provides New insight on New Zealand's extinct megafauna
Quaternary Science Reviews 29:753–762.

https://doi.org/10.1016/j.quascirev.2009.11.022
- Google Scholar
1. Alonso-Llamazares C
2. Pablos A
(2019) Sex estimation from the calcaneus and talus using discriminant function analysis and its possible application in fossil remains
Archaeological and Anthropological Sciences 11:4927–4946.

https://doi.org/10.1007/s12520-019-00855-y
- Google Scholar
1. Alroy J
(2001) A multispecies overkill simulation of the end-Pleistocene megafaunal mass extinction
Science 292:1893–1896.

https://doi.org/10.1126/science.1059342
- PubMed
- Google Scholar
1. Asier L
(2016) Shoulder height, body mass, and shape of proboscideans
Acta Palaeontologica Polonica 61:537–574.

https://doi.org/10.4202/app.00136.2014
- Google Scholar
Data
1. Atlas of Living Australia
(authors) (2020) Dromaius novaehollandiae (Latham, 1790)
Atlas of Living Australia. ID 2714924-4fd5-456e-bb04-d23edbcf888f.

https://bie.ala.org.au/species/urn:lsid:biodiversity.org.au:afd.taxon:c2714924-4fd5-456e-bb04-d23edbcf888f
(2016) Robustness despite uncertainty: regional climate data reveal the dominant role of humans in explaining global extinctions of late Quaternary megafauna
Ecography 39:152–161.

https://doi.org/10.1111/ecog.01566
- Google Scholar
Book
(2008)
Dasyurus maculatus

In: van Dyck S, Strahan R, editors. The Mammals of Australia. New Holland. pp. 60–62.
- Google Scholar
1. Bode M
2. Brennan KEC
(2011) Using population viability analysis to guide research and conservation actions for Australia's threatened malleefowl Leipoa ocellata
Oryx 45:513–521.

https://doi.org/10.1017/S0030605311000688
- Google Scholar
(2012) Robust estimates of extinction time in the geological record
Quaternary Science Reviews 33:14–19.

https://doi.org/10.1016/j.quascirev.2011.11.021
- Google Scholar
(2013) More analytical bite in estimating targets for shark harvest
Marine Ecology Progress Series 488:221–232.

https://doi.org/10.3354/meps10375
- Google Scholar
1. Bradshaw CJA
2. Ehrlich PR
3. Beattie A
4. Ceballos G
5. Crist E
6. Diamond J
7. Dirzo R
8. Ehrlich AH
9. Harte J
10. Harte ME
11. Pyke G
12. Raven PH
13. Ripple WJ
14. Saltré F
15. Turnbull C
16. Wackernagel M
17. Blumstein DT
(2021) Underestimating the challenges of avoiding a ghastly future
Frontiers in Conservation Science 1:9.

https://doi.org/10.3389/fcosc.2020.615419
- Google Scholar
1. Bradshaw CJA
2. Brook BW
(2005) Disease and the devil: density-dependent epidemiological processes explain historical population fluctuations in the Tasmanian devil
Ecography 28:181–190.

https://doi.org/10.1111/j.0906-7590.2005.04088.x
- Google Scholar
(2006) Minimum viable population sizes and global extinction risk are unrelated
Ecology Letters 9:375–382.

https://doi.org/10.1111/j.1461-0248.2006.00883.x
- PubMed
- Google Scholar
(2008) Synergies among extinction drivers under global change
Trends in Ecology & Evolution 23:453–460.

https://doi.org/10.1016/j.tree.2008.03.011
- Google Scholar
1. Brook BW
2. Bowman DMJS
(2004) The uncertain blitzkrieg of Pleistocene megafauna
Journal of Biogeography 31:517–523.

https://doi.org/10.1046/j.1365-2699.2003.01028.x
- Google Scholar
1. Brook BW
2. Bowman DMJS
(2005) One equation fits overkill: why allometry underpins both prehistoric and modern body size-biased extinctions
Population Ecology 47:137–141.

https://doi.org/10.1007/s10144-005-0213-4
- Google Scholar
1. Brook BW
2. Johnson CN
(2006) Selective hunting of juveniles as a cause of the imperceptible overkill of the Australian Pleistocene megafauna
Alcheringa 30:39–48.

https://doi.org/10.1080/03115510609506854
- Google Scholar
(2011) Prey body size and ranking in zooarchaeology: theory, empirical evidence, and applications from the northern Great Basin
American Antiquity 76:403–428.

https://doi.org/10.7183/0002-7316.76.3.403
- Google Scholar
1. Cardillo M
2. Mace GM
3. Jones KE
4. Bielby J
5. Bininda-Emonds OR
6. Sechrest W
7. Orme CD
8. Purvis A
(2005) Multiple causes of high extinction risk in large mammal species
Science 309:1239–1241.

https://doi.org/10.1126/science.1116030
- PubMed
- Google Scholar
Book
1. Caswell H
(2001)
Matrix Population Models: Construction, Analysis, and Interpretation (2nd edn)

Sunderland, USA: Sinauer Associates, Inc.
- Google Scholar
1. Catt DC
(1977) The breeding biology of Bennett's wallaby (Macropus rufogriseus jruticus) in South Canterbury, New Zealand
New Zealand Journal of Zoology 4:401–411.

https://doi.org/10.1080/03014223.1977.9517965
- Google Scholar
Book
1. Caughley G
2. Gunn A
(1996)
Conservation Biology in Theory and Practice

Cambridge, USA: Blackwell Science.
- Google Scholar
(2005) Pleistocene to recent dietary shifts in California condors
PNAS 102:16707–16711.

https://doi.org/10.1073/pnas.0508529102
- PubMed
- Google Scholar
1. Chase JM
2. Blowes SA
3. Knight TM
4. Gerstner K
5. May F
(2020) Ecosystem decay exacerbates biodiversity loss with habitat loss
Nature 584:238–243.

https://doi.org/10.1038/s41586-020-2531-2
- PubMed
- Google Scholar
1. Choquenot D
2. Bowman DMJS
(1998) Marsupial megafauna, Aborigines and the overkill hypothesis: application of predator-prey models to the question of Pleistocene extinction in Australia
Global Ecology and Biogeography Letters 7:167–180.

https://doi.org/10.2307/2997372
- Google Scholar
Book
1. Cockburn A
(1997)
Living slow and dying young: senescence in marsupials

In: Saunders N, Hinds L, editors. Marsupial Biology: Recent Research, New Perspectives. Sydney: University of New South Wales Press. pp. 163–171.
- Google Scholar
Book
1. Corbett L
(1995)
The Dingo in Australia and Asia

Sydney: University of New South Wales Press.
- Google Scholar
1. Cork SJ
2. Dove H
(1989) Lactation in the tammar wallaby (Macropus eugenii). II. Intake of milk components and maternal allocation of energy
Journal of Zoology 219:399–409.

https://doi.org/10.1111/j.1469-7998.1989.tb02588.x
- Google Scholar
(2009) High diversity Pleistocene rainforest dasyurid assemblages with implications for the radiation of the Dasyuridae
Austral Ecology 34:663–669.

https://doi.org/10.1111/j.1442-9993.2009.01972.x
- Google Scholar
(2018) Short-tailed mice with a long fossil record: the genus Leggadina (Rodentia: Muridae) from the Quaternary of Queensland, Australia
PeerJ 6:e5639.

https://doi.org/10.7717/peerj.5639
- PubMed
- Google Scholar
(2017) High mortality and small population size prevent population recovery of a reintroduced mesopredator
Animal Conservation 20:555–563.

https://doi.org/10.1111/acv.12358
- Google Scholar
1. Damuth J
(1981) Population density and body size in mammals
Nature 290:699–700.

https://doi.org/10.1038/290699a0
- Google Scholar
(2009) Multiple ecological pathways to extinction in mammals
PNAS 106:10702–10705.

https://doi.org/10.1073/pnas.0901956106
- PubMed
- Google Scholar
(2007) An analysis of the relationship between metabolism, developmental schedules, and longevity using phylogenetic independent contrasts
Journal of Gerontology: Series A 62:149–160.

https://doi.org/10.1093/gerona/62.2.149
- Google Scholar
1. Dillingham PW
2. Moore JE
3. Fletcher D
4. Cortés E
5. Curtis KA
6. James KC
7. Lewison RL
(2016) Improved estimation of intrinsic growth r_max for long-lived species: integrating matrix models and allometry
Ecological Applications 26:322–333.

https://doi.org/10.1890/14-1990
- PubMed
- Google Scholar
1. Fagan WF
2. Holmes EE
(2006) Quantifying the extinction vortex
Ecology Letters 9:51–60.

https://doi.org/10.1111/j.1461-0248.2005.00845.x
- PubMed
- Google Scholar
(2001) The ecological basis of life history variation in marsupials
Ecology 82:3531–3540.

https://doi.org/10.1890/0012-9658(2001)082[3531:TEBOLH]2.0.CO;2
- Google Scholar
1. Frank M
2. Slaton A
3. Tinta T
4. Capaldi A
(2015) Investigating anthropogenic mammoth extinction with mathematical models
SPORA: A Journal of Biomathematics 1:8–16.

https://doi.org/10.30707/SPORA1.1Frank
- Google Scholar
(2014) Genetics in conservation management: revised recommendations for the 50/500 rules, Red List criteria and population viability analyses
Biological Conservation 170:56–63.

https://doi.org/10.1016/j.biocon.2013.12.036
- Google Scholar
1. Frith HJ
(1959) Breeding of the Mallee-fowl, Leipoa ocellata gould (Megapodiidae)
CSIRO Wildlife Research 4:31–60.

https://doi.org/10.1071/CWR9590031
- Google Scholar
1. Fritz SA
2. Purvis A
(2010) Selectivity in mammalian extinction risk and threat types: a new measure of phylogenetic signal strength in binary traits
Conservation Biology 24:1042–1051.

https://doi.org/10.1111/j.1523-1739.2010.01455.x
- PubMed
- Google Scholar
(2016) Economic utility and nutritional value of the common wombat (Vombatus ursinus): evaluating Australian Aboriginal hunting and butchery patterns
Journal of Archaeological Science: Reports 7:751–763.

https://doi.org/10.1016/j.jasrep.2015.10.037
- Google Scholar
1. Gittleman JL
2. Thompson SD
(1988) Energy allocation in mammalian reproduction
American Zoologist 28:863–875.

https://doi.org/10.1093/icb/28.3.863
- Google Scholar
1. Glen AS
(2008) Population attributes of the spotted-tailed quoll (Dasyurus maculatus) in north-eastern New South Wales
Australian Journal of Zoology 56:137–142.

https://doi.org/10.1071/ZO08025
- Google Scholar
1. Glen AS
2. Dickman CR
(2013) Population viability analysis shows spotted-tailed quolls may be vulnerable to competition
Australian Mammalogy 35:180–183.

https://doi.org/10.1071/AM12045
- Google Scholar
(2013) Which intrinsic traits predict vulnerability to extinction depends on the actual threatening processes
Ecosphere 4:1–16.

https://doi.org/10.1890/ES12-00380.1
- Google Scholar
1. Goosse H
2. Brovkin V
3. Fichefet T
4. Haarsma R
5. Huybrechts P
6. Jongma J
7. Mouchet A
8. Selten F
9. Barriat P-Y
10. Campin J-M
11. Deleersnijder E
12. Driesschaert E
13. Goelzer H
14. Janssens I
15. Loutre M-F
16. Morales Maqueda MA
17. Opsteegh T
18. Mathieu P-P
19. Munhoven G
20. Pettersson EJ
21. Renssen H
22. Roche DM
23. Schaeffer M
24. Tartinville B
25. Timmermann A
26. Weber SL
(2010) Description of the Earth system model of intermediate complexity LOVECLIM version 1.2
Geoscientific Model Development 3:603–633.

https://doi.org/10.5194/gmd-3-603-2010
- Google Scholar
(2019) Widespread male sex bias in mammal fossil and museum collections
PNAS 116:19019–19024.

https://doi.org/10.1073/pnas.1903275116
- PubMed
- Google Scholar
(2016) Is the “Genyornis” egg of a mihirung or another extinct bird from the Australian dreamtime?
Quaternary Science Reviews 133:147–164.

https://doi.org/10.1016/j.quascirev.2015.12.011
- Google Scholar
Book
1. Grzimek B
(1990)
Grzimek's Animal Life Encyclopedia

New York: McGraw-Hill Publishing Company.
- Google Scholar
1. Guiler ER
(1978)
Observations on the Tasmanian devil, Sarcophilus harrisii (Dasyuridae: Marsupiala) at Granville Harbour, 1966-75

Papers and Proceedings of the Royal Society of Tasmania 112:161–188.
- Google Scholar
1. Gurven M
2. Kaplan H
(2007) Longevity among hunter-gatherers: a cross-cultural examination
Population and Development Review 33:321–365.

https://doi.org/10.1111/j.1728-4457.2007.00171.x
- Google Scholar
1. Healy K
2. Guillerme T
3. Finlay S
4. Kane A
5. Kelly SBA
6. McClean D
7. Kelly DJ
8. Donohue I
9. Jackson AL
10. Cooper N
(2014) Ecology and mode-of-life explain lifespan variation in birds and mammals
Proceedings of the Royal Society B: Biological Sciences 281:20140298.

https://doi.org/10.1098/rspb.2014.0298
- Google Scholar
1. Hennemann WW
(1983) Relationship among body mass, metabolic rate and the intrinsic rate of natural increase in mammals
Oecologia 56:104–108.

https://doi.org/10.1007/BF00378224
- PubMed
- Google Scholar
(2006) Population characteristics predict responses in moose body mass to temporal variation in the environment
Journal of Animal Ecology 75:1110–1118.

https://doi.org/10.1111/j.1365-2656.2006.01138.x
- Google Scholar
1. Holz PH
2. Little PB
(1995) Degenerative leukoencephalopathy and myelopathy in dasyurids
Journal of Wildlife Diseases 31:509–513.

https://doi.org/10.7589/0090-3558-31.4.509
- PubMed
- Google Scholar
(2010) Estimates of maximum annual population growth rates (r_m) of mammals and their application in wildlife management
Journal of Applied Ecology 47:507–514.

https://doi.org/10.1111/j.1365-2664.2010.01812.x
- Google Scholar
(2020) Proximal humerus morphology indicates divergent patterns of locomotion in extinct giant kangaroos
Journal of Mammalian Evolution 27:627–647.

https://doi.org/10.1007/s10914-019-09494-5
- Google Scholar
1. Johnson CN
(2002) Determinants of loss of mammal species during the late Quaternary ‘megafauna’ extinctions: life history and ecology, but not body size
Proceedings of the Royal Society B: Biological Sciences 269:2221–2227.

https://doi.org/10.1098/rspb.2002.2130
- Google Scholar
Book
1. Johnson CN
(2006)
Australia's Mammal Extinctions: A 50 000 Year History

Cambridge, United Kingdom: Cambridge University Press.
- Google Scholar
1. Johnson CN
2. Alroy J
3. Beeton NJ
4. Bird MI
5. Brook BW
6. Cooper A
7. Gillespie R
8. Herrando-Pérez S
9. Jacobs Z
10. Miller GH
11. Prideaux GJ
12. Roberts RG
13. Rodríguez-Rey M
14. Saltré F
15. Turney CSM
16. Bradshaw CJA
(2016) What caused extinction of the Pleistocene megafauna of Sahul?
Proceedings of the Royal Society B: Biological Sciences 283:20152399.

https://doi.org/10.1098/rspb.2015.2399
- Google Scholar
1. Johnson CN
2. Prideaux GJ
(2004) Extinctions of herbivorous mammals in the late Pleistocene of Australia in relation to their feeding ecology: no evidence for environmental change as cause of extinction
Austral Ecology 29:553–557.

https://doi.org/10.1111/j.1442-9993.2004.01389.x
- Google Scholar
1. Jones DN
(1988) Hatching success of the Australian brush-turkey Alectura lathami in south-east Queensland
Emu 88:260–263.

https://doi.org/10.1071/MU9880260
- Google Scholar
Book
(1995)
The Megapodes: Megapodiidae, Bird Families of the World

Oxford: Oxford University Press.
- Google Scholar
1. Jones ME
2. Stoddart DM
(1998) Reconstruction of the predatory behaviour of the extinct marsupial thylacine (Thylacinus cynocephalus)
Journal of Zoology 246:239–246.

https://doi.org/10.1111/j.1469-7998.1998.tb00152.x
- Google Scholar
1. Jonzén N
2. Pople T
3. Knape J
4. Sköld M
(2010) Stochastic demography and population dynamics in the red kangaroo Macropus rufus
Journal of Animal Ecology 79:109–116.

https://doi.org/10.1111/j.1365-2656.2009.01601.x
- Google Scholar
1. Juanes F
(1986) Population density and body size in birds
The American Naturalist 128:921–929.

https://doi.org/10.1086/284613
- Google Scholar
(2009) Wild ostrich (Struthio camelus australis) reproduction in Orbata, a nature reserve in Tunisia
Tropical Animal Health and Production 41:1427–1438.

https://doi.org/10.1007/s11250-009-9331-x
- PubMed
- Google Scholar
1. Koch PL
2. Barnosky AD
(2006) Late Quaternary extinctions: state of the debate
Annual Review of Ecology, Evolution, and Systematics 37:215–250.

https://doi.org/10.1146/annurev.ecolsys.34.011802.132415
- Google Scholar
(2009) Demography, disease and the devil: life-history changes in a disease-affected population of Tasmanian devils (Sarcophilus harrisii)
Journal of Animal Ecology 78:427–436.

https://doi.org/10.1111/j.1365-2656.2008.01494.x
- Google Scholar
1. Latham ADM
2. Latham MC
3. Wilmshurst JM
4. Forsyth DM
5. Gormley AM
6. Pech RP
7. Perry GLW
8. Wood JR
(2020) A refined model of body mass and population density in flightless birds reconciles extreme bimodal population estimates for extinct moa
Ecography 43:353–364.

https://doi.org/10.1111/ecog.04917
- Google Scholar
1. Liow LH
2. Fortelius M
3. Bingham E
4. Lintulaakso K
5. Mannila H
6. Flynn L
7. Stenseth NC
(2008) Higher origination and extinction rates in larger mammals
PNAS 105:6097–6102.

https://doi.org/10.1073/pnas.0709763105
- PubMed
- Google Scholar
Preprint
(2020) Sahul’s megafauna were vulnerable to extinction due to their position in the trophic network
bioRxiv.

https://doi.org/10.1101/2021.01.19.427338
- Google Scholar
1. Lorenzen ED
2. Nogués-Bravo D
3. Orlando L
4. Weinstock J
5. Binladen J
6. Marske KA
7. Ugan A
8. Borregaard MK
9. Gilbert MT
10. Nielsen R
11. Ho SY
12. Goebel T
13. Graf KE
14. Byers D
15. Stenderup JT
16. Rasmussen M
17. Campos PF
18. Leonard JA
19. Koepfli KP
20. Froese D
21. Zazula G
22. Stafford TW
23. Aaris-Sørensen K
24. Batra P
25. Haywood AM
26. Singarayer JS
27. Valdes PJ
28. Boeskorov G
29. Burns JA
30. Davydov SP
31. Haile J
32. Jenkins DL
33. Kosintsev P
34. Kuznetsova T
35. Lai X
36. Martin LD
37. McDonald HG
38. Mol D
39. Meldgaard M
40. Munch K
41. Stephan E
42. Sablin M
43. Sommer RS
44. Sipko T
45. Scott E
46. Suchard MA
47. Tikhonov A
48. Willerslev R
49. Wayne RK
50. Cooper A
51. Hofreiter M
52. Sher A
53. Shapiro B
54. Rahbek C
55. Willerslev E
(2011) Species-specific responses of late Quaternary megafauna to climate and humans
Nature 479:359–364.

https://doi.org/10.1038/nature10574
- PubMed
- Google Scholar
1. Lyons SK
2. Miller JH
3. Fraser D
4. Smith FA
5. Boyer A
6. Lindsey E
7. Mychajliw AM
(2016) The changing role of mammal life histories in late Quaternary extinction vulnerability on continents and islands
Biology Letters 12:20160342.

https://doi.org/10.1098/rsbl.2016.0342
- PubMed
- Google Scholar
1. Marshall CR
(1997) Confidence intervals on stratigraphic ranges with nonrandom distributions of fossil horizons
Paleobiology 23:165–173.

https://doi.org/10.1017/S0094837300016766
- Google Scholar
1. McCarthy MA
(1998) Identifying declining and threatened species with museum data
Biological Conservation 83:9–17.

https://doi.org/10.1016/S0006-3207(97)00048-7
- Google Scholar
(2008) Allometric scaling and Bayesian priors for annual survival of birds and mammals
The American Naturalist 172:216–222.

https://doi.org/10.1086/588074
- Google Scholar
Book
1. McIlroy JC
(1996)
Common wombat

In: Strahan R, editors. The Australian Museum Complete Book of Australian Mammals. Sydney: Reed Books. pp. 205–206.
- Google Scholar
Book
1. McIlroy JC
(2008)
Vombatus ursinus

In: van Dyck S, Strahan R, editors. The Mammals of Australia. Sydney: Reed New Holland. pp. 207–208.
- Google Scholar
(2006) Significance of sighting rate in inferring extinction and threat
Conservation Biology 20:562–567.

https://doi.org/10.1111/j.1523-1739.2006.00377.x
- PubMed
- Google Scholar
1. Metcalf JL
2. Turney C
3. Barnett R
4. Martin F
5. Bray SC
6. Vilstrup JT
7. Orlando L
8. Salas-Gismondi R
9. Loponte D
10. Medina M
11. De Nigris M
12. Civalero T
13. Fernández PM
14. Gasco A
15. Duran V
16. Seymour KL
17. Otaola C
18. Gil A
19. Paunero R
20. Prevosti FJ
21. Bradshaw CJA
22. Wheeler JC
23. Borrero L
24. Austin JJ
25. Cooper A
(2016) Synergistic roles of climate warming and human occupation in patagonian megafaunal extinctions during the last deglaciation
Science Advances 2:e1501682.

https://doi.org/10.1126/sciadv.1501682
- PubMed
- Google Scholar
1. Miller GH
2. Fogel ML
3. Magee JW
4. Gagan MK
5. Clarke SJ
6. Johnson BJ
(2005) Ecosystem collapse in Pleistocene Australia and a human role in megafaunal extinction
Science 309:287–290.

https://doi.org/10.1126/science.1111288
- PubMed
- Google Scholar
1. Miller G
2. Magee J
3. Smith M
4. Spooner N
5. Baynes A
6. Lehman S
7. Fogel M
8. Johnston H
9. Williams D
10. Clark P
11. Florian C
12. Holst R
13. DeVogel S
(2016) Human predation contributed to the extinction of the Australian megafaunal bird Genyornis newtoni ∼47 ka
Nature Communications 7:10496.

https://doi.org/10.1038/ncomms10496
- PubMed
- Google Scholar
(2020) Dietary generalism accelerates arrival and persistence of coral-reef fishes in their novel ranges under climate change
Global Change Biology 26:5564–5573.

https://doi.org/10.1111/gcb.15221
- PubMed
- Google Scholar
(2019) Northern quoll persistence is most sensitive to survivorship of juveniles
Wildlife Research 46:165–175.

https://doi.org/10.1071/WR18010
- Google Scholar
1. Morris WF
2. Doak DF
(2004) Buffering of life histories against environmental stochasticity: accounting for a spurious correlation between the variabilities of vital rates and their contributions to fitness
The American Naturalist 163:579–590.

https://doi.org/10.1086/382550
- PubMed
- Google Scholar
1. Munn AJ
2. Dawson TJ
(2001) Thermoregulation in juvenile red kangaroos (Macropus rufus) after pouch exit: higher metabolism and evaporative water requirements
Physiological and Biochemical Zoology 74:917–927.

https://doi.org/10.1086/324568
- PubMed
- Google Scholar
1. Nicol S
2. Andersen NA
(2007) The life history of an egg-laying mammal, the echidna (Tachyglossus aculeatus)
Ecoscience 14:275–285.

https://doi.org/10.2980/1195-6860(2007)14[275:TLHOAE]2.0.CO;2
- Google Scholar
Book
1. Nicol SC
2. Morrow GE
(2012) Sex and seasonality: reproduction in the echidna (Tachyglossus aculeatus)
In: Ruf T, Bieber C, Arnold W, Millesi E, editors. Living in a Seasonal World: Thermoregulatory and Metabolic Adaptations. Berlin, Heidelberg: Springer Berlin Heidelberg. pp. 143–153.

https://doi.org/10.1007/978-3-642-28678-0_13
- Google Scholar
(2008) Climate change, humans, and the extinction of the woolly mammoth
PLOS Biology 6:e79.

https://doi.org/10.1371/journal.pbio.0060079
- PubMed
- Google Scholar
(2004) What are the best correlates of predicted extinction risk?
Biological Conservation 118:513–520.

https://doi.org/10.1016/j.biocon.2003.10.002
- Google Scholar
(2001) Semelparity in a large marsupial
Proceedings of the Royal Society B: Biological Sciences 268:407–411.

https://doi.org/10.1098/rspb.2000.1369
- Google Scholar
Book
1. Oftedal OT
(1985)
Pregnancy and lactation

In: Hudson R. J, White R. G, editors. Bioenergetics of Wild Herbivores. Boca Raton, Florida: CRC Press. pp. 215–238.
- Google Scholar
(1998) Reproductive cessation in female mammals
Nature 392:807–811.

https://doi.org/10.1038/33910
- PubMed
- Google Scholar
Book
1. Pardi MI
2. Smith FA
(2012) Paleoecology in an era of climate change: how the past can provide insights into the future
In: Louys J, editors. Paleontology in Ecology and Conservation. Berlin, Heidelber: Springer Berlin Heidelberg. pp. 93–116.

https://doi.org/10.1007/978-3-642-25038-5_6
- Google Scholar
(2009)
Behavior of emu bird (Dromaius novaehollandiae)

Veterinary World 2:439–440.
- Google Scholar
1. Peters KJ
2. Saltré F
3. Friedrich T
4. Jacobs Z
5. Wood R
6. McDowell M
7. Ulm S
8. Bradshaw CJA
(2019) FosSahul 2.0, an updated database for the late Quaternary fossil records of Sahul
Scientific Data 6:272.

https://doi.org/10.1038/s41597-019-0267-3
- PubMed
- Google Scholar
1. Prowse TAA
2. Johnson CN
3. Lacy RC
4. Bradshaw CJA
5. Pollak JP
6. Watts MJ
7. Brook BW
(2013) No need for disease: testing extinction hypotheses for the thylacine using multi-species metamodels
Journal of Animal Ecology 82:355–364.

https://doi.org/10.1111/1365-2656.12029
- Google Scholar
(2014) An ecological regime shift resulting from disrupted predator-prey interactions in Holocene Australia
Ecology 95:693–702.

https://doi.org/10.1890/13-0746.1
- PubMed
- Google Scholar
(2000) Predicting extinction risk in declining species
Proceedings of the Royal Society B: Biological Sciences 267:1947–1952.

https://doi.org/10.1098/rspb.2000.1234
- Google Scholar
(2003) The frequency and severity of catastrophic die-offs in vertebrates
Animal Conservation 6:109–114.

https://doi.org/10.1017/S1367943003003147
- Google Scholar
(2019) The extraordinary osteology and functional morphology of the limbs in Palorchestidae, a family of strange extinct marsupial giants
PLOS ONE 14:e0221824.

https://doi.org/10.1371/journal.pone.0221824
- PubMed
- Google Scholar
1. Ripple WJ
2. Van Valkenburgh B
(2010) Linking top-down forces to the Pleistocene megafaunal extinctions
BioScience 60:516–526.

https://doi.org/10.1525/bio.2010.60.7.7
- Google Scholar
1. Rismiller PD
2. McKelvey MW
(2000) Frequency of breeding and recruitment in the short-beaked echidna, Tachyglossus aculeatus
Journal of Mammalogy 81:1–17.

https://doi.org/10.1644/1545-1542(2000)081<0001:FOBARI>2.0.CO;2
- Google Scholar
1. Rodríguez-Rey M
2. Herrando-Pérez S
3. Gillespie R
4. Jacobs Z
5. Saltré F
6. Brook BW
7. Prideaux GJ
8. Roberts RG
9. Cooper A
10. Alroy J
11. Miller GH
12. Bird MI
13. Johnson CN
14. Beeton N
15. Turney CSM
16. Bradshaw CJA
(2015) Criteria for assessing the quality of middle Pleistocene to Holocene vertebrate fossil ages
Quaternary Geochronology 30:69–79.

https://doi.org/10.1016/j.quageo.2015.08.002
- Google Scholar
1. Rodríguez-Rey M
2. Herrando-Pérez S
3. Brook BW
4. Saltré F
5. Alroy J
6. Beeton N
7. Bird MI
8. Cooper A
9. Gillespie R
10. Jacobs Z
11. Johnson CN
12. Miller GH
13. Prideaux GJ
14. Roberts RG
15. Turney CS
16. Bradshaw CJA
(2016) A comprehensive database of quality-rated fossil ages for Sahul's Quaternary vertebrates
Scientific Data 3:160053.

https://doi.org/10.1038/sdata.2016.53
- PubMed
- Google Scholar
1. Roger E
2. Laffan SW
3. Ramp D
(2011) Road impacts a tipping point for wildlife populations in threatened landscapes
Population Ecology 53:215–227.

https://doi.org/10.1007/s10144-010-0209-6
- Google Scholar
1. Sales J
(2007) The emu (Dromaius novaehollandiae): a review of its biology and commercial products
Avian and Poultry Biology Reviews 18:1–20.

https://doi.org/10.3184/147020607X245048
- Google Scholar
(2015) Uncertainties in dating constrain model choice for inferring extinction time from fossil records
Quaternary Science Reviews 112:128–137.

https://doi.org/10.1016/j.quascirev.2015.01.022
- Google Scholar
1. Saltré F
2. Rodríguez-Rey M
3. Brook BW
4. Johnson CN
5. Turney CS
6. Alroy J
7. Cooper A
8. Beeton N
9. Bird MI
10. Fordham DA
11. Gillespie R
12. Herrando-Pérez S
13. Jacobs Z
14. Miller GH
15. Nogués-Bravo D
16. Prideaux GJ
17. Roberts RG
18. Bradshaw CJA
(2016) Climate change not to blame for late Quaternary megafauna extinctions in Australia
Nature Communications 7:10511.

https://doi.org/10.1038/ncomms10511
- PubMed
- Google Scholar
1. Saltré F
2. Chadoeuf J
3. Peters KJ
4. McDowell MC
5. Friedrich T
6. Timmermann A
7. Ulm S
8. Bradshaw CJA
(2019) Climate-human interaction associated with southeast Australian megafauna extinction patterns
Nature Communications 10:5311.

https://doi.org/10.1038/s41467-019-13277-0
- PubMed
- Google Scholar
(2014) Global late Quaternary megafauna extinctions linked to humans, not climate change
Proceedings of the Royal Society B: Biological Sciences 281:20133254.

https://doi.org/10.1098/rspb.2013.3254
- Google Scholar
1. Saran KA
2. Parker G
3. Parker R
4. Dickman CR
(2011) Rehabilitation as a conservation tool: a case study using the common wombat
Pacific Conservation Biology 17:310–319.

https://doi.org/10.1071/PC110310
- Google Scholar
(2015) The effects of climate change and land-use change on demographic rates and population viability
Biological Reviews 90:837–853.

https://doi.org/10.1111/brv.12136
- PubMed
- Google Scholar
(2017) Taxonomic review of the late Cenozoic megapodes (Galliformes: Megapodiidae) of Australia
Royal Society Open Science 4:170233.

https://doi.org/10.1098/rsos.170233
- PubMed
- Google Scholar
Book
(2009) Causes and consequences of species extinctions
In: Levin S. A, editors. The Princeton Guide to Ecology. Princeton, USA: Princeton University Press. pp. 514–520.

https://doi.org/10.1515/9781400833023.514
- Google Scholar
(2006) On the Pleistocene extinctions of Alaskan mammoths and horses
PNAS 103:7351–7353.

https://doi.org/10.1073/pnas.0509480103
- PubMed
- Google Scholar
(2019) The limits to population density in birds and mammals
Ecology Letters 22:654–663.

https://doi.org/10.1111/ele.13227
- PubMed
- Google Scholar
Book
1. Strahan R
(1991)
The Australian Museum Complete Book of Australian Mammals

North Ryde, New South Wales: Cornstalk Publications.
- Google Scholar
1. Strauss D
2. Sadler PM
(1989) Classical confidence intervals and Bayesian probability estimates for ends of local taxon ranges
Mathematical Geology 21:411–427.

https://doi.org/10.1007/BF00897326
- Google Scholar
1. Strona G
2. Bradshaw CJA
(2018) Co-extinctions annihilate planetary life during extreme environmental change
Scientific Reports 8:16724.

https://doi.org/10.1038/s41598-018-35068-1
- PubMed
- Google Scholar
1. Timmermann A
(2020) Quantifying the potential causes of Neanderthal extinction: abrupt climate change versus competition and interbreeding
Quaternary Science Reviews 238:106331.

https://doi.org/10.1016/j.quascirev.2020.106331
- Google Scholar
1. Timmermann A
2. Friedrich T
(2016) Late Pleistocene climate drivers of early human migration
Nature 538:92–95.

https://doi.org/10.1038/nature19365
- PubMed
- Google Scholar
1. Tomiya S
(2013) Body size and extinction risk in terrestrial mammals above the species level
The American Naturalist 182:E196–E214.

https://doi.org/10.1086/673489
- Google Scholar
(2010) Pragmatic population viability targets in a rapidly changing world
Biological Conservation 143:28–34.

https://doi.org/10.1016/j.biocon.2009.09.001
- Google Scholar
1. Turbill C
2. Ruf T
(2010) Senescence is more important in the natural lives of long- than short-lived mammals
PLOS ONE 5:e12019.

https://doi.org/10.1371/journal.pone.0012019
- PubMed
- Google Scholar
Book
1. Tyndale-Biscoe CH
(2005)
Life of Marsupials

Clayton, Victoria: CSIRO Publishing.
- Google Scholar
(2016) Combination of humans, climate, and vegetation change triggered late Quaternary megafauna extinction in the Última Esperanza region, southern Patagonia, Chile
Ecography 39:125–140.

https://doi.org/10.1111/ecog.01606
- Google Scholar
1. Wan X
2. Zhang Z
(2017) Climate warming and humans played different roles in triggering late Quaternary extinctions in east and west Eurasia
Proceedings of the Royal Society B: Biological Sciences 284:20162438.

https://doi.org/10.1098/rspb.2016.2438
- Google Scholar
1. Way KC
(1988) Captive management and breeding of the tiger quoll: Dasyurus maculatus
International Zoo Yearbook 27:108–119.

https://doi.org/10.1111/j.1748-1090.1988.tb03203.x
- Google Scholar
1. Weisbecker V
2. Goswami A
(2010) Brain size, life history, and metabolism at the marsupial/placental dichotomy
PNAS 107:16216–16221.

https://doi.org/10.1073/pnas.0906486107
- PubMed
- Google Scholar
(2018) High-quality fossil dates support a synchronous, late Holocene extinction of devils and thylacines in mainland Australia
Biology Letters 14:20170642.

https://doi.org/10.1098/rsbl.2017.0642
- PubMed
- Google Scholar
1. Willis KJ
2. Birks HJ
(2006) What is natural? the need for a long-term perspective in biodiversity conservation
Science 314:1261–1265.

https://doi.org/10.1126/science.1122667
- PubMed
- Google Scholar
Book
1. Wilson DE
2. Mittermeier RA
(2001)
Handbook of the Mammals of the World, 2

Barcelona: Lynx Edicions.
- Google Scholar
(2017) The role of paleoecology in restoration and resource management—the past as a guide to future decision-making: review and example from the Greater Everglades Ecosystem, U.S.A
Frontiers in Ecology and Evolution 5:11.

https://doi.org/10.3389/fevo.2017.00011
- Google Scholar
1. Wroe S
2. Crowther M
3. Dortch J
4. Chong J
(2004) The size of the largest marsupial and why it matters
Proceedings of the Royal Society B: Biological Sciences 271:S34–S36.

https://doi.org/10.1098/rsbl.2003.0095
- Google Scholar

Article and author information

Author details

Corey JA Bradshaw
1. Global Ecology Partuyarta Ngadluku Wardli Kuu, College of Science and Engineering, Flinders University, Tarndanya (Adelaide), Australia
2. ARC Centre of Excellence for Australian Biodiversity and Heritage, Wollongong, Australia
Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Funding acquisition, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing, C.J.A.B conceptualized and wrote the paper, designed the research, wrote the code, sourced the data, and did the analysis

For correspondence
corey.bradshaw@flinders.edu.au

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5328-7741
Christopher N Johnson
1. ARC Centre of Excellence for Australian Biodiversity and Heritage, Wollongong, Australia
2. Dynamics of Eco-Evolutionary Pattern, University of Tasmania, Hobart, Australia
Contribution
Conceptualization, Resources, Methodology, Writing - review and editing, C.J. provided additional analytical approaches and data, and contributed to writing the manuscript

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-9719-3771
John Llewelyn
1. Global Ecology Partuyarta Ngadluku Wardli Kuu, College of Science and Engineering, Flinders University, Tarndanya (Adelaide), Australia
2. ARC Centre of Excellence for Australian Biodiversity and Heritage, Wollongong, Australia
Contribution
Conceptualization, Resources, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5379-5631
Vera Weisbecker
1. ARC Centre of Excellence for Australian Biodiversity and Heritage, Wollongong, Australia
2. College of Science and Engineering, Flinders University, Adelaide, Australia
Contribution
Conceptualization, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2370-4046
Giovanni Strona

Faculty of Biological and Environmental Sciences, University of Helsinki, Helsinki, Finland

Contribution
Conceptualization, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2294-4013
Frédérik Saltré
1. Global Ecology Partuyarta Ngadluku Wardli Kuu, College of Science and Engineering, Flinders University, Tarndanya (Adelaide), Australia
2. ARC Centre of Excellence for Australian Biodiversity and Heritage, Wollongong, Australia
Contribution
Conceptualization, Resources, Software, Formal analysis, Methodology, Writing - review and editing, Conceptualised and assisted writing the paper, and provided climate data

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5040-3911

Funding

Australian Research Council Centre of Excellence for Australian Biodiversity and Heritage (CE170100015)

Corey JA Bradshaw
Christopher N Johnson
Vera Weisbecker

Australian Research Council (DP170103227)

Vera Weisbecker

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This study was supported by the Australian Research Council through a Centre of Excellence grant (CE170100015) to CJAB, CNJ, and VW, and a Discovery Project grant (DP170103227) to VW. We acknowledge the Indigenous Traditional Owners of the land on which Flinders University is built — the Kaurna people of the Adelaide Plains.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

2,183

views
272

downloads
14

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Corey JA Bradshaw
Christopher N Johnson
John Llewelyn
Vera Weisbecker
Giovanni Strona
Frédérik Saltré

(2021)

Relative demographic susceptibility does not explain the extinction chronology of Sahul’s megafauna

eLife 10:e63870.

https://doi.org/10.7554/eLife.63870

Share this article

Cite this article

Description of five dominant mechanisms by which megafauna could have been driven to extinction and the associated seven perturbation scenarios examined.

Relationship between estimated date of species extinction (across the entire continent) and (a) body mass (M, kg) or (b) generation length (G, years) (Scenario LH).

Inset: Increasing extinction risk for birds — quasi-extinction probability: Pr(Eq) — as a function of increasing the mean number of eggs harvested per female per year (Scenario ↓Fe).

Extinction susceptibility of the 21 modelled species.

Histogram of log10 adult body masses (in kg) for the 21 species examined.

Predicted demographic values for each species (equations provided in Materials and methods and this appendix).

Kendall’s rank correlation (τ) among demographic values given in Appendix 2—table 1 across species.

Mass-predicted abundance and biomass for the 21 modelled species.

Relationship between the logarithm of adult female body mass (M, kg) and inter-birth interval (Ib, in days) for 23 extant macropodid herbivores (Fisher et al., 2001).

Relationship between the predicted age (years) at first breeding (αpred) and observed α (αobs) for 23 extant macropodid herbivores (Fisher et al., 2001).

Negative relationship between the logarithm of the maximum rate of intrinsic population growth (log10rm) and the logarithm of the age at primiparity (log10α) for the 21 species examined.

Increasing probabilities of quasi-extinction — quasi-extinction probability: Pr(Eq) — as a function of (a) decreasing juvenile survival (Scenario ↓Sj), and (b) increasing number of juvenile individuals removed year−1 (Scenario ↓indj).

Relationship between estimated date of species extinction (across the entire continent) based on a jack-knifed GRIWM approach (Bradshaw et al., 2012; Saltré et al., 2015) and (a) body mass (kg) or (b) generation length (years) (Scenario LH).

Author details

Corey JA Bradshaw

Contribution

For correspondence

Competing interests

Christopher N Johnson

Contribution

Competing interests

John Llewelyn

Contribution

Competing interests

Vera Weisbecker

Contribution

Competing interests

Giovanni Strona

Contribution

Competing interests

Frédérik Saltré

Contribution

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading

Inset: Increasing extinction risk for birds — quasi-extinction probability: Pr(E_q) — as a function of increasing the mean number of eggs harvested per female per year (Scenario ↓F_e).

Histogram of log₁₀ adult body masses (in kg) for the 21 species examined.

Relationship between the logarithm of adult female body mass (M, kg) and inter-birth interval (I_b, in days) for 23 extant macropodid herbivores (Fisher et al., 2001).

Relationship between the predicted age (years) at first breeding (α_pred) and observed α (α_obs) for 23 extant macropodid herbivores (Fisher et al., 2001).

Negative relationship between the logarithm of the maximum rate of intrinsic population growth (log₁₀r_m) and the logarithm of the age at primiparity (log₁₀α) for the 21 species examined.

Increasing probabilities of quasi-extinction — quasi-extinction probability: Pr(E_q) — as a function of (a) decreasing juvenile survival (Scenario ↓S_j), and (b) increasing number of juvenile individuals removed year⁻¹ (Scenario ↓ind_j).