Climate driven disruption of transitional alpine bumble bee communities

Pollinators at high elevations face multiple threats from climate change including heat stress, failure to phenological match advancing flower resources and competitive pressure from range‐expanding species of lower elevations. We conducted long‐term multi‐site surveys of alpine bumble bees to determine how phenology of range‐stable and range‐expanding species is responding to climate change. We ask whether bumble bee responses generate mismatches with floral resources, and whether these mismatches in turn promote community disruption and potential species replacement. In alpine environments of the central Rocky Mountains, range‐stable and range‐expanding bumble bees exhibit phenological mismatches with flowering host plants due to earlier flowering of preferred resources under warmer spring temperatures. However, workers of range‐stable species are more canalised in their foraging schedules, exploiting a relatively narrow portion of the flowering season. Specifically, range‐stable species show less variance in phenology in response to temporally and spatially changing conditions than range‐expanding ones. Because flowering duration drives the seasonal abundance of floral resources at the landscape scale, we hypothesize that canalisation of phenology in alpine bumble bees could reduce their access to earlier or later season flowers. Warmer conditions are decreasing abundances of range‐stable alpine bumble bees above the timberline, increasing abundance of range‐expanding species, and facilitating a novel and more species‐diverse bumble bee community. However, this trend is not explained by greater phenological mismatch of range‐stable bees. Results suggest that conversion of historic habitats for cold‐adapted alpine bumble bee species into refugia for more heat‐tolerant congeners is disrupting bumble bee communities at high elevations, though the precise mechanisms accounting for these changes are not yet known. If warming continues, we predict that the transient increase in diversity due to colonization by historically low‐elevation species will likely give way to declines of alpine bumble bees in the central Rocky Mountains.

ble bee responses generate mismatches with floral resources, and whether these mismatches in turn promote community disruption and potential species replacement. In alpine environments of the central Rocky Mountains, range-stable and range-expanding bumble bees exhibit phenological mismatches with flowering host plants due to earlier flowering of preferred resources under warmer spring temperatures. However, workers of range-stable species are more canalised in their foraging schedules, exploiting a relatively narrow portion of the flowering season. Specifically, range-stable species show less variance in phenology in response to temporally and spatially changing conditions than range-expanding ones. Because flowering duration drives the seasonal abundance of floral resources at the landscape scale, we hypothesize that canalisation of phenology in alpine bumble bees could reduce their access to earlier or later season flowers. Warmer conditions are decreasing abundances of range-stable alpine bumble bees above the timberline, increasing abundance of range-expanding species, and facilitating a novel and more species-diverse bumble bee community. However, this trend is not explained by greater phenological mismatch of range-stable bees. Results suggest that conversion of historic habitats for cold-adapted alpine bumble bee species into refugia for more heat-tolerant congeners is disrupting bumble bee communities at high elevations, though the precise mechanisms accounting for these changes are not yet known. If warming continues, we predict that the transient increase in diversity due to colonization by historically low-elevation species will likely give way to declines of alpine bumble bees in the central Rocky Mountains.

K E Y W O R D S
alpine, bumble bees, climate change, counter-gradient selection, phenological mismatch, range shift, transitional communities

| INTRODUC TI ON
Mutualisms are fundamental to species persistence and distribution, yet generally involve partners with large differences in mobility, generating potential disruption under conditions promoting range shift. Pollination mutualisms involve associations between largely sedentary plants and more mobile animals. When environmental conditions undergo rapid degradation, such associations are vulnerable to disruption (Aslan et al., 2013;Keeler et al., 2021). Climate change is driving this process via multiple pathways, with unclear outcomes for biodiversity in novel and transitional communities (Gérard et al., 2020;Shah et al., 2020).
While climate change clearly reshuffles local pollinator communities through loss of extant biodiversity coupled with gain of migrants from warmer regions (Williams & Jackson, 2007), what is less understood is how the two processes interact to shape species composition (Gibson-Reinemer et al., 2015;Shah et al., 2020) We address this gap for bumble bee communities in rapidly warming alpine regions.
Phenological synchrony with floral resources is reduced when partners differ in responses to climate cues (match-mismatch theory; Kerby et al., 2013;Ploquin et al., 2013). Conversely, convergence in cue use between partners could maintain phenological matching despite climate variance (Hegland et al., 2009;Inouye, 2020;Roslin et al., 2021). While conceptually straightforward, documenting such convergence is prone to spurious conclusions if forager and flower abundances are monitored together (Forrest & Thomson, 2011).
Linkages between sequential life history stages also may drive temporal mismatches between bumble bees and flowers. Bumble bee queens overwinter, emerging and establishing nests the following spring, with workers maturing several weeks later. For organisms with such complex life histories, climate events at one life cycle stage may affect the timing of critical functions at later stages. For example, in migratory birds, climate cues in the winter habitat influence phenological synchrony with prey in the breeding range (Dunn & Winkler, 2010). Emerging bumble bee queens avoid starvation during prolonged nest searching by resorbing their ovaries (Sarro et al., 2022), potentially delaying worker production and reducing worker synchrony with floral resources.
As alpine environments warm, range-expanding bumble bees from hotter regions may better track advancing flowering schedules than resident species (Soroye et al., 2020;Wadgymar et al., 2015).
Research on developmental timing in amphibians shows that warmer and longer growing seasons at low elevations select for developmental plasticity, while cold-delayed seasons at higher sites favor rapid, early maturation and more canalised developmental schedules (counter-gradient selection, Berven et al., 1979). If bumble bee development has been shaped similarly by past climate regimes, rangeexpanding and range-stable species may diverge in synchrony with floral resources and experience different levels of resource stress.
We consider two pathways for reconfiguration of bumble bee communities under a warming alpine climate coupled with colonization by range-expanding species. First, we test the hypothesis that discrepancy in responses to changing climate cues between bumble bees and flower resources contributes to partners' phenological mismatch; that the magnitude of mismatch in turn reduces bumble bee abundance; and that such losses are amplified by the lack of phenological flexibility in range-stable alpine taxa. Alternatively, we postulate that regardless of phenological mismatch, range-expanding bumble bees from lower elevations are better suited to a warmer alpine environment than coldadapted, range-stable congeners. Both scenarios predict that colonizing bumble bees will replace alpine residents, but by different pathways. To address these hypotheses, we evaluated (1) climate cues affecting flowering time in the host plant-community, (2) cues driving foraging phenology of colonizing versus resident bumble bee species, (3) impacts of climate cues on synchrony of colonizers and residents with floral resources, and (4) relationship of phenological mismatch to abundances of resident and colonizing species.

| Study system
Subalpine bumble bees have moved upward since the mid-1900s and currently nest above treeline (~3500 m) in the Colorado Rocky Mountains (Miller-Struttmann et al., 2015). In our alpine (above treeline) study sites (Colorado, USA), we classified bumble bees as resident or range-stable if pre-1980 records showed them completing their entire life cycle above treeline (B. kirbiellus and B. lapponicus sylvicola, formerly B. balteatus and B. sylvicola, respectively; Byron, 1980). Observations of B. lapponicus likely included the recently discovered cryptic taxon "incognitus" as it overlaps in distribution and is morphologically indistinguishable (Christmas et al., 2021).
Recent warming has occurred at all our study sites, as evidenced by increasing spring and summer temperatures (Miller-Struttmann et al., 2015) and decreasing episodes of extreme cold spring conditions when temperatures averaged below 0°C in May (15% of 20th century records vs. 6% of the 21st; Table S1 and Figure S1; PRISM, Daly et al., 2008;Strachan & Daly, 2017; see below for details).
Current temperature regimes in our alpine field sites are converging on historical norms in the surrounding subalpine regions, consistent with more rapid warming at high elevations (Shah et al., 2020; Figure S1).

| Spatial and temporal characterization of climate
Climate records were obtained from PRISM point interpolation data. Collinearity among climate variables was weak (Table S2), with the variation explained (R 2 ) ranging from 0.0016% to 30% (maximum variance inflation factor of 1.4).
Climate predictors were chosen to reflect critical conditions for different bumble bee life cycle stages and host-plant reproduction.
Mean daily May temperature was selected along with winter precipitation (sum of precipitation or snowfall from January through May), because together these have pronounced impacts on snowmelt timing and the accumulation of growing degree days for plants and insects in alpine regions. Specifically, spring warming coincides with emergence from torpor and colony initiation for bumble bee queens (Szabo & Pengelly, 1973) and onset of new growth in their host plants (Gezon et al., 2016). Summer maximum temperatures and summer precipitation (respectively, mean of monthly diurnal high temperatures and sum of June, July, and August precipitation) were used as they influence accumulation of energy resources in the nest via production and activity of foraging workers (Maebe et al., 2021) and duration of flowering in host plants via soil and plant water TA B L E 1 Sampling design for data addressing the specific aims of this study. Observations of bumble bees (queens [Q], workers [W]) and their host plants were made at several geographically isolated sites and altitudinal zones, over varying time scales to probe effects of climate change on community dynamics. Because not all species were found in all zones, total replication was reduced in some cases

| Climate cues affecting flowering time in the host plant-community
Twenty-three 2 × 10 m plots were established on Pennsylvania Mountain by P.G. Kevan and surveyed weekly from the onset to the end of the flowering season in 1977-1981. Plots were located nonrandomly to sample prevalent alpine habitat types (Byron, 1980).  in the r2glmm package version 0.1.2 (Jaeger, 2017). Here and elsewhere, we tested whether our data met the assumptions of parametric statistics, and if they did not, we used transformations as explained, when necessary.
To ask whether variation in flowering phenology influences the seasonal accumulation of floral resources, we ran a single mixed-effects multiple regression with species and plot as random effects using the lme in the nlme package version 3.1-149 (Pinheiro et al., 2020). We included day of first flowering and day of peak flowering, as independent fixed effects in the model. Cumulative flower density was the dependent variable and was log transformed to meet the assumption of normality. Test statistics were estimated via the anova.lme function in the nlme package version 3.1-149 (Pinheiro et al., 2020), and the partial R 2 for each explanatory variable was extracted via r2beta in the r2glmm package ver- of the first observation of that species/caste, and peak activity as the DOY when the maximum number of foraging individuals were observed in a given elevational zone. We use sightings rather than estimating phenological events via modeling, because model estimates of first occurrence do not provide biologically realistic results for this system. Weibull estimates (Belitz et al., 2020) were, on average, 12.02 days earlier than first sightings (range: 1.06-60.17 day), predicting emergence when mean daily temperatures were well below that required for flight (Heinrich, 1975;Lundberg, 1980 (Pinheiro et al., 2020), and the partial R 2 for each explanatory variable was extracted via r2beta in the r2glmm package version 0.1.2 (Jaeger, 2017). We similarly compared competing models for timing of worker emergence and peak activity (Table S4). Day of queen emergence did not predict worker emergence timing (F 1,39 = 0.080, p = .78) and was removed from the competing models. Mixed effect ANCOVAs tested for fixed effects of range status, summer maximum temperature, and summer precipitation on timing of worker emergence and peak activity, with species as a random effect. For peak worker activity day, we also included as a fixed effect worker emergence day, to test for a carryover effect. When there was an interaction between range status and a climate variable, we ran separate mixedeffect regressions for each range type (i.e., expanding or stable).
To address whether seasonal patterns of forager activity are more canalised for resident than colonizing bumble bees, we compared curves for temporal variation in forager activity over the season. The number of workers of each species on each date was calculated at the regional level by pooling samples from all locations and elevational zones. The activity curves for each species were tested for kurtosis using the PerformanceAnalytics package version 2.0.4 (Peterson et al., 2020), which quantifies the degree to which the peak of a curve is narrower or broader than expected for a normal distribution. Because our sample size was small (n = 2) for resident species, we tested for greater canalisation (more positive kurtotic scores) in resident than colonizing species using a nonparametric one-tailed Mann-Whitney test (wilcox.Test function in the base R package version 4.0.2; R Development Core Team, 2020).
We then explored the degree to which local adaptation to alpine environments has favored canalisation in phenology, by testing for a relationship between the kurtosis of current (2012-2014) worker activity curves and the mean elevation at which a species was collected historically (1966-1969Macior, 1974).  (Table 1).
For workers of each bumble bee species, we estimated the yearly phenological mismatch between foraging activity and flowering of the host-plant community in each elevational zone. To test for mismatch between worker foraging activity and flowering, we calculated the difference between day of peak worker foraging and day of peak flowering. Negative numbers indicate that worker activity peaked before flowering. We tested whether range status and/or climate predicted the phenological mismatch between peak worker activity and peak flowering by comparing competing models constructed via lme in the nlme package version 3.1-149 (Pinheiro et al., 2020).
Predictor variables included range status, summer maximum temperature, summer precipitation, and their interactions, with species as a random effect (Table S5). Competing models were compared as above. For the top performing model, we calculated test statistics using type III sums of squares via the anova.lme function in the nlme package version 3.1-149 (Pinheiro et al., 2020), and the partial R 2 for each explanatory variable was extracted via r2beta in the r2glmm package version 0.1.2 (Jaeger, 2017).

| Relationship of phenological mismatch to abundances of resident and colonizing species
Changing climate factors could influence abundances of range-stable versus range-expanding bumble bees by generating mismatches between seasonal foraging schedules and flowering times. Cumulative abundance of each bumble bee species was estimated as the sum of the weekly counts of all workers of that species weighted by collection effort over the season within an elevational zone. Relative abundance was calculated as the proportional contribution of a species to the bumble bee forager community. We tested whether cumulative and relative abundances of workers were predicted by alternative combinations of range status, mismatch with peak floral resources, summer maximum temperature, summer precipitation, and their interactions (  (Figure 1a).
Cumulative flower density over the season was influenced by flowering phenology. Flower density decreased with later onset of flowering (F 1,764 = 67.32, p < .0001, partial R 2 = 4.9%; Figure S4a) and increased with later peak flowering (F 1,764 = 57.28, p < .0001, partial R 2 = 4.1%; Figure S4b). These results show that patches of plants initiating flowering early and taking longer to reach peak flowering produced more flowers over the season.

| Cues driving foraging phenology of colonizing versus resident bumble bee species
Timing of queen emergence was best described by one model that included range status, May mean temperature, winter precipitation, and one-way interactions of climate factors with range status (Table S4). Neither range status nor mean May temperature influenced queen emergence schedule (Table 3). Instead, emergence of both resident and colonizing bumble bee queens was delayed by heavy winter precipitation (F = 6.60, p = .014; Table 3; Figure 1b), explaining 12.3% of the variation and indicating that snowmelt timing drives queen emergence schedule in these alpine sites. Thus, bumble bee colony initiation is shaped by different climate cues than first flowering of bumble bee host plants.
Worker emergence was best explained by one model that included range status, summer maximum temperature, summer precipitation, and one-way interactions of climate factors with range status (Table S4). Both interaction terms were marginally statistically significant (Table 3), each accounting for ~3% of the variation ( Table 3). Resident bumble bee workers emerged earlier than those of colonizing species (F = 6.84, p = .040, partial R 2 = 7.6%; Figure S5a). However, this trend was weaker in warmer and wetter summers (Table 3; Figure 1c,d). Summer precipitation significantly advanced day of worker emergence, independent of its interaction with range status (F = 7.79, p = .0068, partial R 2 = 8.4%).

TA B L E 2
Parameter estimates for the models best explaining the effects of climate variables on flowering phenology of the bumble bee host-plant community. Day of first flowering and day of peak flowering were best explained by one model (Table S4); therefore, model averaging was not required and chi-square values from the best performing model are reported. Model comparison was not used for cumulative flower density because we tested a single prediction regarding whether flowering phenology influences seasonal flower production. F-values from the mixed effects multiple regression are reported range status and maximum summer temperature best explained variance in the timing of peak worker activity (Table S4). When workers emerged earlier, they reached peak activity earlier (F = 31.23, p < .0001, partial R 2 = 50%; Table 3; Figure S5b; y = 48.44 + 0.81x) indicating that for both resident and colonizing species, peak worker activity is primarily constrained by emergence time.
Variance in the timing of worker accrual was greater for colonizing bumble bee species than for residents ( Figure 2; for unadjusted density curves, see Figure S6). Resident species tended to have lower kurtosis than colonizing species (W = 1.00, p = .071; Figure S7). Colonizers exhibited shallower peak abundance with more worker activity earlier and later in the summer (Figure 2a; Figure S7). Resident worker numbers accumulated and declined more abruptly as the season progressed. Kurtosis in worker activity schedules changed from negative (platykurtotic) to increasingly positive as the mean elevation at which a species was collected historically increased (F 1,6 = 10.32, p = .018, R 2 = 63%; Figure 2b), indicating that greater canalisation in phenology was favored in higher-elevation species.

| Impacts of climate cues on synchrony of colonizers and residents with floral resources
The mismatch between peak worker density and peak flowering was predicted by one model (Table S4) with significant effects of May mean temperature (F = 11.93, p = .0024), summer precipitation (F = 5.06, p = .035), and their interaction (F = 8.68, p = .0077; Table 3; Figure 3). Respectively, each term explained 22%, 14%, and 0.9% of the variance. Overall, warmer May temperatures led to greater phenological mismatch between peak forager and resource abundances regardless of range status ( Figure 3).

| Relationship of phenological mismatch to abundances of resident and colonizing species
Differences in phenological mismatch for colonizing and resident bumble bees did not affect cumulative worker abundances and were not included in the best model (Table S6). Instead, cumulative worker abundance was best explained by a model including statistically significant one-way interactions of range status with summer maximum temperature (F = 9.72, p = .0027) and summer precipitation (F = 8.50, p = .0048; Table 4), indicating that workers of resident and colonizing bumble bee species responded differently to climate factors. Relative abundance also was best explained by a model that included statistically significant interactions of range status with summer maximum temperature (F = 18.54, p = .0001) and summer precipitation (F = 5.92, p = .018; Table 4). Resident worker abundance (both relative and absolute) decreased, whereas colonizing worker abundance increased, with summer maximum temperature, resulting in comparable densities and representation under warmer conditions (Figure 4a,c).

| DISCUSS ION
Under climate change, high-elevation communities may reach novel states due to colonization by species pre-adapted to warmer conditions, as well as loss of resident species lacking the plasticity or adaptive potential to adjust to changing regimes. Predicting outcomes of these combined processes is challenging since warming may have multiple synergistic pathways of impact on organisms.
Our results shed light on ongoing impacts of climate change for high-elevation bumble bee communities of the Rocky Mountains.
Bumble bees expanding their ranges upward from lower elevations appear to cope better with warming than high-elevation residents, although, on average, they are no better matched to phenological

TA B L E 3
Parameter estimates for the model best explaining the effects of climate variables and range status on bumble bee phenology and the phenological mismatch between bumble bee forager activity and host-plant flowering. A single model best explained the timing of each life history event (Table S4) and the phenological mismatch between day of peak worker activity and day peak flowering (Table S5). Therefore, model averaging was not required advancement in host flowering. While neither colonizing nor resident species exhibit mechanisms to synchronize foraging activity with earlier flowering under warm conditions, low-elevation bees have a more generalized strategy of seasonal resource use ( Figure 2) that may allow them to take advantage of the richer floral resources available over the season in locations with longer flowering intervals. Our findings, while based on observational rather than experimental data, concur with the prediction that phenotypic plasticity of generalist species buffers them from the disruptive effects of climate change (Herbertsson et al., 2021;Schweiger et al., 2010;Sponsler et al., 2022). Results support the view that a legacy of strong selection for survival and reproduction in a cold environment constrains the performance of resident alpine bumble bees under warming, and favors colonizing species with more generalized phenology and greater heat tolerance (Martinet et al., 2015;Oyen et al., 2016).
Unique climate drivers cue timing of key life cycle events in alpine flowering plants and their bumble bee pollinators. As reported elsewhere (Iler et al., 2013;Kudo & Cooper, 2019;Menzel et al., 2006), temperature before flowering is a strong driver of flowering phenology. Specifically, temperatures in May, weeks before the onset of the growing season at our high-elevation field sites, most strongly affect the timing of onset and peak flowering in preferred bumble bee food resources ( Figure 1a). Pronounced spatial variance in micro-climate characterizes alpine ecosystems (Kudo, 2020;Ohler et al., 2020) and this, as well as variation among host species in phenological response ( Figure 1a), likely weakened but did not overwhelm these community level effects. For every degree of warming, onset of flowering advanced by 6 days and peak flowering by 5 days (Figure 1a; Figure S3, respectively). Conversely, bumble bee queen emergence in snowy environments is primarily cued by snowmelt schedule rather than temperature (Kudo et al., 2008;Kudo & Cooper, 2019;Pawlikowski et al., 2020;Stemkovski et al., 2020). While winter snowpack and May temperature regimes are correlated (Table S6), the coupling is so weak (R 2 < .05) as to play little role in synchronizing queen and flower appearance. Weak coupling between the snowpack and early spring temperature likely reflects the importance of exposure in alpine habitats, where extreme cold temperatures characterize both wind-swept ridges and snow-packed slopes (Zwinger & Willard, 1996). Byron (1980) postulated that bumble bee emergence at snowmelt optimizes nest site selection by queens of ground nesting alpine species and that nest site availability and quality rather than food limits colony success at high elevations (see also Sarro et al., 2022).
F I G U R E 2 (a) Probability density curves for the seasonal activity of worker bumble bees in colonizing and resident species. Curves were constructed from worker density estimates for each species in each elevational zone and year; therefore, broader peaks reflect more temporal and spatial variation in phenology for each species. Density curves of colonizing species have broader peaks (more negative kurtosis) than those of resident species (W = 1.00, p = .071; Figure S7). Curves were smoothed using a bandwith adjustment of 1.25 times the standard deviation of the smoothing kernel to increase the visibility of the lines, despite overlap among species probability densities. For unadjusted density curves (see Figure S6). (b) Probability density curves exhibit progressively more kurtosis in species found at lower elevations in the past (1966-1969Macior, 1974 Worker emergence schedules were similarly insensitive to summer temperature and precipitation regimes (for temperature, see Pawlikowski et al., 2020), but consistent with historical (mid-20th century) differences in climate between low and high elevations (Table S1). Colonizers tended to emerge later than residents under all but the warmest summer temperatures (Figure 1c) and exhibited greater variation in timing of worker accrual than their alpine congeners ( Figure 2). Such life history differences between low-and highelevation species are predicted under a model of counter-gradient selection on life history traits at high elevations (Berven et al., 1979).
According to this hypothesis, historically shorter growing seasons in alpine environments favor early growth onset and rapid development for poikilothermic organisms, while the longer growing season of lower elevations favors plasticity in the growth phase (Berven et al., 1979;Laugen et al., 2003). The establishment of colonizing species from lower elevations, thus appears to expand the range of life history variation found in extant alpine bumble bee communities. Pawlikowski et al. (2020) also noted substantial interspecific variation in foraging activity schedules among European bumble bees, although underlying selective pressures promoting these differences were unclear.
While our observations cannot pinpoint traits underlying species-level differences between colonizing and resident bumble bees, research probing physiological responses to extreme temperature stress in bumble bees suggests that the two resident alpine species, B. lapponicus and B. kirbiellus are more cold tolerant and heat sensitive than their low-elevation congeners (respectively, Martinet et al., 2015;Oyen et al., 2016). Such phenotypic differences in temperature tolerance could have allowed highelevation bumble bees to begin foraging earlier under historically cooler temperatures (Figure 1c), while incurring a cost in adaptability as temperatures have warmed (e.g., an evolutionary trap, Schlaepfer et al., 2002).
Phenology describes the temporal pattern of resource and consumer abundance, not the abundances themselves. For this reason, we assessed how forager abundances in resident and colonizing species relate to climate change. Results suggest that while mismatches between the timing of peak floral resource availability and bumble bee foraging are exacerbated by warming ( Figure 3a,b), they are less likely to drive changes in bumble bee abundance or community composition than other climate effects (Tables 4 and S6; but see Ogilvie et al., 2017). This conclusion is consistent with observed shifts in traits promoting resource generalization in alpine bumble bees (e.g., shorter proboscis length; Miller-Struttmann et al., 2015). Specifically, the capacity for diet breadth to expand to less-preferred late-season flowers (e.g., Solidago, Sedum) may buffer impacts of phenological mismatches with historically favored resources (Sponsler et al., 2022), but do little to ameliorate other, more detrimental aspects of climate change.
We are aware that gaps in our knowledge could temper these findings. Our sampling schedule for bumble bee forager density precluded observation of reproductives and thus of sensitivity of recruitment to warming. While variation in worker numbers is known to affect reproductive output (Muller & Schmid-Hempel, 1992), direct observations on new queen production would strengthen inferences connecting warming to alpine species decline. Second, monitoring abundance of foragers is not equivalent to monitoring population density in eusocial organisms.
Foraging bumble bees are well known for behavioral flexibility in spatial and phylogenetic associations, potentially biasing density estimates based on forager activity (Crone & Williams, 2016).
Nonetheless, such biases should not influence estimates of rangeexpanding and range-stable species differently and thus are unlikely to account for the increasing representation of subalpine bees in high-elevation communities. Moreover, observations of colonies founded in artificial domiciles on Pennsylvania Mountain indicate that for resident alpine bumble bees the percentage of colonies producing new queens has declined from 56% to 17% over the past 40 years (Byron, 1980;Z. Miller and C. Galen, un-publ. data) consistent with negative impacts of warming in alpine ecosystems. TA B L E 4 Parameter estimates for the models best explaining the effects of climate variables and range status on cumulative and relative worker abundances. Cumulative and relative abundance were each best predicted by one model (Table S6); therefore, model averaging was not necessary

| CON CLUS IONS
As range-stable and range-expanding bumble bee species respond to climate change, novel communities are generated with discouraging implications for long-term persistence of alpine resident species.
Our study suggests that effects of climate change favor colonizing bumble bee species in historically cold environments, to the detriment of resident species. How differences in their responses to climate change affect competition between resident and colonizing species or their pollination services remains to be seen. However, prior work in this system suggests that colonizing taxa are not morphologically equivalent to resident species (Miller-Struttmann et al., 2015) and may behave differently within the community.
Future research exploring the effects of rapid changes in community composition on network structure could clarify the knock-on effects of displacement for community stability and pollination.

ACK N OWLED G M ENTS
We thank the late L. W. Macior and P. A. Byron, along with our colleague and mentor P. G. Kevan, for pioneering research on Rocky Mountain bumble bees that provided the foundation for this work.