Fundamental niche unfilling and potential invasion risk of the slider turtle Trachemys scripta

Occurrence and climatic data

Considering the overlapped distribution of Trachemys scripta subspecies (T.s. scripta. T.s. elegans and T.s. troostii), and the lack of formal discrimination of the turtle’s invasive introductions on the basis of subspecies (Van Dijk, Harding & Hammerson, 2011), we obtained georeferenced records for the three subspecies from the Global Biodiversity Information Facility (GBIF, 2017) and the VertNet (2016) databases. To complement the latter, we also included our own collecting records and those from a thorough revision from the available published literature (data available at https://figshare.com/articles/_/8175158). In order to define T. scripta’s native range that encompasses from eastern and central United States of America to extreme north-eastern Mexico, we used as a basis the range maps proposed by Seidel (2002) about taxonomic observations on extant species and subspecies of Trachemys, updated in Rhodin et al. (2017); next, we used the Freshwater Ecoregions of the World (WWF/TNC, 2013) to define range limits based on a global freshwater regionalization. Data points outside these limits were considered as non-native records; that is where the species has been introduced and is considered invasive worldwide (including North America where it has also been introduced outside its native range). We applied a geographic filter to our full set of occurrence data, considering records further than 5 km distance as independent points (multiple points within a cell indicate abundance is higher on those sites in comparison to cells with fewer occurrences). To that end, we used the Ecospat package in R (Di Cola et al., 2017) to remove duplicated points within a 5 km buffer, thus reducing overfitting and improving predictions (De Oliveira et al., 2014). As a result, we have the most comprehensive data set for the slider turtle, with 2,552 occurrence records, 1,395 from the native distribution range and 1,157 from non-native occurrences around the globe, more than four times the number of records from previous studies.

Climatic data were obtained from the global climate layers of WorldClim (Fick & Hijmans, 2017), which represents a statistical summary of temperature, precipitation and radiation at 5′spatial resolution. We included the six following bioclimatic variables: annual mean temperature (BIO1), mean temperature of the warmest quarter (BIO10), mean temperature of the coldest quarter (BIO11), annual precipitation (BIO12), precipitation of driest month (BIO14), and annual mean radiation (RAD), which were selected on the basis of the slider turtle’s natural history (i.e., physiology), and also to reduced correlated variables (Pearson’s correlation coefficients <0.85). Trachemys slider turtles are found in shallow, slow-moving water that has diverse vegetation and nearby places to bask. During seasonal dryness turtles may wander far seeking water sources, rapidly colonizing any newly available habitat; hence, waterbodies availability is a key feature for the activity, dispersal, distribution, and overall life history and survival of this species (Cox et al., 1998; Minton, 2001; Buhlmann, Tuberville & Gibbons, 2008). Accordingly, we constructed a layer depicting the Euclidian distances to the nearest freshwater body (NearDist) (data available at https://figshare.com/articles/_/8175158). For this, we used the Global Lakes and Wetlands Database (GLWD; WWF, 2017), a combination of best available sources for lakes and wetlands on a global (world) scale. This raster map has a resolution of 30-second that includes polygons of great lakes (area ≥ 50 km²), great reservoirs (capacity ≥ 0.5 km³) and permanent water bodies (area ≥ 0.1 km²). This is the best resolution available and, although smaller water bodies might not be directly detected, nonetheless polygons encompass the maximum extension of different kinds of water bodies (e.g., marshes, flooding areas). Importantly, this distance to the nearest freshwater body is an indirect measure of water availability—not a measure of dispersal—hence it defines a characteristic of the fundamental niche (Anderson, Gutierrez & Romano, 2002; Braun & Phelps, 2016).

We tested our proposed model at a regional scale, for Mexico, based on the same variables described above but with a higher resolution (30 arc sec). We used climatic surfaces for the average monthly climate period 1910–2009 (Cuervo-Robayo et al., 2014), jointly with a layer we built depicting the Euclidian distances to the nearest freshwater body for Mexico (NearDistMex) using the National Hydrographic System catalogue (SGM, 2019).

Comparisons between native and non-native realized niches

We used different approaches to assess similarities and differences between the climatic conditions occupied by T. scripta in its native and non-native ranges, all performed in R 3.0.1 (R Core Team, 2015). First, in order to characterize the environments available worldwide, we used as reference the occurrence records of all freshwater turtles (not only Trachemys) from GBIF (2017), which represent the environments that turtles could occupy. However, because some regions are better sampled for turtles than others and we cannot know if the lack of data is because the environment is not suitable or because they have not been sampled, we constructed a probability layer based on sampling effort by means of assigning a higher probability value where higher occurrences are recorded. We then used the dismo package to create 100,000 random points worldwide weighted by our probability value. Next, based on the bioclimatic variables of those random points and of the 2,556 T. scripta occurrence records (data available at https://figshare.com/articles/_/8175158), we performed a weighted principal components analysis (PCA). Finally, using the first two axes of this PCA we: (i) draw contours estimated with a kernel density function (Broennimann et al., 2012) to delimit the climatic conditions available on the native range, the non-native ranges (in America, Europe, Asia and Australia) and around the world (areas that have not been invaded); and (ii) draw points depicting the estimated realized niche of the native and non-native ranges. This procedure randomizes the position of the kernel density surface of one of the two species within the environmental space available for it, allowing to evaluate whether the realized niches (kernel density surface) occupied in native and non-native ranges are more or less similar to the distribution of similarities under a null model.

As a complementary test, we compared climatic covariance matrices (native versus non-native ranges) with a Common Principal Components Analysis (CPCA; Phillips & Arnold, 1999), which compares two or more matrices considering their eigenvectors and eigenvalues in a hierarchical fashion, to describe their structure in relation to the size, shape and orientation of the matrices. We used the latter to establish principal components that were common between matrices (i.e., progressive differences in shape, orientation and size) and to test hypotheses about equality (identical eigenvectors and eigenvalues, i.e., identical size, shape and orientation), proportionality (equal eigenvectors, but eigenvalues differing in a scalar amount, i.e., same shape and orientation, different but proportional size), and unrelated structure (matrices have dissimilar eigenvectors and eigenvalues). The best solution under the model-building approach of the CPCA of the comparison of climatic covariance matrices is indicated by the minimum value of Akaike information criterion (AIC; Phillips & Arnold, 1999).

The statistical framework proposed by Broennimann et al. (2012) was used to evaluate the assumption of niche conservatism in biological invasions by quantifying the similarity between the native and the non-native realized niches. These comparisons involve a multivariate analysis to calibrate the niche and the occurrence density, performed with Ecospat, and the estimation of niche overlap using Schoener’s D metric (Schoener, 1970), an index ranging from 0 (no overlap) to 1 (total overlap) (Warren, Glor & Turelli, 2008; Broennimann et al., 2012), with statistical significance assessed on the basis of 100 randomizations (α = 0.05). Furthermore, the density of occurrences in environmental space was used to estimate niche expansion (new environmental conditions found in the non-native range), stability (proportion of the native niche conditions found in the non-native one), and unfilling (proportion of the native niche not occupied in the non-native) (Guisan et al., 2014), with Ecospat. Finally, since the previous method only allows the use of two principal components, for a high-dimensional niche overlap analysis we used a technique based on a multivariate kernel density estimation, implemented in the Hypervolume in R (Blonder et al., 2014), that estimates densities based on a Monte Carlo importance sampling approach. To define the six-dimensional hypervolumes (i.e., the six bioclimatic variables described above) for the native and non-native ranges, we chose a threshold that included 100% of the total probability density. We then estimated the overlap (intersection) between the two hypervolumes (native and non-native), and the hypervolume unique to each native or non-native range; results are depicted as pair-plots.

Environmental tolerances, fundamental niche and potential distribution

We performed a literature review (Table S1) to obtain information about physiological optimums and tolerance limits for T. scripta. The environmental variable reported more often and consistently was temperature for the following five features: during eggs incubation, hatching, growth, basking, and activity periods. The proximity to freshwater bodies is another critical limit, a fundamental factor for the turtles’ survival which they depend upon to complete their life cycles, feeding and reproduction (Anderson, Gutierrez & Romano, 2002; Braun & Phelps, 2016); thus, we used the Euclidian distance to the nearest freshwater body (NearDist) as an indirect measure of water availability, as explained above. We defined a 10 km buffer for distance to water availability, considering the general 5–10 km average distance reported where turtles can explore (Cagle, 1944) and, therefore, reach freshwater bodies. Hence, we considered the combined use of these different temperature tolerance limits and freshwater availability as an accurate approximation of the fundamental niche, as bivariate environmental space, for the slider turtle. Accordingly, we draw response curves in R by using beta and logistic functions adjusted to the optimal, maximum and minimum tolerance limits of T. scripta, with Virtualspecies in R (Leroy et al., 2015), which were then used to simulate its potential geographic distribution (see Glossary S1). The tolerance limits used, based on our literature review, were: optimum temperature for eggs incubation (29.5 °C) and daily activity (25.6 °C); ranges for feeding (16.8–32.1 °C), normal activity (10–37 °C), and basking (26.7–31.2 °C); and minimum and maximum tolerances before death (−12.6–42.3 °C). Virtualspecies allows to define response functions for different predictor variables and to combine the responses to obtain a habitat suitability value. This approach increases the ecological realism for generating environmental suitability of virtual species and to create more adequate species distribution models.

Next, to compare the range of the fundamental niche that is actually occupied by T. scripta with the simulated results, kernel density curves based on the temperature variables from all occurrence data (taking into account native and non-native records from which we extracted bioclimatic data) were overlapped with the simulated response curves. Contours and points were also drawn as described above (Broennimann et al., 2012), but specifically considering temperature data from the annual mean temperature (BIO1), the mean temperature of the warmest quarter (BIO10), and the mean temperature of the coldest quarter (BIO11), to encompass both the temperature ranges and the extreme tolerance limits of the five described activities; the variable NearDist was also included. The 100,000 generated worldwide random points were used to delimit the bivariate environmental space available, while the occurrence points were used to characterize the bivariate environmental space occupied (data available at https://figshare.com/articles/_/8175158). This data set was compared with the available conditions worldwide, in native and non-native ranges, as well as with the species’ physiological tolerance limits. In order to evaluate if combining information about these physiological tolerance limits could improve the prediction of T. scripta’s potential invasiveness, we simulated again its potential distribution with Virtualspecies, this time based on the response curves previously estimated and the higher resolution surfaces built for Mexico.

Trachemys scripta unfilled realized and fundamental niches

The approach we followed to define T. scripta’s fundamental niche, based on temperature physiological tolerances and availability of freshwater bodies, both conditions essential for the persistence of its populations, revealed it has a wide environmental range with favourable climatic conditions (within its tolerance limits). This result is consistent with earlier works suggesting that T. scripta has a broad environmental tolerance, particularly regarding temperature ranges (Rödder et al., 2009; Kikillus, Hare & Hartley, 2010; Masin et al., 2014); remarkably, what had not yet been demonstrated is that such broad tolerance does encompass a fragment of the range of temperature optimums it requires for performing key behaviors (egg incubation, hatching, growth, basking and activity periods). Accordingly, we consider that the climatic data we used to obtain the ranges of the fundamental niche actually occupied by T. scripta allowed us to build a robust model encompassing the environmental conditions required for all these behaviors to occur (i.e., the species’ fundamental niche).

Indeed, when compared with the range that T. scripta actually occupies, we evidence that it is not experiencing the entire range of optimal conditions along its native distribution; specifically when compared based on the warmest conditions, where suitable climate is unoccupied along its higher spectrum of tolerances (Fig. 3). The latter might reflect biotic features (absent from the climatic models), namely the fact that the distribution of its congener species, Trachemys cataspila, starts at the southern limit of the slider turtle’s, occupying warmer environmental conditions; these two species do not occur sympatrically, likely due to competition (Seidel, 2002). Interestingly, we also found that although non-native areas include similar temperature ranges comparatively with the native ones, non-native individuals have invaded environments that are both warmer and colder than its native range, as documented by Ficetola, Thuiller & Padoa-Schioppa (2009) and Rodrigues, Coelho & Ribeiro (2018); however, we show that this has happened without deviating from the turtle’s tolerance ranges. This could result from elimination of dispersal barriers and of potential competitors, as shown by studies with different taxa where niche shifts (unfilling and expansion) have been associated with biotic interactions and dispersal limitations, which prevent the species from colonizing the full extent of the available (suitable) conditions (Hargreaves, Samis & Eckert, 2014; Tingley et al., 2014; Tingley et al., 2016; Strubbe, Beauchard & Matthysen, 2015); or, conversely, as factors that facilitate the establishment of invasive species, exhibiting where they can potentially be found (Heikkinen et al., 2007; Tingley, Phillips & Shine, 2011; Giannini et al., 2013). Our finding that T. scripta is also reaching lower temperatures along its invaded distribution could be associated with the trade of the species as a pet (Burger, 2009), where repeated introductions might be forcing the species to occupy areas that are not as optimal (Ficetola, Thuiller & Padoa-Schioppa, 2009; Zhu, Li & Zhao, 2017). For instance, Ficetola, Thuiller & Padoa-Schioppa (2009) recorded T. scripta individuals in northern Europe, but where conditions are likely too harsh for successful breeding.

The evaluation of native and non-native distributions provides key information in order to decipher if observed shifts might simply indicate different portions of the fundamental niche occupied by the species (Araújo & Peterson, 2012), or suggest an adaptation process to new environmental conditions (Hill, Chown & Hoffmann, 2013). When we contrasted T. scripta’s native and non-native realized niches along the environmental space, we found a shift in shape, orientation and size (climatic matrices). In fact, the hypervolume analysis results suggest that the introduction of T. scripta into new areas has facilitated the use of available climates in a different way. For instance, while our findings show that temperature is the most significant variable determining the areas that are occupied by the turtle in its native range, it is precipitation at its non-native range. We acknowledge that the hypervolume method’s assumptions have been criticized (Qiao et al., 2017), nonetheless our results are consistent with the fact that the species always occupies a set of environmental conditions within its physiological tolerance limits and encompassing both analog and non-analog conditions. Furthermore, our results hence support the hypothesis about invasive species not in equilibrium with the environment at their native range (Elith, Kearney & Phillips, 2010; Owens et al., 2013), and ratify the importance of fully considering the species’ capacity to use a wider climate gradient on both space and time (Broennimann & Guisan, 2008; Václavík & Meentemeyer, 2012; Tingley et al., 2014). Accordingly, we show that the slider turtle’s documented niche expansion (Li et al., 2014; Rodrigues et al., 2016) in fact represents a different subset of the species’ fundamental niche.

More complex ways when defining potential invasion risk

Notably, our results show both worldwide suitable conditions (in America, Europe, Asia, Australia) as well as regional areas (Mexico) that have not yet been invaded despite encompassing the range of climatic conditions where T. scripta has successfully invaded (Figs. 4–6). In addition to the relationship between the distribution of organisms and climate availability, there are more complex approaches to define potential invasion risks (Gallien et al., 2010; Qiao et al., 2017), which could also explain why T. scripta does not occupy the total area with suitable environment in its non-native (and neither native) distribution. Václavík & Meentemeyer (2012) have pointed out the need to consider the stage of the invasion to avoid underestimation of habitats at risk of invasion. The history of invasiveness of T. scripta is relatively recent: it has been transported globally since the 1950s as one of the most popular turtle pets (Burger, 2009), facilitating its introduction into new habitats. Considering the long generation times of these turtles, it is highly likely that T. scripta can continue spreading and colonizing new portions of its fundamental niche, which is currently not being fully occupied, allowing it to respond favourably under a climate change scenario (Hellmann et al., 2008; Cosner, 2014; Early & Sax, 2014).

Ecological niche models need to be taken with caution when used to forecast species invasions and their response to environmental change, as it has frequently been shown they may underestimate the potential spread of invasive species (see Parravicini et al., 2015; Merow et al., 2017). Occurrence data from within and outside of a species native range have been amply used in studies evaluating invasion risk and niche change in exotic species (Guisan et al., 2014; Li et al., 2014). Studies have, nonetheless, shown that geographic predictions should be used carefully, especially when species are experiencing niche shifts, because it can lead to underestimation of the potential invasion risk. We believe one crucial condition to evaluate invasion risk is having the most accurate fundamental niche possible of the species in question. Biotic interactions also play an important role during the invasion process and, therefore, are key in defining the range that an invasive species can or cannot occupy. Indeed, although the slider turtle has been able to establish widely, its success has not been ubiquitously identical. One likely reason for the latter is that invasions are species- and case-specific, mainly because it involves complex biotic interactions (i.e., with competitors or predators) not easily detected by environmental modeling approaches. Thus, although still unresolved (see Wisz et al., 2013; Cunningham et al., 2016), it is essential to find means to include biotic interactions when ENM predictions are used to evaluate invasion risk. Also, for an accurate invasiveness assessment, the distinction between colonizing areas where the species can establish and reproduce from those where it can only survive needs be considered; for instance, by comparing predictions with and without those records. In our case, we consider that T. scripta ‘s observed non-native range is an adequate proxy of its invasive potential, given that it can survive for decades in areas outside its breeding requirements (Ficetola, Thuiller & Padoa-Schioppa, 2009) and that the observed non-native ranges do not deviate from its fundamental niche. Finally, spatial heterogeneity in environmental conditions is another challenge to consider when searching for better invasion risk models. Incorporating landscape genetics methods (Balkenhol et al., 2013) could aid in detecting barriers or corridors embedded in the geographic space, ultimately helping improve predictions associated with the expansion range of invasive species (Tingley et al., 2013).