Quaternary climate instability is correlated with patterns of population genetic variability in Bombus huntii

Abstract Climate oscillations have left a significant impact on the patterns of genetic diversity observed in numerous taxa. In this study, we examine the effect of Quaternary climate instability on population genetic variability of a bumble bee pollinator species, Bombus huntii in western North America. Pleistocene and contemporary B. huntii habitat suitability (HS) was estimated with an environmental niche model (ENM) by associating 1,035 locality records with 10 bioclimatic variables. To estimate genetic variability, we genotyped 380 individuals from 33 localities at 13 microsatellite loci. Bayesian inference was used to examine population structure with and without a priori specification of geographic locality. We compared isolation by distance (IBD) and isolation by resistance (IBR) models to examine population differentiation within and among the Bayesian inferred genetic clusters. Furthermore, we tested for the effect of environmental niche stability (ENS) on population genetic diversity with linear regression. As predicted, high‐latitude B. huntii habitats exhibit low ENS when compared to low‐latitude habitats. Two major genetic clusters of B. huntii inhabit western North America: (a) a north genetic cluster predominantly distributed north of 28°N and (b) a south genetic cluster distributed south of 28°N. In the south genetic cluser, both IBD and IBR models are significant. However, in the north genetic cluster, IBD is significant but not IBR. Furthermore, the IBR models suggest that low‐latitude montane populations are surrounded by habitat with low HS, possibly limiting dispersal, and ultimately gene flow between populations. Finally, we detected high genetic diversity across populations in regions that have been climatically unstable since the last glacial maximum (LGM), and low genetic diversity across populations in regions that have been climatically stable since the LGM. Understanding how species have responded to climate change has the potential to inform management and conservation decisions of both ecological and economic concerns.


| INTRODUC TI ON
Geographic instability of ecosystems due to Quaternary climate change has left a lasting imprint on the compostion and diversity of populations and species across the planet (Hewitt, 1996(Hewitt, , 2000. Specifically, climate oscillations since the Pleistocene have facilitated both population divergence and speciation through isolation and recolonization of suitable habitats (Callahan et al., 2013;Carvalho et al., 2017;Galbreath, Hafner, Zamudio, & Agnew, 2010; Gutiérrez-Rodríguez, Barbosa, & Martínez-Solano, 2016;Knowles, 2000). A decrease in the geographic spread of suitable habitat over time may lead to a population range contraction, cascading toward a population bottleneck, genetic drift, and a possible loss of genetic diversity (Pauls et al., 2013). Alternatively, the geographic expansion of suitable habitat over time may facilitate a population expansion, which may also lead to a loss of genetic diversity due to founder effect, as the establishing population is typically made up of a small number of colonizing individuals (Pauls et al., 2013). However, colonization into new suitable habitat may also attract individuals from a diverse pool of populations and result in an increase in population genetic admixture (Ortego, Gugger, & Sork, 2015). Understanding how biodiversity responds to environmental change has the potential to inform effective management decisions for species of ecological and economic concern.
The availability of microsatellite markers and environmental niche modeling techniques provides the opportunity to examine the effects of late Pleistocene climate variability on population genetic variability (Callahan et al., 2013;Gutiérrez-Rodríguez et al., 2016;López-Uribe, Zamudio, Cardoso, & Danforth, 2014). There is converging evidence that species have responded to climate oscillations through either geographic expansion or contraction, depending on their associated life history traits and ecological demands (Callahan et al., 2013;Galbreath et al., 2010;López-Uribe et al., 2014). Small rodents with limited dispersal capacity due to narrow bioclimatic niches exhibit strong population divergence as suggested by highly conserved gene regions during the Pleistocene (Galbreath, Hafner, & Zamudio, 2009). However, animals with great dispersal ability and broad bioclimatic niches have received less attention until recently (Françoso et al., 2016;López-Uribe et al., 2014). Furthermore, many studies have examined species with strong bioclimatic specialization (e.g., montane/ alpine specialist), with few studies examining species that may exhibit a degree of adaptation to regional climate variability. Correlating climate variability with neutral population genetic variability has the potential to build an inference on the role of global environmental change as a driver of population isolation and genetic drift.
Bumble bees (Hymenoptera: Apidae, Bombus) are an appropriate organismal model for studying the effects of climate oscillations on contemporary bioclimatic niches, range dynamics, and populationlevel divergence during the Quaternary. They are large-bodied insects and have the capacity to disperse over several kilometers in search of food, nesting sites, and hibernacula (Jha, 2015;Lozier, Strange, & Koch, 2013;Woodard et al., 2015). Furthermore, bumble bees are densely covered in setae and can fly at low temperatures by warming their wing muscles prior to flight (Heinrich & Esch, 1994). They are dependent on pollen and nectar from flowers to feed developing brood. Thus, dispersal and colonization during the Pleistocene have likely been affected by changes in the distribution of food plants. In western North America, several bumble bee species are associated with temperate, alpine environments that can act as sky islands (Heald, 1951;Lozier et al., 2013). These sky islands and adjacent valleys may have served as refugia for bumble bees during climate oscillations where they may have tracked their preferred habitat across different mountain provinces (Galbreath et al., 2009).
In this study, we examine patterns of contemporary population genetic variability of the Hunt bumble bee, Bombus huntii Greene, 1860 across its geographic range. The latitudinal distribution of B. huntii extends from the southern edge of Canada, south throughout the Intermountain West, and Front Range of the Colorado Rocky Mountains to the Trans-Mexican Volcanic Belt province in southern Mexico (Koch, Strange, & Williams, 2012;Labougle, 1990;Thorp, Horning, & Dunning, 1983). The longitudinal range of the species is primarily bound by the crest of the Sierra Nevada and Cascade Mountains in the west, and the Black Hills of South Dakota in the east ( Figure 1). Although populations of B. huntii have been found east of the montane environment of South Dakota, they have been documented in low abundance relative to surveys of bumble bee communities in Nevada, Utah, Colorado, Montana (Dolan et al., 2016;Koch et al., 2015). Finally, B. huntii and its sibling species B. vosnesenskii diverged from their most recent common ancestor by the early Pliocene (~5 mya) (Cameron, Hines, & Williams, 2007;Hines, 2008). Thus, contemporary population genetic diversity of B. huntii can be investigated using climate scenarios estimated for the Pleistocene.
The primary objective of this study was to examine the role of Quaternary climate oscillations on population genetic variability across B. huntii populations. First, we construct environmental niche models (ENMs) based on contemporary occurrence records and bioclimatic variables. Then, operating under the principle of niche conservatism (Peterson, Sober, & Sanchez-Cordero, 1999), we generated a habitat suitability (HS) map of B. huntii during the last glacial maximum (LGM) (~ca. 22,000 years before present, ybp). We predict that high-latitude populations will exhibit greater niche instability compared to low-latitude populations since the LGM due to the presence of the glaciated regions in northern latitudes (Callahan et al., 2013). Second, we used Bayesian inference to test for differences in genetic structure across contemporary B. huntii populations.
We predict that high-latitude populations will exhibit high genetic diversity as they are distributed across a broad elevation gradient and will be able to disperse throughout the environment and maintain gene flow, whereas low-latitude populations are restricted to high-elevation environments which likely impede gene flow (Koch, Looney, Sheppard, & Strange, 2017). Third, we test for the effect of geographic distance and environmental resistance on contemporary patterns of population differentiation across B. huntii with isolation by distance (IBD) and isolation by resistance models (IBR), respectively. Bumble bees are sensitive to environmental variation; thus, we predict that IBR will be a better predictor than IBD for examining population differentiation (Koch et al., 2017;Lozier et al., 2013). The goal of examining both IBR and IBD is to determine whether B. huntii genetic diversity estimates are best predicted by HS or geographic distance, respectively. Finally, we combined the contemporary and LGM HS maps to estimate an environmental niche stability (ENS) map to test the hypothesis that climate instability predicts contemporary genetic diversity across B. huntii populations. We predict that patterns of genetic diversity can be explained by ENS patterns, with stable areas (low-latitude regions) associated with low genetic diversity and unstable areas (high-latitude regions) associated with high genetic diversity (Callahan et al., 2013).

| Quaternary environmental niche modeling
To estimate the geographic distribution of HS of B. huntii throughout its endemic range at present and during the Pleistocene, ENMs were constructed under the principle of maximum entropy with MAXENT v3.3.1 (Phillips, Dudík, & Schapire, 2004). The software program MAXENT uses presence-only georeferenced locality records and random background points sampled from the study extent to estimate the distribution of the species that is closest to uniform (=maximum entropy) under the suite of independent variables (i.e., bioclimatic variables) supplied to the model (Elith et al., 2011). Georeferenced distribution records were queried in the Global Biodiversity Information Facility database (http://gbif.org) and filtered for unique spatial coordinates ( Figure 1) (Supporting Information Appendix S1). Additional records from a contemporary survey of Mexican bumble bees were included in the final set of georeferenced localities of B. huntii and are provided in Table 1.
Specimen occurrence records that do not agree with historic and contemporary range maps of B. huntii were filtered out of the final dataset (Koch et al., 2012;Thorp et al., 1983;Williams, Thorp, Richardson, & Colla, 2014).
Occurrence records were aggregated with 18 spatially explicit bioclimatic variables representing contemporary conditions F I G U R E 1 Natural history (NH) and current survey records of Bombus huntii throughout western North America. Current survey records = yellow point enclosing black point, NH record = white point. Except for "West/Intermountain West+," all ecoregions presented here are adapted from Olsn et al. (2001). The gross delineation of the "West/Intermountain West+" on this map is modified from maps developed for waterbirds by Ivey and Herzinger (2006). Locality data are presented in decimal latitude and longitude with the WGS1984 coordinate system and the North American Albers Equal Area Conic geographic projection  from the WorldClim v1.4 Bioclim dataset (2.5 arc minutes) (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005 Rather than randomly select a variable for the analyses, we chose to retain the variables that reflect annual and seasonal trends in precipitation and temperature as bumble bees are primarily active during the summer months, but also need to survive winter hibernation (Lozier et al., 2013). ENMs in MAXENT were constructed with default parameters to generate a logistic output (a measure of relative HS), averaged over 100 replicates with a subsampling scheme to evaluate model performance (75% train, 25% test). The Operating under the principle of niche conservatism (Peterson et al., 1999), we predicted the distribution of B. huntii during the TA B L E 1 Location description, sample size (N), and colony assignment of Bombus huntii in this study. Region assignment is inferred using a Bayesian genetic cluster assignment test with 11 microsatellite markers in STRUCTURE v 2.3.4 (Pritchard et al., 2000). Assignment of individuals to a colony is made with the program COLONY 2.0 (Jones & Wang, 2010)  LGM (~22,000 ybp) using the constraints of the contemporary bioclimatic associations of the B. huntii ENM mapped to paleoclimate data available within the WorldClim database (Hijmans et al., 2005).
We used the CCSM4 fully coupled global climate model to estimate the B. huntii LGM environmental niche map. To identify potential climate refugia over the Quaternary, we added contemporary and LGM HS maps together to produce a ENS map. The raster output of the ENS map was standardized by dividing the calculated values by the raster's maximum value, producing values ranging from 0 (low niche stability) to 1 (high niche stability).

| Contemporary taxon sampling
Female worker B. huntii were collected across 33 sampling sites throughout North America from 2008 to 2015 ( Figure 1). This sampling regime captured a major portion of the species' range in western North America (Labougle, 1990;Thorp et al., 1983). Bumble bees were sampled with a diversity of methods including sweep netting, colored bee bowls, and blue vane traps. Given that the aim of our study was to examine population genetic diversity and structure of wild B. huntii, we elected to pool certain sampling sites together if they were <9 km from each other. This decision was made based on known biological properties of dispersal in bumble bees, in that workers of the sister species, B. vosnesenskii, are estimated to disperse approximately 9 km from her nest (Jha & Kremen, 2013).
We used this distance threshold as there are no current data on the landscape genetics of wild B. huntii. Furthermore, barriers to gene flow have been found to be influenced by both land-use change and bioclimatic variability in B. vosnesenskii and B. bifarius (Jha & Kremen, 2013;Lozier et al., 2013). Thus, limiting the distance for pooling sites insures that we did not arbitrarily pool populations together and potentially alter genetic variability estimates.

| DNA extraction and microsatellite genotyping
Within each of the 33 pooled sampling sites, an average of 12.25 (±1.29 SE) female bees were collected (n = 380). DNA was extracted from the mid-leg of each bumble bee using a modified-Chelex100 ™ protocol described in Strange, Knoblett, and Griswold (2009)

| Hardy-Weinberg equilibirium, linkage disequilibirium, and colony assignment
The probability of null alleles was estimated with the software pro- and linkage disequilibrium (LD) across populations and loci were tested with the web-based software program GENEPOP v4.0.10 (Raymond & Rousset, 1995). Sequential Bonferroni corrections were applied to the HWE and LD p-values estimates to minimize type I errors associated with multiple comparisons for both populations and loci (Rice, 1989). We considered the Bonferroni correction test significant at p ≤ 0.05.
As bumble bees are generally monandrous (Estoup et al., 1995), primitively eusocial, and live in annual colonies, it is possible to capture sibling female workers in the wild. To avoid pseudo-replication within sampling locations, full-siblings were first assigned to colonies with COLONY v2.0 (Jones & Wang, 2010). In the colony-assignment exercise, we set the mistyping error rate to 0.05, based on error rates documented in previous studies (Lozier, Strange, Stewart, & Cameron, 2011), and the sex-determination system to "haplodiploid." If full-siblings were detected in the colony-assignment tests (≥95% genotype similarity), we selected only one representative from each family using a coin toss. Postanalysis mistyping error estimates are described in Supporting Information Appendix S2.

| Population genetic structure and landscape genetics
Population genetic structure was examined with an individual assignment Bayesian clustering method in the software program STRUCTURE v 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). We elected to use the admixture model in STRUCTURE to assign genotypes, which assumes that individuals comprise K unknown clusters (i.e., population), to which an individual can be assigned based on their genotype without a priori delineation of populations. We set the model to run with 20,000 burn-in steps and 100,000 samples, with 10 iterations for each K, where K ranged from 1 to 10.
To determine the optimal K, the distributions of the probability of the data LnP(K) and ΔK (as described by Earl & vonHoldt, 2012;Evanno, Regnaut, & Goudet, 2005) were examined with the webbased software program STRUCTURE HARVESTER 0.6.94 (Earl & vonHoldt, 2012). To account for multimodality associated with individual STRUCTURE runs, we averaged an individual's admixture proportions over the 10 replicates for the best K using the full search algorithm in CLUMPP v1.1.2 (Jakobsson & Rosenberg, 2007). We visualized populationassignment by averaging K-assignments within each population with pie charts in ArcGIS v10.1 (ESRI, 2012).
To compliment our STRUCTURE analysis, we also used the program GENELAND v4.05 (Geneland R library) to infer poplation structure . Unlike the STRUCTURE analysis, we included an examination of the effect of geographic location on individual assignment to cluster K in GENELAND. Locality data were transformed from decimal latitude and longitude to UTM with the PBSmapping R library. We set the GENELAND population assignment model to run with the following parameters: 1,000,000 Monte

| Population genetic diversity
We estimated genetic diversity with two different metrics: (a) effective allelic diversity (AD) and (b) expected heterozygosity (H e ) using Nei's genic diversity metric. Given our unequal sampling across field sites, we used rarefaction to estimate genetic diversity by standardizing our estimates to four gene copies per population (Kalinowski, 2005

| Environmental niche comparisons
We used two sample Wilcoxon rank-sum tests to test for significant bioclimatic differences between the two Bayesian-infered genetic clusters (i.e., North vs. South cluster). Unless otherwise indicated, all statistical tests conducted in our study were executed with the R statistical computing platform with the appropriate statistics libraries (R Core Development Team, 2005).

| Quaternary environmental niche modeling
The 18 available bioclimatic variables were reduced to 10 bioclimatic variables after the assessment of colinearity and were in-  (Figures 1 and 2c).

| Hardy-Weinberg equilibirium, linkage disequilibirium, and colony assignment
Of the original 33 sampling sites (n = 380; hereafter refereed to as "populations"), we retained 26 populations (n = 364) for further analyses as they were represented by more than four individuals. In our analyses of the 26 populations, the locus BTERN02 did not amplify in >50% of the specimens genotyped and was removed from any further data processing. MICRO-CHECKER results indicated that the locus BTMS0059 was suspected to have null alleles in 48% of the populations and was also excluded from further analyses.
B124 showed significant deviation from HWE in the Ada (ADA, LGM Contemporary  (Table 1). A summary on the number of colonies detected at each of our sites is available in Supporting Information Appendix S3.

| Population genetic structure and landscape genetics
Based on the admixture ancestry with correlated allele frequencies model in STRUCTURE, the mean log probability of the data was greatest at K = 9 (LnP (9)    for both IBD and IBR models are presented in Table 3. We did not investigate IBD and IBR models using the Geneland genetic cluster assignments as three of the genetic clusters were composed of one or two sampled populations (clusters 2, 3, and 4), and would be an inappropriate sample size for a robust analysis (Figure 4).  (Pritchard et al., 2000). Pie charts represent fractional assignment probability of each sampling site (i.e., population) to each of the K = 2 genetic clusters. Fractional assignment to the North and South clusters is in orange and blue, respectively. Inset bar graphs represent individual STRUCTURE assignment to each of the K = 2 and K = 3 genetic clusters expressed in decreasing latitude (each bar = 1 individual). Three letter site sampling site codes are defined in Table 1 EDM

| Environmental niche comparisons
Pairwise comparisons between the bioclimatic variables found that eight of the ten variables used to construct the ENMs exhibited significant differences between the North and South+ genetic clusters

| D ISCUSS I ON
Our population genetic study of B. huntii found that populations in regions that have been climatically unstable since the LGM (highlatitude environments, e.g., north of 28°N) exhibit high genetic diversity, whereas regions that have been climatically stable since the LGM (low-latitude environments, e.g., south of 28°N) exhibit low genetic diversity (Figure 6b-c). Throughout its northerly distribution, B. huntii population genetic differentiation is significantly correlated with geographic distance (IBD, Figure 6e), but not bioclimatic resistance (IBR, Figure 6f). Nonsignificant IBR suggests that northern populations of B. huntii are not strongly inhibited by environmental heterogeneity, except at long distances between populations (>2,000 km). Populations distributed throughout the north genetic cluster are associated with much colder annual temperatures, a pattern that is certainly correlated with increasing latitude (Figure 3b).
Furthermore, an examination of the elevation profile of B. huntii ENM in the northern part of its range suggests high HS across a F I G U R E 4 Genetic cluster assignment of 26 Bombus huntii populations in western North America inferred with GENELAND v4.05 (Guillot et al., 2012). Cluster 1 = orange, Cluster 2 = green, Cluster 3 = golden yellow, Cluster 4 = gray, and Cluster 5 = blue. Inset bar graphs represent individual GENELAND assignment to each of the K = 5 genetic clusters expressed in decreasing latitude (each bar = 1 individual). Three letter site sampling site codes are defined in Table 1 EDM at distances greater than >500 km (Lozier et al., 2011). The ability of bumble bees to disperse and forage on a diverse array of pollen sources at great distance is a significant factor in facilitating gene flow across populations (Goulson & Stout, 2001;Jha, 2015;Knight et al., 2005;Moreira, Horgan, Murray, & Kakouli-Duarte, 2015).
Bumble bees are dependent on suitable forage and nesting resources, however, dispersal as a mechanism in facilitating gene flow is limited to the reproductive caste, namely queens and males. While some North American bumble bees exhibit low admixture or population genetic differentiation, other, mostly montane species have been found to exhibit significant IBD and IBR (Hines & Williams, 2012;Jha & Kremen, 2013;Lozier et al., 2013). For example, B. bifarius is primarily associated with montane environments throughout its distribution. However, B. bifarius also occur at low-elevation and offshore islands throughout the Pacific Northwest and Central Coast of California (Koch et al., 2012;Williams et al., 2014). These populations have reduced population genetic diversity and exhibit a degree of phenotypic divergence relative to populations distributed across the Colorado Rocky Mountains, the Sierra Nevada, Cascade, Sawtooth, Bighorn, and Uinta Mountains of western North America (Lozier et al., 2013;Lozier et al., 2016). The variability in HS across these montane provinces has been a major barrier to gene flow, revealing significant IBR patterns across B. bifarius populations.  (Table 3).
Comparable significance of the IBR model to the IBD model across the southern genetic cluster implies that adjacent environments of contemporary B. huntii populations are a major barrier to dispersal, and ultimately gene flow. Examination of the ENMs associated with LGM and contemporary HS suggest that high-elevation montane environments south of 28°N are more suitable for B. huntii than lowelevation environments (Figure 2a-b). Furthermore, the patterns of genetic structure uncovered in Mexican B. huntii populations as they relate to their distribution across the different mountain provinces are like the genetic patterns uncovered by Duennes et al. (2017) in the B. epphipiatus species complex. Our study corroborates contemporary investigations into how the complex climate history of the Quaternary, along with the geographic isolation to mountain provinces, has played a significant role in population divergence and cladogensis in Mesoamerican biodiversity.

| Quaternary climate variability and population genetic patterns
Pleistocene climate variability has influenced contemporary patterns of biodiversity observed throughout North America (Hewitt, 2000). Specifically, the oscillation between cooling and warming periods is implicated to have driven the dynamics of species and communities distributed throughout montane environments in plants (Callahan et al., 2013) and mammals (Galbreath et al., 2009(Galbreath et al., , 2010, as well as their underlying patterns of genetic diversity (Hewitt, 1996). The correlation between ENS and latitude, along with increased population genetic diversity with decreasing ENS, suggests that northern populations are part of the leading edge of the species range, whereas the more genetically differentiated southern populations are part of the rear edge of the species range ( Figure 3) (Hampe & Petit, 2005 F I G U R E 6 (a) Scatterplot and linear models examining the relationship between latitude and environmental niche stability (ENS) of Bombus huntii, (b) the relationship between B. huntii effective allelic diversity (AD) and ENS, and (c) the relationship between B. huntii expected heterozygosity (H e ) and ENS. ENS represents the sum of contemporary  and the last glacial maximum (LGM) habitat suitability (0-1) values predicted by a contemporary ENM (see Figure 3). Stability values closer to 1 represent high HS stability, whereas values closer to 0 represent low HS stability (i.e., unstable) geographic distribution (Figure 2a) Figure 2c).
In addition to relative climate stability, stable rear edge populations are associated with heterogeneous topography that are typically restricted to habitat islands (i.e., sky islands) within a matrix of unsuitable environments (Hampe & Petit, 2005). IBD and IBR models of B. huntii populations in southern genetic cluster show that geographic distance and environmental resistance explain a large portion of the variability found in pairwise estimates of genetic differentiation (Figure 5c-d). Finally, reduced within-population genetic diversity pattern and high genetic differentiation, despite close geographic proximity, are characteristic of stable rear edge populations. Both features of genetic variability (diversity and structure) are readily observed in the patterns of allelic diversity and expected heterozygosity uncovered in our study. We found that high-latitude populations (leading edge) are associated with high within-population genetic diversity, and low-latitude populations (rear edge) are associated with reduced within-population genetic diversity (Figure 6b-c).

| CON CLUS ION
A wealth of phylogeographic studies supports the hypothesis that Quaternary climate variability has shaped range-wide patterns of population genetic diversity and divergence that are observed in extant populations and species. In our study, we discover clear geographic and environmental differences in the contemporary B. huntii populations that are predicated on ENS since the LGM.
The major findings of this study include the following: (a) the B. huntii ENS map predicts low-latitude environments to be more climatically stable in comparison with high-latitude environments over the past 22,000 years; (b) there are at least two genetic clusters (north and south) that a population can be assigned to based on a priori Bayesian analysis, with further evidence for genetic subdivision across Mexican populations found across different mountain ranges; (c) IBD and IBR models are significant across populations in the south genetic cluster, whereas only IBD is significant in the north genetic cluster; (d) genetic diversity increases with latitude and, therefore, decreases within increasing ENS; and (e) there are significant bioclimatic differences between the north and south genetic clusters, with southern populations inhabiting mesic environments with low variation in precipitation across seasons.
The results of this study contribute to our understanding of how historic and contemporary climate may have shaped the patterns of genetic variability observed in insects, an understudied group of animals, relative to vertebrate organisms. Bumble bees are a valuable model to study how climate variability affects patterns of genetic variability as they are highly sensitive to recent climate change. Over the past century, bumble bee populations have responded to a warming climate by emerging earlier in the spring, evolved shorter proboscis, and potentially shifted in latitudinal distribution (Bartomeus et al., 2011;Kerr et al., 2015;Miller-Struttman et al., 2015). Documenting contemporary patterns of F I G U R E 7 Boxplot visualization of contemporary  environmental niche space occupied by contemporary North and South Bombus huntii genetic clusters for 10 bioclimatic variables. Significant differences between the North and South+ populations are tested with a Wilcoxon rank-sum test and identified with a * symbol at the midline of each graph (alpha = 0.05) genetic variability has the potential to inform management and conservation policy for both wild and managed bumble bee pollinator populations.

ACK N OWLED G M ENTS
Our study was supported in part by a research grant from CONABIO (#JE016) and SEP-CONACYT (#106043) to Rémy Vandame and from the United States Department of Agriculture (CSREES-NRI 2007-02274) to James P. Strange.

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

AUTH O R CO NTR I B UTI O N S
J.B.K. and J.P.S. conceived the idea; J.B.K., J.P.S., R.V., J.M-R., and P.S., collected the specimens; J.B.K. generated and analyzed the genetic data, and led the writing of the manuscript. J.P.S. and R.V.
provided editorial feedback.

DATA ACCE SS I B I LIT Y
Microsatellite genotypes and collection locations are available on