Improving biological relevance of model projections in response to climate change by considering dispersal amongst lineages in an amphibian

When modelling future or past geographic distributions of a species, attention should be paid to the possible differentiated responses to climate changes between lineages. Dispersal also plays an important role in the capacity of species to track suitable climate, which is particularly relevant for amphibians with limited dispersal ability. In this study, we included different lineages and dispersal distances into species distribution models to make them more biologically relevant in face of climate change scenarios.


| INTRODUC TI ON
Climate change caused by human activities is one of the main drivers of biodiversity decline worldwide (Parmesan, 2006). One of the main effects of climate change is the shift of thermal isoclines upward and poleward, leading to the redistribution of climate conditions on a global scale (Loarie et al., 2009). Species may respond to those changes through phenotypic plasticity and local adaptation (Urban et al., 2014) or may shift their distribution range to track favourable climate conditions (Lenoir et al., 2020). Possibilities and magnitudes of range shifts are determined by the species' propensity and ability to disperse (Travis et al., 2013). Climate-induced dispersal has been reported in many species (Blaustein et al., 2010) and likely contributes to the unprecedented biodiversity redistribution reported in recent studies (Lenoir et al., 2020). However, the evolvability of thermal performances as well as dispersal potential strongly varies at both interspecific and intraspecific levels (Bestion et al., 2015), therefore making accurate prediction of species response to climate change a challenging task.
Niche-based models (also called 'ecological niche models' or 'species distribution models'), relying on correlative relationships between observed species distributions and environmental factors, are powerful tools to predict range shifts driven by future or past climate variations (Araújo et al., 2006). However, those models have historically suffered from conceptual and technical limitations that critically affected the reliability of their predictions (Hernandez et al., 2006). First, relatively few studies explicitly considered species dispersal ability when modelling past and future range shifts, despite the fact that running models with unlimited dispersal overestimates the potential future range, especially for species with moderate to low dispersal ability (Urban et al., 2016). Second, most studies did not consider the intraspecific levels in climate-related phenotypic performances.
Recent studies suggest that populations of the same species may differ in their responses to climate variation (see, for instance, Bestion et al., 2015 andTrochet et al., 2018). Differences in population sensitivity to climate fluctuations likely rely on population-specific reaction norms to thermal and hydric conditions and/or adaptation to local climatic conditions (Gienapp et al., 2007;Urban et al., 2014). Those population-specific responses are expected to be particularly broad in ectotherms whose physiology and fitness components (i.e., survival, reproduction and growth) strongly depend on local/regional thermic and hygrometric variation (Beebee, 1995;Blaustein et al., 2010).
Such variation in population responses to climatic conditions could have a genetic base (Wilson, 2001) depending on population evolutionary history. This is especially true for ectotherm vertebrates from temperate areas that experienced broad climate variation and complex demographic/evolutionary processes during the Pleistocene period (Hewitt, 2000;Zeisset & Beebee, 2008). In many of those species, the occupancy of multiple glacial refuges led to the emergence of glacial lineages that have experienced repeated periods of allopatry, often punctuated by secondary contacts with more or less intense gene flow (Rowe et al., 2006). Those lineages often display partially overlapping climatic niches and contrasted thermal preferences resulting from the effect of climate-driven selection on their respective glacial refugia (Gilbert & Miles, 2017). Therefore, considering dispersal ability and bioclimatic preferences at the intraspecific or lineage level in niche-based models seems a promising way to produce more accurate projections of species distribution in response to past or future climate changes and to guide conservation actions (Fong G. et al., 2015). In this study, we focus on a western Palearctic amphibian -the yellow-bellied toad -B. variegata (sensu lato), which appears to be an excellent candidate on which to test this assumption since it has a relatively limited dispersal ability (most of its dispersal movements have distances of less than 500 m; see Cayuela, Bonnaire, et al., 2018;Cayuela, Rougemont, et al., 2018) and can be separated into different lineages. Indeed, as a result of B. variegata's complex evolutionary history during Pleistocene glaciations, four distinct lineages occupy regions with contrasted climatic conditions (Hofman et al., 2007;Pabijan et al., 2013). Moreover, rainwater fluctuations are known to impact the survival of Bombina variegata and because of its dependency to temporary ponds for its reproduction, extreme drought events increased by global warming are expected to particularly threat the development of this toad (Cayuela et al., 2014). Therefore, we asked the following nested questions: how meaningful is it to run the models while considering intraspecific levels? How has each lineage responded, and how will it respond, to climate change when we consider the dispersal ability of the species? What lineage could suffer the most in the face of future climate change?
First, we fitted correlative niche models for each B. variegata lineage using MaxEnt  and projected both past (LGM and mid-Holocene) and future climate change scenarios.
Second, we compared model projections that had limited and unlimited dispersal abilities. available, in modelling process to catch local variations and obtain more accurate results.

K E Y W O R D S
Bombina variegata, climate refugia, ecological niche modelling, Europe, intraspecific levels, niche differentiation, range-shift 2 | MATERIAL S AND ME THODS

| Model species
Bombina variegata is a pond-breeding amphibian with a widespread distribution in Europe. Over the last century, the species suffered from a sharp decline, especially in Western Europe (Lescure et al., 2011), and is registered in Appendix II of the Bern Convention and in Appendices II and IV of the EU Habitat Directive. This species occupies various habitats, such as forests and meadows, and reproduces mainly in temporary waterbodies filled by rainfall that could result from human activities (Cayuela et al., 2014;Scheele et al., 2014).
Temperature and precipitation are important determinants of postmetamorphic survival and reproduction (Cayuela et al., 2016).
Bombina variegata has a complex biogeography: distinct glacial lineages were previously identified based on mitochondrial DNA (Hofman et al., 2007), including B. v. variegata in western Europe, B. v. variegata in the Carpathian Mountains, B. v. pachypus in Italy and B. v. scabra in the Balkan Peninsula.

| Occurrence data
Occurrence data were obtained from several sources for the period 1979-2016. We gathered data of the different Bombina lineages across Europe, including southern, western and part of eastern Europe (Fijarczyk et al., 2011;Hofman et al., 2007).
Occurrence data of B. v. variegata 'Carpathians', B. v. scabra and B. v. pachypus were all obtained from the Global Biodiversity Information Facility (GBIF, http://www.gbif.org/). The data set for the western lineage of B. v. variegata, which has the largest distribution area, was completed by occurrence data from the association Faune-France, the Écrins National Park, Montagne de Reims Regional Natural Park, naturalistic associations in the region of Centre-Val de Loire, the Swiss Wildlife Mapping Center and some additional data from the northeast of France (Cayuela et al., 2020). Occurrence data were upscaled to correspond at the centroid of a 10 × 10 km grid. This transformation allows us to reduce sample bias by removing over-sampled locations and thus increases the predictive abilities of the model by reducing spatial autocorrelation . Because B. v. scabra and the two lineages of B. v. variegata were not distinguished in the GBIF and International Union for Conservation of Nature (IUCN), we separated them according to the genetic clusters identified by Fijarczyk et al. (2011) andHofman et al. (2007).
After removing suspicious occurrences (i.e., older occurrences and data outside the current species range), the sample size used for modelling was as follows: 1956 occurrence data for B. v. variegata 'western lineage' (V), 123 for B. v. variegata 'Carpathians' (C), 26 for B. v. scabra (S) and 28 for B. v. pachypus (P). These final two lineages occupy the smallest distribution areas, restricted to the Balkan region and to the main part of Italy, whereas the first has the widest distribution area (see Supporting Information Appendix S1).

| Climate data
We selected nine bioclimatic variables (representing temperature and precipitation conditions) known to influence the different development stages of the Bombina toads (see Supporting Information Appendix S2 for a complete description of variables), from egg to adulthood (Cayuela et al., 2014;Kaplan & Phillips, 2006;Reyer & Barandun, 1997). These variables were obtained from 'Climatologies at high resolution for the Earth's land surface area' (CHELSA, (Karger et al., 2017), http://www.chels a-clima te.org/) for the period 1979-2013 (30 arcsec resolution, ~1 km). Climatic layers were averaged at a 10 × 10 km spatial resolution to match the species occurrence grid.
Climate data were extracted at the location of lineage occurrences and used to fit the niche-based model.
We calculated the Pearson correlation coefficient among each of the bioclimatic variables, as well as their variance inflation factor (VIF) in the niche-based model, using the R packages 'usdm' and 'ecospat' (Broennimann et al., 2012;Naimi et al., 2014). Only variables with a correlation coefficient lower than 0.75 and a VIF of less than 10 were retained (Naimi et al., 2014). This method allowed us to avoid collinearity that could bias forecasts made by the niche-based model. Thus, we retained seven variables that were used throughout this study (see Appendix S2).
For past climatic conditions, we used CCSM4 and MIROC-ESM models according to the reviews of Svenning et al. (2011) andVarela et al. (2011). Indeed, these climate models are the most widely used in Europe and show less bias compared to other global circulation models (GCMs) for the LGM (Harrison et al., 2014). Both model outputs were extracted from the CHELSA database (Karger et al., 2017) for the LGM (i.e., 21,000 years ago) and from the WorldClim database (1.4; Fick & Hijmans, 2017) for the mid-Holocene (i.e., 6,000 years ago).
We used two other GCMs for future climate conditions expected during the 21st century: CSIRO-Mk3.6.0 and HadGEM2-ES (Collier et al., 2011;Collins et al., 2008), which we extracted from the CHELSA and WorldClim databases, respectively. These widely used GCMs were selected because they generate robust predictions (Buisson et al., 2010;Farzaneh et al., 2012) and have been recognized for generating divergent predictions (Mitchell et al., 2004). These two GCMs were used to model two future periods, 2041-2060 and 2061-2080, under two Representative Concentration Pathway (RCP) scenarios: RCP2.6 and RCP8.5, which are optimistic and pessimistic scenarios that, respectively, predict a global increase in mean temperature of 0.9-2.3°C to 3.2-5.4°C and changes in global mean precipitation correlated with changes in global mean temperature (Baek et al., 2013;Riahi et al., 2011;van Vuuren et al., 2011). Species redistribution forecasts on several GCMs and climate change scenarios account for future climate condition uncertainties in our work (Buisson et al., 2010).

| Niche similarity between lineages
Prior to modelling, we analysed the bioclimatic niches of the different lineages identified by Fijarczyk et al. (2011) and their possible overlap using the R package 'ecospat'. A similarity test was performed to estimate the climate niche similarity between two lineages by comparing climate conditions (defined by the first two components of a principal component analysis [PCA]) where lineages occur (with replications n = 1,000). We calculated the similarity measure I based on the Hellinger distance to assess niche overlaps (Schoener, 1968;Warren et al., 2008). The I metric ranges from 0 (no overlap) to 1 (perfect overlap), has no biological assumptions concerning the distribution of the two tested species and only considers them as probability distributions (Warren et al., 2008).

| General modelling
We used MaxEnt 3.4.1, one of the most used and robust models (Elith et al., 2011), which also demonstrated good performance in situations with few occurrence data (Elith et al., 2006;Hernandez et al., 2006).
Our models were calibrated according to current bioclimatic conditions from all of Europe (see above), which is the best way to encompass a large range of conditions, limit model extrapolation and identify the overall global niche of the species (Pearson et al., 2002).
We ran 50 replicates with 25% of test data (random seed, bootstrap) and set the other parameters by defaults. Because MaxEnt uses an exponential model and extrapolates projections when environmental data used are outside their present range, we used the clamping method. This method considers variables outside their range as if they were at the limit of their training range and thereby limits extrapolation (Phillips, 2017).
The model performance was assessed with the area under the receiver operating characteristic curve (AUC), the true skill statistic (TSS; usdm package) and the kappa statistic (Araújo et al., 2019). The AUC ranges between 0 and 1, where 0.5 and below indicate that the model has no predictive ability better than random (Hanley & McNeil, 1982). The kappa statistic and TSS range from −1 to 1, where 1 indicates a perfect model and values of less than 0 indicate a performance no better than random predictions (Allouche et al., 2006;Cohen, 1960). We explored models' extrapolations by identifying areas where predictors are outside of their training range with the multivariate environmental similarity surface (MESS) option.
After binarization based on the 'maximum training sensitivity plus specificity' threshold rule (Liu et al., 2013), we calculated the percentage of loss and gain of suitable habitats and projected these results on maps. Changes in areas of suitable habitats of each lineage were combined in one map by model and period.

| Species and intraspecific level
We ran separate models following the methods described above for each lineage and also without considering the intraspecific level (B. variegata sensu lato). Moreover, we ran the models with the most adapted value of regularization multiplier and feature class, according to the AICc, calculated under the R package "ENMeval" (Muscarella et al., 2014). The best regularization multipliers and feature classes were 0.5HP (Hinge plus Product features), 1T (Threshold features), 0.5LQ (Linear plus Quadratic features), 0.5HP (Hinge plus Product features) and 0.5HPQ (Hinge plus Product plus Quadratic features), respectively, for B. v. variegata 'Carpathian', B. v. pachypus, B. v. scabra, B. v. variegata 'western lineage' and B. variegata (i.e., species level).

| Limited versus unlimited dispersal scenarios
We considered the dispersal capacity based on the study from Cayuela, Bonnaire, et al. (2018), who assessed the dispersal abilities of the yellow-bellied toad at different life stages in a forest environment located in eastern France. The median dispersal event distance comprised between 444 and 695 m per year, depending on the age class (Cayuela, Bonnaire, et al., 2018). Based on these results, we considered in this study that B. variegata was able to disperse at a rate of 500 m per year, which corresponds to the distance covered by most of the individuals within a population in a favourable environment.
We took into account the dispersal ability by adding a buffer zone around the current distribution area (obtained from the IUCN) of the considered lineage. The buffer width was calculated by multiplying the number of years that separates the current period from the past or future studied period by 500 m. For each lineage at the species level, future projections were conducted in the current distribution areas with additional buffers of 15 and 30 km for the periods 2041-2060 and 2061-2080, respectively. We also ran models under unlimited dispersal for these periods.
Considering the time scale for past projections, the buffer areas were as large as Europe (equivalent to unlimited dispersal).
Nevertheless, for LGM projections, the continental area has been modified by taking into account the large continental glaciers and changes in emerged land surface due to sea levels' temporal variations (Ehlers & Gibbard, 2007;Lambeck, 2001). We considered continental glaciers to be unsuitable areas. Also, we considered islands to be unreachable and therefore unlikely to be colonized.

| Divergence of bioclimatic niches among the Bombina glacial lineages
The first two components of the PCA explained 69% of the total variance ( Figure 1). The first component represented a water balance gradient ranging from humid/fresh to dry/hot conditions. The second component described specific climate conditions in summer, distinguishing cold/dry and hot/humid locations.
We detected a significant niche similarity between the Italian (B.

| Projections under past climates
Due to the large range of the time period (6,000-21,000), past projections were made without any dispersal limitation (see Materials and Methods).
For models conducted at the species level, the location and percentage of suitable areas were similar to the western lineage (−93.1% to −88.2% for the LGM and −45.1% to −37.8% for the mid-Holocene).
Considering the western lineage of B. v. variegata, only the southwest of France was suitable 21,000 years ago (Figure 3; Table 1).
CCSM4-based projections indicated that the west coast of France, a part of the Apennine Mountains in Italy and a small area of the Balkan Peninsula were also suitable, while MIROC-ESM-based projections showed additional suitable areas in the north of Italy and in the Dinaric Alps (see Appendix S6). As a result, suitable areas were smaller than they are today, with an average of −91% to −68.2%.
During the mid-Holocene, suitable areas appeared towards the north, but these areas were smaller compared to current areas, ranging from −43% to −31.3%. For the Carpathian lineage, models' outputs showed much smaller suitable areas than the current situation, between −100% to −84.2% smaller on average. Those suitable areas were potentially located in the Pannonian Basin and in Germany 21,000 years ago (see Appendix S1). Later, 6,000 years ago, suitable areas were less restricted, with a mean percentage ranging from −80% to −20% in comparison to the current distribution. Therefore, the model showed that this lineage shifted its range towards the centre and the south of the Carpathians (Figure 3; Table 1).
The B. v. pachypus lineage model projections for 21,000 years ago identified larger suitable areas than the lineage's current range, mostly in the Balkan Peninsula, western France and Spain. Models made with the CCSM4 climate scenario resulted in suitable areas 50% larger than the current distribution ( Figure 3; Table 1; see Appendix S6). During the mid-Holocene, models showed a reduction in suitable areas restricted to the south of Italy, resulting in a climatic suitable area of −44.1% to −26.5% compared to current areas. suitable areas were about 12.3%-253% larger than current areas.
During the mid-Holocene, suitable areas were generally smaller than their current distribution by −93.8% to 33.1% (Table 1). Models identified suitable areas mainly in high elevation areas such as the Rhodope but also in Ukraine (see Appendix S6).

| Projections under future climate change
Results correspond to models run at species or lineage level by considering either unlimited or limited dispersal.

| Species level
According to species distribution models (SDMs) made at the spe-    (Figures 4 and 5). In the same way, the western lineage of B. v. variegata will also suffer from a drastic reduction in suitable habitat, and the remaining suitable areas are located in the north or altitude ( Figure 5).

F I G U R E 3
Predicted changes in habitat suitability for Bombina variegata lineages in Europe during 21,000 and 6,000 bp in comparison to the current period. Red indicates areas that were previously unsuitable but are currently suitable, dark blue indicates areas that have been suitable and are still currently suitable and light blue indicates areas that were previously suitable but are now unsuitable. Grey indicates areas that remain unsuitable, taking into account dispersal ability. Emerged land surface due to lower sea levels and major glaciers has been taken into account in the 21,000 bp models. The results of all the models are supported by good performances of the metrics that are broadly used in modelling studies.
In addition, as expected, MESS analysis indicated that reachable areas have a low degree of extrapolation, which is due to our calibration with current bioclimatic data from all of Europe (Appendix S10).
Our study revealed that considering dispersal limitation by expanding the distribution area has a broad influence on model projections of future and past distribution of B. variegata. Indeed, whereas suitable areas will endure in Europe over time, many of them will not be reachable for the yellow-bellied toad.

| Intraspecific level strongly impacts forecasts
The analysis of the bioclimatic niches of the lineages commonly identified in the scientific literature showed enough discrepancies to treat each lineage independently in the models, also suggesting contrasted thermal and hydric preferences between lineages. The integration of a coherent dispersal distance is still a major issue in modelling studies because of the lack of data, but studies that rely on demographic or genetic approaches could help (Cayuela, Rougemont, et al., 2018). Indeed, dispersal distance can be inferred though genetic differentiation between two distinguished populations. Genetic differentiation depends on genetic drift and gene flow, and thus individuals that colonize another distant population could modify its allele frequencies and therefore its genetic differentiation (Broquet & Petit, 2009). However, because genetic differentiation also depends on genetic drift, its use as a dispersal proxy could lead to misleading results. Thus, we recommend jointly the use of genetic and demographic approaches (Cayuela, Rougemont, et al., 2018).
For instance, Cayuela et al. (2020) used a combined demographic and genetic approach to study and compare two populations of B.
v. variegata, one of which exhibited dispersal syndromes. Using this method, they identified that the costs and benefits of dispersal drive the behaviour of spatially structured populations in logging environments, and they were able to distinguish the effect of genetic drift from gene flow. These results bring us to a new point: populations could differ in terms of dispersal behaviour, and those with lower dispersal abilities and weaker genetic structure are expected to be particularly vulnerable to environmental changes (Cayuela et al., 2020).
The integration of the intraspecific differences in dispersal abilities using dispersal kernels would therefore pose a future challenge in obtaining more accurate and realistic SDMs. Nonetheless, integration of intraspecific differences would require precise knowledge of the dispersal behaviour of populations across the range extent of the species studied.
In our study, we considered a dispersal distance of 500 m per year. Even though this distance is biologically relevant, it was assessed in a homogenous landscape, while human infrastructures are known to limit dispersal events (Cayuela, Bonnaire, et al., 2018;Trochet et al., 2019). Thus, at a smaller scale, the integration of the landscape resistance could also allow us to obtain more biologically relevant predictions. According to the negative impact of human infrastructure on animal movements, the width of the dispersal buffers may be too large in some regions. Nevertheless, while human infrastructures have negative effects on dispersal events, logging activities could generate waterbodies (for instance ruts) used for reproduction (Cayuela et al., 2014). Thus, landscape management could counterbalance human impacts in harvested forest that were supposed to provide high quality habitat especially during drought (Scheele et al., 2014).

| CON CLUS ION
We demonstrated that incorporating intraspecific levels into nichebased models is a valuable approach to obtaining more biologically relevant projections in response to environmental changes, especially when the species has lineages with divergent bioclimatic niches across its distribution. Predictions at the intraspecific level were generally more optimistic than those made at the species level, notably by their limited decline in response to future climate change.
Taking into account each lineage separately reflected adaptations to 'regional' climates and allowed us to observe specific responses rather than an average global response. However, when dispersal ability was considered, future projections of the two sister lineages of B. v. variegata, particularly the one from the Carpathians, showed a drastic reduction in suitable areas, which could lead to extinction in less than 70 years. In the face of unprecedentedly high current climate velocity, further studies should focus on obtaining better knowledge about species dispersal. Finally, our approach is particularly valuable for species such as European amphibians that exhibit a complex biogeography inherited by isolation and recolonization events that occurred during Pleistocene glaciation and that resulted in distinct lineages.

ACK N OWLED G EM ENTS
We

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data on which this paper is based are available in