Igor Boyer
email: igor.boyer@univ-tours.fr

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

Keywords: Bombina variegata, climate refugia, ecological niche modelling, Europe, intraspecific levels, niche differentiation, range-shift

Aim: 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.

INTRODUCTION

Climate change caused by human activities is one of the main drivers of biodiversity decline worldwide [START_REF] Parmesan | Ecological and Evolutionary Responses to Recent Climate Change[END_REF]. 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 [START_REF] Loarie | The velocity of climate change[END_REF]. Species may respond to those changes through phenotypic plasticity and local adaptation [START_REF] Urban | Plasticity and genetic adaptation mediate amphibian and reptile responses to climate change[END_REF] or may shift their distribution range to track favourable climate conditions [START_REF] Lenoir | Species better track climate warming in the oceans than on lands[END_REF]. Possibilities and magnitudes of range shifts are determined by the species' propensity and ability to disperse [START_REF] Travis | Dispersal and species' responses to climate change[END_REF]. Climateinduced dispersal has been reported in many species [START_REF] Blaustein | Direct and Indirect Effects of Climate Change on Amphibian Populations[END_REF] and likely contributes to the unprecedented biodiversity redistribution reported in recent studies [START_REF] Lenoir | Species better track climate warming in the oceans than on lands[END_REF]. However, the evolvability of thermal performances as well as dispersal potential strongly varies at both interspecific and intraspecific levels [START_REF] Bestion | Dispersal response to climate change: Scaling down to intraspecific variation[END_REF], 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 [START_REF] Araújo | Climate warming and the decline of amphibians and reptiles in Europe[END_REF]. However, those models have historically suffered from conceptual and technical limitations that critically affected the reliability of their predictions [START_REF] Hernandez | The effect of sample size and species characteristics on performance of different species distribution modeling methods[END_REF]. 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 [START_REF] Urban | Improving the forecast for biodiversity under climate change[END_REF]. 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, [START_REF] Bestion | Dispersal response to climate change: Scaling down to intraspecific variation[END_REF][START_REF] Trochet | Variation of preferred body temperatures along an altitudinal gradient: A multi-species study[END_REF]. 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 [START_REF] Gienapp | Responses to climate change in avian migration time-Microevolution versus phenotypic plasticity[END_REF][START_REF] Urban | Plasticity and genetic adaptation mediate amphibian and reptile responses to climate change[END_REF]. 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 [START_REF] Beebee | Amphibian breeding and climate[END_REF][START_REF] Blaustein | Direct and Indirect Effects of Climate Change on Amphibian Populations[END_REF]. Such variation in population responses to climatic conditions could have a genetic base [START_REF] Wilson | Geographic variation in thermal sensitivity of frog jumping[END_REF] 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 [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF][START_REF] Zeisset | Amphibian phylogeography: A model for understanding historical aspects of species distributions[END_REF]. 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 [START_REF] Rowe | Lusitania revisited: A phylogeographic analysis of the natterjack toad Bufo calamita across its entire biogeographical range[END_REF]. 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 [START_REF] Gilbert | Natural selection on thermal preference, critical thermal maxima and locomotor performance[END_REF].

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 [START_REF] Fong | Amphibian Hotspots and Conservation Priorities in Eastern Cuba Identified by Species Distribution Modeling[END_REF]. 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 [START_REF] Hofman | Phylogeography of the fire-bellied toads Bombina: Independent Pleistocene histories inferred from mitochondrial genomes[END_REF][START_REF] Pabijan | Complete mitochondrial genomes resolve phylogenetic relationships within Bombina (Anura: Bombinatoridae)[END_REF]. 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 [START_REF] Cayuela | To breed or not to breed: Past reproductive status and environmental cues drive current breeding decisions in a long-lived amphibian[END_REF]. 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 [START_REF] Phillips | Opening the black box: An open-source release of Maxent[END_REF] 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.

MATERIALS AND METHODS

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 [START_REF] Lescure | Régression de Bombina variegata (Linné, 1758) en France par l'analyse de sa répartition passée et présente[END_REF], 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 [START_REF] Cayuela | To breed or not to breed: Past reproductive status and environmental cues drive current breeding decisions in a long-lived amphibian[END_REF][START_REF] Scheele | Identifying core habitat before it's too late: The case of Bombina variegata, an internationally endangered amphibian[END_REF]. Temperature and precipitation are important determinants of post-metamorphic survival and reproduction [START_REF] Cayuela | Demographic responses to weather fluctuations are context dependent in a long-lived amphibian[END_REF].

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 [START_REF] Fijarczyk | Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories[END_REF][START_REF] Hofman | Phylogeography of the fire-bellied toads Bombina: Independent Pleistocene histories inferred from mitochondrial genomes[END_REF]. 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 [START_REF] Cayuela | Anthropogenic disturbance drives dispersal syndromes, demography, and gene flow in amphibian populations[END_REF]. 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 [START_REF] Boria | Spatial filtering to reduce sampling bias can improve the performance of ecological niche models[END_REF]. 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 [START_REF] Fijarczyk | Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories[END_REF] and [START_REF] Hofman | Phylogeography of the fire-bellied toads Bombina: Independent Pleistocene histories inferred from mitochondrial genomes[END_REF].

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 [START_REF] Cayuela | To breed or not to breed: Past reproductive status and environmental cues drive current breeding decisions in a long-lived amphibian[END_REF][START_REF] Kaplan | Ecological and developmental context of natural selection: maternal effects and thermally induced plasticity in the frog Bombina orientalis[END_REF][START_REF] Reyer | Reproductive ecology of Bombina variegata: Characterisation of spawning ponds[END_REF]. These variables were obtained from 'Climatologies at high resolution for the Earth's land surface area' (CHELSA, [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF], 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 nichebased 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' [START_REF] Broennimann | Measuring ecological niche overlap from occurrence and spatial environmental data: Measuring niche overlap[END_REF][START_REF] Naimi | Where is positional uncertainty a problem for species distribution modelling?[END_REF]. Only variables with a correlation coefficient lower than 0.75 and a VIF of less than 10 were retained [START_REF] Naimi | Where is positional uncertainty a problem for species distribution modelling?[END_REF]. 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 [START_REF] Svenning | Applications of species distribution modeling to paleobiology[END_REF][START_REF] Varela | Using species distribution models in paleobiogeography: A matter of data, predictors and concepts[END_REF]. 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 [START_REF] Harrison | Climate model benchmarking with glacial and mid-Holocene climates[END_REF]. Both model outputs were extracted from the CHELSA database [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF] for the LGM (i.e., 21 000 years ago) and from the WorldClim database (v1.4; [START_REF] Fick | WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas[END_REF] for the mid-Holocene (i.e., 6000 years ago).

We used two other GCMs for future climate conditions expected during the 21st century: CSIRO-Mk3.6.0 and HadGEM2-ES [START_REF] Collier | The CSIRO-Mk3.6.0 Atmosphere-Ocean GCM: participation in CMIP5 and data publication[END_REF][START_REF] Collins | Evaluation of HadGEM2 model[END_REF], which we extracted from the CHELSA and WorldClim databases, respectively. These widely used GCMs were selected because they generate robust predictions [START_REF] Buisson | Uncertainty in ensemble forecasting of species distribution[END_REF][START_REF] Farzaneh | An appropriate general circulation model (GCM) to investigate climate change impact[END_REF] and have been recognized for generating divergent predictions [START_REF] Mitchell | A comprehensive set of highresolution grids of monthly climate for Europe and the globe: The observed record[END_REF]. 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 [START_REF] Baek | Climate change in the 21st century simulated by HadGEM2-AO under representative concentration pathways[END_REF][START_REF] Riahi | RCP 8.5-A scenario of comparatively high greenhouse gas emissions[END_REF][START_REF] Van Vuuren | RCP2.6: Exploring the possibility to keep global mean temperature increase below 2°C[END_REF]. Species redistribution forecasts on several GCMs and climate change scenarios account for future climate condition uncertainties in our work [START_REF] Buisson | Uncertainty in ensemble forecasting of species distribution[END_REF].

Niche similarity between lineages

Prior to modelling, we analysed the bioclimatic niches of the different lineages identified by [START_REF] Fijarczyk | Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories[END_REF] 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 [START_REF] Schoener | The Anolis Lizards of Bimini: Resource Partitioning in a Complex Fauna[END_REF][START_REF] Warren | Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution[END_REF]. 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 [START_REF] Warren | Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution[END_REF].

Niche-based modelling and forecasts

General modelling. We used MaxEnt v3.4.1, one of the most used and robust models [START_REF] Elith | A statistical explanation of MaxEnt for ecologists: Statistical explanation of MaxEnt[END_REF], which also demonstrated good performance in situations with few occurrence data [START_REF] Elith | Novel methods improve prediction of species' distributions from occurrence data[END_REF][START_REF] Hernandez | The effect of sample size and species characteristics on performance of different species distribution modeling methods[END_REF]. 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 [START_REF] Pearson | SPECIES: A Spatial Evaluation of Climate Impact on the Envelope of Species[END_REF].

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 [START_REF] Phillips | A Brief Tutorial on Maxent[END_REF].

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 [START_REF] Araújo | Standards for distribution models in biodiversity assessments[END_REF]. The AUC ranges between 0 and 1, where 0.5 and below indicate that the model has no predictive ability better than random [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (ROC) curve[END_REF]. 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 [START_REF] Allouche | Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS)[END_REF][START_REF] Cohen | A Coefficient of Agreement for Nominal Scales[END_REF]. 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 [START_REF] Liu | Selecting thresholds for the prediction of species occurrence with presence-only data[END_REF], 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" [START_REF] Muscarella | ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for MAXENT ecological niche models[END_REF]. 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 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 distance comprised between 444 and 695 m per year, depending on the age class [START_REF] Cayuela | Toads and roads: Transport infrastructure severely impacts dispersal regardless of life stage[END_REF]. 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 [START_REF] Ehlers | The extent and chronology of Cenozoic Global Glaciation[END_REF][START_REF] Lambeck | Sea Level Change Through the Last Glacial Cycle[END_REF]. We considered continental glaciers to be unsuitable areas. Also, we considered islands to be unreachable and therefore unlikely to be colonized.

RESULTS

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

Current distribution models

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

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

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

Projections under past climates

Due to the large range of the time period (6000 to 21000), 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 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).

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.

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

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

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

Lineage level

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

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.

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

DISCUSSION

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.

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

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. 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 [START_REF] Fijarczyk | Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories[END_REF] that showed this lineage to be the most genetically divergent. The Italian lineage, sometimes considered a distinct species, Bombina pachypus [START_REF] Canestrelli | Genetic diversity and phylogeography of the Apennine yellow-bellied toad Bombina pachypus, with implications for conservation[END_REF], had a significantly similar bioclimatic niche than the Balkan lineage. Populations of these two lineages are geographically separated, apparently with a quasiabsence 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.

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 [START_REF] Hofman | Phylogeography of the fire-bellied toads Bombina: Independent Pleistocene histories inferred from mitochondrial genomes[END_REF]. 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 [START_REF] Fijarczyk | Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories[END_REF][START_REF] Horreo | Postglacial Colonization of Northern Europe by Reptiles[END_REF]. 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 [START_REF] Canestrelli | Genetic diversity and phylogeography of the Apennine yellow-bellied toad Bombina pachypus, with implications for conservation[END_REF].

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 [START_REF] Gotelli | Climate change, genetic markers and species distribution modelling[END_REF], highlighting the importance of integrating intraspecific levels when data are available. [START_REF] Dowell | Divergent lineages and conserved niches: Using ecological niche modeling to examine the evolutionary patterns of the Nile monitor (Varanus niloticus)[END_REF] 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.

Improving models' relevance with dispersal limitation

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 [START_REF] Araújo | Climate warming and the decline of amphibians and reptiles in Europe[END_REF]. 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 [START_REF] Araújo | Climate warming and the decline of amphibians and reptiles in Europe[END_REF].

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.

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 [START_REF] Cayuela | Demographic and genetic approaches to study dispersal in wild animal populations: A methodological review[END_REF]. 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 [START_REF] Broquet | Molecular Estimation of Dispersal for Ecology and Population Genetics[END_REF]. 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 [START_REF] Cayuela | Demographic and genetic approaches to study dispersal in wild animal populations: A methodological review[END_REF]. For instance, [START_REF] Cayuela | Anthropogenic disturbance drives dispersal syndromes, demography, and gene flow in amphibian populations[END_REF] 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 [START_REF] Cayuela | Anthropogenic disturbance drives dispersal syndromes, demography, and gene flow in amphibian populations[END_REF]. 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 [START_REF] Cayuela | Toads and roads: Transport infrastructure severely impacts dispersal regardless of life stage[END_REF][START_REF] Trochet | Influence of substrate types and morphological traits on movement behavior in a toad and newt species[END_REF]. 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 [START_REF] Cayuela | To breed or not to breed: Past reproductive status and environmental cues drive current breeding decisions in a long-lived amphibian[END_REF]. Thus, landscape management could counterbalance human impacts in harvested forest that were supposed to provide high quality habitat especially during drought [START_REF] Scheele | Identifying core habitat before it's too late: The case of Bombina variegata, an internationally endangered amphibian[END_REF]

CONCLUSION

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. (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 became 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 1

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

FIGURE 2

 2 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

  FIGURE 3 Predicted 541 changes in habitat 542 suitability for Bombina 543 variegata lineages in 544 Europe during 21,000 545 and 6,000 bp in 546 comparison to the 547 current period. Red 548 indicates areas that 549 were previously 550 unsuitable but are 551 currently suitable, 552 dark blue indicates 553 areas that have been 554 suitable and are still 555 currently suitable and 556 light blue indicates 557 areas that were 558 previously suitable but 559 are now unsuitable. 560 Grey indicates areas 561 that remain 562 unsuitable, taking into 563 account dispersal 564 ability. Emerged land 565 surface due to lower 566 sea levels and major 567 glaciers has been 568 taken into account in 569 the 21,000 bp models. 570 Projections are 571 valuable under both 572 limited and unlimited 573 dispersal scenarios. 574 Maps were projected 575 under the WGS84 576 system. 577

FIGURE 5

 5 FIGURE 5 Predicted distribution change in Europe by 2070 (period 2061-2080) for Bombina variegata lineages under dispersal limitation (see details in Figure 4).

  Predicted changes in habitat suitability in Europe for Bombina variegata lineages by 2050

		541		3 Predicted
		542	changes in habitat
		543	suitability for Bombina
		544	variegata lineages in
		545	Europe during 21,000
		546	and 6,000 bp in
		547	comparison to the
		548	current period. Red
		549	indicates areas that
		550	were	previously
		551	unsuitable but are
		552	currently		suitable,
		553	dark blue indicates
		554	areas that have been
		555	suitable and are still
		556	currently suitable and
		557	light blue indicates
		558	areas	that	were
		559	previously suitable but
	FIGURE 4	560	are now unsuitable.
		561	Grey indicates areas
		562	that		remain
		563	unsuitable, taking into
		564	account		dispersal
		565	ability. Emerged land
		566	surface due to lower
		567	sea levels and major
		568	glaciers has been
		569	taken into account in
		570	the 21,000 bp models.
		571	Projections	are
		572	valuable under both
		573	limited and unlimited
		574	dispersal	scenarios.
		575	Maps were projected
		576	under the WGS84
		577	system.
	578
	579
	580

TABLE 1

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

				Percentage
	Lineages	Period	GCM	of suitable
				areas
	B. v. variegata	21000 BP	CCSM4 MIROC-ESM	-100 -84.2
	"Carpathian s"	6000 BP	CCSM4 MIROC-ESM	-22.5 -80
	B. v.	21000 BP	CCSM4 MIROC-ESM	50 1792.4
	pachypus	6000 BP	CCSM4 MIROC-ESM	-26.5 -44.1
	B. v. scabra	21000 BP 6000 BP	CCSM4 MIROC-ESM CCSM4 MIROC-ESM	12.3 253.7 -93.8 33.1
	B. v. variegata "western"	21000 BP 6000 BP	CCSM4 MIROC-ESM CCSM4 MIROC-ESM	-68.2 -91 -31.3 -43
	Species level	21000 BP	CCSM4 MIROC-ESM	-88.2 -93.1
	(Bombina variegata )	6000 BP	CCSM4 MIROC-ESM	-37.8 -45.1

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

supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41).

DATA AVAILABILITY STATEMENT

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


W Bekir

O Messal

L Lambourg

F Dubas

S Harmand

A Lebouc

F Gillon

Multi-Physical Model for Sizing Optimization of Electrical Machines

Keywords: Multi-physical model, optimal sizing, automotive applications, model coupling, electrical machines

this abstract outlines a multi-physical modeling approach, which has been developed and incorporated in an optimization process. This approach allows the conceptual design of an electrical motor for electric hybrid vehicle. It is elaborated to be malleable and to satisfy different precision levels in order to address various sizing problems.

I. INTRODUCTION

The design of an electric machine is a complex and often needs iterative process. Designing requires the use of various physics and expertises. (Semi-)analytical models are wildly used for designing and they can be linked to finite-element analysis during the parameter adjustment and validation stages [START_REF] Felden | Fast design process for a complete machine series applying coupled analytical and numerical simulations[END_REF]. A model is always built for a purpose. A design model should take into consideration different physical aspects and be sufficiently fast for using in an optimization process. The designer should be able to manage multiple models via a fast and malleable design model.

The presented study describes a multi-physical model dedicated to the optimization scaling of a radial-flux interior permanent-magnet synchronous machine for automotive traction use.

II. MODEL COUPLING FOR SIZING

The multi-physical machine design model has been established by coupling a generalized nonlinear adaptive magnetic equivalent circuit (MEC) [START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF], the Loss Surface (LS) hysteresis model [START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF], and an aero-thermal model (SAME) [START_REF] Jandaud | Aero-Thermal Optimization of a Heat Sink using Variable Neighbourhood Search[END_REF]. These models are complex and represent just a part of the whole process needed to meet the design requirements of a traction motor in terms of the modeling of efficiency and heating. The multi-physical model proposed herein is malleable with many physical layers and many accuracy levels (viz., Coarse, Medium, and High) that the designer could adjust for each specific problem.

A flowchart of the developed tool is described in Fig. 1. According to the need, it allows to execute separately different layers. For example, if the optimization problem does not include the temperature effects, there is no need to execute the SAME layer (see.Fig. 1) and it is sufficient to execute the two first ones i.e., the coupling (MEC-LS). Furthermore, each layer accuracy can be different from the others. Thus, in a design process, it is possible to reach the solution using coarse but fast model followed by a more accurate model to reach the optimal solution. This homotopic technic considerably saves calculation time and allows good accuracy. In addition, this coupling provides global data using the local ones (spatial and temporal). For instance, it is possible to study the geometric discretization influence of the MEC on the iron loss distribution in the machine magnetic circuit and thus, its impact on the global value in order to adapt the model to the need.

Three variable types are implemented: geometric variables Gi linked to the dimensions, command variables Ci related to the electrical machine operating point torque/speed and the temperature variables linked to the initial temperature T0i of the electrical machine studied part.

The model outputs are obtained as function of the rotor position and/or the electrical period. The multi-physical design model is incorporated in an optimization process. Examples including to temperature, control, geometry modifications and to the driving cycle will be studied in the extended paper.

Fig. 1 .

 1 Fig. 1. A sizing model decomposed into different layers representing the physical phenomena and with different discretizations allowing a compromise between accuracy and computation time.

ACKNOWLEDGMENT This work was supported by RENAULT-SAS, Guyancourt, France, through the program COCTEL (Conception optimale des chaines de Traction ELectrique) financed by the French Agency for Environnement and Energy Management (ADEME), in the program "Véhicule du future, Programme d'Investissements d'Avenir".