Unraveling the habitat preferences of two closely related bumble bee species in Eastern Europe

Abstract Co‐occurrence of closely related species is often explained through resource partitioning, where key morphological or life‐history traits evolve under strong divergent selection. In bumble bees (genus Bombus), differences in tongue lengths, nest sites, and several life‐history traits are the principal factors in resource partitioning. However, the buff‐tailed and white‐tailed bumble bee (Bombus terrestris and B. lucorum respectively) are very similar in morphology and life history, but their ranges nevertheless partly overlap, raising the question how they are ecologically divergent. What little is known about the environmental factors determining their distributions stems from studies in Central and Western Europe, but even less information is available about their distributions in Eastern Europe, where different subspecies occur. Here, we aimed to disentangle the broad habitat requirements and associated distributions of these species in Romania and Bulgaria. First, we genetically identified sampled individuals from many sites across the study area. We then not only computed species distributions based on presence‐only data, but also expanded on these models using relative abundance data. We found that B. terrestris is a more generalist species than previously thought, but that B. lucorum is restricted to forested areas with colder and wetter climates, which in our study area are primarily found at higher elevations. Both vegetation parameters such as annual mean Leaf Area Index and canopy height, as well as climatic conditions, were important in explaining their distributions. Although our models based on presence‐only data suggest a large overlap in their respective distributions, results on their relative abundance suggest that the two species replace one another across an environmental gradient correlated to elevation. The inclusion of abundance enhances our understanding of the distribution of these species, supporting the emerging recognition of the importance of abundance data in species distribution modeling.


| INTRODUC TI ON
The co-occurrence of closely related species has long puzzled evolutionary ecologists. Closely related species are expected to occupy similar niche space through niche conservatism and thus occur in the same regions, but also to compete more strongly with one another than more distantly related species (Anacker & Strauss, 2016). Strong competition may result in the exclusion of the weaker competitor, as well as in rapid divergent evolution of key life-history traits or phenotypes as a result of resource partitioning (Gause, 1934). In bumble bees (genus Bombus), a classical theory is that species evolved a range of different tongue lengths, allowing them to specialize on different floral resources, and to occur sympatrically as a result (Goulson & Darvill, 2004). Other mechanisms of resource partitioning include differences in nest sites, foraging distances, and the spatial use of habitat (Stanley, Knight, & Stout, 2013;Westphal, Steffan-Dewenter, & Tscharntke, 2006). However, two of the most common bumble bee species in Europe, the buff-tailed bumble bee (Bombus terrestris) and the white-tailed bumble bee (Bombus lucorum), co-occur (pers. obs.; Goulson, 2010;Goulson, Lye, & Darvill, 2008b;Kells, Holland, & Goulson, 2001;Stanley et al., 2013), despite being very similar in their morphology, choice of nest sites, and life-history (Goulson, 2010;Stanley et al., 2013). They even possess the same tongue lengths (Goulson, Hanley, Darvill, Ellis, & Knight, 2005;Stanley et al., 2013) and hence occupy a very similar dietary niche space (Goulson et al., 2008b). This begs the question to what extent their ecological niches overlap and-conversely-how they are divergent. Despite many studies into their ecology and behavior (e.g. Bossert, Gereben-Krenn, Neumayer, Schneller, & Krenn, 2016;Scriven et al., 2015;Stanley et al., 2013;Walther-Hellwig & Frankl, 2000), and broad-scale evidence that their ranges only partially overlap (e.g. Rasmont, Franzen, et al., 2015), their distributions at smaller scales remain equivocal. One reason for this ambiguity may be the fact that these species are morphologically highly variable within species, yet very similar between species (Bossert, 2014;Murray, Fitzpatrick, Brown, & Paxton, 2007;Waters, Darvill, Lye, & Goulson, 2011). As a consequence, they may be difficult to distinguish in the field, depending on where they occur and whether queens, males, or workers are compared.
In mainland Europe, identification can be complex because both species possess a white abdomen (Gammans, Comont, Morgan, & Perkins, 2018;Rasmont, Coppee, Michez, & Meulemeester, 2013), in contrast to those in Great Britain (Murray et al., 2007). Queens and males can be distinguished (Bertsch, Schweer, & Titze, 2004;Goulson, 2010), but workers (especially of some subspecies of B. terrestris) are difficult to discriminate (Williams, 1994). Indeed, in Central Europe, only 45.5% of B. lucorum workers could be correctly identified and distinguished from B. terrestris workers using the most up-to-date morphological key (Wolf, Rohde, & Moritz, 2009). As a consequence, many studies focusing on the ecology or behavior of European bumble bees group these taxonomically well-recognized species together (Bommarco, Lundin, Smith, & Rundlöf, 2011;Carvell, 2002;Goulson et al., 2005;Meek et al., 2002;Pywell et al., 2005;Walther-Hellwig & Frankl, 2000), leading to imprecise information on habitat preferences (Murray et al., 2007;Scriven et al., 2015) and other life-history traits (Stanley et al., 2013). One of the few studies that investigated the ecological preferences of those species separately identified differences in nesting site selection at small spatial scales in Sweden (Svensson, Lagerlöf, & Svensson, 2000) where B. terrestris preferred more open habitat, such as fields in agricultural landscapes, and B. lucorum preferentially built nests close to forest boundaries. Also in Austria, B. lucorum appeared to frequently occur in forested areas (Bossert et al., 2016). At the scale of Europe, B. lucorum occurs at higher latitudes than B. terrestris, suggesting a differentiation based on temperature (Rasmont, Franzen, et al., 2015). Most studies on these species focused on Western and Central Europe, but little attention has been paid to Eastern Europe, where the situation is complicated by the occurrence of two subspecies of B. terrestris, B. t. terrestris and B. t. dalmatinus, that are morphologically variable (Lecocq, Rasmont, Harpke, & Schweiger, 2016;Rasmont et al., 2013). Hence, their distribution patterns in Eastern Europe remain equivocal.
Here, we investigate the distributions and broad habitat characteristics of buff-tailed (B. terrestris) and white-tailed bumble bee (B. lucorum) in Bulgaria and Romania, where they are the two most common bumble bee species. B. lucorum is one of three cryptic species, the other two being B. cryptarum and B. magnus. Here, we only include B. lucorum, because we did not find any individuals of the latter two species, despite extensive sampling efforts. We hypothesized that the differential habitat use at small scales can be scaled up to landscape scale habitat preferences across a spatial extent of hundreds of kilometers. We first genetically identified the species at sites where multiple individuals were sampled, providing a reliable tool for species identification (Bossert et al., 2016;Murray et al., 2007;Stanley et al., 2013;Waters et al., 2011;Williams et al., 2012;Wolf et al., 2009). We subsequently created species distribution models (SDMs), which have been used previously in quantitative ecological studies of bumble bees (Casey, Rebelo, Rotheray, & Goulson, 2015;Herrera, Ploquin, Rodríguez-Pérez, & Obeso, 2014;Kadoya, Ishii, Kikuchi, Suda, & Washitani, 2009;Pradervand, Pellissier, Randin, & Guisan, 2014;Rasmont, Franzen, et al., 2015). They are usually based on presence-only or presence-absence data, with the assumption that the modeled probabilities of occurrence are indicative of abundance (e.g. Dallas & Hastings, 2018). However, recent work across multiple species suggests that these so-called abundance-suitability relationships are often weak (Dallas & Hastings, 2018;Howard et al., 2014;Mi, Huettmann, Sun, & Guo, 2017). For that reason, the collection and use of abundance data to improve the accuracy of species distribution models was highly recommended (Howard et al., 2014), but still not commonly applied to date. To improve our distribution models, and to specifically investigate their co-occurrence and the associated abiotic factors driving variation in abundance patterns, we therefore also collected and modeled the relative abundance of these species.

| Study species and study area
Bumble bees (Bombus sp.) are important pollinators for crops and wild plants, in particular in temperate ecosystems (Corbet, Williams, & Osborne, 1991;Murray et al., 2007). Their body is covered in a dense, colored fur that enables them to be endothermic (Heinrich, 1993), and hence adapt to cold climates, such as alpine and arctic environments. Their distribution extends much further north than that of other bees, and their colonies have been found in the extreme northern regions of the northern hemisphere (Goulson, 2010).
Bombus terrestris and Bombus lucorum are two of the most common bumble bee species in Europe. These species have very similar life cycles and are often found in the same areas. Both species use underground nest sites and often choose already existing holes, previously used by small rodents (Goulson, 2010). They possess similar tongue lengths, and as a result forage primarily on the same short corollas and daisy type of flowers (personal observation and (Goulson, Lye, & Darvill, 2008a). In addition, their workers almost perfectly resemble each other, and only the queens and males can be identified reliably (Murray et al., 2007;Wolf et al., 2009), but field identification remains complicated due to the subtlety of morphological differences.  (2015)). This high habitat heterogeneity represents an interesting area for evaluating habitat preferences and niche differentiation within and among species.

| Sampling
We collected 743 individuals compromising queens and workers of Bombus terrestris and Bombus lucorum over a timeframe of 4 years (2013, 2014, 2015, and 2017), in each between April and July. We sampled 44 locations in total ( Figure 1a and Table 1), which were selected to span a broad range of habitat conditions in both entirely natural and semi-natural or cultivated environments, as well as along environmental gradients (altitude, vegetation, and climate). We visited additional locations where we searched for, but did not find any bumble bees.
These locations were not included as true absences in our species distribution models, but served in computing a sampling bias map (see below). Sampling locations were located at least 20 km apart to rule out the possibility of overlapping foraging ranges (Chapman, Wang, & Bourke, 2003;Westphal et al., 2006) and were visited only once.
At each sampling location, capturing efforts were undertaken by 2-3 people for 1.5-2 hrs between 1 hr after sunrise and 1 hr before sunset. Individuals were collected on suitable forage patches with a radius of 100 m, using an entomological net. Individuals were visually identified as one of the two study species, anesthetized in a killing jar with a 1.5 cm layer of plaster of paris saturated with ethyl acetate, and immediately upon cessation of movement stored in 96% ethanol (Smithers, 1988). After fieldwork, specimens were stored frozen at −20°C in the laboratory at the University of Tübingen.

| Species identification
Because of the previously described difficulties in distinguishing species based on external morphology, we used a 1,043 bp long fragment of the cytochrome c oxidase subunit I (COI) gene, which is known for its relatively fast mutation rate, and is used across many taxa for genetic identification purposes, including bumble bees (Bertsch et al., 2004;Bossert et al., 2016;Hebert, Stoeckle, Zemlak, & Francis, 2004;Murray et al., 2007;Waugh, 2007;Williams et al., 2012). This long fragment is completely overlapping with an ~890 bp region used by Bertsch et al. (2004) to distinguish between the closely related Bombus cryptarum, B. magnus, and B. lucorum and was therefore used here to distinguish between B. lucorum and B. terrestris. DNA was extracted from one or both middle legs using the QIAmp DNA Micro Kit (Qiagen) following the manufacturer's protocol. PCR was performed with primers originally developed for Apis (Tanaka, Roubik, Kato, Liew, & Gunsalam, 2001): forward 5'-ATAATTTTTTTTATAGTTATA-3' and reverse PCRs were performed with a Mastercycler epgradient (Eppendorf) with the following conditions: an initial denaturation step at 95°C for 1 min, followed by 55 cycles of a 3-step process: denaturation for 40 s at 95°C, annealing for 1 min at 45°C, and extension for 2 min at 60°C with a final extension step at 60°C for 4 min. PCR products were visualized using agarose gel electrophoresis to check for the amplification of the fragment. Successfully amplified PCR products were cleaned up using the Promega Wizard SV Gel and PCR Clean-Up System according to the manufacturer's protocol. Cleaned up samples were then sent to LGC Genomics for sequencing with the reverse primer only. Sequences were visualized and edited with Geneious R8 (Biomatters, Ltd.; Kearse et al., 2012). None of the obtained sequences showed any signs of cross-contamination (e.g., double peaks in the chromatograms or ambivalent species identification). We used two methods to assign a species to the sequenced samples. First, we blasted sequences in GenBank (NIH genetic sequences database) and assigned the species with the highest identity (range ~95%-100%) to the corresponding sample (Table S1). In addition, we created a phylogenetic tree ( Figure S1), which included reference sequences for various Bombus species obtained from Genbank ( Table 2). We included reference sequences from various geographic origins, because we expected that genetic intraspecific variability should be smaller than interspecific differences, and thus that a well-supported clustering of our samples with the reference sequences suggests high confidence in the species identification. To construct the phylogenetic tree, we first identified the most likely substitution model in MEGA-X (Kumar, Stecher, Li, Knyaz, & Tamura, 2018). We then created a Bayesian phylogenetic tree in Geneious R8 with the MrBayes module (Huelsenbeck & Ronquist, 2001), using one cold and four heated chains with a temperature of 0.2 and a chain length of 1.1 million and a burn-in of 100,000, five gamma categories, a sampling frequency of 200, unconstrained branch lengths (GammaDir [1,0.1,1,1]), and exponential shape distribution. We used a reference sequence of B. pascuorum as an outgroup (Table 2). Individuals that were included in a monophyletic group with reference sequences of B. terrestris or B. lucorum were considered members of the corresponding target species. Comparisons to the blast results revealed no differences, and all individuals of B. terrestris and B. lucorum were included in subsequent species distribution models.

| Environmental variables
In order to create species distribution models, we used a set of environmental variables at 30 arcsec resolution (Table 3). We F I G U R E 1 Study region with sampling sites, species distribution modeling results based on ensembles of small models, and relative abundance modeling results. (a) Study area with sampling sites indicated in black stars for Bombus terrestris and in white diamonds for B. lucorum. Sites where we searched for bumble bees, but did not find any are indicated in black triangles. The background map is annual mean Leaf Area Index (LAI mean), a measure of vegetation greenness.

| Spatial autocorrelation and sampling bias
Spatial autocorrelation is a major statistical challenge in spatial analyses, causing inflated measures of predictive power and incorrect distribution models (e.g. Guélat, Kéry, & Isaac, 2018;Segurado, Araujo, & Kunin, 2006). There are two main causes for spatial autocorrelation in species distribution modeling.
First, there is often a spatial clustering of sampling sites. Reasons for such clustering are manifold and may be related to the sampling design (for instance ease of access, or issues with the logistics of evenly spaced sampling), or to the biotic and abiotic drivers of species distributions themselves, such as gaps in a species' range due to unsuitable habitat. Many approaches have been proposed to correct for sampling bias, among which subsampling locations to acquire a more even distribution of known presences are optimal in most cases (Fourcade, Engler, Rodder, & Secondi, 2014). As a first step, we therefore removed one of the sites of pairs that were located within 20 km from one another. However, because in our study the number of locations is rather limited, further subsampling would result in an even smaller data set. Hence, in a second step, we instead weighted each location based on the density of known presences within a given radius, which was shown to be a good alternative to subsampling as a correction method (Fourcade et al., 2014;Stolar, Nielsen, & Franklin, 2015). To do so, we created a bias grid in QGIS 3.4.4. (Development Team QGIS, 2017) at 30 arcsec resolution, with each grid cell representing the density of sampling locations within a 50 km radius, and kernel densities following a Gaussian distribution (Balestrieri et al., 2016). We used the inverse of the density to weight each presence and background location (see below), thus downweighting clustered locations. We not only included locations of known presence in this bias grid, but also locations where we searched for bumble bees with similar effort, but did not find any. We restricted these putative absence locations to those that were at least 50 km apart from known presences. We specifically only included these sites in our sampling bias map, and not in our models, because we cannot be sure that these represent true absences. The second cause of spatial autocorrelation in SDMs is the often inherent spatial autocorrelation of habitat conditions, in particular climate variables. In this case, species occurrences are spatially dependent on the underlying habitat variables and thus represent a true association between species presence and local conditions. It is often impossible and undesirable to a priori remove spatial autocorrelation due to spatial dependence. Instead, spatial autocorrelation is expected to be absent in model residuals, regardless of the presence of initial spatial dependence. Models should correctly predict the presence or absence of a species at any given location, independent of its spatial relation to other locations. We thus tested for spatial autocorrelation in the probabilities of occurrence at known presence locations using global Moran's I in the R package "lctools" v.0.2-7. We used four neighbors and tested the significance of correlations with resampling and randomization procedures.

| Presence-only data
To model species distributions based on presence-only data, we used an ensemble method, which has been shown to perform better than any given individual modeling method (e.g. Araújo & New, 2007;Elith & Graham, 2009;Marmion, Luoto, Heikkinen, & Thuiller, 2009 Ripley, 1996). In a recent study comparing 10 different modeling approaches implemented in ecospat and Biomod2, these were shown to be the top performing ones, while keeping computation times manageable (Breiner et al., 2018). We used model tuning to optimize the parameter settings for each model.
We generated input files using presence-only sites and 5,000 background points that were sampled randomly at a minimum distance of 20 km from known presences. To correct for sampling bias, we extracted weights for all locations, which were implemented using the Yweights argument in ecospat.
To evaluate model performance, we computed various evaluation scores and used K-fold cross-validation with subsets of training and testing data. The Boyce index (Hirzel, Lay, Helfer, Randin, & Guisan, 2006) is specifically designed and hence a particularly appropriate evaluation score for presence-only models. It is limited between −1 and 1, with 0 indicating model performance no better than random, and values close to 1 indicating high performance. We used the Boyce index to assess model performance, but also report the area under the receiver operator curve (AUC; Swets, 1988), Cohen's kappa (Cohen, 1960;Hirzel et al., 2006), and the true skill statistic (TSS; Allouche, Tsoar, & Kadmon, 2006). To create training and testing data partitions for K-fold cross-validation, we used spatial blocking. Partitioning the data into spatial blocks has the advantage over random allocation of sites that it is better suited to evaluate model performance in the potential presence of spatial autocorrelation (e.g. Roberts et al., 2017). If a model performs well, it is expected to correctly predict occurrences in both distant as well as nearby locations (Telford & Birks, 2009). We generated spatial blocks of training and testing data with the R package "blockCV" v.

| Relative abundance data
To test whether the relative abundance of B. lucorum compared to B. terrestris is associated with environmental conditions, we used a machine learning approach implemented in the "SuperLearner"  Breiman, 1996;Breiman, 2001;Breiman et al., 1984), and neural networks (ANN; Ripley, 1996). We also ran a very simple mean-of-abundance model as a baseline. We corrected for sampling bias using the weighting method described above, but we also fitted models to uncorrected abundance. We ran models on the full data set, where bagging and randomization were done internally, as well as on a partial data set, where we omitted 20% of the data, which were used as test data for independent cross-validation. For each model, we report its associated risk (a measure of model performance) and coefficient (the weight with which it is included in the ensemble). Response curves were created as described above.
To subsequently create a map of the predicted relative abundance of B. lucorum across the entire study area, we extracted the values for environmental variables at all 30 arcsec gridcells within Bulgaria and Romania. We then used the "predict.SuperLearner" function to project the known relationship between relative abundance and environmental conditions onto the entire landscape.
These values were imported and converted to a raster format in QGIS 3.4.4 (Development Team QGIS, 2017).

| Species identification
The most likely substitution model was the General Time Reversible (GTR) model with gamma distribution, which we implemented to create the phylogenetic tree. We found that 514 individuals clustered with reference sequences of B. terrestris and 220 with those of B. lucorum (Table 1). randomization z = 4.47, and randomization p < .001). We visually inspected the predictive map and compared it to maps of important environmental variables. We found that particularly high probability of occurrence was predicted for sites in Mediterranean Bulgaria, which is consistent with the pattern of seasonality in surface moisture (QSCAT season), the most important variable in predicting the species' distribution. We suspected that the remaining spatial autocorrelation was the result of spatial dependence rather than of sampling bias. We further tested for residual spatial autocorrelation in a second analysis, where we also extracted the predictions for sites where we searched for bumble bees, but did not find any, despite making the same sampling effort. These sites were the same as those used to generate a sampling bias grid and were located at least 50 km from known presences. Although these sites were not included in the models as true absences, we expected that a well-performing model should predict low probability of occurrence for these sites. Note: Five cross-validated models were run based on spatial blocks generated with the R package "blockCV". MAXENT.P is the MAXENT.Phillips model. ENS is the ensemble of small models. p = .096, randomization z = 1.65, and randomization p = .099), and we concluded that sampling bias was sufficiently well corrected for.

| Presence-only models
Interestingly, the most important variables in limiting each species' distribution overlapped between species. The top four variables for B. terrestris were seasonality in surface moisture (QSCAT season), mean temperature of the wettest quarter (Bio 8), mean leaf area index (LAI mean), and temperature seasonality (Bio 4; Table 5). For B. lucorum, these variables comprised mean leaf area index (LAI mean), canopy height, seasonality in surface moisture (QSCAT season), and percent tree cover (Tree ; Table 6), subsequently followed by mean temperature of the wettest quarter (Bio 8). For both species, the ranking of variables by their importance was largely consistent between modeling approaches. The main difference in the response curves between the two species is that those for B. lucorum are generally much steeper than those for B. terrestris, suggesting a stronger influence of the environment on B. lucorum (Figures 2 and 3). This difference is particularly pronounced for the top most important variables that were not overlapping between species, that is, percent tree cover and canopy height.

| Relative abundance models
Models of relative abundance that were corrected for sampling bias performed considerably worse than uncorrected models (corrected models: BART coefficient = 0.885, BART risk = 0.404, GLM coefficient = 0.115, and GLM risk = 0.682). We therefore report results for uncorrected models from here onwards. The only two models included in the ensemble were GAM (coefficient = 0.084, risk = 0.199) and BART (coefficient = 0.916, risk = 0.093). K-fold nested cross-validation with fivefold suggested that the single best model was BART, which performed even slightly better than the ensemble model, yet statistically nonsignificant (Table 7). Because of the high weight of the BART model, we evaluated variable importance based on BART only, providing a robust posterior importance score (Chipman et al., 2010;Hernández, Raftery, Pennington, & Parnell, 2018). The top most important variable was canopy height, subsequently followed by percent tree cover and three temperature variables (Figure 4), which is broadly consistent with the results for the presence-only data. Overall and marginal response curves suggest that B. lucorum is more abundant in more densely vegetated, wet and cool areas ( Figure 5).
To gain further insight in how our two target species differ in their preferred habitat conditions, we visually inspected scatter-   our analyses because of its high VIF, we also created a scatter plot of relative abundance versus elevation. Visual inspection of these plots suggested that B. lucorum does not occur in unforested areas with a canopy height under ~20 m (Figure 6a). Yet, the dichotomy between species is particularly striking for elevation, where B. lucorum is almost completely absent below 600 m, but makes up the majority of the two species at higher elevations ( Figure 6b). Hence, elevation captures the combined influence of correlated environmental variables on limiting the distribution of B. lucorum particularly well.

| D ISCUSS I ON
Here, we modeled the distributions and relative abundance of two cryptic bumble bee species in Bulgaria and Romania from samples that were assigned to one of the species using a long fragment of the COI gene for genetic identification. We demonstrated that even though there is a certain degree of overlap between the ranges of the two species, B. lucorum has a much more restricted  (Bossert et al., 2016;Svensson et al., 2000). Our findings also support observations that B. lucorum occurs at higher elevations (Ban-Calefariu & Sárospataki, 2007;Bossert et al., 2016;Goulson et al., 2008b;Ploquin, Herrera, & Obeso, 2013;Tomozei, 2006), which was suggested to be the result of an adaptation to colder climates (Benton, 2006). Indeed, in northern Europe, B. lucorum generally occurs in colder areas, where it at least partly substitutes B. terrestris (Rasmont, Franzen, et al., 2015). Many environmental variables change along an elevation gradient, and elevation itself is unlikely to determine the distribution of these species, but rather its covariates (Bossert et al., 2016). In our study, mean temperature of the coldest quarter (Bio 11), canopy height, percent tree cover, and mean leaf area index (LAI mean) were particularly highly correlated with elevation (Pearson correlation coefficient > 0.6).
Hence, the distribution of B. lucorum is clearly restricted to the mountainous areas in Bulgaria and Romania (Figure 1c,d), where temperatures are lower, precipitation is higher, and where most of the forest is remaining.
In contrast, our findings for B. terrestris suggest that it is not as restricted to open habitat as previously thought (Bossert et al., 2016;Svensson et al., 2000), but rather is a generalist species, occurring in open as well as more densely vegetated areas. This notion is also apparent in our maps, showing a wide distribution for B. terrestris.
Interestingly, our presence-only species distribution models showed a considerable overlap between the ranges of B. terrestris and B. lucorum, but analyses of their relative abundance evoke a much stronger separation between these species. Although the use of relative abundance does not allow for conclusions regarding the absolute abundance of either one of the species, the large range of relative abundance values for B. lucorum, spanning from 0 to 0.94, suggests that the two species replace one another across an environmental gradient. Thus, the inclusion of abundance enhances our understanding of the distribution of these species based on presence-only models.
Although the conservation status of our study species across the European continent is "least concern" (Rasmont, Roberts,  (Rasmont, Franzen, et al., 2015). The study by Rasmont, Franzen, et al. (2015) provides a great overview of overall Indeed, we found that vegetation characteristics were among the most important variables explaining the distribution and relative abundance of our study species, and it will be difficult to predict how these variables will change in the future, both as a result of climate change, as well as due to changes in forest management.
We did not proceed with an attempt to predict the distribution of B. lucorum onto future climate conditions, because a model based on only current climate conditions failed to even broadly resemble that based on both climate and vegetation variables (not shown).
Moreover, populations are likely adapted to local and regional conditions, and may not respond the same to changing environmental conditions. Our study provides further insight by teasing apart the habitat preferences of these species in southeastern Europe, providing higher resolution range maps that are probably more relevant for the region, where the distribution of B. lucorum is assumed to be rather disjunct. Despite the complexity of predicting future changes in vegetation characteristics, the difference in habitat requirements between these species is expected to have implications for the way they respond to changing climate conditions. Our finding that B. lucorum is rather restricted in its suitable habitat conditions compared to B. terrestris, may suggest that it is more vulnerable to climate change than the latter.
We genetically identified a large number of individuals of two closely related bumble bee species sampled at many sites and modeled their distributions and gained insight into their habitat requirements. We showed that B. terrestris is rather a generalist species, whereas B. lucorum is restricted to colder and wetter climates in forested areas, which in our study area primarily occur at higher elevations. We support the emerging recognition of the importance of abundance data in species distribution modeling. Despite the overlap in occurrence suggested by presence-only data, their relative abundance gradually changes along a major environmental gradient, with one of the species being virtually absent at the extreme ends of this gradient. Our study contributes to the urgent need to fill a major gap of knowledge in the distribution and ecology of these species that can help facilitate the assessment of their conservation status and the development of management plans where necessary.

ACK N OWLED G M ENTS
We are particularly grateful to Frank Breiner, who advised on ESM modeling and implemented new functionality in ecospat in response to our questions. We thank Isabel Fraccaroli for help in the molecular laboratory, and two anonymous reviewers, whose contribution much improved a previous version of this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences are available from Genbank, accession no.