Here, we first use a classical, single‐species PVA to model the fate of the thylacine retrospectively, incorporating information on habitat loss and bounty hunting, and we investigate how including a disease component modifies the model output. We then consider whether a more complete suite of documented, European‐imposed pressures can reasonably account for the thylacine's extinction. To this end, we link multiple PVA models dynamically to develop a predator–herbivore–vegetation metamodel representing the pre‐European trophic system. We then perturb this metamodel to simulate the impact of Europeans, including habitat modification, bounty hunting, macropod (prey) harvesting and the rapid growth of the sheep population that competed with macropods for food. Accepting that important parameters and interspecific interactions cannot be estimated directly, we subjected both models to detailed sensitivity analyses to investigate important assumptions and uncertainties. We hypothesised that the thylacine's extinction could be explained by interactions between known historical stressors, without invoking disease.

Interactions between the thylacine harvest and other European‐imposed pressures might account for the thylacine's extinction, without the need to invoke a hypothetical disease. In particular, the impact of European settlers on the thylacine's macropod prey base was probably important (Paddle 2000 ). European settlers attacked, displaced and transplanted the Aboriginal hunter‐gatherer population and radically altered macropod harvesting regimes through hunting for food and the fur trade (Crowther 1926 ; Barker & Caughley 1990 ). They also introduced non‐native herbivores – by 1951, more than two million sheep inhabited the grasslands favoured by the thylacine's native prey species (Davidson 1938 ; Kirkpatrick 2007 ).

While it is generally accepted that bounty hunting contributed to the thylacine's decline, it has been argued repeatedly that European impacts were insufficient to account for this extinction and that disease must also have been involved (Guiler 1985 ; Paddle 2000 ; IUCN 2011 ). Bulte, Horan & Shogren ( 2003 ) formalised this argument with a deterministic, density‐dependent population model of thylacine‐bounty hunter dynamics, which predicted persistence of the thylacine under the historical bounty harvest. Bounty records show a sudden reduction in annual thylacine kills after 1905, which could be interpreted as the impact of an epidemic (Guiler 1961 , 1985 ; McCallum & Dobson 1995 ; Bowman 2001 ) (see Figure S1 in the electronic supplementary material). However, the only evidence for disease is anecdotal and was recorded long after the likely time of the epidemic (Guiler 1961 ).

Australia has a particularly woeful record of mammal extinctions, having lost 22 mammal species over the last 200 years (Burbidge et al . 2008 ). The most infamous of these was the loss of the wolf‐like thylacine Thylacinus cynocephalus in Tasmania. Before European settlement in 1803, thylacines were sparsely but broadly distributed throughout most of Tasmania (Guiler 1985 ; Paddle 2000 ). In 1886, the Tasmanian government legislated a bounty payment for thylacines, paying bounties for 2 184 thylacine carcasses until the scheme's termination in 1909, although the true harvest number might be up to double this figure (Guiler & Godard 1998 , p. 128). Only a handful of animals were located after the bounty was lifted and, in 1933, the last known thylacine was captured from the wild (Paddle 2000 ).

Recent software development now allows multiple PVA models constructed using the vortex simulation package (Lacy 2000 ) to be integrated within a single ‘metamodel’, thereby extending PVA to incorporate complex, multi‐species systems. Here, we report on the first case study using this approach. A metamodel is a network of component models that share information and represent different, often discipline‐specific elements of a complete system (Miller & Lacy 2003 ; Nyhus et al . 2007 ). PVA‐based metamodels are ideally suited to the investigation of extinction for two reasons. First, they allow the population dynamics of many species to be modelled simultaneously and permit the flexible specification of the functional links between them, while subjecting each component to the stochastic processes that are so important in the ‘extinction vortex’ (Gilpin & Soulé 1986 ). Second, because each component of a metamodel is compartmentalised, critical extinction drivers can be examined by modifying each component separately, or even investigated ‘experimentally’ by removing components from the system.

Mathematical modelling is widely used to assess the extinction risk of threatened species and guide management for recovery. Models of extinction processes must account for stochastic demographic, environmental and genetic processes and are therefore constructed using the synthetic tools of population viability analysis (PVA). PVA is used to evaluate different management interventions, using life‐history information for a single species or population and stochastic computer simulations. Generic PVA software has become increasingly sophisticated and allows incorporation of processes such as harvesting/supplementation regimes and spatially explicit habitat dynamics (Lacy 2000 ; Akçakaya 2009). Nevertheless, classical PVA models continue to operate in a ‘single‐species vacuum’ (Sabo 2008 ) because species interactions are either ignored or oversimplified (e.g. treated as random factors influencing the demographic rates of a focal species). This is a serious deficiency, given that many recent extinctions appear to have been caused by complex interactions among stressors (Brook, Sodhi & Bradshaw 2008 ; IUCN 2011 ), from which the primary extinction drivers can be difficult to isolate.

Given the complexity of these metamodels, we eschewed testing specific scenarios in favour of a broad sensitivity analysis based on Latin hypercube sampling as detailed above (Table 1 ). In addition to parameter sensitivities examined previously for the single‐species approach and detailed above for macropod harvesting by humans, we accounted for uncertainty in macropod densities and numerical responses. Sheep consume approximately 1·4–2·5 times more vegetation than kangaroos (Grigg 2002 ; Munn et al . 2008 ; Munn, Dawson & McLeod 2010 ), so we linked vegetation consumption by sheep to that by kangaroos through a competition coefficient ( c ) controlling this ratio. We varied the thylacine functional response so that per capita predation rates ranged from approximately 0·02–1 macropod killed per day. We assumed that, like wolves, thylacine vital rates were negatively affected when daily prey consumption fell below that required to maintain five times their basal metabolic rate (BMR) (Peterson & Ciucci 2003 ). Using an adult body mass for thylacines of 20 kg (Jones 1997 ) and the procedure of Peterson & Ciucci ( 2003 ), we estimated that 5 times thylacine BMR could be met by 2·15 kg of captured prey daily. As we only modelled a proportion of the thylacine's prey base, however, for sensitivity analysis we assumed the predation rate required for a stable thylacine population ranged from 2·15 kg prey indiv −1 d −1 (i.e. metabolic requirements could be met by consumption of modelled prey only) down to 0·5 kg prey indiv −1 d −1 (i.e. >75% of metabolic requirements could be met by unmodelled prey items). Because the thylacine population size was not set in advance, but emerged dynamically, we burned in metamodels for 500 years prior to European settlement. Simulations were retained if the population mean over this burn‐in ranged between 2 000 and 4 000 individuals, and model output was analysed using BRT as above.

Metamodels were constructed using the command‐centre package metamodel manager, a flexible programme that can link component vortex models and manage the sequence in which these components are invoked as well as the flow of information between them (Miller & Lacy 2003 ; Raboy 2011 ; Bradshaw et al . 2012 ). vortex models can store and manipulate state variables that are accessible by metamodel manager and can be passed on for use and/or manipulation by other models or by mmmacro, a macro (script) interpreter that can implement complex variable operations if required (for a detailed treatment of this approach, see Bradshaw et al . 2012 ). In our metamodels, all demographic processes (e.g. mortality, fertility, inbreeding depression and Allee effects) were implemented by component vortex models, which then stored updated population sizes as state variables at the end of each year simulated. After all component models had completed operations for 1 year, we passed these state variables to mmmacro, which implemented the model of vegetation growth and grazing pressure, calculated predation rates and stored required information as state variables also. Predation of macropods by thylacines was then implemented by component vortex models before demographic processes were executed for the following year.

We also modelled predator–prey dynamics using a fully interactive approach. Each year, total prey abundance was calculated after pooling across both habitats; thylacines consumed prey at a rate determined by their functional response and prey abundances diminished accordingly. We used linear, ratio‐dependent functional responses for thylacines governed by two parameters: the attack rate ( a ) on macropods and a prey‐switching coefficient ( b ) (see Appendix S2 for details). We specified the numerical response of thylacines as a logistic function of their per capita rate of prey consumption, again mediated through mortality‐rate adjustments.

The Forester kangaroo and Bennett's wallaby select open grassland and woody grassland habitats (Tanner & Hocking 2001 ; Le Mar & McArthur 2005 ; Wiggins et al . 2010 ), hereafter collectively termed ‘grassland’. Kirkpatrick ( 2007 ) mapped the pre‐European extent of Tasmanian grassland; after geo‐referencing, we estimated that this habitat covered 7 455 km 2 (or 13%) of the thylacine's original range. The remainder of the thylacine's distribution was dominated by poor‐quality habitats for macropods (Wapstra 1976 ) and is termed ‘forest’ for convenience, although it incorporates a range of nongrassland habitat types (Kirkpatrick & Dickinson 1984 ). We modelled the population dynamics of macropod herbivores in grassland and forest habitats using different approaches. In grasslands, we used a modified version of Caughley's ( 1987 ) interactive model of macropod–vegetation dynamics. In brief, we specified the numerical response of kangaroos and wallabies (mediated through age‐structured vital rates) as a function of vegetation biomass ( V ), which in turn was determined by simulated rainfall, past vegetation biomass and herbivore density. The interactive model allowed us to simulate the impact of sheep on macropod abundances following European settlement. In forest, we assumed an equilibrium density for each macropod species that was modified annually using a rainfall multiplier, thereby enforcing reduced (increased) equilibrium densities in below‐average (above‐average) rainfall years (see Appendix S2 for details).

We constructed a predator–herbivore–vegetation metamodel to recreate the dynamic interaction of thylacines and their macropod prey before the arrival of Europeans. The metamodel architecture consisted of PVA models for individual species constructed in vortex, linked dynamically through functional responses and underpinned by a model of vegetation growth (see Appendix S1 and S2 in the electronic supplementary material for full details of the vortex model parameterisation and metamodel linkages, respectively). Historical accounts indicate that thylacines preyed heavily on two large macropod herbivores, the Forester kangaroo Macropus giganteus tasmaniensis and Bennett's wallaby Macropus rufogriseus (Paddle 2000 ; p. 37, 45, and 79; Cosgrove & Allen 2001 ). Although the thylacine's reliance on large prey is not universally accepted (Jones & Stoddart 1998 ), these accounts are corroborated by aspects of the animal's digestive system: specifically, its alimentary canal was short and its stomach was muscular and capable of great distension (Crisp 1855 ). This suggests a feeding habit involving the rapid input of large amounts of food, with prey kills occurring perhaps once every three days (Guiler 1985 ). We assumed that the thylacine's prey base could be represented by Forester kangaroos and Bennett's wallabies, while accounting for the fact that thylacines undoubtedly consumed smaller prey also, through sensitivity analysis (see below).

For sensitivity analysis, we used the R package lhs (Carnell 2009 ) to construct a 9‐dimension Latin hypercube with 10 000 divisions (Table 1 ) by sampling from (i) three binary factors representing the reproductive mode of thylacines and the inclusion/exclusion of inbreeding depression/Allee effects; and (ii) wide, but plausible, ranges for six continuous parameters (including parameters governing simulated harvest numbers and disease effects). To account for uncertainty in thylacine vital rates, and because the strength of density dependence is often critically important for PVA models (Henle, Sarre & Wiegand 2004 ), we tested a range of maximum and minimum population growth rates permissible by the deterministic, mortality‐rate function. A single PVA simulation was done for each hypercube division and scored as either a ‘success’ (i.e. thylacines were driven extinct) or ‘failure’. To evaluate the relative importance of parameter inputs to simulation success, and to simplify the results of this sensitivity analysis for presentation, we analysed this output using boosted regression trees (BRT) (Elith, Leathwick & Hastie 2008 ) implemented via functions in the R package dismo (Hijmans et al . 2011 ). All BRT models were fitted with the routine gbm.step , using a binomial error and logit link, a learning rate of 0·01 and a bag fraction of 0·75. We simplified BRT models by examining metrics of relative influence for each factor/parameter and comparing misclassification rates between models based on 10‐fold cross‐validation.

All single‐species PVA models ran from 1753 ad (50 years before European settlement) to 1935 (2 years after the last thylacine died in captivity). We conducted 1 000 simulations of each PVA model. We defined extinctions as simulations in which the thylacine population dropped to 50 individuals or fewer and generated cumulative probability of extinction curves with 95% confidence intervals calculated using the Kolmogorov–Smirnov test statistic (Sokal & Rohlf 1981 ).

We then finalised PVA models for each starting population size according to either a baseline or pessimistic scenario. The baseline scenario assumed a polygynous breeding system and that the bounty records report the true numbers of thylacines harvested. The pessimistic scenario assumed a short‐term, monogamous breeding system and twice the reported thylacine harvest (Paddle 2000 ). The latter scenario also included inbreeding depression, assuming 3·14 lethal equivalents per individual of which 50% was due to recessive lethal alleles (Miller & Lacy 2005 ), and a simple demographic Allee effect that, once thylacine density fell below 10% of initial density, reduced the mean percentage of breeding females linearly to a minimum of zero at a population density of zero. Both scenarios included moderate environmental variation (5% standard deviation) in vital rates, in addition to the demographic stochasticity inherent in the vortex model structure.

We constructed age‐structured, stochastic demographic models for thylacines using(Lacy). We derived demographic parameters from published values or obtained them from similar species (full details are provided in Appendix S1 of the electronic supplementary material). Briefly, we derived age‐structured fertility rates for thylacines from data for the Tasmanian Devil, a related marsupial carnivore (also Order Dasyuromorphia) that similarly has four teats (Lachish, McCallum & Jones). Initial, age‐structured mortality rates were derived from a cursorial predator with similar longevity, the spotted hyena (Watts & Holekamp). The deterministic component of the models employed a density‐dependent mortality function that allowed improvements in thylacine vital rates at low densities. We estimated a maximum population growth rate () for thylacines of 0·3405 (Appendix S1). We assumed that mortality rates for a given densityand age, denoted by μ, were reduced along a logistic curve from an upper limit when density is high (μ) to a minimum when density is low (μ) according to the formula:wherecontrols the‐intercept (i.e. μ) andcontrols the slope of the mortality curve. We set (μ) by proportionally reducing the initial age‐specific mortality rates until the deterministic population growth rate () equalled our estimate of. We fit the slope parameterso that the deterministic population growth ratewas zero when thylacine density equalled the initial density.

Guiler & Godard ( 1998 ) estimated the pre‐European thylacine population at 2000–4000 individuals after considering habitat‐specific differences in thylacine abundance, based on home range area estimates derived from the Woolnorth sheep station records and from the spatial pattern of government bounty payments. We also used Kelt & Van Vuren's ( 2001 ) allometric relationship between body mass and home range area for carnivorous mammals, and assumed an adult body mass for thylacines of 20 kg (Jones 1997 ), to estimate a home range area for thylacines of 25·7 km 2 . Assuming nonoverlapping home ranges, this equates to 2 184 individuals across the thylacine's distribution, which is close to the lower limit of Guiler & Godard's ( 1998 ) estimated range. However, because thylacine home ranges could have overlapped, we tested PVA models with five starting population sizes spanning their full range of 2 000–4 000 individuals.

We estimated an initial range of 56 051 km 2 for thylacines after georeferencing Paddle's ( 2000 ) reconstruction of the species' pre‐European distribution. We estimated progressive habitat loss from maps documenting Tasmanian land alienation (land grants) by Europeans (Davies 1965 ). To afford single‐species models the greatest opportunity of recreating the thylacine's extinction, we assumed all alienated land was lost from the modelled habitat area (Table 1 ), resulting in a range reduction for thylacines of 46% by 1935. We took annual harvest numbers for thylacines over 1886–1909 directly from the government bounty records (Guiler 1985 ). Disease‐induced increases in mortality are difficult to quantify in the wild but probably range up to 40% for common mammalian diseases such as rabies and distemper (Vucetich & Creel 1999 ). We simulated sensitivity to a disease outbreak in the thylacine population from 1906 to 1909 by increasing annual, age‐structured mortality rates over this period by up to 40%, thereby emulating plausible disease‐induced reductions in fitness. Paddle ( 2000 , pp. 202–204) suggested, based on records of sick thylacines in zoos, that epidemic disease could have affected the population for a longer period; however, it is impossible to extrapolate diseases that occurred in captivity to the Tasmanian wilderness. We therefore initiated simulated disease outbreaks in 1906, the first year in which government bounty payments decreased dramatically, and terminated them before 1910, the year by which thylacines are deemed to have been functionally extinct (Figure S1).

A BRT model including all metamodel parameters as predictors (Table 2 ) and five splits yielded a misclassification rate of 17·5%. Parameters governing macropod prey dynamics and the strength of the human harvest exerted little influence on the probability of simulation success. Again, a simplified model including the top five predictors and no interactions did not compromise the misclassification rate (an increase of 2·0% only) (Fig. 4 a). Again the thylacine starting population size, maximum population growth rate and harvest multiplier were influential predictors, as was the sheep competition coefficient. Even when only moderate competition from sheep was assumed ( c = 2), metamodels that did not include disease were able to simulate the thylacine extinction as readily as single‐species models including a virulent disease (cf. Figs 2 b & 4 b).

Example simulation for a metamodel with a pre‐European thylacine population of 3 985 individuals (averaged over the 500 years burn‐in). Population size for (a) thylacines, (b) the kangaroo Macropus giganteus (Mg) and the wallaby M. rufogriseus (Mr) in the native grassland habitat, (c) macropods in the forest habitat and (d) macropods and introduced sheep (Shp) in the alienated grassland habitat. In this simulation, the thylacine population goes extinct in 1901. All population sizes are in thousands of individuals, except that sheep abundance has been reduced by a factor of ten to facilitate presentation. The vegetation component of the metamodel is not shown.

Metamodels that accounted for the effects of Europeans on the thylacine's prey base simulated extinction across the full range of starting population sizes (Fig. 3 ). From 10 000 metamodel simulations, one for each hypercube division, we obtained 3 295 simulations for which the thylacine population size averaged between 2 000 and 4000 individuals over the 500 year burn‐in period. After 1803, macropod abundance in the grassland habitat was reduced rapidly by competition from the growing sheep population (Fig. 3 b). This effect, in combination with concurrent habitat loss, caused the availability of the thylacine's prey to decline. By the first year of the bounty scheme in 1886, thylacine abundance had on average declined by 52% in the metamodels, compared with only 38% in the single‐species models, despite the former employing a more conservative scenario of habitat loss.

Boosted regression tree (BRT) summary of the sensitivity analysis for single‐species models. (a) Relative influence of model parameters on the probability of simulating the thylacine extinction, for the selected BRT model including five predictors only. (b) Plots of interactions between the most influential predictors of model success (thylacine r max , pre‐European population size and the harvest multiplier used to calculate the actual thylacine harvest implemented for each simulation). Contour lines illustrate a 50% probability of simulating the thylacine extinction successfully, for five different starting population sizes, for scenarios that exclude or include simulation of a virulent disease epidemic (40% mortality rate increase over 1906–1909).

Using the sensitivity analysis model output, a BRT model including all predictors and five splits (i.e. fourth‐order interactions) predicted the success or failure of simulations with an average misclassification rate of only 8·9% (as measured using 10‐fold cross‐validation). A simplified model including only the top five predictors and single splits (i.e. no interactions) performed similarly well (misclassification rate of 9·0%), so we used this model to summarise the sensitivity analysis results (Fig. 2 ). As indicated by the baseline and pessimistic scenarios tested (Fig. 1 ), single‐species PVA models were unable to recreate the thylacine's extinction unless a high human harvest, small starting population size or low maximum population growth rate was assumed, even when the effects of a virulent disease were simulated over 1906–1909 (Fig. 2 ).

‘Classical’, single‐species PVA results for thylacines (). Mean thylacine population size (calculated from 1 000 simulations) for five starting (pre‐European) population sizes, and cumulative probability of extinction curves, for the baseline scenario excluding (a, b) and including (c, d) a simulated disease (20% mortality rate increase over 1906–1906), and the pessimistic scenario excluding (e, f) and including (g, h) a virulent disease component (40% mortality rate increase). Vertical guidelines indicate (i) European settlement in 1803, (ii) initiation of the government bounty payment for thylacines in 1886 and (iii) termination of the bounty scheme in 1909. Cumulative probability of extinction curves show 95% confidence intervals at selected times for a starting thylacine population size of 2 000 individuals (see 2 for details). PVA models with different starting population sizes are indicated with the following line formats: 2 000 (grey, solid), 2 500 (black, dashed), 3 000 (black, solid), 3 500 (grey, dashed) and 4 000 (black, dotted).

For the baseline, single‐species scenario, extinctions occurred only for the lower starting population sizes of 2 000 and 2 500 individuals (26·6 and 0·4% of simulations went extinct, respectively) (Fig. 1 a, b). At the higher starting population sizes, populations declined during the bounty period but rebounded to a lower new equilibrium population size, due to progressive alienation of habitat (Fig. 1 a). Inclusion of moderate disease‐induced increases (20%) in mortality rates over the period 1906–1909 had little effect (Fig. 1 c, d). For the pessimistic scenario, which included an annual thylacine harvest equal to twice that reported by the bounty records, extinction probabilities were raised and ranged from 100 to 20·3% for starting population sizes of 2 000 and 4 000, respectively (Fig. 1 e, f). Simulation of a virulent disease epidemic (mortality rates increased by 40% over 1906–1909) raised extinction probabilities further for larger population sizes (e.g. 78·1% probability of extinction for 4 000 individuals initially) (Fig. 1 g, h).

Discussion

In the absence of hard evidence, two indirect arguments are commonly proposed in support of a disease explanation for the thylacine extinction: (i) known European impacts cannot account for the extinction, leading to the deduction that another contributing phenomenon must have been the ultimate cause (Bulte, Horan & Shogren 2003) and (ii) government bounty records indicate that the thylacine population was sustainably harvested until 1905, after which the rapid decline in claimed bounties is most likely explained by a disease epidemic (Guiler 1985). We conclude that both these arguments are unnecessary. The thylacine extinction (including a rapid population decline) can be explained by the thylacine harvest interacting with other European‐imposed stressors.

If we assume that the government bounty records report the true thylacine harvest, single‐species PVA models can only simulate the thylacine's disappearance for relatively small, pre‐European thylacine population sizes (Figs 1, 2). This contrasts with the results of Bulte, Horan & Shogren‘s (2003) model of thylacine‐bounty hunter dynamics, which predicted persistence. However, the Bulte model was deterministic, unstructured and neglected important components of the extinction process that we include (e.g. inbreeding effects, demographic and environmental stochasticity) (Gilpin & Soulé 1986). Inclusion of a simulated disease epidemic (through increased mortality over the period 1906–1909) naturally raised the probability of extinction, but not dramatically so (Figs 1, 2). Of course, if we had simulated a thylacine disease that inflicted close to 100% mortality like devil facial tumour disease (McCallum 2008), then more simulated extinctions would result. Importantly though, even when a disease was not included, single‐species models simulated a thylacine population crash in the early 1900s, from which some populations could not recover (Fig. 1). This decline can be explained by the peak in bounty payments recorded at the turn of the century. Uncertainties regarding the pre‐European thylacine abundance notwithstanding, thylacine harvest rate and r max were more important determinants of extinction probability than the disease parameters used (Fig. 2). Starting population sizes are always critical components of PVA models, and we assumed a commonly cited range of 2 000–4 000 thylacines initially (Guiler & Godard 1998). Owen (2003, p. 26) suggests a slightly higher estimate (5 000) but provides no supporting evidence for this figure. Nevertheless, if we assumed higher pre‐European abundances, the ability of these PVA models to recreate the thylacine's extinction would be reduced even further.

Incorporating the thylacine PVA model within a metamodel that linked thylacine abundance to prey availability improved our capacity to simulate the thylacine's demise. This different outcome is explained by the metamodel's ability to integrate the thylacine's numerical response with spatial variation in prey abundance and the pattern of land alienation by Europeans. Even when we assumed moderate competition between sheep and kangaroos, the metamodels simulated the thylacine extinction more readily than single‐species models including a disease epidemic. Of course, we assumed that thylacines were not major predators of sheep, and it is abundantly clear that the portrayal of the animal as a ‘sheep‐killer’ by a small fraction of graziers was a gross exaggeration. Paddle (2000, p. 137) could only find six published accounts that actually detail thylacines preying on sheep, the last in 1871 (15 years prior to the government bounty scheme). Further, none of the official reports of the Tasmanian sheep industry list the thylacine as even a minor threat (Paddle 2000, p. 143). While we cannot discount the possibility that some thylacines near grazing lands might have preyed on sheep, we consider it unlikely that sheep sustained the thylacine population as a whole in any meaningful way. The government bounty records indicate that, by the late 1800s at least, thylacines were mostly distributed in areas removed from the prime grazing lands (Davies 1965; Guiler 1985).

The influence of introduced herbivores on native Australian herbivores, including their potential contribution to species extinctions, has been hotly debated (e.g. Edwards, Croft & Dawson 1996; Fisher, Blomberg & Owens 2003). Most recent estimates suggest that sheep consume pasture at approximately 2·5 times the rate of kangaroos (Munn et al. 2008; Munn, Dawson & McLeod 2010). However, the negative impacts of sheep have been hard to demonstrate experimentally, partially due to logistical problems with conducting these experiments (Edwards, Croft & Dawson 1996). We have shown that competition from introduced sheep provides one threatening process that could have contributed to the thylacine's extinction. However, our metamodels probably underestimate the true impact of Europeans on the thylacine population. Europeans also introduced other herbivores (particularly cattle) and transformed large tracts of open habitat for agricultural purposes (Davies 1965). Further, dogs (both domestic and feral) rapidly spread throughout Tasmania after European settlement and might well have adversely affected thylacines via direct predation and/or interference competition (Boyce 2006). One interesting feature of our metamodels was that macropod harvesting by humans, treated as a proportional harvest that increased over 1803–1935, had little effect on simulation output. Instead, the pre‐European thylacine abundance and sheep competition coefficient were more influential parameters because they strongly affected thylacine abundance prior to the bounty scheme in 1886. There is ample evidence that thylacines were in decline before this; by the late 1870s, for example, Tasmanian naturalists accepted that the species was extremely rare (Paddle 2000, p. 136).

Given that the thylacine's demise can be explained by known European impacts, the case for a disease epidemic remains weak. Thylacines and dasyurid mammals reportedly disappeared from one Tasmanian area in around 1910, possibly due to a ‘distemper‐like’ disease (Guiler 1985). To our knowledge, however, there are no unambiguous cases of diseased thylacines from the wild. Melbourne Zoo lost sixteen thylacines to an unknown disease over the period 1901–1903 (Paddle 2000), but it is impossible to extrapolate these events to the island of Tasmania, particularly when comparable losses were not experienced in any other zoo. The disease explanation appears to be an ad hoc hypothesis, born of bewilderment at witnessing the rapid decline of a formerly ubiquitous species, underestimation of negative European impacts and poor comprehension of the power of synergies among anthropogenic stressors to elevate extinction risk (Brook, Sodhi & Bradshaw 2008).

It is likely that many extinctions result from disrupted interactions within complex biological systems (Lennartsson 2002). Single‐species modelling approaches can be unsuitable for studying such processes for a number of reasons. First, these models conflate environmental variation with predator–prey cycles, leading to overly pessimistic predictions of extinction risk (Brook 2000; Sabo 2008). Second, single‐species approaches typically necessitate carrying capacity adjustments for a focal species in response to external stressors. For example, Bulte, Horan & Shogren (2003) assumed that habitat loss and prey destruction by Europeans could be represented by a linear reduction in thylacine carrying capacity of 3% per year. This is an arbitrary figure and difficult to justify. Finally, without simultaneous modelling of the relevant system components, proper consideration cannot be given to the causes of past extinctions or the relative importance of threatening processes.

vortex software provides an individual‐based model of the extinction process and has been used extensively over the last two decades for conservation assessment and planning (Lacy 2000). We have extended the PVA concept to incorporate dynamic species interactions by integrating vortex population models within a single metamodel using the command‐centre package metamodel manager (Raboy 2011; Bradshaw et al. 2012). This allows complex systems models to be flexibly constructed from linked component models, while affording the advantages of stochastic PVA‐based approaches for the evaluation of both historical extinctions and contemporary threatened species. The thylacine extinction provides an ideal case study because, while difficult to explain by single factors, we know that the species was confronted with multiple extinction drivers. We find that the thylacine's demise can be explained by known negative European impacts and suggest that their interacting effects were powerful enough that, even without a disease epidemic, the species was committed to extinction.