Temperature and Prey Species Richness Drive the Broad-Scale Distribution of a Generalist Predator

: The ongoing climate change and the unprecedented rate of biodiversity loss render the need to accurately project future species distributional patterns more critical than ever. Mounting evidence suggests that not only abiotic factors, but also biotic interactions drive broad-scale distributional patterns. Here, we explored the effect of predator-prey interaction on the predator distribution, using as target species the widespread and generalist grass snake ( Natrix natrix ). We used ensemble Species Distribution Modeling (SDM) to build a model only with abiotic variables (abiotic model) and a biotic one including prey species richness. Then we projected the future grass snake distribution using a modest emission scenario assuming an unhindered and no dispersal scenario. The two models performed equally well, with temperature and prey species richness emerging as the top drivers of species distribution in the abiotic and biotic models, respectively. In the future, a severe range contraction is anticipated in the case of no dispersal, a likely possibility as reptiles are poor dispersers. If the species can disperse freely, an improbable scenario due to habitat loss and fragmentation, it will lose part of its contemporary distribution, but it will expand northwards.


Introduction
Extended human-induced changes have given rise to a more "fragile" world which is getting hotter, with changing precipitation and more frequent extreme weather events [1], e.g., there was an increase of 0.7 • C in average temperature in the last 100 years [2], and a further increase of at least 1.5 • C is anticipated by 2052 [3]. Humans feel the climate change effect through everyday weather [4] and its socio-economic impacts [5,6], while climate change is related to global migration [7,8]. The rate of biodiversity loss is alarming and species have or will be summoned to cope with new climatic conditions by formed new communities. In this context, species shift their behavior [9], phenology [10][11][12], or their range [13], and these changes might allow them to persist and evade extinction. Range shifts include changes in species distributions in different geographic dimensions [14][15][16], with some species suffering range contraction and others exhibiting range expansion [17][18][19]. As the biosphere changes rapidly and in various directions [20], ecologists and biogeographers strive to project the future of species and communities under climate change in order to ensure effective biodiversity conservation and human well-being [21].
Species Distribution Modeling (SDM) [22] is a useful tool to predict climate-change impacts on species distributions (among others [17,23]). SDMs correlate spatially explicit species occurrences with sets of predictors to estimate current distributions and can be projected to other scenarios over time and space [24]. The majority of studies use environmental variables as a set of predictors. This is based on the assumption that the strength of biotic interactions fades at large spatial scales and climate is the preponderant driver of species distributions [24,25]. However, a growing number of studies show that the inclusion of biotic interactions such as competition, facilitation, mutualism and predation into SDMs increase predictive accuracy of the models [26][27][28][29][30][31][32][33][34], suggesting that biotic interactions should be brought into SDMs for more accurate predictions of climatechange impacts [35,36], though caution is required in the interpretation of the observed patterns [37].
Predation is a key biotic interaction shaping structure, dynamics and functioning of ecosystems [38][39][40][41][42], while it contributes to the maintenance of biodiversity [43]. Past global climate change events resulted in species extinctions, shifted species distributions and altered predator-prey interactions, favoring the prevalence of generalists [41]. Ongoing climate change can affect the distribution, population density, behavior, morphology or physiology of both the predator and its prey [42,44,45], inducing shifts in predator-prey dynamics [46], such as disrupted or novel interactions [47] and changes in their strength and frequency [41,48]. However, the first prerequisite for predators and their prey to interact is to overlap in space and in time. Thus, climate change-induced shifts in species distributions might result in spatial or temporal mismatch [49] or overlap into new regions affecting the predator-prey interaction. Given that predator-prey interactions affect the distribution of both the predator and the prey across spatial scales [36], and the strong relationship between a predator and its prey diversity [43], it is critical to take into consideration predator-prey interactions to accurately predict their future distributions under climate change. For instance, Aragón and Sánchez-Fernández [29] found that reptile distribution and richness had a significant effect on the distribution of bird predator in Spain, while Gherghel et al. [31] found that the SDMs predicting sea kraits' distributional patterns with biotic variables at broad spatial scales were better than the ones including only abiotic variables.
In this study, we explored the role of prey species richness in estimating contemporary distribution and projecting the future distribution of the grass snake Natrix natrix. Reptiles as ectotherms depend on the local temperature to regulate their body temperature and thus, are considered more susceptible to climate change and especially to climate warming [50][51][52]. The grass snake is a trophic generalist and a widespread species ranging from north Africa across Europe up to Scandinavia and east to Lake Baikal, Russia [53]. It benefits from anthropogenic habitats [54] and utilizes anthropogenic heat sources for nesting [55]. Here, using grass snakes we explore the effects of prey species richness and dispersal on the current and future species distribution, as estimated by Species Distribution Modeling of the widespread grass snake.

Species Data
We compiled species occurrence data from the Global Biodiversity Information Facility (GBIF) [56] and Sillero et al. [57] as centroids of the atlas grid cells for the grass snake (Natrix natrix) which is widespread in Europe. Natrix natrix taxonomy was recently revised and the species was divided in N. astreptophora found in Spain, Portugal and North Africa, N. helvetica in France, Italy and Great Britain and N. natrix which is distributed east of the Rhine and in the rest of Europe [58][59][60]. Here, we used data of the distribution of Natrix natrix as was newly defined in Europe and thus, occurrences obtained from Sillero et al. [57] were filtered to include species occurrences east of the Rhine and Italy. The grass snake is considered a generalist species. It feeds on amphibians, mostly anurans [61][62][63], but its diet varies with geographical location, habitat or competitor's presence and may include also small mammals, other reptiles and fish [62,64]. Among its most common prey species are Bufo bufo, Rana temporaria and Lissotriton vulgaris [61][62][63]. Based on this literature, we compiled occurrence data from GBIF [56] for these three species and six more species that were referenced in at least one study and represent various classes (amphibians, other reptiles, small mammals) in Europe. Occurrence data for reptilesamphibians and mammals were supplemented with occurrences from Sillero et al. [57] and Mitchell-Jones et al. [65] respectively as centroids of the atlas grid cells to ensure a more accurate depiction of the species' actual distributions. We included in our analysis only georeferenced observations (machine or human) from GBIF and we further tested for invalid records (e.g., countries' centroids, museums and GBIF headquarters) in the occurrence data with the CoordinateCleaner package [66]. Furthermore, we excluded occurrences with uncertainty higher than the radius of the study's grid pixels (GBIF field "uncertainty in meters" ≥ 25 km) or outside the species' range polygons as obtained by the International Union for Conservation of Nature (IUCN) [67] to limit the inclusion of environmental conditions irrelevant to the species' niche. The comparison with the IUCN range maps was not performed for Natrix natrix as IUCN does not provide spatial data for that species. We accounted for sampling bias by selecting only one random occurrence per pixel [68]. The number of occurrences used in the SDMs, hereafter presences, varied highly per species (see Figure S1).

Abiotic Variables
We selected as environmental predictors variables related to climate, elevation and water availability to describe current and future environmental conditions (Table S1). All predictors covered the study area which was set between 34 • 55 -81 • 45 N and 12 • 0 W-68 • 57 E and were upscaled to a 50 × 50 km grid cell resolution. Specifically, we obtained 19 bioclimatic variables, averaged over the period 1970-2000, provided by the WorldClim 2.1 database [69]. Elevation data were retrieved from the European Digital Elevation Model (EU-DEM) [70]. Furthermore, we estimated the percentage coverage by lakes and the river length included per pixel using data from the CCM River and Catchment Database [71]. Multicollinearity among bioclimatic variables was initially explored by estimating Spearman's rho correlation coefficient (Table S2). We excluded from the SDM climate predictors that exhibited collinearity, since we examined species with very different physiologies (e.g., endotherms, ectotherms) and we preferred a uniform modeling approach across species [72]. Variance Inflation Factor analysis performed by the usdm package [73] showed strong collinearity between variables, and we included only bioclimatic variables scoring VIF < 5 [74][75][76] in any subsequent analysis (Table S1).
The future bioclimatic variables used in the SDMs were produced by WorldClim 2.1 [69] using downscaled climate projections. These projections were based on the global climate model developed by the Institute Pierre-Simon Laplace (IPSL-CM6A-LR) [77] as a contribution to the Coupled Model Intercomparison Project Phase 6 (CMIP6) [78]. We selected the data for the Shared Socio-Economic Pathway (SSP) 2-4.5 [79] as it represents a modest scenario ("the Middle of the Road" [80]). This scenario predicts a CO 2 emissions pick at approximately 44 Gt/year by 2040, followed by a gradual reduction to almost 27 Gt/year by 2080. SSP 2-4.5 estimates an upward trajectory of increase in mean global temperature that will reach 2.5 • C by 2080. We focused on two datasets of bioclimatic variables spanning over two time periods: (i) 2040-2060, hereafter 2060 and (ii) 2060-2080, hereafter 2080.

Biotic Variables
As a proxy for the trophic interactions of Natrix natrix with its prey we used the prey species richness in the biotic model following Gherghel et al. [31]. Specifically, the potential distribution of each prey species was modeled in the same way as for Natrix natrix for current and future conditions using only abiotic variables (see 2.4. Ensemble species distribution modeling for the methodology followed). Future prey species richness was estimated separately under the unhindered dispersal and no dispersal scenarios; that is, potential future prey distributions were used intact in the first case and limited within their current potential distributions in the latter. We estimated current and future prey species richness as the sum of the prey species' potential distributions per pixel, which was then used as a predictor in the biotic model of the grass snake.

Ensemble Species Distribution Modeling
The potential distribution of each species under current and future conditions was estimated with Ensemble Species Distribution modeling using the biomod2 package (version 3.4.6) [81]. Ensemble SDM, by combining predictions from different modeling ap- proaches, can handle prediction variability and results in increased performance [82][83][84][85][86]. We built two models following the same approach: an abiotic model including only abiotic variables, and a biotic model in which the prey species richness was also used as a predictor variable. Prior to model fitting, we selected 1000 pseudo-absences (PAs) from the background data following the surface range envelope strategy (SRE) [87]. The SRE defines the most common environmental conditions where the species occurs and the PAs are drawn randomly outside this envelope. We fitted four algorithms generally classified either as regression-based or machine learning approaches, namely Generalized Linear Models (GLM), Generalized Additive Models (GAM), Flexible Discriminant Analysis (FDA) and Generalized Boosted Models or Boosted Regression Trees (GBM) with the settings shown in Table S3. Each algorithm was run ten times with 70% of the presence-absence data randomly selected, while the other 30% was set aside for cross-validation. The produced models were assessed using the True Skill Statistic (TSS) criterion that ranges between −1 and 1 (perfect score) and reflects the ability of a model to predict correct presences, i.e., sensitivity, and absences, i.e., specificity [88]. We also estimated the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) plot which is a threshold-independent metric that ranges between 0 and 1 with 0.5 indicating that the model discriminates presences randomly and 1 perfectly [89].
The ensemble model for each species was built as a total consensus model using all models and cross-validation runs that satisfy a predetermined evaluation threshold of TSS ≥ 0.7 or, in case of all having TSS < 0.7, with 10% of model runs with the highest TSS [88]. The ensemble model predicted the probability of the presence of each species as the weighted mean of probabilities for the current and future environmental conditions. Specifically, the probability of the presence of each model run complying with the aforementioned criterion contributed to the ensemble probability of presence according to its TSS score. The continuous probability of the presence of the ensemble model was finally transformed to binary using a cut-off value that maximized TSS in order to estimate species vulnerability [90,91]. An overview of the estimation method of grass snake distribution is presented in Figure 1. We estimated the coefficient of variation of the ensemble models and the Multivariate Environmental Suitability Surface (MESS) from dismo package (version 1.3-3) [92] to examine the uncertainty related to the extrapolation outside the range of current environmental conditions used to fit the models. Finally, we evaluated the performance of the SDMs for N. natrix with null models as proposed by Raes and ter Steege [93] to explore whether they differed significantly from random expectations. Specifically, we created 99 null model repetitions for each model type (abiotic and biotic) by randomly selecting points from the study area equal to the number of presences and pseudo-absences used in the real models. We created a frequency histogram with the AUC scores of the null models and assessed the 95% CI upper limit of these values, i.e., selected the 95th AUC value of the ranked scores. All null models ran with the same settings as the real models. The real models were considered significantly better than random if their evaluation scores were above the corresponding 95% CI AUC value of the null models.

Species' Vulnerability
We adopted commonly used indices to quantify species' vulnerability to future environmental conditions (among others [94,95]). Specifically, we estimated the relative change in area of potential distribution (change in suitable area, CSA) between current and future time periods as: and the relative loss of area (loss of suitable area, LSA) as: where current SA is the number of pixels with currently suitable habitat, future SA is the number of pixels with suitable habitat in the future, and future SA U current SA is the number of pixels with suitable habitat in the present that persists in the future. All data preparation and analysis were performed in R version 4.0.3 [96].
Diversity 2021, 13, x 5 of 16 Figure 1. Summary plot of the workflow followed to estimate potential distribution of Natrix natrix by Species Distribution Modeling. The future distribution of the species was projected using future environmental data in the abiotic model, and environmental data and projected prey species richness in the biotic model, assuming unhindered and no dispersal of Natrix natrix and its prey species.

Species' Vulnerability
We adopted commonly used indices to quantify species' vulnerability to future environmental conditions (among others [94,95]). Specifically, we estimated the relative change in area of potential distribution (change in suitable area, CSA) between current and future time periods as: and the relative loss of area (loss of suitable area, LSA) as: where current SA is the number of pixels with currently suitable habitat, future SA is the number of pixels with suitable habitat in the future, and future SA U current SA is the number of pixels with suitable habitat in the present that persists in the future. All data preparation and analysis were performed in R version 4.0.3 [96].

Performance of the Grass Snake and Prey Species Models and Variable Contribution
The abiotic and the biotic model for the estimation of the grass snake's current potential distribution were similar in terms of quality as measured by the TSS, AUC and specificity (Table 1). However, the biotic model outperformed slightly the abiotic in sensitivity ( Table 1). The ensemble models and the individual models predicting the distribution of prey species performed adequately according to TSS and AUC scores (Table S4, Figure S2). Model uncertainties fluctuated among species (the coefficient of variation ranged between 2.15 and 66.47, Tables 1 and S4) and in most cases, except for Bufo bufo, uncertainty increased in the Iberian Peninsula and northern Europe ( Figure S3). The bio- Figure 1. Summary plot of the workflow followed to estimate potential distribution of Natrix natrix by Species Distribution Modeling. The future distribution of the species was projected using future environmental data in the abiotic model, and environmental data and projected prey species richness in the biotic model, assuming unhindered and no dispersal of Natrix natrix and its prey species.

Performance of the Grass Snake and Prey Species Models and Variable Contribution
The abiotic and the biotic model for the estimation of the grass snake's current potential distribution were similar in terms of quality as measured by the TSS, AUC and specificity (Table 1). However, the biotic model outperformed slightly the abiotic in sensitivity ( Table 1). The ensemble models and the individual models predicting the distribution of prey species performed adequately according to TSS and AUC scores (Table S4, Figure S2). Model uncertainties fluctuated among species (the coefficient of variation ranged between 2.15 and 66.47, Tables 1 and S4) and in most cases, except for Bufo bufo, uncertainty increased in the Iberian Peninsula and northern Europe ( Figure S3). The biotic model for the grass snake exhibited lower uncertainty than the abiotic both as a mean value and across the study areas (Table 1, Figure S3). The MESS maps for Natrix natrix occurrences showed that the extrapolation of model predictions was within the range used to calibrate the model, suggesting reliable predictions in most of the study area ( Figure S4). Only in a few areas, mainly at areas of high elevation like the Alps mountains and areas with high aridity like the Iberian Peninsula, extrapolation exceeded this range ( Figure S4). Environmental differences between model training and extrapolation were smaller for the abiotic model than the biotic in central and southern Europe ( Figure S4). The AUC value of both ensemble models for Natrix natrix was higher than the 95% CI upper limit of the AUC values of the null models suggesting that both models differed significantly from random expectations compared to null models (Table 1). ) for the specified cut-off value, the coefficient of variation of the ensemble model prediction averaged over the area of prediction, the 95% CI upper limit of AUC values by the null models, range size (number of pixels), the number of occurrences separately for GBIF and the European atlas [4] originally in the datasets, the number of presences kept in the modeling procedure and the variables' importance (%) for the abiotic and biotic ensemble model for N. natrix rescaled to sum up to 100. The three most important variables are highlighted in bold. The annual mean temperature was the main predictor of grass snake distribution in the abiotic model (~53%) according to variable importance scores, with the species probability of presence showing a hump-shaped relationship with temperature (Table 1, Figure S5). On the other hand, in the biotic model, prey species richness had the greatest effect on the estimated grass snake distribution (~45%) with the probability of species presence increasing with increasing prey species richness (Table 1, Figure S5). In both models, isothermality contribution on predicting the potential distribution of grass snake was relatively high (~32% and~26% in the abiotic and biotic model, respectively, approximately hump-shaped relationship) (Table 1, Figure S5). The cumulative contribution of temperature, prey species richness and isothermality was~85%. In both models, the contribution of any of the remaining variables was smaller than 10% (Table 1). Species response to all variables was similar across model algorithms (FDA, GAM, GBM, GLM) and in both models ( Figure S5). Regarding prey species, temperature and isothermality were the strongest predictors of their distribution in the majority of the species (Table S5).

Potential Distribution of Grass Snake and Its Prey under Current and Future Conditions
The abiotic and biotic models predicted similar potential distribution of grass snakes under current conditions, with the spatial overlap of their estimations being approximately 91% ( Table 2). The current potential distribution of the species according to the two models is presented in Figure 2a, and its future projected distribution under different dispersal scenarios in Figure 2b-e. The abiotic model predicted approximately 3% of additional pixels suitable for grass snake and the biotic 6% ( Table 2). The inclusion of prey species richness resulted in a slight increase in the potential distribution of the grass snake (Table 1), with its potential distribution spanning a bit further in the south and east of Europe compared to the abiotic model estimation (Figure 2a). Both models predicted presence of the species in Italy and France as predictor variables did not include genetic information based on which N. natrix was divided in N. natrix and N. helvetica [60,97]. Furthermore, the discontinuities of the species distribution in the Balkan Peninsula are possibly impelled by the higher elevation in these areas and the associated decrease of mean annual temperature. The two models exhibited similar patterns of spatial overlap with prey species richness (Figure 3a). The highest overlap was observed in central Europe, the northern Balkans, and Sweden, but at least one prey species is present in almost all of Europe ( Figures S6a and S7 for the predicted distribution of each species individually). The spatial overlap of the grass snake's potential distribution with its prey under future environmental conditions is presented in Figure 3b-e. Table 2. Spatial overlap (%) of Natrix natrix potential distribution for current conditions and future projections as predicted by both models and dispersal scenarios and the respective mismatch for each model.       In the no dispersal scenario, the two models predict similar future potential distribution of the grass snake ( Table 2). The species distribution in central Europe becomes more discontinuous with time according to both models (Figure 2d,e). However, the grass snake is anticipated to lose parts of its current distribution in the east, according to the biotic model, and in the south, according to the abiotic model (Figure 4d,e). Furthermore, the biotic model estimates larger losses than the abiotic, i.e., almost 26% and 42.5% of its current distribution by 2060 (Figures 2d and 4c) and 2080 (Figures 2e and 4d) respectively. In this case, prey species were also subjected to the no dispersal limitation and are expected to suffer range contraction. In the future, the grass snake overlap with its prey differs between models (Figure 3d,e). The abiotic model predicts losses of grass snake distribution in the south, where prey species richness is anticipated to be higher compared to eastern Europe (Figures 2d,e and S6c,e), and thus, lower co-occurrence of predator with high prey species richness is expected (Figure 3d,e). On the other hand, the biotic model predicts losses of both predator and prey species in the same areas, i.e., in the east (Figures 2d,e and S6c,e), resulting in an increased overlap of the predator with at least four prey species (Figure 3d,e).

Abiotic
In the unhindered dispersal scenario, there was a slightly lower disagreement in the predictions of the future potential distribution of grass snake provided by the abiotic and biotic model than the no dispersal scenario ( Table 2). In this case, the grass snake still counts losses south and east depending on the model but expands further to the north according to both models (~3% increase and 3% decrease in range size by 2080 for the abiotic and biotic model respectively; Figures 2c and 4b). In the biotic model, the species exhibits patchy distribution in the east, while discontinuities are predicted by the abiotic model in the south, analogous to the no dispersal scenario (Figure 2b,c). Apart from the northward expansion expected for the predator and its prey in the unhindered dispersal scenario (Figures 2b,c and S6b,d), prey species richness distribution within the predator distribution was similar in both dispersal scenarios (Figure 3b-e).

Discussion
Species distributions at broad scales are considered to be driven by abiotic factors, where biotic interactions, that manifest at local scales and diminish at broad scales, add only noise to SDMs [25]. However, a series of studies have demonstrated that biotic interactions affect species distributions across spatial scales [26,27,[31][32][33]. Our results showed that the Species Distribution Model of Natrix natrix, including prey species richness, had similar performance with the one including only abiotic variables. Both models performed well, i.e., differentiated strongly species presences and absences.
The similar performance of biotic and abiotic SDM might be due to the type of interaction explored here, i.e., predator-prey, and the feeding ecology of the grass snake, i.e., a generalist snake. The effect of biotic interactions across scales depends on the type of interaction, and a positive interaction is more probable to scale up than a negative interaction [98,99]. Predator-prey is a consumer-resource interaction, meaning that is positive for the predator and negative for the prey, and if it scales up it depends on the net effect of the interaction [98]. Furthermore, the role of biotic interactions on species distribution has been related to the degree of a species' specialization [26,100,101]. The grass snake's trophic generalization might explain why the inclusion of prey species richness did not improve the SDM. Studies that reported that predator-prey interaction is imprinted on the species distribution at broader spatial scales focused primarily on trophic specialists, e.g. [29,31]. For instance, Gherghel et al. [31], who proposed the incorporation of prey species richness into SDMs, used as target species sea kraits, i.e., snakes that are trophic specialists, and found that the inclusion of prey species richness improved the model's quality only slightly.
Temperature was the top driver of the distributional pattern of grass snake in the abiotic model, while isothermality also played a role. Reptiles depend strongly on the local temperature to regulate their body temperature [102], with other aspects such as metabolic rate, sex determination, and nesting depending on temperature as well. Temperature is a strong limiting factor in the northern part of their distribution (cold physiological limit) [103,104] and other factors, e.g., biotic interactions drive the southern limit of their distribution [105,106]. It seems that the presence of the grass snake is favored up to a mean annual temperature ranging approximately 4-11 °C (this depends on the

Discussion
Species distributions at broad scales are considered to be driven by abiotic factors, where biotic interactions, that manifest at local scales and diminish at broad scales, add only noise to SDMs [25]. However, a series of studies have demonstrated that biotic interactions affect species distributions across spatial scales [26,27,[31][32][33]. Our results showed that the Species Distribution Model of Natrix natrix, including prey species richness, had similar performance with the one including only abiotic variables. Both models performed well, i.e., differentiated strongly species presences and absences.
The similar performance of biotic and abiotic SDM might be due to the type of interaction explored here, i.e., predator-prey, and the feeding ecology of the grass snake, i.e., a generalist snake. The effect of biotic interactions across scales depends on the type of interaction, and a positive interaction is more probable to scale up than a negative interaction [98,99]. Predator-prey is a consumer-resource interaction, meaning that is positive for the predator and negative for the prey, and if it scales up it depends on the net effect of the interaction [98]. Furthermore, the role of biotic interactions on species distribution has been related to the degree of a species' specialization [26,100,101]. The grass snake's trophic generalization might explain why the inclusion of prey species richness did not improve the SDM. Studies that reported that predator-prey interaction is imprinted on the species distribution at broader spatial scales focused primarily on trophic specialists, e.g., [29,31]. For instance, Gherghel et al. [31], who proposed the incorporation of prey species richness into SDMs, used as target species sea kraits, i.e., snakes that are trophic specialists, and found that the inclusion of prey species richness improved the model's quality only slightly.
Temperature was the top driver of the distributional pattern of grass snake in the abiotic model, while isothermality also played a role. Reptiles depend strongly on the local temperature to regulate their body temperature [102], with other aspects such as metabolic rate, sex determination, and nesting depending on temperature as well. Temperature is a strong limiting factor in the northern part of their distribution (cold physiological limit) [103,104] and other factors, e.g., biotic interactions drive the southern limit of their distribution [105,106]. It seems that the presence of the grass snake is favored up to a mean annual temperature ranging approximately 4-11 • C (this depends on the model used to estimate its distribution), and also from smaller levels of isothermality, perhaps in order to retain activity. Therefore, the grass snake can survive in a wide array of thermal conditions, as it is also reflected by its wide distribution, i.e., from warmer conditions, like those prevailing in the Mediterranean, to cold, e.g., Scandinavia. Although the biotic SDM did not perform better than the abiotic, the prey species richness was the most important driver of grass snake distribution, contributing approximately 45% to the SDM. The distribution of the majority of prey species and the predator (according to the abiotic model) is driven by similar predictors (temperature and isothermality). These species seem to thrive in the same areas dictated by their thermal biology, and this perhaps can explain the similar performance of abiotic and biotic SDMs, i.e., although the species interact, the similar environmental niches might partially mask the imprint of biotic interaction on the broad-scale distribution of this generalist predator.
Both models predicted shifts of the grass snake's distribution under the climate change scenario assuming unlimited dispersal, and losses assuming no dispersal. Under the no dispersal scenario, the species is projected to suffer range contraction by 2080, while under the unhindered dispersal scenario, the losses of current distribution will be compensated by expansion to the north. According to the abiotic model with no dispersal, range contraction is inevitable in the face of climate change with the species losing a large part of its current distribution in the Balkan Peninsula by 2060 and approximately 37% of its current distribution by 2080. However, if the species can disperse freely, it is expected to expand into new areas northwards and eastwards; in the latter case, losing simultaneously to the south until 2080. Note that the unhindered dispersal scenario used here is in congruence with the limited dispersal assumed for amphibians by Thuiller et al. [107], i.e., 20 km/year, and a suggestion made by Smith and Green [108], i.e., 10 km/year for amphibian dispersal ability, although reptiles are considered poor dispersers [17,109] with reported maximum dispersal of 3 km to our knowledge [110]. In the case of a biotic model, grass snake is anticipated to retain the southern parts of its distribution, but to become extinct in the east. Assuming no dispersal, the species will suffer a range of contraction (similar to the abiotic model), but it will expand northwards if it can disperse. The two models exhibited differences regarding the distribution in the south which were sharpened by 2080. Specifically, the biotic model predicted a distribution with a further southern edge than that of the abiotic model. Perhaps, although future climatic conditions exceed the physiological limits of the species there, the presence of some prey mitigates them. This is in congruence with the finding that the southern range limits of amphibians and reptiles are determined by species interactions (see Cunningham et al. [105] and references therein).
Snakes and reptiles, in general, are considered to globally decline [51,109] due to climate change and induced changes (e.g., habitat loss, prey availability), while regional declines in grass snakes have been reported [111]. Even assuming that the species can disperse unhindered based on their dispersal ability, habitat loss and fragmentation and land-use changes, which are among the main drivers of the global decline of amphibians and reptiles [112], might render the expanded suitable space inaccessible [113]. Therefore, the dispersal ability of species and connectivity will crucially affect the future distribution of the species [17,114]. It seems that climate warming along with the ongoing habitat fragmentation [115] minimizes the reptiles' chance to move to areas with suitable environmental conditions [17,110]. This is vividly depicted by the range contraction projected by the SDMs assuming no dispersal. Furthermore, global warming might favor reptiles by increasing the available suitable area for them but can negatively impact phenology-related aspects resulting in higher metabolic rates and skewed sex ratio due to their temperaturedependent sex determination [50,115,116]. Although adaptations to these changes have "rescued" reptile species from historic climate warming events, a timely adaptation is not the most probable scenario [50,[117][118][119].

Conclusions
We found a similarly good performance of the biotic and the abiotic Species Distribution Model predicting the distribution of the generalist Natrix natrix. However, the prey species richness was the top driver of grass snake distribution, highlighting the role of the prey availability in shaping species distribution. Temperature and isothermality played a significant role in shaping the grass snake distribution in all cases. The abiotic and the biotic model predicted a northward shift of the species distribution under climate change assuming unlimited dispersal, but severe range contractions in the case of no dispersal. Our results highlighted that biotic interactions and dispersal ability are game-changers in future projections of species distributions and should be included in SDMs. However, they differ widely among taxonomic groups and regions and cannot be collectively modeled for different taxonomic groups. Thus, more taxon and region-specific studies are needed if we want to predict a species' response to climate warming and develop effective conservation strategies.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13040169/s1, Table S1: The predictors considered for the species distribution modeling. Table S2: Spearman's rho correlation coefficients between bioclimatic variables. Table S3: The models used for the ensemble modeling and the associated parameter settings. Table S4: The total number of occurrences and separately for GBIF and the European atlases originally in the datasets, the number of presences kept in the modeling procedure, True Skill Statistic (TSS), Area Under the Curve (AUC), sensitivity and specificity scores for the presented cut-off value for each prey species and the coefficient of variation of the ensemble model prediction for each species averaged over the area of prediction. Table S5: Variable importance (%) for the prediction of each prey species probability of presence in the abiotic ensemble species distribution model., Figure S1: Presences used in the modeling process for each species. Figure S2: Scores of True Skill Statistic (TSS) and Area Under the Curve (AUC) for each individual model, prey species and model type (abiotic and biotic) for Natrix natrix. Figure S3: Coefficients of variation (standard deviation/mean) of the predictions for the ensemble species distribution model of each species and model type (abiotic and biotic) for Natrix natrix. Figure S4: Multivariate environmental suitability surface maps (MESS) of the predictors included in the abiotic and the biotic model for current environmental conditions for Natrix natrix. Figure S5: Predicted response curves from the different algorithms that were used in the individual and ensemble models for Natrix natrix. Figure S6: Predicted prey species richness under current environmental conditions and by 2060 and 2080 respectively according to unhindered dispersal scenario and no dispersal scenario. Figure   Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://www.gbif.org/; http://na2re.ismai.pt/; https://www.european-mammals.org/ php/mapmaker.php; https://www.iucnredlist.org/; https://worldclim.org/; https://www.eea. europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem; https://ec.europa. eu/jrc/en/publication/eur-scientific-and-technical-research-reports/ccm-river-and-catchment-databaseversion-20-analysis-tools.