www.sciencemag.org/content/349/6255/1541/suppl/DC1

Supplementary Material for Functional mismatch in a bumble bee pollination mutualism under climate change Nicole E. Miller-Struttmann,* Jennifer C. Geib, James D. Franklin, Peter G. Kevan, Ricardo M. Holdo, Diane Ebert-May, Austin M. Lynn, Jessica A. Kettenbach, Elizabeth Hedrick, Candace Galen *Corresponding author. E-mail: [email protected] Published 25 September 2015, Science 349, 1541 (2015) DOI: 10.1126/science.aab0868 This PDF file includes: Materials and Methods Figs. S1 to S4 Tables S1 to S9

2

Materials and Methods Change in bumble bee tongue length and impact of covariance with body size Using archived and field-collected worker bumble bees, we tested for temporal changes in the tongue length and intertegular space (the distance between wing bases) of two alpine resident species Bombus balteatus and B. sylvicola. We used intertegular space to measure body size, since it is highly correlated with body mass (32) and is not sensitive to pollen load or specimen storage time. The three alpine field sites surveyed in this study correspond to locations of historic collections of B. balteatus made by L.W. Macior (18) on Mount Evans in Clear Creek County (Mount Evans Wilderness Area, Arapaho National Forest, 39°35.033'N, 105°38.307'W) and on Niwot Ridge in Boulder County (Niwot Ridge Long Term Ecological Research Site, 40°3.567'N, 105°37.000'W), and by P. A. Byron (17), R. C. Plowright and others on Pennsylvania Mountain in Park County (39°15.803'N, 106°8.564'W). Collections of B. sylvicola were made concurrently at two of these sites (Niwot Ridge and Pennsylvania Mountain). From 1966-1969, Macior (18) surveyed and identified more than 8,000 bees from alpine habitats (>3500m) at Mount Evans and Niwot Ridge. Bumble bees collected by Macior from 1966-1969 were deposited in the University of Colorado Entomology Museum. While only the county of origin was listed for each specimen and exact locations were not provided, 97.5% and 97.8% of individuals of B. balteatus and B. sylvicola, respectively, were collected in alpine habitats (18). In 2012-2014, we resurveyed bumble bee communities at four altitudes on Niwot Ridge (3500m, 3600m, 3675m and 3750m) and five altitudes on Mount Evans (3500m, 3600m, 3700m, 3800m and 3900m) in accordance with Macior (18). We collected, identified and determined the caste of all bumble bees encountered during weekly, one hour surveys along altitudinal transects at each site. On Pennsylvania Mountain, Byron (17) surveyed and identified bumble bee workers weekly during the flowering season along transects at two altitudes (3505m in willow krummholz habitat and 3960m on the summit) from 1977-1979. In addition to Byron’s specimens, collections were made as part of a UNESCO Man-and-Biosphere project between 1977-1980. In 2008, 2011 and 2013, we collected worker bumble bees along transects at 3659m and 3730m, near locations of Byron’s original surveys. While these transects did not overlap completely with Byron’s collection sites (17), average distance separating current and historic sites (<1km) is within the flight distance of both bumble bee species (19). Those individuals collected in 2013 on Pennsylvania Mountain were compared with historic specimens. To match the timing of historical collections at each field site and minimize possible negative demographic effects on colony production, lethal collections were conducted later in the flowering season (mid-July to August) when foraging workers were at peak abundance. The glossa of each archived individual collected from past populations on Pennsylvania Mountain (1977-1980) and present populations at all sites (2008, 2011-2014) was removed and measured. Historical specimens collected from Niwot Ridge and Mount Evans (1966-1969) could not be dissected. Therefore, we rehydrated specimens for 24 hours, extended the mouthparts, and measured the glossa. To correct for any potential bias in measurement method, we measured specimens collected from Niwot Ridge and Mount Evans in 2014 using both methods, and analyzed the relationship between them

3 via linear regression. We then used this relationship (y = 0.18 + 0.974x; R2 = 0.731, t1,61 = 12.87, P < 0.0001) to correct the tongue length of historic specimens from these sites. To match collection effort (number of specimens) in the past, we randomly sub-sampled individuals from the present, larger collections. This and all other analyses were conducted in the R Statistical Software (33), unless otherwise noted. We tested for a temporal change in tongue length of each species via analysis of covariance (ANCOVA) with location and time (i.e., past and present) as fixed effects (lm in the stats package (33) and Anova in the car package (34)). Because tongue length is correlated with body size, we used intertegular space as a covariate. The proportion of variance explained by each independent variable was calculated as partial eta squared (35; etaSquared in the lsr package (36)). Similarly, we conducted a two-way analysis of variance (ANOVA) to test for a change in body size (intertegular space) over time with location and time as fixed effects. Using bees collected during 2012-2014 at Mount Evans and Niwot Ridge, we tested for seasonal and interannual changes in tongue length. If bees vary in tongue length over the season, independent of body size, then use of intertegular space as a covariate would not fully account for temporal plasticity and the comparison between past and present specimens could be sensitive to date of collection. We tested for unmeasured sources of seasonal plasticity in tongue length via ANCOVA with collection day of year as a continuous effect, location as a categorical fixed effect and body size as a covariate (lm in the stats package (33)) (fig. S4). Similarly, we tested for inter-annual plasticity via ANCOVA with year and location as categorical variables and intertegular space as a covariate (fig. S4). A significant day of year effect or year effect would indicate the unmeasured sources of temporal plasticity. Changes in flower tube depth as a potential driver of tongue length evolution Using archived plant specimens and two historical datasets, we tested the prediction that directional changes in flower depth have driven the evolutionary decreases in tongue length of alpine bumble bees. Changes in flower tube depth within host plant species were measured using published records and archived specimens collected from or near Niwot Ridge and Mt Evans. Macior (18) reported corolla tube depth for eleven plant species, including many that are important bumble bee food resources, based on collections made from 1966-1969 on Mount Evans and Niwot Ridge. In 2012 and 2013, we collected, pressed and dried flowering plants for six of these species (Castilleja occidentalis, Mertensia lanceolata (synonym M. viridis (37)), Polemonium viscosum, Trifolium dasyphyllum, T. nanum and T. parryi) from the same localities. We then rehydrated the flowers to make replicate measurements of the same floral traits (18). Tube depth was measured as follows: for C. occidentalis as the distance from the base of the flower to longest point of fusion in the corolla; for M. lanceolata as the distance from the base of the corolla to point where the lobes separate from each other; and for Trifolium spp. as the length of the longest petal. To ensure that past and present measurements were comparable, we also measured pressed flowers of each species from herbarium specimens collected at Mount Evans, Niwot Ridge or within a 15 km radius of each site between 1960 and 1982. Measurements of herbarium specimens matched published records for all but one plant species (T. parryi; table S7). For this species, we compared the modern specimens to herbarium specimens instead of Macior’s (18) records. To test the prediction that flowers are shallower now than historically, we

4 conducted a multivariate analysis of variance (MANOVA) with tube depth as the response variable (manova in the stats package (33)). If significant temporal change in flower tube depth was supported, univariate analyses (ANOVA) were conducted for each species (summary.aov in the stats package (33)). To examine whether declines in abundance of long-tubed floral resources favor bumble bees with shorter tongues, we surveyed bumble bee food plants that exhibit interspecific variation in flower tube depth (18) at sites where density of flowers or relative abundance of host plants was recorded historically: Pennsylvania Mountain and Niwot Ridge. Five taxa were monitored at both locations (C. occidentalis, M. lanceolata, T. dasyphyllum, T. nanum, and T. parryi). One additional species, Polemonium viscosum, was monitored at Pennsylvania Mountain, and, at Niwot Ridge, six additional species were sampled (Chionophila jamesii, Pedicularis sudetica, Phlox condensata, Primula parryi, Rhodiola rhodantha and Silene acaulis). Survey methods varied between locations to match historical data. We tested for a relationship between species mean flower depth and percentage change in relative abundance over time based on vegetative cover (Niwot Ridge) or flower density at peak bloom (PFD; Pennsylvania Mountain) via linear regression (lm in the stats package (33)). In 1971, author Ebert-May established thirty 1 × 10 m plots across the range of vegetation zones occurring at 3528 m on Niwot Ridge (38). Percent vegetative cover of each plant species was estimated in 0.10 × 1 m strips at 1m intervals in each plot and resampled in 1991, 2001, and 2011. Species’ vegetative cover in 1971 and 2011 was summed across all plots to estimate total vegetative cover in years that correspond with timing of bumble bee surveys at Niwot Ridge (18). Flower production was not measured; we use vegetative cover as an indication of total plant production. Vegetative cover predicts flower production in many alpine plants (39); however, changes in vegetative cover should be interpreted with caution in terms of implications for flower production. On Pennsylvania Mountain, author P.G. Kevan surveyed flower abundances weekly from 1977-1980 in 23, 2 × 10 m plots spanning five habitat types from 3600-3969 m (17). Using aerial photographs and maps, we relocated historic plots. We ensured that each relocated plot included the historically surveyed area by adding a buffer strip around the re-located “edge” of the plot. Including this buffer increased the original 2 m dimension to 10 m and expanded the 10 m edge to a range of 24-34 m. To account for the larger sampling scale in re-located plots, statistical comparison of flower density over time was based on the average flower density in all possible independent 2 × 10 m subplots within the larger plots surveyed in 2012-2013. Open flowers or inflorescences were counted weekly from the onset to the end of the flowering season. As inflorescence structure varies among plant species, we used the average number of open flowers for a given species m-2 at peak flowering (peak flower density, PFD) as a common metric of floral resource production. When necessary, field counts of inflorescences were converted to flowers using counts from a sub-sample of individuals in each plot (average of 30 plants per species per plot, 2012-2013); published counts from Pennsylvania Mountain (40; T. dasyphyllum and T. parryi) and for one species (C. occidentalis) flower counts from Niwot Ridge and Mount Evans (2012-2013).

5 Competition from immigrating subalpine bumble bees We resampled historically surveyed bumble bee communities (17, 18) to address the prediction that competition from immigrating subalpine bumble bees has driven tongue length evolution in resident alpine congeners. Using all bumble bee species sampled at our field sites during historic and modern surveys, we tested for a change in community composition via χ2 test for independence (chisq.test in the stats package (33)). Abbreviations for bumble bee species are as follows: B. hun, B. huntii; B. ruf, B. rufocinctus; B.mix, B. mixtus; B.occ, B. occidentalis; B.fri, B. frigidus; B.bif, B. bifarius; B.sylv, B. sylvicola; B.mel, B. melanopygus; B.cen, B. centralis; B.fla, B. flavifrons; B.bal, B. balteatus; B.nev, B. nevadensis; B.app, B. appositus; and B.mor, B. morrisoni. Because Macior’s original field notebooks have been lost and (unlike his archived specimens) his published records (18) do not indicate exact location or date of collection, bumble bees observed on Niwot Ridge and Mount Evans (hereafter referred to as the Front Range) were pooled across sites and years (16). Changes in the interspecific distribution of tongue length among alpine bumble bees To address how changes in bumble bee species composition over time have altered the distribution of interspecific variation in tongue length at our alpine field sites, we used past surveys of bee visitation (Front Range: 1966-1969 (18); Pennsylvania Mountain: 1977-1979 (17)) and present collections of foragers from the same sites (Front Range: 2012-14 and Pennsylvania Mountain: 2008, 2011 and 2013) to sample relative abundances of bumble bee species. We focused solely on workers to avoid the larger demographic impact of removing queens. For past surveys of Front Range bees, we calculated the proportion of individuals in each species that were workers using queen to worker visitation ratios on alpine plants as reported in the source data (18). We simulated 1000 communities for each time period and assigned tongue length to individuals based on the historical mean and standard deviation of their species from the Front Range, excluding two species (B. huntii and B. rufocinctus) which were not measured (18). These species were rare, accounting for 0% and 1.6% of individuals collected in the past and present, respectively. The population size, N, of each simulated community was set to the number of individuals collected in present surveys, as these N are lower than the historical values. Next, we characterized the multi-modal distribution of each simulated community using mixture models (normalmixEM in the mixtools package (41)), which describe the distribution as a combination of two subpopulations (e.g., subpopulation 1, comprised of short-tongued bumble bees and subpopulation 2, of longer-tongued bumble bees), each with a normal distribution. Preliminary analysis of the tongue length distribution in communities of past and present bumble bees indicated bimodal distributions with peaks at around 6 mm (subpopulation 1) and 9 mm (subpopulation 2). Using these estimates as starting conditions for mixture models, we characterized the tongue length distribution for each simulated community and recorded the mean and variance in tongue length, and proportion of individuals in each peak (phenotypic subpopulation). We tested for a temporal change in density function after excluding any simulations that did not resolve to two peaks (less than 1%) via MANOVA (manova in the stats package (33)). If the multivariate test indicated a significant change in the shape of the density function, temporal changes in the mean and variance for subpopulations 1 and 2 and the proportion

6 of individuals in subpopulation 1 were tested via univariate ANOVA (summary.aov in the stats package (33)). Shift in realized foraging breadth We compared resource use (host species richness) of past and present bumble bee communities in the Front Range. While we resampled the majority of sites originally surveyed by Macior (for sampling details see “Change in bumble bee tongue length and impact of covariance with body size”) in 2012 and 2013, we did not resample his three highest-altitude sites on Mount Evans, because high tourist activity at these sites could disrupt bumble bee populations. Because bumble bee visitation rates could not be separated by altitude (unlike community composition) (16), additional observations at 4100 m and 4200 m were conducted in 2014 to ensure complete coverage of the altitudinal range sampled by Macior (18). To compare past and present foraging breadth (Levin’s niche breadth: B = 1/∑pj2, where pj is the proportion of individuals visiting plant species j; 42) of alpine bumble bees, we calculated Z scores. Resampling (n = 1000 iterations) of historic bumble bee visitation records was constrained to the sample size of present collections to account for differences in sampling effort (n = 4099 and 519 visits respectively). We tested for a change in the distribution of flower tube depth for plants visited by B. balteatus and B. sylvicola. Competition from subalpine immigrants could lead to visitation of a narrower range of flower tube phenotypes by alpine resident bumble bees. We simulated the distributions of flower tube depth for plants visited by each alpine bumble bee species in the past (1966-1969; 9 and 10 plant species, respectively) and present (2012-2014; 13 and 18 plant species, respectively) (n = 1000 iterations), with bumble bee population size set to the number of individuals collected in current surveys. For each flower visited, tube depth was simulated using the species mean and standard deviation from published records or rehydrated, pressed specimens (table S8). Tube depth of 11 species (table S8) follows (18) and one species (Primula parryii) follows (43). Tube depth was measured as follows: for all species in the Asteraceae (i.e., Cirsium scopularum, Erigeron simplex, Packera cana, Solidago simplex, Tetraneuris grandiflora, Tonestus pygmaeus) and Gentianella amarelle, as the distance from the base of the corolla to the crease where the fused petals part; for Castilleja puberula, as the distance from the base of the flower to longest point of fusion in the corolla; and for Eriogonum flavum as the distance from the base of the corolla to the end of the fused part of the petals. Several plant species (i.e., Arenaria fendleri, Rhodiola rhodantha, Sedum lanceolatum) have open flowers with exposed nectaries rather than flower tubes. We excluded visits to these species from temporal comparison of flower tube depth but note their frequency at each time interval (table S3). Mean and variance in tube depth within the assemblage of host plants visited were calculated for resampled distributions. We tested for temporal change in tube depth phenotype (mean and variance) via MANOVA (manova in the stats package (33)). Following a significant temporal change in phenotype distribution, changes in mean and variance were tested separately via univariate ANOVA (summary.aov in the stats package (33)). Declining resource availability under climate change

7 We examined the response of foraging bumble bees (B. balteatus) to changes in flower density in the field to test for an effect of resource scarcity on foraging niche breadth. We then compared present to past flower abundance on Pennsylvania Mountain where detailed records of flowering were made from 1977-1980 for six bumble bee-pollinated species in habitats from 3600-4000 m (17) to determine if flower densities have declined. Last, we asked whether observed expansion of foraging breadth in response to resource scarcity is predicted from a mathematical model of pollinator foraging based on optimal foraging theory (26). Response of bumble bee foraging breadth to changes in flower density Bombus balteatus queens and workers were followed while foraging in the field to test whether they are more likely to incorporate new resources into their diet when preferred host plants are scarce. On Pennsylvania Mountain (3640 m), we monitored movement of foraging bumble bees from plant to plant noting species visited, distance flown, and distance to the four nearest conspecific neighbors of each visited plant. Average nearest neighbor distance was converted to flower density (44) using species-specific measurements of flowers per inflorescence as described above (“Changes in flower tube depth as drivers of tongue length evolution”). We tested whether flower density of the first species visited in sequence predicted the change in floral tube depth (Δt) from one visit to the next. For intraspecific moves, Δt = 0, while on moves to shorter-tubed species, Δt >0. If bumble bees are more likely to incorporate shallower flowers in their diet as resources become rare, then Δt should increase in sparse patches. We used linear regression of Δt on flower density to test this hypothesis (SAS 9.2 PROC REG (45)). Flower density was log-transformed prior to the analysis to correct for departures from normality. We tested for temporal change in summer low temperatures at Mount Evans and Pennsylvania Mountain using PRISM interpolation (46). We examined the relationship between summer low temperature and flower density on Pennsylvania Mountain using four bumble bee host plant species that together accounted for 57% of flowers visited by B. sylvicola and B. balteatus historically (17). As a standardized metric for flower density of species with different inflorescence architecture, we used mean flower density at peak bloom (flowers m-2, PFD). We calculated mean PFD per species per plot (n = 23) along the altitudinal gradient (from timberline to the summit) in each of six years between 1977 and 2014. These mean PFD values were subjected to quadratic regression on mean minimum summer temperatures (from June 1-Aug 31; estimated using PRISM interpolation (46)). We log-transformed PFD to meet assumptions of normality of parametric statistics. To account for potential interspecific variation in flowering response, we included species as a random effect in the regression model (ANCOVA; SAS 9.2 PROC Mixed (45)). Differences among species were non-significant (Z = 0.91, P > 0.183). The equation for the best fit line was given by: ln (PFD) = 2.71(temp) – 0.425(temp2) – 2.29. The linear term of the model (F1,y = 3.44, P = 0.081) suggested a trend for increasing flower production up to a threshold of 3.250C while the negative quadratic term (F1,17 = 4.73, P < 0.044) indicated decreasing flower production above this point.

8 To ask whether floral resources have declined at the landscape scale in response to warming over the time span of this study, we used a digital elevation model from the U.S. Geological Survey (47) in ArcGIS (48) to define polygons referenced to habitats below timberline (subalpine, < 3600m; excluded), uninhabitable areas (cliffs and steep talus slopes greater than 30o in incline; mining disturbances; also excluded), and each of five alpine habitat types: krummholz, slope, swale, false summit and summit. Where altitudinal span of habitats overlapped (e.g., swale and false summit), boundaries were mapped using waypoints obtained with high resolution GPS (SX-Blue II). We compared past (1977-1980) and present (2012-13) flower densities of six bumble bee host plants (which in combination accounted for 88% of historical bumble bee diet) on Pennsylvania Mountain in similarly snowy (1979-1980 and 2013) and dry (1977 and 2012) years. At the plot scale, the change in PFD was analyzed using a mixed model ANOVA. We treated time (past vs. present), habitat, and precipitation (high vs. low winter snowpack) as fixed effects, and plot and plant species as random effects (ANOVA, SAS 9.2 PROC Mixed (45)). For sampling details see above, “Changes in flower tube depth as drivers of tongue length evolution”. We combined measurements of PFD in plots of each habitat with habitat area (table S5) to compare the loss in flower production in low altitude habitats to the gain in flower production at the summit. For a bumble bee foraging from a habitat set containing habitats i and j, the contribution of flowers in each habitat to the bee’s resource intake is defined as the product of flower density (PFD) in that habitat and total surface area (s) of that habitat. Assuming that all available surface area supports plants then, PFDi × si = flower production in i PFDj × sj = flower production in j PFDtot = PFDi + PFDj If the two habitats are of equal size, a decline in PFD in one habitat would simply be balanced by an equivalent increase in PFD in the other. PFD decreased by 2.61 flowers m-2 in the krummholz and 3.51 flowers m-2 on the slopes between 1977-80 and 2012-13. On the summit, PFD increased over the same time span by 3.0 flowers m-2. However, together the krummholz and slope habitats are more than an order of magnitude larger than the summit (table S5), resulting in a net 60% decline in flowers over time at the landscape scale. Modeling energetic returns of pollinator foraging breadth under variable flower density We modified a published mathematical model (26) based on optimal-foraging theory (49) to predict the impact of declining floral density on foraging efficiency (rewards gained per unit time) of different resource choice strategies. The model describes a simple, twospecies system where a generalist pollinator visits all available resources (i.e., both deep and shallow flowers), whereas a specialist pollinator forages only on deep flowers. Because short-tongued bees exhibit greater foraging generalization than long-tongue bees (16, 30), we assumed that generalist short-tongued bees will visit both deep and shallow flowers, whereas specialist long-tongued bees will only visit deep flowers.

9 In this model, total foraging time per flower is defined as the sum of the time required to travel to each flower (tf, hitherto referred to as flight time) and handling time per flower (th). If the time required to travel to each flower increases linearly with the distance between flowers, then flight time (tf,g) of the generalist species is given by

t f ,g =

1

σ δ

where σ is flight speed and δ is flower density. For a specialist visiting only deep flowers, flight time is dependent on the proportion of available resources that are deep flowers. Specifically,

t f ,s =

1

σ ρδ

where ρ is the proportion of flowers with deep tubes. Assuming that handling time varies with flower tube depth (10, 50, 51), then total foraging time (ttot,g) for a generalist is given by

 1  + ρτ g ,d + (1 − ρ )τ g ,s  ttot ,g = N  σ δ  where N is the number of flowers visited, and τg,d and τg,s are the handling times for a generalist feeding on a deep flower and shallow flower, respectively. The generalist visits both flower types, and deep flowers are encountered with probability ρ. Total foraging time of a specialist is given by  1  + τ s ,d  ttot ,s = N   σ ρδ 

where τs,d is the handling time for a specialist feeding on a deep flower. The relative advantage (RA) of being a generalist is defined as the ratio of total foraging time of the specialist to that of the generalist. Specifically,

1 t RA = tot ,s = ttot ,g

σ ρδ 1

σ δ

+ τ s ,d

+ ρτ g ,d + (1 − ρ )τ g ,s

When the total foraging time of the specialist is greater than that of the generalist (RA > 1), then the generalist has an energetic advantage and is favored. Conversely when the total foraging time of the specialist is less than that of the generalist (RA < 1), the specialist has the energetic advantage. We parameterized the model (table S9) with past and present resource abundances (PFD) on Pennsylvania Mountain (Fig 3B; see above, “Response of bumble bee foraging breadth to changes in flower density” for methods) to explore how the relative advantage of being a generalist has changed over time. For search speed of bumble bee foragers we

10 followed (52-54); and for handling times of long and short-tongued bumble bees we used data from (50). This model makes several simplifying assumptions that warrant further discussion. First, handling of flowers and travel between them are assumed to have equal energetic costs. However, the energetic cost of flight is likely greater than handling of flowers (55). Model outcomes for effects of flower density on the generalist advantage are conservative in this regard. When energetic costs of handling flowers and flight are equal, the specialist foraging strategy is more likely to be favored, since the advantage of being a specialist is driven by increased handling time efficiency (table S9). Second, the relative advantage of being a generalist is based only on handling and flight times, and all flowers are assumed to be equally rewarding. Instead, plants often vary in nectar rewards (e.g., 56, 57), which may shift the conditions under which the specialist or generalist foraging strategy is favored. Exploring the impacts of reward variation on bumble bee behavior is beyond the scope of this study. Nonetheless, the simplified model presented here provides broad predictions concerning the conditions under which generalization is favored. Finally, the model assumes that pollinators visit all flowers they encounter. While this is not always the case (e.g., 58), relaxing this assumption only modestly altered the effects of density on visitation patterns in the original model (26).

11

Fig. S1 Temporal change in body size and the relationship between body size and tongue length. Intertegular space (mean +/- S.E.) is a measure of body size for bumble bee workers (table S1). Bombus balteatus (A) from Mount Evans, Niwot Ridge and Pennsylvania Mountain (Penn. Mtn.), and B. sylvicola (B) from Niwot Ridge and Pennsylvania Mountain (Penn. Mtn.). Means of groups denoted by asterisks differ significantly at P < 0.05. Means denoted by ns do not differ at P < 0.05. Tongue length is correlated with body size for both species (C and D; respectively; y = 2.35 + 0.965x, Pearson’s r = 0.50, t1,28 = 3.04, P = 0.005; y = e0.213 + 0.341x, Pearson’s r = 0.60, t1,70 = 6.24, P < 0.0001; shaded symbols indicate past specimens and open symbols indicate present ones) (cor.test in the stats package (33)). However, body size alone is insufficient to explain declines in tongue length (table S1).

12

Fig. S2 Shifts in flower tube depth (A; mean +/- S.E.) for bumble bee host plants in the Front Range and the relationship of tube depth to the temporal change in abundances of flowers (B) or plants (C). Changes in tube depth vary in direction and magnitude among species despite a significant effect of time across all species (A, MANOVA, F6,13 = 9.42, P = 0.0004). Declines in tube depth are not evident for plant species visited most heavily by bumblebees (e.g., TRNA and TRPA). Similarly, the change in (B) mean flower density at peak bloom (PFD) on Pennsylvania Mountain (1977-80 vs. 2012-13) and (C) vegetative cover at Niwot Ridge (1971 vs. 2011) is not related to flower tube depth (respectively; R2 = 0.227, t1,4 = 1.21, P = 0.294; R2 = 0.0004, t1,9 = -0.062, P = 0.952). Asterisks denote significant differences in mean tube depth over time (* P < 0.05, ** P < 0.001, *** P < 0.0001). Four letter codes denote plant species: CAOC, Castilleja occidentalis; MELA, Mertensia lanceolata; POVI, Polemonium viscosum; TRDA, Trifolium dasyphyllum; TRNA, T. nanum; and TRPA, T. parryi.

13

Fig. S3 Recent warming at Mount Evans and Pennsylvania Mountain and the relationship between mean minimum summer temperature and flower density (PFD) on Pennsylvania Mountain. Average minimum summer temperatures (AMST) have increased ~ 2°C since the mid-1960s at Mount Evans (A; R2 = 0.383, t1,52 = 5.68, P < 0.0001) and Pennsylvania Mountain (B; R2 = 0.341, t1,52 = 5.20, P < 0.0001; red symbols denote years with AMST above 3.25⁰C, predicting reduced PFD; blue symbols denote years with AMST below 3.25⁰C) as estimated using PRISM interpolation (46). The relationship of AMST to PFD (flower density at peak bloom, flrs m-2) averaged across 23 sampling plots on Pennsylvania Mountain for four bumble bee host species sampled in six years between 1977-2014 (quadratic R2 = 0.23, t1,17 = -2.18, P = 0.041) shows PFD declining above 3.25⁰C (C). Unshaded symbols represent recent observations (2012-14) and shaded symbols denote past observations (1977,78-80).

14

Fig. S4 Intra-annual and inter-annual variation in tongue length for Bombus balteatus and B. sylvicola collected in the Front Range (Mount Evans and Niwot Ridge). Bars represent least squares mean estimates that account for variation in body size, and error bars denote standard errors. After removing variation due to body size, tongue length did not vary significantly for either species with Day of Year (A and B; respectively; F1,18 = 1.43, P = 0.25; F1,38 = 0.70, P = 0.41) or year (C and D; respectively; F2,16 = 0.77, P = 0.48; F2,36 = 0.64, P = 0.53). Tongue length also did not vary among locations for either species (respectively; F1,18 = 0.20, P = 0.66; F1,38 = 0.029 P = 0.87; data not shown).

15 Table S1. ANCOVA and ANOVA for temporal changes in tongue length and intertegular space (respectively) of two alpine bumble bee species (Bombus balteatus and B. sylvicola). Location (Mount Evans, Niwot Ridge and Pennsylvania Mountain for B. balteatus, and Niwot Ridge and Pennsylvania Mountain for B. sylvicola) and time (past, 1966-80 vs. present, 2012-14) were fixed categorical effects. Intertegular space was included as a covariate in the analysis of tongue length, to account for the potential relationship between tongue length and body size. The partial eta squared (η2p; 35) describes the variance explained by each independent variable. Tongue length of B. sylvicola was log transformed to meet the assumption of normality. Results from Type II ANOVA are reported when the interactions are not significant (tongue length for B. balteatus and intertegular space for B. sylvicola) (59). Main Effect

B. balteatus P

df

F

Tongue length Intertegular space Location Time Location × Time Error

1 2 1 2 23

5.06 0.35 17.02 1.51

Intertegular space Location Time Location × Time Error

2 1 2 96

1.34 1.02 8.61

B. sylvicola P η2p

η2p

df

F

0.034 0.71 0.0004 0.24

0.18 0.03 0.43 0.12

1 1 1 2 67

5.80 10.71 46.14 8.13

0.019 0.0017 0.0001 0.0058

0.08 0.14 0.41 0.11

0.26 0.31 0.0004

0.03 0.01 0.15

1 1 2 76

75.73 0.0001 29.01 0.0001 1.97 0.17

0.37 0.21 0.12

16 Table S2. Temporal change in parameters of the tongue length density function for bumble bee communities in the Front Range (individuals from Niwot Ridge and Mount Evans pooled; 1966-69 vs. 2012-14) and on Pennsylvania Mountain (1977-80 vs. 2008, 2011, 2013). Results for each parameter are from univariate tests following a significant effect of time across all parameters as tested via MANOVA (Pillai approximate F1,1994 = 22,860, P < 0.0001; Pillai approximate F1,1984 = 16,635, P < 0.0001; for the Front Range and Pennsylvania Mountain, respectively). The proportion of individuals in subpopulation 2 (belonging to long-tongued bees) was not included in the MANOVA, because it is not independent of the proportion in population 1. Asterisks denote significant differences between group means (* P < 0.05, ** P < 0.001, *** P < 0.001). Past Mean S.D.

Present Mean S.D.

Univariate ANOVA df F

0.568 5.86 0.334

1 1 1

94,618*** 797.0*** 1477.5***

0.295 0.020 9.30 0.097 0.706 0.153

--1 1 1998

--303.9*** 3667.4***

0.017 0.049 0.042

0.377 0.032 5.75 0.074 0.256 0.083

1 1 1

85.6*** 165.5*** 700.29***

0.017 0.041 0.365

0.623 0.032 8.54 0.108 1.15 0.158

1 1 1 1988

--50,092*** 22,209***

Front Range (Mount Evans and Niwot Ridge) Short-tongue subpopulation (Peak 1) Proportion of individuals 0.432 0.020 Mean tongue length (mm) 5.81 0.041 Variance in tongue length 0.399 0.039 Long-tongue subpopulation (Peak 2) Proportion of individuals 0.568 0.020 Mean tongue length (mm) 9.36 0.039 Variance in tongue length 0.405 0.038 Error Pennsylvania Mountain Short-tongue subpopulation (Peak 1) Proportion of individuals 0.366 Mean tongue length (mm) 5.79 Variance in tongue length 0.334 Long-tongue subpopulation (Peak 2) Proportion of individuals 0.634 Mean tongue length (mm) 9.36 Variance in tongue length 0.384 Error

0.020 0.039 0.037

Table S3. Temporal change in tube depth phenotype of flowers visited by Bombus balteatus and B. sylvicola in the Front Range (individuals from Niwot Ridge and Mount Evans pooled, 16; 1966-69 vs. 2012-14). Results for each parameter are from univariate tests following a significant effect of time as tested via MANOVA (Pillai approximate F1,1997 = 7554, P < 0.0001; Pillai approximate F1,1997 = 64,851, P < 0.0001; for B. balteatus and B. sylvicola, respectively). Plants with exposed nectaries rather than flower tubes were excluded from the analysis. The proportion of visits to these tubeless flowers increased over time from 0.37% to 0.63% for B. balteatus and 9.9% to 24.7% for B. sylvicola. Asterisks denote significant differences between groups (* P < 0.05, ** P < 0.001, *** P < 0.001). Past Mean S.D.

Present Mean S.D.

Univariate ANOVA df F

B. balteatus Mean flower tube depth (mm) Variance in flower tube depth Error

14.14 2.36

0.11 1.04

13.86 0.22 6.84 1.50

1 1 1998

1029*** 15,072***

B. sylvicola Mean flower tube depth (mm) Variance in flower tube depth Error

12.62 20.65

0.31 3.16

7.99 0.36 28.58 1.86

1 1 1998

129,569*** 17,213***

Table S4. Regression to test the effect of flower density on the difference in average flower tube depth (Δt) between sequentially visited plants chosen by free-foraging Bombus balteatus. Density refers to that of the species visited initially and was log-transformed prior to regression to correct for deviations from normality. Bees were observed visiting tundra plants of Polemonium viscosum, Trifolium parryi and Mertensia lanceolata on Pennsylvania Mountain. For the best fit line: Δt = 6.95 – 1.87 log flower density, R2 = 0.127. Effect Flower density Error

df 1 194

F 29.39

P 0.0001

Table S5. Surface area and altitude for each habitat on Pennsylvania Mountain for which flowering abundance was recorded in the past (1977-1980) and present (2012-2013). Surface area, altitudinal midpoint and altitudinal range were estimated for each habitat polygon created using a digital elevation model from the U.S. Geological Survey (47) and waypoints to delimit habitat boundaries in ArcGIS (48). Habitat Krummholz Slope Swale False summit Summit

Surface area (km2) 0.60 1.35 0.41 0.88 0.10

Altitude (m) Range Midpoint 3600-3650 3625 3650-3800 3725 3757-3938 3847 3800-3938 3869 3938- 3970 3954

Table S6. ANOVA for change in mean flower density at peak bloom (PFD) of six plant species that accounted in the past for 88% of flower visitation by bumble bees on Pennsylvania Mountain (17). Time (past, 1977-80 vs. present, 2012-14), habitat, and precipitation (winter snowpack high: 69-96 percentile vs. low: < 22 percentile over years from 18952013, based on cumulative January-May precipitation) were fixed categorical effects, and plant species and plot were random effects. PFD was log-transformed prior to the analysis to correct for departure from normality. Effect Time Habitat Winter snowpack Time × habitat Time × winter snowpack Habitat × winter snowpack Time × habitat × winter snowpack Error

df 1 4 1 4 1 4 4 385

F 2.38 1.21 1.97 5.55 12.13 0.80 0.28

P 0.12 0.34 0.16 0.0002 0.0006 0.53 0.89

Table S7. Flower tube depth for past populations of bumble bee host plants in the Front Range of Colorado as measured by Macior (1966-69) and as re-measured using herbarium specimens from the region (1960-80). Macior reported only means and standard deviations (18); therefore, we sampled from a normal distribution with the mean and standard deviation given in the published record to estimate tube depth. We then compared measurements from archived specimens to the resampled tube depth measurements. Asterisks denote significant differences between group means (* P < 0.05). Host Plant Species

Tube depth from published record (mm) n Mean S.D. Castilleja occidentalis 25 16.2 1.5 Mertensia lanceolata 25 7.8 0.4 Polemonium viscosum 25 21.8 0.8 Trifolium dasyphyllum 25 11.5 0.6 T. nanum 25 14.5 0.8 T. parryi 25 13.9 0.6

Resampled tube depth from published record (mm) n Mean S.D. 2 15.8 1.3 20 8.0 0.5 8 21.8 0.6 6 11.6 0.4 11 14.2 0.7 10 14.1 1.0

Tube depth from herbarium specimens (mm) n Mean S.D. 2 15.9 1.7 20 7.6 1.8 8 21.5 1.9 6 12.4 1.2 11 15.0 1.8 10 15.4 1.3

Resampled vs. herbarium specimens df F 1 0.010 1 2.06 1 0.66 1 1.99 1 2.16 1 7.04*

22 Table S8. Mean and standard deviation of flower tube depth for plants visited by bumble bees on Niwot Ridge and Mount Evans. Asterisks denote species for which tube depth was measured from rehydrated, pressed specimens as described in the text (“Shift in realized foraging breadth”). For Cirsium scopulorum, field-collected specimens were augmented with archived specimens from Niwot Ridge and Mount Evans. Tube depth of Primula parryi (denoted by a dagger) on Niwot Ridge follows (43). All other measurements are from (18). Host Plant Species Polemonium viscosum Chionophila jamesii Castilleja occidentalis Penstemon whippleanus Trifolium nanum Trifolium parryi Erysimum nivale Castilleja puberula* Cirsium scopulorum* Primula parryi† Pedicularis sudetica Trifolium dasyphyllum Mertensia lanceolata Gentianella amarella* Silene acaulis Tonestus pygmaeus* Packera cana* Tetraneuris grandiflora* Solidago simplex* Erigeron simplex* Eriogonum flavum*

Tube depth (mm) n Mean S.D. 25 21.8 0.8 25 16.6 0.6 25 16.2 1.5 25 15.7 0.7 25 14.5 0.8 25 13.9 0.6 25 13.5 0.6 15 13.1 2.4 4 13 0.5 50 12 2.0 25 11.5 0.5 25 11.5 0.6 25 7.8 0.4 3 6.6 2.9 25 6.6 1.0 18 5.7 0.6 14 5.2 0.7 7 4.9 0.3 19 3.4 0.6 23 2.7 0.3 4 1.8 0.6

23 Table S9. Parameters used to model relative efficiency of generalist and specialist foraging strategies. Range in flower density was based on records from Pennsylvania Mountain (17). In flight search speeds are taken from (52), (53), and (54); and handling time from (50). Parameter Δ Σ Ρ τs,d τg,d τg,s

Description Flower density Search speed in flight Proportion of flowers with deep tubes Specialist handling time on deep flowers Generalist handling time on deep flowers Generalist handling time on shallow flowers

Value 0.1-10 0.2, 0.5, 1, 2, 5 0.1-1 2.04 3.32 2.08

Units flrs/m2 m/s -s/flr s/flr s/flr