Invasive risk assessment and expansion of the realized niche of the Oriental Garden Lizard Calotes versicolor species complex (Daudin, 1802)

Correlative species distribution modelling (SDM) can be a useful tool to quantify a species’ realized niche and to predict its potential distribution for non-native ranges. The agamid lizard Calotes versicolor s.l. belongs to the most widely distributed reptile taxa worldwide. In the past, C. versicolor s.l. has been introduced to several countries, including regions in the Oriental, the Neotropical and the Afrotropical realms, where strong negative impact on the local fauna is assumed. Due to the complicated taxonomy and the existence of several cryptic species, which are covered by this taxon, we used C. versicolor sensu lato and its four subtaxa ( C. versicolor sensu stricto, C. irawadi , C. vultuosus , C. farooqi ) as target species to (1) compute correlative SDMs for C. versicolor s.l. and its subtaxa and project them across the globe to highlight climatically suitable areas of risk for future invasion and (2) based on the ecological niche concept, we investigate if the species complex expanded its realized climatic niche during the invasion process. We use two different SDM approaches, namely n-dimensional hypervolumes and Maxent. N-dimensional hypervolumes are a non-hierarchically ranked approach, which is a useful tool to investigate the expansion in the realized niche, while Maxent, a hierarchically ranked model, is used to focus on potentially suitable areas for future invasion. We calculated two final models for C. versicolor s.l., one based on records from the native range and one based on records from the native and invaded range, as well as one model for each subtaxon. Our results show a geographic expansion into novel climatic conditions as well as an expansion in the realized niche. Our results reveal that C. versicolor s.l. is currently inhabiting 13% of its potential range but could find suitable climatic conditions on a global surface area between 14,025,100 km 2 and 53,142,600 km 2 . Our predictions reveal large areas of highly suitable climatic conditions for the Oriental, Australian, Afrotropical and Neotropical realms, whereas only small regions of the Palearctic and Nearctic realms provide moderately suitable conditions. Further, some localities, especially those with a high amount of human traffic like ports or airports, might act as multiplicators and might therefore be a stepping stone into further areas.


Introduction
Anthropogenic habitat alterations, climate change, and increased human-mediated transport, alter the distribution of organisms, which leads to novel invasions and the expansion of established populations of non-indigenous species (Fordham et al. 2012, Kearney et al. 2008, Todd et al. 2008, Shabani et al. 2020. As a result, non-native taxa spread worldwide (Sala et al. 2000, Shabani et al. 2020. Thus, identifying the limits of an organisms' distribution is particularly important for invasive risk assessments (Kearney et al. 2008, Jimenez-Valverde et al. 2011.
Species distribution modelling (SDM) can be a first step in invasive risk assessments. For this task, the concept of the ecological niche of a species is essential. The ecological niche comprises abiotic (fundamental) and biotic factors (interactions among organisms like competition, parasitism, mutualism, predator-prey-relationship, etc.). The fundamental niche comprises all abiotic factors that are necessary for species' survival and reproduction (Grinnell 1917, Hutchinson 1978, Soberón & Peterson 2005, Soberón 2007). The intersection among biotic and abiotic factors and geographic accessibility determines the species' realized niche (Grinnell 1917, Hutchinson 1978, Soberón & Peterson 2005, Soberón 2007). Climate is one of the main parameters determining the ecological niche and the potential distribution of taxa (Soberón 2007, Thuiller et al. 2004). Therefore, ectothermic organisms such as reptiles are ideal subjects for SDM assessments as they have a strong dependency on climatic conditions. To predict the potential invasive risk of species, correlative SDM approaches linking geographic occurrence data with environmental predictors, can be used (Merow et al. 2013, Guisan & Thuiller 2005. In the past, the realized niche of a species was assumed to be conservative among space and time (Peterson 2011). However, subsequent studies provide a more complex pattern with some degree of flexibility. Shifts in the realized niche were found for plants (Broennimann et al. 2007) and animals (i.e. Medley 2010, Rödder et al. 2017, including reptiles (Nania et al. 2020).
The Oriental Garden Lizard, Calotes versicolor (Daudin, 1802) complex is among the rather widespread agamid lizards. However, the taxon comprises several cryptic species (Zug et al. 2006). The type locality of C. versicolor is unknown and the holotype was lost (Auffenberg & Rehman 1993). A neotype was designated from Pondicherry University Campus, Kalapet, Pondicherry, India (Gowande et al. 2016), which was invalidated one year later related to the complicated taxonomic status of this species complex (Chaitanya et al. 2017). Despite the taxons' wide-ranging distribution across large parts of Asia, in historical times most taxonomists did not recognize other populations to be a different species. This phenomenon may be explained by the seemingly uniform appearance of this taxon and the ease of labelling an already existing scientific name to collected or observed specimens (Zug et al. 2006).
However, some researchers also noticed subtle differences among populations (i.e. Kästle et al. 2013, Radder 2006. As a result, Auffenberg and Rehman (1995) described a Pakistani population as the only currently recognized subspecies, C. v. farooqi, next to the nominal form (Uetz et al. 2021). Since the last two decades, some authors started to particularly shed light into the species' complex phylogenetic relationships (Huang et al. 2013), distributional patterns (Liu et al. 2021) and elevated some populations to distinct species level (Zug et al. 2006). However, the majority of populations is still assigned to C. versicolor sensu lato (below abbreviated s.l.). Recently, Gowande et al. (2021) published a phylogeny based on the mitochondrial 16S and CO1 genes and supported this by morphological data to disentangle the Indian taxa of the C. versicolor complex. According to the latter authors, the species complex contains at least four distinct lineages: (1) C. versicolor s.str. from South India and parts of the Southeast coast, (2) C. irawadi s.l. Zug et al., 2006 from Northeast India across Southeast Asia, (3) C. vultuosus (Harlan, 1825) (resurrected synonym of C. versicolor) from the eastern Indo-Gangetic plains, and (4) C. farooqi Auffenberg & Rehman, 1993 (former subspecies of C. versicolor) from Pakistan and adjacent regions. However, it must be mentioned that at least C. irawadi is still known to be a species complex and also C. farooqi includes another distinct phylogenetic lineage and further surveys are required to clarify these issues. Furthermore, Gowande et al. (2021) found that C. calotes (Linnaeus, 1758) is a phylogenetic sister taxon to C. farooqi and is therefore also included in the C. versicolor complex.
C. versicolor s.l. naturally occurs across large parts of the Oriental realm from eastern Iran (Anderson 1999, Mobaraki et al. 2013, and Afghanistan (Clark et al. 1969), eastwards to China and southwards to the Malay Peninsula (Smith 1935, Auffenberg & Rehman 1993, Zug et al. 2006). In the past, this agamid lizard has been introduced to several countries worldwide where it could establish reproducing populations. These areas include regions in the Oriental (i.e. Das et al. 2008, Pili 2019, the Neotropical (Enge & Krysko 2004) and the Afrotropical realms (i.e. Sandera 2009, Spawls et al. 2018, see Table 1 for an overview). It has a semi-arboreal, heliophile (= "sun loving") lifestyle and prefers open, often anthropogenic or disturbed landscapes like open forest, grassy savanna, scrubland, wasteland, gardens, cemeteries, parks and areas along roadsides (Diong et al. 1994, Erdelen 1984, Stuart 1999, Khan & Mahmoud 2003, Blanchard 2000, Permalnaick et al. 1993. The taxon avoids dense forest with closed canopy (Permalnaick et al. 1993, Diong et al. 1994, Pawar 1999. C. versicolor s.l. occurs in the lowlands (Stuart 1999) but also in much higher altitudes (Kästle et al. 2013, Nanhoe &Ouboter 1987, Srinivasulu et al. 2014, Vassilieva et al. 2016), even up to 3,200 m (Gautam et al. 2020. Studies investigating the lizard's impact on the native biodiversity are limited to very few localities and mostly rely on observations that have not been tested. In Singapore, it was found that C. versicolor s.l. to some extent displaced the native agamid lizard Bronchocela cristatella (Kuhl, 1820) (Diong et. al. 1994). Further, Mauremootoo et al. (2003) stated that, on Mauritius, C. versicolor s.l. competes with endemic geckos by preying on native invertebrates. Vinson (1968) speculated that the absence and decrease in population size of endemic stick insects (Insecta: Phasmatodea) on Mauritius and Réunion is related to the presence of this non-native lizard. It was shown that this lizard is an opportunistic generalist that preys on several kinds of invertebrates as well as small vertebrates and plant material (Matyot 2004). Furthermore, C. versicolor s.l. may act as vector to spread several parasites on native species (Matyot 2004). Concluding these observations and facts provide strong evidence that C. versicolor s.l. may have strong negative impact on the native biodiversity.
The existing literature provides evidence that C. versicolor s.l. could spread by a combination of multiple, natural and manmade, factors and the pathways of introduction seem various. Most likely, C. versicolor s.l. started a natural invasion of the Malay Peninsula in the late 19 th or early 20 th century, as the taxon was not listed in the comprehensive first herpetological literature for this region (Cantor 1847). Later, Boulenger (1912) documented the lizard from the northern portion of the Malay Peninsula, where it was very abundant, while it could not be observed in the southern parts. In the 1980s, C. versicolor s.l. reached Singapore, where it has been established (Lim & Chou 1990, Lim & Lim 1992, Chou 1994, Diong et al. 1994. It was suspected that the species might have benefited from the direct railroad from Thailand and northern Malaysia to Singapore (Tan et al. 2007). In the early 20 th century, a single specimen of C. versicolor s.l. was observed in the very North of Sumatra (De Rooij 1915), from where the lizard started to spread across the region (Brongersma 1931). Das et al. (2008) suggested an introduction from the western coast of the Malay Peninsula, which is located c. 65 km northeast of Sumatra at the closest point of the Strait of Malacca. Recently, new records of the lizard were also documented for Iran (Damadi et al. 2018). For the Philippines, the introduction of C. versicolor s.l. on Luzon as stowaway is assumed as the international airport is in proximity (Emerson 2013). In Oman, the lizard was first recorded in 1982. Its introduction most likely occurred unintentionally as stowaway of shipments of ornamental plants or other nature products (Grossmann & Kowalski 2019).
For the Seychelles, two distinct introduction events were suggested. The first one on Mahé in 1985-86, which has been unsuccessful, and the second one on Ste Anne in 2003, which succeeded in a reproducing and dispersing population (Matyot 2004). Most likely C. versicolor s.l. arrived at the Seychelles by unintentional transport by humans, which is also the case for the Mascarenes (Réunion), where the lizard is guessed to be a stowaway in a shipment of sugarcane cuttings in year 1865 from Java, Indonesia (Permalnaick et al. 1993, Staub 1993. However, the presence of the taxon is not mentioned for Java yet (Uetz et al. 2021).
In Florida (USA) it is known that several C. versicolor s.l. escaped in 1978 from a local reptile dealer. These animals may most likely derive of a shipment from Pakistan (Enge & Krysko 2004, Krysko et al. 2011, where C. versicolor s.l. is the only member of the genus (Auffenberg & Rehman 1993). In Florida, until 2004 the taxon already invaded an area of a maximum extend of 11.7 km east-west and 9.3 km north-south direction (Enge & Krysko 2004, Krysko et al. 2011. Considering this geographic range expansion into non-native regions in combination with negative impact on local biodiversity, a first step for global invasive risk assessment is strongly required for this species complex. Thus, the aims of this study are two-folded: (1) we calculated correlative SDMs for C. versicolor s.l. and projected them across the globe to highlight climatically suitable areas of risk for future invasion and (2) based on the ecological niche concept, we investigated if the species complex expanded its realized climatic niche during the invasion process. To assess the invasive risk of C. versicolor s.l., we computed 8625 Maxent models for each taxon (C. versicolor s.l., C. versicolor s.str., C. irawadi s.l., C. vultuosus, C. farooqi) including 23 regularization multipliers, and 15 feature class combinations, each combination 25 times to find the best fitting model to highlight areas of high climatic suitability. To investigate an expansion in the realized niche, we used an n-dimensional hypervolume approach, which uses support vector machines to assess if geographic occurrences and their linked climatic conditions are in or out of the species' realized niche volume.

Data preparation
Considering the high degree of cryptic diversity and the lack of knowledge on the origin of the introduced specimens, we calculated SDMs for C. versicolor s.l. and added SDMs for the four subtaxa as sensitivity analysis. To predict the invasive risk of C. versicolor s.l., we obtained 1705 occurrence records from GBIF (https://doi.org/10.15468/dl.axawsq) covering the native and invasive range of the taxon based on preserved specimens only. Further, we georeferenced six occurrence records from the literature (see Table 1) using GEOLocate (https://www.geo-locate.org/) to increase the number of records for the invasive range. Furthermore, we used 93 genetically confirmed records from Gowande et al. (2021) to increase the number of records. After checking the dataset for outliers in ArcMap (ESRI 2020)(see distribution in introduction), we corrected the dataset for potential sampling bias (see systematic sampling in Fourcade et al. 2014) and spatial autocorrelation using the raster package (Hijmans 2016) and the ecospat.mantel.correlogram() function in the ecospat package (Di Cola et al. 2017) for R (R Core Team 2021). For systematic sampling, occurrence records in a defined radius are subsampled to reveal a regular distribution of records in geographic space. Subsampling reduces spatial aggregation of records but does not correct low sampling efforts in some areas. This approach underestimates highly suitable areas, where the number of sampled occurrence records reflects the real abundance of a species (Fourcade et al. 2014). Considering the ecospat. mantel.correlogram(), we subsampled the dataset and deleted neighbouring occurrence points in a radius of 50 km to reduce the mantel r value to zero. Finally, we used 195 occurrence records from the native range and 209 records from the native and invaded range for modelling (see Supplemental Material S1).
As Gowande et al. (2021) recently split C. versicolor s.l. into four distinct subtaxa, we conducted a sensitivity analysis and also calculated SDMs for each of those subtaxa. For this purpose, we split the mentioned occurrence records and assigned them to one of the four species based on geographic proximity. We only used occurrence records that could clearly be assigned to one of the four species and removed uncertain ones that were found in transition zones. Finally, we used 27 occurrences for C. versicolor s.str., 97 for C. irawadi s.l., 56 for C. vultuosus and 15 records for C. farooqi. Despite, Gowande et al. (2021) found C. calotes to be a phylogenetic sister species to C. farooqi, we did not include the taxon into our sensitivity analysis as it is morphologically easily distinguishable from C. versicolor s.l., it has never been synonymised with the latter taxon and it seems not to play a role as non-native species yet (Uetz et al. 2022).
We obtained 19 bioclimatic variables from Worldclim 2.0, containing climatic conditions from 1970-2000, at a resolution of 2.5 arc minutes (Fick & Hijmans 2017). We restricted the environmental variables to those reflecting thermal energy and water availability for an ectothermic organism (i.e. minimal, maximal and mean values at the species' records and rejected derived variables [mean diurnal range (bio2), isothermality (bio3), temperature seasonality (bio4), precipitation seasonality (bio15)]. As correlative SDMs are sensitive to multi-co-linearity of predicting variables, we computed variance inflation factor (VIF) between all possible pairs of predictors across the species' background area (see below for explanation) using the vifstep() function of the usdm package for R (Naimi et al. 2014). We restricted the predictor pairs to only one variable, whenever VIF exceeds 10.
For the Maxent models of C. versicolor s.l., the final variables comprised temperature annual range (bio7), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), precipitation of wettest month (bio13), precipitation of driest month (bio14), precipitation of warmest quarter (bio18) and precipitation of coldest quarter (bio19). For the four subtaxa, which were used as sensitivity analysis, a subset of the 19 worldclim variables was used (see results for details).

Species distribution modelling
For SDM calculation and model fitting, we combined two SDM approaches, namely the hypervolume package (Blonder & Harris 2019) for R ver. 4.0.5 (R Core Team 2021) and Maxent ver. 3.4.4 (Phillips et al. 2006). While the hypervolume approach gives a broad visualization of the realized niche based on nonhierarchically ranked variables, Maxent is an approach providing predictions based on hierarchically ranked environmental variables. Using the hypervolume package and Maxent, we computed each two final models for C. versicolor s.l., one based on records from the native range and one based on records from the native and invaded range. Furthermore, we ran a sensitivity analysis to consider the recent split of C. versicolor s.l. into four distinct subtaxa (Gowande et al. 2021). For each subtaxon, we computed Maxent as well as hypervolume models based on records from the native range only. Unfortunately, most records from the invaded range could not be assigned to one of the four species and therefore we could not include the invaded range as well. Additionally, the R-packages raster (Hijmans 2016), dismo (Hijmans et al. 2017) and ENMeval (Muscarella et al. 2014) were used for data processing and SDM calculation.

Realized niche expansion -Hypervolumes
To analyse the expansion in the realized niche of C. versicolor s.l., n-dimensional hypervolumes were calculated based on the same two datasets used for the Maxent models, one based on records from the native range (SDM Hyper_nat ) and the other one from the native and invaded range (SDM Hyper_nat_inv ). As mentioned in the introduction, C. versicolor s.l. could spread across several non-native areas and we hypothesise that an expansion in the geographic space also expands the species' realized niche in the environmental space. As we only have a small record number (14 records) from the invaded range and the time of introduction in some cases is uncertain, we could not split the datasets into smaller subsets. For hypervolume construction, we used the hypervolume_svm-algorithm available in the hypervolume package for R building a one-class support vector machine (Blonder & Harris 2019). This one-class support vector machine (SVM) transforms the input data into n-dimensional non-linear space in which the data points can be optimally separated from background by a single hyperplane. The hyperplane is then transformed backwards into the original space, which delineates an adaptive grid of random points next to the original data points. Finally, the SVM predicts if each of these points is in or out of the hypervolume (Blonder & Harris 2019). The SVM is implemented with a radial basis function (RBF), from which only two parameters, svm.nu and svm.gamma, can be defined by the user. The svm.nu-parameter defines the lower bound of support vectors and an upper bound on the fraction of training errors. Lower values reveal in a tighter wrapping of the shape to the data. The gamma-parameter tunes how far the influence of a single training point reaches and it can be seen as the inverse radius of influence of points chosen by the model as support vectors (Blonder & Harris 2019). For model construction, we used the default parameters of svm.nu = 0.01 and svm.gamma = 0.5.
Before hypervolumes were constructed, the predicting environmental variables were masked using an overall climatic envelope model (bioclim)model (Hijmans et al. 2017) to reduce calculation time as climatically unsuitable areas were excluded. The maximum possible hypervolume-model is thus nested in the bioclim-model (Nania et al. 2020). For hypervolume construction, we conducted a principal component analysis (PCA) to correctly delineate the shape of the hypervolumes to principal components (PCs) with an eigenvalue >= 1. This is necessary, as the hypervolume_svm-algorithm requires an orthogonal space for the correct delineation of the hypervolume shape. For both final hypervolumes and the sensitivity analysis, the same bioclim-model and PCA was conducted from the full record set based on occurrences from the native and invaded range together as this allows a better comparison of the models. Furthermore, it was necessary to use the full set of bioclimatic variables (excluding the derived ones as described above) as the variable selection revealed different final variables for each taxon, which would not allow comparison based on hypervolumes in environmental space. For hypervolume construction, the full record set was subset to perform one hypervolume for the native record set (SDM Hyper_nat ) and one for the native and invasive record set (SDM Hyper_nat_inv ). For the sensitivity analysis, we used the respective records of the four subtaxa. Hypervolumes were evaluated by TSS (Allouche et al. 2006).
The hypervolume_variable_importance() function was applied over each 100 replicates for all hypervolumes to compute the PCs contribution to the hypervolumes volume, respectively. Pairwise overlap statistics applying Jaccard and Sorensen index to the plotted hypervolumes for C. versicolor s.l. and the four subtaxa were calculated with the respective function of hypervolume package (Blonder & Harris 2019) to quantify the realized niche expansion in environmental space. Further, all hypervolumes were detected to contain holes. For this, a convex shape of the respective hypervolume was assumed and this idealistic shape was compared with the original hypervolume [functions: expectation_convex() and hypervolume_holes()], As the detection of holes for higher dimensions is very time-consuming, we set the number of maximum points to 10,000. We averaged the results of the hypervolume_holes() function across 100 replicates as the function uses a certain number of random points to detect holes. Afterwards, we used the hypervolume_segment() function to cluster the closest occurrence records and hypervolume_prune() function to prune small holes in the hypervolume, which may occur from too low sampling effort (Blonder 2016). Unfortunately, this function may also prune holes that occur naturally in ecological niche space. Finally, the hypervolumes were projected into geographic space using the hypervolume_project() function (Blonder & Harris 2019).
For better visualization and comparison between the native and invasive climate space, we plotted the density distributions of the environmental predictors, which were used for PCA, and the revealed PCs using the sm.density.compare() function of the sm package for R (Bowman & Azzalini 2018). Furthermore, we used the area() function of the raster package (Hijmans 2016) to calculate the predicted areas for all hypervolumes. Due to the inaccuracy of the function the results were rounded to the nearest 100 km 2 .

Invasive risk assessment -Maxent models
Maxent is one of the most frequently used correlative SDM software (Srivastava et al. 2019). The algorithm uses the principle of maximum entropy to predict a species' potential distribution from presence-pseudoabsence data and environmental variables. Presence data are confirmed occurrence records of the target species, while the pseudoabsence data are sampled in the species' background area, where the environmental conditions might allow the target species' presence (Phillips et al. 2006). As environmental background, we selected a 500 km radius around each occurrence record to provide a seamless cover across the species' distributional range. The selected background area resembles quite well the IUCN map for this taxon (Wogan et al. 2021).
Using Maxent's raw output, we calculated corrected Akaike Information Criterion [AICc; (Warren & Seifert 2011)] for each 25 replicates of each combination of model settings (Figure 1). The best combination of settings was selected according to the lowest average AICc over 25 replicates and an AUC Test above 0.7 [AUC = Area under the ROC curve (Lobo et al. 2008, Phillips & Dudík 2008, Elith & Graham 2009]. AUC was calculated by using a bootstrap approach with an 80% data split for model training and 20% used for model testing. The best fitting model was replicated 100 times using again a bootstrap approach with an 80:20 split for model training and testing to calculate the AUC values. Further, True Skills Statistic (TSS) was calculated as second evaluation metric (Somodi et al. 2017, Allouche et al. 2006). The final model was averaged across its 100 replicates. MESS (Multivariate Environmental Similarity Surfaces) maps were provided to highlight areas where the model must be extrapolated and therefore the prediction should be interpreted with caution. We visualised the SDM using Maxents cloglog output format, which provides an estimate between 0 (very unsuitable habitat) and 1 (very suitable habitat). We calculated lowest 5% thresholds for all models based on their respective average output, from which we extracted the Maxent cloglog scores for each occurrence record used for modelling. This threshold was set as threshold for presence-absence and areas below these 5% thresholds were removed [compare to Phillips et al. (2006)]. Finally, we used the area() function of the raster package (Hijmans 2016) to calculate the surface areas used as background and the predicted areas of all models. For the predictions, we calculated the surface areas after the lowest 5% thresholds were applied. Due to the inaccuracy of the function the results were rounded to the nearest 100 km 2 .
As sensitivity analysis, we calculated each a final model for C. versicolor s.str., C. irawadi s.l., C. vultuosus and C. farooqi following the same procedure described above. Subsequently, we combined the four prediction maps to a single map. For this, we calculated the median of each grid cell across the raster stack, which comprised the four prediction maps. As combined MESS area we used the areas where all four subtaxa showed MESS areas.

Realized niche expansion -Hypervolumes
The performed bioclim-model predicts and excludes larger unsuitable areas in northern North America, the Andes, the Pole regions, northern and south-eastern Africa, northern, central and western Asia and north-western and southern Australia. The conducted PCA reveals four PCs with an eigenvalue >=1 (Table 2). PC1 is mainly correlated (factor loadings >0.70 or <-0.70) with annual mean temperature (bio1), minimum temperature of coldest month (bio6), mean temperature of driest (bio9), wettest (bio8), coldest (bio11) and warmest quarter (bio10). PC2 is mainly correlated with annual precipitation (bio12), maximum temperature of warmest month (bio5), precipitation of wettest month (bio13) and quarter (bio16). PC3 is mainly correlated with precipitation of driest month (bio14) and quarter (bio17). PC4 only shows low correlations but is strongest correlated with precipitation of warmest quarter (bio18). For all models, the PCs' variable importance over 100 replicates to the models' volume reveals the greatest importance for PC1 and PC2 (Table 2). Most models show a moderate performance (>0), however, the hypervolumes of C. versicolor s.str. (27 records) and C. farooqi (15 records) perform poor due to the  (Table 3). Further, in all hypervolumes several holes are detected. The volumes, point densities and numbers of random points of the holes strongly differ among the models (Table 3).
The density distribution plots can give a more detailed explanation of a potential realized niche expansion. We consider the realized niche expansion when novel environmental conditions are met in the invaded range, which were not found in the native range. For the four PCs, the density plots reveal major expansions in the realized niche of C. versicolor s.l. for all PCs (Figure 2). Major expansions in the realized niche can also be found for temperature annual range (bio7), precipitation of driest month (bio14), warmest (bio18) and coldest quarter (bio19), while only minor expansions are found for the remaining variables ( Figure 3). Considering the four subtaxa, the density distribution plots reveal strong differences among the taxa and among C. versicolor s.l. due to the different geographic distributions and different climatic conditions (Supplementary Material S2). The pairwise overlap statistics for SDM Hyper_nat and SDM Hyper_nat_inv reveal 1.00 ± 0.00 a large unique fraction for the latter model, which also indicates an expansion in the realized climate niche due to the invaded ranges. Further, the low Jaccard and Sorensen Indices as well as the unique fractions of the hypervolumes for the four subtaxa show that the respective subtaxa only represent small subsets of the whole climatic niche of C. versicolor s.l. (Table 4). The predictions of SDM Hyper_nat and SDM Hyper_nat_inv reveal surface areas of 30,166,500 km 2 and 35,361,100 km 2 , respectively, which correspond to a gain of 17%. The hypervolume models for the four subtaxa reveal much smaller surface areas (Table 5).   according to AICc reveals moderate to high AUC values for all SDMs. We selected the model with the lowest AICc for further processing (Table 6).
All final models also perform moderate to high AUC and TSS values across their 100 replicates (Table 6). For SDM Maxent_nat , temperature annual range (bio7) contributes highest to the models performance, followed by precipitation of wettest (bio13) and driest month (bio14). For SDM Maxent_nat_inv , the variables contribute similar to the model with temperature annual range (bio7) followed by precipitation of wettest month (bio13) and mean temperature of driest quarter (bio9), whereas the contribution for the models of the four subtaxa differs from the latter two models (Table 7).
The used background areas reveal surface areas of 6,932,200 km 2 (SDM Maxent_nat ) and 6,941,200 km 2 (SDM Maxent_nat_inv ), respectively, while the predictions (5% thresholds applied) result in an area of 40,250,800 km 2 and 53,142,600 km 2 . This corresponds to gains of 481% and 666%. The sensitivity analysis (combined map of the four subtaxa; Fig 3) reveals a surface area of 14,025,100 km 2 , which is 35% of the predicted area of SDM Maxent_nat . For the four subtaxa, the predicted surface areas strongly differ among each other (Table 6).

Species complex vs. subtaxa
In contrast to the Maxent models, the hypervolumes reveal similar but rougher predictions as they provide presence-absence maps without probabilities (Figure 2, Supplementary Material S2). For the invasive risk assessment, we therefore focus on the Maxent models. The two Maxent models for C. versicolor s.l., SDM Maxent_nat and SDM Maxent_nat_inv , reveal similar global predictions and show highly suitable environmental conditions across large parts of the Oriental realm ( Figure 3). The native distribution of C. versicolor s.l. largely covers this region, where both models predict highly suitable environmental conditions. Habitat suitability decreases towards the Himalaya and towards drier regions (i.e. northeast India, Afghanistan, Pakistan) where the number of occurrence records is also lower or no records are confirmed (Figure 3). According to our models, the Pacific and Indian Ocean islands of the Oriental realm towards the Wallace line are climatically moderately to highly suitable including the already invaded regions [Borneo, Celebes (= Sulawesi), Philippines, Maldives, Singapore, Andaman  Table 5. Predicted surface areas of hypervolume and Maxent models. For Maxent models, the surface areas of the background and the predictions (5% thresholds applied) are provided. The percentage surface area in respect to the model for C. versicolor s.l. is given in brackets.
[native] Table 6. Results of the five best Maxent ensembles according to AICc. Models for C. versicolor s.l., based on records from the native range (SDM Maxent_nat ) and the native and invaded range (SDM Maxent_nat_inv ), and a sensitivity analysis based on the four subtaxa (C. versicolor s.str., C. irawadi s.l., C. vultuosus, C. farooqi): Regularization multiplier, feature classes, number of parameters, AICc, AUC Training and AUC Test . The model used for further processing is highlighted in bold font. Both Maxent models reveal that very large parts of the Afrotropics (Figure 3) are climatically suitable for C. versicolor s.l.. Especially, the regions along the equator and the coastal regions as well as the Indian Ocean islands, including all already invaded regions (Mascarene archipelago, Seychelles, Kenya, Oman), provide moderately to highly suitable conditions. Habitat suitability decreases towards the Sahel zone and the deserts of the Arabian Peninsula where only the very coastal regions show lowly to moderately suitable environmental conditions. Furthermore, large MESS areas are found across the Sahara, the Arabian Peninsula and the coastal regions of southern Africa. SDM Maxent_nat_inv also predicts less climatically suitable area than SDM Maxent_nat . For the Neotropics, our Maxent SDMs predict high climatic suitability across most parts of this biogeographic region including the Caribbean, Central America, northern and central South America (Figure 3). Southern Florida, where C. versicolor s.l. was already confirmed, also shows lowly to moderately suitable environmental conditions for both SDMs. MESS areas are found in parts of the Andes and some southern parts of mainland South America. The MESS areas are decreasing for SDM Maxent_nat_inv .

Regularization Features nParameters
Considering the Maxent models as well as the sensitivity analysis, the Palearctic and Nearctic realms are found to be almost entirely climatically unsuitable or low habitat suitability is predicted for both Maxent models (Figure 3). For the Palearctic realm, the only moderately suitable areas are found along some coastal regions and islands of the Mediterranean Sea and the Atlantic coast of the Iberian Peninsula. Furthermore, in Asia the transition zone to the Paleotropics is also partially predicted as lowly to moderately suitable. Large MESS areas are shown in North Africa, the Middle East and North to East Asia. For the Nearctic realm, only the transition zone to the Neotropics shows some regions with lowly to moderately suitable environmental conditions. Large MESS areas are found in the North towards the North Pole.
Compared to the models for C. versicolor s.l., the averaged map of the sensitivity analysis ( Figure 3) reveals less suitable area in the Oriental realm (i.e. southern India, Sri Lanka, large parts of mainland Southeast Asia and Indonesia) and the Australian realm (i.e. coastal region of northern, eastern and southern Australia), Large MESS areas are found on Borneo, New Guinea and several adjacent islands. Large parts of mainland Africa across the equator as well as Madagascar are climatically suitable for most of the four subtaxa. Further, several coastal areas i.e. of southern Africa and the Arabian Peninsula are suitable as well. Large MESS areas are found across the Sahara and the deserts of the Arabian Peninsula. Large parts of Central America, the Caribbean, northern and central South America are climatically suitable for most of the four subtaxa. Furthermore, small parts of southern South America are suitable as well. Large fragmented MESS areas are found across the Pacific coast of South America. The predictions for each subtaxon used for the sensitivity analysis differ strongly among each other (Supplementary Material S2). Particularly, the models for the species with low sample sizes reveal large MESS areas and it must be mentioned that more accurate predictions require a more species-specific tuning and a higher sample size.

Discussion
The geographic expansion into novel climatic conditions reveals an expansion of the realized climatic niche for C. versicolor s.l. which is supported by the hypervolumes and density plots. The results of the invasive risk assessment predict large areas of suitable climatic conditions for most parts of the tropics and subtropics for the species complex and the four subtaxa. Furthermore, the models for the species complex reveal larger surface area predictions and also partially more realistic predictions as the sensitivity analysis.

Introduction pathways
The unintentional introduction as stowaway from ornamental plants, timber or other nature products, seems to be the most frequent introduction pathway into non-native oversea regions for this lizard (see introduction for details). While a combination of deforestation and manmade transport might also boost the species' expansion in its native range. A comparison with maritime bottlenecks of container traffic revealed that several core and secondary trading routes are intersecting with the species' range (https:// porteconomicsmanagement.org/pemp/contents/ part1/interoceanic-passages/connectivity-patternworlds-major-maritime-bottlenecks/) and might be likely introduction pathways into the Afrotropical, Australian and further regions of the Oriental realms, where our models predict suitable climatic conditions for large areas. Our invasive risk assessment reveals that coastal regions often provide more climatically suitable conditions, also if main parts of the country are predicted as unsuitable (i.e. Oman). Considering that C. versicolor s.l. prefers open and often anthropogenic landscapes, coastal ports with high amount of traffic could serve as stepping stones and might allow to disperse easily into novel and even far distant regions. The Strait of Malacca is one of the most important trading routes worldwide and the main passage between the Indian and Pacific Oceans. Specifically, the port of Singapore, which is among the most important trans-shipment hubs, might be an important stepping stone for C. versicolor s.l. to reach novel regions. For Kenya it is suggested that the species derives from the introduced populations of Réunion and Mauritius (Sandera 2009). Also for other reptiles such as Lepidodactylus lugubris it was discussed that its invasion success is linked to maritime bottlenecks (Nania et al. 2020).
Our results reveal that the already invaded regions provide suitable climate for all four subtaxa and based on this, we cannot assign a certain species to an invaded locality. However, due to geographic proximity, it can be assumed that the non-native populations on Borneo, Sulawesi, the Philippines and the Andaman Islands most likely belong to C. irawadi s.l., which should also be the case for the populations in Malaysia and Singapore, if they truly derive from natural invasions. According to Enge & Krysko (2004) and Krysko et al. (2011), the population in Florida (USA) derives from animals from Pakistan and if so, they would belong to C. farooqi s.l. Unfortunately, genetic data on the non-native populations are scarce and only the populations on the Maldives were assigned to C. versicolor s.str. (Gowande et al. 2021).
Next to stowaway, the pet trade is one of the main pathways of introduction for non-native reptiles and amphibians (Kraus 2009). However, only the population in Florida (USA) seems to derive from escaped animals from the pet trade (Enge & Krysko 2004, Krysko et al. 2011). According to own observations (PG), C. versicolor s.l. is rarely kept long-terms in captivity as the species is very sensitive and often subject to poor husbandry conditions. This might be the reason that no other population derives from this source and the risk of invasion due to the pet trade seems low.
Despite C. versicolor s.l. has a huge native range, surprisingly only few occurrences are known from outside this area. Reasons for this might be that the taxon is just neglected by biologists and labelled as "unspectacular lizard". Urban environments, where the species is found, might often be understudied as they seem uninteresting for most herpetologists. The case of the African Clawed Frog (Xenopus laevis) in Portugal shows that exotic species can be undetected for several decades, even in densely urbanised landscapes (Sousa et al. 2018). Alternatively, the species might arrive in novel and climatically suitable environments much more often than it is observed but it cannot establish populations, how it was observed on Mahé in 1985-86 (Matyot 2004). Reasons for this might be complex and often include biotic interactions (i.e. competition with other species) but also stochastic effects (i.e. low number of founder individuals; Kraus 2009).

Realized niche expansion
Our hypervolumes, the pairwise overlap statistics and the density distributions reveal an expansion of the realized niche in geographic as well as in environmental space during the invasion process. For the novel invaded regions, C. versicolor s.l. meets a lower temperature annual range and less precipitation in the driest month, warmest and coldest quarter compared to the native range ( Figure 3). For all four PCs, the density distribution plots reveal a major expansion of the realized niche as well.
A closer view on the Maxent sample predictions (Supplementary Material S1 and S2) shows that the records from the invaded range in average reveal higher suitability (mean ± SD [range]: 0.55 ± 0.22 [0.19 -1.00]) compared to the native records (0.71 ± 0.18 [0.33 -0.96]), which also indicates that the invaded localities are not on the edge of the species' realized climatic niche but in contrast provide highly suitable climatic conditions. However, particularly for the Maxent models, some predictions of climatically low suitability regions might be an artefact derived from the 500 km buffer around each record used for modelling. This buffer was necessary to obtain seamless cover across the native range of C. versicolor s.l. and to provide Maxent contrasting environmental data (absencepseudoabsence points used for modelling). Further, it allows a better comparison among the different subtaxa. However, some occurrence points close to high mountain ranges (i.e. Himalaya) led to the inclusion of environmental gradients that are outside of the species' ecological niche space and might led to erroneous predictions (i.e. temperate rainforests in Canada).
In both realized niches (SDM Hyper_nat and SDM Hyper_nat_ inv ), several holes are detected. Holes in hypervolumes, or modelled ecological niches, respectively, might be explained by the hypothesis that niche evolution does not always proceed through simple shifts. Genetic or developmental constraints or natural selection might prevent the presence of certain phenotypes, leading to the presence of holes in trait or potential niche space (Blonder 2016). Furthermore, holes might be the result of a genetically diverse fitness landscape. Populations of very different phenotypes can be connected by a small number of intermediate genotypes (Blonder 2016). This hypothesis is supported by the fact that C. versicolor s.l. and even its subtaxa are a complex of several cryptic species (see introduction for details). Further, the most holes might most likely be an artefact from undersampled niche space as C. versicolor s.l. is a very abundant and widespread taxon. Unfortunately, data on real absences of the species were not available. Furthermore, the holes might also derive from too restrictive settings used for hypervolume construction.

Conclusion
In summary, much of the tropical and subtropical regions seem climatically suitable for the C. versicolor species complex. Considering this, as well as the species' dispersal history and evidence for negative impacts on the native biodiversity, this lizard should be further monitored. According to our basal risk assessment (SDM Maxent_nat_inv ), the species is currently inhabiting 13% of its potential range but could find suitable climatic conditions on a global surface area between 14,025,100 km 2 (sensitivity analysis) and 53,142,600 km 2 .
Regions with a high degree of endemism (i.e. islands like Mascarenes, Madagascar) might be particularly sensitive to an invasion by this lizard and the subsequent negative impacts (i.e. Duttaphrynus melanostictus on Madagascar; Moore et al. 2015). Furthermore, the introduction of C. versicolor s.l. to regions of other genetic lineages of this species complex might most likely lead to hybridizations and might threaten close relatives that are endemic to small areas, which is the case for many Calotes species. Therefore, strategies to hamper and prevent further spread of C. versicolor s.l. are urgently needed. Strategies to prevent stowaways, especially in regions with a high amount of human traffic as ports or airports, are also required. Lastly, a further systematic revision of this species complex will result in different predictions for each species and will be subsets of our current prediction.

Author Contributions
PG did the main work of model construction and the main part of writing. NWCT and DR contributed in model construction and gave advice. All authors read the MS several times and contributed to the final version.

Data Accessibility
The data that support the findings of this article are available on request from the corresponding author (PG).

Supplemental Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Table S1. Maxent sample predictions and coordinates for SDM Maxent_nat , SDM Maxent_nat_inv and the models for C. versicolor s.str., C. irawadi, C. vultuosus and C. farooqi as well as the mean values ± SD [min -max] across the sample predictions for each model Appendix S2. Hypervolume and Maxent models as well as the relevant density distribution plots for the four subtaxa C. versicolor s.str., C. irawadi, C. vultuosus and C. farooqi and plot of the sample predictions for SDM Maxent_nat_inv