# Ecological niche modelling and coalescent simulations to explore the recent geographical range history of five widespread bumblebee species in Europe

Simon Dellicour, Chedly Kastally, Sara Varela, Denis Michez, Pierre Rasmont, Patrick Mardulyn, Thomas Lecocq

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

## To cite this version:

Simon Dellicour, Chedly Kastally, Sara Varela, Denis Michez, Pierre Rasmont, et al.. Ecological niche modelling and coalescent simulations to explore the recent geographical range history of five widespread bumblebee species in Europe. Journal of Biogeography, 2017, 44 (1), pp.39-50. 10.1111/jbi.12748. hal-01575588

HAL Id: hal-01575588

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

Submitted on 30 Aug 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ecological niche modelling and coalescent simulations to explore the recent geographical range history of five widespread bumblebee species in Europe

Simon Dellicour<sup>a</sup>, Chedly Kastally<sup>b</sup>, Sara Varela<sup>c,d</sup>, Denis Michez<sup>e</sup>, Pierre Rasmont<sup>e</sup>, Patrick Mardulyn<sup>b</sup>, Thomas Lecocq<sup>e,f</sup>.

a : Department of Microbiology and Immunology, Rega Institute for Medical Research, Clinical and Epidemiological Virology, KU Leuven—University of Leuven, Minderbroedersstaat 10, 3000 Leuven, Belgium

b : Evolutionary Biology and Ecology, Université Libre de Bruxelles, avenue Franklin Delano Roosevelt 50, 1050 Bruxelles, Belgium

c : Departamento de Ciencias de la Vida, Edificio de Ciencias, Campus Externo, Universidad de Alcalá, 28805 Alcalá de Henares Madrid, Spain

d : Museum für Naturkunde, Leibniz Institute for Evolution and Biodiversity Science, Invalidenstr. 43, 10115 Berlin, Germany

e : Laboratoire de Zoologie, Research institute of Biosciences, University of Mons, Place du Parc 23, 7000 Mons, Belgium

f : Research Unit Animal and Functionalities of Animal Products (URAFPA), University of Lorraine – INRA, 2 Avenue de la Forêt de Haye, BP 172, 54505 Vandoeuvre-Lès-Nancy, France

Published in *Journal of Biogeography*. Available online: https://doi.org/10.1111/jbi.12748

## Keywords

Bombus, coalescent simulations, comparative phylogeography, DNA sequences, ecological niche modelling, glacial refugia, last glaciation.

## Abstract

Aim Studying the changes in species ranges during the last glaciation event is an important step towards the understanding of the observed patterns of intra-specific genetic variability. We focused on bumblebees, an interesting biological model to address these questions because cold-adapted species are likely to have experienced different geographical range histories during the last glacial period compared to more commonly studied, strictly temperate, species. We investigated and compared historical hypotheses regarding the geographical range of five common and co-distributed West Palaearctic bumblebee species.

Location Europe, West Palaearctic.

Methods For each species, we inferred present and past (Last Glacial Maximum) distributions from species occurrence records, and present and past climatic data, using the ecological niche modelling (ENM) approach implemented in MAXENT. Based on genetic data previously obtained from the sequencing of three gene fragments (mitochondrial locus COI and two nuclear loci EF-1a and PEPCK), we then compared global and local patterns of genetic variation using several summary statistics as well as a visual mapping of genetic variation. Finally, we used a spatially explicit model of DNA sequence coalescence to test and compare four evolutionary scenarios derived from ENM results and patterns of genetic diversity.

Results Ecological niche modelling results based on climatic data clearly suggested a range continuum in Europe during the last glaciation. Yet, the related evolutionary scenario involving such continuum was less supported than alternative scenarios involving a more fragmented distribution. Indeed, for the three
out of five species for which genetic data allowed discriminating among tested scenarios, the scenario
that included a fragmented range during the last glaciation was identified as the most likely.

Main conclusions Although ENM suggested that bumblebees would have maintained a range continuum across Europe during the last glaciation, coalescent simulations tended to refute the persistence of a large range continuum for these species during this period. This suggests that even for cold-adapted species, the cooling periods have significantly shrunk and fragmented their respective ranges.

# 1. Introduction

The demographic history of many species is strongly impacted by Pleistocene climatic oscillations (Andersen & Borns, 1994), especially in Europe (Avise, 2000; Hewitt, 2004). It is generally acknowledged that during glacial episodes, the ranges of temperate species are reduced (e.g. Varela et al., 2010), as they cannot survive in extremely cold regions, in particular, those covered by an ice sheet. As a result, during glaciations temperate species are confined to small regions with more favourable climates, which are referred to as glacial refugia (Zagwijn, 1982; Stewart et al., 2010). At the end of each glacial cycle, populations expand from their glacial refugia to recolonize at least part of their previous range (Taberlet et al., 1998; Varela et al., 2010). While this general view allowed data interpretation and historical inference in many previous phylogeographical studies, a more complex and realistic representation of species range changes during the Pleistocene has emerged more recently (e.g. Stewart et al., 2010). This representation involves the potential presence of additional cryptic glacial refugia in central and northern Europe, especially for cold-tolerant species, and now requires to be tested using various biological models. Among potential organisms of interest, bumblebees (Hymenoptera, Apidae, genus Bombus) are an excellent system because (1) they live in some of the coldest areas inhabited by insects and (2) have been able to respond relatively quickly to climate changes, recolonizing the depopulated areas after each glacial cycle during the last 3 Myr (Hines, 2008). Thus, Bombus species represent a relevant biological model to investigate the geographical range history of cold-adapted species during the last glaciation. Although the early biogeographical history (from Eocene-Oligocene to Miocene) is well studied (Williams, 1985; Hines, 2008; Wappler et al., 2012; Dehon et al., 2014), their recent history (from Pliocene to Holocene) has received comparatively little attention. Some studies have investigated patterns of genetic variation of West Palaearctic bumblebees at the intra-species level (e.g. Estoup et al., 1996; Widmer & Schmid-Hempel, 1999; Lecocq et al., 2013a,b; Dellicour et al., 2015), but so far, a comparative phylogeographical analysis of several widespread and common West Palaearctic bumblebees is still lacking.

Hypotheses testing based on coalescent simulations of genetic data (Knowles, 2009; Hickerson et al., 2010) has become an important tool in phylogeographical studies to investigate past demographic history (e.g. Mardulyn et al., 2009; Ray et al., 2010a; Veeramah et al., 2012; Metcalf et al., 2014). In practice, as an infinite number of scenarios can be formulated, a potential limitation of this approach is the

necessity to select, a priori, a restricted (and manageable) number of clearly defined alternative evolutionary scenarios to be tested. Prior knowledge about the history of a species can be used to build these scenarios (e.g. Knowles, 2001; Fagundes et al., 2007; Laurent et al., 2011). Yet, in many cases, such prior information is absent. One way to circumvent this problem is to perform preliminary ecological niche modelling (ENM) analyses to restrict and/or inform the set of geographical range scenarios to test (Knowles et al., 2007; Carstens & Richards, 2007; Knowles & Alvarado-Serrano, 2010). Ecological niche modelling methods aim to model species distribution by estimating the relationship between species records at sites and the environmental and/or spatial characteristics of those sites (Franklin, 2009; Elith et al., 2011). In bumblebees, it has been suggested that determining habitat suitability, which is analogous to the ENM approach, can be used as a baseline approach for understanding species-specific responses to climate change (Herrera et al., 2014; Lecocq et al., 2015a; Rasmont et al., 2015). Coalescence-based genetic models have been combined with ENM in numerous studies (Carstens & Richards, 2007; Knowles & Carstens, 2007; Shepard & Burbrink, 2008, 2009; Solomon et al., 2008; Bagley et al., 2013; Dellicour et al., 2014a). More recently, there have been studies (e.g. Knowles & Alvarado-Serrano, 2010; Dellicour et al., 2014a) that have also combined the ENM with spatially explicit models of coalescence (Ray et al., 2010b; Dellicour et al., 2014b). In such models, the species range is explicitly drawn on a map, which allows implementing more detailed geographical information than in a more classical approach where different portions of the species range are modelled a priori as panmictic populations.

The overall aim of this study was to investigate whether the geographical range of bumblebees experienced similar changes to those of many exclusively temperate species, which have been characterized by a highly fragmented range during cold episodes of the Pleistocene. We employed two complementary approaches: ENM and genetic variation analyses to investigate and compare the recent geographical range history of five common and co-distributed Bombus species. In particular, we (1) used an ENM method to infer the potential distribution of these species today and during the last glaciation and (2) mapped the geographical distribution of DNA sequence variation for three loci. Alternative hypotheses for the changes in species ranges were then formulated based on these preliminary results, and evaluated via computer simulations of population evolution using a spatially explicit model of coalescence.
### 2.1. Selected study system and genetic data

We selected five common and co-distributed (West) Palaearctic Bombus species: B. (Megabombus) hortorum, B. (Melanobombus) lapidarius, B. (Thoracobombus) pascuorum, B. (Pyrobombus) pratorum and B. (Bombus) terrestris, for which genetic data were already available. Fragments of the mitochondrial gene COI (cytochrome oxidase I) and two protein-coding nuclear genes, EF-1a (elongation factor-1 alpha, F2 copy) and PEPCK (phosphoenolpyruvate carboxykinase), were previously sequenced (total of ~2750 base pairs) for individuals sampled across the distribution of each species by Lecocq et al. (2013a, 2015b) for B. terrestris, Lecocq et al. (2013b, 2015c) for B. lapidarius, Dellicour et al. (2015) for B. hortorum and B. pratorum, and Lecocq et al. (2015d) for B. pascuorum (see Fig. 1 for the sampling locations for each species). Here, we focused on the West Palaearctic portion of their ranges as it corresponded to (1) the co-distribution area of the five species and (2) the area that was thoroughly sampled for each species.

### 2.2. Inferring the past and present species distributions

We estimated the current and Last Glacial Maximum (LGM) West Palaearctic distributions of several Bombus species by using the maximum entropy method implemented in Max-ent 3.3.3 (Phillips et al., 2006; Phillips & Dudík, 2008). Maxent software calculates the environmental suitability of the species, and thus a potential geographical distribution of the species, by minimizing the relative entropy between two probability densities: one estimated from the species occurrence data and one estimated from the landscape (Elith et al., 2011). Here, we estimated each model parameters using the relationship between the present occurrences of the species, which was filtered to eliminate the environmental biases (see below), and some climatic data available for the West Palaearctic area.

We downloaded the 19 bioclimatic variables (Bio1-Bio19) for the current time (1950-2000) at a 2.5 arc-minutes (c. 5 km) resolution, from the WorldClim database (World-Clim 1.4; Hijmans et al., 2005). We cropped the variables to the European range, longitude limits from -25 to 90° and latitude limits from 30 to 80°. We calculated the correlations between these bioclimatic variables using the Pearson correlation coefficient to avoid the use of significantly correlated variables (where the Pearson coefficient > 0.9). By doing so, we selected 13 variables: Bio1 (annual mean temperature), Bio2 (mean diurnal range), Bio3 (isothermality), Bio5 (maximum temperature of warmest month), Bio6 (minimum temperature of coldest month), Bio7 (temperature annual range), Bio8 (mean temperature of wettest quarter), Bio9 (mean temperature of driest quarter), Bio15 (precipitation seasonality), Bio16 (precipitation of wettest quarter), Bio17 (precipitation of driest quarter), Bio18 (precipitation of warmest quarter), Bio19 (precipitation of coldest quarter). We discarded Bio4 (temperature seasonality), Bio10 (mean temperature of warmest quarter), Bio11 (mean temperature of coldest quarter), Bio12 (annual precipitation), Bio13 (precipitation of wettest month) and Bio14 (precipitation of driest month).

We selected the current species occurrences to calibrate our models from the BDFGM database (Base de données fauniques Gembloux-Mons): 16,916 records for B. hortorum, 21,356 for B. lapidarius, 32,115 for B. pascuorum, 17,720 for B. pratorum and 21,262 for B. terrestris. These records mostly correspond to exact locations with a resolution of ~10 km or higher, except for some old observations and/or observation in the eastern part of the distribution for which the coordinates refer to the nearest locality (i.e. town or village). As there are fewer occurrence data for the eastern part of the species range, we decided to keep these records in order to avoid further increasing the sampling bias. Indeed, one potential problem with these data sets is the large difference in sampling effort across the species distribution (e.g. Western Europe versus Eastern Europe). A fundamental assumption in Maxent is that the entire area of interest has been systematically sampled. However, in practice, occurrence records are spatially biased towards better-surveyed areas (Kramer-Schadt et al., 2013). Furthermore, Yackulic et al. (2013) have estimated that 87% of the studies that used Maxent between 2008 and 2012 are likely to suffer from sampling bias. This is important because sampling biases can strongly affect the performance of ecological niche models (Varela et al., 2014a). The problem resides in the fact that field records can be aggregated, leading to the over-sampling (and respectively under-sampling) of some specific ecological conditions. It has been demonstrated that removing occurrence data (field records) that are aggregated in the environmental space can improve ENM performance (Varela et al., 2014b). To avoid this problem, we have employed a climatic filter (Varela et al., 2014b) to select occurrence data for calibrating our models. Specifically, we have selected four representative bioclimatic variables among the 13 variables listed above to filter the species occurrences: maximum temperature of the warmest month (Bio5), minimum temperature of the coldest month (Bio6), mean temperature of the warmest quarter (Bio10), and the annual precipitation (Bio12). This resulted in the selection of 133 records for B. hortorum, 102 for B. lapidarius, 134 for B. pascuorum, 128 for B. pratorum and 129 for B. terrestris. The geographical distribution of the species occurrence records used for ENM is displayed in Fig. 1.

We performed ten replicates for each analysis in Maxent. We used the default convergence threshold (10-5), 10,000 iterations and a 'random seed' to generate a random partition of our records into training (90%) and test (10%) records. These random partitions were used to test the model. Model results were evaluated using the AUC (area under the curve) - area under the ROC (receiver operating characteristic) curve - commonly used to assess the Maxent estimation performance (Lobo et al., 2008). All our models had a high discriminative power for the training data sets (AUC > 0.89) and they were also able to predict the testing points (AUC > 0.82; see Table S1 in Supporting Information for further details). Subsequently, we projected the models into two different LGM climatic scenarios, CCSM (Community Climate System Model) and MIROC (Model for Interdisciplinary Research on Climate), to estimate the geographical ranges of the species during this cold climatic event. The bioclimatic variables from the LGM climatic scenarios of the two models (CCSM and MIROC) were based on the models from PMIP2 database (Paleoclimate Modelling Intercomparison Project Phase II, Braconnot et al., 2007) and we also downloaded them from the WorldClim database (Hijmans et al., 2005).
## 2.3. Comparing and mapping intra-specific genetic variation

Patterns of genetic variability were investigated at two geographical levels: (1) at a global scale by considering genetic diversity and differentiation (i.e. phylogeographical signal) at the level of the entire range for each species, and (2) at a local scale by mapping spatial variation in genetic diversity within populations and in genetic differentiation among populations (see below). The following statistics were estimated using the toolbox SPADS 1.0 (Dellicour & Mardulyn, 2014) to compare the genetic variability in each species at a global scale: nucleotide diversity $p$ (Nei & Li, 1979) over the entire range of each species and phylogeographical structure within species as measured by NST-GST (Pons & Petit, 1996), separately for each locus. These statistics were selected because they take into account differences in both allelic frequencies and genetic distances (mismatches between DNA sequences). The statistical significance of the difference between NST and GST was evaluated using 10,000 random permutations of haplotypes in the original data sets.

To study genetic variation at the local scale, we used the R functions GDisPAL and GDivPAL (Dellicour & Mardulyn, 2014) to generate graphical representations of genetic variation across species ranges. These functions generated graphs displaying geographical variation in genetic diversity and genetic differentiation by extending a method initially developed by Miller (2005) and based on an interpolation procedure (inverse distance-weighted interpolation; Watson & Philips, 1985; Watson, 1992). For genetic diversity, interpolation was based on diversity values directly estimated at each sampling locality, but for genetic differentiation, interpolation was based on distance values assigned at mid-points of each edge of a connectivity network built between the sampling localities (i.e. a Delaunay triangulation). For each species, we generated two differentiation surfaces based on different inter-individual distances and three surfaces based on distinct genetic diversity indices. All diversity and distance measures used are independent of sample size (which is different between species) and were estimated with SPADS. The two inter-individual distances were (1) an inter-individual distance based on allelic frequencies as defined by Miller (2005), and (2) an inter-individual distance based on distance DNA sequence mismatches averaged over all loci (i.e. average P-distances across loci, see the SPADS manual for further details). The three diversity statistics were (1) allelic richness $AR$ of each population (estimation of the expected number of different haplotypes in a subsample of $n$ sequences for a given population, $n$ being the smallest number of sequences obtained for a sampling site; El Mousadik & Petit, 1996), (2) nucleotide diversity $p$ of each population (average number of nucleotide differences per site between any two DNA sequences chosen randomly in a given population; Nei & Li, 1979) and (3) relative nucleotide diversity $pR$ of each population (nucleotide diversity within this population divided by the nucleotide diversity within the group formed by all other populations; Mardulyn et al., 2009). Note that for the allelic richness estimation, as $n = 2$ for all five species, interpolation surfaces based on this diversity measure were directly comparable. All surfaces were generated using a distance weighting parameter $a = 1, 2$ and $5$, and differentiation surfaces were both based on a Delaunay triangulation connectivity network. Furthermore, we performed the distance interpolations using residual distances derived from the linear

regression of genetic versus geographical distances (Manni et al., 2004). Geographical distances were measured as great circle geographical distances (i.e. distances on the surface of the earth in kilometers).

## 2.4. Spatially explicit coalescent simulations

While the estimates of past species distributions for the LGM obtained in this study with Maxent mainly suggest a continuous distribution, or at least strong connectivity among potential refuge areas, our mapping of genetic diversity for the five species, along with genetic data from previous studies (Lecocq et al., 2013b; Dellicour et al., 2015), suggests relatively strong population structure across their ranges, more compatible with a long history of fragmented distributions in Europe (see Results). Given this apparent contradiction between ENM and genetic variation results, we designed and compared four evolutionary scenarios (Fig. 2), by varying the level of fragmentation of the species range during the last glaciation (110–12 ka). The first scenario (S1), directly inspired by the 1VIAXENT results, postulates a largely continuous range during the last glaciation. The second scenario (S2) was identical to S1, but for the isolation of the Anatolia-Caucasian region during the LGM. A strong past isolation of this region was indeed already suggested by results from Dellicour et al. (2015; for B. hortorum and B. pratorum). Scenarios S3 and S4 were designed to model a much more fragmented range during the last glaciation with distinct isolated refuge areas. For this period, scenario S3 depicts four isolated regions corresponding to the three classically defined European glacial refuges (Iberian Peninsula, Italian Peninsula and the Balkans; e.g. Taberlet et al., 1998), plus the Anatolia-Caucasian region. Scenario S4 was similarly defined, but included an additional refuge located in Central Europe, whose existence was suggested for B. lapidarius in a previous study (Lecocq et al., 2013b).

DNA sequences were simulated along these evolutionary scenarios using a spatially explicit coalescence model as implemented in PHYLOGEOSIM 1.0 (Dellicour et al., 2014b). In this model, the geographical structure of the studied range at any point in time is defined by a two-dimensional grid in which each cell is considered as a panmictic population. Cells accessible to individuals are specified, and each sampled sequence is attributed to one cell on the grid displaying the current species range ($t = 0$). The simulation begins at the current time $t = 0$ and proceeds backward in time until all gene copies (DNA sequences) for the locus considered have coalesced (Dellicour et al., 2014b). At each generation, a given gene copy has the opportunity to move to adjacent cells (migration) and to coalesce with another gene copy located in the same cell. At the end of the simulation, a genealogy is reconstructed based on the recorded coalescence events (Dellicour et al., 2014b). The total number of mutations added on the genealogy is specified a priori and inferred from the empirical data. Note that in the framework implemented in PHYLOGEOSIM, a forward pre-simulation is performed to estimate demographic parameters such as cell population effective size and migration rates [See figure 4 in Dellicour et al. (1014b) for a representation of the general workflow].

To define spatially explicit scenarios in PHYLOGEOSIM, matrices of maximal cell population effective sizes were built according to each of the four distinct scenarios by superim-
posing a grid of 130 9 75 cells to a current map of Europe or to a map of the same area at the LGM. As we do not have access to prior information regarding the inter-cells migration rate and cell population effective size, we tested a large range of parameter values. For each scenario, we tested all combinations of the following parameter values: (1) forward migration rate between adjacent cells on the grid (10-3, 10-4 and 10-5), (2) maximal effective population sizes for grid cells that sequences are allowed colonizing (1000; 10,000 and 100,000 for the two nuclear loci and 333; 3333 and 33,333 for the COI mitochondrial locus, as haploid mitochondria are expected to have an effective population size three times smaller than that of the nuclear genome in haplodiploid Hymenoptera organisms; Hedrick & Parker, 1997) and (3) reproduction rate within a cell (2 and 5). Note that some matrices of maximal effective population sizes include some 'connecting' cells (green cells in Fig. 2) with a maximal cell population effective size systematically 10 times lower than that of the other accessible cells (orange cells in Fig. 2). As displayed in Fig. 2, these cells allowed potential connections between islands and mainland, and their lower effective size was controlled for restricted migrations between these areas. We performed a total of 100 independent backward simulations for each individual scenario. To take into account stochastic variation associated with the forward simulation, we reiterated a forward simulation for every ten backward simulations to estimate the backward simulation parameters (backward migration rates and effective population sizes), yielding a total of 10 forward simulations. The sequence lengths and a fixed number of mutations were specified for each simulated locus, matching the parameters of our real DNA sequence alignments.

The historical scenarios we have tested were all divided into four separate geographical layers (Fig. 2). The first layer (going backward in time) corresponded to the current interglacial period (from now to 12 ka), the second to the last glacial period (from 12 to 110 ka), the third to the previous interglacial period (from 110 to 130 ka), and the fourth to the penultimate glacial period (from 130 ka to the time to the most recent common ancestor). With the exception of B. terrestris, for

which more than one generation per year was recorded on a portion of its distribution (Rasmont et al., 2008), the considered Bombus species are all univoltine. We therefore assumed that 1 year approximately corresponded to one generation. Using this framework, we designed the four scenarios introduced above. The comparison of these four scenarios aimed at investigating connectivity among regions during the last glacial period and, among others, to evaluate the possibility of a large continuous range, as suggested by ENM results.

The simulated data sets were compared to real data sets, to evaluate the relative probabilities of the proposed scenarios, through calculation of different summary statistics. The comparisons were based, for the five species, on summary statistics measuring genetic diversity and genetic differentiation: (1) nucleotide diversity (p) (Nei & Li, 1979) estimated within each defined group of sampled populations and (2) the three AMOVA (analysis of molecular variance; Excoffier et al., 1992) statistics ΦSC, ΦST and ΦCT estimated from the same group partition. By definition, ΦSC measures the proportion of variation among populations within groups, ΦST the proportion of variation among populations and ΦCT the proportion of variation among groups. For estimating these summary statistics, we thus defined several groups of sampled populations: Iberian Peninsula, Italy, Balkans, Anatolia-Caucasus, Central Europe and a sixth group gathering the remaining sampled populations. Summary statistics were automatically computed by PHYLOGEOSIM on simulated data sets and were computed by SPADS on real data sets. We combined all summary statistics computed for one data set into a single v2 statistic (Dellicour et al., 2014a). For a given scenario and a given set of parameters, we were thus able to generate a distribution of v2 statistics, to which we compared the v2 statistic estimated on the corresponding real data set. This comparison yielded a nonparametric test returning a P-value equal to the proportion of simulated values greater than or equal to the observed value. Finally, P-values inferred for each locus were combined using the method of Fisher (1948) to provide a unique and global probability value for each scenario and set of simulation parameter values.

# 3. Results

## 3.1. Inference of current and past geographical distributions

All our MAXENT models have a high discriminative power for the training data sets (AUC > 0.89) and they are also able to predict the testing points (AUC > 0.82; see Appendix S1 for all the AUC values). The MAXENT predicted ranges are different for the five species, but also according to the climatic model used for the LGM (Fig. 3): CCSM and MIROC. Yet, all inferred distributions suggest a continuous range during the last glaciation, with the notable exception of the Caucasian region, which appears relatively isolated on LGM predictions based on CCSM for B. hortorum, B. lapi-darius, B. pascuorum and B. pratorum. Therefore, except for the Anatolia-Caucasian region, ENM analyses performed with MAXENT do not tend to support a hypothesis of fragmented range for these species in the past (except for B. hor-torum where the LGM distribution estimate is based on MIROC). It is noteworthy that high probabilities of occurrence are inferred even at the positions of the three main ice sheets (i.e. Alps, British and Fennoscandian ice sheets; see Toucanne et al., 2010), which cover a large portion of Europe during the LGM.

## 3.2. Comparison of intra-specific genetic variation

Overall nucleotide diversity and phylogeographical signal estimated within each species are reported in Table 1, and surfaces generated with distance weighting parameter a = 1, 2 and 5 are displayed in Appendix S2 (in Figs S1–S3 respectively). Comparison of genetic diversity measured over the entire range of each species reveals divergent levels of nucleotide diversity among species: B. hortorum presents the smallest nucleotide diversity for all genetic markers, while the highest COI and nuclear nucleotide diversities are respectively found in B. terrestris and B. pascuorum (Table 1). B. hortorum is associated with the highest phylogeographical structure estimated for COI, whereas significant phylogeographical structures for all three markers are only detected in B. lapidarius and B. pascuorum (NST-GST; Table 1). Although we observe differences among interpolation surfaces (see Figs S1–S3 in Appendix S2), it is difficult to highlight meaningful overall trends from the comparison among species of patterns of genetic diversity, or population differentiation, at the local level. Bombus pratorum is associated with the highest
genetic diversity regardless of the diversity index used. Bombus hortorum displays a low relative nucleotide diversity (pR) for most localities but the distribution of all four other species shows at least some areas of higher relative nucleotide diversity. The comparison of nucleotide diversity (p) and allelic richness (AR) does not reveal any clear difference between B. hortorum and the other species. Finally, differentiation surfaces based on the two inter-individual measures show comparable levels among adjacent populations for all species.

### 3.3. Spatially explicit coalescent simulations

Results of the comparisons between real and simulated datasets are summarized in Table 2 (see also Tables S1–S5 in Appendix S3 for results detailed locus by locus). For B. hortorum, except for one particular set of simulation parameters (cell population effective size = 333/1000; forward

migration rate = 0.0001 and reproduction rate = 2 or 5), scenario S3 (i.e. four isolated regions) is systematically associated with the highest combined P-value when comparing data sets simulated under the same cell population effective size and forward migration rate values. Combined P-values obtained for B. lapidarius are all smaller than 0.02 and in most cases even smaller than 0.01. As a consequence, all tested scenarios appear incompatible with the data available for this species. On the contrary, for B. pratorum, all scenarios are at least once associated with a P-value higher than 0.05 (with similar values ranging from 0.05 to 0.13). The highest P-value (P = 0.16 and 0.07) obtained for the other species, B. ter-restris and B. pascuorum, respectively, is again associated with scenario S3, although for the same simulation parameters, a P-value > 0.05 is also found for scenario S2 in the case of B. terrestris (P = 0.08). Note that in this context, each P-value represents the estimated probability that the simulated scenario has generated the observed pattern of genetic variation.

## 4. Discussion

### 4.1. Recent geographical range history of West Palaearctic bumblebee species

Bumblebees are endothermic organisms (Heinrich, 1979), consequently we expect, as the ENM predictions suggest, that the colder climate of the last glaciation had only a minor impact on their distribution, which may have remained mostly continuous at the time. Indeed, the five studied species are polylectic (i.e. host-plant generalist) and currently display a widespread range, occurring from the south (Mediterranean area) to the north of Europe (B. hortorum, B. pascuorum, and B. pratorum reach the Arctic Ocean). Thus, they are likely to have a large adaptive potential to eco-climatic changes, in particular to global cooling. However, our results based on coalescent simulations contradict the hypothesis of a large connected range occurring at the LGM for three out of five species (i.e. B. hortorum, B. pascuorum and B. terrestris). Indeed, scenario S3 (classical hypothesis of a fragmented range with isolated regions in the Iberian Peninsula, Italy, Balkans and Anatolia) often identified for other organisms, including some vertebrates and insects (e.g. Hewitt, 2004; Stewart et al., 2010), was associated with the highest P-value for these Bombus species. The phenotypic differentiation observed within these species (i.e. several allopatric subspecies described, especially in the southern peninsulas and Anatolia; Rasmont et al., 2008; Lecocq et al., 2015d) further supports a history of intraspecific fragmentation.

For the two other species, B. lapidarius and B. pratorum, comparisons between real and simulated data sets could not discriminate among the four tested scenarios. For B. lapidarius, none of the tested sets of simulation parameters led to simulated data sets associated with P-values higher than 0.02. This suggests that the set of scenarios here tested were not sufficiently representative of the true geographical range history. In the case of B. pratorum, the relatively similar P-values associated with all four scenarios makes it difficult to draw any strong conclusions. Consequently, this suggests that the available genetic variation signal for this species is not sufficiently informative. Yet, results obtained for B. lapidarius and B. pratorum do not contradict the convergent results obtained for the three other species.

### 4.2. Ecological niche modelling limitations

Our results based on coalescent simulations contradict those of ENM analyses for three out of five species. A similar observation has been observed in Dellicour et al. (2014a), highlighting the limitations of the ENM approach for estimating current and past species ranges. From a practical point of view, a main issue in ENM based range estimates is the availability of past environmental data. Indeed, taking only temperature and humidity information into account, as is the case in this study, could lead to erroneous estimates. For example, ENM here suggests the potential occurrence of a species over areas that were presumably covered by ice sheets (e.g. Alps, British and Fennoscandian ice sheets; Toucanne et al., 2010) during the LGM (Fig. 3). Algorithms used to model the climatic requirements of the species, like Max-ent, allows fitting linear relationships with the temperature and precipitation variables. However, having no lower limits of tolerance for the species implies that the estimated range maps overlap extreme cold climatic conditions like, for example, areas covered by ice sheets. Beside climatic requirements, such ENM methods also do not take into account potential range limitation induced by biotic interactions (Schweiger et al., 2008; Hanspach et al., 2014) or related to limited dispersal capacities (reviewed in Rasmont et al., 2015). Such factors could also lead to potential overestimations of ecological niches.

In recent years, the ENM approach has been applied to a large number of studies (Radosavljevic & Anderson, 2014). Among alternative methods available, Maxent has been the most widely used (Fourcade et al., 2014) because it can easily and rapidly provide detailed results on the current and past occurrences of a target species. Yet, as discussed above, it is important to note that a Maxent-estimated distribution represents the potential range rather than the realized range of a species. Although increasingly used, ecological niche models based on climatic data alone to generate hypotheses over the evolution of species ranges (mainly from the LGM to the current) should be interpreted with caution.
<table><caption>Table 1. Overall nucleotide diversity p and phylogeographical signal as measured by NST-GST estimated for five species of Bombus in the West Palaearctic within each species.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Nucleotide diversity π</th><th colspan="3">NST-GST</th></tr><tr><th>COI</th><th>EF-1α</th><th>PEPCK</th><th>COI</th><th>EF-1α</th><th>PEPCK</th></tr></thead><tbody><tr><td><i>B. hortorum</i></td><td>0.00226</td><td>0.00000</td><td>0.00021</td><td>0.546*</td><td>—</td><td>0.000</td></tr><tr><td><i>B. lapidarius</i></td><td>0.00483</td><td>0.00131</td><td>0.00091</td><td>0.221*</td><td>0.125*</td><td>0.067*</td></tr><tr><td><i>B. pascuorum</i></td><td>0.00510</td><td>0.00240</td><td>0.00679</td><td>0.273*</td><td>0.196*</td><td>0.186*</td></tr><tr><td><i>B. pratorum</i></td><td>0.00318</td><td>0.00143</td><td>0.00462</td><td>0.214*</td><td>0.061</td><td>0.150*</td></tr><tr><td><i>B. terrestris</i></td><td>0.01450</td><td>0.00000</td><td>0.00133</td><td>0.123*</td><td>—</td><td>0.107</td></tr></tbody></table>

<footnote>(*) indicates significant NST-GST value (P-value < 0.05). Results for Bombus hortorum and B. pratorum come from Dellicour et al. (2015).</footnote>
<table><caption>Table 2. Combined P-values based on 100 replicates per set of parameters and obtained from the comparison between real and simulated data sets</caption><thead><tr><th rowspan="2">Simulation parameters<br>Ne (mitoch./nuclear)</th><th rowspan="2">fm</th><th rowspan="2">Scena<br>rio</th><th colspan="2">B. hortorum</th><th colspan="2">B. lapidarius</th><th colspan="2">B. pascuorum</th><th colspan="2">B. pratorum</th><th colspan="2">B. terrestris</th></tr><tr><th>t<sub>R</sub> = 2</th><th>t<sub>R</sub> = 5</th><th>t<sub>R</sub> = 2</th><th>t<sub>R</sub> = 5</th><th>t<sub>R</sub> = 2</th><th>t<sub>R</sub> = 5</th><th>t<sub>R</sub> = 2</th><th>t<sub>R</sub> = 5</th><th>t<sub>R</sub> = 2</th><th>t<sub>R</sub> = 5</th></tr></thead><tbody><tr><td>333/1000</td><td>0.0001</td><td>S1</td><td><strong>0.06</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.0001</td><td>S2</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.0001</td><td>S3</td><td>&lt; 0.01</td><td><strong>0.05</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.0001</td><td>S4</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.001</td><td>S1</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.001</td><td>S2</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.001</td><td>S3</td><td>&lt; 0.01</td><td><strong>0.05</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>333/1000</td><td>0.001</td><td>S4</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.00001</td><td>S1</td><td>0.02</td><td>0.04</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.00001</td><td>S2</td><td>0.01</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.00001</td><td>S3</td><td><strong>0.06</strong></td><td>0.04</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.00001</td><td>S4</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.0001</td><td>S1</td><td>0.02</td><td>&lt; 0.01</td><td><strong>0.09</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.01</td><td>0.02</td><td><strong>0.10</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.0001</td><td>S2</td><td>0.04</td><td><strong>0.09</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.07</strong></td><td>0.04</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.0001</td><td>S3</td><td>0.02</td><td><strong>0.19</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.07</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.0001</td><td>S4</td><td><strong>0.05</strong></td><td>0.02</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.01</td><td>&lt; 0.01</td><td>0.02</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.001</td><td>S1</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.11</strong></td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.001</td><td>S2</td><td>&lt; 0.01</td><td><strong>0.05</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.13</strong></td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.001</td><td>S3</td><td><strong>0.07</strong></td><td>0.04</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.10</strong></td><td>&lt; 0.01</td></tr><tr><td>3333/10,000</td><td>0.001</td><td>S4</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.05</strong></td><td>0.04</td></tr><tr><td>33,333/100,000</td><td>0.00001</td><td>S1</td><td>0.04</td><td>&lt; 0.01</td><td>0.01</td><td>&lt; 0.01</td><td>0.04</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.00001</td><td>S2</td><td>0.04</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.00001</td><td>S3</td><td><strong>0.09</strong></td><td><strong>0.06</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.03</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.00001</td><td>S4</td><td>0.01</td><td>0.02</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.02</td><td>0.04</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.0001</td><td>S1</td><td>&lt; 0.01</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.04</td><td><strong>0.09</strong></td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.0001</td><td>S2</td><td>0.04</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.08</strong></td><td>0.02</td></tr><tr><td>33,333/100,000</td><td>0.0001</td><td>S3</td><td><strong>0.07</strong></td><td><strong>0.15</strong></td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.16</strong></td><td>0.03</td></tr><tr><td>33,333/100,000</td><td>0.0001</td><td>S4</td><td>0.02</td><td>0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.001</td><td>S1</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.001</td><td>S2</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.001</td><td>S3</td><td>0.04</td><td>0.03</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td><strong>0.08</strong></td><td>0.03</td><td>&lt; 0.01</td></tr><tr><td>33,333/100,000</td><td>0.001</td><td>S4</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr></tbody></table>

Combined P-values based on 100 replicates per set of parameters and obtained from the comparison between real and simulated data sets. Ne refers to the maximal effective population sizes, fm to the forward migration rate and t<sub>R</sub> to the reproduction rate. Summary statistics used for the comparison between real and simulated data sets: (1) the nucleotide diversity Π of each defined group of sampled populations (see the text) and (2) the three AMOVA statistics (Φ<sub>SC</sub>, Φ<sub>ST</sub> and Φ<sub>CT</sub>) estimated from this same group partition. See Tables S2-6 for the P-values obtained for each locus taken separately. Values in bold correspond to combined P-values > 0.05
<figure><img src="image_3.png" /><figcaption>Figure 1. West Palaearctic maps with Bombus hortorum, B. lapidarius, B. pascuorum, B. pratorum and B. terrestris occurrence data used in ecological niche modelling analyses (above) and maps of sampling localities for which DNA sequence data (one mitochondrial and two nuclear loci) were available and analysed in this study.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Four alternative evolutionary scenarios simulated with PHYLOGEOSIM 1.0 for Bombus hortorum (B. hor). Each scenario includes chronological layers (the oldest on the left) that show potential presence of the species for each grid cell on the map. Orange cells are those accessible to gene copies during the simulation. Green cells are connecting cells (allowing connection between separated regions), also accessible to gene copies, but with a maximal effective size systematically 10 times lower than the effective size of orange cells. The time interval (in years, going backward in time) corresponding to each layer is indicated. (b) Four layers specific to each of the other species: B. lapidarius (B. lap), B. pascuorum (B. pas), B. pratorum (B. pra) and B. terrestris (B. ter). Scenarios associated with these four species can be reconstructed simply by replacing layers '*B. hor' by the corresponding other layer.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3. Results of ecological niche modelling: inferred European distributions of five Bombus species obtained with MAXENT and based on current climatic data ('current times') or reconstructed Last Glacial Maximum (LGM) palaeoclimatic data ('CCSM' and 'MIROC'). Probabilities of occurrence are displayed using different colours.</figcaption></figure>

5. References

Andersen, B.G. & Borns, H.W., Jr (1994) The Ice Age world. Scandinavian University Press, Oslo.

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

Bagley, J.C., Sandel, M., Travis, J., Lozano-Vilano, M.D.L. & Johnson, J.B. (2013) Paleoclimatic modeling and phylo-geography of least killifish, Heterandria formosa: Insights into Pleistocene expansion-contraction dynamics and evolutionary history of North American Coastal Plain freshwater biota. BMC Evolutionary Biology, 13, 223.

Braconnot, P., Otto-Bliesner, B., Harrison, S. et al. (2007) Results of PMIP2 coupled simulations of the Mid-Holocene and last glacial maximum -Part 1: Experiments and large-scale features. Climate of the Past, 3, 261-277.

Carstens, B.C. & Richards, C.L. (2007) Integrating coalescent and ecological niche modeling in comparative phylogeography. Evolution, 61, 1439-1454.

Dehon, M., Michez, D., Nel, A., Engel, M.S. & De Meule-meester, T. (2014) Wing shape of four new bee fossils (Hymenoptera: Anthophila) provides insights to bee evolution. PLoS ONE, 9, e108865.

Dellicour, S. & Mardulyn, P. (2014) SPADS 1.0: A toolbox to perform spatial analyses on DNA sequence data sets. Molecular Ecology Resources, 14, 647-651.

Dellicour, S., Fearnley, S., Lombal, A., Heidl, S., Dahlhoff, E.P., Rank, N.E. & Mardulyn, P. (2014a) Inferring the past and present connectivity across the range of a North American leaf beetle: combining ecological niche modeling and a geographically explicit model of coalescence. Evolution, 68, 2371-2385.

Dellicour, S., Kastally, C., Hardy, O.J. & Mardulyn, P. (2014b) Comparing phylogeographic hypotheses by simulating DNA sequences under a spatially explicit model of coalescence. Molecular Biology and Evolution, 31, 3359-3372.

Dellicour, S., Michez, D. & Mardulyn, P. (2015) Comparative phylogeography of five bumblebees: impact of range fragmentation, range size and diet specialization. Biological Journal of the Linnean Society, 116, 926-939.

El Mousadik, A. & Petit, R.J. (1996) High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) Skeels] endemic to Morocco. Theoretical and Applied Genetics, 92, 832-839.

Elith, J., Phillips, S.J., Hastie, T., Duda, R., M., Chee, Y.E. & Yates, C.J. (2011) A statistical explanation of MaxEnt for ecologists. Diversity and Distributions, 17, 43-57.

Estoup, A., Solignac, M., Cornuet, J.M., Goudet, J. & Scholl, A. (1996) Genetic differentiation of continental and island populations of Bombus terrestris (Hymenoptera: Apidae) in Europe. Molecular Ecology, 5, 19-31.

Excoffier, L., Smouse, P.E. & Quattro, J.M. (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics, 131, 479-491.

Fagundes, N.J.R., Ray, N., Beaumont, M., Neuenschwander, S., Salzano, F.M., Bonatto, S.L. & Excoffier, L. (2007) Statistical evaluation of alternative models of human evolution. Proceedings of the National Academy of Sciences of the USA, 104, 17614-17619.

Fisher, R.A. (1948) Questions and answers #14. The American Statistician, 2, 30-31.

Fourcade, Y., Engler, J.O., Reodder, D. & Secondi, J. (2014) Mapping species distributions with MAXENT using a geographically biased sample of presence data: a performance assessment of methods for correcting sampling bias. PLoS ONE, 9, e97122.

Franklin, J. (2009) Mapping species distributions: spatial inference and prediction. Cambridge University Press, Cambridge, UK.

Hanspach, J., Schweiger, O., Kühn, I., Plattner, M., Pearman, P.B., Zimmermann, N.E. & Settele, J. (2014) Host plant availability potentially limits butterfly distributions under cold environmental conditions. Ecography, 37, 301-308.

Hedrick, P.W. & Parker, J.D. (1997) Evolutionary genetics and genetic variation of haplotype and X-linked genes. Annual Review of Ecology and Systematics,

28, 55-83.

Heinrich, B. (1979) Bumblebee economics. Harvard University Press, Cambridge, MA.

Herrera, J.M., Ploquin, E.F., Rodríguez-Pérez, J. & Obeso, J.R. (2014) Determining habitat suitability for bumblebees in a mountain system: a baseline approach for testing the impact of climate change on the occurrence and abundance of species. Journal of Biogeography, 41, 700-712.

Hewitt, G.M. (2004) Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society B: Biological Sciences, 359, 183-195.

Hickerson, M.J., Carstens, B.C., Cavender-Bares, J., Crandall, K.A., Graham, C.H., Johnson, J.B., Rissler, L., Victoriano, P.F. & Yoder, A.D. (2010) Phylogeography's past, present, and future: 10 years after Avise, 2000. Molecular Phylo-netics and Evolution, 54, 291-301.

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.

Hines, H.M. (2008) Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: Bombus). Systematic Biology, 57, 58-75.

Knowles, L.L. (2001) Did the pleistocene glaciations promote divergence? Tests of explicit refugial models in montane grasshoppers. Molecular Ecology, 10, 691-701.

Knowles, L.L. (2009) Statistical phylogeography. Annual Review of Ecology Evolution and Systematics, 40, 593-61

Knowles, L.L. & Alvarado-Serrano, D.F. (2010) Exploring the population genetic consequences of the colonization process with spatio-temporally explicit models: insights from coupled ecological, demographic and genetic models in montane grasshoppers. Molecular Ecology, 19, 3727-3745.

Knowles, L.L. & Carstens, B.C. (2007) Estimating a geographically explicit model of population divergence. Evolution, 61, 477-493.

Knowles, L.L., Carstens, B.C. & Keat, M. (2007) Coupling genetic and ecological-niche models to examine how past population distributions contribute to divergence. Current Biology, 17, 940-946.

Kramer-Schadt, S., Niedballa, J., Pilgrim, J.D. et al. (2013) The importance of correcting for sampling bias in MaxEnt species distribution models. Diversity and Distributions, 19, 1366-1379

Laurent, S.J.Y., Werzner, A., Excoffier, L. & Stephan, W. (2011) Approximate Bayesian analysis of Drosophila melanogaster polymorphism data reveals a recent colonization of Southeast Asia. Molecular Biology and Evolution, 28, 2041-2051.

Lecocq, T., Vereecken, N.J., Michez, D., Dellicour, S., Lhomme, P., Valterovaia, I., Rasplus, J.Y. & Rasmont, P. (2013a) Patterns of genetic and reproductive traits differ from mainland vs. Corsican populations of bumblebees. PLoS One, 8, e65642.

Lecocq, T., Dellicour, S., Michez, D., Lhomme, P., Vander-Planck, M., Valterovaia, I., Rasplus, J.Y. & Rasmont, P. (2013b) Scent of a break-up: phylogeography and reproductive trait divergences in the red-tailed bumblebee (Bombus lapidarius). BMC Evolutionary Biology, 13, 263.

Lecocq, T., Rasmont, P., Harpke, A. & Schweiger, O. (2015a) Improving international trade regulation by considering intraspecific variation for invasion risk assessment of commercially traded species: the Bombus terrestris case. Conservation Letters. doi: 10.1111/conl.12215.

Lecocq, T., Brasero, N., De Meulemeester, T., Michez, D., Dellicour, S., Lhomme, P., de Jonghe, R., Valterovaia, I., Urbanovaia, K. & Rasmont, P. (2015b) An integrative taxonomic approach to assess the status of Corsican bumblebees: implications for conservation. Animal Conservation, 18, 236-248.

Lecocq, T., Dellicour, S., Michez, D., Dehon, M., Dewulf, A., De Meulemeester, T., Brasero, N., Valterovaia, I., Rasplus, J.-Y. & Rasmont, P. (2015c) Methods for species delimitation in bumblebees (Hymenoptera, Apidae, Bombus): towards an integrative
approach. Zoologica Scripta, 44, 281-297.

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

Lobo, J.M., Jiménez-valverde, A. & Real, R. (2008) AUC: A misleading measure of the performances of predictive distribution models. Global Ecology and Biogeography, 17, 145-151.

Manni, F., Guerard, E. & Heyer, E. (2004) Geographic patterns of (genetic, morphologic, linguistic) variation: How barriers can be detected by using Monmonier's algorithm. Human Biology, 76, 173-190.

Mardulyn, P., Mikhailov, Y.E. & Pasteels, J.M. (2009) Testing phylogeographic hypotheses in a Euro-Siberian cold-adapted leaf beetle with coalescent simulations. Evolution, 63, 2717-2729.

Metcalf, J.L., Prost, S., Nogues-Bravo, D., Dechaine, E.G., Anderson, C., Batra, P., Araujo, M.B., Cooper, A. & Gural-nick, R.P. (2014) Integrating multiple lines of evidence into historical biogeography hypothesis testing: a Bison bison case study. Proceedings of the Royal Society B: Biological Sciences, 281, 20132782.

Miller, M.P. (2005) Alleles In Space (AIS): computer software for the joint analysis of interindividual spatial and genetic information. Journal of Heredity, 96, 722-724.

Nei, M. & Li, W.H. (1979) Mathematical model for studying genetic variation in terms of restriction endonucleases. Proceedings of the National Academy of Sciences of the USA, 76, 5269-5273.

Phillips, S.J. & Dudík, M. (2008) Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 31, 161-175.

Phillips, S.J., Anderson, R.P. & Schapire, R.E. (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190, 231-259.

Pons, O. & Petit, R.J. (1996) Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics, 144, 1237-1245.

Radosavljevic, A. & Anderson, R.P. (2014) Making better Maxent models of species distributions: complexity, overfitting and evaluation. Journal of Biogeography, 41, 629-643.

Rasmont, P., Coppée, A., Michez, D. & De Meulemeester, T. (2008) An overview of the Bombus terrestris (L. 1758) subspecies (Hymenoptera: Apidae). Annales de La Société Entomologique De France, 44, 243-250.

Rasmont, P., Franzén, M., Lecocq, T. et al. (2015) Climatic risk and distribution atlas of European bumblebees. Biodiversity and Ecosystem Risk Assessment, 10, 1-236.

Ray, N., Wegmann, D., Fagundes, N.J.R., Wang, S., Ruiz-Linares, A. & Excoffier, L. (2010a) A statistical evaluation of models for the initial settlement of the American con- tinent emphasizes the importance of gene flow with Asia. Molecular Biology and Evolution, 27, 337-345.

Ray, N., Currat, M., Foll, M. & Excoffier, L. (2010b) SPLATCHE2: a spatially explicit simulation framework for complex demography, genetic admixture and recombination. Bioinformatics, 26, 2993-2994.

Schweiger, O., Settele, J., Kudrna, O., Klotz, S. & Kěuhn, I. (2008) Climate change can cause spatial mismatch of trophically interacting species. Ecology, 89, 3472-3479.

Shepard, D.B. & Burbrink, F.T. (2008) Lineage diversification and historical demography of a sky island salamander, Plethodon ouachitae, from the Interior Highlands. Molecular Ecology, 17, 5315-5335.

Shepard, D.B. & Burbrink, F.T. (2009) Phylogeographic and demographic effects of Pleistocene climatic fluctuations in a montane salamander, Plethodon fourchensis. Molecular Ecology, 18, 2243-2262.

Solomon, S.E., Bacci, M., Jr, Martins, J., Jr, Vinha, G.G. & Mueller, U.G. (2008) Paleodistributions and comparative molecular phylogeography of leafcutter ants (Atta spp.) provide new insight into the origins of Amazonian diversity. PLoS ONE, 3, e2738.

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.

Taberlet, P., Fumagalli, L., Wust-Saucy, A.G. & Cosson, J.F. (1998) Comparative phylogeography and postglacial colonization routes in Europe. Molecular Ecology, 7, 453-464.

Toucanne, S., Zaragosi, S., Bourillet, J.F., Marieu, V., Cremer, M., Kageyama, M., Van Vliet-Lanoë, B., Eynaud, F., Turon, J.L. & Gibbard, P.L. (2010) The first estimation of Fleuve Manche palaeoriver discharge during the last deglaciation: evidence for Fennoscandian ice sheet meltwater flow in the English Channel ca 20-18 ka ago. Earth and Planetary Science Letters, 290, 459-473.

Varela, S., Lobo, J.M., Rodríguez, J. & Batra, P. (2010) Were the Late Pleistocene climatic changes responsible for the disappearance of the European spotted hyena populations? Hindcasting a species geographic distribution across time. Quaternary Science Reviews, 29, 2027-2035.

Varela, S., Mateo, R.G., García-Valdés, R. & Fernández-González, F. (2014a) Macroecology and ecoinformatics: evaluating the accuracy of the ecological niche models calibrated with species occurrence data with biases and/or errors. Ecosistemas, 23, 46-53.

Varela, S., Anderson, R.P., García-Valdés, R. & Fernández-González, F. (2014b) Environmental filters reduce the effects of sampling bias and improve predictions of ecological niche models. Ecography, 37, 1084-1091.

Veeramah, K.R., Wegmann, D., Woerner, A., Mendez, F.L., Watkins, J.C., Destro-Bisol, G., Soodyall, H., Louie, L. & Hammer, M.F. (2012) An early divergence of KhoeSan ancestors from those of other modern humans is supported by an ABC-based analysis of autosomal resequencing data. Molecular Biology and Evolution, 29, 617-630.

Watson, D.F. & Philips, G.M. (1985) A refinement of inverse distance weighted interpolation. Geo-processing, 2, 315-327.

Watson, D.F. (1992) Contouring: A Guide to the Analysis and Display of Spatial Data. Pergamon Press, New York, New York.

Wappler, T., De Meulemeester, T., Murat Aytekin, A., Michez, D. & Engel, M.S. (2012) Geometric morphometric analysis of a new Miocene bumble bee from the Randeck

Maar of southwestern Germany (Hymenoptera: Apidae). Systematic Entomology, 37, 784-792.

Widmer, A. & Schmid-Hempel, P. (1999) The population genetic structure of a large temperate pollinator species, Bombus pascuorum (Scopoli) (Hymenoptera: Apidae). Molecular Ecology, 8, 387-398.

Williams, P.H. (1985) A preliminary cladistic investigation of relationships among the bumble bees (Hymenoptera, Apidae). Systematic Entomology, 10, 239-255.

Yackulic, C.B., Chandler, R., Zipkin, E.F., Royle, J.A., Nichols, J.D., Campbell Grant, E.H. & Veran, S. (2013) Presence-only modelling using MAXENT: When can we trust the inferences? Methods in Ecology and Evolution, 4, 236-243.

Zagwijn, W.H. (1982) Migration of vegetation during the Quaternary in Europe. Courier Forschungsinstitut Senckenberg, 153, 9-20.

# 6. Acknowledgements

We thank Paul Williams and two anonymous referees for providing helpful comments and advices. We are also grateful to Jayna Raghwani for her comments and careful proof-reading of the final version of our manuscript. Computational resources were provided by the High Performance Computing Centre co-funded by ULB and VUB (HPC cluster 'Hydra') and by the Consortium des Équipements de Calcul Intensif (CÉCI, HPC cluster 'HMEM'), funded by the Fonds de la Recherche Scientifique (F.R.S.-FNRS, convention 2.5020.11). S.D., C.K. and T.L. were supported by a grant from the Belgian Fonds pour la Recherche dans l'Industrie et l'Agriculture (FRIA). S.D. was also supported

by the Wiener-Anspach Foundation and is now a post-doctoral research fellow funded by the Fonds Wetenschappelijk Onderzoek (FWO, Flanders, Belgium). T.L. was post-doc in the BELBEES project (BR/132/A1/BELBEES) and S.V. was funded by the project 'Support of establishment, development and mobility of quality research teams at the Charles University' (CZ.1.07/2.3.00/30.0022), of the Education for Competitiveness Operational Program (ECOP) of the Czech Republic and the European Science Foundation.