Dynamics and predicted distribution of an irrupting ‘sleeper’ population: fallow deer in Tasmania

Sleeper populations of non-native species can remain at low abundance for decades before irrupting. For over a century, fallow deer (Dama dama) in the island state of Tasmania, Australia, remained at low abundance and close to the region in which they were released. Recently, there are indications the population has increased in abundance and distribution. Here, we spatially quantify the population change using a time series of annual spotlight counts from 1985 to 2019 (up to 172 transects annually, totalling of 5756 transect counts). Next, we predict the potential for further range expansion, using global occurrences to characterise the species’ climatic niche, and remote-camera surveys (3225 camera sites) to model fine-grained habitat suitability. Spotlight counts of fallow deer increased by 11.5% annually, resulting in a 40-fold increase from 1985 to 2019. The core distribution increased 2.9-fold during this 35-year period, and now spans c. 27% of Tasmania’s land area. Satellite populations have established in locations where farmed deer have escaped or been released, suggesting that humans have facilitated range expansion via new introduction events. Based on climate and habitat suitability, our models predict that 56% of Tasmania is suitable under the current climate. This suggests range expansion is likely to continue unless the population is actively managed, which could include the eradication of satellite populations and containment of core populations. This case study cautions that despite over a century of slow population growth, sleeper populations of non-native species can abruptly increase.

Abstract Sleeper populations of non-native species can remain at low abundance for decades before irrupting. For over a century, fallow deer (Dama dama) in the island state of Tasmania, Australia, remained at low abundance and close to the region in which they were released. Recently, there are indications the population has increased in abundance and distribution. Here, we spatially quantify the population change using a time series of annual spotlight counts from 1985 to 2019 (up to 172 transects annually, totalling of 5756 transect counts). Next, we predict the potential for further range expansion, using global occurrences to characterise the species' climatic niche, and remote-camera surveys (3225 camera sites) to model fine-grained habitat suitability. Spotlight counts of fallow deer increased by 11.5% annually, resulting in a 40-fold increase from 1985 to 2019. The core distribution increased 2.9-fold during this 35-year period, and now spans c. 27% of Tasmania's land area. Satellite populations have established in locations where farmed deer have escaped or been released, suggesting that humans have facilitated range expansion via new introduction events. Based on climate and habitat suitability, our models predict that 56% of Tasmania is suitable under the current climate. This suggests range expansion is likely to continue unless the population is actively managed, which could include the eradication of satellite populations and containment of core populations. This case study cautions that despite over a century of slow population growth, sleeper populations of non-native species can abruptly increase.

Introduction
Sleeper populations of introduced species can persist at low abundance for decades before being triggered by an environmental factor to become abundant and problematic (Spear et al. 2021). Deer are among the most successful biological invaders globally (Clout, Russell 2008;Davis et al. 2016), and six species established wild populations in Australia in the mid-1800s (Bentley 1998;Davis et al. 2016). As is common among biotic invasions (Crooks 2005), these deer populations initially remained small and localised in ranges reflecting the idiosyncratic pattern of release rather than optimal habitat (Caley et al. 2011). During the past four decades, however, deer populations across Australia have increased substantially and now pose threats to conservation and agriculture (Davis et al. 2016). These recent increases (e.g., Forsyth et al. 2018) are characteristic of the irruptive dynamics that can occur when non-native species (Spear et al. 2021), and ungulates specifically (Caughley 1970;Forsyth, Caley 2006), are introduced to new areas, released from predation, or favoured by environmental change.
Fallow deer (Dama dama) were introduced to the large island state of Tasmania, Australia (c 68,400 km 2 ) in 1829 (Chapman, Chapman 1980). During the next c. 150 years, the Tasmanian population of fallow deer remained at low abundance within c. 60 km of the initial release regions in the Midlands and eastern regions of Tasmania. Based on field observations but not systematic surveys, Wapstra (1973) conservatively estimated the population to be 7-8000 in the 1970s. Since then abundance has evidently increased, with an aerial survey of eastern Tasmania in 2019 estimating a population size of 53,000 deer (Lethbridge et al. 2020a). This survey was conducted late in the 2019 hunting season, during which c. 30,000 deer were killed (Lethbridge et al. 2020a), suggesting the population numbered at least c. 83,000 at its 2019 peak. Based on estimated climate suitability, linked with a conservative estimate of carrying capacity, Potts et al. (2015) suggest that Tasmania could support over a million fallow deer. There is evidence of range expansion (Commonwealth of Australia 2014), with deer reported in the eastern Central Plateau region of the Tasmanian Wilderness World Heritage Area (TWWHA; Driessen et al. 2020;Locke 2007). Due to the high natural and cultural value of the TWWHA and concerns about the damaging impact of deer, predicting their expansion into this area, and other currently unoccupied areas of Tasmania, is a management priority (Driessen et al. 2020).
Climate suitability has been a fundamental determinant of the spread of deer populations throughout Australia (Forsyth et al. 2004). Predicting habitat suitability for invading species is challenging because they violate the implicit assumption of correlative species distribution models (SDMs) that the species' range is at equilibrium with its environment (Gallien et al. 2012;Guisan, Thuiller 2005). Ignoring the nonequilibrium of invading species would usually result in under-predictions of habitat suitability (Guisan, Thuiller 2005). One way to circumvent this problem is to model suitability in areas where a species is at or close to equilibrium, and then project modelled suitability to areas of potential expansion (Guisan, Thuiller 2005).
Here, we seek to quantify changes in the abundance and distribution of fallow deer in Tasmania and predict the potential for continued population and range expansion. We begin by quantifying population trends and changes in the distribution of fallow deer in Tasmania over the last 35 years. To do this, we use a time series of annual spotlighting surveys of up to 172 transects across Tasmania from 1985 to 2019, totalling 5756 transect surveys. Next, we evaluate whether there is habitat in Tasmania that is suitable for fallow deer but not yet occupied. To do this, we characterise the climatic niche of fallow deer using c. 69,000 global occurrences from regions where populations are closer to equilibrium (including the species' native range in Europe). Finally, to characterise fine-scale habitat suitability, we use a dataset from 3225 camera traps across Tasmania. The habitat suitability modelling improves upon that of Potts et al. (2015) by (i) incorporating many more occurrences, (ii) developing models of suitability for populations at equilibrium and then projecting those models to areas of potential expansion, and (iii) incorporating recent camera surveys to delineate the current distribution and to model fine-scale habitat suitability.

Study area
Tasmania is Australia's largest island (68,400 km 2 ). Fifty one percent of land is reserved under some form of conservation classification (DPIPWE 2020), mainly in the western half of the state. The west of the island has rugged and mountainous topography (Fig. 1D), high annual rainfall (Fig. 1C), and vegetation dominated by dense wet sclerophyll forests, rainforests, alpine treeless vegetation and low-productivity buttongrass moorlands (Fig. 1E). Large parts of western and central Tasmania are considered wilderness, with very few people, roads and limited commercial development. The eastern half of the island and the northern coastal areas have gentler topography (Fig. 1D), moderate rainfall ( Fig. 1C) and drier forest types, making these areas more suitable for farming (Figs. 1E) and producing habitats that are more typical of the mixed forest-grassland preferences of fallow deer (Apollonio et al. 1998;Chapman, Chapman 1980).

Data sources (i) Spotlighting counts
The Tasmanian State Government has conducted annual spotlight surveys along 172 transects across northern, eastern, and central Tasmania from 1985 to 2019. The surveys were established to monitor harvested herbivore populations, while also recording all sightings of free-ranging mammal species, including fallow deer . The same transects are surveyed each year; however, the number of transects increased from 128 in 1985 to 172 in 1998 (Table S1). The transects do not sample large areas of western Tasmania, largely due to inaccessibility, and also because wildlife are not harvested in those areas Driessen, Mallick 2003). Each transect comprises a 10-km section of road and adjacent land. Transects are surveyed at night with the aid of a spotlight from a vehicle driven at 20 km/hr (for details, see . To ensure detection probability did not vary through time because of differences in the spotlight device, all surveys used a 100-W halogen spotlight. Transects are surveyed once each year during the summer months. This ensures comparability between years but prevents the use of statistical techniques that separate the detection and state processes using repeat surveys within a year. We consider the count of fallow deer per transect as an index of abundance, but we acknowledge that the spotlight counts reflect the dual processes of detection and abundance. Because we are unable to explicitly model both the detection and abundance processes with these data, we must keep in mind that any changes in the abundance index cannot solely be assigned to changes in density. Nevertheless, these spotlight surveys have shown population declines of a similar magnitude to that estimated by concurrent mark-recapture surveys for Tasmanian devils (Sarcophilus harrisii) (Cunningham et al. 2021) and linetransect methods for common wombats (Vombatus ursinus) (Carver et al. 2021), so we expect the spotlight surveys are a reasonable measure of population change.
We used camera traps to model the current distribution and habitat preferences of fallow deer in Tasmania. We compiled a database of our own camera deployments between 2014 and 2020, totalling 3225 camera sites, and 315,120 camera-nights (Table S2). Cameras were deployed by different research teams using a variety of camera types (Table S2). Importantly, there were no geographical or temporal biases to the use of particular camera types, with all being deployed Tasmania-wide. For each camera, we recorded the number of unique detection events of fallow deer, the duration of camera deployment, and environmental covariates (described below). We defined a detection event as unique if at least 30 min separated consecutive occurrences of fallow deer on the same camera.
(iii) Global occurrences of fallow deer To characterise climatic conditions suitable for fallow deer, we collated occurrence records from regions where fallow deer are likely to be near to equilibrium. We used citizen-science data from Europe, where fallow deer occur as both a native and long-established introduced species, and the Australian mainland, where fallow deer have been widely introduced across a broader range of environmental conditions than in Tasmania. We compiled occurrences in Europe using the Global Biodiversity Information Facility (GBIF) database, and in Australia using the Atlas of Living Australia. Both databases are comprised of species occurrences (presence-only) and derived from a similar range of sources, including scientific studies, museum collections and citizen-science observations. Broadly speaking, there are analogous climates between regions of south-eastern Australia and Europe (Peel et al. 2007). Tasmania's climate is analogous to the United Kingdom, Ireland and parts of France, whereas parts of pseudo-absences (orange) used to create the model of climate suitability. In Europe (A), pseudo-absences were distributed according to a proxy of survey bias, which is why eastern Europe has few pseudo-absences. Maps of some relevant bioclimatic variables across Tasmania, including elevation (C), which ranges from 0 to 1617 m, and mean annual rainfall (D), which ranges from 490 to 3100 mm, with blue representing lower values for both (Australian Bureau of Meteorology). (E) The four main vegetation types in Tasmania, with buttongrass moorlands shown in brown, wet eucalypt forest and rainforest in blue, agricultural or grassland in tan, and dry eucalypt forest in green (Department of Primary Industries Parks Water and Environment 2020) Victoria and New South Wales have climates similar to Italy, France and Spain (Peel et al. 2007).
For Europe, we downloaded 58,381 presence-only occurrences of fallow deer from GBIF (2020). We selected only human-observed records (i.e., excluding fossil records). To reduce the influence of intensively sampled areas, we randomly thinned records so that at least 5 km separated occurrences, leaving 4726 occurrences (Fig. 1A). Substantial regional biases in survey effort were apparent across Europe, with many fewer mammals reported in eastern than in western Europe (Fig. 1A). Ignoring this survey bias would increase the likelihood of misclassifying areas with low survey effort as unsuitable for fallow deer (Phillips et al. 2009). Thus, as recommended (Phillips et al. 2009), we selected pseudo-absences according to the same survey bias that generated the occurrences. We did this selection using reported occurrences (excluding fossils) of any mammal (excluding fallow deer) in the GBIF database as a proxy of relative survey effort across Europe. From these mammal occurrences, we randomly selected 4726 records as pseudo-absences, ensuring they were not in the immediate neighbourhood of known presences (i.e., within 10 km). Distributing pseudo-absences in this way helped to distinguish heavily sampled areas where fallow deer are probably absent (e.g., northern Norway; Fig. 1A) from areas where there is low survey effort and that could contain suitable habitat (e.g., eastern Europe; Fig. 1A).
For Australian occurrences of fallow deer, we downloaded 10,925 presence-only occurrences from the Atlas of Living Australia (Atlas of Living Australia 26 August 2020). We removed three unconfirmed and probably erroneous records from the South Australian desert and randomly thinned records so that occurrences were separated by at least 5 km, leaving 490 occurrences (Fig. 1B). Unlike Europe, fallow deer have not had the opportunity to disperse to all suitable habitat in Australia. Thus, we used a different approach to distribute pseudo-absences, aiming to distribute them across the area that we expect fallow deer have had the opportunity to occupy. We also included areas around failed introduction attempts that might indicate unsuitable habitat; most notably the Cobourg Peninsula, Northern Territory, where 52 fallow deer were introduced in 1912, but did not establish (Chapman, Chapman 1980). In total, we randomly distributed 490 pseudo-absences at least 5 km apart within a 200 km buffer around occurrences and failed introductions on the Australian mainland, and within a 90% minimum convex polygon in Tasmania (where deer populations are clearly expanding, and therefore not at equilibrium at the current range boundary). As for the European data, we ensured that pseudo-absence locations were not within 10 km of a known occurrence.

Statistical analysis (i) Population trends from 1985 to 2019 and current distribution
First, we fitted three plausible models to describe the Tasmania-wide trend in the mean number of deer detected each year from 1985 to 2019. We modelled mean deer detections in response to (i) a linear effect of survey year, (ii) a quadratic effect of survey year, and (iii) an exponential effect of survey year (i.e., using the log-transformed mean deer detections as the response variable). The exponential model follows the hypothesis that the apparent increase in the population might be the result of exponential population growth (i.e., a constant growth rate). We fitted these models using the 1m function in R.
Next, we modelled spatial changes in the spotlight counts using a spatiotemporally autoregressive model (Gómez-Rubio 2020). This approach correlates deer detections in space and time, with the goal of characterising changes in the actual distribution of fallow deer. We fitted the autoregressive model using integrated nested Laplace approximation (INLA) (Illian et al. 2013), which is a computationally efficient option for Bayesian inference from spatial data. In INLA, spatial dependence between observations is modelled using a spatially continuous correlation process known as a Gaussian random field (Bachl et al. 2019). Continuous space is approximated by the solution to a stochastic partial differential equation (SPDE), which is based on a discretised 'mesh' of abutting triangles (for details, see Lindgren et al. 2011). To avoid problematic boundary effects (Bakka et al. 2018), we extended the mesh beyond the boundaries of Tasmania, with maximum interior edge lengths of 10 km and maximum exterior edge lengths of 30 km. We first attempted to model a spatiotemporal random field using first-and second-order autoregressive models. These models had convergence problems, likely because of substantial year-to-year variation in the numbers of deer counted on some spotlight transects. Thus, to smooth some of the year-to-year variability, we grouped surveys into five-year periods, and then modelled the spatiotemporal random field using a first-order autoregressive model (Gómez-Rubio 2020), which correlated observations in consecutive five-year periods. See Gómez-Rubio (2020) for more details of spatiotemporal random fields. To account for overdispersion, we modelled the count of deer per transect using a negative binomial distribution (Lindén, Mäntyniemi 2011;Ver Hoef, Boveng 2007). We visualised the results by producing maps of predicted deer spotlight counts in each five-year period.
Finally, we estimated the recent distribution of fallow deer by simultaneously incorporating the camera data and recent spotlight surveys (2015-2019) into a joint-likelihood model. Joint-likelihood models combine multiple sources of data into single integrated models that estimate one or more shared parameters (Isaac et al. 2020;Miller et al. 2019). In this case, we used the joint model to estimate a shared spatial random field (explained in previous paragraph) that describes the distribution of fallow deer. The joint model was comprised of two sub-models, one for each data source, and took the form of where b 1 and b 2 denote intercepts, and random field denotes a shared spatial random field. We modelled both data sources using the negative binomial distribution to account for overdispersion in the count data (Lindén, Mäntyniemi 2011;Ver Hoef, Boveng 2007). From this model, we predicted a map of fallow deer distribution. We note that due to imperfect detection, this map is inherently a conservative estimate of the current distribution. We fitted the spatiotemporal and joint-likelihood models using v2.1.13 of the inlabru R package (Bachl et al. 2019; R Core Team 2019), an interface for the R-INLA package (Bakka et al. 2018;Rue et al. 2009).
(ii) Climate suitability We constructed a species distribution model using the combined Australian and European occurrence data. We modelled deer occurrences in response to explanatory variables using the WorldClim global climate dataset v2.1 (Fick, Hijmans 2017), which provides global rasters of 19 climatic and environmental variables. Using rasters with a cell size of 2.5 min, we selected five uncorrelated (Pearson's r \|0.7|; Dormann et al. 2013) predictor variables that we expected would influence the distribution of fallow deer: (i) maximum temperature of the warmest month ('maxTemp'); (ii) minimum temperature of the coolest month ('minTemp'); (iii) average annual precipitation ('precip'); (iv) the seasonality of precipitation, measured by the coefficient of variation ('precipCV'), and (v) terrain ruggedness index ('ruggedness') calculated from elevation (tri function of the spa-tialEco package; Riley 1999);.
We constructed the species distribution model as an ensemble of generalised additive, random forest and boosted regression tree models. We fitted the model using five-fold cross validation with 10 replicates for each of the three model types, totalling 150 model runs (defaults for all other settings). We fitted the SDM using v1.0-89 of the sdm package (Naimi, Araújo 2016) in R v4.0.4 (R Core Team 2020). We projected an ensemble map of climate suitability, with projections from individual models weighted according to out-of-sample area under the receiver operating-characteristic curve (AUC) (Marmion et al. 2009).
(iii) Fine-scale habitat suitability Because fallow deer have probably not yet dispersed to all suitable habitat, we modelled fine-grained habitat suitability using camera data from the core distributional range in Tasmania where they are likely to be closer to equilibrium. We defined the core range using a 90% minimum convex polygon around the locations where fallow deer have been observed on camera. There were 674 camera sites within the corerange polygon, of which we withheld 20% as a test dataset for model validation. We projected the model of habitat suitability in the core range across all of Tasmania.
We compiled five explanatory variables that we expect would influence fallow deer abundance and distribution: (i) the percentage of grassy vegetation types (%grassyVeg); (ii) the percentage of forest cover (%forest); (iii) the percentage of alpine treeless vegetation (excluding alpine grasslands); (iv) terrain ruggedness index and (v) year of survey. For covariate details, see Table S3. We expected that grassy vegetation types would positively influence fallow deer abundance (King, Forsyth 2020;Nugent 1990;Putman et al. 1993;Spitzer et al. 2020). We expected that highly rugged terrain would negatively affect deer abundance, and that alpine treeless vegetation-a lowproductivity vegetation group occurring in harsh climates-might negatively affect deer abundance. We excluded %forest from model fitting because it was negatively correlated with %grassyVeg (Pearson's r = 0.76). Because these variables could have effects at different spatial scales (Connor et al. 2018), we created them using a buffer of three different sizes around the camera locations: 250 m, 500 m and 1000 m. Because cameras were deployed for different time periods, we used the natural logarithm of survey duration as an offset to scale deer detections according to the length of survey.
We used generalised additive mixed models (GAMM) to investigate the effect of habitat on the daily count of unique camera detections of deer. To account for overdispersion, we fitted this model using the negative binomial distribution (Lindén, Mäntyniemi 2011; Ver Hoef, Boveng 2007). We fitted all combinations of smooth effects of the explanatory variables at different spatial scales, excluding variables that had correlations of r [|0.7| from the same model. Thus, different spatial scales of the same vegetation type were not included in the same model. To prevent overfitting, we a priori restricted smoothing parameters to a maximum of five knots. Because each camera had a single count of the number of deer detections over the survey period, we did not include a camera-level random effect (i.e., observation-level random effect) because that would not impose any relationship between nearby cameras. However, to account for the non-independence of nearby cameras, we grouped cameras into clusters (Fig S1) using hierarchical clustering (hclust function in R) with a cut distance of 10 km (cutree function). We included cluster ID as a random effect in all models. We fitted the GAMM using the mgcv package (Wood 2017) in R (R Core Team 2020).
We ranked models using small-sample corrected Akaike information criterion (AIC c ) (Burnham, Anderson 2002). Because several models had similarly strong support, we produced a model-averaged map of habitat suitability across Tasmania by weighting all models within five DAIC c according to their model weights. To produce the habitat suitability map, we used explanatory vegetation variables at the spatial scale indicated by each model. For instance, if %grassyVeg at 1000 m was in the model, the raster layer for projection was the percentage cover of grassy vegetation types within 2000 9 2000 m pixels (i.e., twice the buffer width because the buffer goes in all directions).
(iv) Combining climate and habitat suitability We expected that climate would influence the broad-scale distribution of fallow deer, and that deer would discriminate suitable habitat mostly from within climatically suitable areas. Nevertheless, it is possible that deer could use habitat that our models suggested was in climatically unsuitable areas. We therefore used a weighted product model (Tofallis 2014) to combine the maps of climate and habitat suitability. We first scaled the maps of climate and habitat suitability by dividing pixel values by the maximum value of each map. This constrained suitability to the range 0-1 but did not force any values to zero, which is advantageous because it preserves the proportionality of the data (Tofallis 2014). Next, we combined these scaled maps of climate (C) and habitat suitability (H) using a weighted product model where each map is raised to the power of its weight and then the maps are multiplied together. We determined the weights for each map by selecting the combination of weights, constrained to sum to one, that maximised the AUC statistic for the 20% of camera data that was withheld from model fitting. This resulted in a combined map that was calculated as C 0.78 9 H 0.22 . Finally, we divided the resulting map by the maximum value, constraining relative suitability values to a maximum of one.

Population trends
The annual spotlight surveys revealed a steady and large increase in the spotlight counts of fallow deer from 1985 to 2019 ( Fig. 2A). The Tasmania-wide trend in mean deer detections was best fit by an exponential effect of year (adjusted R 2 = 0.71; coefficients ± SE: Intercept = 218.2 ± 23.63, byear = 0.109 ± 0.012), relative to quadratic (adj. R 2 = 0.61) and linear (adj. R 2 = 0.55) effects of year. The exponential model suggests the spotlight counts have increased at an approximately constant rate of 11.5% (SE: 1.2%) per year over the last 35 years, representing a 40-fold increase from 1985 to 2019 ( Fig. 2A). The exponential model, however, did not explain a large increase followed by a decline in deer detections in the mid-2000s ( Fig. 2A). This pattern coincided with a period of above-average rainfall followed by a period of well-below-average rainfall across much of the species' Tasmanian range (Fig S2). The number of deer reported as being killed under crop protection and hunting permits has increased substantially in the last five years (Fig. 2C), but this offtake has evidently not prevented population growth.
The spatiotemporal model of the spotlight surveys shows that fallow deer distribution has steadily expanded (Fig. 2D), with a 2.9-fold increase in the model-estimated area of 'core' distribution since [1985][1986][1987][1988][1989] (Fig. 2B). The joint likelihood model, which simultaneously modelled the spotlight and camera datasets, conservatively estimates that fallow deer now occupy c. 27% of Tasmania's land area (Fig. 3). Much of the expansion in deer distribution has occurred in areas where there are anecdotal reports of farm escapes or intentional releases (Fig. 2D).

Climate suitability
The climate suitability models revealed that large areas of Tasmania that are currently unoccupied by fallow deer have climates that could support the species (Fig, 4E). Maximum temperature of the warmest month and minimum temperature of the coldest month had the highest relative importance (Fig. 4C). The concave response curves for the temperature variables and average annual precipitation suggest upper and lower constraints of climate suitability. For instance, maximum temperatures of less than * 18°C in the warmest month, and  Lethbridge et al. 2020b). The black polygon shows the boundary of the Tasmanian Wilderness World Heritage Area minimum temperatures less than * 7°C in the coolest month, resulted in an abrupt reduction in climate suitability (Fig. 4D). Highly rugged terrain had a negative effect on fallow deer occurrences (Fig. 4D). All models performed better than random, with out-of-sample AUC values of 0.89, 0.75 and 0.76 for the random forest, GAM and boosted regression tree, respectively. Based on an ensemble of these models, weighted by out-of-sample AUC (Marmion et al. 2009), much of Tasmania has a suitable climate for fallow deer, with highest suitability in the east and the north and the two largest off-shore islands (Fig. 4E). Due to its high annual rainfall the southwest of Tasmania had relatively low climate suitability for fallow deer, and the Central Plateau and highlands in the northeast had low climate suitability due to low temperatures.

Fine-grained habitat suitability
The fine-grained habitat suitability model based on camera trap detections of fallow deer revealed large areas of potentially suitable but as-yet unoccupied habitat in northern, eastern and southern Tasmania (Fig. 5B). All top-performing models included a positive effect (Fig. 5C) of %grassyVeg within 1000 m of a camera (variable importance = 1; Table S4). There was strong support for topographic ruggedness negatively influencing (Fig. 5C) deer detections, with similar support for this variable at spatial scales of 500 m and 1000 m (cumulative importance of the different spatial scales of ruggedness = 0.93; Table S4). Survey year had a positive effect on deer detections (Fig. 5C). Given the strong positive effect of %grassyVeg, most predicted suitable habitat (Fig. 5B) is in the agricultural regions of Tasmania or areas with native grasslands and eucalypt woodlands with grassy understories. Importantly, though, habitats without a major grass component still seem able to support low deer densities (i.e., positive y-intercept; Fig. 5), as do rugged areas. This suggests that although much of western Tasmania has low habitat suitability, it may be capable of supporting low densities of deer.

Combined climate-habitat suitability
We combined the maps of climate and habitat suitability using weights that maximised the AUC of withheld test data, indicating a weight of 0.78 for habitat suitability and 0.22 for climate suitability (Fig. 6C). The combined map of climate-habitat suitability (Fig. 6D) suggests that 56% of Tasmania is suitable habitat (Fig. 6E, categorized based on a threshold of 0.083, as determined by maximising the true skill statistic; Allouche et al. 2006). If fully occupied (Fig. 6E), this would represent a doubling of the current range extent (Fig. 3). Fig. 3 The estimated distribution of fallow deer. We estimated the distribution using a joint-likelihood model that simultaneously modelled the spotlight (2015-2019) and camera datasets. Dots show cameras and crosses show spotlight transects where deer have been observed. This figure shows the jointly estimated spatial random field, with the contour denoting the area occupied by deer, as determined using a 0.025 omission threshold. This map is inherently a conservative estimate of distribution due to imperfect detection

Discussion
Tasmania's population of fallow deer is increasing in abundance and distribution. The range increase has been concentrated in areas currently disjunct from the main distribution, in places where there have been known escapes from farms and intentional but unlawful releases of deer. Coupled with a localised distribution during the previous c. 150 years, this suggests that humans have assisted the recent colonisation of new areas. Our models conservatively estimate that deer now occupy 27% of Tasmania. Based on global occurrences of fallow deer, much of Tasmania is climatically suitable under the current climate, and based on our extensive camera dataset, there also appears to be a large amount of unoccupied suitable habitat. Collectively, the combined climatehabitat suitability suggests 56% of Tasmania, primarily in the east and north of the state, is suitable for Fig. 4 Climate suitability of Europe, Australia and Tasmania for fallow deer. The climate suitability of fallow deer, fitted to occurrence data (black dots) from Europe (A), mainland Australia (B) and Tasmania (E). The relative variable importance (C) and response curves (D) calculated across all model runs and weighted according to out-of-sample AUC, with error bars denoting the 95% confidence interval across the different model runs fallow deer, which would represent more than a doubling of the current range area.

Population trends
Non-native sleeper populations can remain at low abundance for decades before irrupting (Crooks 2005;Spear et al. 2021). Until recent decades, fallow deer remained at relatively low abundance in Tasmania (\ 1% of estimated carrying capacity; Potts et al. 2015). However, over the last 35 years, spotlight counts of deer have increased at * 11.5% annually, suggesting the population is now in the increasing phase of the irruptive dynamics that can occur when ungulates are introduced to new environments lacking predators (Caughley 1970;Forsyth, Caley 2006;Spear et al. 2021). The 11.5% annual increase in spotlight counts cannot be assumed to be directly equivalent to growth rate because we were unable to explicitly separate abundance from the detection process (MacKenzie, Kendall 2002). However, we are not aware of any large changes through time that would increase detection probability, and given that each transect was repeatedly surveyed using a standardised protocol, we expect no systematic temporal bias in detectability. Further, this rate of increase matches the 11% annual growth of a Mediterranean population of fallow deer (Focardi et al. 1996). It remains unclear why the Tasmanian deer population was relatively stable for over a century preceding the current increasing phase, but this is common for exotic species, which often exhibit time lags and exist at low numbers for decades before increasing rapidly (Crooks 2005;Spear et al. 2021). Such irruptions need not occur as a result of a sudden change or perturbation; for instance, a gradual shift in a climatic driver can lead to abrupt population dynamics if the population is controlled by non-linear or threshold responses (Ratajczak et al. 2018;Spear et al. 2021). Wapstra (1973) suggested that poisoning of deer prior to the 1980s might have kept the population small. It is Although we did not survey King Island (the north-western island) with remote cameras, the island is known to have a wild deer population. (B) Habitat suitability map of fallow deer using data from the equilibrium range (black polygon), and then projected across all of Tasmania. (C) The effect of the variables in the best-performing GAM also possible that the change from relative population stability to irruption was triggered by a change in environmental conditions (e.g., rainfall; Fig S2) or changes to resource availability, such as increased irrigation (Romanin et al. 2015) and a decline in stocking numbers of wool sheep with which deer would compete for forage (Hassall & Associates Pty Ltd 2006;Kirkpatrick, Bridle 2007).
Until the mid-1970s-150 years after the initial introduction -Tasmania's fallow deer population had remained within c. 60 km of the region where deer were originally released (Wapstra 1973). This suggests range extension over the first 150 years occurred at an average of c. 0.4 km/year, which is consistent with slow dispersal of the species elsewhere in their invasive range. In New Zealand, distribution expanded at 0.8 km/year at Wanganui (Caughley 1963), 0.6 km/ year at South Kaipara (Lane 1982 cited by King & Forsyth 2020), and at 0.8 km/year at Paparoa (Clarke 1976 cited by King and Forsyth 2020). In Nebraska, fallow deer extended their range more quickly at 5.6 km/year (Packard 1955). During the last two Fig. 6 Combined climate and habitat suitability of Tasmania for fallow deer. We combined maps of climate (A) and habitat (B) suitability using a weighted product of the two maps. (C) We determined relative weights by selecting the combinations of weights that produced the highest AUC value for a portion of the remote-camera data that was withheld from model fitting, suggesting a weighting of 0.22 for climate and 0.78 for habitat. (D) The combined climate-suitability map, which when converted to a binary map E of habitat suitability (threshold determined by maximising the true skill score), suggests 56% of Tasmania is suitable for fallow deer decades, much of the range expansion in Tasmania has occurred in areas where deer have escaped from farms or were intentionally released (Fig. 1D), suggesting humans have played a key role in facilitating the recent range expansion.

Climate and habitat suitability
Our modelling showed that maximum and minimum temperatures had a substantial influence on climate suitability. Consequently, the cold, northern part of the Tasmanian Wilderness World Heritage Area (TWWHA), an area of high natural and cultural value, seems relatively unsuitable for fallow deer under the current climate. Nevertheless, climate change is expected to produce more favourable growing conditions in this area (Williamson et al. 2014). The response curves from our models (Fig. 4D) suggest that warming temperatures should make this region more climatically suitable for fallow deer. Climate change is also expected to increase the extent and severity of bushfires, which might create forests that are more open and herbaceous (Fairman et al. 2016) and therefore more suitable for fallow deer. Even lowdensity populations of deer in habitats like the central highlands could have significant impacts because the vegetation there is slow growing (Holz et al. 2015). Our models predict that low-density populations could occur over much of the northern TWWHA under current conditions. Recent large and high-severity bushfires in the central highlands (Jones 2019) should be exploited as natural experiments to reveal the effect of deer on ecosystem structure and function during periods of bushfire recovery.
Our model of habitat suitability suggested that fallow deer prefer flat areas with grassy vegetation types, probably because of a dietary preference for grasses (Elliott, Barrett 1985;Nugent 1990;Spitzer et al. 2020). Flatter, open areas are also favourable for gregarious behaviour and sight-based anti-predator responses of ungulates (Hopewell et al. 2005;Underwood 1982). The east-central part of the island, where deer first established and are currently most abundant (Fig. 3), is relatively flat and was originally dominated by native grasslands and open woodlands with grassy understories (Fensham 1989). This favourable vegetation type likely facilitated the establishment of the population. Much of the native grassy vegetation types in the midlands have been replaced by introduced grasses, which are possibly an even more palatable food source for fallow deer (Nugent 1990).
Some patches of high montane grasslands and shrublands in the northern and north-eastern TWWHA appear to be suitable habitat for fallow deer (Fig. 5B). These are the most likely locations for currently undetected populations and for the establishment of new populations, making them prime locations for systematic monitoring of incursions of deer into the TWWHA. For instance, deer have been detected in the last decade at Skullbone Plains on the Five Rivers Reserve (Tasmanian Land Conservancy 2019) and in the Walls of Jerusalem National Park (Fig. 6) (personal obs. Rob Buck c. 2010; personal obs. C Cunningham 2021), although not detected in the latter area during camera surveys in 2020 (M Driessen 2020 unpublished data). This area has large expanses of native grasslands and herbfields alongside relict stands of endemic conifers (Department of Primary Industries Parks Water and Environment 2020). The habitat suitability map indicates this area has well-connected suitable habitat, suggesting that the most likely route of dispersal into the northern TWWHA is from higher density populations to the east and south. The map also highlights the most likely dispersal corridors in other areas; for instance, the Fingal Valley (Fig. 6) is the most direct corridor to the east coast. We present rasters of the suitability and distribution maps in the Supplementary Material that may be useful in planning management to prevent expansion of populations into areas that are currently unoccupied but apparently suitable.
Although we characterised fine-scale habitat preferences in the high-density area of Tasmania, densities across much of this area are considerably lower (mean = 2.7/km 2 ; Lethbridge et al. 2020b) than those in some similar environments on the Australian mainland (e.g., 40.9/km 2 , 95% CI: 29.9-56.1; A.J. Bengsen & D.M. Forsyth, unpubl.). Consequently, deer might not yet occupy lower-quality habitat in Tasmania, which means the fine-scale map of habitat suitability likely underestimates the extent of suitable habitat. We combined the maps of climate and habitat suitability in a way that maximised out-ofsample AUC, which we expect will yield a more realistic potential distribution of fallow deer in Tasmania.
Given the deer population is likely to continue to expand into unoccupied, suitable habitat, the next important research step is to develop a spatially explicit demographic-dispersal model. This model could take the form of a grid-based spatial population viability analysis (e.g., Fordham et al. 2021;Visintin et al. 2020) that combines information about vital rates, dispersal capacity, initial population size and distribution (e.g., Fig. 2), and habitat suitability (e.g., Fig. 6). The model could be parametrised using a pattern-oriented modelling approach, which provides a systematic, data-oriented way of tuning complex simulation models by matching the simulation results with independent real-world data, known as 'targets' (Grimm, Railsback 2012;Grimm et al. 2005). In our case, Approximate Bayesian Computation (Csilléry et al. 2010) could be used to fine-tune the model parameters based only on the simulations that produce results matching the 'targets' for site occupancy, distribution and abundance (e.g., Fig. 2). Once parameterised, the model could forecast how deer distribution is likely to change under a suite of different management scenarios and compare counterfactual scenarios that investigate the long-term effects of deer escapes and unlawful releases. This approach could also test hypotheses relating to the long lag-times in range expansion, inflexion points and triggers of population growth.

Management suggestions
In light of our findings, we make several management suggestions to reduce the threats posed by expanding deer populations to environmental values and the agricultural sector, at a time when total eradication is not supported by Government. First, satellite populations could be assessed and prioritised for eradication. Obvious candidates for eradication are populations far from the traditional range, like Bruny Island, King Island, the Arthur Pieman Conservation Area (Fig. 6), and those in south-east Tasmania. Decades of invasion ecology show that early eradication is the most costeffective, long-term management strategy (Simberloff 2008;Taylor, Hastings 2004;Vander Zanden et al. 2010). The eradication of fallow deer on Australia's third largest island, Kangaroo Island (4405 km 2 ), demonstrates that isolated populations of fallow deer can be eradicated, but that community support and well-resourced management agencies are essential (Masters et al. 2018;Simberloff 2008). Eradication of a large population of red deer was also achieved on Secretary Island (81.4 km 2 ), New Zealand (Macdonald et al. 2019).
Next, we suggest that deer be prevented from expanding into protected areas, especially the Tasmanian Wilderness World Heritage Area. Habitats and pathways of likely spread could be identified and prioritised for monitoring and management. Deer subpopulations most likely to contribute to further spread could be identified and prioritised for management. We suggest that deer be contained to the 'traditional range' at a population size that allows landholders to achieve the outcomes desired from their land. Demographic models are needed to develop strategies to achieve pre-determined outcomes.
Once management strategies are implemented, the population should be monitored to determine the effectiveness of these strategies. Our habitat suitability maps could inform the establishment of a systematic surveillance network of camera traps in areas with and without deer. Coupled with the long-term spotlighting dataset analysed in this paper, a systematic camera network would be useful for monitoring the deer population, as well as detecting incursions into the TWWHA, escapes from farms, and unlawful releases.

Conclusion
After more than 100 years of relatively low and stable abundance and restricted distribution, the Tasmanian population of fallow deer is now increasing substantially in distribution and abundance. Our case study highlights that ungulate species, and invasive species more generally, can have 'sleeper populations' with long lag-times before the population irrupts (Caughley 1970;Crooks 2005;Forsyth, Caley 2006;Spear et al. 2021). The trajectory of the Tasmanian deer population, together with climate and habitat suitability modelling, indicate that the population is likely to expand further, including into areas of high conservation value. Demographic scenario modelling, which can be trained and validated using our population expansion and habitat suitability data, is needed to inform where expansion and population irruptions are most likely. Such modelling can help inform policy, frame management responses, and educate diverse stakeholders. Management interventions to prevent further spread of fallow deer in Tasmania could involve the eradication of small, isolated populations, containment of the core population, and prevention of further escapes and translocations.