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

Igor Boyer, Hugo Cayuela, Romain Bertrand, Francis Isselin-nondedeu

► To cite this version:

Igor Boyer, Hugo Cayuela, Romain Bertrand, Francis Isselin-nondedeu. Improving biological relevance of model projections in response to climate change by considering dispersal amongst lineages in an amphibian. Journal of Biogeography, 2021, 48 (3), pp.561-576. 10.1111/jbi.14019. hal-03213497

HAL Id: hal-03213497

https://amu.hal.science/hal-03213497

Submitted on 3 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

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

Igor Boyer ¹, Hugo Cayuela ², Romain Bertrand ³, Francis Isselin-Nondedeu ¹, ⁴, ⁵

Correspondence

Igor Boyer, UMR 7324 – CNRS CITERES, 33 allée Ferdinand de Lesseps, 37200 TOURS, France.

Email: igor.boyer@univ-tours.fr

[Authors' affiliations]

¹UMR 7324 – CNRS CITERES, 37200 TOURS France

²Institut de Biologie Intégrative et des Systèmes (IBIS), Université Laval, Québec, QC, Canada

³Laboratoire Evolution & Diversité Biologique (EDB), UMR5174 Université de Toulouse Midi-Pyrénées, CNRS, IRD, UPS, 118 route de Narbonne, Bat 4R1, 31062 Toulouse Cedex 9, France

⁴IMBE, UMR Université Aix-Marseille Avignon, 7263-CNRS, 237-IRD IRPNC (Ingénierie de la Restauration des Patrimoines Naturels et Culturels), France

⁵Ecole Polytechnique de l'Université de Tours, Département d'Aménagement et d'Environnement, 33 allée Ferdinand de Lesseps, Tours, France

ORCID

Igor Boyer : https://orcid.org/0000-0002-2870-9391

Hugo Cayuela : https://orcid.org/0000-0002-3529-0736

Romain Bertrand : https://orcid.org/0000-0003-2386-9174

Francis Isselin-Nondedeu : https://orcid.org/0000-0002-9513-4524

ACKNOWLEDGMENTS

We thank Matthieu Willmes from the Direction Régionale de l'Environnement, de l'Aménagement et du Logement Centre-Val de Loire (DREAL), the Office National des Forêts (ONF), Faune-France and all of the naturalistic associations and parks which contributed to the project by sharing the data of occurrence of the species. The funding was supported by the Region Centre-Val de Loire

20

21

22

23

24

25

26

27

28

29

30
31 (France) under the acronym MODELISE. RB's work was supported by the TULIP Laboratory of
32 Excellence (ANR-10-LABX-41).

33 CONFLICT OF INTEREST STATEMENT

34 The authors declare they have no conflict of interest.

35

36

37 ABSTRACT

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

44 **Location:** Europe

45 **Taxon:** *Bombina variegata* and its lineages

46 **Methods:** Using MaxEnt, we fitted correlative niche models for *B. variegata* lineages and also for
47 the species level (*B. variegata sensu lato*) that we projected under Last Glacial Maximum (LGM)
48 and future climate scenarios (RCP2.6, RCP8.5). A comparison of projections was conducted
49 considering both unlimited and limited dispersal abilities.

50 **Results:** We found that the *B. variegata* lineages differed in their bioclimatic niches. In general,  
51 models run without discriminating the lineages showed reduced suitable areas compared to  
52 models run at each lineage, especially those for *Bombina variegata pachypus* and *Bombina  
53 variegata scabra* (southern lineages). Suitable areas identified for the LGM match with the climatic  
54 refugia identified in phylogeographic studies. Projections for the mid-Holocene showed the  
55 increase of suitable areas for the lineages *B. variegata variegata* "Carpathians" and "Western"  
56 while they decreased for the southern lineages. In the future (2040-2080), only models under  
57 limited dispersal predict the extinction of the Carpathian lineage and a reduction in suitable areas  
58 for the other lineages, regardless of periods and scenarios.

59 **Main conclusions:**

60 We demonstrated important shifts in habitat suitability in Europe suggesting range-shifts about
61 hundreds of kilometres in response to past climate changes. However, rapid and unprecedented
62 changes in suitability are expected in the future due to the accelerated climate warming during
63 the 21st. According to their limited dispersal abilities, the persistence of all *B. variegata* lineages
64 over time appears uncertain. We further recommend to integrate the intraspecific levels, when
65 available, in modelling process to catch local variations and obtain more accurate results.
66 KW: *Bombina variegata*, climate refugia, ecological niche modelling, Europe, intraspecific levels,
67 niche differentiation, range-shift

68 1. INTRODUCTION

69 Climate change caused by human activities is one of the main drivers of biodiversity decline
70 worldwide (Parmesan, 2006). One of the main effects of climate change is the shift of thermal
71 isoclines upward and poleward, leading to the redistribution of climate conditions on a global
72 scale (Loarie et al., 2009). Species may respond to those changes through phenotypic plasticity
73 and local adaptation (Urban et al., 2014) or may shift their distribution range to track favourable
74 climate conditions (Lenoir et al., 2020). Possibilities and magnitudes of range shifts are
75 determined by the species' propensity and ability to disperse (Travis et al., 2013). Climate-
76 induced dispersal has been reported in many species (Blaustein et al., 2010) and likely contributes
77 to the unprecedented biodiversity redistribution reported in recent studies (Lenoir et al., 2020).
78 However, the evolvability of thermal performances as well as dispersal potential strongly varies
79 at both interspecific and intraspecific levels (Bestion et al., 2015), therefore making accurate
80 prediction of species response to climate change a challenging task.

81 Niche-based models (also called ‘ecological niche models’ or ‘species distribution models’), relying
82 on correlative relationships between observed species distributions and environmental factors,
83 are powerful tools to predict range shifts driven by future or past climate variations (Araújo et al.,
84 2006). However, those models have historically suffered from conceptual and technical
85 limitations that critically affected the reliability of their predictions (Hernandez et al., 2006). First,
86 relatively few studies explicitly considered species dispersal ability when modelling past and
87 future range shifts, despite the fact that running models with unlimited dispersal overestimates
88 the potential future range, especially for species with moderate to low dispersal ability (Urban et
89 al., 2016). Second, most studies did not consider the intraspecific levels in climate-related
90 phenotypic performances.

91 Recent studies suggest that populations of the same species may differ in their responses to
92 climate variation (see, for instance, Bestion et al., 2015 and Trochet et al., 2018). Differences in
93 population sensitivity to climate fluctuations likely rely on population-specific reaction norms to
94 thermal and hydric conditions and/or adaptation to local climatic conditions (Gienapp et al., 2007;
95 Urban et al., 2014). Those population-specific responses are expected to be particularly broad in
96 ectotherms whose physiology and fitness components (i.e., survival, reproduction and growth)
97 strongly depend on local/regional thermic and hygrometric variation (Beebee, 1995; Blaustein et
98 al., 2010). Such variation in population responses to climatic conditions could have a genetic base
99 (Wilson, 2001) depending on population evolutionary history. This is especially true for
100 ectotherm vertebrates from temperate areas that experienced broad climate variation and
101 complex demographic/evolutionary processes during the Pleistocene period (Hewitt, 2000;
102 Zeisset & Beebee, 2008). In many of those species, the occupancy of multiple glacial refuges led to
103 the emergence of glacial lineages that have experienced repeated periods of allopatry, often
104 punctuated by secondary contacts with more or less intense gene flow (Rowe et al., 2006). Those
105
106
107

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).

108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

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 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?

125
126
127

First, we fitted correlative niche models for each *B. variegata* lineage using MaxEnt (Phillips *et al.*, 2017) 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.

128

129

## 2. MATERIALS AND METHODS

130

### 2.1. Model species

131

132
133
134
135
136
137
138

*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 post-metamorphic survival and reproduction (Cayuela et al., 2016).

139
140
141
142

*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.
143

144

## 2.2. Occurrence data

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

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 GBIF.org (accessed May 2019) GBIF Occurrence Download (DOI: https://doi.org/10.15468/dl.ndtync for *B. v. pachypus*; DOI: https://doi.org/10.15468/dl.86wc6h for the other lineages). 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 x 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 (Boria et al., 2014). 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) and Hofman et al. (2007).

162

163

164

165

166

167

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).

168

169

## 2.3. Climate data

170

171

172

173

174

175

176

177

178

179

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.chelsa-climate.org/) for the period 1979–2013 (30 arcsec resolution, ~1 km). Climatic layers were averaged at a 10 x 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.
180

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).

181

182

183

184

185

186

For past climatic conditions, we used CCSM4 and MIROC-ESM models according to the reviews of Svenning et al. (2011) and Varela 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 (v1.4; Fick & Hijmans, 2017) for the mid-Holocene (i.e., 6000 years ago).

187

188

189

190

191

192

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).

193

194

195

196

197

198

199

200

201

202

203

204

205

## 2.4. Niche similarity between lineages

206

207

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* = 1000). 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).

208

209

210

211

212

213

214

215

216

217

## 2.5. Niche-based modelling and forecasts
218

219

220

221

222

223

224

General modelling. We used MaxEnt v3.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).

225

226

227

228

229

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).

230

231

232

233

234

235

236

237

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.

238

239

240

241

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.

242

243

244

245

246

247

248

249

250

251

*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).

252

253

254

255

*Limited vs. unlimited dispersal scenarios.* We considered the dispersal capacity based on the study from Cayuela 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
256

257

258

259

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.

260

261

262

263

264

265

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.

266

267

268

269

270

271

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.

272

273

## 3. RESULTS

274

### 3.1 Divergence of bioclimatic niches among the *Bombina* glacial lineages

275

276

277

278

279

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.

280

281

282

283

284

285

286

287

288

289

We detected a significant niche similarity between the Italian (*B. v. pachypus*) and the Balkan (*B. v. scabra*) lineages, suggesting they occur in similar climate conditions ($I = 0.693$; see Supporting Information Appendix S3). Other climate niche comparisons among lineages were not significant, showing a high level of intraspecific niche differentiation within the *B. variegata* species. The Carpathian lineage had the smallest bioclimatic niche and always presented the lowest scores of niche similarity with the other three lineages (Figure 1). The maximal niche similarity for this lineage was observed with its sister lineage, western *B. v. variegata* ($I = 0.252$). Whereas the western lineage of *B. v. variegata* had the largest niche, which overlapped greatly with the other lineages, niche similarity was low ($I \min/\max = 0.252/0.428$) as a large part of the climate niche of the western lineage is not shared with any other lineages.

290

291

### 3.2 Current distribution models
292

293

Models at both the species and lineage levels performed well according to the mean values of the AUC, TSS and kappa, which ranged from 0.937 (± 0.036) to 0.961 (± 0.044), 0.760 (± 0.016) to 0.923 (± 0.076) and 0.412 (± 0.101) to 0.781 (± 0.015), respectively (see Supporting Information Appendix S4).

294

295

296

297

According to the models' outputs, the most suitable areas were located in the central Carpathian Mountains for the Carpathian lineage of *B. v. variegata*, in the Balkan Mountains for *B. v. scabra*, in southern Italy for *B. v. pachypus* and in western Europe for the western lineage of *B. v. variegata* (Figure 2a). Models at the species level identified suitable areas that correspond almost exclusively to those identified for the western lineages (i.e., in western Europe) (Figure 2b).

298

299

300

301

302

The bioclimatic variables that contributed the most to the models were 'Precipitation of the Wettest Quarter' (Bio16: 33.4% and 39.8%) for the Carpathian and Italian lineages, 'Mean Temperature of the Driest Quarter' (Bio9: 36.5%) for *B. v. scabra* and 'Precipitation of the Driest Quarter' (Bio17) for the western lineage of *B. v. variegata* and also for the species-level model (53.1% and 52.4%, respectively) (see Supporting Information Appendix S5).

303

304

305

306

307

308

### 3.3. Projections under past climates

309

310

Due to the large range of the time period (6000 to 21000), past projections were made without any dispersal limitation (see Materials and Methods).

311

312

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).

313

314

315

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 Supporting Information 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 toward 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, 6000 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 toward the centre and the south of the Carpathians (Figure 3; Table 1).

316

317

318

319

320

321

322

323

324

325

326

327

328
329

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.

330

331

332

333

334

335

For the Balkan lineage, *B. v. scabra*, a large part of the Balkans appeared suitable, including the Dinaric Alps and Rhodope Mountains during the LGM (Figure 3, see Appendix S1). Results based on climate conditions inferred from the CCSM4 model are more restricted than those based on the MIROC-ESM model (see Appendix S6). Thus, suitable areas were about 12.3% to 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).

336

337

338

339

340

341

342

343

### 3.4. Projections under future climate change

344

345

Results correspond to models run at species or lineage level by considering either unlimited or limited dispersal.

346

347

348

#### Species level

349

According to species distribution models (SDMs) made at the species level, *B. variegata* is expected to undergo a reduction in suitable areas similar to the western lineage of *B. v. variegata*, ranging from -99.6% to -54.1% in 2061–2080 under dispersal limitation. Thus, percentages of predicted suitable areas under dispersal limitation are lower than those obtained for the southern lineages (*B. v. pachypus* and *B. v. scabra*).

350

351

352

353

354

Under unlimited dispersal, forecasts indicated that *B. variegata* will likely not expend its range, and whereas the western lineage of *B. v. variegata* could gain 37.6% more suitable areas in 2061–2080 according to the HadGEM2-ES GCM coupled with the RCP2.6, *B. variegata* (*sensu lato*) could only gain 5.5% more suitable areas. Nevertheless, losses in suitable areas are less significant under unlimited dispersal, ranging from -93.3% to 5.5% (Table 2; see Supporting Information Appendix S7).

355

356

357

358

359

360

361

#### Lineage level

362

The models showed that the two lineages of *B. v. variegata* had the lowest increase of suitable areas, with a mean percentage increase of -99.6% to 170.2% and -39.1% to 37.6%, compared to current areas, in 2061–2080 for the Carpathians and western lineage under unlimited dispersal,

363

364
365
366
367

respectively (all scenarios and periods). Under dispersal limitations, forecasts were pessimistic, with suitable areas ranging from -100% to 44.5% compared to current areas (all scenarios, periods and lineages).

368
369
370
371
372
373
374
375

According to projections based on the RCP8.5 under unlimited dispersal, potential suitable areas were expected to be 2.3 to more than 20 times larger than current areas for *B. v. pachypus* and *B. v. scabra* (Table 2; see Supporting Information Appendix S8, S9). Even if suitable areas were expected to increase for these lineages according to this scenario under unlimited dispersal, the two lineages of *B. v. variegata* should undergo reduction in suitable areas from -97.7% to -59% in 2061-2080. Under dispersal limitations, forecasts obtained for the RCP8.5 were always pessimistic compared to those obtained under unlimited dispersal, even for *B. v. pachypus* and *B. v. scabra*, with a mean percentage of suitable areas between -63.2% and 1.5% in 2061-2080.

376
377
378
379
380
381
382

Thus, SDMs often predicted less loss of suitable areas under unlimited dispersal even if some predictions still indicated an almost complete disappearance of suitable areas (Table 2). Finally, the Carpathian lineage appears to be the most threatened, and only a few areas remain suitable in altitude according to the optimistic climate scenario, while the other models project an almost complete loss of suitable habitat in 2061-2080 under dispersal limitations (Figure 4, 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).

383
384

385

## 4. DISCUSSION

386
387
388
389
390

For this paper, we developed models that simulated the changes in suitable areas for *B. variegata* in response to past and future climate changes by testing the different lineages independently and taking into account the dispersal capacity. This study is the first forecast of the future geographic distribution of bioclimatic suitable areas for this vulnerable amphibian. This is also the first that explored the animal's potential distribution during the LGM and mid-Holocene period.

391
392
393
394

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 (Supporting Information Appendix S10).

395
396
397
398

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.

399

400

### 4.1. Intraspecific level strongly impacts forecasts
401

402

403

404

405

406

407

408

409

410

411

412

413

414

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. According to the niche similarity measures, each lineage had relatively little overlap, and we distinguished two main groups: the southern group including *B. v. pachypus* and *B. v. scabra* and the northern group including the other lineages. With the smallest overlap scores, the Carpathian lineage appears to be the most different, with populations living in a particular climatic environment. Our results are consistent with genetic studies (Fijarczyk et al., 2011) that showed this lineage to be the most genetically divergent. The Italian lineage, sometimes considered a distinct species, *Bombina pachypus* (Canestrelli et al., 2006), had a significantly similar bioclimatic niche than the Balkan lineage. Populations of these two lineages are geographically separated, apparently with a quasi-absence of populations between their current range areas. Thus, the proximity of these two niches is likely due to similarities in the bioclimatic conditions of the Balkan and Italian peninsulas.

415

416

417

418

419

420

421

422

423

424

Models' outputs were different if we considered the species or the lineage level, and simulations for the LGM gave contrasting results, depending on the lineages. According to species models for past projections, no suitable area was present near the Carpathians or in the Balkan Peninsula, even if these areas had served as climatic refugia during the LGM (Hofman et al., 2007). For the Carpathian's lineage, models have identified suitable areas in lowland regions, in the Pannonian Basin. The localization of these suitable areas is consistent with phylogeographic studies and supports the fact that this region served as a climatic refuge for ectotherms during the LGM (Fijarczyk et al., 2011; Horreo & Fitze, 2018). For *B. v. pachypus*, suitable areas were located in lowlands between the Tyrrhenian Sea and the Apennine Mountains and in southern Italy, as suggested by genetic studies during the LGM (Canestrelli et al., 2006).

425

426

427

428

429

430

431

432

433

434

435

436

437

438

Current models at species level indicated that Carpathians and the Balkan and Italian Peninsulas were not suitable to host *B. variegata* populations. By contrast, when we considered the intraspecific level, we showed that the suitability of these areas was higher than expected according to the species model. These results support the fact that distinct populations of the same species may react differently to climate and habitat changes (Gotelli & Stanton-Geddes, 2015), highlighting the importance of integrating intraspecific levels when data are available. Dowell and Hekkala (2016) have shown that the predicted suitable areas differed depending on the lineages of *Varanus niloticus* and therefore should not face the same future threats. In our case, forecasts for the species models were similar to those obtained for the western lineage of *B. v. variegata* (see Appendix S10). One plausible explanation is that the western lineage of *B. v. variegata*, which has the largest distribution and accounts for 91.7% of the occurrences in Europe, largely influenced the model, hiding the diversity within the species models. Consequently, the suitable areas identified in the Italian and Balkan Peninsulas for the southern lineages have been considerably reduced, leading to pessimistic predictions.

439
440

## 4.2. Improving models' relevance with dispersal limitation

441

442

443

444

445

446

447

448

449

450

451

452

453

Dispersal is one of the keys that allows species to migrate and to modify their range in response to environmental changes. Several studies have discussed modelling species distributions using different dispersal scenarios. Two 'dispersal assumptions' are particularly widespread in modelling studies: unlimited dispersal and no dispersal. In the case of unlimited dispersal, it is recognized that forecasts tend to over-predict suitable areas and are unrealistic because the areas predicted as suitable are, in general, unreachable by the species (Araújo et al., 2006). With regard to the no dispersal scenario, forecasts often under-predict suitable areas because they are limited by the current distribution area of the species studied. Even if this latter dispersal assumption appears more realistic than that of unlimited dispersal, it suffers from severe limitations, mainly because many species have the ability to track their bioclimatic niches to some extent, emphasizing the importance of integrating a biologically relevant dispersal distance (Araújo et al., 2006).

454

455

456

457

458

459

460

461

462

463

464

Under unlimited dispersal, our models often predicted an increase in new suitable areas for all lineages, mainly predicted to be located in north or at altitude, out of the lineages' current range. However, under limited dispersal scenarios, despite predicted new suitable areas, results showed a global reduction of suitable areas and sometimes an almost total disappearance of them. Consequently, three of the four lineages (*B. v. scabra*, Carpathians and western *B. v. variegata*) are predicted to lose all of their suitable areas in the future, according to at least one model. The discrepancy between the results under limited and unlimited dispersal could be explained by the dispersal abilities of the yellow-bellied toad. Indeed, even if suitable areas remain in Europe over time, many of them will likely be unreachable, according to our dispersal buffer. Thus, only a fraction of these suitable areas will likely be reached, mainly in altitude or to the north of the lineages' current range.

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

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
480

481

482

483

484

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.

485

486

487

488

489

490

491

492

493

494

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)

495

496

5. CONCLUSION

497

498

499

500

501

502

503

504

505

506

507

508

509

510

We demonstrated that incorporating intraspecific levels into niche-based 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.

511

512

513

514

515

516
<figure><img src="image_2.png" /><figcaption>FIGURE 1 Bioclimatic niche similarities amongst the four *Bombina variegata* lineages. The horizontal axis in the correlation circle represents the first component (explaining 46.29% of the total variation), whereas the vertical axis represents the second component (22.81% of the total variation). Bio2 = Mean Diurnal Range, Bio5 = Maximum Temperature of Warmest Month, Bio6 = Minimum Temperature of Coldest Month, Bio9 = Mean Temperature of Driest Quarter, Bio16 = Precipitation of Wettest Quarter, Bio17 = Precipitation of Driest Quarter, Bio18 = Precipitation of Warmest Quarter. Grey shading represents the density of the occurrences by cell, with darker shading indicating higher density of occurrences. Solid and dashed contour lines represent 100% and 50% of the available environment, respectively. Green colour represents niche overlap between two lineages, whereas yellow and lavender represent non-overlapping niches. PC1 and PC2 are the first two axes of principal component analysis (PCA). Similarity I indices between two lineages are mentioned in parentheses above the representation of related bioclimatic niches. A significant value denotes a pair of species that are ecologically similar (*, p ≤ 0.05).</figcaption></figure>
<figure><img src="image_3.png" /></figure>

FIGURE 2 (a) Habitat suitability of each lineage of Bombina variegata in Europe, performed separately, within their distribution areas. Models of each lineage were fitted separately, and their outputs were clipped by the current distribution limits of each lineage and pooled in the same map. The letters correspond to the lineages and locate their distribution areas: C = *Bombina variegata variegata* 'Carpathians'; P = *Bombina variegata pachypus*; S = *Bombina variegata scabra*; V = *Bombina variegata variegata* 'western'. The dotted white line separates the distribution areas of B. v. scabra from the Carpathian and western lineages of B. v. variegata. (b) Habitat suitability of Bombina variegata in Europe (species level). Maps were projected under the WGS84 system.
<figure><img src="image_4.png" /></figure>

FIGURE 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. Projections are valuable under both limited and unlimited dispersal scenarios. Maps were projected under the WGS84 system.

577

578

579

580
<figure><img src="image_5.png" /></figure>

FIGURE 4 Predicted changes in habitat suitability in Europe for *Bombina variegata* lineages by 2050 (period 2041–2060) under dispersal limitation. Because predictions for each lineage were made in different areas, they were pooled in the same map. The colours show the differences from the current distributions: red shows lost suitable areas (i.e., areas that will become unsuitable), dark blue indicates suitable areas that are predicted to persist, light blue indicates new suitable areas and grey indicates areas that remain unsuitable, taking into account dispersal ability. A dashed line separates the predictions for *Bombina variegata scabra* (in the south) from those of the Carpathian and western lineages of *Bombina variegata variegata*. Maps were projected under the WGS84 system.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 Predicted distribution change in Europe by 2070 (period 2061–2080) for *Bombina variegata* lineages under dispersal limitation (see details in Figure 4).</figcaption></figure>

<table><caption>TABLE 1 Predicted change of suitable areas in the past compared to the current period under unlimited dispersal. Species models did not take into account intraspecific diversity (i.e., *Bombina variegata* sensu lato).</caption><thead><tr><th>Lineages</th><th>Period</th><th>GCM</th><th>Percentage of suitable areas</th></tr></thead><tbody><tr><td rowspan="2"><i>B. v. variegata</i><br>"Carpathian s"</td><td rowspan="2">21000 BP</td><td>CCSM4</td><td>-100</td></tr><tr><td>MIROC-ESM</td><td>-84.2</td></tr><tr><td rowspan="2"><i>B. v. pachypus</i></td><td rowspan="2">6000 BP</td><td>CCSM4</td><td>-22.5</td></tr><tr><td>MIROC-ESM</td><td>-80</td></tr><tr><td rowspan="2"><i>B. v. scabra</i></td><td rowspan="2">21000 BP</td><td>CCSM4</td><td>50</td></tr><tr><td>MIROC-ESM</td><td>1792.4</td></tr><tr><td rowspan="2"><i>B. v. variegata</i><br>"western"</td><td rowspan="2">6000 BP</td><td>CCSM4</td><td>-26.5</td></tr><tr><td>MIROC-ESM</td><td>-44.1</td></tr><tr><td rowspan="2"><i>B. v. scabra</i></td><td rowspan="2">21000 BP</td><td>CCSM4</td><td>12.3</td></tr><tr><td>MIROC-ESM</td><td>253.7</td></tr><tr><td rowspan="2"><i>B. v. variegata</i><br>"western"</td><td rowspan="2">6000 BP</td><td>CCSM4</td><td>-93.8</td></tr><tr><td>MIROC-ESM</td><td>33.1</td></tr><tr><td rowspan="2"><i>B. v. variegata</i><br>"western"</td><td rowspan="2">21000 BP</td><td>CCSM4</td><td>-68.2</td></tr><tr><td>MIROC-ESM</td><td>-91</td></tr><tr><td rowspan="2"><i>B. v. variegata</i><br>"western"</td><td rowspan="2">6000 BP</td><td>CCSM4</td><td>-31.3</td></tr><tr><td>MIROC-ESM</td><td>-43</td></tr><tr><td rowspan="2"><i>B. v. variegata</i><br>"western"</td><td rowspan="2">21000 BP</td><td>CCSM4</td><td>-88.2</td></tr><tr><td>MIROC-ESM</td><td>-93.1</td></tr><tr><td rowspan="2"><i>B. v. variegata</i><br>"western"</td><td rowspan="2">6000 BP</td><td>CCSM4</td><td>-37.8</td></tr><tr><td>MIROC-ESM</td><td>-45.1</td></tr></tbody></table>
<table><caption>TABLE 2 Predicted change of suitable areas in the future compared to the current period under limited and unlimited dispersal. Species models did not take into account intraspecific diversity (i.e., *Bombina variegata* sensu lato).</caption><thead><tr><th rowspan="2">Lineages</th><th rowspan="2">Period</th><th rowspan="2">RCP</th><th rowspan="2">GCM</th><th colspan="2">Percentage of suitable areas compared to current (%)</th></tr><tr><th>Limited dispersal</th><th>Unlimited dispersal</th></tr></thead><tbody><tr><td rowspan="6">B. v. variegata "Carpathians"</td><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-77.4<br>-100</td><td>92.7<br>-99.6</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-89.4<br>-100</td><td>-50.7<br>-99.2</td></tr><tr><td rowspan="2">2070</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-62.9<br>-100</td><td>170.2<br>-97.5</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-100<br>-100</td><td>-76.3<br>-73</td></tr><tr><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-16.4<br>-65.8</td><td>543.2<br>257</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-6.1<br>-77.2</td><td>1688.6<br>411.4</td></tr><tr><td rowspan="6">B. v. pachypus</td><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-28.7<br>-41.4</td><td>1051.7<br>310.4</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>1.5<br>-43.8</td><td>2103<br>1342.6</td></tr><tr><td rowspan="2">2070</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-6.6<br>-97.9</td><td>151.8<br>33</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-19.7<br>-91.3</td><td>234.8<br>121</td></tr><tr><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>44.5<br>-100</td><td>176.9<br>-58.8</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-14.5<br>-63.2</td><td>535.5<br>359.3</td></tr><tr><td rowspan="6">B. v. scabra</td><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-68.9<br>-85</td><td>-39.1<br>35.3</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-92<br>-93.9</td><td>-70.2<br>-3.7</td></tr><tr><td rowspan="2">2070</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-56.7<br>-84.4</td><td>6.9<br>37.6</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-100<br>-97.3</td><td>-97.7<br>-59</td></tr><tr><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-62.4<br>-86.8</td><td>-15.1<br>1.7</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-89.5<br>-94.8</td><td>-51.2<br>-30.5</td></tr><tr><td rowspan="6">Species level (Bombina variegata)</td><td rowspan="2">2050</td><td>RCP2.6</td><td>CSIRO<br>HadGEM2-ES</td><td>-54.1<br>-83.5</td><td>5.5<br>-0.9</td></tr><tr><td>RCP8.5</td><td>CSIRO<br>HadGEM2-ES</td><td>-99.6<br>-96.8</td><td>-93.3<br>-46.2</td></tr></tbody></table>
635 Allouche, O., Tsoar, A., & Kadmon, R. (2006). Assessing the accuracy of species distribution models:

636 Prevalence, kappa and the true skill statistic (TSS). *Journal of Applied Ecology*, **43**(6), 1223-1232.

637 https://doi.org/10.1111/j.1365-2664.2006.01214.x

638 Araújo, M. B., Thuiller, W., & Pearson, R. G. (2006). Climate warming and the decline of amphibians

639 and reptiles in Europe. *Journal of Biogeography*, **33**(10), 1712-1728.

640 https://doi.org/10.1111/j.1365-2699.2006.01482.x

641 Araújo, M. B., Anderson, R. P., Márcia Barbosa, A., Beale, C. M., Dormann, C. F., Early, R., ... Rahbek,

642 C. (2019). Standards for distribution models in biodiversity assessments. *Science Advances*, **5**(1),

643 eaat4858. https://doi.org/10.1126/sciadv.aat4858

644 Baek, H.-J., Lee, J., Lee, H.-S., Hyun, Y.-K., Cho, C., Kwon, W.-T., ... Byun, Y.-H. (2013). Climate change

645 in the 21st century simulated by HadGEM2-AO under representative concentration pathways.

646 Asia-Pacific Journal of Atmospheric Sciences, **49**(5), 603-618. https://doi.org/10.1007/s13143-

647 013-0053-7

648 Beebee, T. J. C. (1995). Amphibian breeding and climate. *Nature*, **374**(6519), 219-220.

649 https://doi.org/10.1038/374219a0

650 Bestion, E., Clobert, J., & Cote, J. (2015). Dispersal response to climate change: Scaling down to

651 intraspecific variation. *Ecology Letters*, **18**(11), 1226-1233. https://doi.org/10.1111/ele.12502

652 Blaustein, A. R., Walls, S. C., Bancroft, B. A., Lawler, J. J., Searle, C. L., & Gervasi, S. S. (2010). Direct

653 and Indirect Effects of Climate Change on Amphibian Populations. *Diversity*, **2**(2), 281-313.

654 https://doi.org/10.3390/d2020281

655 Boria, R. A., Olson, L. E., Goodman, S. M., & Anderson, R. P. (2014). Spatial filtering to reduce

656 sampling bias can improve the performance of ecological niche models. *Ecological Modelling*, **275**,

657 73-77. https://doi.org/10.1016/j.ecolmodel.2013.12.012

658 Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., Petitpierre, B., Pellissier, L., Yoccoz, N. G., ...

659 Guisan, A. (2012). Measuring ecological niche overlap from occurrence and spatial environmental

660 data: Measuring niche overlap. *Global Ecology and Biogeography*, **21**(4), 481-497.

661 https://doi.org/10.1111/j.1466-8238.2011.00698.x

662 Broquet, T., & Petit, E. J. (2009). Molecular Estimation of Dispersal for Ecology and Population

663 Genetics. *Annual Review of Ecology, Evolution, and Systematics*, **40**(1), 193-216.

664 https://doi.org/10.1146/annurev.ecolsys.110308.120324

665 Buisson, L., Thuiller, W., Casajus, N., Lek, S., & Grenouillet, G. (2010). Uncertainty in ensemble

666 forecasting of species distribution. *Global Change Biology*, **16**(4), 1145-1157.

667 https://doi.org/10.1111/j.1365-2486.2009.02000.x
668

669

670

671

Canestrelli, D., Cimmaruta, R., Costantini, V., & Nascetti, G. (2006). Genetic diversity and phylogeography of the Apennine yellow-bellied toad *Bombina pachypus*, with implications for conservation. *Molecular Ecology*, 15(12), 3741-3754. https://doi.org/10.1111/j.1365-294X.2006.03055.x

672

673

674

Cayuela, H., Arsovski, D., Thirion, J.-M., Bonnaire, E., Pichenot, J., Boitaud, S., ... Besnard, A. (2016). Demographic responses to weather fluctuations are context dependent in a long-lived amphibian. *Global Change Biology*, 22(8), 2676-2687. https://doi.org/10.1111/gcb.13290

675

676

677

678

Cayuela, H., Besnard, A., Bonnaire, E., Perret, H., Rivoalen, J., Miaud, C., & Joly, P. (2014). To breed or not to breed: Past reproductive status and environmental cues drive current breeding decisions in a long-lived amphibian. *Oecologia*, 176(1), 107-116. https://doi.org/10.1007/s00442-014-3003-x

679

680

681

Cayuela, H., Besnard, A., Cote, J., Laporte, M., Bonnaire, E., Pichenot, J., ... Léna, J.-P. (2020). Anthropogenic disturbance drives dispersal syndromes, demography, and gene flow in amphibian populations. *Ecological monographs*, 90(2). https://doi.org/10.5061/dryad.q83bk3jdz

682

683

Cayuela, H., Bonnaire, E., & Besnard, A. (2018). Toads and roads: Transport infrastructure severely impacts dispersal regardless of life stage. *BioRxiv*. https://doi.org/10.1101/349639

684

685

686

687

Cayuela, H., Rougemont, Q., Prunier, J. G., Moore, J.-S., Clobert, J., Besnard, A., & Bernatchez, L. (2018). Demographic and genetic approaches to study dispersal in wild animal populations: A methodological review. *Molecular Ecology*, 27(20), 3976-4010. https://doi.org/10.1111/mec.14848

688

689

Cohen, J. (1960). A Coefficient of Agreement for Nominal Scales. *Educational and Psychological Measurement*, 20(1), 37-46. https://doi.org/10.1177/001316446002000104

690

691

692

Collier, M. A., Jeffrey, S. J., Rotstayn, L. D., Wong, K. K.-H., Dravitzki, S. M., Moeseneder, C., ... Atif, M. (2011). The CSIRO-Mk3.6.0 Atmosphere-Ocean GCM: participation in CMIP5 and data publication. *International Congress on Modelling and Simulation-MODSIM*, 7.

693

694

Collins, W. J., Bellouin, N., Doutriaux-Boucher, M., Gedney, N., Hinton, T., Jones, C. D., ... Kim, J. (2008). Evaluation of HadGEM2 model. *Hadley Cent. Tech. Note*, 47.

695

696

697

Dowell, S. A., & Hekkala, E. R. (2016). Divergent lineages and conserved niches: Using ecological niche modeling to examine the evolutionary patterns of the Nile monitor (*Varanus niloticus*). *Evolutionary Ecology*, 30(3), 471-485. https://doi.org/10.1007/s10682-016-9818-7

698

699

Ehlers, J., & Gibbard, P. L. (2007). The extent and chronology of Cenozoic Global Glaciation. *Quaternary International*, 164-165, 6-20. https://doi.org/10.1016/j.quaint.2006.10.008

700

701

702

Elith, J., H. Graham, C., P. Anderson, R., Dudík, M., Ferrier, S., Guisan, A., ... E. Zimmermann, N. (2006). Novel methods improve prediction of species' distributions from occurrence data. *Ecography*, 29(2), 129-151. https://doi.org/10.1111/j.2006.0906-7590.04596.x
703

Elith, J., Phillips, S. J., Hastie, T., Dudík, M., Chee, Y. E., & Yates, C. J. (2011). A statistical explanation of MaxEnt for ecologists: Statistical explanation of MaxEnt. *Diversity and Distributions*, *17*(1), 43–57. https://doi.org/10.1111/j.1472-4642.2010.00725.x

704

705

706

Farzaneh, M. R., Eslamian, S., Samadi, S. Z., & Akbarpour, A. (2012). An appropriate general circulation model (GCM) to investigate climate change impact. *International Journal of Hydrology* *Science and Technology*, *2*(1), 34. https://doi.org/10.1504/IJHST.2012.045938

707

708

709

Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology*, *37*(12), 4302–4315. https://doi.org/10.1002/joc.5086

710

711

712

Fijarczyk, A., Nadachowska, K., Hofman, S., Litvinchuk, S. N., Babik, W., Stuglik, M., ... Szymura, J. M. (2011). Nuclear and mitochondrial phylogeography of the European fire-bellied *toads Bombina bombina* and *Bombina variegata* supports their independent histories. *Molecular Ecology*, *20*(16), 3381–3398. https://doi.org/10.1111/j.1365-294X.2011.05175.x

713

714

715

716

Fong G., A., Viña Dávila, N., & López-Iborra, G. M. (2015). Amphibian Hotspots and Conservation Priorities in Eastern Cuba Identified by Species Distribution Modeling. *Biotropica*, *47*(1), 119–127. https://doi.org/10.1111/btp.12186

717

718

719

Gienapp, P., Leimu, R., & Merilä, J. (2007). Responses to climate change in avian migration time— Microevolution versus phenotypic plasticity. *Climate Research*, *35*, 25–35. https://doi.org/10.3354/cr00712

720

721

722

Gilbert, A. L., & Miles, D. B. (2017). Natural selection on thermal preference, critical thermal maxima and locomotor performance. *Proceedings of the Royal Society B: Biological Sciences*, *284*(1860), 20170536. https://doi.org/10.1098/rspb.2017.0536

723

724

725

Gotelli, N. J., & Stanton-Geddes, J. (2015). Climate change, genetic markers and species distribution modelling. *Journal of Biogeography*, *42*(9), 1577–1585. https://doi.org/10.1111/jbi.12562

726

727

Hanley, J. A., & McNeil, B. J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. *Radiology*, *143*(1), 29–36. https://doi.org/10.1148/radiology.143.1.7063747

728

729

730

Harrison, S. P., Bartlein, P. J., Brewer, S., Prentice, I. C., Boyd, M., Hessler, I., ... Willis, K. (2014). Climate model benchmarking with glacial and mid-Holocene climates. *Climate Dynamics*, *43*(3–4), 671–688. https://doi.org/10.1007/s00382-013-1922-6

731

732

733

Hernandez, P. A., Graham, C. H., Master, L. L., & Albert, D. L. (2006). The effect of sample size and species characteristics on performance of different species distribution modeling methods. *Ecography*, *29*(5), 773–785. https://doi.org/10.1111/j.0906-7590.2006.04700.x

734

735

736

Hewitt, G. (2000). The genetic legacy of the Quaternary ice ages. *Nature*, *405*(6789), 907–913. https://doi.org/10.1038/35016000

737
738

Hill, J. K., Griffiths, H. M., & Thomas, C. D. (2011). Climate Change and Evolutionary Adaptations at

739

Species' Range Margins. *Annual Review of Entomology*, 56(1), 143–159.

740

https://doi.org/10.1146/annurev-ento-120709-144746

741

Hofman, S., Spolsky, C., Uzzell, T., Cogălniceanu, D., Babik, W., & Szymura, J. M. (2007).

742

Phylogeography of the fire-bellied toads Bombina: Independent Pleistocene histories inferred

743

from mitochondrial genomes. *Molecular Ecology*, 16(11), 2301–2316.

744

https://doi.org/10.1111/j.1365-294X.2007.03309.x

745

Horreo, J. L., & Fitze, P. S. (2018). Postglacial Colonization of Northern Europe by Reptiles. In P.

746

Pontarotti (Ed.), *Origin and Evolution of Biodiversity* (pp. 197–214). Springer International

747

Publishing. https://doi.org/10.1007/978-3-319-95954-2_12

748

Kaplan, R. H., & Phillips, P. C. (2006). Ecological and developmental context of natural selection:

749

maternal effects and thermally induced plasticity in the frog *Bombina orientalis*. *Evolution*, 60(1),

750

142–156. https://doi.org/10.1111/j.0014-3820.2006.tb01089.x

751

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., ... Kessler, M. (2017).

752

Climatologies at high resolution for the earth's land surface areas. *Scientific Data*, 4, 170122.

753

https://doi.org/10.1038/sdata.2017.122

754

Lambeck, K. (2001). Sea Level Change Through the Last Glacial Cycle. *Science*, 292(5517), 679–

755

686. https://doi.org/10.1126/science.1059549

756

Lenoir, J., Bertrand, R., Comte, L., Bourgeaud, L., Hattab, T., Murienne, J., & Grenouillet, G. (2020).

757

Species better track climate warming in the oceans than on lands. *Nature Ecology and Evolution*

758

4, 1044–1059. https://doi.org/10.1038/s41559-020-1198-2

759

Lescure, J., Pichenot, J., & Cochard, P.-O. (2011). *Régression de Bombina variegata (Linné, 1758) en*

760

*France par l'analyse de sa répartition passée et présente*. 38.

761

Liu, C., White, M., & Newell, G. (2013). Selecting thresholds for the prediction of species occurrence

762

with presence-only data. *Journal of Biogeography*, 40(4), 778–789.

763

https://doi.org/10.1111/jbi.12058

764

Loarie, S. R., Duffy, P. B., Hamilton, H., Asner, G. P., Field, C. B., & Ackerly, D. D. (2009). The velocity

765

of climate change. *Nature*, 462(7276), 1052–1055. https://doi.org/10.1038/nature08649

766

Mitchell, T. D., Carter, T. R., Jones, P. D., Hulme, M., & New, M. (2004). A comprehensive set of high-

767

resolution grids of monthly climate for Europe and the globe: The observed record (1901-2000)

768

and 16 scenarios (2001-2100). *Tyndall Working Paper*, 30.

769

Muscarella, R., Galante, P. J., Soley-Guardia, M., Boria, R. A., Kass, J. M., Uriarte, M., & Anderson, R.

770

P. (2014). ENMeval: An R package for conducting spatially independent evaluations and

771

estimating optimal model complexity for MAXENT ecological niche models. *Methods in Ecology and Evolution*, 5(11), 1198–1205. https://doi.org/10.1111/2041-210X.12261

772
773

Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K., & Toxopeus, A. G. (2014). Where is positional uncertainty a problem for species distribution modelling? *Ecography*, 37(2), 191–203.

774

775

https://doi.org/10.1111/j.1600-0587.2013.00205.x

776

Pabijan, M., Wandycz, A., Hofman, S., Węcek, K., Piwczyński, M., & Szymura, J. M. (2013). Complete mitochondrial genomes resolve phylogenetic relationships within Bombina (Anura: Bombinatoridae). *Molecular Phylogenetics and Evolution*, 69(1), 63–74.

777

778

779

https://doi.org/10.1016/j.ympev.2013.05.007

780

Parmesan, C. (2006). *Ecological and Evolutionary Responses to Recent Climate Change*. *Annual Review of Ecology, Evolution, and Systematics*, 37(1), 637–669.

781

782

https://doi.org/10.1146/annurev.ecolsys.37.091305.110100

783

Pearson, R. G., Dawson, T. P., Berry, P. M., & Harrison, P. A. (2002). SPECIES: A Spatial Evaluation of Climate Impact on the Envelope of Species. *Ecological Modelling*, 154(3), 289–300.

784

785

https://doi.org/10.1016/S0304-3800(02)00056-X

786

Phillips, S. J. (2017). *A Brief Tutorial on Maxent*. 39.

787

Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E., & Blair, M. E. (2017). Opening the black box: An open-source release of Maxent. *Ecography*, 40(7), 887–893.

788

789

https://doi.org/10.1111/ecog.03049

790

Reyer, H.-U., & Barandun, J. (1997). Reproductive ecology of Bombina variegata: Characterisation of spawning ponds. *Amphibia-Reptilia*, 18(2), 143–154.

791

792

https://doi.org/10.1163/156853897X00035

793

Riahi, K., Rao, S., Krey, V., Cho, C., Chirkov, V., Fischer, G., ... Rafaj, P. (2011). RCP 8.5—A scenario of comparatively high greenhouse gas emissions. *Climatic Change*, 109(1–2), 33–57.

794

795

https://doi.org/10.1007/s10584-011-0149-y

796

Rowe, G., Harris, D. J., & Beebee, T. J. C. (2006). Lusitania revisited: A phylogeographic analysis of the natterjack toad Bufo calamita across its entire biogeographical range. *Molecular Phylogenetics and Evolution*, 39(2), 335–346. https://doi.org/10.1016/j.ympev.2005.08.021

797

798

799

Scheele, B. C., Boyd, C. E., Fischer, J., Fletcher, A. W., Hanspach, J., & Hartel, T. (2014). Identifying core habitat before it's too late: The case of Bombina variegata, an internationally endangered amphibian. *Biodiversity and Conservation*, 23(3), 775–780. https://doi.org/10.1007/s10531-014-0624-7

800

801

802

803

Schoener, T. W. (1968). The Anolis Lizards of Bimini: Resource Partitioning in a Complex Fauna.

804

*Ecology*, 49(4), 704–726. https://doi.org/10.2307/1935534

805

Svenning, J.-C., Fløjgaard, C., Marske, K. A., Nógues-Bravo, D., & Normand, S. (2011). Applications of species distribution modeling to paleobiology. *Quaternary Science Reviews*, 30(21–22), 2930–2947. https://doi.org/10.1016/j.quascirev.2011.06.012

806

807
808

Travis, J. M. J., Delgado, M., Bocedi, G., Baguette, M., Bartoń, K., Bonte, D., ... Bullock, J. M. (2013).

809

Dispersal and species' responses to climate change. *Oikos*, **122**(11), 1532-1540.

810

https://doi.org/10.1111/j.1600-0706.2013.00399.x

811

Trochet, A., Dupoué, A., Souchet, J., Bertrand, R., Deluen, M., Murarasu, S., ... Aubret, F. (2018).

812

Variation of preferred body temperatures along an altitudinal gradient: A multi-species study.

813

*Journal of Thermal Biology*, **77**, 38-44. https://doi.org/10.1016/j.jtherbio.2018.08.002

814

Trochet, A., Le Chevalier, H., Calvez, O., Ribéron, A., Bertrand, R., & Blanchet, S. (2019). Influence

815

of substrate types and morphological traits on movement behavior in a toad and newt species.

816

*PeerJ*, **6**, e6053. https://doi.org/10.7717/peerj.6053

817

Urban, M. C., Bocedi, G., Hendry, A. P., Mihoub, J.-B., Peer, G., Singer, A., ... Travis, J. M. J. (2016).

818

Improving the forecast for biodiversity under climate change. *Science*, **353**(6304), aad8466-

819

aad8466. https://doi.org/10.1126/science.aad8466

820

Urban, M. C., Richardson, J. L., & Freidenfelds, N. A. (2014). Plasticity and genetic adaptation

821

mediate amphibian and reptile responses to climate change. *Evolutionary Applications*, **7**(1), 88-

822

103. https://doi.org/10.1111/eva.12114

823

van Vuuren, D. P., Stehfest, E., den Elzen, M. G. J., Kram, T., van Vliet, J., Deetman, S., ... van Ruijven,

824

B. (2011). RCP2.6: Exploring the possibility to keep global mean temperature increase below 2°C.

825

*Climatic Change*, **109**(1-2), 95-116. https://doi.org/10.1007/s10584-011-0152-3

826

Varela, S., Lobo, J. M., & Hortal, J. (2011). Using species distribution models in paleobiogeography:

827

A matter of data, predictors and concepts. *Palaeogeography, Palaeoclimatology, Palaeoecology*,

828

**310**(3-4), 451-463. https://doi.org/10.1016/j.palaeo.2011.07.021

829

Warren, D. L., Glor, R. E., & Turelli, M. (2008). Environmental niche equivalency versus

830

conservatism: quantitative approaches to niche evolution. *Evolution*, **62**(11), 2868-2883.

831

https://doi.org/10.1111/j.1558-5646.2008.00482.x

832

Wilson, R. S. (2001). Geographic variation in thermal sensitivity of frog jumping. *Journal of*

833

*Experimental Biology*, **204**(24), 10.

834

Zeisset, I., & Beebee, T. J. C. (2008). Amphibian phylogeography: A model for understanding

835

historical aspects of species distributions. *Heredity*, **101**(2), 109-119.

836

https://doi.org/10.1038/hdy.2008.30

837

838

839

840

841
842

843

844

845 BIOSKETCH

846 Igor Boyer is a PhD student in Ecology at the Cltés, TERritoires, Environnement et Sociétés (CITERES) lab of the National Center for Scientific Research (CNRS). He is particularly interested in the consequences of global changes on future distribution of European herpetofauna. He asks questions on species biogeography assessing the environmental factors that shape their distributions and therefore the potential consequences of future changes using modeling approaches at different geographical scales.

852 Author contributions: I.B., F.I.N. and H.C. conceived the ideas; I.B. performed the analysis, R.B., F.I.N. co-analysed some of the data; and I.B. wrote, R.B., H.C. and F.I.N reviewed the paper. F.I.N. managed the whole project.

855

856 DATA AVAILABILITY STATEMENT

857 The data on which this paper is based are available in the Figshare digital repository available at
858 https://doi.org/10.6084/m9 figshare.12982016.v2.