Spread risk assessment of invasive axis deer using bioclimatic niche models

Ecological niche models are useful tools for predicting and managing biological invasions. However, their reliable implementation and derived applications could be biased by niche shifts among native and invaded areas and by non-equilibrium dynamics. We calculated bioclimatic niche overlap between native and invaded ranges, and tuned species distribution models trained on both the native range (native model), and the native and invasive ranges (complete model) for the axis deer (Axis axis) to assess its invasion risk potential worldwide. Our results suggest that the axis deer shifted its realized niche in some invaded areas where most of the climatic niche available was contained within the environmental conditions from the native range. We found signals of non-equilibrium dynamics between native and invaded niches due to several unfilled areas in invaded regions. The native and complete models are similar in their projection for the equatorial region but differ significantly in the tropical regions. The first fails to predict areas where the species was successfully established. This difference reinforces the idea of using all the occurrences for predicting possible invaded areas, and that invasive species can somehow shift their niche. The complete model predicted that the most climatically suitable regions for the species worldwide were in tropical and subtropical climates, particularly the Neotropics and Afrotropics. While the species has been already established in the Neotropics, it hasn’t in the latter. Therefore, our work emphasized the importance of a nuanced understanding of the potential distribution of invasive species to inform more effective management measures.

methods rely on fundamental assumptions such as niche conservatism (i.e., niche equivalence), the tendency of a species to maintain ancestral ecological requirements (Wiens and Graham 2005;Soberón 2007;, and distributional equilibrium within the current climate (Elith et al. 2006). Under these hypotheses, invasive species tend to occupy a niche equivalent to the realized niche within their native range during the invasion process Liu et al. 2020). ENMs made with climatic predictors estimate an approximation of the realized niche since observed occurrences are readily constrained by biotic interactions, dispersal limitation, and habitat availability and quality (Soberón 2007;Guisan et al. 2014).
Projections of predicted environmental suitability onto geographic space (i.e., ENM) for invaded regions are also challenged by recurrent evidence of niche shifts after species introductions (Guisan et al. 2014). Niche shifts can be caused by changes in the species' realized niche due to differences in biotic interactions, the effect of novel habitat types, including their structure and availability (DeWalt et al. 2004;Broennimann et al. 2007), phenotypic plasticity (Dlugosch and Parker 2008;Davidson et al. 2011), or an evolutionary shift in the fundamental niche (James and Abbott 2005). Niche shifts after invasion have been documented in several plant and animal species and remain an important and controversial topic for the reliability of inferences based on ENMs (Sales et al. 2021). Part of the disagreement is likely interpretative, yet some authors have identified issues concerning study design, methods, and data types (West-Eberhard 2008;Barbet-Massin et al. 2018;Pili et al. 2020;Liu et al. 2020).
The Cervidae family has the highest proportion of successful invaders among mammals (Clout and Russell 2008). Negative economic and environmental effects have been associated with these introductions worldwide (Conover 1997;Reimoser 2003;Dolman and Wäber 2008;Davis et al. 2016). The axis deer (Axis axis Erxleben 1777) is the most successful invader among the deer species worldwide, representing a concern both for economic and environmental reasons, due to possible competition with livestock and as a vector of disease (Riemann et al. 1979;Richardson and Demarais 1992;Mertins et al. 2011;Debárbora et al. 2012;Davis et al. 2016), and also due to interspecific competition with native fauna (Ali 2004;Faas and Weckerly 2010;Mohanty et al. 2016), respectively. The species is native to the Indian subcontinent (Grubb 2005). However, the axis deer was introduced for hunting purposes -with variable rates of success -in several biogeographical regions and it is currently expanding in many areas (Nentwig 2007;Clout and Russell 2008;Sponchiado et al. 2011).
Understanding the ecological needs of invasive species is fundamental for their efficient management (Stohlgren and Schnase 2006;Tobin 2018). In the initial phase of the invasion, management should aim at limiting the spread potential (Simberloff 2009). To create a more reliable ENM, it is important to assess whether the species is able to occur in bioclimatic areas outside of those present within the native range or not (Elith et al. 2010;Pili et al. 2020). In this study, we aimed at determining whether the realized climatic niches for the axis deer differ between native and introduced ranges by means of fitting ENMs. We also estimated the extent of the potential climatically suitable areas worldwide to help wildlife managers in the decision-making regarding the control of this species. For this, we tested the reciprocal hypotheses of niche shift and conservatism to have a better understanding of the potential areas of invasion. Being a generalist mammal with documented successful introductions in various parts of the world, we expected that the species could occupy climatically different areas from those found in its native distribution.

Studied species
The axis deer is a mid-sized deer with a shoulder height of 90 cm on average. The species is predominantly a grazer, but also consumes leaves, flowers, fruits, mushrooms, and bark (Sankar 1994;Raman et al. 1996;Sankar and Acharya 2004;Moe and Wegge 2011). It has also been seen feeding on garbage deposits in peri-urban areas (Raman et al. 1996), which hints at its wide trophic niche and adaptation potential.
The axis deer is native to India, Sri Lanka, Nepal, Bhutan, and Bangladesh, where it inhabits transitional areas between forest and grassland (Grubb 2005). For hunting purposes, the first records of introductions outside the native range took place Spread risk assessment of invasive axis deer using bioclimatic niche models in Australia, in a park near Sydney in 1800 (Vos et al. 1956;Long 2005). Within a decade, individuals were already established and dispersing. From these individuals, new introductions took place in various regions of Queensland in 1872, where populations can still be found today (Vos et al. 1956;Long 2005). In the United States, the axis deer was introduced in Texas in 1932, in Hawaii in 1867, in Florida in 1930, and in California in 1948(Vos et al. 1956Long 2005). In 1979, the axis deer occurred in 20 counties in Texas, with an estimated population of over 7800 individuals (Lever 1985). In Florida, however, populations have not further expanded, remaining restricted to the counties to which they were initially introduced, while in California, due to hunting control initiatives, there are fewer than 350 individuals (Long 2005). In the Hawaiian Islands, axis deer populations expanded to more than 1700 individuals on Lania and 3000 individuals on Molokai (Long 2005). There are also records of introductions in France and England, but the populations did not persist, probably due to the low number of introduced individuals (Long 2005). For Mexico the species has been reported to be present in some protected areas (Álvarez-Romero et al. 2008).
Most of the recent records of axis deer outside its native range come from South America. In Uruguay for example, the first individuals were introduced in 1930 in the department of Colonia, where they became established and subsequently expanded across the country (Pereira-Garbero et al. 2013). In Argentina, they were introduced in the provinces of La Pampa in 1906 and Buenos Aires in 1930, and later in Santa Fe, Neuquén, and Río Negro, where the species is established in wild populations and expanded its distribution throughout the country (Novillo and Ojeda 2008). In Brazil, the species was first recorded in the southernmost region of the country, the state of Rio Grande do Sul, in 2009, in the Espinilho State Park. The area is close to the border with Uruguay and Argentina, where the incursion into Brazil presumably originated (Sponchiado et al. 2011). Since the first record, several other sightings have been reported, suggesting imminent establishment and expansion throughout southern Brazil. Currently, it is on the list of exotic species of Rio Grande do Sul (Portaria Sema nº 79/2013), and its monitoring is aided by citizens through an application (InvasorasRS) developed by the environment secretariat and launched in late 2018.

Species occurrence
We extracted 1,134 georeferenced records within the species' native range from the Global Biodiversity Information Facility (992 records, https:// www. gbif. org) and India Biodiversity Portal (142 records, https:// india biodi versi ty. org). An IUCN range map was used to restrict the countries where the species is considered native. With regards to the invaded regions, we extracted 492 records from GBIF and 353 records from gray and scientific literature, camera trapping, and sightings reported by experts.
We removed duplicates and checked all occurrences individually to guarantee that there were no location errors, excluding occurrences located in water bodies, urban perimeters, or those that did not match descriptions in the databases. We used only records of confirmed invasive populations to reduce bias, checking available literature and information in databases. To minimize the effects of sampling bias, we randomly removed occurrences located less than 5 km apart (Phillips et al. 2009; Moe and Wegge 2011) using the spThin package (Aiello-Lammens et al. 2015) in the software R v4.1.1. Following these procedures, we ended up with 328 records from the native range, 234 from South America, 157 from Texas, 24 from Australia, and 13 from Hawaii ( Fig. 1; Supplementary Material).

Bioclimatic variables
We downloaded 19 bioclimatic variables for the globe from WorldClim v2 at 2.5 arc-minutes resolution (Fick and Hijmans 2017), from which we extracted bioclimatic values for axis deer occurrences from both native and invaded ranges. We aimed to screen broadly for climatically suitable areas across the world's main geographic regions. To avoid multicollinearity and overfitting the models, we used the usdm package (Naimi et al. 2014) to select variables with variance inflation factor (VIF) < 10 (Guisan and Thuiller 2005;Vittinghoff et al. 2012), while also only retaining variables with |r| < 0.75 (Pearson's correlation coefficient) (Butler et al. 2012;Silva et al. 2016;Mohammadi et al. 2019;Herrando-Moraira et al. 2019). Upon this routine, we retained nine variables for ENM and niche overlap analyses (Table 1; Supplementary Material).

Niche overlap
We quantified niche overlap between the niches occupied in native and invaded ranges in the multivariate climatic space derived from a principal component analysis (PCA) of the nine selected variables following Broennimann et al. (2012). We selected the environmental backgrounds for each area (i.e., calibration area) from a minimum convex polygon (MCP) with a 20 km radius buffer centered at the occurrence locations. We extracted the values of the climatic variables to construct the available environmental space represented by the two principal axes that accounted for the maximum bioclimatic variation within the species' native and invaded ranges (Silva et al. 2016;Herrando-Moraira et al. 2019). We calculated Schoener's D to estimate the niche overlap (Schoener 1970) and performed tests of niche equivalency and similarity between native and invaded ranges to check for non-random patterns (Warren et al. 2008), to evaluate lower niche equivalencies/similarities than expected by chance. To assess statistical significance, we ran both tests with 10,000 iterations. To understand the nature of possible niche shifts in space, we tested for stability, expansion, and unfilling patterns. We employed niche stability as a measure of the climatic  Precipitation of coldest quarter X Spread risk assessment of invasive axis deer using bioclimatic niche models conditions of the regions where the species inhabits that are shared between the two ranges. Niche expansion describes new climate conditions occupied by the species in the invaded area, whereas niche unfilling refers to those climate conditions occupied in the native range that are not occupied in the invaded area. We performed all niche overlap analyses in the ecospat R package (Di Cola et al. 2017).

Niche modeling
We created two separate bioclimatic suitability models for the species, one using only native occurrence records (hereafter "Native Model"), and the other one using both native and invaded range occurrences (hereafter "Complete Model"). We fitted both models using the presence-background algorithm Maxent v.3.4.4 (Phillips et al. 2004(Phillips et al. , 2017, which has shown good performance across diverse datasets (Elith et al. 2006;Heikkinen et al. 2006;Wisz et al. 2008;Renner and Warton 2013;Valavi et al. 2021).
We assessed different combinations of parameters aiming at the most parsimonious models using the ENMeval v2.0 package (Kass et al. 2021). We tried values of the regularization multiplier from one to five at increments of 0.5 and all feature class transformations, those being linear (L), quadratic (Q), hinge (H), product (P), and threshold (T). Regularization multipliers penalize model complexity with increasing values, and feature classes add to the complexity of the model; thus, tuning both hyperparameters can help you find optimal model settings (Phillips et al. 2006;Phillips and Dudík 2008). We ranked models based on their ΔAICc, which represents the difference in Akaike's information criterion (AICc) units between the best and subsequent models. Models with a low rank according to AICc are therefore the ones with the best performance (Warren and Seifert 2011). AICc trades off explanatory power with model complexity; e.g., most parsimonious models explaining more variation have the lowest AICc values (Warren and Seifert 2011; Ficetola et al. 2014). Additionally, we calculated the area under the curve (AUC) on withheld data using 10-fold random cross-validation and considered it alongside AIC to determine model performance.
The AUC is a threshold-independent performance measure that reflects the probability that a randomly chosen presence site will rank above a randomly chosen background site (Phillips et al. 2016). Values near one indicate high discriminatory ability, whereas values of 0.5 (or less) indicate discrimination no better than random (Elith et al. 2006). However, the efficiency of this metric has been much debated (Lobo et al. 2008, Jiménez-Valverde et al. 2012. We parametrized the optimal, final Native and Complete models according to the parameters indicated by the previous selection with the following adjustments. For both models, we used 10,000 background points, extracted from the same MCP from PCA analyses. For the Native Model, only the MCP area from the native region was used, whereas, for the Complete Model, all MCP were used as background. The whole background points data set represented the calibration area, hypothesized to be the available area and bioclimatic conditions that were accessible to the species over evolutionary time (Barve et al. 2011), assuming that the species is at equilibrium with the environment at least in the areas between the occurrence points. We chose to use minimal convex polygons due to the density and distribution of the points, the dispersal capacity of the species, its home range size, and the success of that approach in studies of niche shift with other species (Soberón and Peterson 2005;Warren et al. 2008;Barve et al. 2011;Herrando-Moraira et al. 2019).
We increased the maximum iterations to 5000, providing the model with adequate time for convergence. For model projection onto geographical space, we used the option "fade by clamping," which reduces prediction at each point in projections by the difference between clamped and non-clamped outputs at that point (Phillips et al. 2009), where clamped means that the values of the projection raster are constrained between the maximum and minimum values associated with the training data (Elith et al. 2010).
To evaluate the final model predictions, we calculated the true skill statistic (TSS), a measure of model performance that corrects for the dependence on the prevalence of the data being modeled, using a threshold value. The TSS is the measure of choice for presence-absence predictions and ranges from − 1 to 1 (Allouche et al. 2006). Values close to 1 indicate good prediction and values equal to or smaller than zero are not better than random predictions (Allouche et al. 2006). For this analysis, we set maximum sensitivity plus specificity as the threshold as recommended by Liu et al. (2013Liu et al. ( , 2015. We also calculated the variable importance using the permutation importance routine provided by Maxent (Phillips et al. 2017).
For a better understanding of possible expansion routes, we imported the model prediction raster (cloglog transformation; Phillips et al. 2017) into QGIS 3.16.10 and regrouped habitat suitability scores in five classes as follows: unsuitable (0-0.2); barely suitable (0.2-0.4); suitable (0.4-0.6); highly suitable (0.6-0.8); and optimal habitat suitability (0.8-1.0). Further, we used the maximum sensitivity and specificity values (Native Model − 0.347, Complete Model − 0.411) as the threshold to convert the prediction into binary projections (Liu et al. 2013(Liu et al. , 2015. In addition, we overlayed model projections with the world's biomes (Dinerstein et al. 2017) and analyzed which of these would present a greater portion of suitable bioclimatic areas for the species.

Results
The first two axes of the PCA of environmental variables cumulatively explained 57.2% of the variation associated with the observed global distribution for the axis deer (Fig. 2D). The niche equivalency analysis revealed significant differences between the bioclimatic niche occupied by the axis deer in the native versus that of the invaded range in Texas, Australia, and South America (Table 2; Fig. 2a). Niche similarity yielded non-significant results for all invaded ranges (Table 2; Fig. 2b-c).
Regarding niche characteristics, in Texas and South America, following the invasion, the axis deer shifted from the realized niche envelope of the native range, with low stability and substantial expansion signal (Table 2). Australia showed intermediate niche stability and expansion, and Hawaii shows high stability and a low expansion. All invaded areas showed a high availability of native climatic niche considering Fig. 2 Climatic niche representation for the axis deer (Axis axis) from niche overlap analysis. The color gradients indicate the density of occurrences for the species (Green-native range; Red-Hawaii; Yellow-Australia; Black-South America; Blue-Texas). The solid lines in a indicate 95% of the occurrence, while lines in b depict 95% of the available environment from the selected background for each area. In c, 95% of the occurrence density is illustrated with thin lines, and 100% of available climatic background with thick lines. Vectors onto the biplot d show the variable contributions in direction and magnitude with respect to the two main axes from the PCA ordination; explained variation (%) for each axis is given in parenthesis the occupied climates in native areas (i.e., high niche unfilling).
The selection of parameter combinations for the Maxent models yielded the feature classes LQH and regularization multiplier 1.5 for the Native Model, and feature classes LQPH with regularization multiplier 1 for the Complete Model. Both final models showed good predictive capacity (i.e., Native Model: AUC v = 0.85, TSS = 0.67; Complete Model: AUC v = 0.87; TSS = 0.68). The Native Model predicted the most suitable regions for the axis deer in tropical climates, including large suitable areas in the center of South America and sub-Saharan Africa, a large portion in Central America, and some regions in Southeastern Asia and Northern Australia (Fig. 3). The Complete Model predicted the most suitable regions in tropical and subtropical climates around the equator and the southern hemisphere, including large suitable areas across central Andes mountains, north, northern, and southeast South America, central sub-Saharan Africa, Madagascar, Caribbean islands in Central America and southeastern Asia, and southeast North America (Fig. 3). The most important variables for the Native Model were isothermality (40.1%), precipitation seasonality (18.3%), and precipitation of the coldest quarter (14.9%). In turn, for the Complete Model, the most important variables were precipitation of driest month (32.2%), precipitation seasonality (16.2%), precipitation of wettest month (16.2%), precipitation of wettest month (13.2%), and precipitation of coldest quarter (11.8% - Table 3).
Concerning biomes, the two of them that contained the greatest amount of projected areas, independently of the model, are tropical and subtropical moist broadleaf forests and tropical and subtropical grasslands, savannas and shrublands (Fig. 4). The Native Model included a small area of the deserts and xeric shrublands biome in Mexico and the native region while the Complete Model presents areas in temperate broadleaf and mixed forests and temperate grasslands, savannas and shrublands, and small portions in temperate conifer forests and montane grasslands and shrublands. Both models overlapped with flooded grasslands and savannas; tropical and subtropical dry broadleaf forests; tropical and subtropical grasslands, savannas and shrublands; and tropical and subtropical moist broadleaf forests.

Discussion
Our modeling of the bioclimatic niche and global potential distribution of the axis deer demonstrated a shift of the species' realized bioclimatic niche between native and invaded areas, matching previous studies (Holt et al. 2005). We showed that vast regions -particularly in South and Central America and Central Africa -were climatically suitable for invasion based on the climatic (realized) niche from regions where it was successfully introduced.
We found a low overlap between the climatic niche of native and invaded regions. In Hawaii, despite appearing to have a greater overlap in the modeled niche (Fig. 1a), the relative low overlap occurs because the modeled niche for Hawaii and that of the native range have overlap mostly in areas with low density for the native range (Fig. 2a). In Texas and South America regions, on the other Table 2 Estimated metrics of bioclimatic niche dynamics comparing the native range with that of known invaded areas around the world for the axis deer (Axis axis) Interpretation of niche parameters: Niche overlap using density grid (D; 0 = no overlap, 1 = complete overlap), equivalency/similarity using density grid (p < 0.05, indicating that niches are less equivalent/similar than expected by chance), expansion (0 = in equilibrium, 1 = fully in expansion), stability (0 = not stabilized, 1 = totally stabilized), unfilling (0 = totally filled, 1 = totally unfilled).  . 3 Global climatic suitability estimates and potential invadable areas for the axis deer (Axis axis). On the left, predictions are based on models calibrated solely with occurrences within the species' native range (i.e., Native Model). On the right, models were calibrated with four invaded locations plus the native range (i.e., Complete Model). On the upper panels, the climatic suitability is divided into five categories, whereas the bottom panels display geographic ranges from binarizing predictions using the maximum sensitivity plus specificity as suitability threshold hand, following the invasion, the axis deer shifted most of the realized niche envelope and exhibited substantial niche expansion, indicating a divergence in the realized niche (Fig. 2a). Still, part of the new climatic conditions from the occupied areas are present in the native region even if the species has no records within those areas (Fig. 2b-c). This means that the axis deer could be unable to occupy areas similar to those in its native region due to other factors, such as the absence of key resources (e.g., water bodies), biotic interactions, and/or anthropogenic pressures (e.g., hunting and habitat degradation) (Madhusudan 2004;Soberón and Peterson 2005;Soberón 2007;Dave 2008;Karanth et al. 2010;. Another possibility is that, in the native area, habitat fragmentation and degradation may have prevented the dispersion and maintenance of populations in different regions where the species was historically present; thus, nowadays, there are no records of occupation within these bioclimatic zones (Vance 1984).
The non-equilibrium between native and invasive niches of the axis deer was suggested by many unfilled areas in the invasive range. Proportionally, large areas climatically similar to the native region are not occupied by the axis deer in these invaded regions. Niche unfilling occurs when some environmental conditions occupied within the native niche are available, but unoccupied outside the native range (Guisan et al. 2014). These unoccupied areas could result from recent colonization and ongoing dispersal (Guisan et al. 2014). However, this is not clear, and more studies specifically focused on these regions are needed.
In regions where the species occupies a distinct climatic niche (e.g., a portion of Texas and South America), phenotypic plasticity and novel genotypes could be part of the answer (Richmond et al. 2005;Dlugosch and Parker 2008;Davidson et al. 2011). However, this is just a hypothesis, and more studies in this area are needed. Another plausible explanation could be the continuous dispersal of individuals from Fig. 4 Overlay of the world's terrestrial biomes with the binarized predictions (i.e., maximum sensitivity plus specificity as suitability threshold) of geographic distribution for axis deer (Axis axis) from models calibrated with the native range (upper panel), and with invaded and native ranges (lower panel). The paler coloration represents biome areas outside of the predicted ranges for the species climatically similar areas to the expanded climatic areas. Thus, in such a case, expanded areas would be maintained due to the number of immigrants since the propagule pressure can be fundamental for the maintenance of invasive populations and ensure the successful invasion of new climatic areas (Simberloff 2009;Bradie and Leung 2015;Hudgins et al. 2017). Currently, there is not enough evidence to assess the relative importance of these possible explanations and grant further investigation.
Our worldwide projections suggest a shift in the species' realized niche. Although the two models showed good predictive performance, both presented important differences when compared to the already established invasive populations. The Native Model was unable to predict areas where currently invasive populations occur (i.e., parts of South America and the United States) but correctly predicted parts of Australia. On the other hand, the Complete Model predicted suitability in some parts of England, where the species was introduced but has never spread, and in some parts of Mexico. This highlights the importance of thinking about what data to use to model species invasions in a new environment or at a more restricted scale (Broennimann and Guisan 2008;Elith et al. 2010;Fernández and Hamilton 2015). Focusing only on the native distribution, both models diverged in their predictions. Although there was some overlap between them, the Native model showed more suitable areas in Central and West India, while the Complete Model focused more on East India. Although historically the axis deer was considered present all across India, that area might overestimate the actual area of occupancy of the species; for example, currently, it is found in greater density in protected and forested areas used as a refuge (Dave 2008). While both model predictions included part of these areas, the degree of suitability changed across models. We believe that the Native Model was more reflective of the current distribution of the species while the Complete Model shows the potential distribution of the species, as it considers a greater diversity of environments.
Regarding potential invasion areas, both models predicted that the major regions suitable for the axis deer are in tropical and subtropical climates. For the Complete Model, in the Americas, the axis deer was introduced mainly in areas with high bioclimatic suitability, suggesting a high risk of further range expansion from those areas. The areas with the greatest suitability in Argentina, southern Brazil, and Uruguay are consistent with the current distribution of the species and reports of new sightings (Long 2005;Novillo and Ojeda 2008;Rosa et al. 2020). Based on our suitability estimates, we found that the species' bioclimatic niche mostly coincides with the climates found within two biomes -the tropical and subtropical moist broadleaf forests and tropical and subtropical grasslands, savannas and shrublands. Roughly, we predicted two coarse major centers of high climatic suitability for the axis deer in the Neotropics. One corresponds to Southern Brazil, Northern Argentina, and Uruguay; the other includes Northeastern Brazil near the coast, and north and northeast South America. The former center is already invaded and shows signs of expansion (Sponchiado et al. 2011). Given our models, further spread of the species over the Neotropics is plausible, considering that potential corridors across regions with lower suitability exist (i.e., forested areas in Central-West Brazil), and there are already some records of the species in the southern regions (Sponchiado et al. 2011). Since the axis deer is a game species, there is also the risk of intentional introduction in the surrounding areas to influence hunting-related laws. Illegal introductions are thus concerning because propagule pressure by multiple introductions is a major cause of the success of introduced species (Simberloff 2009;Lockwood et al. 2013).
Other highly suitable regions for the axis deer, according to the Complete Model, included regions of Central Africa, Southeast Asia, and some regions in Central America and North America. Concerning Central Africa, we have not found any studies showing introductions of the axis deer. There are many game species native to this region, and there has not been an interest in introducing new species for this purpose (Lindsey et al. 2007). In Mexico where the species is reported to be present (Álvarez-Romero et al. 2008), the Native Model predicted the areas with the highest habitat suitability well; in contrast, the Complete Model covers only a portion of these reported areas. We have not found further articles reporting whether populations expanded, were managed, or did not survive, and therefore cannot assess which of the two models best predicts the possible distribution for the area. In the USA, the Native Model didn't predict bioclimatic suitability in some areas where the species is already established. The Complete Model, on the other hand, presents a mosaic of areas with low and high bioclimatic suitability. The variation may explain why the species does not present high densities and has not spread through southern Florida. At the same time, the highest densities are found in Texas, an area in which the model predicted high suitability (Long 2005).
In Southeast Asia, regions of medium to high probabilities were found along the border of the axis deer's original distribution. Competition with Axis porcinus could be a potential explanation for the absence of A. axis in these regions, considering that the limits of their distributions match and that both species share several ecological traits (Timmins et al. 2015).
In Australia, predictions of the Complete Model differed from those proposed by Davis et al. (2016), which yielded greater suitability values for the northern region. Our estimates were much more restricted to the eastern region of the country, areas in which the species is already present. However, the Native Model made similar predictions to that of Davis et al. (2016), albeit more restricted to the coast. Using the most recent occurrence records available in public databases, the Native Model and the model proposed by Davis et al. (2016) do not match the density of occurrences, therefore the Complete Model would be the most accurate. From a biome perspective, the Native Model does not predict the occurrence in the Temperate Broadleaf and Mixed Forests areas, unlike the Complete Model. This difference could be explained by looking at the results found between the niche analyses of the native region versus the Australian region. There are many climatically similar areas between the regions. A larger survey in search of the real distribution of the species in Australia may help to better clarify its potential areas of distribution, as there are currently few records mainly in three states (Queensland, New South Wales, and Victoria).
Although we only used climatic variables in our modeling, we acknowledge that other explanatory factors also influence the distribution of the species, since not all climatically suitable areas possess the other conditions to be occupied (e.g. biotic interactions). Climatic variables address the abiotic components of the niche, so, to advance in the local or regional scale predictors associated with species, biotic limitations and dispersal characteristics must be incorporated into the models. Variables reflecting land use, distance from water bodies, accessibility, and dispersibility at finer scales could help refine model predictions. All modeling is dependent upon the nature of the input data, and this could explain discrepancies with other studies. In addition, we know that the values used for limits in the creation of binary models as well as the evaluation method can be debated and will lead to modification in the projected boundaries. Given the data currently available, we think the estimates we present are reasonable range estimates that can hopefully be improved with continued monitoring.

Conclusion
Biological invasions are major threats to biodiversity. Many vertebrate species are intentionally introduced, and cervids have been particularly successful invaders. The axis deer is a priority species for management since it shows a high success rate in the documented cases of introduction and is currently expanding. In this work, we show that there are climatically susceptible areas for the axis to invade in several continents where it has been either already introduced or remains absent. In addition, we highlighted that the species readily occupied climatic niches beyond those found within the native area, thereby reinforcing the importance of early warning systems and management actions to prevent dispersal and potential settlement of the axis deer into new areas.