- Open Access
Unraveling climate influences on the distribution of the parapatric newts Lissotriton vulgaris meridionalis and L. italicus
© The Author(s). 2017
- Received: 4 July 2017
- Accepted: 27 October 2017
- Published: 12 December 2017
Climate is often considered as a key ecological factor limiting the capability of expansion of most species and the extent of suitable habitats. In this contribution, we implement Species Distribution Models (SDMs) to study two parapatric amphibians, Lissotriton vulgaris meridionalis and L. italicus, investigating if and how climate has influenced their present and past (Last Glacial Maximum and Holocene) distributions. A database of 901 GPS presence records was generated for the two newts. SDMs were built through Boosted Regression Trees and Maxent, using the Worldclim bioclimatic variables as predictors.
Precipitation-linked variables and the temperature annual range strongly influence the current occurrence patterns of the two Lissotriton species analyzed. The two newts show opposite responses to the most contributing variables, such as BIO7 (temperature annual range), BIO12 (annual precipitation), BIO17 (precipitation of the driest quarter) and BIO19 (precipitation of the coldest quarter). The hypothesis of climate influencing the distributions of these species is also supported by the fact that the co-occurrences within the sympatric area fall in localities characterized by intermediate values of these predictors. Projections to the Last Glacial Maximum and Holocene scenarios provided a coherent representation of climate influences on the past distributions of the target species. Computation of pairwise variables interactions and the discriminant analysis allowed a deeper interpretation of SDMs’ outputs. Further, we propose a multivariate environmental dissimilarity index (MEDI), derived through a transformation of the multivariate environmental similarity surface (MESS), to deal with extrapolation-linked uncertainties in model projections to past climate. Finally, the niche equivalency and niche similarity tests confirmed the link between SDMs outputs and actual differences in the ecological niches of the two species.
The different responses of the two species to climatic factors have significantly contributed to shape their current distribution, through contractions, expansions and shifts over time, allowing to maintain two wide allopatric areas with an area of sympatry in Central Italy. Moreover, our SDMs hindcasting shows many concordances with previous phylogeographic studies carried out on the same species, thus corroborating the scenarios of potential distribution during the Last Glacial Maximum and the Holocene emerging from the models obtained.
- Last glacial maximum
- Multivariate environmental dissimilarity index
- Niche divergence
- Species distribution models
The interpretation of species distribution patterns is crucial to understand biogeographical, ecological and conservation aspects of biodiversity [1–6]. The presence of certain species in a territory can be read as the result of many factors interacting in space and time, such as historical distributions, paleoclimatic events, large and fine-scale fragmentation, biotic interactions, niche width and dispersal ability (e.g., [7–10]). Thus, modelling the distribution of species within discrete biogeographical units, especially when dealing with low-dispersal ability species [11, 12], requires a proper implementation of both Species Distribution Models (SDMs) and phylogeographic analyses based on genetic evidence [13–15]. Species’ distribution is usually constrained by biotic interactions, dispersal capability and geographic accessibility (the “B” and “M” of a BAM diagram ) or by abiotic factors (the “A” part). Within the SDMs approach, wide areas with high predicted suitability/probability of presence with respect to a certain species but currently unoccupied suggest the existence of some unmeasured variables constraining the species’ distribution . Dispersal, for instance, could be limited by geographical barriers, competition between two species with similar ecological niches, or tolerance limits to some environmental factors . As an example, cycles of glaciations during the Pleistocene, because of their limiting effects on species’ dispersal possibilities and tolerance to climatic stressors, largely contributed to shape European fauna , leading to the extinction of some species and influencing the distribution of other taxa [20–22].
In this contribution, we try to understand if and how the current distributions of the two parapatric newts Lissotriton vulgaris meridionalis (Boulenger, 1882) and L. italicus (Peracca, 1898), have been influenced by climatic conditions. These two urodeles live in lentic environments, have no strict habitat preferences but differ in some phenological traits [23, 24], and are clearly reproductively isolated [25, 26]. Our analysis is focused on the climatic conditions, particularly those related to temperature and precipitation patterns because of the great influence of these parameters on amphibian life-history traits , occurring within the current range of these two amphibians, in order to identify possible variables contributing to limit their capability of range expansion. Particular attention was given to a sympatric zone, because of the great importance accorded to these areas in many fields of theoretical and applied ecology. In fact, understanding the patterns of species’ overlapping occurrences represents a challenging topic of ecological research, ranging from the reconstruction of the species’ evolutionary history  to ecological modelling [29, 30], especially when applied to the study of the relationships between past biogeographical patterns and climate in closely related taxa [9, 12, 31]. We built SDMs for both the target species under the current conditions using two machine-learning techniques, Boosted Regression Trees (BRT) and Maxent. Then, the resulting models were also projected to past climatic conditions [21, 32] in order to infer historical habitat suitability and, consequently, to hypothesize how the climatic oscillations during the Last Glacial Maximum and the Holocene have influenced the distribution of these two Lissotriton species until they reached the current status. Finally, we investigated also possible statistically significant divergence between the environmental niches of the two target species.
Target species and study area
Presence records in the study area were drawn from the database reported in Iannella , with a total of 401 and 500 occurrences for L. italicus and L. v. meridionalis, respectively; 25 records, all falling into the sympatric area, are referred to syntopic localities shared by both species. The geographic coordinates of the presence records are reported in the Additional file 1.
Our study area comprehends the whole Apennine area and the Apulian Province, since the latter hosts two L. italicus haplogroups shared with some other L. italicus populations of Central and Southern Apennines . We did not consider, however, the portion of L. v. meridionalis’ range corresponding to the Padano-Venetian plain (sensu Canestrelli et al. ) and to the Prealps, since our aim was to focus on the climatic conditions of the Mediterranean biogeographical region where the parapatry between our two target species occurs. Moreover, we excluded the L. v. meridionalis’ populations located in north-eastern Italy because of the introgression events with L. v. graecus’ populations from northern Balkans . However, since restrictions of the environmental range in SDMs may lead to artefacts in response curves and spurious projections , we investigated if the set of environmental conditions characterizing the portion of the L. v. meridionalis’ range not included in our study area was noticeably different from the one used to calibrate our model. 170 presence and 1000 background points were randomly sampled from each of the two portions of L. v. meridionalis’ range (included and excluded); the values of the input bioclimatic variables at each point were extracted using the ‘Extract Values to Points’ tool in Arcmap 10.0 and the Pearson coefficient (r) between the resulting matrices was computed in R through the ‘cor.test’ function. Presence records for L. v. meridionalis outside the study area considered for this correlation test were drawn from .
Model building and GIS analysis
In order to build the SDMs for the target species, Boosted Regression Trees (BRT) and Maxent were used. Both are among the modelling techniques which perform best, based on statistics referring to discrimination, calibration and correlation between observations and predictions [38, 39]. Contrarily to the high number of contributions based on Maxent analyses , there are not as many studies implying BRT to address biogeographical questions, probably because it requires a more substantial effort to parameterize and build the model. Nevertheless, BRT was chosen because this technique permits to accurately model complex responses of the target species to the predictors [39, 41], and to investigate possible synergistic effects of the considered variables through the assessment of pairwise interactions. The BRT models were built in R 3.3.2 (R Core Team, 2016) using the package ‘gbm’ 2.1.1  and the ‘gbm.step’ function provided by J. Elith and J.R. Leathwick , while the Maxent models were built through the software MAXENT 3.3.3 k . Notwithstanding BRT is a “tolerant” method with respect to correlated variables , a correlation matrix considering all the variables selected as possible predictors was built (see Additional file 2), in order to exclude the ones showing Pearson | r | > 0.85 , thus avoiding any multicollinearity side effect. Bioclimatic variables of current, Mid-Holocene (MOL, ~ 6000 years ago), and Last Glacial Maximum (LGM, ~ 22,000 years ago) periods and three topographic variables were considered as candidate predictors. Nineteen bioclimatic variables were downloaded from the WorldClim.org website, version 1.4 , at 30 arc-seconds resolution, when available. LGM bioclimatic variables, available only at 2.5 min resolution, were rescaled to 30 arc-seconds in ArcMap 10.0 . For both the Mid-Holocene and the Last Glacial Maximum, two Global Climate Models (GCMs) were considered, namely the CCSM4  and the MIROC-ESM . The details on the bioclimatic variables are reported in Additional file 3. With regard to the topographic variables, altitude (ALT) was obtained through the digital elevation model of the study area, downloaded from the geo-portal of the Italian Ministry of the Environment (http://www.pcn.minambiente.it). The other two topographic variables, aspect (ASPECT) and slope (SLOPE) were derived from ALT through their respective functions in ArcMap 10.0.
A two-step pseudo-absence selection method with preliminary environmental profiling was adopted [49, 50] to generate the pseudo-absences used to build the BRT models. The environmental profiling was performed by building a Bioclim  model in R through the package ‘dismo’ 1.1–1 , for each species, considering as predictors the bioclimatic variables which resulted to be non-correlated from the above-mentioned correlation matrix. Bioclim continuous output was converted in five discrete classes in GIS environment (ArcMap, Natural breaks) and the pseudo-absences were generated by selecting points at random within the polygons corresponding to the three classes of lowest modelled habitat suitability, so as to characterize the pseudo-absences as possible true absences. The use of pseudo-absences generated after a prior identification of climatically unsuitable areas allowed us to properly model the potential distribution of the two Lissotriton species considered [49, 53]. The number of the pseudo-absences generated for each species through this procedure was the same as the number of presence records, since a ratio between presences and pseudo-absences close to one seems to assure an optimal performance of BRT models . The BRT models were built with the following setting of parameters for each species: 10-fold cross-validation; ‘tree complexity’ = 5; ‘bag fraction’ = 0.5 and ‘learning rate’ = 0.001. Maxent, instead, was parameterized as follows: ‘Auto features’; 10-fold cross-validation; 10,000 background points; 5000 iterations. Additional SDMs were built for the current scenario, through both BRT and Maxent, considering as predictors both the non-correlated bioclimatic variables and the topographic ones, in order to assess if the latter outclass the climate-linked variables in terms of contribution to the models.
In order to point out possible spatial autocorrelation of presence data before building the SDMs , a Moran’s I test on a 5 × 5 km grid, derived from the 10 × 10 km UTM grid, was performed for each species, with each cell containing spatial information of presence data. The raster values used in the statistical analysis were extracted through the ‘Extract Values to Points’ tool in ArcMap 10.0.
All spatial input data (i.e. presence records and rasters of the candidate predictors) were clipped to the extent of the whole study area and projected in WGS84 reference system.
BRT and Maxent models obtained were evaluated in their discrimination power by means of the AUC (i.e. Area Under Curve) of the ROC (Receiver Operating Characteristic) curve , a threshold-independent statistic which assesses the capability of a SDM to discriminate between presences and (pseudo-)absences [38, 43]. Secondly, calibration within the BRT models was evaluated considering the deviance score, calculated through the ‘gbm.step’ function, which represents a measure of loss in predictive performance due to suboptimal models . Both for the BRT and the Maxent models, the relative contribution of each predictor was assessed, and the partial dependence plots for the most contributing ones were successively produced. Moreover, within the BRT models, the pairwise interaction score of each possible pair of predictors was computed in order to individuate any synergistic effect between the most influential bioclimatic variables. Finally, a forward stepwise discriminant function analysis  was performed to derive functions discriminating the respective allopatric and sympatric areas for L. italicus and L. v. meridionalis, using a selected group of predictors. No data standardization or normalization were performed for these variables.
Statistical analyses and graphics were performed using the package NCSS version 11 for Windows.
Model projection to past scenarios
SDMs obtained for the two target species under the current conditions were projected to the LGM and MOL scenarios using both the CCSM4 and the MIROC-ESM paleoclimatic reconstructions. A multivariate environmental similarity surface (MESS)  was computed in order to quantify the degree of extrapolation (i.e. prediction to environmental conditions that differ from those represented within the training data ) that the obtained SDMs would require when projected to the considered past scenarios . This analysis was performed using the ‘mess’ function implemented in the ‘dismo’ package. For each combination of species-scenario-GCM, the percent extent of the study area showing MESS values lower than −20 was calculated in Arcmap 10.0.
SDMs’ projections resulting from the two different GCMs were combined firstly by averaging them (‘Mean’ consensus approach ); then, in order to reduce the influence of possible spurious projections due to model extrapolation [58, 60, 61], a weighted average was also computed by assigning to the projection resulting from a GCM a weight being inversely proportional to the degree of “environmental novelty” of the GCM’s hindcasted climate with respect to the current conditions across the presence-background points. In this way, the more the SDM projection considering a certain GCM requires extrapolation, the less it contributes to the combined projection, so that this latter would be driven mostly by the projections having less intrinsic uncertainty.
In order to plainly formulate the above-mentioned weighted average (i.e. avoiding negative weights which could drive the values of the averaged projection outside the range 0–1), a transformation of the multivariate environmental similarity (MES)  to a positively-valued multivariate environmental dissimilarity index (MEDI) is here proposed.
Hence, MEDI values higher than 100 (corresponding to MES < 0) will indicate sites with one or more variables outside the range of values characterizing the reference points (i.e. non-analog climate ), while MEDI values lower than 100 (MES > 0) will indicate sites with environmental conditions within the range of the reference points, and increasingly common across the reference sites as MEDI tends to 0.
We propose the use of this weighted average, based on the inverse of MEDI, because it is concurrently intuitive and proper with respect to the aim of down-weighting extrapolation. Indeed, for example, in case we attributed weights based on the inverse squared MES values, as a common transformation to avoid negative weights, analogue and non-analogue environments with the same MES absolute value (but with positive and negative sign, respectively) would have the same weight.
To test for possible statistically significant differences between the environmental niches of the target species, two different methods were used. The first one is based on the comparison of the observed value of niche overlap, estimated by calculating the Schoener’s D metric  on SDMs-derived geographical projections of the suitability/probability of occurrence for the two target species, to a null distribution of overlap D values computed on simulated niches built through randomization procedures [64–66]. This approach was implemented through the ENMTools software . The second method is based on a preliminary analysis by means of an ordination technique aimed at reducing the set of input variables to two principal component axes defining a gridded environmental space inside which the kernel-smoothed density of occurrence of the two target species is calculated taking into account the density of the different combinations of environmental conditions available to each species . Also in this case, the observed overlap D value between the two species in the gridded environmental space is compared to a null distribution of D values computed between simulated niches built through randomization procedures . For this latter method the R package ‘ecospat’ , using the ‘PCA-env’  as ordination technique, was used.
For both the approaches, we first tested for niche equivalency and then extended the analysis to the less restrictive hypothesis of niche similarity , comparing the observed values of niche overlap to the 95th percentile density of the simulated values forming the null distribution [64, 65]. The random points used to perform the background test in ENMTools were generated through the Sampling Design Tool in ArcMap 10.0 .
The ecospat.plot.niche.dyn function was used to create plots of each species’ occurrence density with respect to the single input variables .
The 77.6% (388) of the 500 occurrences analyzed for Lissotriton v. meridionalis falls in the allopatric area, while for L. italicus the percentage is only of 43.1% (173/401). The sympatric area, located in Central and part of Southern Apennines, corresponds to the 16.4% of the total area of L. italicus, with 228 occurrences, and to the 14.8% of the total area of L. v. meridionalis, with 112 occurrences (Fig. 1).
Presence data showed no spatial autocorrelation for both species (Moran’s I = 0.186, p-value = 0.348, z-score = 0.939 for L. italicus and Moran’s I = 0.161, p-value = 0.426, z-score = 0.796 for L. v. meridionalis).
On the basis of the obtained correlation matrix (see Additional file 2), we selected eleven out of the nineteen Worldclim bioclimatic variables, namely BIO2, BIO4, BIO5, BIO6, BIO7, BIO8, BIO9, BIO12, BIO15, BIO17 and BIO19.
The correlation analysis between the matrices representing environmental conditions within the portions of L. v. meridionalis’ range falling inside and outside the study area resulted in r = 0.990 (p < 0.001). Thus, it can be reasonably assumed that the SDMs built for L. v. meridionalis on the considered study area are properly informed about the full set of climatic conditions relevant to the species, so that the Padano-Venetian plain and Prealps can be excluded from the study area without weakening model calibration.
The following results refer primarily to the models built through the BRT technique. However, we have also reported in an additional table and in an additional figure (see Additional file 4) the results of the Maxent analyses, which were performed in order to verify the overall consensus of predictions between different model classes . Both the techniques provided very similar responses with respect to the two target species considered.
Models of current distribution obtained through the BRT technique showed high discrimination for both species, with a cross-validated AUC = 0.892 ± 0.008 for L. italicus and a cross-validated AUC = 0.854 ± 0.011 for L. v. meridionalis. Models built through Maxent achieved similar, though lower, scores, with cross-validated AUC = 0.851 ± 0.009 for L. italicus and cross-validated AUC = 0.794 ± 0.013 for L. v. meridionalis. With respect to the calibration performance shown by the BRT models, as measured by the deviance scores, similar patterns emerged for both Lissotriton species: the model built for L. italicus showed an initial mean total deviance of 1.386, and the final model, resulting from 12,100 trees, yielded a mean residual deviance of 0.468; for L. v. meridionalis, instead, starting from the same mean total deviance of 1.386, the final model, built on 10,050 trees, yielded a mean residual deviance of 0.615.
Relative contribution of the six most influential predictors for each target species
Lissotriton vulgaris meridionalis
Relative contribution (%)
Relative contribution (%)
We built models for the current scenario, using the previously selected bioclimatic variables (except for BIO5, which was highly correlated with ALT, see Additional file 2) and the topographic ones, in order to assess the importance of the latter. These models resulted in a clear predominance of the climate-related predictors; in fact, ASPECT and SLOPE showed low contribution values for both species. With regard to ALT, even though it resulted as the second highest contributor for L. italicus (ALT = 11.2%, 1st predictor: BIO19 = 15.3%), and the third highest contributor for L. v. meridionalis (ALT = 9.3%, 1st predictor: BIO17 = 20.9%, 2nd predictor: BIO2 = 10.5%), the corresponding marginal response curves showed a similar and slightly multi-modal trend for both species, probably because of the correlation between this variable and temperature- and precipitation-related ones.
With respect to the current climatic conditions across the training points, the MESS analysis returned a considerable degree of extrapolation only for the LGM scenarios, particularly for the MIROC-ESM model (Additional file 5). Nonetheless, the presence of non-analog climate within the LGM scenarios does not seem to noticeably affect the projections of the BRT models built on the current conditions when it is accounted for by penalizing the models requiring extrapolation in our weighted average formula, as emerges from the plot in Additional file 5.
Discriminant Stepwise Analysis considering the six most contributing variables within the BRT models obtained
F to enter
Classification matrix resulting from the Discriminant Stepwise Analysis
L. italicus/L. v. meridionalis
L. v. meridionalis
L. italicus/L. v. meridionalis
L. v. meridionalis
Both the niche equivalency tests, performed through ENMTools and Ecospat, confirmed significant differences between the environmental niches of the two species. In ENMTools, the observed niche overlap value between the BRT models obtained for the two Lissotriton species under the current climatic conditions was D = 0.371, falling far outside the 95th percentile of the null distribution (D 5% = 0.885, D 95% = 0.914); in Ecospat, instead, the observed niche overlap value was D = 0.463 while the null distribution of simulated niche overlaps ranges between D 5% = 0.591 and D 95% = 0.700, leading to the rejection of the null hypothesis (p < 0.001).
The predicted suitability along the latitudinal gradient spanning the study area, for both L. v. meridionalis and L. italicus, is shown in Fig. 7b.
Environmental conditions have deeply affected the distributional patterns of many animal groups in southern European peninsulas during the Quaternary (e.g., [70–72]). For instance, climate played an important role in amphibians’ biogeography, as recently observed by Reino et al.  for some parapatric anurans.
In agreement with this hypothesis, our comprehensive analysis of the SDMs obtained for Lissotriton v. meridionalis and L. italicus clearly shows that these species occur preferentially in sites hosting specific bioclimatic conditions, with precipitation-linked variables (i.e. BIO17 and BIO19) and annual temperature variation (BIO7) strongly influencing their past and current distribution. Discriminant analysis and the outputs from the niche equivalency and niche similarity tests furtherly confirm that the environmental niches of these two newts have different peculiarities but are not strongly divergent, consistently with the presence of a sympatric area hosting conditions suitable for both the taxa. The discrepancies between the results obtained from the two niche similarity tests are probably due to the different techniques that the two methods use to select and weight the variables on which niche overlap is then computed  (SDM-based for ENMtools and PCA for Ecospat).
The opposite trends shown for the two urodeles by the marginal response curves of the most contributing variables suggest that the SDMs obtained properly distinguished the different sets of bioclimatic conditions favoring or not the presence of each Lissotriton species in allopatric or sympatric areas. Precipitation of the driest quarter (BIO17) and precipitation of the coldest quarter (BIO19) are not only the most contributing variables for L. v. meridionalis and L. italicus respectively, but they also interact each other (see Fig. 3a) and with other bioclimatic variables (see Fig. 3b and c), increasing synergistically the predicted suitability in a certain range. The first and the second interactions shown in Fig. 3a and b reflect some phenological traits of L. italicus: this species can mate and lay eggs, depending on microclimatic conditions, during the whole year . With water available even during cold months (i.e. high values of BIO19) and relatively stable annual temperatures (low values of BIO7), L. italicus can reproduce during autumn and winter. The low values of precipitation during the driest quarter (i.e. BIO17) associated with relatively high values of predicted suitability (see Fig. 3) are consistent with the reproductive phenology of this species, whose embryo and larval development takes, on average, half of the time with respect to L. v. meridionalis , so that L. italicus individuals generally reach the adult phase before the driest months. Furthermore, L. italicus easily aestivate during extremely dry periods . On the contrary, L. v. meridionalis is strongly dependent on the precipitation of driest months, because of the longer time span needed for its embryo and larval development, and also because it usually lives in shallow water bodies (up to 50–60 cm of depth), which may rapidly drain during dry periods . The coldest temperatures of the year (BIO6) weakly interact with BIO17 within the BRT model obtained for L. v. meridionalis probably because this newt, which can hibernate in water during winter, increasing the probability of a successful reproduction during spring , may not survive the harsh winter conditions which occur at higher latitudes with respect to the area inhabited by L. italicus. Contrarily to these precipitation-and-temperature-linked variables, topographic predictors were shown to be less informative with respect to the distributional patterns of the two target species.
The potential altitudinal range of the sympatric area located in Central Italy spans from the sea level to almost 3000 m a.s.l., thus including from Mediterranean to alpine habitats. However, the real altitudinal range of the two newts analyzed, based on the available data, goes from 1 to 1650 m a.s.l.. In this area and in this altitudinal interval, highly contributing variables for both the two newts show intermediate values (see Fig. 3), making possible the coexistence of L. v. meridionalis and L. italicus, coherently with a scenario of niche differentiation driven by the climate.
Particularly interesting is also the evidence of the isolated “patch” with high predicted suitability for L. italicus located in the central portion of Northern Apennines (see Fig. 2); it emerges from both the BRT and the Maxent models even though the species has never been detected there. This potentially suitable area may not have been colonized both because of the existence of a climatically unsuitable interposed area and/or because of the presence of L. v. meridionalis determining possible competitive displacement.
The projections of the BRT models to the two past scenarios were made under the assumption of temporal niche conservatism , since the temporal scale of our hindcasting is relatively short. The MESS analysis confirmed that our projections to past scenarios required some degree of extrapolation only for the LGM. Nonetheless, the transformation of the MES to the MEDI permitted us to verify that the patterns of predicted suitability show minimal changes when model extrapolation is penalized in averaging projections from different GCMs.
Model hindcasting provided some interesting insights into the influence that climatic oscillations during the Late Pleistocene and the Holocene had on the distribution of the two Lissotriton species analyzed. The diffused area with medium-to-high values of predicted suitability for L. v. meridionalis during the Last Glacial Maximum is linked to the precipitation trends in the Central Apennines. Even if it may seem counter-intuitive, since it is commonly thought that water should not have been fully available in that period due to retention within glaciers, precipitation during the LGM increased in peninsular Italy, because of a regional southerly atmospheric circulation which brought humidity from the central region of the Mediterranean Sea towards the Italian peninsula . These climatic conditions could have allowed L. v. meridionalis to move further southwards than Central Apennines, as suggested by some areas in the southern peninsula showing medium suitability values. The past presence of this species at lower latitudes with respect to its actual occurrence area is also confirmed by some fossil records in the Apulian region . On the contrary, the increasing values of the precipitation of the driest quarter within the central and south-western areas of Apennines during the LGM may have determined a shift in the distribution of L. italicus, which is predicted by our models to prefer low-to-medium values of precipitation in the driest quarter, towards the drier Adriatic region, which would have been also at lower risk of competition with L. v. meridionalis since this latter suffers from dry conditions (see above). This probably has favoured the colonization of the Apulian province, with high values of the temperature annual range hindering a northwards Adriatic colonization. This scenario agrees with what hypothesized by Canestrelli et al.  for the differentiation timing of the Apulian L. italicus clade.
During the Middle Holocene, the increase of the temperature annual range within the central and south-eastern portion of Apennines determined unfavorable conditions for L. italicus, contracting its range to the southern part of the LGM Tyrrhenian corridor, where the medium-to-high values of the precipitation of coldest quarter maintained suitable conditions. The absence of records of this species at latitudes higher than the current sympatric zone supports this scenario. Similarly, the climatic changes during the transition from the Last Glacial Maximum to the interglacial Holocene period seem to have determined the contraction and fragmentation of the potential range for L. v. meridionalis within the study area, which is now centered in the central and northern Apennines. The reason of the medium-to-low predicted suitability for both species in the MOL scenario can be found in a general decrease of precipitations in the study area [75, 77, 78], with slight increases only in the north-eastern area of Apennines.
Our results regarding L. v. meridionalis distribution during the Last Glacial Maximum also agree with the species’ evolutionary history throughout the Pleistocene inferred by Maura et al. . Their rejection of a scenario of southern refugia with postglacial northwards re-colonization is corroborated by our models. In fact, multiple and separated areas showing medium-to-highly suitable areas can be found in Central and Northern Apennines, supporting the hypothesis of multiple northern glacial refugia.
The comparison of hindcasted SDMs with results from previous studies based on genetic evidence confirmed to be a good strategy to shed light on some of the factors affecting the distributional shifts of the target species over time [32, 80, 81]. Trends of contraction and expansion emerging from our work generally agree with the scenarios of Late Pleistocene glacial refugia well supported in literature [12, 71, 82, 83], especially for the Italian peninsula (e.g., [11, 22, 35, 84]).
The strong influence of temperature annual range and precipitation-linked variables on the distribution of Lissotriton v. meridionalis and L. italicus clearly emerges from our SDMs. The effect of climate on the target species’ distribution is further supported by the results emerging from discriminant analysis and niche equivalency and similarity tests. Moreover, the coherence of the SDMs’ hindcasting with respect to the Pleistocene evolutionary scenarios inferred in previous studies for the two species analyzed suggests that the proper implementation of SDMs and the comparison of their outputs with evidences emerging from molecular-based phylogeographical analyses permit to better understand the complex interactions within the different biogeographical processes shaping the distribution of species. Nonetheless, despite the noticeable agreement between molecular phylogeographic reconstructions and patterns of past predicted suitability resulting from the SDMs obtained, it is still possible that some climatic variables individuated as less contributing in the models built for the current scenario may have indeed played a role in shaping the past distributions of the two Lissotriton species analyzed.
We are very grateful to Daniele Salvi for his precious comments and to Simone Fattorini for his constructive review of a first version of the manuscript. We are also thankful to the editor and to the two anonymous reviewers who helped us in improving the clearness and quality of the manuscript.
Mattia Iannella (PhD) is a post-doctoral researcher; his field of research is focused on amphibians and reptiles’ ecology and conservation, applying both population monitoring techniques on field and GIS and SDMs analysis.
Francesco Cerasoli is a PhD student in Environmental Sciences at the University of L’Aquila, with major interest in the application of modelling techniques and phylogeographical analyses to conservation, biogeographic and evolutionary issues.
Maurizio Biondi is a professor of Zoology at the University of L’Aquila, active in mathematical modelling applied to phylogenetic reconstructions, biogeographical research and ecological bio-characterization. He is also specializing in the systematics and evolution of Coleoptera Chrysomelidae, with a special interest in the host plant-insect relationships.
Availability of data and materials
The dataset supporting the conclusions of this article is included within the additional files. In order to avoid the risk of illegal withdrawal by poachers, who could use published data, the coordinates provided for the presence points of the target species are not at high GPS resolution.
MI, FC and MB conceived the idea; MI collected and gathered the data; MI, FC and MB analyzed the data; FC conceived the idea of the MEDI and the relative application to SDMs projections; MI, FC and MB wrote the paper. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Spellerberg IF, Sawyer JW. An introduction to applied biogeography. Cambridge: Cambridge University Press; 1999.Google Scholar
- Lomolino MV, Heaney LR. Frontiers of biogeography: new directions in the geography of nature. MA: Sinauer Associates Sunderland; 2004.Google Scholar
- Whittaker RJ, Araújo MB, Jepson P, Ladle RJ, Watson JE, Willis KJ. Conservation biogeography: assessment and prospect. Divers Distrib. 2005;11:3–23.View ArticleGoogle Scholar
- Cayuela L, Golicher D, Newton A, Kolb M, de Alburquerque F, Arets E, Alkemade J, Pérez A. Species distribution modeling in the tropics: problems, potentialities, and the role of biological data for effective species conservation. Tropical Conservation Science. 2009;2:319–52.View ArticleGoogle Scholar
- KMPMdB F, SFdB F, RCd P, Beisiegel B, Breitenmoser C. Species distribution modeling for conservation purposes. Nat Conserv. 2012;10:214–20.View ArticleGoogle Scholar
- Brown JL, Cameron A, Yoder AD, Vences MA. Necessarily complex model to explain the biogeography of the amphibians and reptiles of Madagascar. Nat Commun. 2014;5:5046.View ArticlePubMedGoogle Scholar
- Gaston KJ. Geographic range limits: achieving synthesis. Proc R Soc Lond B Biol Sci. 2009;276:1395–406.View ArticleGoogle Scholar
- Stewart JR, Lister AM, Barnes I, Dalén L. Refugia revisited: individualistic responses of species in space and time. Proc R Soc Lond B Biol Sci. 2010;277:661–71.View ArticleGoogle Scholar
- Pelletier TA, Carstens BC. Comparing range evolution in two western Plethodon salamanders: glacial refugia, competition, ecological niches, and spatial sorting. J Biogeogr. 2016;43:2237–49.View ArticleGoogle Scholar
- Baselga A, Lobo JM, Svenning JC, Araujo MB. Global patterns in the shape of species geographical ranges reveal range determinants. J Biogeogr. 2012;39:760–71.View ArticleGoogle Scholar
- Canestrelli D, Aloise G, Cecchetti S, Nascetti G. Birth of a hotspot of intraspecific genetic diversity: notes from the underground. Mol Ecol. 2010;19:5432–51.View ArticlePubMedGoogle Scholar
- Reino L, Ferreira M, Martínez-Solano Í, Segurado P, Xu C, Márcia Barbosa A. Favourable areas for co-occurrence of parapatric species: niche conservatism and niche divergence in Iberian tree frogs and midwife toads. J Biogeogr. 2017;44:88–98.View ArticleGoogle Scholar
- Richards CL, Carstens BC, Lacey Knowles L. Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. J Biogeogr. 2007;34:1833–45.View ArticleGoogle Scholar
- Vega R, Fløjgaard C, Lira-Noriega A, Nakazawa Y, Svenning JC, Searle JB. Northern glacial refugia for the pygmy shrew Sorex Minutus in Europe revealed by phylogeographic analyses and species distribution modelling. Ecography. 2010;33:260–71.Google Scholar
- Schorr G, Holstein N, Pearman P, Guisan A, Kadereit J. Integrating species distribution models (SDMs) and phylogeography for two species of alpine Primula. Ecology and evolution. 2012;2:1260–77.View ArticlePubMedPubMed CentralGoogle Scholar
- Soberon J, Peterson AT. Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodiv Inform. 2005;2:1–10.View ArticleGoogle Scholar
- Godsoe WI. can't define the niche but I know it when I see it: a formal link between statistical theory and the ecological niche. Oikos. 2010;119:53–60.View ArticleGoogle Scholar
- Heads M. The relationship between biogeography and ecology: envelopes, models, predictions. Biol J Linn Soc. 2015;115:456–68.View ArticleGoogle Scholar
- Hofreiter M, Stewart J. Ecological change, range fluctuations and population dynamics during the Pleistocene. Curr Biol. 2009;19:R584–94.View ArticlePubMedGoogle Scholar
- Araújo MB, Nogués-Bravo D, Diniz-Filho JAF, Haywood AM, Valdes PJ, Rahbek C. Quaternary climate changes explain diversity among reptiles and amphibians. Ecography. 2008;31:8–15.View ArticleGoogle Scholar
- Nogués-Bravo D, Ohlemüller R, Batra P, Araújo MB. Climate predictors of late quaternary extinctions. Evolution. 2010;64:2442–9.PubMedGoogle Scholar
- La Greca M. Vicende paleogeografiche e componenti della fauna italiana. La fauna in Italia Touring Editore, Milano e Ministero dell’Ambiente e della Tutela del Territorio, Roma 2002.Google Scholar
- Razzetti E, Lapini L, Bernini F: Anfibi d'Italia. In Quad Conserv Nat. Volume 29. Edited by Lanza B, Nistri A, Vanni S: Ministero dell'Ambiente e della Tutela del Territorio e del Mare; Istituto Superiore per la protezione la ricerca ambientale; 2009: 456 + 451 CD bilingue.[Editori GG (Series Editor).Google Scholar
- Scillitani G, Tripepi S: Anfibi d'Italia. In Quad Conserv Nat. Volume 29. Edited by Lanza B, Nistri A, Vanni S: Ministero dell'Ambiente e della Tutela del Territorio e del Mare; Istituto Superiore per la protezione la ricerca ambientale; 2009: 456 + 451 CD bilingue.[Editori GG (Series Editor).Google Scholar
- Mancino G. Anomalie meiotiche nella spermatogenesi dell'ibrido Triturus italicus♀ x Triturus vulgaris♂. Italian Journal of Zoology. 1961;28:691–701.Google Scholar
- Macgregor HC, Sessions SK, Arntzen J. An integrative analysis of phylogenetic relationships among newts of the genus Triturus (family Salamandridae), using comparative biochemistry, cytogenetics and reproductive interactions. J Evol Biol. 1990;3:329–73.View ArticleGoogle Scholar
- Ficetola GF, Maiorano L. Contrasting effects of temperature and precipitation change on amphibian phenology, abundance and performance. Oecologia. 2016;181:683–93.View ArticlePubMedGoogle Scholar
- Coyne JA, Orr HA. Speciation. Sunderland, MA: Sinauer Associates, Inc; 2004.Google Scholar
- Harris DJ. Inferring species interactions from co-occurrence data with Markov networks. Ecology. 2016;97:3308–14.View ArticlePubMedGoogle Scholar
- Cazelles K, Araújo MB, Mouquet N, Gravel DA. Theory for species co-occurrence in interaction networks. Theor Ecol. 2016;9:39–48.View ArticleGoogle Scholar
- Pauls SU, Theissinger K, Ujvarosi L, Balint M, Haase P. Patterns of population structure in two closely related, partially sympatric caddisflies in Eastern Europe: historic introgression, limited dispersal, and cryptic diversity1. J North Am Benthol Soc. 2009;28:517–36.View ArticleGoogle Scholar
- Svenning J-C, Fløjgaard C, Marske KA, Nógues-Bravo D, Normand S. Applications of species distribution modeling to paleobiology. Quat Sci Rev. 2011;30:2930–47.View ArticleGoogle Scholar
- Lanza B, Nistri A, Vanni S: Anfibi d'Italia. Ministero dell'Ambiente e della Tutela del Territorio e del Mare; Istituto Superiore per la protezione la ricerca ambientale; 2009.Google Scholar
- Iannella M: Central Apennines batrachofauna: status of knowledge, chronogeonemy and gap analysis, aiming to an informed wildlife management. University of L'Aquila, clinical medicine, Public Health, Life and Environment Sciences Dept - Environmental Sciences Sect 2015.Google Scholar
- Canestrelli D, Sacco F, Nascetti G. On glacial refugia, genetic diversity, and microevolutionary processes: deep phylogeographical structure in the endemic newt Lissotriton Italicus. Biol J Linn Soc. 2012;105:42–55.View ArticleGoogle Scholar
- Canestrelli D, Salvi D, Maura M, Bologna MA, Nascetti G. One species, three Pleistocene evolutionary histories: phylogeography of the Italian crested newt, Triturus Carnifex. PLoS One. 2012;7:e41754.View ArticlePubMedPubMed CentralGoogle Scholar
- Pabijan M, Zieliński P, Dudek K, Chloupek M, Sotiropoulos K, Liana M, Babik W. The dissection of a Pleistocene refugium: phylogeography of the smooth newt, Lissotriton Vulgaris, in the Balkans. J Biogeogr. 2015;42:671–83.View ArticleGoogle Scholar
- Elith J, Graham HC, Anderson PR, Dudík M, Ferrier S, Guisan A, Hijmans JR, Huettmann F, Leathwick RJ, Lehmann A, et al. novel methods improve prediction of species’ distributions from occurrence data. Ecography. 2006;29:129–51.View ArticleGoogle Scholar
- Elith J, Graham CH. Do they? How do they? WHY do they differ? On finding reasons for differing performances of species distribution models. Ecography. 2009;32:66–77.View ArticleGoogle Scholar
- Merow C, Smith MJ, Silander JAA. Practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography. 2013;36:1058–69.View ArticleGoogle Scholar
- Elith J, Leathwick JR, Hastie TA. Working guide to boosted regression trees. J Anim Ecol. 2008;77:802–13.View ArticlePubMedGoogle Scholar
- Ridgeway G. Generalized boosted regression models. The Comprehensive R Archive Network. 2015. https://CRAN.R-project.org/package=gbm.
- Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190:231–59.View ArticleGoogle Scholar
- Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, Marquéz JRG, Gruber B, Lafourcade B, Leitão PJ. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography. 2013;36:27–46.View ArticleGoogle Scholar
- Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.View ArticleGoogle Scholar
- ESRI. ArcGIS Desktop: Release 10. Redlands: Environmental Systems Research Institute; 2010.Google Scholar
- Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, Lawrence DM, Neale RB, Rasch PJ, Vertenstein M. The community climate system model version 4. J Clim. 2011;24:4973–91.View ArticleGoogle Scholar
- Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, Nozawa T, Kawase H, Abe M, Yokohata T. MIROC-ESM 2010: model description and basic results of CMIP5-20c3m experiments. Geosci Model Dev. 2011;4:845.View ArticleGoogle Scholar
- Chefaoui RM, Lobo JM. Assessing the effects of pseudo-absences on predictive distribution model performance. Ecol Model. 2008;210:478–86.View ArticleGoogle Scholar
- Senay SD, Worner SP, Ikeda T. Novel three-step pseudo-absence selection technique for improved species distribution modelling. PLoS One. 2013;8:e71218.View ArticlePubMedPubMed CentralGoogle Scholar
- Booth TH, Nix HA, Busby JR, Hutchinson MF. BIOCLIM: The first species distribution modelling package, its early applications and relevance to most current MAXENT studies. Divers Distrib. 2014;20:1–9.Google Scholar
- Hijmans RJ, Phillips S, Leathwick J, Elith J. dismo: Species Distribution Modeling. The Comprehensive R Archive Network. 2016. https://cran.r-project.org/package=dismo.
- Jiménez-Valverde A, Lobo JM, Hortal J. Not as good as they seem: the importance of concepts in species distribution modelling. Divers Distrib. 2008;14:885–90.View ArticleGoogle Scholar
- Barbet-Massin M, Jiguet F, Albert CH, Thuiller W. Selecting pseudo-absences for species distribution models: how, where and how many? Methods Ecol Evol. 2012;3:327–38.View ArticleGoogle Scholar
- Elith J, Leathwick JR. Species distribution models: ecological explanation and prediction across space and time. Annu Rev Ecol Evol Syst. 2009;40:677–97.View ArticleGoogle Scholar
- Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology. 1982;143:29–36.View ArticlePubMedGoogle Scholar
- Tabachnick BG, Fidell LS. Using multivariate statistics. 4th ed. Needham Heights: Allyn and Bacon; 2001.Google Scholar
- Elith J, Kearney M, Phillips S. The art of modelling range-shifting species. Methods Ecol Evol. 2010;1:330–42.View ArticleGoogle Scholar
- Nogués-Bravo D. Predicting the past distribution of species climatic niches. Glob Ecol Biogeogr. 2009;18:521–31.View ArticleGoogle Scholar
- Marmion M, Parviainen M, Luoto M, Heikkinen RK, Thuiller W. Evaluation of consensus methods in predictive species distribution modelling. Divers Distrib. 2009;15:59–69.View ArticleGoogle Scholar
- Zurell D, Elith J, Schröder B. Predicting to new environments: tools for visualizing model behaviour and impacts on mapped distributions. Divers Distrib. 2012;18:628–34.View ArticleGoogle Scholar
- Fitzpatrick MC, Hargrove WW. The projection of species distribution models and the problem of non-analog climate. Biodivers Conserv. 2009;18:2255–61.View ArticleGoogle Scholar
- Schoener TW. The Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology. 1968;49:704–26.View ArticleGoogle Scholar
- Broennimann O, Fitzpatrick MC, Pearman PB, Petitpierre B, Pellissier L, Yoccoz NG, Thuiller W, Fortin M-J, Randin C, Zimmermann NE, et al. Measuring ecological niche overlap from occurrence and spatial environmental data. Glob Ecol Biogeogr. 2012;21:481–97.View ArticleGoogle Scholar
- Warren DL, Glor RE, Turelli M. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution. 2008;62:2868–83.View ArticlePubMedGoogle Scholar
- McCormack JE, Zellmer AJ, Knowles LL. Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation?: insights from tests with niche models. Evolution. 2010;64:1231–44.PubMedGoogle Scholar
- Warren DL, Glor RE, Turelli M. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography. 2010;33:607–11.View ArticleGoogle Scholar
- Di Cola V, Broennimann O, Petitpierre B, Breiner FT, D'Amen M, Randin C, Engler R, Pottier J, Pio D, Dubuis A, et al. ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography. 2017;40:774–87.View ArticleGoogle Scholar
- Buja K, Menza C. Sampling design tool for ArcGIS: instruction manual [for ESRI ArcGIS 10.0 service pack 3 or higher]. (NOAA ed. SIlver Spring, MD; 2013.Google Scholar
- Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405:907–13.View ArticlePubMedGoogle Scholar
- Hewitt GM. Quaternary phylogeography: the roots of hybrid zones. Genetica. 2011;139:617–38.View ArticlePubMedGoogle Scholar
- Schmitt T, Varga Z. Extra-Mediterranean refugia: the rule and not the exception? Front Zool. 2012;9:22.View ArticlePubMedPubMed CentralGoogle Scholar
- Zurell D, Berger U, Cabral JS, Jeltsch F, Meynard CN, Münkemüller T, Nehrbass N, Pagel J, Reineking B, Schröder B. The virtual ecologist approach: simulating data and observers. Oikos. 2010;119:622–35.View ArticleGoogle Scholar
- Wiens JJ, Ackerly DD, Allen AP, Anacker BL, Buckley LB, Cornell HV, Damschen EI, Jonathan Davies T, Grytnes JA, Harrison SP. Niche conservatism as an emerging principle in ecology and conservation biology. Ecol Lett. 2010;13:1310–24.View ArticlePubMedGoogle Scholar
- Giraudi C. Middle Pleistocene to Holocene glaciations in the Italian Apennines. In: Ehlers J, Gibbard PL, editors. Quaternary glaciations-extent and chronology: part I: Europe, vol. 2. Amsterdam: Elsevier; 2011. p. 211–9.View ArticleGoogle Scholar
- Database of lower vertebrate: fossil fishes, amphibian, reptiles (fosFARbase). [http://www.wahre-staerke.com/; Archived by WebCite® at http://www.webcitation.org/6eNe3PsuJ].
- Davis BA, Brewer S, Stevenson AC, Guiot J. The temperature of Europe during the Holocene reconstructed from pollen data. Quat Sci Rev. 2003;22:1701–16.View ArticleGoogle Scholar
- Combourieu-Nebout N, Peyron O, Bout-Roumazeilles V, Goring S, Dormoy I, Joannin S, Sadori L, Siani G, Magny M. Holocene vegetation and climate changes in the central Mediterranean inferred from a high-resolution marine pollen record (Adriatic Sea). Clim Past. 2013;9:2023.View ArticleGoogle Scholar
- Maura M, Salvi D, Bologna MA, Nascetti G, Canestrelli D. Northern richness and cryptic refugia: Phylogeography of the Italian smooth newt Lissotriton Vulgaris Meridionalis. Biol J Linn Soc. 2014;113:590–603.View ArticleGoogle Scholar
- Waltari E, Hijmans RJ, Peterson AT, Nyári ÁS, Perkins SL, Guralnick RP. Locating Pleistocene refugia: comparing phylogeographic and ecological niche model predictions. PLoS One. 2007;2:e563.View ArticlePubMedPubMed CentralGoogle Scholar
- Wielstra B, Crnobrnja-Isailović J, Litvinchuk SN, Reijnen BT, Skidmore AK, Sotiropoulos K, Toxopeus AG, Tzankov N, Vukov T, Arntzen JW. Tracing glacial refugia of Triturus newts based on mitochondrial DNA phylogeography and species distribution modeling. Front Zool. 2013;10:13.View ArticlePubMedPubMed CentralGoogle Scholar
- Gonçalves H, Martìnez-Solano I, Pereira R, Carvalho B, Garcìa-Parìs M, Ferrand N. High levels of population subdivision in a morphologically conserved Mediterranean toad (Alytes Cisternasii) result from recent, multiple refugia: evidence from mtDNA, microsatellites and nuclear genealogies. Mol Ecol. 2009;18:5143–60.View ArticlePubMedGoogle Scholar
- Previšić A, Walton C, Kučinić M, Mitrikeski PT, Kerovec M. Pleistocene divergence of Dinaric Drusus endemics (Trichoptera, Limnephilidae) in multiple microrefugia within the Balkan peninsula. Mol Ecol. 2009;18:634–47.View ArticlePubMedGoogle Scholar
- Salvi D, Harris DJ, Kaliontzopoulou A, Carretero MA, Pinho C. Persistence across Pleistocene ice ages in Mediterranean and extra-Mediterranean refugia: phylogeographic insights from the common wall lizard. BMC Evol Biol. 2013;13:147.View ArticlePubMedPubMed CentralGoogle Scholar
- Minelli A, Ruffo S, Vigna TA. Le province faunistiche italiane. Checklist e distribuzione della fauna italiana Memorie del Museo Civico di Storia Naturale di Verona, 2a serie. Sezione Scienze della Vita. 2006;16:37–9.Google Scholar
- Biondi M, Urbani F, D’Alessandro P. Endemism patterns in the Italian leaf beetle fauna (Coleoptera, Chrysomelidae). ZooKeys. 2013;332:177–205.View ArticleGoogle Scholar