Invasive species threaten biological diversity throughout the world. Understanding the dynamics of their spread is critical to mitigating this threat. In Australia, efforts are underway to control the invasive cane toad (Chaunus [Bufo] marinus). Range models based on their native bioclimatic envelope suggest that the cane toad is nearing the end of its invasion phase. However, such models assume a conserved niche between native and invaded regions and the absence of evolution to novel habitats. Here, we develop a dynamically updated statistical model to predict the growing extent of cane toad range based on their current distribution in Australia. Results demonstrate that Australian cane toads may already have the ability to spread across an area that almost doubles their current range and that triples projections based on their native distribution. Most of the expansion in suitable habitat area has occurred in the last decade and in regions characterized by high temperatures. Increasing use of extreme habitats may indicate that novel ecological conditions have facilitated a broader realized niche or that toad populations at the invasion front have evolved greater tolerance to extreme abiotic conditions. Rapid evolution to novel habitats combined with ecological release from native enemies may explain why some species become highly successful global invaders. Predicting species ranges following invasion or climate change may often require dynamically updated range models that incorporate a broader realization of niches in the absence of natural enemies and evolution in response to novel habitats.
Invasive species pose a major threat to global biodiversity (Elton 1958; Sala et al. 2000). Efforts to control invasive species can be enhanced by developing an accurate model of an invader's suitable habitat in the introduced region (Welk 2004; Butin et al. 2005; Mau-Crimmins et al. 2006). Such models provide critical information about where to allocate limited surveillance, control and conservation efforts (Wiens & Graham 2005). Usually invasive range models are based on the conditions found in a species' native range; therefore, they assume a conserved niche between native and introduced habitats (Sutherst et al. 1995; Wiens & Graham 2005). However, unique conditions, adaptation, phenotypic plasticity or random sampling of native genetic variation may reduce the accuracy of such models (Thomas et al. 2001; Lee 2002; Kolbe et al. 2004; Yeh & Price 2004; Butin et al. 2005; Facon et al. 2006; Mau-Crimmins et al. 2006). Despite the potential for violations to the assumption of a conserved niche in exotic habitats, few models predict the potential future distribution of an invasive species based on their actual distribution in invaded regions (Drake & Bossenbroek 2004; Welk 2004; Mau-Crimmins et al. 2006).
A massive invasion of Australia is currently underway by the cane toad, a global invader of incomparable success (IUCN 2001). Following its introduction in 1935 to control sugar cane pests, the cane toad has expanded its range to approximately 1.2 million km2 (figure 1a). Along the way, the toxic and often competitively superior cane toad has initiated declines in dozens of native species (Phillips & Shine 2004; Murray & Hose 2005). Continued expansion is likely to elicit further declines in the endemic fauna of Australia (Lever 2001; Murray & Hose 2005). In an earlier (1995) modelling effort, information on rainfall, humidity and temperature in its native range were used to construct a bioclimatic probabilistic index of toad population persistence that circumscribed their ancestral range limits in Central and South America (Sutherst et al. 1995). The probability of toad colonization was then calculated for point locations throughout Australia on the basis of these factors. Given a 50% probability of establishment, this prior model predicted the colonization of 0.7 million km2 of northern and eastern Australia (see electronic supplementary material for further details of range calculations). This prediction has already been essentially realized, and in many cases exceeded. Yet, expansion rates continue to increase rather than to decrease (Phillips et al. 2006). Evidence for post-introduction evolution (Phillips et al. 2006) and the recent colonization of habitat previously thought unsuitable (Sutherst et al. 1995) indicates an untapped potential for future range expansion and associated harm to endemic fauna.
Here, we forecast the cane toad's potential invasion range in Australia based on a similar set of bioclimatic variables as those used in the earlier model (Sutherst et al. 1995). However, instead of the native range, we base our calculations upon the cane toad's current (2005) distribution in invaded regions of Australia as the basis for a bioclimatic envelope. To this dataset, we also add anthropogenic factors thought to be aiding its spread (Brown et al. 2006). This effort represents the first statistically derived model of Australian cane toad's potential range based on its invasive range. In a departure from the previous methods, we apply a dynamically updated modelling approach that makes use of the data available over the course of the toad's invasion to inform range predictions. This approach illustrates that predictions of invasive distributions change through time and highlights a recent episode of rapidly expanding niche breadth.
2. Material and methods
(a) Toad data collection
We used cane toad (Chaunus [Bufo] marinus) distribution data assembled from the Queensland Museum and several published sources (Phillips & Shine 2004) and more recent records from Estoup et al. (2004), the Northern Territory Department of Natural Resources, Environment and the Arts and from New South Wales National Parks and Wildlife Service and Forests New South Wales. We aggregated toad data into minute-by-minute grid cells (area∼2.7 km2), the finest resolution possible given the positional accuracy. Toad presences and absences were assigned to each grid cell if their point locations fell within the cell's boundaries. Coincident presences and absences were recorded as a presence. This resolution resulted in a spatial database of 1366 presences and 486 absences that were each referenced by minute grid cell within continental Australia. The area occupied by the current cane toad range, and that predicted by an earlier model based on the cane toad's native distribution, was calculated using an adaptive kernel density estimator (Worton 1989; see electronic supplementary material for details).
(b) Model derivation
We limited our analysis to a small set of predictor variables selected using prior knowledge of cane toad physiological tolerances, habitat requirements, distribution limits in their native range and anthropogenic elements potentially linked with their spread (Zug & Zug 1979; Floyd 1985; Sutherst et al. 1995; Lever 2001; Estoup et al. 2004). We expected that annual temperature (minimum and maximum and their squared terms), annual precipitation, minimum moisture index, mean annual evaporation, elevation, topographical heterogeneity, road cover and urban land use would shape the distribution of cane toads in Australia (see electronic supplementary material for details of predictor variable estimation).
We used logistic regression to model the relationship between predictor variables and toad occurrences in minute-by-minute grid cells. We then applied the Akaike information criterion (AIC) to calculate a weighted log-likelihood distance between candidate models and data. AIC was deemed appropriate because no evidence was found for overdispersion (c=0.61) and the model had a high data : predictor ratio (Burnham & Anderson 2002). We calculated AIC for all subsets of variables in a program written in Matlab v. 7. Since several models performed relatively well with respect to their AIC values, we chose to employ a model-averaging approach (Burnham & Anderson 2002). We calculated the average regression coefficients from the 20 models that comprised the class of models aggregating 95% of summed Akaike weights. These coefficients were calculated as the average of all regression coefficients, substituting zeros for those models for which the particular variable was excluded and weighted by the Akaike weights recalculated for the model subset (Burnham & Anderson 2002). Selected models generally included some combination of 9 or 10 of the 11 variables. Proportion paved roads (found in 35% of models) and precipitation×topographical variation (60%) were selected the least often. The overall model performed well as assessed by regression statistics (N=1852; model Χ2 p<0.0001), the area under the curve of the receiver operating characteristic plot (A=0.92) and the total percentage of correct presence/absence assignments (86.7%). In addition, model predictions were insensitive to the effects of spatial autocorrelation (see electronic supplementary material).
The estimated regression model then was applied to the variables in each minute grid cell on the continent to derive a suitability index for Australian cane toads. A threshold value above which a toad presence is predicted was derived parametrically from the maximum-likelihood classification of the binormal distribution of predictions for toad occurrences. To accomplish this aim, data were converted to normality using the logit transformation (Robert et al. 1991). This threshold was then weighted by the ratio of presences to absences in the dataset while assuming an equal cost to predicting a false positive or negative. Such thresholds derived from data-based methods tend to perform better than subjective thresholds in classifications (Liu et al. 2005). We then converted the habitat suitability map into projections of cane toad presence or absence by indicating the potential for invasion when the calculated suitability index met or exceeded the threshold.
(c) Dynamic updating of model predictions
Here, we evaluated the area of predicted suitable habitat in relation to the time of data collection. We were limited to the period following 1974, because this is the first date in which systematic surveys reported absences. We applied the same methods as for the full model to each dataset containing toad presences and absences prior to 1975, 1985 and 1995. Models created using these temporal subsets of the data included a recalculation of regression parameters from the 95% set of best AIC models and the application of a threshold adjusted to the distribution of prediction values (range 0.518–0.527).
We then used a classification tree to differentiate toad presences recorded before and after 1994 based on model parameters. Classification trees offer a flexible solution to separating complex data into a minimal set of groups based on objective criteria (De'ath & Fabricius 2000). We then used 10-fold cross-validation to pick the tree size that was within 1 s.d. of the minimum classification error, as obtained in the subroutine RPART in S-plus v. 7 (Insightful Corp., Seattle, WA; Therneau & Atkinson 1997).
(a) An increase in suitable invasive range
Our model suggests that, unless interrupted, cane toads eventually will occupy over 2 million km2 of Australia (figure 1b). Moreover, cane toads are predicted to colonize 76% of the Australian coastline. This projection greatly increases the predicted invasion area compared with a prior model (Sutherst et al. 1995) and includes for the first time substantial regions of Western Australia, South Australia and Victoria. If this range is realized, cane toads would approximately double their current range and triple projections from a model parametrized on their native distribution (Sutherst et al. 1995). Suitable habitat for cane toad colonization is predicted in regions with minimum annual temperatures of less than 5.0°C (e.g. parts of southeastern Australia) and in regions with maximum annual temperatures of greater than 37.0°C (e.g. parts of the Northern Territory). These findings suggest that cane toads are increasingly occupying zones outside their native bioclimatic envelope.
In addition, we found a positive link between cane toads and proportional urban area (N=1852; odds ratio=18.1; P=3.1×10−11). Cane toads appear to be most successful in open habitats associated with human disturbances, such as roadsides and suburban developments (Lever 2001). The statistical model supported this observation and indicated that most large coastal Australian cities are susceptible to cane toad colonization and hence make ideal targets for control efforts.
(b) The dynamics of invasive range expansion
Evidence suggests that the recent colonization of hot savannah regions in the Northern Territory has been associated with an expanded climatic envelope. To understand this expansion and its effect on model predictions in greater detail, we split the data into nested subsets of information available prior to 1975, 1985 or 1995. We then derived models for each subset using the same methods as for the complete dataset. Models generated with different temporal resolutions of data suggested that net suitable toad habitat changed little until a substantial increase after 1994 (figure 2a).
To determine what aspects of the toad biology may have changed during this expansion, we used a classification tree of bioclimatic variables to split toad presences into homogeneous groups based on whether they were recorded prior to or after 1994. This analysis resulted in a tree with two nodes, demonstrating that recent expansions have occurred in areas where the minimum and maximum annual temperatures exceed 12.1 and 37.0°C, respectively (figure 3). In line with this result, the 95th percentile maximum temperatures associated with all toad presences increased from 36.1 to 38.3°C from 1974 to 2004 (figure 2b). However, the 95th percentile minimum annual temperature of colonized habitats showed no significant change over the same period.
(a) Increasing invasibility of extreme regions
Our results suggest an expansion of suitable habitat area for invasive cane toads in Australia. We predict that cane toads now have the potential to inhabit over 2 million km2 of the continent. This estimate includes three-quarters of Australia's coast, a region where most of the continent's human population and biological diversity are concentrated. If the cane toad's advance continues, this prolific and problematic species is likely to cause further harm to Australia's unique wildlife and economy.
Our model differed from prior work in a variety of ways, including its higher spatial resolution, statistical derivation and the inclusion of anthropogenic variables. However, the predicted distribution of cane toads and the total area of suitable habitat were relatively insensitive to the removal of anthropogenic factors (see electronic supplementary material). The most important variables in our model, including maximum and minimum temperatures, precipitation, moisture index and evaporation, were similar to those used to build a predictive model of Australian toad range based on the species' native range (Sutherst et al. 1995). We also applied approximately the same criterion (greater than or equal to 50% probability of colonization) to determine the eventual range size predicted by native and introduced populations. Hence, we believe that the major difference between our model and past work, and the primary reason for the observed increase in the predicted area of suitable toad habitat, is that our model does not depend on data from the cane toad's native range or its native physiology in that range. Instead, the invasive range model incorporates changes in the toad's bioclimatic niche that have arisen in the exotic habitat.
Our invasive range model predicts a broader distribution of toads in both colder and hotter climates in Australia when compared with those based on native range (Sutherst et al. 1995). Australian cane toads now occupy regions where the minimum monthly temperature falls below 5.0°C and the maximum monthly temperature climbs above 37.0°C. Hence, Australian cane toads now are predicted to be able to inhabit a broad region of cooler climates in southern Australia (figure 1b). At the same time, the toads' recent and successful invasion into high-temperature habitats of the Northern Territory suggests that cane toads are capable of surviving in extensive areas of northern and western Australia. Since we based our climate data on monthly means, daily thermal fluctuations will probably be more extreme and might expose toads to their lower or upper physiological temperature limits in these regions (Floyd 1983). In addition, in hot regions, increases in potential water loss correlated with high temperatures are likely to elevate the potential for lethal desiccation stress (Zug & Zug 1979). Although broad-scale climatic variation will not always equate with the range of available thermal microhabitats on the ground, our study suggests that toads increasingly have colonized areas where thermal and desiccation stresses had been expected to prevent their successful establishment.
(b) Invasion dynamics
By updating our model with time-limited data subsets, we constructed a time-series of invasive range predictions. From this series, we discovered a remarkable increase in the predicted area of suitable habitat that has occurred over the last 10 years (figure 2). This expansion in suitable habitat followed a longer-term rise in the 95th percentile maximum temperature of colonized toad habitat from 36.1 to 38.3°C. Moreover, colonized regions after 1994 were more likely to occur where maximum annual temperatures rose above 37.0°C, according to a classification tree analysis. This pattern cannot be explained by a simple time lag following introduction before cane toads encountered extreme temperatures, because comparably hot and dry regions found in interior Australia halted westward invasions prior to their expansion into the Northern Territory. In contrast to maximum temperatures, 95th percentile minimum annual temperatures did not change significantly over the same period.
These results suggest two points. First, the dynamics of invasion into hot and cold regions of Australia appears to be quite different. Although toads inhabit cooler regions of Australia than expected based on their native range, this cold tolerance has not expanded greatly over the last 30 years. In contrast, cane toads rapidly expanded their range into hot regions in the last decade. This rapid expansion was preceded by the colonization of climatic regions close to the toad's known upper physiological limits of maximum temperature. This pattern is in agreement with the empirical data demonstrating that cane toads from the warmer northern invasion front are expanding their range at increasing rates (Phillips et al. 2006). Moreover, these recent toad expansions suggest that model predictions may be conservative for regions characterized by high maximum annual temperatures.
Clearly, the breadth of predictors associated with an expanding invasion range can be expected to increase until the invasive taxon reaches the edge of its bioclimatic envelope in a new region (Wiens & Graham 2005; Facon et al. 2006). However, several observations suggest that the cane toad's niche has expanded in Australia. (i) Range expansions have been restricted to specific geographical regions (i.e. northern expansion front) rather than to broad regions with similar bioclimatic conditions. (ii) Australian toad range boundaries now lie outside those predicted by their native range. (iii) Toad range expansions have accelerated rather than decelerated with time.
Both ecological and evolutionary explanations can be offered to explain niche expansions in introduced species (Lee 2002; Torchin et al. 2003; Holt et al. 2005). Ecologists have argued that a species' native niche, as defined by its ancestral range climate envelope, may not equate with its realized niche in an exotic habitat because the negative demographic effects of native species are absent (Holt et al. 2005). The successful establishment of introduced species has often been attributed to their ecological release from native enemies and competitors (Diamond 1970; Keane & Crawley 2002; Torchin et al. 2003). Ecological release could result in a higher tolerance of extreme abiotic conditions in novel habitats if both biotic interactions and abiotic conditions combine to limit fitness in the ancestral habitat (Holt et al. 2005). For example, abundant ticks are thought to increase cane toad susceptibility to dehydration and starvation in their native range (Zug & Zug 1979). The cane toad left most of its parasites, pathogens and predators behind, when it was introduced to exotic habitats (Speare 1990). Hence, ecological release then may explain part of the cane toad's success in hot and cold areas of Australia. However, evidence for ecological release remains equivocal (Keane & Crawley 2002; Colautti et al. 2004) and little is known about how biotic interactions and abiotic conditions jointly affect cane toad demography in Australia.
We suggest that an evolutionary, in addition to an ecological, explanation may be necessary to account for the changing distribution of Australian toads. Cane toads increasingly occupy regions in Australia with extreme temperatures, with daily temperatures potentially falling outside experimentally determined physiological limits (Floyd 1983). The number of cane toads in Australia probably has surpassed the population genetic threshold at which mutation rate constrains emergence of new genetic variation (Butin et al. 2005). The observed lag in range expansion followed by explosive growth (Phillips et al. 2006) corresponds with predicted patterns of expansion following niche evolution (Holt et al. 2005). Most importantly, recent evidence suggests that toads from the invasion front in the Northern Territory have evolved enhanced dispersal ability (Phillips et al. 2006). Adequate genetic variation also may be expected for other traits, including those involved in overcoming physiological limits. Toads may be evolving to move not only faster, but also further into extreme abiotic habitats in Australia.
Invasion dynamics usually is viewed as constituting the ecological processes of arrival, establishment and spread (Lockwood et al. 2005). These steps generally ignore the possibility of adaptive evolution. Mounting evidence shows that invaders often evolve rapidly in introduced habitats (Williams & Moore 1989; Reznick & Ghalambor 2001; Lee 2002; Gilchrist & Huey 2004; Butin et al. 2005; Huey et al. 2005). Despite an initial genetic bottleneck, genetic variation can arise in invasive populations through hybrid origins, recombination or mutation following accretion of a sufficiently large number of individuals (Kolbe et al. 2004; Butin et al. 2005). From a theoretical perspective, high gene flow could curb further expansion due to maladaptation in peripheral populations (Kirkpatrick & Barton 1997; Case & Taper 2000; Holt 2003). However, theoretical predictions also suggest that shallow ecological gradients or genetic changes can foster range expansion, despite the constraints imposed by moderate gene flow (Kirkpatrick & Barton 1997).
Therefore, invasion dynamics may include a fourth stage (after arrival, establishment and spread): once adaptation and niche expansion are sufficient, adaptive genetic variation accrues in an expanding invasive population (Lee 2002). Ecological considerations are expected to dominate early on in a successful invasion as the newly introduced species becomes established in a habitat that matches its ancestral niche. Following establishment of a stable population source—a beachhead of sorts—an increasing number of propagules can disperse into suboptimal habitats with little impact on the demography of the source population (Holt et al. 2005). In the meantime, the rapidly growing core population may restore the genetic variation lost through the initial bottleneck, via recombination and mutation (Butin et al. 2005). Together, these factors can increase the probability that advantageous adaptations arise in marginal habitats such that range expansion can occur through niche evolution. Demographic or genetic traits that allow species to reach the evolutionary invasion stage may explain why some taxa become globally invasive while others remain infrequent intruders. For example, cane toads are explosive breeders and often are found at extremely high densities in their introduced habitats (Lever 2001). These demographic characteristics both support high rates of population increase and influence the rate at which successful mutations arise in new habitats.
From a practical standpoint, forecasts of species distributions, such as those applied to invasive species or to species distributions following climate change, should incorporate the possibility that the niches of invasive species may not be conserved in novel habitats. This may occur due to either ecological release or evolutionary dynamics. More accurate forecasts can be accomplished using flexible models that are based on the introduced, in addition to the native, range and that are dynamically updated during a species' range expansion. Unlike the standard approach of using native range to predict the potential distribution of invasive species in novel habitats, models based on the invasive range include the potential for niche expansion through either ecological or evolutionary means (Lee 2002; Kolbe et al. 2004; Butin et al. 2005; Holt et al. 2005; Facon et al. 2006). Range predictions through time can indicate areas of model uncertainty, highlight specific variables to which organisms may be evolving higher tolerances and delineate particular geographical regions where trait divergence is suspected.
Controlling the further expansion of cane toads in Australia remains a daunting task. Future range expansions could devastate Australia's endemic species. However, effective quarantine might still limit the range of this invasive taxon. Model predictions can facilitate attempts to curb further expansion, by focusing monitoring and eradication in cities, at habitat bottlenecks, and between isolated regions of suitable habitat in western parts of Australia. However, if the continuing acceleration of toad expansion reflects continuing adaptation to a new environment (Phillips et al. 2006), then even our current forecast may prove conservative.
The research was supported by a grant from the Yale Center for International and Area Studies, and data accumulation by the Australian Research Council. Manuscript preparation was completed while MU was a postdoctoral fellow at the National Center for Ecological Analysis and Synthesis, a centre funded by NSF (grant no. DEB-0553768), the University of California, Santa Barbara, and the State of California. We thank S. Webster, F. Lemckert and A. Estoup for providing additional data. J. Brownstein provided advice on modelling aspects. M. Baskett, K. Freidenburg, J. MacLellan and M. Smith provided helpful comments.