Predicting The Invasion Risk of Non-Native Reptiles as Pets in The Middle East

Wildlife trade for non-native pets is an important and increasing driver of biodiversity loss and often compromises the standards required for protection. However, the growing interest in non-native pets has posed the issue of invasive non-native species to wildlife managers and conservationists. Instituting effective policies regarding non-native species requires a thorough understanding of the potential range of the species in new environments. In this study, we used an ensemble of ten species distribution models to predict the potential distribution for 23 commonly traded species of reptiles across the Middle East. We used ten modeling techniques implemented in the Biomod2 package and ensemble forecasts. Final models contained fourteen environmental variables, including climatic, topographic, and land cover/land use variables. Our results indicate that all Middle Eastern countries included suitable habitats for at least six species, except Qatar, Kuwait and Bahrain for which the models did not predict any suitable habitats. Our study showed that Lebanon, Palestine, Turkey, and Israel face the highest risk of biological invasion by the species on the whole. Also, the results showed that Centrochelys sulcata, Chamaeleo calyptratus, and Trachemys scripta posed the highest risk of spreading in this area. Information on which species pose a greater danger as invaders and the possible impacts of their introduction will be a valuable contribution to the development of conservation plans and policies.


Introduction
Over the past two decades, non-native pet trade has become the main source of non-native reptile and amphibian species worldwide (Capinha et al., 2017;Stringham & Lockwood, 2018). The aesthetic and entertainment value of some species has made them coveted non-native pets (Reed, 2005) and created a market around capturing, breeding and selling them (Bush et al., 2014) while docility in captivity, large size and rarity has contributed to the popularity of other species as pets (e.g. pythons and boas) (Luiselli et al., 2012;Reed, 2005). As a result, over half a million live reptiles are traded annually across the globe (Karesh et al., 2005). "Non-native pets" refers to species without a history of domestication by humans which are sold as pets; these species are caught from their natural habitats or grown in human-made environments. These animals are not sold with the goal of being released into the wild. However, non-native species are frequently observed as free-living, indicating that large numbers of these species eventually nd their way into natural ecosystems (Hulme, 2015). The released individuals might establish populations in the wild, which could in turn lead to negative consequences for native species and the environment (Stringham & Lockwood, 2018).
Invasive reptiles can disturb ecosystems through predation, herbivory, competition, genetic hybridization (Kraus, 2015), or introduction of non-native pathogens (Burridge et al., 2000;Nowak, 2010). The introduction of non-native species can lead to disturbance in local ecosystems and cause decline or change in local populations (Holbrook & Chesnes, 2011); extinction of native species as a consequence of the introduction of non-native species has also been reported (Savidge, 1987). The detrimental consequences of invasive reptiles are not limited to ecosystems and can extend into economic and social realms as well (Bomford et al., 2005). Wildlife managers should also take note that many of the illegally traded reptile species are venomous Knowledge about the invasion range of these species plays a vital role in understanding the ecology of invasive species and in creating a sustainable conservation and management plan. Predictive models, such as species distribution models (SDMs), have become an essential tool for researchers to predict and assess the distribution of non-native species by quantifying species-environment relationships (Guisan and Thuiller, 2005;Elith and Leathwick, 2009). SDM rst became species in Australia (Guisan and Thuiller, 2005). In the 1990s the approach experienced a renaissance caused by simultaneous developments in computer-and statistical sciences (Elith and Leathwick, 2009 Jha and Jha, 2021). There are now many methods used for distribution modelling, varying in how they approach model selection, how they de ne tted functions and interactions, whether they can handle imperfect detection and sampling biases and so on (Franklin, 2010;Guisan et al., 2017). Predictive outcomes across SDM methods are known to be variable, and the choice of modelling method can signi cantly affect model predictive performance. However, no one method is consistently superior in performance across species, regions and applications (Segurado & Araujo, 2004). This makes it di cult to choose which method (individual model hereafter) to use, prompting the idea of combining predictive outputs from different models in a so-called ensemble (Araújo & New, 2007). The underlying philosophy of ensemble modelling is that each individual model carries both some true "signal" about the relationships the model is aiming to capture, and some "noise" created by errors and uncertainties in the data and the model structure. Ensembles combine models with the intention of obtaining better separation of the signal from noise (Araújo & New, 2007;Dormann et al., 2018). These types of ensembles have been shown in many cases to have superior predictive performance to individual models (Seni & Elder, 2010).
In this study, an ensemble framework of ten SDMs was used to investigate the potential distribution of the 23 most commonly traded non-native reptile species in the Middle East. Thus, the aims of this study were as follows: (1) to identify potential distributions of the 23 most commonly traded non-native reptile species in the Middle East; (2) to determine the important environmental variables that shape the variations in the potential geographic range of nonnative species; and (3) identify countries in the Middle East that are at risk of biological invasion by these 23 species.

Studied species and presence points
In this study, 23 species which met two criteria were selected: 1) the species had to be widely traded globally and 2) an adequate number of presence points had to be accessible for the species (more than 50 presence points). We initially found 89 species that were traded on social media, uno cial platforms, and on the black market, but only 23 species of the initial set had enough presence points to be used for modeling. We compiled presence points for the 23 species from the Global Biodiversity Information Facility (GBIF), VertNet, iNaturalist, and Berkeley Ecoinformatics Engine. GBIF is an international platform aiming to provide publicly-available data on life on Earth. GBIF focuses on providing templates, standards, and tools for sharing information between institutions and individuals. Data available on GBIF come from a number of sources including eld observations by experts and citizen scientists, museum records, and surveys. VertNet is an NSF-funded project created through the combination of four taxon-speci c repositories (FishNet, MaNIS, HerpNET, ORNIS) to make biodiversity data available online. iNaturalist, jointly funded by California Academy of Sciences and National Geographic is a crowdsourced species identi cation system and a platform to record presence points. Berkeley Ecoinformatics Engine allows access to several biodiversity information repositories (https://holos.berkeley.edu/LearnMore/datasources/) through its R package. Each database was accessed through their respective packages in R ("rgbif", "rvertnet", "rinat", and "ecoengine" respectively). Queries were made using the scienti c name of the species, and georeferenced presence points from January 1, 1998 to January 1, 2018 were extracted. In this study, we used data from the last 20 years to minimize the error in the data (Alhajeri and Fourcade, 2019).
Here, we used the datasets for both native and non-native geographical ranges of the specie. In order to obtain more accurate results in modelling, presence points were compared with the native range of the species. In the next step, only the none-native range in countries where established populations of the species had been reported was included (Cordier et al., 2020). In cases where the same data point was recorded more than once, duplicates were omitted in R.
After removing duplicates, presence points were used to model the potential distribution of the species. Moran's correlograms were created in Spatial Analysis in Macroecology (SAM; Rangel et al., 2006) to determine spatial autocorrelation for the species. The signi cance of Moran's I was tested using a randomization test with 9,999 Monte Carlo permutations, adjusted for multiple testing. In instances where spatial autocorrelation was detected, we restricted the testing and training data as follows: rst, a distance threshold was set using distance lags which showed positive spatial autocorrelation (from 10 to 25 km); in the next step, the distance between data pairs was compared with the threshold, and pairs which had a distance smaller than the threshold were placed in the same partition (Parolo et al., 2008). Given that we only had access to species presence points, a number of pseudo-absences were randomly generated (Elith et al., 2011). To evaluate human impact on the reptile, we calculated the human footprint index using data on settlements, land transformation, accessibility and infrastructure.
To exclude highly-correlated inputs, all pairs of variables were checked for correlation using Pearson's correlation coe cient (r) in SDMToolbox (Brown, 2014) for ArcGIS. r > 0.70 and r < -0.75 were considered as the threshold values, and variables with stronger correlation were excluded (Kalboussi & Achour, 2018), resulting in a nal set of 14 variables (including: 9 climatic, 2 topographic, and 3 land cover/land use variables) ( Table 1). Only algorithms that met the 0.7 quality threshold were included in the ensemble calculations; those that do not meet this standard were discarded (Thuiller et al., 2009(Thuiller et al., , 2014. To compensate for variation between the number of pseudo-absences (PA), replicates, and model runs based on which algorithm is used, we grouped the 10 algorithms into three PA subsets to allow for optimal ensemble model performance Groups 1, 2, and 3 (Table 3 . Group results were compared for accuracy using the true skill statistic (TSS; scaled -1 to +1, where performance ≤ 0 means the model output is no better than random and > 0 means the proposed model successfully distinguishes between suitable and unsuitable habitat). Models were also evaluated for sensitivity and speci city. Sensitivity is the accuracy rate for true positive outcomes (i.e., the probability that the model correctly predicted presence). Speci city is the accuracy rate for true negative outcomes (i.e., the probability that the model correctly predicted absence) (Allouche et al., 2006). Finally, we used a quality threshold to accept models and TSS, sensitivity and speci city values to compare and select among alternative ensemble SDMs (Table 4).

Modelling techniques and ensemble forecasting
Variable importance is determined in Biomod2 package using permutation as follows: after models are calibrated, a baseline prediction is generated. In the next step, the values of one variable are randomized and a new prediction is generated using the randomized values of the randomized variable and the unchanged values of other variables. The correlation between the baseline prediction and the new prediction (using one randomized variable) indicates the relative importance of the variable that has been randomized (Thuiller et al., 2009). This procedure does not depend on the modelling algorithms employed.

Results
Evaluation of modeling results based on the TSS, sensitivity and speci city values showed that the ensemble SDMs performed better than individual models (Table 4).
Results of variable importance tests showed that the importance of environmental variables varied for different species. The three most important environment variables that determine species potential geographic distributions are shown in Table 2. The most important environmental variables in predicting potential geographic distribution for most of the species were temperature seasonality, land cover, and NDVI.
The potential distribution maps of the ensemble model for each species are presented in Fig 1. The current overlap between suitable habitats and the Middle East is shown in Considering the area of suitable habitats for the studied species, this study showed that Lebanon, Palestine, Turkey, and Israel face the highest risk of biological invasion by the species on the whole. We ranked countries in the Middle East in the order of highest to lowest risk of biological invasion (Fig 2a). We also ranked the studied species based on their risk of biological invasion; accordingly, African Spurred Tortoise, Veiled Chameleon, and Yellow Bellied Slider had the highest risk of spreading in this area (Fig 2b).

Discussion And Conclusion
The Our results indicate that all Middle Eastern countries included suitable habitats for at least six species, except Qatar, Kuwait and Bahrain for which the models did not predict any suitable habitats. Our study showed that Lebanon, Palestine, Turkey, and Israel face the highest risk of biological invasion by the species on the whole (Fig. 2a). Previous studies in these countries have not paid much attention to the invasion and trade of non-native reptiles, which makes the management and conservation of these species di cult or even impossible. Also, the results showed that African Spurred Tortoise, Veiled Chameleon, and Yellow Bellied Slider posed the highest risk of spreading in this area (Fig. 2b). Although these species have limited global range, their ability to spread in other parts of the word should be taken into consideration.
The Veiled Chameleon is natively found in Saudi Arabia and Yemen, mostly in plateaus and grasslands (Showler, 1995;Schmidt, 2001 (Newbery, 1984). Due to its extensive global spread, the species has been categorized as one of the 100 worst invasive species (ISSG, 2021). The African Spurred Tortoise ranks second among terrestrial turtles in terms of size, and naturally occurs in the western parts of the Sahel region. Although it is a threatened species and its population has been in decline through its range, little information is available about the causes of this decline. Researchers have pointed to competition with cattle, wild res, and pet trade as the possible drivers of the decline in its population (Petrozzi et al., 2018). Despite the species having been reported from other countries, it is not clear whether it has been able to establish populations outside its native range.
Wildlife managers should pay attention to species potential of invasion when making decisions about the import of non-native species into their countries.
Although the model did not indicate any areas with a potential distribution for Amazon Tree Boa and Burmese Python, managers should exercise caution when making decisions about these species as well, and refer to previous studies or conduct more detailed research to further assess the risks associated with non-native species.
SDMs take into account only the bioclimatic realized niche of species, meaning that the model assumes the species is only able to establish in areas with a similar climate to its native range (Kearney 2006). It is obvious that biotic and abiotic factors and their interactions can alter a species' realized niche (Broennimann et al., 2007). This leads to the generated niche predictions in the new environments being only rough approximations. Consequently, a mismatch between bioclimatic conditions in the native and introduced range can lead models to underestimate the extent of suitable areas. Further studies focusing on each species and utilizing a wider range of data with better accuracy, including biotic variables, abiotic variables, and their interactions are needed to make more nuanced judgments about each species. Such studies would also bene t from a comparison with areas that have the studied species as invasive species.
In the meantime, managers should not interpret the absence of suitable habitats for some of the studies species as a license for unrestricted import or introduction. Thus, the ndings of this study are mainly helpful as a basis for imposing restrictions on species for which a suitable habitat was found; our ndings should not be seen as justi cation for allowing the introduction of species or loosening regulations. However, SDMs can still serve as the rst step in understanding the biological invasion history of species (Silva Rocha et al., 2015).   CITES Appendices: I (it lists species that are the most endangered among CITES-listed animals and plants), and II (it lists species that are not necessarily now threatened with extinction but that may become so unless trade is closely controlled).  (Lek and Guegan, 1999) Surface range envelope (SRE) (Busby, 1991) Generalized additive models (GAM) (Guisan et al., 2002) Generalized linear models (GLM) (Guisan et al., 2002) Maximum      Potential distribution map of the studied species in the world (red: suitable habitats, gray: unsuitable habitats)