Conserv Genet (2013) 14:1099–1110 DOI 10.1007/s10592-013-0498-3

RESEARCH ARTICLE

Landscape heterogeneity predicts gene flow in a widespread polymorphic bumble bee, Bombus bifarius (Hymenoptera: Apidae) Jeffrey D. Lozier • James P. Strange Jonathan B. Koch

•

Received: 8 March 2013 / Accepted: 18 May 2013 / Published online: 26 May 2013 Ó Springer Science+Business Media Dordrecht 2013

Abstract Bombus bifarius is a widespread bumble bee that occurs in montane regions of western North America. This species has several major color pattern polymorphisms and shows evidence of genetic structuring among regional populations, and the taxonomic status of regional populations has repeatedly been debated. We test whether observed structure is evidence for discrete gene flow barriers that might indicate isolation or instead reflects clinal variation associated with spatially limited dispersal in a complex landscape. We first consider color pattern variation and identify geographical patterns of B. bifarius color variation using cluster analysis. We then use climate data and a comprehensive set of B. bifarius natural history records with an existing genetic data set to model the distribution of environmentally suitable habitat in western North America and predict pathways of potential gene flow using circuit theory. Resistance distances among populations that incorporate environmental suitability information predict patterns of genetic structure much better than geographic distance or Bayesian clustering alone. Results Electronic supplementary material The online version of this article (doi:10.1007/s10592-013-0498-3) contains supplementary material, which is available to authorized users. J. D. Lozier (&) Department of Biological Sciences, University of Alabama, Tuscaloosa, AL 35407, USA e-mail: [email protected] J. P. Strange  J. B. Koch United States Department of Agriculture-Agricultural Research Service, Pollinating Insects Research Unit, Utah State University, Logan, UT 84322, USA J. B. Koch Department of Biology, Utah State University, Logan, UT 84321, USA

suggest that there may not be barriers to gene flow warranting further taxonomic considerations, but rather that the arrangement of suitable habitat at broad scales limits dispersal sufficiently to explain observed levels of population differentiation in B. bifarius. Keywords Landscape genetics  Isolation by distance  Isolation by resistance  Environmental niche model  Microsatellites  Circuit theory  Color polymorphism

Introduction Understanding the processes governing genetic diversity and differentiation of populations is a central focus of population genetics. Gene flow is the major force that maintains connectivity among populations, and many studies seek to identify relationships between space and genetic variation with the aim of identifying barriers to gene flow (Manel et al. 2003; Storfer et al. 2007; Guillot et al. 2009). Classic approaches involve testing for correlations between pairwise genetic (e.g. FST or other measures) and geographic distance separating individuals or populations assuming models of isolation by distance (IBD) (Slatkin 1993; Rousset 1997). A significant IBD signature, for example, suggests that a species is at genetic equilibrium with spatially restricted dispersal, while lack thereof can suggest that long-distance dispersal is sufficiently common to overcome IBD signatures or that populations are not at equilibrium (Slatkin 1993; Hutchison and Templeton 1999; Guillot et al. 2009). Model-based approaches such as Structure or Geneland (e.g., Falush et al. 2003; Guillot 2005) can also be used to assign genotypes to potentially unknown source populations, and may reveal more discrete patterns of population structure

123

1100

that indicate barriers to gene flow, as might be found for recently isolated species. For organisms that occur in regions with substantial environmental heterogeneity, populations will not be continuously distributed throughout the landscape, but will occur primarily in areas with suitable ecological conditions. Likewise, dispersal will not be random with respect to two-dimensional space, but will occur with higher probability along paths of high habitat suitability. Identifying how such heterogeneity affects inferences of gene flow from genetic data is one of the primary goals of landscape genetics (Manel et al. 2003; Storfer et al. 2007; Storfer et al. 2010). Isolation by resistance (IBR) modeling incorporates spatial data to create a ‘‘resistance’’ surface representing the probability with which a species can disperse across a landscape. IBR models often outperform simple IBD analyses because they can better predict actual pathways of dispersal among populations (McRae 2006; McRae and Beier 2007) and identify environmental factors shaping these paths (Cushman et al. 2006). Accurate inference of connectivity has clear implications for conservation, providing a foundation for predicting future gene flow as the spatial arrangement of suitable habitat is altered by climate and land use changes (Cushman et al. 2006; Wilson et al. 2005; McRae and Beier 2007; Storfer et al. 2010). Similarly, because IBR can be used to distinguish the signatures of apparent isolation from one of equilibrium gene flow, such methods can improve detection of more discrete population structure that might represent subspecies or other evolutionarily significant units, versus patterns that can largely be attributed to range dimensions (McRae and Beier 2007). Bombus bifarius is a widespread North American (NAm) bumble bee that ranges from Alaska to New Mexico and from sea level to [3,000 m in elevation. B. bifarius is among the more phenotypically variable species in NAm, with abdominal pile color ranging from bright red-banded forms in the southern Rocky Mountains to black-banded forms in central populations, and red-banded forms reemerging in the northwestern portion of the species range (Fig. 1). Such variability led Stephen (1957, p. 139) to note that B. bifarius ‘‘has undoubtedly caused more consternation to bumble bee taxonomists than any other species in western NAm,’’ and various taxonomic epithets have been ascribed to the different B. bifarius color morphs (Stephen 1957) (Fig. 1). Color pattern divergence is a potential signal pointing to cryptic specific or subspecific barriers to gene flow in insects (e.g., Jiggins et al. 2001;), and this appears to be the case in some bumble bees (Rasmont et al. 2008; Duennes et al. 2012; Williams et al. 2012a) but not necessarily in others (Carolan et al. 2012). Besides this striking color variation, B. bifarius is also distinct among the NAm bumble bees in

123

Conserv Genet (2013) 14:1099–1110

Fig. 1 Map of B. bifarius data from the western U.S. used for genetic analyses. Gray scale shading indicates the predicted Maxent distribution (logistic suitability) used in the Circuitscape analysis. Symbols indicate membership of sampled populations to Geneland clusters (K = 4; see Supplementary Fig. 2 for probability surfaces). Intercluster FST values estimated by Geneland are provided in the inset. For analyses with K = 3, Geneland clusters the weakly differentiated diamond and triangle populations into a single cluster. The major color pattern groups with varietal names used in the text are also shown

terms of population genetic structure. Most NAm species examined to date show little evidence for regional gene flow barriers, whereas B. bifarius exhibits evidence for more regional clustering into southeastern, central, and northwestern populations (Lozier et al. 2011). Genetic clustering is broadly consistent with the distribution of major color patterns (Stephen 1957), although groupings are not perfect. Overall, the spatial trends in color and genetic structure raise questions regarding the possibility of historical isolation of regional populations, warranting further study. Montane species such as B. bifarius are likely to be particularly affected by range dimensions (Cushman et al. 2006; Schwartz et al. 2009; Wasserman et al. 2010). Because the habitat mosaic associated with elevation changes may limit dispersal pathways, discrepancies between distance and resistance may be especially strong, altering our perceptions of population connectivity. Environmental niche models (ENMs) offer a promising tool to quantify these complex species distributions (Elith and Leathwick 2009). By modeling the relationship between species observations and environmental characteristics associated with those localities, ENMs generate a spatial prediction of environmental suitability that can be applied to estimate gene flow routes or barriers (Rissler and Apodaca 2007; McRae et al. 2008; Abella´n et al. 2011;

Conserv Genet (2013) 14:1099–1110

Chan et al. 2011). ENM predictions can also identify the environmental variables contributing most to a species’ distribution (Elith et al. 2010; Warren 2012) and, when integrated with population genetic data, might have the greatest impacts on gene flow. A comprehensive understanding of B. bifarius population structure serves several purposes. First, this species is a potential candidate for development as a commercial pollinator in western NAm (Strange 2010), and the nature of genetic structure in B. bifarius, including the need for delimitation of species or other evolutionarily significant units, will provide data for assessing risks associated with domestic transport for pollination services (Goulson 2010; Duennes et al. 2012; Williams et al. 2012b). Second, other montane western NAm bumble bees have suffered serious declines in recent decades (Cameron et al. 2011). Earlier analyses for B. bifarius show that populations at higher elevations tend to show greater genetic isolation (Lozier et al. 2011), and although this species is currently common, dispersal could be challenged if pathways of suitable habitat are altered in the future (Wilson et al. 2005). By understanding how environmental variables, including climate, shapes dispersal routes in a species that remains abundant, we may gain insights into factors that could contribute to such declines in other species. To more explicitly model the forces shaping gene flow in B. bifarius, we apply novel landscape genetic approaches to previously generated microsatellite data that focused on describing general population genetic patterns for B. bifarius and other NAm Bombus (Lozier et al. 2011). First we utilize natural history collection specimens to test whether previously suggested subspecific B. bifarius groups can be recovered with a quantitative analysis of color pattern variation. We then apply an existing genetic data set toward testing the hypothesis that environmental heterogeneity, rather than pure IBD or discrete barriers to gene flow, shapes patterns of genetic differentiation among populations of this polymorphic species. Our analysis provides the first attempt to link ENMs with population genetics in bumble bees, and offers a straightforward approach for studying gene flow and testing species-level taxonomic hypotheses in other montane Bombus.

1101

(Supplementary Table 1). Pinned specimens reflecting the geographic range of our genetic analysis were selected from a recent study (Cameron et al. 2011). We scored ten body regions associated with setae color variability in bumble bees (Supplementary Fig. 1; from Williams 2007) and are similar to those historically used to qualitatively diagnose B. bifarius subspecies (Stephen 1957; Thorp et al. 1983). Setae color values were captured with a digital photograph of the dorsal region of each specimen using a Keyence VHx-500f camera-mounted dissecting scope. Magnification was kept constant and photographs of each specimen were taken within 2 h of each other in the late evening to reduce variability from ambient sunlight. Setae colors were sampled from a single pixel in each body region using the RGB color model (red, green, blue), which were later transformed to a HSL color model (hue, saturation, lightness), which has proven valuable to studies of insect color evolution and mating systems (Davis et al. 2007). We excluded a small number of poor-quality specimens with matted setae that produced extreme outlier color values. To determine our modeling approach, we visualized the distribution of the HSL color model data for each body region with histograms and box plots with R 2.15.1 (R Core Team 2012). Therefore, each body region was associated with three measurements of color properties for 30 data points per specimen. H was normally distributed, but S and L were positively skewed. Non-normality of the data was confirmed using a Generalized Shapiro–Wilk test (MVW = 0.9652, P = 0.0353) and we used a non-parametric K-means cluster analysis to identify distinct color phenotypes in two-dimensional space. To determine the optimal K for our analysis we calculated the within group sum of squares for all values between one and the maximum number of geographic regions (eight or eleven). We also conducted a Monte Carlo Cluster Analysis to determine the optimal K for the final K-means cluster analysis. HSL values were transformed from RGB values using the color library in Python 2.7.3 (http://www.python.org). All color pattern statistical analyses were executed in R with the car, psych, cluster, mclust, and adegenet libraries. Genetic data

Methods Analysis of color pattern To quantify pile color, we examined setae for 10 individuals from a set of western NAm populations representing regions included in genetic analyses (Fig. 1) as well as from a number of geographically intermediate populations from which we did not have sufficient genetic data

We take our molecular data from an earlier study of population structure in NAm Bombus (Lozier et al. 2011). Briefly, these data consist of coarsely sampled populations from throughout the B. bifarius range in the U.S.A., including sites in Alaska, California, Colorado, Idaho, Montana, Nevada, Oregon, Utah, Washington, and Wyoming. Specimens were genotyped for nine microsatellite markers (doi:10.5061/dryad.d403s) that showed no

123

1102

deviations from Hardy–Weinberg equilibrium, and fullsibs were excluded. Here we exclude sites with fewer than 10 sampled individuals to minimize noise from sampling error. We also excluded a single Alaska site because of the large distance from all other samples. The final data set consisted of 26 populations and 447 individuals (Fig. 1). Environmental niche model (ENM) To quantify landscape heterogeneity we generated an ENM using the principle of maximum entropy (Elith et al. 2010). We first obtained spatially referenced natural history collection records for B. bifarius as described in Cameron et al. (2011). After exclusion of duplicate records and specimens falling outside of our study extent (approximately 50.0°N, 30°N, 125°W, 104°W), we retained 1,537 data points for analysis (Fig. 2). Spatially explicit environmental variables for contemporary conditions were taken from the WorldClim V1.4 (Hijmans et al. 2005) BIOCLIM data set (2.5 arc-minute, or *5 km, resolution) and clipped to our study area using ArcMap 10.0 (ESRI, Redmond CA). We used Maxent v3.3.3 (Phillips et al. 2006) to generate the ENM. Maxent uses presence-only locality data and random background points sampled from a study area to estimate the species distribution that is closest to uniform, subject to limited information on the true distribution and environmental conditions (Elith et al. 2010). We elected to run models using a set of eight variables that captured annual and seasonal trends in temperature and precipitation likely to be relevant to bumble bees, which are only active during part of the year, but

Fig. 2 Map of 1,537 western North American B. bifarius natural history collection records (see Cameron et al. 2011) used for environmental niche modeling in Maxent

123

Conserv Genet (2013) 14:1099–1110

must also survive winter hibernation. Thus, montane bumble bees require suitable year round temperatures, but also sufficient yearly precipitation to provide both snowpack and rain for floral resources (Aldridge et al. 2011) (Table 1). We ran a final Maxent analysis with default parameters to generate a logistic output, which we refer to as a relative measure of environmental suitability, averaged over 10 cross-validated replicates (see Elith et al. 2010 for detailed discussion of the logistic output). We performed jackknife analysis to examine individual contributions of each environmental variable to the model. Estimating resistance and distance among populations We developed estimates of population connectivity by combining Maxent models of habitat suitability with the circuit theory approach implemented in Circuitscape V3.5 (Shah and McRae 2008). Circuitscape uses electrical circuit theory to generate IBR models of gene flow. The resistance distance between two points reflects the likelihood of potential gene flow, integrating over all possible alternative pathways of dispersal. When spatial heterogeneity in habitat quality is an important predictor of species dispersal, this approach should improve the relationship with genetic distance over IBD analyses (McRae and Beier 2007). We used the logistic output raster from Maxent as a conductance (the inverse of resistance) layer for Circuitscape, where areas of high predicted suitability values specify greater conductance, and thus the potential for greater gene flow; water was specified as ‘‘no data.’’ We ran Circuitscape in pairwise mode with the four-neighbor cell connection scheme calculated using average resistances. Because of the resolution of environmental data used for distribution models, the two coastal islands (San Juan Island and Orcas Island, WA) were connected to the mainland by 1–2 cell-wide corridors, but otherwise surrounded by completely unsuitable habitat that excluded gene flow. Given our interest in improving inferences about coarse regional-scale population structure, we considered such narrow corridors to be a sufficient representation of restricted over-water dispersal, rather than attempting to assign specific resistance values to water. As seen in results (below) this seems adequate at this regional scale. To compare our IBR model to more standard IBD approaches, we estimated resistance distances using a flat two-dimensional landscape with the same dimensions as the Maxent logistic output raster but where all cells were given conductance values of 1.0. Results using this flat landscape scenario closely mirrored patterns from analysis of ln-transformed great circle distances (not shown), and because it shares spatial properties with the raster-based resistance estimates, we focus our results on the former for more straightforward comparison of IBD and IBR.

Conserv Genet (2013) 14:1099–1110

1103

Table 1 Analysis of environmental variable importance for the Maxent environmental niche model used to estimate heterogeneity of habitat suitability for B. bifarius Jackknife analysis Percent contribution Annual precipitation

Training gain withoutb

Training gain with onlya

Test gain with onlyc

Test gain withoutd

52.11

0.39

0.58

0.40

0.61

Max temperature of warmest month

7.03

0.30

0.61

0.31

0.64

Mean temperature of wettest quarter

19.76

0.28

0.60

0.29

0.64

Precipitation of driest month

5.13

0.25

0.58

0.26

0.61

Annual mean temperature

4.13

0.22

0.62

0.23

0.65

Precipitation of wettest month

3.67

0.19

0.59

0.21

0.61

Mean temperature of driest quarter

2.19

0.11

0.61

0.12

0.64

Min temperature of coldest month

5.98

0.06

0.61

0.07

0.64

a

Training gain achieved with each climate variable in isolation

b

Training gain achieved excluding each variable in turn Test gain achieved with each climate variable in isolation

c d

Test gain achieved excluding each variable in turn

Genetic structure To estimate potential locations of discrete population structure in the data we applied the Bayesian clustering algorithm employed in Geneland 4.0.0 (Guillot 2005), implemented in R, which explicitly takes into account the spatial location of sampling sites and estimates the optimal number of population clusters. Site coordinates were provided for each individual, with a random adjustment of [-0.001 - 0.001] decimal degrees. We ran five replicates for two million iterations (thinning = 200) with the number of possible clusters K ranging from 2 to 5 and checked for consistency. We processed a final run on a landscape of 100 9 100 cells and with a burn-in of 2,000 iterations, fixing the maximum K to four, the value with the highest posterior density from all preliminary runs (e.g., Leache´ 2011). Given the geographic patterns evident in both color pattern and Geneland clusters (see below), we treat these clusters as a proxy to test the alternate hypothesis of barriers to gene flow between B. bifarius subspecies, without proposing an explicit mechanism for these barriers. Pairwise FST was estimated using Genepop 4.1.2 (Rousset 2008) and transformed to FST/(1-FST). Note that FST is highly correlated with other measures of genetic distance/ differentiation (e.g. Jost 2008) in B. bifarius. We used the R package ecodist (Goslee and Urban 2007) to perform Mantel tests and partial Mantel tests comparing the relationship between each of the landscape distances (geographic = IBD and resistance = IBR) and the FST/(1-FST) matrix. To examine the potential effects of discrete population structure, we also performed Mantel tests coding population pairs as derived from the same (0) or different (1) genetic clusters as identified with Geneland and applying partial Mantel tests

to partial out effects of IBD or IBR. Although Mantel tests have been criticized (Raufaste and Rousset 2001), several studies have demonstrated the continued utility of this approach for distinguishing IBD and IBR phenomena (Cushman and Landguth 2010; Legendre and Fortin 2010; Landguth et al. 2011). Models were compared by visual examination of distance scatter plots, Mantel r-values, and Bonferroni adjusted P-values for tests of positive correlation between matrices, estimated with 10,000 permutations of the response matrix in ecodist (Goslee and Urban 2007; Legendre and Fortin 2010). Because of concerns that island populations could be driving results (based on visual results of scatterplots), we performed Mantel tests for IBD and IBR that excluded island populations. We also tested whether resistance-based estimates of population connectivity better explained patterns of genetic structure than distance-based estimates using the hierarchal Bayesian F-model (Foll and Gaggiotti 2006; Gaggiotti and Foll 2010), implemented in GESTE 2.0 (Foll and Gaggiotti 2006). GESTE estimates population-specific FST values that indicate the degree of drift relative to the metapopulation as a whole, allowing for local differences in population size and migration rates (Gaggiotti and Foll 2010). We tested the effects of population isolation on FST using GESTE’s generalized linear modeling approach. Following Foll and Gaggiotti (2006) and Balkenhol et al. (2009), we estimated connectivity for each population j, Sj, according to j X expðbdij Þ Sj ¼ i¼1 i6¼j

where dij represents the distance/resistance between a population pair, and ß indicates the effect of distance on

123

1104

migration probability. We specified ß = 1, although preliminary analyses showed that, qualitatively, results were insensitive to the choice of ß. Sj estimates were normalized, GESTE was run with default parameters using multiple replicates to ensure consistency, and models were evaluated based on relative posterior probabilities.

Results Color pattern We performed an initial color pattern analysis focusing on populations from regions available for genetic analyses. The optimal color pattern cluster size was estimated at K = 5 based on the within group sum of squares calculation. However, three cluster centroids substantially overlapped upon visualization of the data, all of which reflected subtle variations of the black-banded form of B. bifarius, perhaps unsurprising given the broad geographic range of this form (see also Population Genetics results below). We thus simplified the model by rerunning the cluster analysis for K = 3, which collapsed the overlapping clusters (Fig. 3a). This reduced cluster comprises black-banded individuals (Fig. 1, 3a) found in the western and northern part of the B. bifarius range. The second cluster comprises individuals from the southeastern parts of the B. bifarius range, reflecting the red-banded from. The third cluster is specifically composed of individuals from the San Juan Islands, Washington population (Fig. 1, 3a). We found no evidence for misclassification of phenotypes among regional populations in this analysis (Fig. 3a). Overall, these findings are consistent with previous designations applied to B. bifarius: B. bifarius nearcticus, B. bifarius bifarius and B. vancouverensis (here, considered a ‘variety’), respectively, (Viereck et al. 1904; Stephen 1957), suggesting this analysis can recover differences among major morphological groups. The optimal cluster size for the second analysis including populations from geographically intermediate regions was estimated at K = 6 (Fig. 3b, Supplementary Table 1). Similar to the first analysis (Fig. 3a), three clusters representing black-banded B. bifarius nearcticus individuals overlap substantially. Adjacent to the blackbanded clusters, three more clusters overlapped, representing red-banded B. bifarius bifarius, and also included the San Juan Islands population (‘vancouverensis’ variety), both of which have red abdominal setae. Even so, pigmentation of the San Juan Island population remains somewhat distinct (Fig. 3b). Overall, there is a greater degree of intermediate color variation than observed in Fig. 3a when intermediate populations were included, although the three typical forms are still the most extreme

123

Conserv Genet (2013) 14:1099–1110

phenotypes, and the clouds representing red ‘bifarius’ and black ‘nearcticus’ forms remain largely distinguishable. Ecological niche model Average AUC for the cross validated B. bifarius distribution model was 0.820 (s.d. = 0.002) for training and 0.810 (s.d. = 0.017) for test points, and the overall pattern of high-suitability regions reflects previous coarse-grained distribution maps for this species (e.g. Thorp et al. 1983). The Maxent analysis revealed annual precipitation to be a strong limiting factor for B. bifarius, with this variable having the strongest influence on the final model as indicated by percent contribution (52.1 %), highest gain when used in isolation, and greatest decrease in gain when excluded (Table 1). Annual precipitation also overwhelmingly provided the greatest contribution to preliminary models run using all 19 BIOCLIM variables. Overall, the model predicts suitable habitat to be strongly associated with forested mountain ranges in the western US, with drier intervening basins and deserts constituting relatively low suitability habitat through which dispersal might be limited. Population genetics The clustering approach implemented in Geneland consistently detected four genetic clusters (Supplementary Fig. 2), versus the three most obvious clusters that were observed in earlier Structure analyses (Lozier et al. 2011). Results suggest an optimal partitioning of populations in southeastern, northern, southwestern, and northwestern populations (Fig. 1). If clustering is restricted to K = 3 in Geneland, however, results are comparable to Structure (Lozier et al. 2011) and the geographic distribution of three major color forms, with southeastern and northwestern sites separated from more central sites (diamond ? triangle populations; Fig. 1). The new Geneland cluster results from the split of populations from California and Oregon, similar to the subclusters observed within the ‘nearcticus’ form with the color data. However, FST indicates that this cluster is very weakly differentiated from other central populations (Fig. 1), and we evaluate both the K = 3 and K = 4 models. Just as color pattern variation became more complex with inclusion of certain geographically intermediate populations, landscape genetic analyses revealed signatures consistent with gene flow among regions. For all analyses there was a significant increase in genetic distance with increasing distance between populations, however, there was a large amount of scatter associated with IBD plots relative to the IBR plot that could be in part attributed to within-and between-Geneland cluster population pairs (Fig. 4; Lozier et al. 2011). There was a notable improvement in Mantel test

Conserv Genet (2013) 14:1099–1110 Fig. 3 a Principal components analysis (K = 3) of HSL color pattern data taken from geographic regions used in our genetic analysis and representing the typical extreme color forms for B. bifarius (N = 76; see Supplementary Table 1). The first two axes explain 33.92 and 12.28 % of the variation, respectively. b Principal components analysis (K = 6) including intermediate geographic populations (N = 103). Axes 1 and 2 explain 30.14 and 10.64 % of the variation, respectively. Clusters reflect 95 % data concentration ellipses determined with the scatterplot function of the car library in R

1105 6

A

“vancouverensis” forms

San Juan Islands

N. Idaho

S. Oregon

N. Wyoming

N. Washington

S. Wyoming Colorado

E. Oregon

4

“nearcticus” forms

2

0 -7

-5

-3

-1

1

3

5

7

-2

“bifarius” forms

-4

-6 San Juan Islands S. Oregon

6

B

N. Washington E. Oregon 4

N. Idaho N. Wyoming S. Wyoming Colorado Utah New Mexico Nevada

2

0 -8

-6

-4

-2

0

2

4

6

8

-2

-4

-6

r-values when the ENM-based distribution of B. bifarius was taken into account in the IBR analysis (r = 0.835, P \ 0.001, 95 % CI: 0.810–0.853) versus the IBD analysis (r = 0.473, P \ 0.001; 95 % CI: 0.417–0.519). The IBR model had the highest r-value observed (Table 2), with a strong linear relationship obtained between genetic and resistance distance (Fig. 4b). We did not explicitly model dispersal over water by assigning specific resistance values, but the resulting plots showed that incorporating information about restricted dispersal associated with islands via narrow habitat corridors made substantial improvements to the relationship between genetic and spatial distances. We thus

wanted to test whether the island populations alone were responsible for the larger r-values in the IBR analysis. Significant differences in Mantel r-values were also observed when island populations were excluded (r = 0.550, 95 % CI: 0.499-0.597 for IBR vs. r = 0.342, 95 % CI: 0.292-0.409 for IBD), demonstrating that the improvement from incorporating ENM-based resistances also applies to environmentally heterogeneous mainland landscapes. We further tested the role of connectivity versus distance among mainland populations using GESTE’s F-model assessment of population structure and generalized linear modeling framework. Results from Mantel tests are supported by the

123

1106

Conserv Genet (2013) 14:1099–1110

FST /(1-FST )

Table 2 Results from Mantel and partial Mantel tests of effects of distance, resistance, or discrete population structure on genetic differentiation [FST/(1-FST)] among 26 B. bifarius populations

Geographic Distance

Modela

Mantel r-valueb

FST/(1–FST) * resistance

0.835b

FST/(1–FST) * distance

0.473b

FST/(1–FST) * resistance ? distance FST/(1–FST) * clusterK4

0.845b 0.527b

FST/(1–FST) * clusterK4 ? distance

0.385b

FST/(1–FST) * clusterK4 ? resistance

0.156c

FST/(1–FST) * clusterK3

0.681b

FST/(1–FST) * clusterK3 ? distance

0.581b

FST/(1–FST) * clusterK3 ? resistance

0.225c

a For partial Mantel tests, the variable being partialed out follows the ‘?’ symbol b Bonferroni corrected P-values \ 0.005

FST /(1-FST )

c

N.S. not significant

Table 3 GESTE F-model analysis of resistance vs. distance for mainland populations, demonstrating the impact of landscape resistance on genetic structure over distance alone Model

Resistance Distance Fig. 4 Isolation by distance (a) and isolation by resistance (b) plots. In (a) geographic distance was calculated in Circuitscape using a flat conductance raster (Mantel test r = 0.473) and in (b) pairwise resistance (Mantel test r = 0.835) was calculated using the logistic Maxent output (Figure 1) as a conductance raster. In both plots, black points represent population pairs from the same Geneland cluster, and open points indicate pairs from different clusters (K = 4)

F-model analysis of FST and population connectivity (Table 3). The model including a constant and the resistance-based connectivity measure alone was more strongly supported (posterior probability = 0.753) than those including distance-based connectivity (posterior probability = 0.136) or both measures of connectivity (posterior probability = 0.075). As expected under IBR, populations isolated by greater levels of resistance had higher FST values (a2 = -0.642), with connectivity providing a reasonably good explanation for the observed differentiation (r2 = 0.384) (Foll and Gaggiotti 2006). We also tested whether regional genetic clustering may reflect isolation of populations not explained by spatial distance or landscape resistance alone. We used partial Mantel tests to examine the relationship between genetic distance and Geneland cluster assignment for each population pair (same cluster = 0 vs. different cluster = 1) and partialling out resistance or spatial distances. We predict that if discrete gene

123

Posterior probability

Constant

0.035

Constant ? distance-based connectivity Constant ? resistance-based connectivity

0.136 0.753

Constant ? distance-based ? resistance-based connectivity

0.075

Parameter (factor) of model 3

Posterior mode [95 % HPDI]

a1 (constant) a2 (resistance-based connectivity)

-4.26 [-4.70; -3.89] -0.642 [-1.08; -0.253]

r2

0.384 [0.150; 0.962]

a represents regression coefficient for factors, r2 measures deviation from the regression (Foll and Gaggiotti 2006)

flow barriers exist among clusters, then partial Mantel test r-values should remain significantly positive, whereas if distance/resistance accounts for most of the genetic structure, there will be no significant correlation. A simple Mantel test showed a highly significant effect of cluster membership alone for both K = 3 and 4 (r = 0.527 and 0.681, respectively). When IBD was taken into account, a significant relationship between genetic distance and clustering remained, albeit with reduced r-values (r = 0.385 and 0.581; P \ 0.005). However, partialling out resistance distances substantially reduced the clustering effect (r = 0.156 and 0.225; P [ 0.05 after Bonferroni correction). Together, these results suggest that limited, but not complete, isolation of populations associated with habitat heterogeneity offers a reasonable explanation for the B. bifarius data.

Conserv Genet (2013) 14:1099–1110

Discussion The use of range dimensions and spatial data on environmental heterogeneity to estimate gene flow corridors among populations can often improve inferences about population genetic processes (Cushman et al. 2006; McRae and Beier 2007; Goulson et al. 2010; Apodaca et al. 2012; Devitt et al. 2013). We show that ENMs and isolation by resistance analysis can explain much of the genetic differentiation among B. bifarius populations, providing additional insights over the use of geographic distance and clustering methods alone (Lozier et al. 2011). In our experience, the majority of B. bifarius specimens from the field or natural history collections can be readily assigned subjectively to one of three morphological groups based on their geographic provenance (Fig. 1), including the two widely accepted subspecies B. bifarius bifarius (red-banded) and B. bifarius nearcticus (black-banded)(Thorp et al. 1983, Koch et al. 2012), and a third phenotype that likely reflects the vancouverensis form described by Cresson (1878) and Viereck et al. (1904) but later grouped with B. bifarius bifarius (Stephen 1957). Given these color phenotypes, and comparable patterns of regional genetic clustering (Fig. 1, 3 Supplementary Fig. 2; Lozier et al. 2011), we were thus interested in testing whether reproductive isolation could be responsible for structuring diversity in B. bifarius, in contrast to a null model of migration-drift equilibrium (McRae and Beier 2007; Frantz et al. 2009). Color pattern morphometrics for populations used in genetic analyses initially revealed five distinct phenotypic clusters in B. bifarius, which were subsequently collapsed to three groups in two-dimensional space as three of the initial five were superimposed. These final three groups corresponded directly to the above morphological groups, with the nearcticus black-banded form comprising the three overlapping clusters, suggesting greater variation within this form. As suggested by this sub-clustering, the subsequent inclusion of geographically intermediate populations indicates greater complexity, with a number of intermediate pigment clusters obscuring the much simpler three color-group scenario. Such patterns support the correlation of genetic distance with resistance that suggests that gene flow in B. bifarius is largely shaped by rangewide variation in environmental suitability, with populations connected by narrow bands of suitable habitat showing larger effective separation than those in more homogenous regions (Fig. 4). Our results do not necessarily indicate that there have never been important periods of historical isolation that have shaped variation in B. bifarius, but the linear IBR relationship and diversity of intermediate coloration suggests that extant populations are governed to a large extent by equilibrium processes.

1107

Equilibrium does not imply an absence of ecologically important limits to dispersal across the B. bifarius range, however. Incorporating the ENM-based B. bifarius distribution into genetic analyses demonstrates that both overly dry habitats and bodies of water are important for restricting gene flow. Linear features and climatic gradients often constitute dispersal barriers for terrestrial organisms (Storfer et al. 2010), but it does not appear that these factors completely impede gene flow, evidenced by the relatively low r-values for partial Mantel tests including both cluster membership and resistance. We suspect the clustering observed in Geneland and Structure can be largely attributed to the ‘‘clines versus cluster’’ dilemma (Frantz et al. 2009; Guillot et al. 2009), where populations along a spatial gradient can be assigned to unique genetic clusters in the absence of true barriers. The fairly weak clustering observed from structure (Lozier et al. 2011) suggested the potential for ongoing gene flow among geographic regions, but insufficient genetic data or incomplete lineage sorting might conceivably affect such analyses in large recently isolated populations (Latch et al. 2006; Orozco-terWengel et al. 2011). The present results provide further evidence that ongoing gene flow is an important part of the genetic structure of this species, and that species-level taxonomic revisions are likely unwarranted without additional evidence. Such ongoing gene flow may also explain the intermediate phenotypes falling between color pattern cluster ellipses (Fig. 3). In contrast, B. ephippiatus, another new world species that exhibits substantial color polymorphism, comprises a set of distinct lineages associated with both geographic isolation and color pattern (Duennes et al. 2012). This study adds to the mounting evidence for a complex relationship between population structure and color pattern variation, and suggests that insights into phylogeography and species delimitation of Bombus with respect to color pattern variability will benefit from in depth analyses on a case-by-case basis. The lack of evidence for strongly isolated populations certainly raises questions regarding the maintenance of color polymorphisms in B. bifarius. The convergence of multiple co-distributed species on similar color patterns in western NAm (Stephen 1957; Thorp et al. 1983) strongly suggests an adaptive function, such as Mu¨llerian mimicry or thermoregulation (Williams 2007). At this point it is difficult to address the specific adaptive mechanisms driving the spatial distribution of these polymorphisms in B. bifarius. Additional research should focus on the extent to which the distribution of color patterns reflects local adaptation versus primarily neutral processes. The study of the San Juan Island vancouverensis phenotype could be particularly interesting, as despite proximity to the blackbanded nearcticus color form, these B. bifarius more closely resemble the eastern red-banded bifarius phenotype

123

1108

(Fig. 1, 3). Genome-wide analyses may enable detection of genetic variants associated more closely with color than with space, which might represent candidate genome regions underlying adaptive genetic structure. The relatively homogenous ranges of many bumble bees may explain why, in contrast to B. bifarius, little genetic structure has been identified in several other NAm Bombus species (Lozier et al. 2011). For example, ENMs for three eastern species (B. impatiens, B. bimaculatus, and B. pensylvanicus) predict large continuous distributions with few environmental barriers to dispersal (Cameron et al. 2011). Isolation on islands appears important in eastern species, but there was no evidence for a relationship between geographic and genetic distances such as that observed for B. bifarius (Lozier et al. 2011). It may thus be the case that substantial habitat heterogeneity is required to produce detectable genetic structure in bumble bees collected at broad geographic scales (Wasserman et al. 2010). Similar patterns are apparent in Europe, where it is often necessary to examine populations isolated by islands or mountain ranges to identify strong intraspecific differentiation (Estoup et al. 1996; Widmer et al. 1998; Widmer and SchmidHempel 1999; Darvill et al. 2010; Goulson et al. 2011). This is not to say that sampling at the appropriate scale and incorporating relevant environmental data will not prove useful for studying patterns of dispersal at finer scales in other bumble bees. Our study focused on broad scale effects of climate on gene flow in B. bifarius, but lacks the necessary resolution to identify how local landscape features may impact dispersal. More intensive local studies with appropriately scaled landscape data would be worthwhile (e.g. Murphy et al. 2010; Wasserman et al. 2010). For example, Goulson et al. (2011) applied Circuitscape models incorporating bathymetric data to infer patterns of isolation and drift associated with sea level changes in B. hortorum populations. Fine-scale analyses incorporating other land use variables in a causal-modeling framework (e.g. Cushman et al. 2006; Wasserman et al. 2010) may provide valuable data on how dispersal in bumble bees is affected by particular components of urban or agricultural landscape (Jha and Kremen 2013). Efforts with a focus on spatial relationships among individuals, rather than populations, may be particularly valuable (Knight et al. 2009; Goulson et al. 2010). From a conservation perspective, our results also suggest that research into the effects of increased habitat isolation associated with climate change is warranted. Although B. bifarius is not currently in obvious decline (Cameron et al. 2011), montane species are predicted to become increasingly isolated with climate change as populations track suitable conditions upslope (Hill et al. 2010; McCain and Colwell 2011). Given that gene flow in B. bifarius can at least in part be predicted by the distribution

123

Conserv Genet (2013) 14:1099–1110

of climatically suitable habitat, environmental change may result in increased genetic isolation (Wilson et al. 2005). In terms of environmental drivers of gene flow, it is interesting that precipitation consistently contributed most to the B. bifarius ENM as opposed to temperature, which was our a priori expectation for this primarily montane species. It may be that overly wet conditions limit opportunities for foraging (Williams 1986), or contribute to flooding of suitable nest sites for principally underground nesting species like B. bifarius (Sakagami 1976; Harder 1986; Williams and Osborne 2009). At the opposite extreme, limitations from overly dry conditions may reflect biotic interactions, as sufficient yearly precipitation, including snow, is likely to be crucial for providing floral resources required throughout the colony lifespan (e.g., Aldridge et al. 2011). Further research into the effects of variable precipitation on bee distributions is warranted, particularly in light of potential climatic changes that could increasingly disrupt the synchrony of plant-pollinator interactions.

Conclusions Taking range dimensions and variation in environmental suitability into account is likely to improve inferences about population structure at the range-wide scale in bumble bees that occur in montane environments. B. bifarius exhibits striking color pattern differentiation and regional genetic clustering that has led to hypotheses about sub-specific or specific taxonomic status of populations. However, accounting for environmental suitability across the landscape suggests that ongoing gene flow along restricted spatial dimensions can explain observed genetic patterns. The ENM approach used here requires extensive spatially referenced observation data for the taxa of interest, but the increasing availability of digitized natural history collections should alleviate this limitation for many Bombus species (Koch and Strange 2009; Cameron et al. 2011; Bartomeus et al. 2011). Likewise, the increasing availability of high-resolution land use and environmental data will likely make studies of dispersal or gene flow possible at finer scales. Understanding how gene flow is limited by physical and environmental factors at both broad and fine spatial scales will be important for predicting how these important pollinators may respond to changing climate and land use practices in the coming century. Acknowledgments We thank the many curators and institutions listed in Cameron et al. (2011) for access to B. bifarius specimens used for our ENM. We thank J. Knoblett for assistance in generating molecular data, and S. Jha and reviewers for comments on an earlier version of this manuscript. The data for this research was supported in part by a grant from the United States Department of Agriculture (CSREES-NRI 2007-02274)

Conserv Genet (2013) 14:1099–1110

References Abella´n P, Arribas P, Svenning J-C (2011) Geological habitat template overrides late quaternary climate change as a determinant of range dynamics and phylogeography in some habitatspecialist water beetles. J Biogeogr 39:970–983 Aldridge G, Inouye DW, Forrest JRK, Barr WA, Miller-Rushing AJ (2011) Emergence of a mid-season period of low floral resources in a montane meadow ecosystem associated with climate change. J Ecol 99:905–913 Apodaca JJ, Rissler LJ, Godwin JC (2012) Population structure and gene flow in a heavily disturbed habitat: implications for the management of the imperilled Red Hills salamander (Phaeognathus hubrichti). Cons Gen 13:913–923 Balkenhol N, Waits LP, Dezzani RJ (2009) Statistical approaches in landscape genetics: an evaluation of methods for linking landscape and genetic data. Ecography 32:818–830 Bartomeus I, Ascher JS, Wagner D, Danforth BN, Colla S, Kornbluth S, Winfree R (2011) Climate-associated phenological advances in bee pollinators and bee-pollinated plants. P Natl Acad SciBiol 108:20654–20659 Cameron SA, Lozier JD, Strange JP, Koch JB, Cordes N, Solter LF, Griswold TL (2011) Patterns of widespread decline in North American bumble bees. P Natl Acad Sci-Biol 108:662–667 ´ , Crossley J, Schmidt H, Carolan JC, Murray TE, Fitzpatrick U Cederberg B, McNally L, Paxton RJ, Williams PH, Brown MJF (2012) Colour patterns do not diagnose species: quantitative evaluation of a DNA barcoded cryptic bumblebee complex. PLoS One 7:e29251 Chan LM, Brown JL, Yoder AD (2011) Integrating statistical genetic and geospatial methods brings new power to phylogeography. Mol phylogenet evol 59:523–537 Cushman SA, Landguth EL (2010) Spurious correlations and inference in landscape genetics. Mol Ecol 19:3592–3602 Cushman SA, McKelvey KS, Hayden J, Schwartz MK (2006) Gene flow in complex landscapes: testing multiple hypotheses with causal modeling. Am Nat 168:486–499 Darvill B, O’Connor S, Lye GC, Waters J, Lepais O, Goulson D (2010) Cryptic differences in dispersal lead to differential sensitivity to habitat fragmentation in two bumblebee species. Mol Ecol 19:53–63 Davis AK, Cope N, Smith A, Solensky MJ (2007) Wing color predicts future mating success in male Monarch butterflies. Ann Entomol Soc Am 100:339–344 Devitt TJ, Devitt SEC, Hollingsworth BD, McGuire JA, Moritz C (2013) Montane refugia predict population genetic structure in the Large-blotched Ensatina salamander. Mol Ecol online early. doi:10.1111/mec.12196 Duennes MA, Lozier JD, Hines HM, Cameron SA (2012) Geographical patterns of genetic divergence in the widespread Mesoamerican bumble bee Bombus ephippiatus (Hymenoptera: Apidae). Mol Phyl Evol 64:219–231 Elith J, Leathwick JR (2009) Species distribution models: ecological explanation and prediction across space and time. Ann Rev Ecol Evol S 40:677–697 Elith J, Phillips SJ, Hastie T, Dudı´k M, Chee YE, Yates CJ (2010) A statistical explanation of MaxEnt for ecologists. Divers Distrib 17:43–57 Estoup A, Soignac M, Cornuet JM, Goudet J, Scholl A (1996) Genetic differentiation of continental and island populations of Bombus terrestris (hymenoptera: apidae) in Europe. Mol Ecol 5:19–31 Falush D, Stephens M, Pritchard JK (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164:1567–1587

1109 Foll M, Gaggiotti O (2006) Identifying the environmental factors that determine the genetic structure of populations. Genetics 174: 875–891 Frantz AC, Cellina S, Krier A, Schley L, Burke T (2009) Using spatial Bayesian methods to determine the genetic structure of a continuously distributed population: clusters or isolation by distance? J Appl Ecol 46:493–505 Gaggiotti OE, Foll M (2010) Quantifying population structure using the F-model. Mol Ecol Res 10:821–830 Goslee SC, Urban DL (2007) The ecodist package for dissimilaritybased analysis of ecological data. J Stat Softw 22:1–19 Goulson D (2010) Impacts of non-native bumblebees in Western Europe and North America. Appl Entomol Zool 45:7–12 Goulson D, Lepais O, O’Connor S, Osborne JL, Sanderson RA, Cussans J, Goffe L et al (2010) Effects of land use at a landscape scale on bumblebee nest density and survival. J Appl Ecol 47:1207–1215 Goulson D, Kaden JC, Lepais O, Lye GC, Darvill B (2011) Population structure, dispersal and colonization history of the garden bumblebee Bombus hortorum in the Western Isles of Scotland. Conserv Genet 12:867–879 Guillot G (2005) Geneland: a computer package for landscape genetics. Mol Ecol Notes 5:712–715 Guillot G, Leblois R, Coulon A, Frantz AC (2009) Statistical methods in spatial genetics. Mol Ecol 18:4734–4756 Harder LD (1986) Influences on the density and dispersion of bumble bee nests (Hymenoptera: Apidae). Ecography 9:99–103 Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. Int J Climatol 25:1965–1978 Hill JK, Griffiths HM, Thomas CD (2010) Climate change and evolutionary adaptations at species’ range margins. Ann Rev Entomol 56:143–159 Hutchison DW, Templeton AR (1999) Correlation of pairwise genetic and geographic distance measures: Inferring the relative influences of gene flow and drift on the distribution of genetic variability. Evolution 53:1898 Jha S, Kremen C (2013) Urban land use limits regional bumble bee gene flow. Mol Ecol 22:2483–2495 In press Jiggins CD, Naisbit RE, Coe RL, Mallet J (2001) Reproductive isolation caused by colour pattern mimicry. Nature 411:302–305 Jost L (2008) GST and its relatives do not measure differentiation. Mol Ecol 17:4015–4026 Knight ME, Osborne JL, Sanderson RA, Hale RJ, Martin AP, Goulson D (2009) Bumblebee nest density and the scale of available forage in arable landscapes. Insect Conserv Diver 2:116–124 Koch JB, Strange JP (2009) Constructing a species database and historic range maps for North American bumble bees (Bombus sensu stricto Latreille) to inform conservation decisions. Uludag Bee Journal 9:97–108 Koch J, Strange JP, Williams P (2012) Bumble bees of the Western United States. USDA Forest Service Research Notes Landguth EL, Fedy BC, Oyler-McCance SJ, Garey AL, Emel SL, Mumma M, Wagner HH, Fortin M-J, Cushman SA (2011) Effects of sample size, number of markers, and allelic richness on the detection of spatial genetic pattern. Mol Ecol Res 12:276–284 Latch EK, Dharmarajan G, Glaubitz JC, Rhodes OE (2006) Relative performance of Bayesian clustering software for inferring population substructure and individual assignment at low levels of population differentiation. Cons Gen 7:295–302 Leache´ AD (2011) Multi-locus estimates of population structure and migration in a fence lizard hybrid zone. PLoS ONE 6:e25827. doi:10.1371/journal.pone.0025827

123

1110 Legendre P, Fortin M-J (2010) Comparison of the Mantel test and alternative approaches for detecting complex multivariate relationships in the spatial analysis of genetic data. Mol Ecol Res 10:831–844 Lozier JD, Strange JP, Stewart IJ, Cameron SA (2011) Patterns of range-wide genetic variation in six North American bumble bee (Apidae: Bombus) species. Mol Ecol 23:4870–4888 Manel S, Schwartz MK, Luikart G, Taberlet P (2003) Landscape genetics: combining landscape ecology and population genetics. Trends in Ecol Evol 18:189–197 McCain CM, Colwell RK (2011) Assessing the threat to montane biodiversity from discordant shifts in temperature and precipitation in a changing climate. Ecol Lett 14:1236–1245 McRae BH (2006) Isolation by resistance. Evolution 60:1551–1561 McRae BH, Beier P (2007) Circuit theory predicts gene flow in plant and animal populations. P Natl Acad Sci-Biol 104:19885–19890 McRae BH, Dickson BG, Keitt TH, Shah VB (2008) Using circuit theory to model connectivity in ecology, evolution, and conservation. Ecology 89:2712–2724 Murphy MA, Evans JS, Storfer A (2010) Quantifying Bufo boreas connectivity in Yellowstone National Park with landscape genetics. Ecology 91:252–261 Orozco-terWengel P, Corander J, Schlo¨tterer C (2011) Genealogical lineage sorting leads to significant, but incorrect Bayesian multilocus inference of population structure. Mol Ecol 20:1108–1121 Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecol Model 190:231–259 R Core Team (2012) R A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/ Rasmont P, Coppe´e A, Michez D, De Meulemeester T (2008) An overview of the Bombus terrestris (L. 1758) subspecies (Hymenoptera: Apidae). Ann Soc Entomol Fr 44:243–250 Raufaste N, Rousset F (2001) Are partial Mantel tests adequate? Evolution 55:1703–1705 Rissler L, Apodaca J (2007) Adding more ecology into species delimitation: ecological niche models and phylogeography help define cryptic species in the black salamander (Aneides flavipunctatus). Syst Biol 56:924–942 Rousset F (1997) Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145: 1219–1228 Rousset F (2008) genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Res 8:103–106 Sakagami SF (1976) Specific differences in the bionomic characters of bumblebees. A comparative review. Jour Fac Sci Hokkaido Univ Seri VI, Zool 20:390–447 Schwartz MK, Copeland JP, Anderson NJ, Squires JR, Inman RM, Mckelvey KS, Pilgrim KL, Waits LP, Cushman SA (2009) Wolverine Gene Flow across a Narrow Climatic Niche. Ecology 90:3222–3232 Shah VB McRae BH (2008) Circuitscape: a tool for landscape ecology. Proceedings of the 7th Python in Science Conference (SciPy 2008):62–65

123

Conserv Genet (2013) 14:1099–1110 Slatkin M (1993) Isolation by distance in equilibrium and nonequilibrium populations. Evolution 47:264–279 Stephen WP (1957) Bumble Bees of Western America (Hymenoptera: Apoidea). Technical Bulletin 40, Agricultural Experiment Station, Oregon State College, Corvallis, OR, USA Storfer A, Murphy MA, Evans JS, Goldberg CS, Robinson S, Spear SF, Dezzani R, Delmelle E, Vierling L, Waits LP (2007) Putting the ‘‘landscape’’ in landscape genetics. Heredity 98:128–142 Storfer A, Murphy MA, Spear SF, Holderegger R, Waits LP (2010) Landscape genetics: where are we now? Mol Ecol 19:3496–3514 Strange JP (2010) Nest initiation in three North American bumble bees (Bombus): gyne number and presence of honey bee workers influence establishment success and colony size. J Insect Sci 10:1–11 Thorp RW, Horning DSJ, Dunning LL (1983) Bumble bees and cuckoo bumble bees of California. Bulletin of the California Insect Survey 23. University of California Press, Berkeley CA Viereck HL, Cockerell TDA, Titus ESG, Crawford JC Jr, Swenk MH (1904) Synopsis of bees of Oregon. Washington and Vancouver. Can Entomol 36(4):93–100 Warren DL (2012) In defense of ‘‘niche modeling’’. Trends Ecol Evol 27:497–500 Wasserman TN, Cushman SA, Schwartz MK, Wallin DO (2010) Spatial scaling and multi-model inference in landscape genetics: Martes americana in northern Idaho. Landscape Ecol 25:1601–1612 Widmer A, Schmid-Hempel P (1999) The population genetic structure of a large temperate pollinator species, Bombus pascuorum (Scopoli) (Hymenoptera: Apidae). Mol Ecol 8:387–398 Widmer A, Schmid-Hempel P, Estoup A, Scholl A (1998) Population genetic structure and colonization history of Bombus terrestris s.l. (Hymenoptera: Apidae) from the Canary Islands and Madeira. Heredity 81:563–572 Williams PH (1986) Environmental change and the distributions of British bumble bees (Bombus Latr.). Bee World 67:50–61 Williams P (2007) The distribution of bumblebee colour patterns worldwide: possible significance for thermoregulation, crypsis, and warning mimicry. Biol J Linn Soc 92:97–118 Williams PH, Osborne JL (2009) Bumblebee vulnerability and conservation world-wide. Apidologie 40:367–387 Williams PH, Brown MJF, Carolan JC, An J, Goulson D et al (2012a) Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae). Syst Biodivers 10:21–56 Williams PH, An J, Brown MJF, Carolan JC, Goulson D, Huang J, Ito M (2012b) Cryptic bumblebee species: consequences for conservation and the trade in greenhouse pollinators. PLoS ONE 7:e32992 Wilson RJ, Gutie´rrez D, Gutie´rrez J, Martı´nez D, Agudo R, Monserrat VJ (2005) Changes to the elevational limits and extent of species ranges associated with climate change. Ecol Lett 8(11):1138–1146