Specialization of plant–pollinator interactions increases with temperature at Mt. Kilimanjaro

Abstract Aim Species differ in their degree of specialization when interacting with other species, with significant consequences for the function and robustness of ecosystems. In order to better estimate such consequences, we need to improve our understanding of the spatial patterns and drivers of specialization in interaction networks. Methods Here, we used the extensive environmental gradient of Mt. Kilimanjaro (Tanzania, East Africa) to study patterns and drivers of specialization, and robustness of plant–pollinator interactions against simulated species extinction with standardized sampling methods. We studied specialization, network robustness and other network indices of 67 quantitative plant–pollinator networks consisting of 268 observational hours and 4,380 plant–pollinator interactions along a 3.4 km elevational gradient. Using path analysis, we tested whether resource availability, pollinator richness, visitation rates, temperature, and/or area explain average specialization in pollinator communities. We further linked pollinator specialization to different pollinator taxa, and species traits, that is, proboscis length, body size, and species elevational ranges. Results We found that specialization decreased with increasing elevation at different levels of biological organization. Among all variables, mean annual temperature was the best predictor of average specialization in pollinator communities. Specialization differed between pollinator taxa, but was not related to pollinator traits. Network robustness against simulated species extinctions of both plants and pollinators was lowest in the most specialized interaction networks, that is, in the lowlands. Conclusions Our study uncovers patterns in plant–pollinator specialization along elevational gradients. Mean annual temperature was closely linked to pollinator specialization. Energetic constraints, caused by short activity timeframes in cold highlands, may force ectothermic species to broaden their dietary spectrum. Alternatively or in addition, accelerated evolutionary rates might facilitate the establishment of specialization under warm climates. Despite the mechanisms behind the patterns have yet to be fully resolved, our data suggest that temperature shifts in the course of climate change may destabilize pollination networks by affecting network architecture.

Plant-pollinator interactions belong to the most frequently studied mutualistic interactions in terrestrial ecosystems (Waser & Ollerton, 2006). The networks share typical topological features, such as high degrees of nestedness, arising from the tendency of specialists to interact with generalists, which tend to interact among each other . Also skewed distributions of links per species resulting from the dominance of few generalists among plenty of species that only interact occasionally (Jordano, Bascompte, & Olesen, 2002), and dependence asymmetry, that is, the differences in mutual dependencies of interacting species (Bascompte, 2006;Blüthgen, Menzel, Hovestadt, Fiala, & Blüthgen, 2007), are network commonalities.
Due to its impact on pollination success and network stability, an interesting feature of plant-pollinator interactions is the species and network specialization. From a plant's perspective, higher specialization on specific pollinators may promote reproductive success and increase genetic diversity, because better morphological adaptations between the pollinator and the reproductive parts of the plant can increase the amount of transferred pollen (Waser & Ollerton, 2006).
In contrast, higher generalization may decrease the dependence on specific pollen vectors and stabilize pollination over broad temporal and spatial scales (Brosi, 2016). From a pollinator's perspective, specialization is a way to reduce interspecific competition and foraging costs, as switching between different search images and handling a diverse range of flower types can be costly (Chittka & Thomson, 1997). On the other hand, specialized pollinators risk investing energy in additional flight time and ignoring lucrative floral resources nearby, which may outweigh the benefits of specialization in energy-restricted habitats and destabilize food safety in environments with high spatial and temporal resource turnover (Waser & Ollerton, 2006). Furthermore, the adaptive potential of specialists to optimize foraging dependent on macronutrient requirements (Vaudo, Patch, Mortensen, Tooker, & Grozinger, 2016) is limited. These trade-offs for plants and pollinators point either toward a strong selection pressure on the degree of floral specialization or, alternatively, require a high plasticity to allow for adaptive foraging, depending on the ecological context (Miller-Struttmann & Galen, 2014;Spiesman & Gratton, 2016).
Factors that may shape the specialization of pollinator communities are, inter alia, the abundance and richness of floral resources, the interaction strengths among pollinators, climatic variables, and the area of available habitat. First, the abundance of floral resources determines average foraging distances and the net energy gain of pollinators (Carvell et al., 2012). Under the assumption that the net energy gain of foraging flights on average decreases with decreasing abundances of flowering plants, specialization may be favored in habitats with high abundances of flowering plants (Kunin & Iwasa, 1996). Flower richness, in contrast, sets the limits for resource partitioning. Specialization should be more likely to occur in habitats offering a variety of resources. Second, a temporarily reduced number of interactions with an otherwise common pollinator species has been shown to broaden the food choice of another pollinator (Brosi & Briggs, 2013). Frequent interspecific interactions, resulting from, for example, high pollinator richness or high visitation rates, may permanently restrict diet breadth promoting species coexistence in the long term (Goulson, Lye, & Darvill, 2008). Third, climate shapes plant and pollinator richness, composition and phenology and has been directly linked to network properties, including specialization (Petanidou et al., 2018).
Temperature determines the costs of foraging flights in ectothermic pollinators (Kovac, Stabentheiner, & Brodschneider, 2015) and may thus modulate resource usage strategies in a way that species broaden their dietary spectrum in energy-limited habitats (Miller-Struttmann & Galen, 2014). Restricted foraging times due to persistent mist and/or temperatures below a threshold in which foraging is possible, should equally result in more generalized foraging. Furthermore, true specialization might establish faster and K E Y W O R D S altitudinal gradient, climate change, ecological network, functional traits, generalization, mutualistic interactions, network specialization index (H 2 ′), pollination, robustness, specialization more often under warm climates, as evolutionary rates accelerate with temperature (Allen et al., 2006;Lin et al., 2019). Finally, habitat area may influence the mean degree of specialization in mutualistic networks (Sugiura, 2010). Specialists typically depend on larger habitat areas than generalists (Bommarco et al., 2010) and might thus become less abundant in high elevations, were habitat area is significantly reduced. Yet, the relative importance of resources, the interactions among pollinators, climate, and area in structuring plant-pollinator interactions remains unclear.
The search for factors that explain specialization is aggravated by the fact that global change can modulate the architecture of species interaction networks (Tylianakis, Didham, Bascompte, & Wardle, 2008). The transformation of natural habitats into arable land and the introduction of invasive species changes species composition drastically and in a short time, requiring permanent adaptations to new interaction partners. Species loss in one trophic level may cause secondary species extinctions in the other level, thereby reducing network robustness, that is, the capacity of a network to buffer such secondary extinctions (Memmott, Waser, & Price, 2004).
Generalization and nestedness may generally increase network robustness, because species have alternative interaction partners, suggesting changing sensitivity of networks to species loss along environmental gradients.
It is assumed that specialization in plant-pollinator networks is linked to functional traits, which restrict species flexibility to switch between different interaction partners (Dehling, Jordano, Schaefer, Böhning-Gaese, & Schleuning, 2016;Stang, Klinkhamer, Waser, Stang, & Meijden, 2009). Broad-scale correlations between specialization and species traits provide important information about such trait-based feedback on specialization, but have hardly been studied on a community level in insects (Albrecht et al., 2018;Lara-Romero et al., 2019). Bees and syrphid flies, for example, differ in their requirements for floral resources. In bees-but not in syrphid fliesthe whole offspring depends on the pollen selection (Praz, Müller, & Dorn, 2008), which might increase bees' selectiveness. Similarly, morphological traits like, for example, proboscis length could restrict the number of potential interaction partners (Ibanez, 2012).
Physiological and energetic constraints are suggested to shape the mean and the variance of species traits along elevational gradients (Classen, Steffan-Dewenter, Kindeketa, & Peters, 2017;Hoiss, Krauss, Potts, Roberts, & Steffan-Dewenter, 2012), indicating that morphological barriers restricting the choice of interaction partners can change with increasing elevation.
Here, we analyzed patterns of specialization of plants and pollinators at different levels of organization in natural and disturbed habitats along a 3.4 km elevational gradient on Mt. Kilimanjaro (Tanzania, East Africa). First, we tested the hypothesis that the structure of plant-pollinator networks changes along elevational gradients. In particular, we hypothesized that species and plantpollinator networks are more specialized in the warm lowlands than in the cool highlands, as energy restrictions should favor generalization in higher elevations. Second, we aimed to explain changes in network structure by a set of major factors that are assumed to have a positive effect on specialization. Using path analysis, we separated the direct and indirect effects of resource availability, pollinator richness and visitation rates, temperature, and habitat area on the community mean of pollinator specialization. Third, we tested whether changes in the specialization of pollinator species are related to taxonomic identity, morphological traits or species elevational ranges. Finally, we explored whether changes in the architecture of plant-pollinator networks lead to a higher sensitivity of some elevational zones to simulated species extinctions.  Table S1). The average distance between study sites was 22.6 ± SD 13.1 km; only two sites were nearer than 2 km (1,920 m), which is still above the average foraging ranges of most pollinators. Along this elevational gradient mean annual temperature varies between 3.1 and 24.0°C, and mean annual precipitation ranges between 590 and 2,740 mm with maximal precipitation around 2,200 m a.s.l. (Hemp, 2006; Figure S1.3).

| Plant-pollinator interactions, species identification, and pollinator traits
In total, we conducted 80 four-hour transect walks on the selected sites, summing up to 320 observational hours. Transect walks were conducted over the course of two consecutive years (2011,2012), covering different seasons of the year. We recorded plant-pollinator interactions between 07.30 and 17.00 hr on days when the weather was sunny or moderately cloudy. In case of rain, mist or heavy wind, we interrupted transect walks and continued it as soon as the weather was suitable again-but at the latest within the next 2 days. Due to logistic constraints and often unsuitable climatic conditions like rain or dense fog at high elevations, the number of transect walks was not homogeneously distributed among sites but ranged between one and eight (Table S1.1). We addressed this by using networks deriving from individual transect walks as sampling units within a mixed-effects model framework (see Section 2.5). This approach ensures that all species contributing to one network co-occurred in space and time.
Additionally, it reduces the susceptibility of network metrics to errors in species identification, because morphospecies were separated only within but not across networks. During each transect walk, we moved slowly and without fixed corridors through the vegetation of each site and recorded each interaction in which a pollinator touched reproductive parts of herbaceous plant species or bushes. If pollinators visited different flowers of the same plant individual before catching, we counted it as single interaction. Note that flower visitors are termed "pollinators" here, although their contribution to the pollination success is unknown. In 95% of all considered interactions, we either identified pollinator species in the field (Apis mellifera), or caught them with sweep nets for further identification by experienced taxonomists. Escaped pollinators, that is, 5% of considered interactions, were also recorded and separated with a conservative approach within single networks (see Appendix S1 for detailed information). Exemplars of interacting plant species, including both herbs and shrubs, were collected or photographed and identified by the botanist AH on a species level. Nomenclature follows the Flora of Tropical East Africa (FTEA, 1952(FTEA, -2012. We restricted our analyses to pollinator taxa, which we could sort on a species or morphospecies level (46% and 54% of all specimens, respectively). Most major groups of pollinators were included, that is, all Hymenoptera: Apoidea: Apiformes ("bees"), the paraphyletic group of nonbee aculeates and symphyta ("wasps"), and Diptera: Syrphidae ("syrphid flies"). Butterflies were excluded from analyses because only relatively few interactions were observed, and the voucher sampling success of the few specimens was poor.
In addition, we excluded nonsyrphid Diptera from analyses, because reliable species delimitation based on outer morphology was not feasible for this group (see also Table S1.4). We further restricted analyses to networks with a minimum of five interactions, but for most networks we could sample a much higher number of interactions (mean number of observed interactions ± SD = 64.9 ± 69.4).
This filtering resulted in 67 networks and led to the exclusion of one network collected in the disturbed Ocotea forest, thereby reducing the number of sites from 19 to 18.
We measured pollinator's proboscis length and head widths using a stereo microscope with calibrated ocular micrometer and a precision of 0.01 mm. Trait matching between proboscis length and corolla tube lengths might restrict the number of interaction partners and has been linked to species specialization before (Miller-Struttmann et al., 2015). Head width, used as a proxy for body size (Branquart & Hemptinne, 2000), is related to energy requirements and foraging ranges (Greenleaf, Williams, Winfree, & Kremen, 2007) and can thus be related to species specialization.
To calculate trait means, we measured the traits of all available, but not more than 10 individuals per species and study site (three individual measurements for syrphid flies) and averaged those values per species. We assessed elevational ranges of pollinator species by subtracting the minimum from the maximum elevation of occurrence. See Appendix S1 for more details on trait measurements and range estimations.

| Resources, pollinator richness, climate, area, and land use intensity
After each transect walk, we counted total flower abundance and flower richness within 10 4 × 5 m rectangles and used the sum of all flower heads and the total number of species as estimates for flower abundance and flower richness per transect walk. Replicated pan trap sampling across seasons was used to estimate networkindependent species richness of pollinator species per site ( Figure   S1.2a). Pan traps were installed in the same years when transect walks were conducted. Sampling effort was equal on all study sites here. Species richness of the pan trap sampling was correlated with the total species richness of pollinators per site collected with nets (Pearson correlation, r = .7, p < .001). Visitation rates were calculated by dividing the number of observed interactions per transect walk (as a measure for insect activity) by flower abundance. Temperature was recorded on each site using temperature loggers (Appelhans et al., 2016). We calculated two temperature measures: For the mean annual temperature (MAT), we averaged all measurements per study site. For the mean temperature during each transect walk (ACT = "actual temperature"), we averaged all measurements during each 4-hr transect walk. Mean annual precipitation (MAP) was interpolated for every study site using a kriging approach on the basis of 15-year-long data records from a network of about 70 rain gauges on Mount Kilimanjaro (Appelhans et al., 2016).
Study sites differed in their land use intensity. To account for that, we used a quantitative composite index of human land use, hereafter termed LUI, which was designed in earlier studies based on data collected on 60 study sites (Peters et al., 2019).
In a nutshell, we averaged standardized estimates of (a) annual plant biomass removal and (b) agricultural inputs (irrigation, fertilization, insecticides, fungicide, and herbicides) and quantified (c) differences of the vegetation structure to the natural vegetation (quantified in terms of canopy closure, canopy height, vegetation heterogeneity) on each study site. Changes in structural characteristics with elevation are partly natural and independent of land use intensification (e.g., different canopy cover in savannah and forest) such that raw data of vegetation structure would not be an informative indicator of land use. Therefore, we calculated the mean Euclidian dissimilarity of vegetation structure measures of the respective study site to the average vegetation structure in natural habitats at the same elevational level. We further quantified (d) landscape composition 1.5 km around each study site-that is, the proportion of areas with agriculturally managed habitats.
We standardized all components (a-d), before averaging them to the final index. More details on this index and other variables (resources, pollinator richness, and temperature) are given in Appendix S1 and Peters et al. (2019).
Elevational belt area was extracted from a digital elevation model of Mt. Kilimanjaro with a resolution of 30 m and an extension from 37.00074 to 37.75602 E and 3.507533 to 2.750183 S (Appendix S1: Figure S1.2b). We calculated the area within a range of 100 m below and 100 m above the respective study site.
For nestedness, we chose the "weighted nestedness overlap and decreasing fills" metric for quantitative networks (Almeida-Neto & Ulrich, 2011). Matrix size equals the product of the number of plant and pollinator species included in each network. Dependence asymmetry ranges between −1 and 1; positive values indicate higher dependences in the higher trophic level (Bascompte, 2006;Blüthgen et al., 2007). High values of nestedness indicate high tendencies of specialists to interact with generalists, which tend to interact among each other, while values around zero indicate the opposite. As absolute values of nestedness partly depend on network size, we standardized them by comparing observed values with results of null models (Dormann, Fründ, Blüthgen, & Gruber, 2009; Appendix S2).
We quantified the degree of specialization at different levels of biological organization: species specialization was calculated for each pollinator and plant species by the d' index, also implemented in the R package "bipartite." The d' index ranges between zero and one, indicating maximal generalization and maximal specialization, respectively. It describes to which extent the observed interaction frequencies of plant and pollinator species deviate from expected frequencies based on random pattern of interactions considering the total frequency of interactions of each partner available (Blüthgen, Fründ, Vázquez, & Menzel, 2008). Compared with alternative specialization indices, d' is relatively robust to observation effort, that is, the specialization of rarely observed species is not overestimated (Poisot, Canard, Mouquet, & Hochberg, 2012).
We averaged d' both on a community and on a species level. The community mean of d' is the average specialization of all plants/pollinators contributing to one network (i.e., all interactions sampled during one transect walk), whereas the species mean d' is the average specialization of single pollinator species across different sites (see Appendix S2 for details). As also the dispersion of d' in a community might be of ecological relevance, we extracted the standard deviation of d' in communities and divided it by the respective community mean (=coefficient of variance [CV]).
For entire networks, we calculated Blüthgens' network specialization index H 2 ′ using the H2fun function implemented in the "bipartite" package (Blüthgen et al., 2006;Dormann et al., 2009). H 2 ′ describes complementary specialization on a network level and has been shown to be robust against sampling intensity and network size, making it a useful tool for the comparison of networks across multiple habitats. Five networks had to be excluded from analyses because we recorded only one pollinator or one plant species (three cases) during the 4-hr walks, or because the distribution of interactions in the network matrix did not allow the generation of more than one random (shuffled) network, as required for the calculation of H 2 ′ (two cases). To confirm that specialization patterns (d′ and H 2 ′) are not driven by network properties like the observed number of interactions or matrix size, we additionally compared those indices with indices derived from null models (Table S2.2a).
Finally, we estimated network robustness against pollinator and plant extinctions for each network using the robustness function of the "bipartite" package. Network robustness equals the area below a secondary extinction curve. This was derived from the stepwise, random removal of species from one trophic level-by setting all entries of this species to zero-, and counting the number of secondarily extinct species from the other trophic level, that is, species with no interactions remaining. As also network robustness may partly depend on network size, we again standardized this metric through comparisons with null models, as described for nestedness in Appendix S2.

| Statistical analyses
We analyzed how (log-transformed) matrix size, dependence asymmetry, nestedness, H 2 ′, community mean of d′ of pollinators and plants and their CV, as well as standardized network robustness (against simulated pollinator and plant extinction) changed along the elevation gradient. We fitted linear mixed-effects (lme) models with elevation as single predictor variable and added study site as a random term, to control for repeated measurements on the same sites. R 2 values were obtained using the Nakagawa & Schielzeth approach (Nakagawa & Schielzeth, 2013) implemented in the R package "r2glmm" (Jaeger, 2016). We conducted several sensitivity analyses to test the stability of detected patterns: (a) As LUI may significantly influence patterns of biotic variables along the elevation gradient (Peters et al., 2019), and because LUI is correlated to elevation (Pearson's r = −.57; Table S1.3), we examined if LUI is potentially a better predictor of network metrics than elevation. We did this by comparing models with elevation as a predictor variable with models that included LUI as single predictor variable using the Akaike information criterion (AIC; Table 1). Since the sample size was low compared with the number of estimated parameters, we employed the AIC C with a second-order bias correction quantifying model support. (b) Network metrics calculated from single sampling days might be strongly influenced by actual weather conditions and might not be representative of the "average" interaction patterns.
To verify the robustness of our results in this respect, we additionally calculated the mean of each metric per study site and analyzed the elevational pattern with ordinary linear models (N = 18; Table   S2.2b). (c) We further tested whether the comparatively poorly tor-driven specialization). Specialization might be a strategy to avoid competition among pollinators (Brosi & Briggs, 2013). Interaction frequencies will increase either with pollinator richness and/or with visitation rates, that is the ratio between the abundance or activ-  Note: All models were fitted with elevation as fixed factor and study site as a random term. ΔAIC C gives AIC C differences of the presented model to a model that includes LUI as single fixed factor. A negative ΔAIC C indicates that the LUI model performed better than the model with elevation; |ΔAIC C | ≤ 2 indicates, that the two models were similarly supported by the data.
Abbreviations: df, degrees of freedom; G, number of study sites; N, number of networks included in analysis; R 2 , semipartial R 2 for the fixed effect; SE, standard error.
least three different sites (43 species). We used generalized linear mixed-effects models with Gaussian error distribution, MAT as single fixed factor and species and site as crossed random effects (Bates, Mächler, Bolker, & Walker, 2015). As this analysis might be especially error-prone toward species misidentification, we run the same model also on a subdataset that only included species with proper species names (20 species).
Differences in pollinator species specialization of different pollinator groups, that is, Hymenoptera (bees and wasps) versus Diptera (syrphid flies), as well as the association of species d' with proboscis length, head width, and elevational range size were investigated across all pollinators and within pollinator groups with linear mixed-effects models. We partly controlled for nonindependence between species by integrating taxonomic information as nested random terms in the models, that is, "genus" for differences between pollinator groups, and "order/family/genus" for differences between species traits. We restricted trait analyses to pollinators that we observed at least three times (n = 87).
We tested if network robustness against plant or pollinator extinction is influenced by H 2 ′, standardized nestedness, MAT and LUI, which were included as additive explaining variables using lme models. Study site was added as random term here again to account for the hierarchical structure of the data. For model simplification, nonsignificant terms (p > .05) were successively removed from the model.

| Elevational patterns in plantpollinator networks
Networks (N = 67) included between 1 and 32 pollinator species and between 1 and 14 plant species. The relative proportions of bees, hoverflies, and wasps in the pollinating community did not significantly change along the elevational gradient (p > .05).
Matrix size ranged between 2 and 448 (mean: 78.9 ± 85.4), and significantly decreased along the elevational gradient ( The community mean of pollinator and plant specialization (d′), as well as network specialization (H 2 ′) decreased with elevation ( Figure 1a-c, Table 1). Comparisons with null models indicated that these declines of specialization are not driven by network size or interaction frequencies alone (Table S2. Table S2.2b). Second, significant declines of network and mean pollinator (but not plant) specialization with elevation were not only depending on the inclusion of high elevation sites but also detectable for the more intensively sampled study sites below 2,000 m a.s.l. (p < .05).
Third, the decline in H 2 ' was also detected when five networks per sites were lumped to one joint network, indicating that the pattern was robust to differences in sampling effort (lm, N = 10, t = −2.631, p < .03).
Finally, the decline of specialization in networks, as well as in pollinator and plant communities, was stable when analyzing only sites with low canopy cover (≤median canopy cover), indicating that vegetation structure did not influence the results (all p < .05).
While the d′ of pollinators on average decreased with elevation Interestingly, a decline in pollinator specialization was also detected at the intraspecific level as the level of specialization of pollinator species decreased along the temperature gradient, that is, species that showed specialized foraging behavior in warm areas, tended to forage more generally in cooler regions (lmer, df = 39.76, t = 4.746, p < .001). This trend was also found when restricting analysis to species with proper species name (lmer, df = 57.95, t = 2.76, p = .008). However, the proportion of explained variance was low in both cases (R 2 = .13 and R 2 = .09, respectively).

| Drivers of specialization in pollinator communities
Mean annual temperature, pollinator richness, flower richness, and habitat area were significantly related to the drop in the community mean of pollinator specialization, when analyzing these variables in separate linear mixed-effects models (Table S2.3). When we added these variables in a joint mixed-effects model and evaluated the support for the full model and all nested models, models including only MAT or MAT and flower richness received the highest level of support (Table S2.3). The best path model revealed that MAT was the p-values > .05 for C indicate that the specific causal structure reflects the data properly. ACT, actual temperature; LUI, land use intensity, MAP, mean annual precipitation; MAT, mean annual temperature; area = habitat area (100 m above and 100 m below the respective study site). All variables were z-transformed prior to analyses. Statistical details are given in Table S2.3

| Elevational patterns in plant-pollinator networks
Here, patterns ranged from higher pollinator specialization in higher (Schleuning et al., 2012), or lower latitudes respectively (Trøjelsgaard & Olesen, 2013), to no trend at all (Ollerton & Cranmer, 2002). In view of the fact that temperature has a big impact on specialization, we suggest that besides the difficulty to standardize studies properly with different sampling strategies for latitudinal meta-analyses, the dominance of endothermic pollinators (e.g., hummingbirds, bats) toward the tropics might partly explain these discrepancies. For example, in hummingbirds species richness, contemporary precipitation and quaternary climate change velocity, but not MAT, predicted specialization (Dalsgaard et al., 2011). Depending on the ratio of endothermic and ectothermic pollinators considered in meta-analyses, the role of temperature in structuring specialization patterns might become weaker or stronger causing different patterns along latitude.
With this we excluded some flower visitors, which are known to contribute to pollination like nonsyrphid Diptera (Orford, Vaughan, & Memmott, 2015). These Diptera were more abundant in higher than in lower elevations abundant in high than in low elevations.
However, as they are known to be as specialized as syrphid flies (Orford et al., 2015), and, in accordance with syrphid flies, less specialized than Hymenoptera (Benadi, Hovestadt, Poethke, & Blüthgen, 2014, and our results), the exclusion of this group in our analyses unlikely affected the general direction of the reported patterns of specialization.

| Drivers of specialization
Along the slopes of Mt. Kilimanjaro, MAT was the best predictor of mean pollinator specialization. Even within pollinator species, individuals tended to forage more generalized in the cold highlands than in the warm lowlands. The mechanisms behind the link of temperature and specialization remain elusive. However, one major hypothesis is that is the increase of metabolic costs in foraging flights under cool temperatures (Kovac et al., 2015;Stabentheiner, Vollmann, Kovac, & Crailsheim, 2003), triggers generalized foraging. Even within pollinator species, individuals tended to forage more generalized in the highlands than in the lowlands, supporting the concept of energetic constraints. Long-tongued bumble bees species also showed comparable shifts in intraspecific foraging behavior in alpine habitats (Miller-Struttmann & Galen, 2014). The authors suggested that restricted seasonal lengths drive generalized resource usage of bumble bees in cold habitats, which is plausible in temperate, but not in cold tropical mountain habitats that lack distinct seasonality. We suggest that in tropical highlands energy-limitations are caused directly by temperature, as temperature controls the foraging costs for ectothermic pollinators (Kovac et al., 2015;Stabentheiner et al., 2003). Furthermore, in the alpine and subalpine zone, we observed that foraging of bees was restricted to very few warm and cloud-free hours of a day, while in the lowlands bee foraging took place all day long. Alternatively or in parallel, higher evolutionary rates under warm climates might favor the evolution of specialists. For instance, bumble bee species within the same subgenera have been shown to evolve faster in low elevations than in the highlands (Lin et al., 2019). This parallels with higher bee richness in the lowlands (Classen et al., 2015;Hoiss et al., 2012) and a higher chance for the evolution of true specialists. In addition, we detected a potential indirect effect of MAT on specialization via flower richness. Flower resource diversity, and a sufficient amount of flowers per resource type, might be a premise for pollinator specialization, as specialized pollinators need to acquire enough food.
Against expectations and in contrast to other studies (Dalsgaard et al., 2011;Schleuning et al., 2012), pollinator species richness and visitation rates were poor predictors of pollinator specialization. This may indicate that physiological and energetic constraints set by temperature rather than the interactions among pollinators (e.g., competitive interactions), shape insect pollinator specialization at Mt. Kilimanjaro. This is also reflected by a previously reported decline of intraspecific variation in morphological traits of wild bees with increasing elevation (Classen et al., 2017): if reduced competition was the major force shaping these functional traits of bees in high elevations, the opposite pattern, that is, "character release," would be expected.
Despite the high proportion of explained variance in pollinator specialization in the path model (R 2 cond. = .79), and even though we verified our results by a set of robustness analyses, it must be considered that our study is based on correlation analyses and that Extinction Extinction example, Helichrysum species, which are accessible by more or other, more generalized, pollinator species than tubular flowers.
Additional, more detailed analyses of plant-pollinator interactions and their traits in the field in combination with true experiments will be necessary to verify the results of our correlative study.

| Species traits
We showed that Hymenoptera were on average more specialized than Diptera, which is in agreement with reports from the literature (Benadi et al., 2014;Weiner, Werner, Linsenmair, & Blüthgen, 2011) and which can be explained by the fact that bee larval survival depends on the pollen source (Praz et al., 2008). Morphological traits were not related to the degree of species specialization, indicating that traits of one trophic level are not informative to predict specialization (Dalsgaard et al., 2018). Trait matching between trophic levels, in contrast, can sharpen our understanding about how species traits structure species interactions and network architecture (Albrecht et al., 2018;Bender et al., 2018;Dehling et al., 2016). Interestingly, we found no link between species specialization and species elevational ranges. Thus, specialized foraging does not necessarily restrict pollinators' ability to inhabit new habitats, which is in line with our finding that even within species the degree of specialization can be adapted to changing abiotic conditions.

| Robustness
Network robustness against pollinator extinction was generally higher than robustness against plant extinction. This is in line with simulation studies under climate change scenarios . Importantly, network robustness increased with elevation. In other words, network sensitivity to species loss was highest where habitat transformation and land use intensification, as main driving forces of species loss are most prevalent (Nogués-Bravo, Araújo, Romdal, & Rahbek, 2008). Robustness against species extinction was negatively correlated with network specialization, but not with nestedness, as shown in other studies (Bastolla et al., 2009;Rohr, Saavedra, & Bascompte, 2014). Network generalization decreases the co-dependence of interaction partners and is thus an important driver of network robustness. We assessed robustness against random species extinctions of one trophic level, as we currently had no concrete indication for other extinction orders than random.
Especially in human-modified landscapes, abundant and well-connected species are frequent visitors of crops, suffering from pesticide usage, soil cultivation or the bee-keeping associated spread of pathogens. They might nowadays face extinction risks that are comparable with the ones of less-abundant and more specified species. Alternatively, trait-based extinction orders can drastically alter the robustness pattern along elevational gradients and related pollination services in unpredictable ways (Larsen, Williams, & Kremen, 2005). Generally, the robustness metric needs to be handled with care. First, species may adapt to a certain degree to other interaction partners, when their previous partner goes extinct. Our finding that specified species from the lowlands tended to forage more generalized in the highlands, points toward such an adaptive capacity, which will relax the co-dependence of certain species and strengthen network robustness (Kaiser-Bunbury, Muff, Memmott, Müller, & Caflisch, 2010). Second, robustness was calculated for networks that derived from short-term observations. However, many pollinators show floral fidelity, making observed networks more specialized and thus less robust than networks observed over a longer period (Petanidou et al., 2008;Spiesman & Gratton, 2016). Finally, this approach ignores that the loss of one species can change the behavior or another species from the same trophic level, with consequences for interaction partners and network stability (Brosi & Briggs, 2013).

| CON CLUS IONS
Plant-pollinator network architecture strongly responded to changing climate along the slopes of Mt. Kilimanjaro. We identified MAT as the main driving force for specialization, which in turn affected network robustness against species extinction. Rising metabolic costs for foraging flights in cool environments might explain the detected decrease of (ectothermic) pollinator specialization in the highlands.
We expect that the temperature-dependence of metabolic costs affects also the structure of other species interactions including ectothermic organisms (e.g., parasitoid-host, plant-herbivores), with consequences for diversity, ecosystem functions and stability (Plowman et al., 2017).
Although the mechanisms behind the temperature-specialization relationship remain speculative, the outstanding role of temperature in structuring plant-pollinator interactions is alarming: global temperatures are predicted to increase by up to 0.7°C in the next two decades (compared to 1985-2005IPCC, 2013). So far, climate change was expected to disrupt plant-pollinator interactions by causing spatial and phenological mismatches between interaction partners (Hegland, Nielsen, Lázaro, Bjerknes, & Totland, 2009), or by increased frequency of extreme weather events (Hoiss et al., 2015), with negative consequences for interaction resilience and fitness of plants and pollinators (Forrest, 2015;Schenk, Krauss, & Holzschuh, 2018). The direct impact of temperature on the specialization of species, which has also been reported from smallscale climatic gradients for network specialization (Petanidou et al., 2018), imposes additional challenges to species interactions.
The loss of interaction partners may improve pollination quality on small spatial and temporal scales, as less heterospecific pollen will be transported (Brosi, 2016). This, however, may drastically reduce pollination insurance in the long term (Brosi, 2016).
We conclude that rising temperatures in the course of climate change will destabilize species interactions along entire elevational gradients, thereby exerting additional pressure on species, which already live close to their maximum thermal capacity (Colwell, Brehm, Cardelus, Gilman, & Longino, 2008).