A novel tool to assess the effect of intraspecific spatial niche variation on species distribution shifts under climate change

Aim: Niche-based models often ignore spatial variation in the climatic niche of a species across its occupied range and the related variation in the response to changing climate conditions. This assumption may lead to inaccurate predictions of species distribution shifts under climate change. Models have been developed to address this issue, but most of them depend upon prior knowledge on evolutionary lineages, phenotypic traits or ecological processes underlying local adaptation or adaptive plasticity. As such information is often lacking, these models are not frequently used to predict distribution shifts for many species. This limits our ability to explore general patterns of change across species. Innovation: Here, we propose a modelling framework that can be applied across a large sample of species to assess their distribution shifts under future climate while exploring the effect of intraspecific spatial variation in the response to climate conditions. The proposed approach does not require a detailed understanding of the processes underlying such variation. The geographical distribution of a species is split into spatial subsets along the gradient of occupied climate conditions. These sub-sets are considered as proxies for intraspecific spatial niche variation. Local models are built with each subset and their predictions are assembled across the study area


| INTRODUC TI ON
Niche-based models are developed to identify climate conditions suitable for a species by estimating correlations between the presence of the species and local climate conditions.Niche-based models are then used to explore potential shifts in species distributions under future climate conditions (Elith, Kearney, & Phillips, 2010;Franklin, 2013).The reliability of these statistical models has been questioned due to their limited ability to integrate key ecological and evolutionary processes affecting changes in species distributions (Franklin, 2013;Schweiger et al., 2012;Thuiller et al., 2013).For instance, the most frequently used approaches consider a species as a fixed entity in the modelling process and assume that all the populations respond to changing climate conditions in a similar way throughout the occupied geographical range (Ikeda et al., 2017;Moran, Hartig, & Bell, 2016;Pearman, Guisan, Broennimann, & Randin, 2008).
However, it has been shown that there can be important geographical intraspecific variation in the sensitivity and response to changing climate conditions (Mills et al., 2017;Nice et al., 2019;Peterson, Doak, & Morris, 2019).Intraspecific variation in climatic requirements and tolerance may sometimes be as large as interspecific variation (e.g.Díaz-Almeyda et al., 2017;Malyshev et al., 2016) and arises from local adaptation or adaptive plasticity along the gradient of climate conditions (Macdonald, Llewelyn, & Phillips, 2018;Moran et al., 2016;Vargas et al., 2017).Standard niche-based models may, therefore, incorrectly assess the full range of suitable climate conditions for the species and may produce inaccurate or even wrong predictions of distribution shifts (Hällfors et al., 2016;Nice et al., 2019;Peterson et al., 2019;Smith, Godsoe, Rodríguez-Sánchez, Wang, & Warren, 2019;Valladares et al., 2014).
Incorporating intraspecific evolutionary and ecological mechanisms in process-based models that deal with distribution shifts has already been proposed (Chevin, Lande, & Mace, 2010;Cotto et al., 2017;Peterson et al., 2019).Yet, the development of such models is still restricted due to the limited availability of such quantitative information at large spatial scales (Briscoe et al., 2019;Morin & Thuiller, 2009).Alternative approaches that integrate basic information about mechanisms into niche-based models have recently been developed (Diniz-Filho et al., 2019;Gotelli & Stanton-Geddes, 2015).With some of these approaches, the distribution range of a species is first split into distinct entities (hereafter 'partitions') based either on evolutionary lineages or on phenotypic traits that are assumed to correlate with different responses to local climate conditions (Hällfors et al., 2016;Lecocq, Harpke, Rasmont, & Schweiger, 2019;Smith et al., 2019).
Local models are built separately for each partition and their predictions are combined across the entire range of the species (D'Amen, Zimmermann, & Pearman, 2013;Ikeda et al., 2017;Pearman, D'Amen, Graham, Thuiller, & Zimmermann, 2010).These studies indicate that such local approaches can significantly alter model outcomes (Smith et al., 2019;Theodoridis, Patsiou, Randin, & Conti, 2018).However, detailed information on relevant intraspecific differentiation is available for few species only.Hence, these methods are currently hardly applicable to extensive datasets of species (Moran et al., 2016).From these studies, it is therefore difficult to derive general recommendations on the type of species or conditions under which niche-based models should account for intraspecific spatial variation in sensitivity and response to climate conditions (Hällfors et al., 2016).
The distribution of several species can be partitioned simultaneously according to pre-defined ecological regions (Estrada-Peña & Thuiller, 2008) or regular-shaped geographical entities (Osborne & Suárez-Seoane, 2002).Then local models can be built to capture the response of the species to climate conditions at the level of each partition.These approaches do not require prior knowledge of the processes, lineages or traits that (are assumed to) affect different responses to local climate conditions.Yet, such pre-established partitions are not necessarily distributed well along the gradient of climate conditions used by the species.Climate conditions play, however, an important role in shaping local adaptation or in structuring the spatial distribution of ecological traits that affect intraspecific variation in the response of organisms to changing climate conditions (Amburgey et al., 2018;Ley et al., 2018;Macdonald et al., 2018;Moran et al., 2016;Nice et al., 2019).Local approaches applicable to a large sample of under different dispersal assumptions.Using European butterflies as an example, we show that this approach can be used to explore the uncertainty about predicted distribution shifts arising from intraspecific spatial variation in sensitivity and response to changing climate conditions.
Main conclusions: Our modelling approach is not intended to replace advanced modelling methods based on species-specific knowledge of ecological and evolutionary processes, but it is useful as an exploratory tool to detect species for which detailed information on intraspecific responses to climate conditions is likely to make a difference for prediction of future distribution shifts.

K E Y W O R D S
butterflies, climate change, data partitioning, dispersal, ecological niche models, local adaptation, model aggregation, niche dynamics, range shift, species distribution models species are therefore needed that split the individual distribution of each of them into partitions to capture the diversity of idiosyncratic climate conditions used by the species throughout their ranges.
Here, we developed such a modelling framework that accommodates potential spatial variation in the response of species to changing climate conditions and that can be applied for a large number of species and within different dispersal constraints to predict their potential future distribution shifts.The proposed method does not require any detailed understanding of the processes underlying intraspecific variation, as it only relies on species records and climate data to partition the geographical distribution of the species and then uses these partitions as proxies for such variation.We compared the predicted distribution shifts obtained from our proposed method to the predictions derived from a standard niche-based modelling approach assuming a homogeneous response of the species to climate conditions.This comparison was carried out for a large set of butterfly species coping with varying climate conditions throughout their range in Europe.With this example, we tested the usefulness of our proposed framework to explore the uncertainty around predicted species distribution shifts that arises from intraspecific variation in sensitivity and response to changing climate conditions.

| Species data
We used a large database of European butterflies because the geographical ranges of these species have been shown to shift in response to climate change (Devictor et al., 2012;Parmesan et al., 1999) and their populations vary in sensitivity to temperature and precipitation and show heterogeneous responses to climatic variation across their ranges (Mills et al., 2017;Nice et al., 2019).
Butterfly distribution data were compiled from the 'Mapping European Butterflies' project (Kudrna, 2002;Kudrna et al., 2011;Settele et al., 2008).Data from 1980-2000 were extracted for all species in Europe (N = 427) using the European Terrestrial Reference System Lambert Azimuthal Equal Area (ETRS89-LAEA) at 50-km resolution.Species with fewer than 30 occupied grid cells in the study area (n = 142) were discarded from further analyses.To ensure that we used the whole diversity of European butterfly distribution patterns when applying our proposed method, butterfly species were classified into three categories reflecting the geographical extent of their distribution range (see Supporting Information Appendix S1 and Figure S1.1):(a) narrow-range species, (b) medium-range species and (c) widerange species.We randomly selected 40 species within each category (Supporting Information Table S1.1) as a trade-off between covering a significant proportion of species and limiting computational time.

| Climatic data
Climate data were obtained from the Climatic Research Unit (CRU) climatological database with 10-minute spatial resolution pixels (Mitchell, Carter, Jones, Hulme, & New, 2004).We derived the three most commonly used climate variables reflecting the main constraints on butterfly growth and survival (Titeux et al., 2017): annual daily temperature sum above 5°C (growing degree-days), mean temperature of the coldest month and annual water balance (as estimated in Hickler et al., 2009).To rescale the climate variables to the resolution of the species data, we intersected the 10-minute spatial resolution climate pixels of the CRU dataset with the 50-km grid cells and we calculated for each of the grid cells the area-weighted average of the climate variables among the intersecting pixels (Martin, Van Dyck, Dendoncker, & Titeux, 2013).

| Partitioning of species distributions
Species coping with diverse climate conditions across their ranges are likely to show a high level of spatial variation in sensitivity, tolerance and response to these conditions (Fitzpatrick & Keller, 2015).
The geographical distribution of each species (i.e. the occupied grid cells) was split into a number of partitions to mirror this intraspecific spatial variation.Partitions were produced so that each of them included populations experiencing idiosyncratic climate conditions within the full range of conditions occupied by the species across the study area.The Ward's minimum variance agglomerative clustering method implemented in the R package 'const.clust'(Legendre, 2011) was applied to create the partitions.This algorithm used the multivariate dissimilarity matrix of Euclidean distances between grid cells based on the three standardized climate variables reflecting climate conditions during the period 1971-2000 (Legendre & Legendre, 2012).
We implemented four different partitioning procedures that differed in input data and spatial constraints used in the clustering algorithm: (a) generic and unconstrained partitioning, (b) generic and constrained partitioning, (c) specific and unconstrained partitioning, and (d) specific and constrained partitioning.With the 'generic partitioning' procedures, we disassembled the whole set of grid cells at once to produce climatic strata irrespective of the distribution of the species (Figure 1).To partition the individual distribution of a species, the occupied grid cells were then intersected with the obtained climatic strata (Figure 2a,b).With the 'specific partitioning' procedures, the distribution of each species was partitioned in a species-specific way using only the set of grid cells occupied by the species as input data in the clustering algorithm (Figure 2c,d).The generic and specific partitioning procedures were applied with spatial constraints in the clustering algorithm ('constrained partitioning', Figures 1b and 2b,d) and without these constraints ('unconstrained partitioning', Figures 1a and 2a,c).With spatial constraints, the grid cells included in each partition were forced to be spatially contiguous in addition to sharing similar climate conditions.
For each procedure and for each species, the number of partitions was set as to minimize the cross-validated residual error estimated based on 100 cross-validation iterations (Terrier, Girardin, Périé, Legendre, & Bergeron, 2013).Automated reallocation was applied to avoid the creation of partitions with an insufficient number of occupied grid cells to build the models (Hernandez, Graham, Master, & Albert, 2006): each partition with fewer than 15 occupied grid cells was merged in an iterative manner with its climatically closest partition according to the Mahalanobis distance between the different partitions (Legendre & Legendre, 2012).

| Niche-based models
The main idea behind our proposed modelling approach was to combine across the study area the outcomes of a set of niche-based models (hereafter 'local models') built separately with the different partitions of occupied grid cells.As we implemented four procedures to partition the distribution of each species (Figure 2), four different sets of local models were built.For each of these sets separately, the outcomes of individual models were assembled at the European scale.This local modelling approach aimed at capturing the full range of climate conditions used by the species throughout the study area.For the purpose of comparison (see Supporting Information Appendix S2: Figure S2.2), we also implemented a second approach using the entire set of grid cells occupied by the species to build a single model (hereafter 'global model') irrespective of the partitions created above.This standard modelling approach assumed and intended to capture a homogeneous response of the species to climate conditions across their ranges.
We used Maxent version 3.3.3e(Phillips & Dudík, 2008) to build the local and global models and estimate variation in climatic suitability for each species across the study area.Maxent has been shown to perform well in comparison with other presence-only modelling methods, especially with few species records (Hernandez et al., 2006).To limit model overfitting when dealing with small sample sizes, we used only linear and quadratic features of the climate variables (Hernandez et al., 2006).Each model was built 10 times for each species based on a random sample of 70% of occupied grid cells for model training and 30% for model evaluation.The random selection was applied to the set of occupied cells included in each partition for the local models and to the whole set of occupied grid cells for the global models.The whole study area was used as background data (number of grid cells: N = 2,643) for training the models (Oney, Reineking, O'Neill, & Kreyling, 2013) as described in detail in Supporting Information Appendix S2.
The predictions of local and global models were used to estimate the current (2000) and future (2020, 2050 and 2080 under BAMBU, GRAS and SEDG scenarios) distributions of the species.The continuous logistic predictions reflected the climatic suitability for the species (or for some populations of the species in the case of local models) within each grid cell.These predictions were transformed into binary values, that is, suitable or unsuitable climate conditions, using the maximum sum of sensitivity and specificity in the training dataset as a threshold (Liu, White, & Newell, 2013).To estimate the current distribution of a species, we considered that suitable and unsuitable climate conditions were occupied and unoccupied by the species, respectively.For the global models, grid cells where climate conditions were currently estimated as suitable were therefore used to predict the presence of the species in 2000.For the F I G U R E 1 European climatic strata obtained from the clustering algorithm applied to the whole set of grid cells in Europe (a) without spatial constraints (12 strata) and (b) with spatial constraints (10 strata).Grid cells were clustered based on the same three climate variables as the ones used to build the niche-based models.These European climatic strata were then used as a basis to partition the grid cells occupied by the different butterfly species with the 'generic partitioning procedure' (see Figure 2a,b).Colours were selected arbitrarily to differentiate between the climatic strata local models, we produced an assembled prediction for each of the four partitioning procedures separately.Here, predicted presences from the different local models were combined across the study area assuming that a species was present in a grid cell when at least one local model predicted that the cell was suitable (D'Amen et al., 2013; see Supporting Information Appendix S2).Before this assemblage, we discarded the predictions from any local model considered as potentially unreliable, that is, when it predicted the absence of the species in more than 30% of grid cells from the test dataset where the species was actually present (Radosavljevic & Anderson, 2014).
We applied a similar approach to predict the future distribution of the species, but here we assumed that the species (or the populations of the species) was (were) either not limited when colonizing grid cells with suitable climate conditions (hereafter 'full dispersal') or unable to colonize such newly suitable areas (hereafter 'no dispersal') (Bateman, Murphy, Reside, Mokany, & VanDerWal, 2013;Valladares et al., 2014).Modelling performance was evaluated based on the predictions of the global model and on the assembled predictions of the local models.We used the omission rate (proportion of occupied grid cells in the test dataset with predicted species absence) and the current predicted area of occupancy (PAO 2000 -proportion of grid cells with predicted species presence relative to the total number of grid cells in the study area) as measures of modelling performance (Phillips & Dudík, 2008).

| Comparison between local and global models
We examined whether the modelling performance and the predictions of current and future species distributions differed between F I G U R E 2 Partitions of grid cells occupied by the small tortoiseshell (Aglais urticae) obtained from the four partitioning procedures and used to build the local models for that example species.In the 'generic' procedure, the whole set of grid cells occupied by the species was overlaid upon the European climatic strata (Figure 1a For each species under current (2000) and future (2020, 2050 and 2080) climate conditions and according to the different scenarios and dispersal hypotheses, we estimated the extent to which the assembled predictions of the four sets of local models deviated from the predictions of the global model.We calculated the ratio between the predicted areas of occupancy (PAOs) obtained from the local and from the global models.We then expressed this ratio on a logarithmic scale to document dissimilarities in predictions between our proposed method and the standard approach in a symmetrical way.
Variation in this PAO log-ratio was related to a number of categorical variables and their interactions in ANOVA analyses carried out for current and future predictions separately.For current predictions, the categorical variables were partitioning procedure (i.e. the four types of partitions used to build the local models) and category of species range.For future predictions under the 'full dispersal' and 'no dispersal' hypotheses, we also included scenario (SEDG, BAMBU, GRAS) and period (2020,2050,2080).
Finally, we plotted the response curves for the different species to assess how the climatic suitability estimated by both types of models varied along the gradient of current climate conditions (see Supporting Information Appendix S3) and we mapped the current (2000) and future (2080) distribution of selected species and species richness as predicted by both types of models (see Supporting Information Appendix S4).

| Partitioning of species distributions
The generic unconstrained and constrained partitioning procedures created 12 and 10 European climatic strata in Europe, respectively (Figure 1a,b).The grid cells occupied by each species were then intersected with these strata to split the distribution of each species into partitions sharing similar climate conditions (Figure 2a,b).In contrast, the specific unconstrained and constrained partitioning procedures produced partitions of occupied grid cells that were specific to each species (Figure 2c,d).The number of partitions increased significantly from narrow-to wide-range species (F 2;117 = 341.8;p < .001)and differed significantly between the four partitioning procedures (F 3;116 = 41.1;p < .001).A higher number of partitions were created with the specific partitioning procedures than with the generic procedures (Table 1).

| Modelling performance
Two local models (out of 2,899 in total) were identified as not sufficiently robust and were removed before assembling the predicted distributions of the species.Omission rate and current predicted area of occupancy (PAO 2000 ) varied significantly among models and categories of species range (Table 2).Omission rate was significantly larger for predictions based on global models than for those obtained from the assemblages of local models (Figure 3a).PAO 2000 was smaller for predictions derived from the global model than for those based on the assemblages of local models (Figure 3b).The difference in omission rate and in PAO 2000 between the two types of models increased gradually from narrow-to wide-range species (Figure 3, Table 2).
Post-hoc analyses for multiple comparisons with Bonferroni correction showed that the differences in omission rate and PAO 2000 between the global model and each assemblage of local models were significant for each combination of partitioning procedure and category of species range, except in the case of PAO 2000 for narrow-range species (Figure 3).Omission rate and PAO 2000 were also similar across the different assemblages of local models derived from the four partitioning procedures (p > .05for all pairwise comparisons).

| Estimated climatic suitability
For medium-and, particularly, wide-range species, the two types of models produced estimations of currently suitable conditions that TA B L E 1 Minimum, mean (with standard deviation) and maximum number of partitions of grid cells occupied by the butterfly species resulting from the different partitioning procedures for each category of species range in Europe  3).In contrast, local models covered these parts of the distributions more thoroughly and, therefore,

| Predicted shifts in species distributions
The dissimilarity between the predictions of the two types of models (PAO log-ratio) varied strongly between the different categories of species range and also with time (Table 3, Figure 4).There was a significant but smaller effect of the climate change scenario, with the most severe scenario (GRAS) producing the largest differences between the two model types.PAO log-ratio differed among the four assemblages of local models and this difference was particularly marked for the 'no dispersal' hypothesis, indicating that the assembled predictions obtained from the four partitioning procedures deviated from the predictions of the global model in different ways, especially when the species was assumed to be unable to track suitable conditions.
Variation in PAO log-ratio from 2000 to 2080 was larger for widely distributed species than for geographically restricted species (Figure 4).PAO log-ratio was close to 0 for narrow-range species in 2000 and did not vary markedly with time under any of the two dispersal hypotheses, indicating a general agreement between the two types of models across time.Higher PAO log-ratio in medium-range species indicated that local models predicted a larger PAO than global models, but this difference decreased with time, especially under the 'no dispersal' hypothesis.In wide-range species, PAO log-ratio slightly increased with time (especially between 2050 and 2080) under the 'full dispersal' hypothesis, but strongly decreased under the 'no dispersal' hypothesis to become close to 0 and even lower with some partitioning procedures, indicating that the local models predicted a smaller future area of occupancy than the global models.
Large between-model differences were observed in future spatial predictions towards the marginal parts of the distributions of

| D ISCUSS I ON
Our analytical framework disassembles the distribution of species along the occupied gradient of climate conditions throughout their distribution range, and reassembles the predictions obtained from individual models based on each of these climatically distinct partitions.This approach assumes that the response of species to climate conditions may vary spatially.Using the partitions as proxies for such intraspecific variation, our method aimed at capturing the whole variation of climate conditions used by the species across their distribution ranges.Such local models contrast with standard approaches that assume a homogeneous species response to climate conditions across the range.The proposed approach does not rely on any detailed understanding of the ecological and evolutionary processes underlying intraspecific spatial niche variation.Hence, it can be applied to a large sample of species and, as we will explain below, it can serve as a basis for exploring the level of uncertainty arising from intraspecific variation when predicting future species distributions.
Overall, the predictions obtained from our proposed method deviated significantly from those based on a standard approach.A first important difference between the two types of models resulted from the fact that assembling the outcomes of local models covered more systematically the diversity of idiosyncratic climate conditions used by each species throughout its distribution range.This difference gradually increased with the geographical extent of the range of the species (see also Pearman et al., 2010).Local and global modelling approaches predicted similar areas of suitable climate conditions for species associated with geographically restricted ranges, whereas significantly different predictions were obtained for widely distributed species occupying a variety of climate conditions.The geographical distributions of wide-range species were split into a larger number of partitions than those of narrow-range species.For that reason, the assemblage of local models for wide-range species captured the diversity of occupied climate conditions more thoroughly than a single model built across the entire study area.
The response curves estimated by the two types of models showed that local models tended to place a stronger emphasis than global models on climate conditions used by widespread species at the marginal areas of their distributions far from the core of the

TA B L E 3
Results of the analysis of variance for linear model (ANOVA) testing for differences in predicted area of occupancy (PAO) log-ratio between the different assemblages of local models derived from the four partitioning procedures, the categories of butterfly species range in Europe and the interaction between the two factors for the current and future predictions.PAO log-ratio reflected similarity between the predictions from the global model and those from the assemblages of local models.For future predictions, the scenario, the period and all resulting interactions between factors were also included in the analyses ranges.Some of the largest differences in climatic suitability between the two types of models were predicted on latitudinal or altitudinal extremes of the gradient of climate conditions in Europe.
Standard modelling approaches may at least partly neglect marginal populations of widespread species when assessing the response curves and estimating spatial variation in climatic suitability.This leads to concern because there is growing evidence for the biological significance of specific population responses at the distribution margins (Diniz-Filho et al., 2019;Pironon, Villellas, Morris, Doak, & García, 2015).The local approach implemented here may, however, overemphasize the importance of these marginal populations.We assumed that a grid cell was climatically suitable for a species if at least one of the local models predicted so (D'Amen et al., 2013).
Here, the same importance was given to the predictions derived from local models built with a few occupied grid cells in marginal climate conditions as to the predictions obtained from models generated with a larger number of occupied grid cells near the core of the distribution.Therefore, the assemblage of local models may potentially overemphasize the presence of widespread species in marginally used climate conditions.Yet, we opted for this approach to represent as much as possible the variety of local climate conditions The generic constrained partitioning procedure was used here for the local modelling approach.The observed frequency distribution and the response curves of the species are shown along the gradient of current climate conditions in Europe, that is, the first axis of a principal component analysis (PCA) based on the three climate variables used in the models, divided into 10 bins each including the same number of grid cells.The observed frequency distributions shown as histograms represent the proportion of grid cells within each bin along the PCA axis that are occupied by the species either in Europe (light grey histograms) or within each individual partition used to build the local models (dark grey histograms).The response curves of the species represent how the proportion of grid cells predicted as suitable for the species (or for some populations of the species in the case of local models) varies along the PCA axis.These curves reflect the predicted climatic suitability for the species (or its local populations) along the gradient of current climate conditions.Supporting Information Appendix S3 provides detailed information on the estimation of these response curves that further work is needed here and we opted for quantifying and mapping the spatial agreement and mismatch between the climatically suitable areas predicted from the two types of models as a way to assess the level of uncertainty associated with intraspecific variation (Figure 5 and Supporting Information Figures S4.6-S4.14).
The second important difference between the two types of models is related to the hypothesis made when converting future climatically suitable areas into predicted presences.When assuming that the local populations of a species will not be able to colonize newly suitable areas, its geographical distribution was predicted to become increasingly fragmented even at the core of the range.Under this 'no dispersal' hypothesis, our proposed approach reflected a situation in which populations are locally adapted to idiosyncratic climate conditions and will go extinct if local conditions become unsuitable because they are not able to move and colonize areas that will become suitable.
Here, the predictions of distribution shifts for widespread species under climate change turned out to be more pessimistic than predictions obtained from the standard approach (PAO log-ratio lower than 0).Our results are consistent with those of Valladares et al. (2014) showing that species ranges may shrink considerably in the future due to the interaction between population differentiation and limited dispersal.Although dispersal constraints have often been integrated into niche-based modelling approaches, the potential interaction with intraspecific variation in climatic requirements has been less studied and warrants further research (Valladares et al., 2014).
The predicted distributions of species associated with very small geographical ranges were similar according to both types of modelling approaches.With the partitioning procedures, the distribution of each narrow-range species was split into a small number of partitions (Table 1) -sometimes it remained a single unit such as for the Frigga fritillary (Boloria frigga; Supporting Information Figure S4.8).For such species, the assembled predictions from a single (or a few) local model(s) were very similar to those of the global model (see Supporting Information Appendix S2 for detailed explanation).Our proposed local approach resulted in slightly lower omission rates and similar predicted areas of occupancy for narrow-range species.Hence, the approach did not systematically put a disproportionate emphasis on the climatic margins of the distribution of any species.It potentially may have done so only for widespread species coping with a large variety of climate conditions across their range and for which global models run the risk of neglecting the conditions that are marginally used when estimating a homogeneous response.Future distribution shifts predicted by the two types of models are also comparable for narrow-range species because applying the dispersal constraints at the intraspecific level -that is, for a small number of local populations -was roughly equivalent to applying these constraints at the specific level.
Instead of advocating the use of one approach instead of the other, we encourage the use of both types of models to explore the range of uncertainties in model development and predictions of future distribution shifts that arise from intraspecific spatial niche variation.Such model comparison can be done as a first-step approach for a large number of species (Hällfors et al., 2016) and serve as a basis for identifying cases where standard modelling approaches are unlikely to capture the diversity of climate conditions experienced by the species across their range.Local idiosyncrasies of climate conditions may explain variation in the functional response of butterflies to a changing climate through physiological or behavioural mechanisms (Nice et al., 2019).Our approach will not be able to capture these mechanisms, but it could be useful to detect and prioritize species for which experimental evidence on such key processes would be required and where advanced modelling methods based on these processes should be applied to better understand and predict how local adaptation, adaptive plasticity and dispersal abilities would interact with changing climate conditions to affect distribution shifts (Diniz-Filho et al., 2019;Peterson et al., 2019).
Most of the previous studies that incorporated population differentiation into the models suggested that intraspecific variation in response to climate conditions could buffer the future impacts of climate change (Oney et al., 2013;Pearman et al., 2010).Although our predictions at the margins of the distribution of widespread species are consistent with these suggestions, the integration of dispersal constraints also showed an opposite pattern if the different populations of the species are locally adapted but not able to disperse (Valladares et al., 2014).Therefore, the application of standard modelling approaches assuming that all populations of a species will respond in the same manner to changing climate conditions bears unappreciated risks for conservation decision-makers as they may lead to too pessimistic or too optimistic predictions on future species persistence.We provide a tool to filter species for which further information should be collected to document processes underlying potential local adaptation, adaptive plasticity and dispersal abilities.It will contribute to increasing our ability to predict future shifts in species distributions and their resilience under climate change (Peterson et al., 2019;Smith et al., 2019).

ACK N OWLED G M ENTS
Most of the European distributional data were compiled by Dr Kudrna, especially within the framework of Kudrna (2002) and Kudrna et al. (2011) Under the 'full dispersal' hypothesis, suitable grid cells in 2020, 2050 and 2080 were therefore considered as future predicted presences irrespective of the predictions in 2000.Under the 'no dispersal' hypothesis, suitable cells in 2020, 2050 and 2080 were considered as predicted presence only if conditions were also suitable in 2000.This dispersal constraint was applied at the specific level (based on the predictions of the global models) or at the intraspecific level (based on the predictions of each local model individually).
,b) to obtain the different partitions (a) without spatial constraints and (b) with spatial constraints.In the 'specific' procedures, the clustering algorithm was directly applied to the set of grid cells occupied by the focal species (c) without spatial constraints and (d) with spatial constraints.Colours in (a) and (b) match those of Figure 1a,b, respectively, and colours in (c) and (d) were selected arbitrarily to differentiate between the partitions the global model and the four assemblages of local models.We also tested whether the differences were related to the geographical extent of the range of each species.We carried out an analysis of variance for linear models (ANOVA) to test for the effect of the following two categorical variables (and their interactions) on omission rate and predicted area of occupancy (PAO 2000 ): model (global model and each of the four assemblages of local models) and category of species range (narrow-, medium-and wide-range species).
predicted a larger area of occupancy under current climate conditions (PAO 2000 ).For widely distributed species, the average response curves between the two types of models often closely matched the observed frequency distribution along the gradient of climate conditions (see Supporting Information Appendix S3).

F
Modelling performance measures based on the global model (dark grey) and on each of the four assemblages of local models (light grey): (a) omission rate in the test dataset and (b) predicted area of occupancy in 2000 (PAO 2000 -percentage relative to the geographical extent of Europe).Modelling performance is shown using boxplots for the random sample of butterfly species within each category of species range (narrow-, medium-and wide-range species) (see Supporting Information Appendix S1).Asterisks indicate the level of significance for the differences in modelling performance between the global model and each assemblage of local models for the different categories of species range (see that are potentially suitable for a species throughout its range.With the difference between the predictions of the two types of models, we assess the level of uncertainty arising from intraspecific variation in response to climate conditions.Interestingly, our results for most widespread species show that the average response curves reflecting the predicted climatic suitability between the global and local models closely match the observed frequency distribution of the species along the gradient of climate conditions.It is important to note, however, that this intermediate predicted suitability results from averaging the response curves of the global and local models; it is not based on the combined predictions of the two types of models at the level of each individual grid cell.Combining the predictions of these two conceptually different types of models and mapping the resulting predicted distribution of the species is warranted to evaluate if this hybrid approach would produce a better balance between omission rate and predicted area of occupancy than either of the two individual approaches.Implementing such an approach is, however, not straightforward, because the predictions derived from these two types of models are not directly comparable and easy to combine with each other (see Supporting Information Appendix S3).We acknowledge F I G U R E 4 Degree of agreement between the global model and each of the four assemblages of local models [predicted area of occupancy (PAO) log-ratio] on their current predictions for 2000 (white) and on their future predictions for 2020 (light grey), 2050 (medium grey) and 2080 (dark grey) under the most severe climate change scenario ('Growth Applied Strategy' -GRAS) in Europe and according to the (a) 'full dispersal' hypothesis or the (b) 'no dispersal' hypothesis.PAO log-ratio is the ratio between the predicted areas of occupancy obtained from the local models (PAO local) and from the global model (PAO global) expressed on a logarithmic scale: values higher (lower) than 0 indicate that local models predicted a larger (smaller) area of occupancy than global models.PAO log-ratios are shown using boxplots for the random sample of species within each category of butterfly species range (narrow-, medium-and wide-range species) (see Supporting Information Appendix S1) F I G U R E 5 Observed (black) and predicted butterfly species distributions in 2000 and in 2020, 2050 and 2080 ['Growth Applied Strategy' (GRAS) scenario] under the 'full dispersal' and 'no dispersal' hypotheses.Areas where the species was predicted to be present by the two types of models are shown in green and presences predicted exclusively by the local models or the global model are shown in yellow and blue, respectively.The different steps of the modelling approaches to estimate suitable climate conditions are shown for a wide-range species (light grey background), the small tortoiseshell (Aglais urticae).The same procedure was applied to narrow-range species (dark grey background) and medium-range species (intermediate grey background) (see Supporting Information Appendix S4 for other examples).
Results of the analysis of variance for linear model (ANOVA) testing for differences in omission rate and current predicted area of occupancy across Europe (PAO 2000 ) between models (global model and each assemblage of local models), categories of butterfly species range and their interaction TA B L E 2

Table 2
; ANOVA and post hoc Bonferroni tests for multiple pairwise comparisons; ***p < .001.*p < .05)wide-range species under the two dispersal hypotheses (Figure 5, , and their use for scientific analyses is agreed on the basis of a contract between Dr Kudrna and the Helmholtz Centre for Environmental Research -UFZ, Germany.We used the scenarios developed in the 'ALARM' project (Framework Programme FP6 of the European Commission, grant number GOCE-CT-2003-506675). Financial support was obtained from the Luxembourg National Research Fund (FNR-AFR PHD-09-121), Fédération Wallonie-Bruxelles and Université catholique de Louvain (ARC-grant 17/22-086) and European Commission (contract no.308454; FP7-ENV-2012).This is publication BRC355 of the Biodiversity Research Centre (Earth and Life Institute, UCLouvain, Belgium).