# Integrating intraspecific differentiation in species distribution models: Consequences on projections of current and future climatically suitable areas of species

Thomas Lecocq, Alexander Harpke, Pierre Rasmont, Oliver Schweiger

► To cite this version:

Thomas Lecocq, Alexander Harpke, Pierre Rasmont, Oliver Schweiger. Integrating intraspecific differentiation in species distribution models: Consequences on projections of current and future climatically suitable areas of species. Diversity and Distributions, inPress, 25 (7), pp.1088-1100. 10.1111/ddi.12916 . hal-02120201

HAL Id: hal-02120201

https://hal.univ-lorraine.fr/hal-02120201

Submitted on 5 May 2019

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Integrating intraspecific differentiation in species distribution models: Consequences on projections of current and future climatically suitable areas of species

Thomas Lecocq<fnref n="1" />,<fnref n="2" /> ID

Alexander Harpke<fnref n="3" />

Pierre Rasmont<fnref n="2" />

Oliver Schweiger<fnref n="3" /> ID

**Correspondence**

Thomas Lecocq, Université de Lorraine, Inra, URAFPA, Nancy, France.
Email: thomas.lecocq@univ-lorraine.fr

**Funding information**

European Community's Seventh Framework Programme, Grant/Award Number: FP7/2007-2013

Editor: Janine Bolliger

**Abstract**

**Aim:** Conventional species distribution models (SDMs) usually focus on the species level but disregard intraspecific variability. Phylogeographic structure and evolutionary significant units (ESU) have been proposed as pragmatic proxies to incorporate intraspecific differentiation in SDMs. Nevertheless, the efficiency of using these proxies in SDMs has been poorly investigated. We analysed how the projections of current and future climatically suitable areas can be affected when using ESU-based or lineage-based models compared to a species-level model.

**Location:** West-Palaearctic region.

**Methods:** As examples, we used three bumblebee species (133,787 observations). We assessed potential climatic niche differentiation between species level, lineages and ESUs, by determining the niche position and niche breadth for each classification level by outlying mean index analyses. Subsequently, we developed SDMs for each species and classification level separately using boosted regression trees prior to a comparison of their performances. Finally, we used the alternative models to project the extent of climatically suitable areas in 2070.

**Results:** We found that in spite of highly similar overall model accuracy, integrating intraspecific variability significantly increases model sensitivity (i.e., better predicting presences) while decreased model specificity (i.e., over-predicting the range). Consequently, both predictions of current and projections of future suitable conditions differed among the three approaches.

**Main conclusions:** We showed that although integrating lineage or ESU information did not improve the accuracy of conventional species-level SDMs, it led to considerably different conclusions. As SDM-based climatic risk assessments are increasingly used to help and improve conservation plans, divergences and limitations of each modelling approach should be taken into account for developing efficient biodiversity management strategies. Lineage and ESU-based SDMs offer the advantage to draw attention to species in which allopatric populations could display physiologically

<footnote n="1">¹Université de Lorraine, Inra, URAFPA, Nancy, France</footnote>

<footnote n="2">²Laboratory of Zoology, Research Institute of Biosciences, University of Mons, Mons, Belgium</footnote>

<footnote n="3">³Department of Community Ecology, Helmholtz Centre for Environmental Research - UFZ, Halle, Germany</footnote>
different responses to climate change when they lead to different results than spe-
cies-based models.

KEYWORDS

bumblebee, climate change, evolutionary significant unit, intraspecific variability, model performance, species distribution model

**1** | INTRODUCTION

Species distribution modelling (SDM) is among the most widely
used methods in macro-ecology and conservation ecology (Mainali
et al., 2015; Zurell et al., 2016). Several methods are available but
most studies use a correlative approach where the species niche is
modelled by relating distribution data to environmental conditions
(Guisan & Zimmermann, 2000). SDMs are typically applied to map
species ranges and to project their changes under future global
change scenarios (Guisan & Zimmermann, 2000; Rasmont et al.,
2015). However, most of the SDM studies consider species as a unit,
while intraspecific differentiation and potential niche divergences
have largely been ignored (D'Amen, Zimmermann, & Pearman,
2013; Pearman, D'Amen, Graham, Thuiller, & Zimmermann, 2010;
Valladares et al., 2014). This can bias models especially when smaller
populations with potential local adaptations and individualistic re-
sponses are less-well represented by species-level models (D'Amen
et al., 2013; Lecocq, Rasmont, Harpke, & Schweiger, 2016; Pearman
et al., 2010). This is of particular concern in the context of climate
change and respective mitigation options, especially when the re-
sulting models underestimate or overestimate the climatic tolerance
or sensitivity of species (Benito Garzón, Alía, Robson, & Zavala,
2011; Lecocq, Rasmont et al., 2016).

One way to lower such uncertainties of SDMs could be to inte-
grate intraspecific niche divergences in the models (Valladares et al.,
2014). However, this requires a priori definition of subunits differ-
entiated within a species of concern (Lecocq, Rasmont et al., 2016).
Defining such subunits on the basis of features shaped by adaptation
to regional climatic conditions (e.g., higher heat stress resistance for
some populations; Martinet, Lecocq, Smet, & Rasmont, 2015) is the
best solution (Alexander & Edwards, 2010; Lecocq, Rasmont et al.,
2016). However, comprehensive information on such features is sel-
dom available for the majority of species. A more pragmatic approach
may use proxies for the detection of intraspecific differentiation
(Lecocq, Brasero, Meulemeester et al., 2015b) that could be related
to different climatic requirements of populations (Lynch, 1996). This
relies on the assumption that intraspecific differentiation triggered
by demographic history and/or local adaptations shapes population
specificities in genetic, morphology, physiology, behaviour, or ecol-
ogy that can impact the environmental requirements (Hewitt, 2004;
Lecocq, Rasmont et al., 2016; Lynch, 1996).

Considering the phylogeographic structure has been sug-
gested as an effective way to integrate intraspecific differentia-
tion into SDMs (i.e., lineage-based models) (D'Amen et al., 2013;

Marcer, Méndez-Vigo, Alonso-Blanco, & Picó, 2016; Pearman et al., 2010). The phylogeographic structure is usually estimated with neutral markers which are impacted by the heritage of the populations resulting from the effects of neutral evolutionary forces (i.e., genetic drift, mutation, migration, hybridization; Avise, 2000). This allows defining differentiated allopatric lineages with individualistic evolutionary history (Avise, 2000) which has possibly shaped lineage-specific climatic requirements (Marcer et al., 2016). However, the definition of units below the species level based on solely neutral marker remains controversial because (a) DNA sequences analysed are often chosen arbitrarily (Cruaud et al., 2014) and (b) ecologically differentiated populations are not always characterized by the accumulation of many genetic differences (Ferguson, 2002; Salvato et al., 2002). These limitations can be overcome by combining different lines of evidence such as phylogeographic structure based on multi-markers along with phenotypic traits and/or ecological features within the context of the integrative taxonomy (Schlick-Steiner et al., 2010). This approach allows detecting conspicuous allopatric differentiated populations (Lecocq, Brasero, Meulemeester et al., 2015b) that are closely related to evolutionarily significant units (ESUs). Initially, the ESUs were defined as populations with current geographic separation, neutral marker genetic differentiation or locally adapted phenotypic traits (Conner & Hartl, 2004). However, this initial definition has been criticized because (a) current geographic disjunction can be triggered by very recent demographic events making adaptive differences between allopatric populations unlikely (e.g., see examples of recent local extinctions of bumblebees in Rasmont et al., 2015), (b) ESUs defined solely using neutral markers ignore adaptive differences (Frankham, Ballou, & Briscoe, 2010), (c) trait divergences considered as locally adapted phenotypic traits can be shaped by phenotypic plasticity (i.e., the ability of a genotype to produce different phenotypes under different environmental or developmental conditions; Valladares et al., 2014). Therefore, we here use a more restrictive definition: ESUs are geographically isolated conspicuous groups differentiated in genetic and phenotypic traits (Frankham et al., 2010). These ESUs are used as an effective shortcut for estimating patterns of intraspecific diversity in conservation biology (Frankham et al., 2010). Since the integrative taxonomy approach is increasingly applied (e.g., arthropods: Hendrixson, Guice, & Bond, 2015; vertebrates: Costa & Amorim, 2014), corresponding information on ESUs is constantly growing for an increasing number of species. This could pave the way to a massive development of ESU-based SDMs and an improved
projection of tomorrow's biodiversity. Nevertheless, the efficiency of ESU-based SDMs compared to that of lineage-based SDMs or specie-level SDMs has been less-well investigated to date.

In this paper, we analyse how the projections of current and future climatically suitable areas can be affected by using ESU-based or lineage-based SDMs compared to a species-level model. As examples, we use three widespread bumblebee species (*Bombus lapidarius*, *B. pascuorum* and *B. terrestris*) with comparable distributions across Europe, for which phylogeography and ESU delimitations are available (Lecocq, Brasero, Martinet, Valterovà, & Rasmont, 2015; Lecocq, Coppée et al., 2016; Lecocq, Dellicour et al., 2015).

## 2 | METHODS

### 2.1 | Studied species and geographic range

We focused on three (West)-Palaearctic bumblebee species: *Bombus lapidarius*, *B. pascuorum* and *B. terrestris*. These species are some of the most abundant bumblebees where they occur (Rasmont et al., 2015). They are sympatric in most of their range and not threatened (i.e., IUCN Red List status: Least Concern with stable populations) (Rasmont et al., 2015). The occurrence data came from the database "Base de données fauniques de Gembloux-Mons" (Rasmont et al., 2015). We extracted all species occurrences between 1950 and 2014 in the West-Palaearctic region but because of the rather low sampling effort and unavailability of lineage/ESU statuses, we had to exclude Russia, Belorussia, Kazakhstan and Ukraine. We obtained 35,252 observations for

*B. lapidarius*, 63,565 observations for *B. pascuorum* and 34,970 observations for *B. terrestris* and called the resulting datasets "raw observation datasets." We aggregated these data to presence/absence data at 10 arc-minutes resolution grid (i.e., = 0.17° or ~14.5 km at 40°N) to account for differences in local sampling effort and to obtain reliable absence data. Additionally, we also used a lower resolution (i.e., 30 arc-minutes) to further taking into account the sampling effort differences. To further increase the reliability of absence data, we considered only grid cells as valid for absence data when at least one bumblebee individual from one of the 69 *Bombus* species had been observed between 1950 and 2014 (based on 575,421 observations in the database). As (a) most of observations from the database "Base de données fauniques de Gembloux-Mons" were based on studies investigating the whole bumblebee fauna from a particular geographic area and (b) the three studied species were very abundant, we considered that such grid cells without *B. lapidarius*, *B. pascuorum* or *B. terrestris* observations were most likely true absences for these species.

### 2.2 | Phylogeographic and ESU framework

We chose these three species because their phylogeographic structure and the ESU delimitations were available (Lecocq, Brasero, Martinet et al., 2015; Lecocq, Coppée et al., 2016; Lecocq, Dellicour et al., 2015) (Table 1). These lineages/ESUs occurred in quite distinct biogeographic regions (Supporting Information Figure S1) and differed considerably in morphology, genetic, behaviour, hair length or phenology (Lecocq, Brasero, Martinet et al., 2015; Lecocq, Brasero, Meulemeester et al., 2015b; Lecocq, Coppée et al., 2016; Lecocq,

<table><caption>**TABLE 1** The component phylogeographic lineages and ESUs for three bumblebee species according to the literature. The numbers in brackets refer to the number of observations. Taxo. nomen. refers to the corresponding taxonomic names of the lineage and ESU (if available)</caption><thead><tr><th>Species</th><th>Lineages</th><th>ESUs</th><th>Taxo. nomen.</th><th>References</th></tr></thead><tbody><tr><td rowspan="6">B. lapidarius<br>(33078)</td><td rowspan="3">Main lineage (32125)</td><td>Lapidarius (30811)</td><td>B. l. lapidarius</td><td>Lecocq, Dellicour et al. (2015)</td></tr><tr><td>Atlanticus (29)</td><td>B. l. atlanticus</td><td>Lecocq, Dellicour et al. (2015)</td></tr><tr><td>Decipiens (1285)</td><td>B. l. decipiens</td><td>Lecocq, Dellicour et al. (2015)</td></tr><tr><td>South-eastern lineage (774)</td><td>SE Lapidarius (774)</td><td>-</td><td>Lecocq, Dellicour et al. (2015)</td></tr><tr><td>South Italian lineage (179)</td><td>SI Decipiens (179)</td><td>-</td><td>Lecocq, Dellicour et al. (2015)</td></tr><tr><td>South-western lineage (3995)</td><td>Main unit (55622)</td><td>several subspecies</td><td>Lecocq, Brasero, Meulemeester et al. (2015b)</td></tr><tr><td rowspan="6">B. pascuorum<br>(59774)</td><td rowspan="3">Main lineage (55779)</td><td>Siciliensis (53)</td><td>B. p. siciliensis</td><td>Lecocq, Brasero, Meulemeester et al. (2015b)</td></tr><tr><td>Rehbinderi (104)</td><td>B. p. rehbinderi</td><td>Lecocq, Brasero, Meulemeester et al. (2015b)</td></tr><tr><td>SW unit (3995)</td><td>several subspecies</td><td>Lecocq, Brasero, Meulemeester et al. (2015b)</td></tr><tr><td rowspan="7">Terrestris (31802)</td><td rowspan="7">Terrestris (31802)</td><td>Africanus (187)</td><td>B. t. africanus</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>Audax (11411)</td><td>B. t. audax</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>Calabricus (75)</td><td>B. t. calabricus</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>Canariensis (80)</td><td>B. t. canariensis</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>Dalmatinus (4061)</td><td>B. t. dalmatinus</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>Lusitanicus (2098)</td><td>B. t. lusitanicus</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>Sassaricus (87)</td><td>B. t. sassaricus</td><td>Lecocq, Coppée et al. (2016)</td></tr><tr><td>B. terrestris<br>(31802)</td><td>Terrestris (13803)</td><td>B. t. terrestris</td><td>Lecocq, Coppée et al. (2016)</td></tr></tbody></table>
Déllicour et al., 2015; Lecocq, Rasmont et al., 2016). These specifici-
ties make differences in the climatic requirements between lineages/
ESUs likely for all species. The differences in climatic requirements
had been demonstrated in *B. terrestris* with changes in hair length,
diapause condition and phenology (de Jonghe, 1986; Lecocq,
Rasmont et al., 2016; Velthuis & van Doorn, 2006). We used infra-
specific nomenclature according to the International Commission on
Zoological Nomenclature. When infra-specific nomenclature did not
correspond to phylogeographic lineage/ESU or was non-existent,
we named lineages and ESU according to their geographic ranges
(Table 1).

We attributed the lineage/ESU status to the observation points according to the following criteria. For the data points corresponding to the individuals used in the baseline phylogeographic and ESU delimitation studies (about 5% of the dataset), we used the status underlined by these studies (Lecocq, Brasero, Martinet et al., 2015; Lecocq, Coppée et al., 2016; Lecocq, Dellicour et al., 2015). When the observation was identified to subspecies level (i.e., data were available from previous studies that identified the subspecies of the observations; about 25% of the remaining dataset) and this subspecies status corresponded to a particular lineage/ESU exclusively, we attributed the corresponding lineage/ESU to the observation (see matches between lineages/ESUs and subspecies in Table 1). For the remaining data points, we assigned lineage/ESU status based on the geographic region of the observation when it was in an area where only one lineage/ESU occurs. Doubtful observations (i.e., occurrences in areas where several lineages/ESUs overlap and none of the abovementioned criteria were applicable) were removed from the datasets (about 10% of the dataset). To avoid a potential bias by considering grid cells with doubtful observations as absences, we excluded these grid cells entirely for model development. The resulting dataset was called "used observation datasets" and was used for niche determinations and for SDMs (Table 1).

**2.3 | Climatic data**

We obtained 19 bioclimatic variables from the WORLDCLIM database version 1.4 (downloaded 01 October 2017) (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) for current (1950–2000) climatic conditions at a 10 arc-minutes resolution grid. For the analyses based on 30 arc-minutes resolution, we aggregated the bioclimatic layers to the same resolution as the bumblebee observation dataset (a) by merging the 10 arc-minutes and the 30 arc-minutes resolution grid in ArcGis 10.2 (www.esri.com) and (b) by calculating the mean of each bioclimatic variable within each of the 30 arc-minutes grid cell in R 3.2.5. (R Core Team, 2017). To avoid statistical problems due to high levels of collinearity among bioclimatic variables, we selected the following relevant and least correlated variables (r < 0.6) using a complete linkage cluster analysis based on a Spearman's correlation matrix calculated for the 19 bioclimatic variables: annual mean temperature, temperature seasonality, maximum temperature of the warmest month, precipitation of the wettest month and precipitation seasonality.

**2.4 | Climatic niche comparisons**

We assessed potential climatic niche differentiation between species level, lineages and ESUs, by determining the niche position and niche breadth for each classification level via outlying mean index (OMI) analyses (Dolédec, Chessel, & Gimaret-Carpentier, 2000). This measure describes the average relative position of species in a multi-variate environmental space by measuring the distance between the mean habitat conditions used by a target group and the mean habitat conditions across the entire study area. The higher the OMI value for a particular group of organisms (i.e., species, lineage, or ESU) is, the more marginal and atypical the group's niche is compared to mean habitat conditions observed in the studied region. In contrast, low (i.e., close to zero) OMI values indicate groups occurring in typical habitats of a region. Significance of OMI was tested with a random-ization test based on 1,000 permutations. We estimated niche breadth per classification level and species by using the tolerance index which measures the amplitude in the distribution of each tar-get group along a climatic gradient (Dolédec et al., 2000; Thuiller, Lavorel, Midgley, Lavergne, & Rebelo, 2004). Low tolerance values indicate climate specialists, while high tolerance values indicate cli-mate generalists. We first ran a principal component analysis (PCA) on the total climatic niche available in the study area based on the five selected bioclimatic variables for each species. Then, the cli-matic space occupied by each group of each classification level was projected to the climatic space available. Further, we used a multi-variate analysis of variance (MANOVA) to test for significant differ-ences in the mean climatic niche values of all classification levels in (R Core Team, 2017). Finally, we performed a pairwise MANOVA to identify niche differentiation among all pairs of classification levels. The p-values were adjusted for multiple testing using Benjamini and Hochberg's correction.

**2.5 | SDM development and evaluation**

We developed SDMs for each species and classification level sepa-
rately using boosted regression trees (BRTs) (Elith, Leathwick, &
Hastie, 2008) in R (R-package gbm version 2.1.; Ridgeway, 2013)
assuming a binomial error structure and using a learning rate of
0.001, a tree complexity of three (lower complexities resulted in
worse models while higher complexities performed similarly), and
a bag fraction of 0.75. To ensure an overall prevalence of 0.5 and
thus comparable resulting occurrence probabilities, we weighted
the absences by the ratio of the number of presences on the num-
ber of absences (Maggini, Lehmann, Zimmermann, & Guisan, 2006).
The BRT analyses converged to the best fitting model for all clas-
sification levels for all test case species except two ESUs with a very
small geographic ranges at 30 arc-minute resolution: ESU *Siciliensis*
(*B. pascuorum*) and ESU *Atlanticus* (*B. lapidarius*). To compare the con-
sequences of an ESU-based model with other models at this resolu-
tion, we removed the observations corresponding to these units in
the observation data for each kind of model (same procedure than
for doubtful observations during lineage/ESU status attribution).
We merged individual lineage/ESU models (so-called summed models) by keeping the highest probability obtained among models for a particular grid cell. To transform the resulting occurrence probabilities into presence/absence maps, we identified a threshold by maximizing the true skill statistic (TSS) (Allouche, Tsoar, & Kadmon, 2006).

Spatial autocorrelation in model residuals is known to result in biased parameter estimates and the inflation of type I errors in SDM (Crase, Liedloff, Vesk, Fukuda, & Wintle, 2014). Therefore, we used correlograms plot Moran's I (i.e., a measure of spatial autocorrelation; Moran, 1950) between grid cells as a function of the distance between them (Kissling & Carl, 2007) to compare autocorrelation in the model residuals in R (R-package ncf version 1.2-6; Bjornstad & Cai, 2018). Distance classes of geographic distances were created following Sturge's rule. The mean distance per distance class was used to draw correlograms. Moran's I statistics of individual distance classes were tested for significance by a resampling procedure (1,000 permutations). The global significance of the entire correlograms was evaluated using the progressive Bonferroni technique by dividing the significance level (p-value < 0.05) according to the rank of each subsequent distance class. A value of zero of Moran's I for within model residuals indicates an absence of spatial autocorrelation while a significant deviation from zero suggests that the model is not adequately accounting for spatial autocorrelation.

We evaluated the models using (a) 10-fold cross-validated area under the receiver characteristic curve (AUC) scores (i.e., receiver operating characteristic; Hanley & McNeil, 1982), (b) sensitivity (the proportion of observed presences that are predicted as such), (c) specificity (the proportion of observed absences that are predicted as such) and (iv) TSS (sensitivity + specificity - 1). We first ran the models with a 10-fold cross-validation to get performance metrics and then fitted the model again using the whole set of points to build the final models. We assessed significance for differences in all indices (sensitivity, specificity and TSS) among the summed models for each classification level by (a) randomly resampling 90% of occurrences, (b) developing the SDM based on these resampled datasets and (c) calculating sensitivity/specificity/TSS based on the "raw observation datasets." We repeated this procedure 50 times prior to assess potential sensitivity/specificity/TSS index differentiation between classification levels with a global repeated-measures analysis of variance and subsequent pairwise paired Student's t tests with p-values adjusted for multiple testing using Benjamini and Hochberg's correction. All analyses were conducted in R using the packages dismo version 1.4 (Hijmans, Phillips, Leathwick, & Elith, 2016) and GBM version 2.1. (Ridgeway, 2013).

To compare prediction and projection ability of the three different approaches (species-based model, lineage-based models and ESU-based models), we overlaid the results for the individual lineage-based and ESU-based models and considered the highest value for occurrence probability. These probabilities were again transformed into presence/absence data by maximizing TSS, and performance of the summed models was evaluated using sensitivity, specificity and TSS. We based these assessments on the "raw observation dataset"

to take into account all observations for each species (including observations for which we were not able to determine the lineage/ESU status).

## 2.6 | Future projections

We aimed at examining the impact of potential differences among the three SDM approaches on future projections of climatically suitable areas. We used the respective climate values obtained for the five global circulation models (GCMs; i.e., CCSM4, HadGEM2-AO, MIROC-ESM, MPI-ESM-LR and NorESM1-M) and according to the representative concentration pathways 8.5 scenario (Moss et al., 2008). We used only one climate change scenario as we aimed only at comparing the general consequences when using the three different SDM approaches and not at performing for a comprehensive assessment of potential future changes in the three species. We obtained these data from the WorldClim database (Hijmans et al., 2005) version 1.4 (downloaded 01 October 2017) at a resolution of 10 arc-minutes for the year 2070. The summed lineage-based and summed ESU-based occurrence predictions for 2070 were obtained using the same procedure as used for the current conditions. For each GCM, a particular grid cell is considered as climatically suitable in the future when its occurrence probability in 2070 is higher than the threshold maximizing the TSS obtained in distribution models. We applied an ensemble approach to identify levels of agreement among the projections of the different GCMs (Araújo & New, 2007): We considered a particular grid cell as climatically suitable in the future when at least three different GCMs considered it as such. Finally, we estimated the proportion of (a) the geographic range that could be gained or lost in 2070 and (b) the current predicted geographic range that could remain climatically suitable in 2070 for each classification level.

# 3 | RESULTS

As the results were quite similar at the two different resolutions, we here present only the results at 10 arc-minutes resolution. The results at 30 arc-minutes resolution are shown in Supporting Information (Figures S2 and S4), and only major differences with the other resolution are given in the main text.

## 3.1 | Climatic niche differentiation

The first three axes of the PCA based on the overall climatic conditions in the study area accounted for 95% of the total variability (1st axis: 41%; 2nd axis: 33%; and 3rd axis: 21%). The bioclimatic variables annual mean temperature, maximum temperature of the warmest month and precipitation seasonality were highly correlated with the first axis (Pearson correlation coefficients = 0.70, 0.97 and 0.54, respectively) while temperature seasonality and precipitation of the wettest month highly loaded on the second axis (Pearson correlation coefficients = -0.92 and 0.70, respectively).
Our analyses indicated climatic niche differentiation within each of the three species. All species, lineages and ESUs showed significant OMI values indicating distinct climatic requirements compared to overall conditions within the study area (Supporting Information Table S1). However, the high residual tolerance values for most of the species/lineages/ESUs (Supporting Information Table S1) suggested that in addition to their strong dependence on climatic conditions other factors may also affect their distribution. The OMI analyses indicated differentiation of the mean climatic niche positions within the species (Figure 1, Supporting Information Figure S2). The differences were confirmed by the global MANOVA (B. lapidarius: Pillai's trace = 0.18; B. pascuorum: Pillai's trace = 0.09; B. terrestris: Pillai's trace = 0.33; p-value < 0.01 for all species), and pairwise analyses showed that almost all lineages and ESUs within a species differed significantly

in their mean climatic niche positions (all *p*-value < 0.01) except for two ESUs within *B. terrestris* occurring in the south of Italy and in Sardinia, respectively (*Calabricus* vs. *Sassaricus*, *p*-value = 0.12 at 30 arc-minutes resolution). According to the different niche characteristics, also the relative importance of each chosen bioclimatic variable differed among lineages and ESUs and the species level model (Figure 1 and Supporting Information Figure S2).

## 3.2 | SDM performance comparison

The correlograms of the Moran's I statistics did not indicate significant autocorrelation (Supporting Information Figure S3). Model performance was good to excellent for the species-level models and each of the lineage- and ESU-based models (Table 2 and Supporting Information Figure S2). The comparisons of performances between

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

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

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

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

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

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

**FIGURE 1** Outlying mean index (OMI) analyses of each classification level (species, lineage and evolutionary significant units) for three bumblebee species. These analyses are based on the datasets at the resolution of 10 arc-minutes. Spaces occupied by studied species along the first and second axes of the OMI analysis are given on the left with the projection of mean position (centre of red crosses) and its variance (branches of the red crosses) occupied by each group of each classification level. Canonical weights of the five bioclimatic variables are given on the right (Bio1 = annual mean temperature, Bio4 = temperature seasonality, Bio5 = maximum temperature of the warmest month, Bio13 = precipitation of the wettest month and Bio15 = precipitation seasonality)
alternative models showed quite similar results at each resolution (Tables 2 and Supporting Information Table S2). AUC values were higher for the species-level models than for other models (except for *B. pascuorum*; Tables 2 and Supporting Information Table S2). Overall, TSS values showed a lower prediction ability of summed ESU-based models compared to the species-level models (Tables 2 and Supporting Information Tables S2-S4). TSS values of lineage-based models were as good as species-level models for *B. pascuorum* but were worse for *B. lapidarius* (Tables 2 and Supporting Information Table S2). In general, summing up lineage-based and ESU-based

models led to larger areas considered as climatically suitable compared to species-based models (Figure 2 and Supporting Information Figure S4; see individual SDM for each lineage/ESU in Supporting Information Figure S5). This was indicated by lower model specificity but correspondingly higher sensitivity for lineage-based and ESU-based models (Tables 2, Supporting Information Tables S2 and S4).

Summed lineage-based and summed ESU-based models did not only predict larger areas as climatically suitable compared to species-based models but also the unsuitable areas differed from species-based models (Figure 2). Using species-based models as a

<table><caption>TABLE 2 Evaluation of alternative models for the three bumblebee species</caption><thead><tr><th>Species</th><th>Models</th><th>Lineages/ESUs</th><th>AUC</th><th>Sensitivity</th><th>Specificity</th><th>TSS</th></tr></thead><tbody><tr><td rowspan="12">B. lapidarius</td><td>Species-level model</td><td>-</td><td>0.86</td><td>0.80</td><td>0.79</td><td>0.58</td></tr><tr><td rowspan="4">Lineage-based model</td><td>Main lineage</td><td>0.87</td><td>0.80</td><td>0.80</td><td>0.60</td></tr><tr><td>South-eastern lineage</td><td>0.89</td><td>0.87</td><td>0.87</td><td>0.74</td></tr><tr><td>South Italian lineage</td><td>0.98</td><td>1.00</td><td>0.99</td><td>0.99</td></tr><tr><td>Summed lineage-based model</td><td>-</td><td>0.79</td><td>0.84</td><td>0.65</td><td>0.49</td></tr><tr><td rowspan="6">ESU-based model</td><td>Lapidarius</td><td>0.87</td><td>0.81</td><td>0.81</td><td>0.62</td></tr><tr><td>Atlanticus</td><td>0.70</td><td>1.00</td><td>0.99</td><td>0.99</td></tr><tr><td>Decipiens</td><td>0.97</td><td>0.97</td><td>0.97</td><td>0.93</td></tr><tr><td>SE Lapidarius</td><td>0.88</td><td>0.87</td><td>0.87</td><td>0.74</td></tr><tr><td>SI Decipiens</td><td>0.97</td><td>1.00</td><td>0.99</td><td>0.99</td></tr><tr><td>Summed ESU-based model</td><td>-</td><td>0.80</td><td>0.87</td><td>0.62</td><td>0.49</td></tr><tr><td rowspan="12">B. pascuorum</td><td>Species-level model</td><td>-</td><td>0.81</td><td>0.74</td><td>0.74</td><td>0.48</td></tr><tr><td rowspan="2">Lineage-based model</td><td>Main lineage</td><td>0.83</td><td>0.75</td><td>0.75</td><td>0.51</td></tr><tr><td>South-western lineage</td><td>0.98</td><td>0.96</td><td>0.96</td><td>0.92</td></tr><tr><td>Summed lineage-based model</td><td>-</td><td>0.80</td><td>0.79</td><td>0.69</td><td>0.48</td></tr><tr><td rowspan="5">ESU-based model</td><td>Main unit</td><td>0.83</td><td>0.76</td><td>0.75</td><td>0.51</td></tr><tr><td>Siciliensis</td><td>0.70</td><td>1.00</td><td>0.99</td><td>0.99</td></tr><tr><td>Rehbinderi</td><td>0.88</td><td>0.93</td><td>0.93</td><td>0.85</td></tr><tr><td>SW unit</td><td>0.98</td><td>0.96</td><td>0.96</td><td>0.92</td></tr><tr><td>Summed ESU-based model</td><td>-</td><td>0.81</td><td>0.82</td><td>0.61</td><td>0.43</td></tr><tr><td rowspan="12">B. terrestris</td><td>Species-level model</td><td>-</td><td>0.86</td><td>0.77</td><td>0.78</td><td>0.55</td></tr><tr><td rowspan="11">ESU-based model</td><td>Africanus</td><td>0.99</td><td>0.99</td><td>0.98</td><td>0.97</td></tr><tr><td>Audax</td><td>0.99</td><td>0.96</td><td>0.96</td><td>0.92</td></tr><tr><td>Calabricus</td><td>0.98</td><td>1.00</td><td>0.99</td><td>0.99</td></tr><tr><td>Canariensis</td><td>0.99</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><td>Dalmatinus</td><td>0.91</td><td>0.87</td><td>0.87</td><td>0.74</td></tr><tr><td>Lusitanicus</td><td>0.96</td><td>0.93</td><td>0.93</td><td>0.86</td></tr><tr><td>Sassaricus</td><td>0.99</td><td>1.00</td><td>0.99</td><td>0.99</td></tr><tr><td>Terrestris</td><td>0.92</td><td>0.86</td><td>0.86</td><td>0.71</td></tr><tr><td>Summed ESU-based model</td><td>-</td><td>0.78</td><td>0.94</td><td>0.59</td><td>0.53</td></tr></tbody></table>

Note. NA = not analysed.

TSS is true skill statistics. Results are displayed at the resolutions of 10 arc-minutes. AUC is 10-fold cross-validated area under the receiver characteristic curve.
baseline (baseline predicted distribution [BPD] = 100% for each species), the predicted suitable range by summed lineage-based models was 157% of the BPD for *B. lapidarius* but included only 99% of the BPD and 118% for *B. pascuorum* but included only 96% respectively of the BPD. Summed ESU-based models for *B. lapidarius* predicted 165% of the BPD including 99% of the BPD, for *B. pascuorum* 157% of the BPD including 97% of the BPD and for *B. terrestris* 124% of the BPD but including only 90% of the BPD (Figure 2). For all species, lineage- and ESU-based SDMs predicted also a larger geographic extent for Eastern and Southern Europe (Figure 2 and Supporting Information Figure S4).

## 3.3 | Future projections

As performances of each SDM were similar between the two resolutions, we investigated future projections at 10 arc-minutes resolution. Climate change under the representative concentration pathways 8.5 scenario was projected to have considerable effects on future distributions of all three species (Figure 3 and Supporting Information Figure S6). However, the severity of this impact depended on the classification level used to assess species responses (Table 3; Figure 3). In particular, species-level models were more

pessimistic than lineage-based or ESU-based models for all species in terms of relative loss of suitable area until 2070 (Table 3).

# 4 | DISCUSSION

Our analyses were based on occurrence data from several previous field surveys. We argue that these datasets can be considered reliable for SDM analyses (see Methods). However, we cannot completely rule out that field surveys failed to detect or disregarded (i.e., some surveys may have only targeted a particular species and therefore did not record others) our species of concern. Therefore, climatic niche and modelled current and future distribution of the three bumblebee species should be considered with caution.

## 4.1 | Climatic niche differentiation within species

Our analyses show that the climatic niches of most of the lineages and ESUs are distinct (Figure 1; Supporting Information Table S1), as also observed in mammals (D'Amen et al., 2013) or other insects (Homburg, Brandt, Drees, & Assmann, 2014). One could argue that this is an expected result since the environment is spatially

<figure><img src="image_8.png" /><figcaption>**FIGURE 2** Results of species distribution models for three bumblebee species. Potential climatically suitable area according to species-based models, summed lineage-based models and summed ESU-based models for three bumblebee species. Green areas indicate model predictions. Analyses are based on datasets at 10 arc-minutes resolution. The thresholds used to transform the occurrence probabilities into presence/absence map are those maximizing the true skill statistic: 0.56/0.7/0.72 for *Bombus lapidarius*, 0.5/0.61/0.65 for *Bombus pascuorum* and 0.6/0.6/0.67 for *Bombus terrestris* (for species-level model/for summed lineage-based model/for summed ESU-based model)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**FIGURE 3** Predictions of suitable climatic areas in 2070 for three bumblebee species according to alternative species distribution models based on the consensus of five general circulation models. Changes in suitable climatic conditions projected for 2070 according to the representative concentration pathways 8.5 scenario for species-based models, summed lineage-based models and summed ESU-based models of three bumblebee species: orange, conserved area; red, lost area; green, gained area. Analyses are based on datasets at 10 arc-seconds resolution. The thresholds used to transform the occurrence probabilities into presence/absence map are those maximizing the true skill statistic: 0.56/0.7/0.72 for *Bombus lapidarius*, 0.5/0.61/0.65 for *Bombus pascuorum* and 0.6/0.6/0.67 for *Bombus terrestris* (for species-level model/for summed lineage-based model/for summed ESU-based model)</figcaption></figure>

autocorrelated and the intraspecific units have allopatric distributions. Therefore, the results could be explained by geographic range differences shaped by biogeographic history or limited dispersal ability rather than by actual niche divergence. However, a large body of work documents differences in climatic requirements between allopatric populations (e.g., Addo-Bediako, Chown, & Gaston, 2000; Sunday, Bates, & Dulvy, 2011), including bumblebees (mainly *B. terrestris*) through bioassays (de Jonghe, 1986; Oyen, Giri, & Dillon, 2016; Velthuis & van Doorn, 2006). Therefore, we assume that there are indeed niche differentiations for the analysed bumblebee species, but further assessments, for example, via bioassays (e.g., Martinet, Lecocq et al., 2015), are needed.

## 4.2 | Does integrating ESU and lineage information improve species distribution models?

The existence of spatial autocorrelation may affect the performance metrics (e.g., Crase et al., 2014) which we used to compare species-level, lineage-based and ESU-based models. Although the Moran's I statistics were high for some distance classes, they were

not significant. Therefore, we consider that our results were poorly impacted by potential issues triggered by the presence of spatial autocorrelation in model residuals.

In our case, integrating lineage or ESU information did not improve overall measures of accuracy like TSS or AUC, but model sensitivity and specificity changed (ESU and lineage-based SDMs vs. species-based SDM). While prediction of presence data was improved (i.e., increased sensitivity), decreased specificity led to constant over-prediction for all three species ranges. However, these results must be qualified. On one hand, the specificity of ESU and lineage-based SDMs could be underestimated in our analyses because of the lack of observations in some areas: ESU and lineage-based SDMs predicted larger extents of suitable climatic conditions in Eastern Europe where sampling intensity was generally lower (Figure 2) although the species most likely actually occur there (see Rasmont et al., 2015). On the other hand, the higher sensitivity of lineage/ESU-based SDMs could be a result of the over-prediction of these models compared to the species-based SDMs. Indeed, the exclusion of records from Russia, Belorussia, Kazakhstan and Ukraine for SDM developments can bias the models to under-predict range
<table><caption>**TABLE 3** Summary of conserved, lost and gained areas in 2070 for different classification levels of three bumblebee species according to the representative concentration pathways 8.5 scenario and the consensus of five general circulation models under the assumption of unlimited dispersal. Conserved means area projected as remaining suitable. Lost means area projected as turning unsuitable. Gained means novel area projected to become suitable. Total gives the extent of suitable area in 2070. The percentage reflects the future surface of suitable climatic areas relative to the current prediction</caption><thead><tr><th>Species</th><th>Models</th><th>Conserved (%)</th><th>Lost (%)</th><th>Gained (%)</th><th>Total (%)</th><th>Global trend</th></tr></thead><tbody><tr><td rowspan="3">B. lapidarius</td><td>Species-level model</td><td>29</td><td>71</td><td>9</td><td>38</td><td>Displacement and decrease in global suitable area (~62%)</td></tr><tr><td>Summed lineage-based model</td><td>53</td><td>47</td><td>17</td><td>70</td><td>Displacement and decrease in global suitable area (~30%)</td></tr><tr><td>Summed ESU-based model</td><td>51</td><td>49</td><td>16</td><td>66</td><td>Displacement and decrease in global suitable area (~34%)</td></tr><tr><td rowspan="3">B. pascuorum</td><td>Species-level model</td><td>65</td><td>35</td><td>16</td><td>80</td><td>Displacement and decrease in global suitable area (~20%)</td></tr><tr><td>Summed lineage-based model</td><td>54</td><td>46</td><td>14</td><td>68</td><td>Displacement and decrease in global suitable area (~31%)</td></tr><tr><td>Summed ESU-based model</td><td>75</td><td>25</td><td>36</td><td>111</td><td>Displacement and increase in global suitable area (+11%)</td></tr><tr><td rowspan="2">B. terrestris</td><td>Species-level model</td><td>70</td><td>29</td><td>8</td><td>78</td><td>Displacement and decrease in global suitable area (~22%)</td></tr><tr><td>Summed ESU-based model</td><td>77</td><td>23</td><td>13</td><td>89</td><td>Displacement and decrease in global suitable area (~11%)</td></tr></tbody></table>

suitability in Eastern Europe. Therefore, the prediction by lineage/ESU-based SDMs that these areas are climatically suitable for the species (thus increasing the sensitivity) could be an artefact of the over-prediction of these models rather than a performance increase. Overall, our current results do not highlight that integrating ESU or lineage information increases the SDM performance under current climatic conditions. Therefore, none of the alternative SDM approaches appears to be superior under current conditions.

climatic changes can be the source of future species expansion, similarly to the population movements fostered by past climatic oscillations (Stewart, Lister, Barnes, & Dalén, 2010). Therefore, integrating intraspecific differentiation into models (i.e., summed lineage/ESU-based SDMs) may highlight potential resilience of particularly differentiated units acting as potential buffer against adverse effects of climate change (Benito Garzón et al., 2011; Oney, Reineking, O'Neill, & Kreyling, 2013).

## 4.3 | Integrating intraspecific differentiation for climatic risk assessment

Our projections underline different climatically suitable areas in 2070 for all species but with pronounced differences depending on the classification level used for SDMs (Figure 3). Indeed, the extent of change and the proportion of conserved, lost and gained areas depended on the kind of model used (Table 3): (a) The losses are critical for *B. lapidarius* and *B. terrestris* with the species-level models but are much less severe for lineage/ESU-based SDMs; (b) the losses are smaller at species-level model for *B. pascuorum* but the species range predicted by this model is much smaller than those predicted by other models. Overall, for these three species with comparable distribution and overall niche requirements (a) predicted gains are much larger with ESU-based SDMs, especially to Northern Europe which can be considered as reliable given current observations of range expansions in Scandinavia (Martinet, Rasmont et al., 2015b; but see Kerr et al., 2015) and (b) predicted losses are low with ESU-based SDMs when the species includes more warm/dry-adapted units (e.g., *B. lapidarius* and *B. terrestris* vs. *B. pascuorum*). This could be biologically relevant since even small units able to survive under

## 4.4 | Limitations and relevance of the approach

The main problem associated with integrating intraspecific in SDMs is the required split of species into smaller units that are geographically separated from each other. This could result in two alternative issues. First, the resulting projections are likely to suffer even more from the mismatch between fundamental and realized niches than whole-species projections. Indeed, SDMs are based on observations to estimate the realized niche although a subspecific taxon/unit range can be restricted by limited dispersal or geographic barriers rather than by climatic requirements. As the development of individual lineage/ESU models requires to remove all other lineages/ESUs and to consider areas inhabited by other units as unoccupied by the target unit, the climatic niche could be underestimated. Subsequently, this may *per se* cause errors when projecting SDMs into new areas, since suitable areas may be excluded although being physiologically suitable for the target taxon (Pulliam, 2000). This issue could happen for each independent SDM of each intraspecific unit leading to dramatic consequences on the summed SDM based on all of these independent models. Second, the combination of several SDMs of small intraspecific units (i.e., which is based on
less observations) can lead to an overestimation of species range because smaller observation datasets can decrease the SDM precision resulting in broader climatic parameters inference (Stockwell & Peterson, 2002). In our study, the quite similar performances and the larger species range predictions of ESU and lineage-based SDMs do not confirm the first concern but could reflect the second potential issue.

Another potential issue is the relevance of the a priori classification of differentiated populations. Indeed, intraspecific units defined on the basis of genetic markers or ESU delimitation may not display physiologically different responses to climate change. For instance, a previous study showed that tolerance to heat is largely conserved among closely related species (Araújo et al., 2013). This would suggest that heat tolerance may also be conserved within species. However, several empirical studies (including studies on bumblebee species) underlined differentiations in genes coding for heat tolerance (Du, Li, Zhang, Meng, & Zhang, 2014) or in heat resistance (Homburg et al., 2014; Martinet, Lecocq et al., 2015; Oyen et al., 2016; Slatyer, Nash, & Hoffmann, 2016; Sorte, Jones, & Miller, 2011) among conspecific populations or ESUs. Therefore, such assessments require experiments to determine the tolerances to all climatic parameters of differentiated populations. However, such experimental assessments are still not available for most of the species.

We show for the example of the three selected bumblebee species that integrating lineage or ESU information does not necessarily improve overall accuracy of conventional species-level SDMs but it can lead to considerably different conclusions about future climatic risks. From a conservation point of view, although all models indicate displacement and decrease in global suitable area in the future, the divergent predictions can deeply impact the planning of conservation action, especially in the current context of spatially prioritizing conservation actions (Margules & Pressey, 2000). As SDM-based climatic risk predictions are increasingly used to produce guidelines to (non-)governmental agencies and to assess conservation plans (Carroll, 2010), divergences and limitations of each modelling approach should be taken into account for developing efficient biodiversity management strategies. Lineage and ESU-based SDMs offer the advantage to drag attention to species in which allopatric populations could display physiologically different responses to climate change when they lead to different results than species-based models. Therefore, they pave the way to further assessments of ecoclimatic tolerance of particular species where intraspecific differentiation could lead to an unexpected response to climate change.

## ACKNOWLEDGEMENTS

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement no 244090, STEP Project (Status and Trends of European Pollinators, www.step-project.net) (Potts et al., 2011).

## CONFLICT OF INTERESTS

We have no competing interests.

## DATA ACCESSIBILITY

The occurrence data came from the database "Base de données fauniques Gembloux-Mons" and were published in Rasmont et al. (2015).

## ORCID

Thomas Lecocq https://orcid.org/0000-0002-4947-0332

Oliver Schweiger https://orcid.org/0000-0001-8779-2335

## REFERENCES

* Addo-Bediako, A., Chown, S. L., & Gaston, K. J. (2000). Thermal tolerance, climatic variability and latitude. *Proceedings of the Royal Society B: Biological Sciences*, 267, 739-745. https://doi.org/10.1098/rspb.2000.1065

* Alexander, J. M., & Edwards, P. J. (2010). Limits to the niche and range margins of alien species. *Oikos*, 119, 1377-1386. https://doi.org/10.1111/j.1600-0706.2009.17977.x

* Allouche, O., Tsoar, A., & Kadmon, R. (2006). Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). *Journal of Applied Ecology*, 43, 1223-1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x

* Araújo, M. B., Ferri-Yáñez, F., Bozinovic, F., Marquet, P. A., Valladares, F., & Chown, S. L. (2013). Heat freezes niche evolution. *Ecology Letters*, 16, 1206-1219. https://doi.org/10.1111/ele.12155

* Araújo, M. B., & New, M. (2007). Ensemble forecasting of species distributions. *Trends in Ecology and Evolution*, 22, 42-47. https://doi.org/10.1016/j.tree.2006.09.010

* Avise, J. C. (2000). *Phylogeography: The history and formation of species*. Cambridge, MA: Harvard University Press.

* Benito Garzón, M., Alia, R., Robson, T. M., & Zavala, M. A. (2011). Intraspecific variability and plasticity influence potential tree species distributions under climate change. *Global Ecology and Biogeography*, 20, 766-778. https://doi.org/10.1111/j.1466-8238.2010.00646.x

* Bjornstad, O. N., & Cai, J. (2018). *ncf: Spatial covariance functions*. R package version 1.2-6.

* Carroll, C. (2010). Role of climatic niche models in focal-species-based conservation planning: Assessing potential effects of climate change on Northern Spotted Owl in the Pacific Northwest, USA. *Biological Conservation*, 143, 1432-1437. https://doi.org/10.1016/j.biocon.2010.03.018

* Conner, J. K., & Hartl, D. L. (2004). A primer of ecological genetics. Sunderland, MA: Sinauer Associates.

* Costa, W. J. E. M., & Amorim, P. F. (2014). Integrative taxonomy and conservation of seasonal killifishes, *Xenurolebias* (Teleostei: Rivulidae), and the Brazilian Atlantic Forest. *Systematics and Biodiversity*, 12, 350-365.

* Crase, B., Liedloff, A., Vesk, P. A., Fukuda, Y., & Wintle, B. A. (2014). Incorporating spatial autocorrelation into species distribution models alters forecasts of climate-mediated range shifts. *Global Change Biology*, 20, 2566-2579. https://doi.org/10.1111/gcb.12598

* Cruaud, A., Gautier, M., Galan, M., Foucaud, J., Sauné, L., Genson, G., ... Rasplus, J.-Y. (2014). Empirical assessment of RAD sequencing for interspecific phylogeny. *Molecular Biology and Evolution*, 31, 1272-1274. https://doi.org/10.1093/molbev/msu063
D'Amen, M., Zimmermann, N. E., & Pearman, P. B. (2013). Conservation of phylogeographic lineages under climate change. *Global Ecology and Biogeography*, **22**, 93–104. https://doi.org/10.1111/j.1466-8238.2012.00774.x

de Jonghe, R. (1986). Crossing experiments with *Bombus terrestris terrestris* (LINNAEUS, 1758) and *Bombus terrestris xanthopus* KRICHBAUMER, 1870 and some notes on diapause and nosemose (Hymenoptera: Apoidea). *Phegea*, **14**, 19–23.

Dolédec, S., Chessel, D., & Gimaret-Carpentier, C. (2000). Niche separation in community analysis: A new method. *Ecology*, **81**, 2914–2927. https://doi.org/10.2307/177351

Du, X., Li, L., Zhang, S., Meng, F., & Zhang, G. (2014). SNP identification by transcriptome sequencing and candidate gene-based association analysis for heat tolerance in the bay scallop *Argopecten irradians*. *PLoS One*, **9**, e104960.

Elith, J., Leathwick, J. R., & Hastie, T. (2008). A working guide to boosted regression trees. *The Journal of Animal Ecology*, **77**, 802–813. https://doi.org/10.1111/j.1365-2656.2008.01390.x

Ferguson, J. W. H. (2002). On the use of genetic divergence for identifying species. *Biological Journal of the Linnean Society*, **75**, 509–516. https://doi.org/10.1046/j.1095-8312.2002.00042.x

Frankham, R., Ballou, J. D., & Briscoe, D. A. (2010). *Introduction to conservation genetics*. Cambridge, UK: Cambridge University Press.

Guisan, A., & Zimmermann, N. E. (2000). Predictive habitat distribution models in ecology. *Ecological Modelling*, **135**, 147–186. https://doi.org/10.1016/S0304-3800(00)00354-9

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

Hendrixson, B. E., Guice, A. V., & Bond, J. E. (2015). Integrative species delimitation and conservation of tarantulas (Araneae, Mygalomorphae, Theraphosidae) from a North American biodiversity hotspot. *Insect Conservation and Diversity*, **8**, 120–131. https://doi.org/10.1111/icad.12089

Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the Quaternary. *Philosophical Transactions of the Royal Society B-Biological Sciences*, **359**, 183–195. https://doi.org/10.1098/rstb.2003.1388

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., & Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. *International Journal of Climatology*, **25**, 1965–1978. https://doi.org/10.1002/joc.1276

Hijmans, R. J., Phillips, S., Leathwick, J., & Elith, J. (2016). *Package "dismo"*.

Homburg, K., Brandt, P., Drees, C., & Assmann, T. (2014). Evolutionarily significant units in a flightless ground beetle show different climate niches and high extinction risk due to climate change. *Journal of Insect Conservation*, **18**, 781–790. https://doi.org/10.1007/s10841-014-9685-x

Kerr, J. T., Pindar, A., Galpern, P., Packer, L., Potts, S. G., Roberts, S. M., ... Pantoja, A. (2015). Climate change impacts on bumblebees converge across continents. *Science*, **349**, 177–180. https://doi.org/10.1126/science.aaa7031

Kissling, W. D., & Carl, G. (2007). Spatial autocorrelation and the selection of simultaneous autoregressive models. *Global Ecology and Biogeography*, **17**, 59–71. https://doi.org/10.1111/j.1466-8238.2007.00334.x

Lecocq, T., Brasero, N., Martinet, B., Valterová, I., & Rasmont, P. (2015). Highly polytypic taxon complex: Interspecific and intraspecific integrative taxonomic assessment of the widespread pollinator *Bombus pascuorum* Scopoli 1763 (Hymenoptera: Apidae). *Systematic Entomology*, **40**, 881–888.

Lecocq, T., Brasero, N., Meulemeester, T., De, M. D., Dellicour, S., Lhomme, P., ... Rasmont, P. (2015b). An integrative taxonomic approach to assess the status of Corsican bumblebees: Implications for conservation. *Animal Conservation*, **18**, 236–248.

Lecocq, T., Coppée, A., Michez, D., Brasero, N., Rasplus, J.-Y., Valterová, I., & Rasmont, P. (2016). The alien's identity: Consequences of

taxonomic status for the international bumblebee trade regulations. *Biological Conservation*, **195**, 169–176.

Lecocq, T., Dellicour, S., Michez, D., Dehon, M., Dewulf, A., De Meulemeester, T., ... Rasmont, P. (2015). Methods for species delimitation in bumblebees (Hymenoptera, Apidae, *Bombus*): Towards an integrative approach. *Zoologica Scripta*, **44**, 281–297.

Lecocq, T., Rasmont, P., Harpke, A., & Schweiger, O. (2016). Improving international trade regulation by considering intraspecific variation for invasion risk assessment of commercially traded species: The *Bombus terrestris* case. *Conservation Letters*, **9**, 281–289.

Lynch, M. (1996). A quantitative-genetic perspective on conservation issues. In J. C. Avise, & J. L. Hamrick (Eds.), *Conservation genetics. Case histories from nature* (pp. 471–501). New York, NY: Chapman and Hall Ltd.

Maggini, R., Lehmann, A., Zimmermann, N. E., & Guisan, A. (2006). Improving generalized regression analysis for the spatial prediction of forest communities. *Journal of Biogeography*, **33**, 1729–1749. https://doi.org/10.1111/j.1365-2699.2006.01465.x

Mainali, K. P., Warren, D. L., Dhileepan, K., McConnachie, A., Strathie, L., Hassan, G., ... Parmesan, C. (2015). Projecting future expansion of invasive species: Comparing and improving methodologies for species distribution modeling. *Global Change Biology*, **21**, 4464–4480. https://doi.org/10.1111/gcb.13038

Marcer, A., Méndez-Vigo, B., Alonso-Blanco, C., & Picó, F. X. (2016). Tackling intraspecific genetic structure in distribution models better reflects species geographical range. *Ecology and Evolution*, **6**, 2084–2097. https://doi.org/10.1002/ece3.2010

Margules, C. R., & Pressey, R. L. (2000). Systematic conservation planning. *Nature*, **405**, 243–253. https://doi.org/10.1038/35012251

Martinet, B., Lecocq, T., Smet, J., & Rasmont, P. (2015). A protocol to assess insect resistance to heat waves, applied to bumblebees (*Bombus Latreille*, 1802). *PLOS One*, **10**, e0118591.

Martinet, B., Rasmont, P., Cederberg, B., Evrard, D., Ødegaard, F., Paukkunen, J., & Lecocq, T. (2015b). Forward to the north: Two Euro-Mediterranean bumblebee species now cross the Arctic Circle. *Annales De La Société Entomologique De France (N.S.)*, **51**, 303–309.

Moran, P. (1950). Notes on continuous stochastic phenomena. *Biometrika*, **37**, 17–23. https://doi.org/10.1093/biomet/37.1.2.17

Moss, R., Babiker, M., Brinkman, S., Calvo, E., Carter, T., Edmonds, J., ... Zurek, M. (2008). *Towards new scenarios for analysis of emissions, climate change, impacts, and response strategies*.

Oney, B., Reineking, B., O'Neill, G., & Kreyling, J. (2013). Intraspecific variation buffers projected climate change impacts on *Pinus contorta*. *Ecology and Evolution*, **3**, 437–449.

Oyen, K. J., Giri, S., & Dillon, M. E. (2016). Altitudinal variation in bumble bee (*Bombus*) critical thermal limits. *Journal of Thermal Biology*, **59**, 52–57. https://doi.org/10.1016/j.jtherbio.2016.04.015

Pearman, P. B., D'Amen, M., Graham, C. H., Thuiller, W., & Zimmermann, N. E. (2010). Within-taxon niche structure: Niche conservatism, divergence and predicted effects of climate change. *Ecography*, **33**, 990–1003. https://doi.org/10.1111/j.1600-0587.2010.06443.x

Potts, S. G., Biesmeijer, J. C., Bommarco, R., Felicioli, A., Fischer, M., Jokinen, P., ... Schweiger, O. (2011). Developing European conservation and mitigation tools for pollination services: Approaches of the STEP (Status and Trends of European Pollinators) project. *Journal of Apicultural Research*, **50**, 152–164. https://doi.org/10.3896/IBRA.1.50.2.07

Pulliam, H. R. (2000). On the relationship between niche and distribution. *Ecology Letters*, **3**, 349–361. https://doi.org/10.1046/j.1461-0248.2000.00143.x

R Core Team (2017). *R: A language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing.

Rasmont, P., Franzen, M., Lecocq, T., Harpke, A., Roberts, S., Biesmeijer, K., ... Schweiger, O. (2015). Climatic risk and distribution atlas of
European bumblebees. *BioRisk*, 10, 1-236. https://doi.org/10.3897/biorisk.10.4749

Ridgeway, G. (2013). *gbm: Generalized boosted regression models*. *R Package Version*, 2, 1.

Salvato, P., Battisti, A., Concato, S., Masutti, L., Patarnello, T., & Zane, L. (2002). Genetic differentiation in the winter pine processionary moth (*Thaumetopoea pityocampa* - wilkinsoni complex), inferred by AFLP and mitochondrial DNA markers. *Molecular Ecology*, *11*, 2435-2444. https://doi.org/10.1046/j.1365-294X.2002.01631.x

Schlick-Steiner, B. C., Steiner, F. M., Seifert, B., Stauffer, C., Christian, E., & Crozier, R. H. (2010). Integrative taxonomy: A multisource approach to exploring biodiversity. *Annual Review of Entomology*, *55*, 421-438. https://doi.org/10.1146/annurev-ento-112408-085432

Slatyer, R. A., Nash, M. A., & Hoffmann, A. A. (2016). Scale-dependent thermal tolerance variation in Australian mountain grasshoppers. *Ecography*, *39*, 572-582. https://doi.org/10.1111/ecog.01616

Sorte, C. J. B., Jones, S. J., & Miller, L. P. (2011). Geographic variation in temperature tolerance as an indicator of potential population responses to climate change. *Journal of Experimental Marine Biology and Ecology*, *400*, 209-217. https://doi.org/10.1016/j.jembe.2011.02.009

Stewart, J. R., Lister, A. M., Barnes, I., & Dalén, L. (2010). Refugia revisited: Individualistic responses of species in space and time. *Proceedings of the Royal Society B: Biological Sciences*, *277*, 661-671. https://doi.org/10.1098/rspb.2009.1272

Stockwell, D. R., & Peterson, A. T. (2002). Effects of sample size on accuracy of species distribution models. *Ecological Modelling*, *148*, 1-13. https://doi.org/10.1016/S0304-3800(01)00388-X

Sunday, J. M., Bates, A. E., & Dulvy, N. K. (2011). Global analysis of thermal tolerance and latitude in ectotherms. *Proceedings of the Royal Society B: Biological Sciences*, *278*, 1823-1830. https://doi.org/10.1098/rspb.2010.1295

Thuiller, W., Lavorel, S., Midgley, G., Lavergne, S., & Rebelo, T. (2004). Relating plant traits and species distributions along bioclimatic gradients for 88 Leucadendron taxa. *Ecology*, *85*, 1688-1699. https://doi.org/10.1890/03-0148

Valladares, F., Matesanz, S., Guilhaumon, F., Araújo, M. B., Balaguer, L., Benito-Garzón, M., ... Zavala, M. A. (2014). The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. *Ecology Letters*, *17*, 1351-1364. https://doi.org/10.1111/ele.12348

Velthuis, H. H. W., & van Doorn, A. (2006). A century of advances in bumblebee domestication and the economic and environmental aspects of its commercialization for pollination. *Apidologie*, *37*, 421-451. https://doi.org/10.1051/apido:2006019

Zurell, D., Thuiller, W., Pagel, J., Cabral, J. S., Münkemüller, T., Gravel, D., ... Zimmermann, N. E. (2016). Benchmarking novel approaches for modelling species range dynamics. *Global Change Biology*, *22*, 2651-2664.

## BIOSKETCH

Thomas Lecocq is Associate Professor (senior lecturer) in animal biology and ecology. His current research focuses upon the integration of intraspecific diversity in conservation biology and domestication programs.

Author contributions: T.L., P.R. and O.S. conceived the original idea. T.L., A.H. and O.S. performed statistical analyses. T.L. led the writing. All authors discussed the results and contributed to the final manuscript, revising it critically for content.

## SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section at the end of the article.

**How to cite this article:** Lecocq T, Harpke A, Rasmont P, Schweiger O. Integrating intraspecific differentiation in species distribution models: Consequences on projections of current and future climatically suitable areas of species. *Divers Distrib.* 2019;00:1-13. https://doi.org/10.1111/ddi.12916