Reaching new environments through illegal trade: evidence of a widely traded turtle in Colombia

A major threat to biodiversity is illegal trade, with many unwanted wildlife pets released into exotic environments outside their native distribution. Therefore, many potential invasive species have established in new ecosystems. Ecological niche modeling (ENM) has been used to predict and compare the environmental conditions of natural and exotic population in many groups. We used ENM to compare the climatic niche between natural and exotic areas of Trachemys venusta callirostris, one of the most traded turtles in Colombia. We generated a niche model using the MaxEnt algorithm through the R package kuenm to test several parametrizations and four sets of fresh water environmental predictors. Models were calibrated in the native distribution and projected to non-native zones in Colombia to identify suitable areas for the species. Further, we use a niche similarity test to compare native and exotic environmental space. We found few suitable areas within the projected zone even when using extrapolation; there was a greater suitability in the Magdalena River basin than in the Cauca River basin. Low similarity was detected between the niche comparison of native and exotic areas, suggesting that exotic populations have reached different environmental conditions than the native zone through ilegal trade. Although there was low extrapolation in the exotic area, the models projected ideal conditions in localities with new records for this turtle. The generalist strategies for feeding, thermoregulation, and reproduction in changing conditions may help this T. v. callirostris establish in new ecosystems, and with no current knowledge on dynamics between this exotic species and local fauna, its effects on aquatic communities are unpredictable.


Introduction
Illegal wildlife trade is a major threat to biodiversity (Zhang et al. 2020;Morton et al. 2021) and is a profitable activity for some countries, especially when the final goal is to use the animals as pets (Yi-Ming et al. 2000;Zimmerman 2003). Although national and international regulations have been established to protect biodiversity, measures are insufficient to control trading (Lee et al. 2005;Spangenberg 2007). In many cases, wildlife buyers cannot keep the acquired animals and end up releasing them in places outside their native range. As a result, many potential exotic and invasive species have successfully established in new locations around the world (Johnson and Padilla 1996;Sakai et al. 2001;Rödder et al. 2009;Banha et al. 2017). The presence of a new species in an ecosystem influences the direct and indirect competition for resources with native species (Sakai et al. 2001;Brown et al. 2002;Braks et al. 2004). Therefore, detecting potential areas out of a species´ known native distribution is a main objective for ecosystems conservation (Sakai et al. 2001;Allendorf and Lundquist 2003;Mehta et al. 2007).
Ecological niche modeling (ENM) has been widely used to estimate the potential ranges of invasive or exotic species in new localities (Rödder et al. 2009;Rodrigues et al. 2016;Khosravi et al. 2019;Liu et al. 2019). Further, niche analysis allows us to understand potential species response to new environments by comparing native and invasive populations (Rödder et al. 2009;Broennimann et al. 2012;Fernández and Hamilton 2015;Di Cola et al. 2017). These tools provide valuable information to promote ecological and social strategies to mitigate or control the invasion of alien species (Peterson and Vieglais 2001;Clout and Williams 2009;Rodrigues et al. 2016).
It is suspected that new locations of this slider turtle in Colombia are due to the release of unwanted pets in natural areas (Adames-Jiménez et al. 2018). Further, T. v. callirostris is the most traded turtle in the country, accounting for more than 50% of the seizures by Regional Environmental Corporations RECs (Bonilla et al. 2012;Arroyave et al. 2014;Bock et al. 2015;Suárez Giorgi 2017). Therefore, this turtle occurs in new suitable regions where its population may increase competition for resources with other species (Sakai et al. 2001;Ficetola et al. 2009;Díaz-Paniagua et al. 2011). This species could shift from exotic to a potential invasive species in some regions as has occurred with other traded slider turtles (Ficetola et al. 2009;Rödder et al. 2009;Banha et al. 2017). In this study we tested if the species illegal trade allows to reach new environment al conditions in exotic areas using one of the most traded turtles in Colombia, Trachemys venusta callirostris, as a study model. We expected that the climatic conditions of introduced populations of T. v callirostris would differ statistically from the native niche using an ENM approach and we identified potential suitable regions for establishment of the turtle.

Occurrence and environmental data
Occurrences of T. v. callirostris were obtained through the Global Biodiversity Information Facility database (GBIF 2021) and personal field observations. We split the records of the species into native and exotic areas based on previous systematic, taxonomic, and distribution studies (see Fritz et al. 2012;Vargas-Ramírez et al. 2017;Páez et al. 2022). Then, we cleaned the data based on Cobos et al. (2018), removing duplicated occurrences, then applying a spatial thin of one kilometer to have only one record in each raster pixel, and deleted potential identification results based on previous taxonomic revisions (e.g., Páez et al. 2022). Finally, we move those occurrences out the pixels calculating the Euclidean distance of each record to the surrounding pixels, then moving manually to the closest pixel, using Arcgis (ESRI 2022).
Because T. v. callirostris is restricted to freshwater habitats we obtained environmental variables from Domisch database (Domisch et al. 2015) with a spatial resolution of 30 s. These variables are more related to the ecology of freshwater species and present relevant implications within model construction and niche comparison in fresh water ecosystems (see Nori and Rojas-Soto 2019). For selection of climatic variables, we first created the accessibility area (M; Soberón and Peterson 2005) for the native distribution and a projection zone that represents the introduced populations. Both areas were estimated using level 6 HydroShed basins (Lehner et al. 2021). Then, we performed a Pearson correlation test using the variables from the native distribution area and removed those with correlation coefficients below −0.8 and above 0.8 (Echeverry-Cárdenas et al. 2021). However, because in some cases multiple highly correlated variables can be important for the species, we created four environmental sets to test during the niche modeling process ( Table 2 in Supplementary material).

Ecological niche modeling
We used native occurrences (80% for training and 20% for model evaluation of each one) to model the ecological niche using the Maxent algorithm (Elith et al. 2011) in the kuenm R package (Cobos et al. 2019) and projected each model in geography for the native and exotic area using all extrapolation types to compare the results. The parametrization of the model did not include "hinge" and "threshold" features because these can lead to a complex function response (Elith et al. 2011). Further, we tested different levels of complexity for model construction as Warren and Seifert (2011) suggested, shifting the regularization multiplier values from 0.1 to 1 at 0.1 intervals and then 2, 3, 4, and 5. We selected the best model for both native and exotic based on the evaluation metrics performed in Kuenm (AUC ratio, Omission rate, and AIC) and used the ten-percentile threshold to transform the cloglog Maxent's output into presence/absence maps. If several models fulfilled the evaluation criteria, we estimated a final model from each climatic set using the median of the replicates obtained for the native and exotic zones for all extrapolation methods. Finally, we applied the Mobility-Oriented Parity test (MOP; Owens et al. 2013) for the models with clamping extrapolation, which represent areas of potential extrapolation risk.

Niche similarity test
To test if T. v. callirostris could be reaching new environments in the exotic area, we performed a niche similarity test using the ecospat R package (Broennimann et al. 2012;Di cola et al. 2017). This test summarizes the environmental predictors in two principal components and perform a kernel density function for native and exotic occurrences. Then estimates the overlap or similarity between both (native and exotic population) when the native niche is randomly introduced in the exotic area, then a high environmental similarity value trend to one while low similarity to zero (Broennimann et al. 2012). If the observed value (obtained by comparing climates from populations from both zones) falls within the null models, then the differences are just explained by chance, because environments of the exotic zone are too different from the native one (Di cola et al. 2017). We performed niche similarity test between both areas for each environmental data set selected by kuenm evaluation metrics to obtain the best models and their corresponding potential distribution.

Results
We obtained 60 occurrences of the species for the native distribution and 32 for the exotic localities ( Table 1 in Supplementary material), of which four presented evidence of reproduction (Fig. 4). We obtained four ecological niche models with different parametrization that fulfilled the evaluation criteria from kuenm (Table 3 in Supplementary material). From these four models, two used the set three and two the set four of environmental predictors. Both consensus models (using the two environmental dataset predictors) showed slight differences in the number of predicted pixels (as presence) and set four had more pixels in native and exotic zones even when using different extrapolation methods (Fig. 4). Predictions of potential distribution in the native area showed the expected presence in the geographical basins of Caribe and Magdalena where most of the occurrences have been recorded (Bock et al. 2010(Bock et al. , 2012Fig. 2).
The projection of the model in the exotic area showed few potential presence zones even using free extrapolation and clamping. We found a similar Cesar-Magdalena, c Sucre, and d Eastern Antioquia populations. Black points represent species occurrences in its natural area used for niche modeling number of predicted pixels in the exotic zone using both predictors' data sets when no extrapolation was performed. However, by applying clamping or free extrapolation, we found that set 4 predicted more pixels than set 3 (Fig. 3). Our models predicted suitable areas in the nearest zones of the Magdalena Basin, near urban centers in Caldas and Tolima. However, we identified few suitable areas in the Cauca basin where two new localities have recently been found with evidence of reproduction. Further, important locations such as Bogotá (Cundinamarca) and Cali (Valle del Cauca) where adults have been detected, but no juveniles, presented few potential zones (Fig. 4).
MOP results showed a greater similarity (~1) among areas of occurrence in localities of Caldas, Cundinamarca, Quindío, and Valle del Cauca, but a low similarity of extrapolation areas occurred in localities of western Antioquia (Fig. 1 in supplementary material).
The niche similarity tests (Fig. 5) using both sets of environmental predictors suggest a great variability in the environmental space inhabited by native and exotic populations (D = 0.01; I = 0.03 for both environmental predictor sets). Further, both results were not statistically significant (p-value > 0.05). However, niches similarity is low, and the difference can be explained by random chance (null models) following the environments situated between native and invasive areas.

Discussion
The reduced projected suitable conditions in the exotic area can be explained by the low similarity between native and exotic distribution niches (Fig. 4a-c), showing a potential niche shifting process between T. v. callirostris native and exotic populations. Niche differences may be highly associated with the establishment of populations at upper elevations (1300-2500 masl) in Quindío and Bogotá (Adames-Jiménez et al. 2018) where lower mean temperature zones occur than those found in the native distribution below ~ 500 masl (Rueda-Almonacid et al. 2007;Bock et al. 2012). This suggest that T. v. callirostris presents a wider fundamental niche (physiological tolerances) than the environmental space occupied in its native distribution area, allowing the species to reach different climatic conditions in the exotic area. These findings have also been reported for other invasive species like Trachemys scripta (e.g. Espindola et al. 2019). Then, it is possible that some behaviors in T. v. callirostris individuals, e.g., longer time performing aerial basking for thermoregulation, are used for establishment in relatively higher latitudes and altitudes, as has been well documented in the invasions of the Emydid turtle Trachemys scripta (Polo-Cavia et al. 2012Lambert et al. 2013;Zhang et al. 2020).
Between the years 2005 and 2009 T. v. callirostris was the most traded turtle (Bermudez et al. 2014) and is one of the most confiscated reptiles by regional environmental corporations RECs (Suárez 2017; Castro et al. 2022). With a higher trade between different Colombian regions, there is a greater chance to inhabit new localities as found in the central area of the country (exotic area), which occurs by intentional release as unwanted pets in natural areas (Adames-Jiménez et al. 2018). However, for the establishment of populations in new areas as recorded in this study in the departments of Valle del Cauca and mainly in colder environments such as Bogota (Cundinamarca department), ecological and physiological Fig. 3 Number of predicted pixels as presences in the exotic zone using all extrapolation types. Black bars correspond to prediction of consensus models for set 4, and gray bars to set 3. E = free extrapolation, EC = clamping extrapolation, and NE = no extrapolation requirements must be considered (Allendorf, and Lundquist 2003;Spear 2018). For example, several studies with the invasive Emydid turtle Trachemys scripta elegans have identified that the generalist strategies (e.g., a wide diet and high tolerance to organic accumulation process in water bodies) of this species allow it an advantage in a new environment among other turtle species (Mauermann 1995;Ficetola et al. 2009;Díaz-Paniagua et al. 2011;Polo-Cavia et al. 2014;Spear 2018). Similar traits have been reported for T. v. callirostris like varied diet (consuming algae, seeds, flowers, fishes, mollusks, and arthropods) and the capability of reproducing and nesting in highly degraded habitats (De Vivero and De La Ossa 2018;Restrepo et al. 2006;Leguízamo-Pardo and Bonilla Gomez 2014). These generalist habits represent the potential of the species to respond to changing environmental conditions.
Our results showed the impact of illegal trade on the distribution of T. v. callirostris populations, showing that even in the same country an adjacent exotic area environmental differences can occur in short distances or scales. With no knowledge about the ecological dynamics of T. v. callirostris in the exotic area, e.g., interspecific competition, rates of survival, and population density, it is unclear if the turtle is behaving as an invasive species. Indeed, our results show the potential of the reptile to establish in exotic environments, and there is a need to determine whether its presence affects biotic and abiotic conditions of waterbodies. Further, future analysis from a niche modeling perspective and conservation management would help to understand which environmental Fig. 4 Results of the suitable areas (color lines) for establishment of the turtle in the exotic area according to A = Set 3 and B = Set 4. Gray dots = no reproduction records; green star = reproduction records. Predicted pixels are represented by CE: only from clamping extrapolation; E: only from free extrapolation; NE + EC: common pixels between No Extrapo-lation and Clamping; E + NE: common pixels between No Extrapolation and free extrapolation; E + CE: common pixels from clamping and free extrapolation; E + NE + CE: common pixels predicted using all extrapolation types. Shaded areas with light purple = Cauca basin, with light brown = Magdalena basin conditions allows the establishment success of exotic populations.
We emphasize the need of relating the information of confiscations and relocations of individuals with all the RECs of Colombia, with special efforts in monitoring the new recorded localities from this study in the departments of Valle del Cauca and Cundinamarca. These departments are areas where adults have been inhabiting natural areas but no reproduction has been identified yet. There is also a need to improve management strategies and determine routes of traffic and potential new areas for establishment of the turtle. Further, monitoring plans to estimate the viability of the exotic populations is needed, as well as new research adding climate change scenarios to predict suitable areas in future.
Funding Open Access funding provided by Colombia Consortium.

Declarations
Conflict of interest No conflicts of interest are presented in this manuscript. All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.
Data availability All data generated or analyzed during this study are included in this article [and its supplementary information files].

Fig. 5
Niche similarity test between natural and exotic population of T. v. callirostris for set 3 (a and c) and set 4 (b and d).
A and B represents the native and exotic environmental space is represented in gray and red, respectively, while the environ-mental overlap is represented in blue. Dashed and solid lines correspond to the 50 and 100%, respectively, of the environmental space within native and exotic areas. C and D, shows the significance test for similarity test for both set predictors Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.