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

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.

INTRODUCTION

Interspecific species interactions are known to limit the diversity and distribution of species (Bairey, Kelsic, & Kishony, 2016; Bascompte, 2006; Bastolla et al., 2009; Chan, Shih, Chang, Shen, & Chen, 2019; Wisz et al., 2013), to promote species evolution (Ramos & Schiestl, 2019), and to determine ecosystem functions (Brosi & Briggs, 2013; Garibaldi et al., 2013). Much progress was made in understanding the structure and dynamics of species interaction networks (Bastolla et al., 2009; Blüthgen, Menzel, & Blüthgen, 2006; Olesen, Bascompte, Dupont, & Jordano, 2007; Petanidou, Kallimanis, Tzanopoulos, Sgardelis, & Pantis, 2008). Nevertheless, knowledge about the spatial patterns and drivers of network properties remains surprisingly ambiguous (Morris, Gripenberg, Lewis, & Roslin, 2014; Novotny, 2006; Ollerton & Cranmer, 2002; Schleuning et al., 2012; Song, Rohr, & Saavedra, 2017; Lara‐Romero et al., 2019). Elevational gradients offer the opportunity to study the architecture of species interaction networks along broad climatic gradients at feasible spatial and temporal scales (Hoiss, Krauss, & Steffan‐Dewenter, 2015; Ramos‐Jiliberto et al., 2010; Lara‐Romero et al., 2019). Importantly, on mountains, network metrics, and potential underlying drivers can be measured in a standardized and thus informative manner. This rare information is essential for understanding the evolution and coexistence of species communities and for predicting the functionality of ecosystems under global change (Tylianakis, Laliberté, Nielsen, & Bascompte, 2010).

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 (Bascompte & Jordano, 2007). 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 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 flies—the 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 plant–pollinator 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.

MATERIALS AND METHODS

Study sites

We selected 19 100 × 100 m study sites on the southern slopes of Mt. Kilimanjaro, spanning an elevational gradient from 993 m above sea level (m a.s.l.) up to 4,390 m a.s.l. Study sites covered the major natural and anthropogenic habitat types of Mt. Kilimanjaro: colline savanna and maize fields (990–1,020 m a.s.l.), lower montane forest, agroforestry systems (Chagga home gardens), grasslands, coffee plantations (1,260–1,920 m a.s.l.), montane undisturbed and disturbed by former logging Ocotea forest (2,120–2,470 m a.s.l.), upper montane undisturbed and fire‐disturbed Podocarpus forest (2,850–2,990 m a.s.l.), subalpine Erica forest (3,880 m a.s.l.), and alpine Helichrysum vegetation (3,880–4,390 m a.s.l.; 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–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 network‐independent 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.

Network indices

We used the R package “bipartite” (Dormann, 2011; Dormann, Gruber, & Fründ, 2008) to calculate matrix size, dependence (or “interaction strength”) asymmetry and nestedness for all networks. 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 AICC 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 sampled high‐elevated habitats with generally low interaction rates are responsible for the strong elevational patterns detected in specialization, by restricting analyses to all study sites below 2,000 m a.s.l. (d) To illuminate whether the differences in the number of repeated measurements per site influences the elevational pattern of H 2′, and to test if higher sampling effort per study site leads to different patterns of network metrics along the elevation gradient, we lumped a minimum of five replicated transect walks per site to one joint network (possible for 10 sites, Table S1.1), extracted H 2′ and analyzed the pattern along elevation with an ordinary linear model (N = 10). (e) To evaluate the impact of changes in vegetation structure on elevational patterns, we analyzed the elevational trends in network and species specialization for a subset of all study sites in which the canopy cover is smaller or equal to the median canopy cover (i.e., all “open” habitats).

Table 1

Outputs of linear mixed‐effects models, showing the changes in network metrics along the elevational gradient on Mt. Kilimanjaro

Network metric	Predictor	N	G	SD rand.	Estimate	SE	df	t‐Value	p‐Value	R 2	ΔAICc
log (matrix size)	Intercept	67	18	.69	5.12	.42	49	12.28
	Elevation				−8.3E–04	1.7E–04	16	−4.86	<.001	.37	−8.74
Dependence asymmetry	Intercept	67	18	.10	.25	.08	49	3.16
	Elevation				−9.24E–05	3.675E–05	16	−2.52	.023	.11	0.99
Std. nestedness	Intercept	61	16	.52	−4.17	.56	45	−7.45
	Elevation				8.9E–04	2.9E–04	14	3.08	.008	.15	2.09
Mean d′ pollinators	Intercept	66	17	.09	.64	.05	49	11.98
	Elevation				−1.5E–04	2.3E–05	15	−6.60	<.001	.53	13.11
Mean d′ plants	Intercept	66	17	.16	.79	.09	49	8.67
	Elevation				−1.7E–04	3.8E–05	15	−4.57	<.001	.36	8.60
H 2′	Intercept	62	17	.11	.89	.08	45	10.88
	Elevation				−1.9E–04	3.9E–05	15	−4.77	<.001	.32	7.34
Std. robustness (against pollinator extinction)	Intercept	64	17	.54	−3.06	.66	47	−4.67
	Elevation				7.6E–04	3.4E–04	15	2.23	.041	.08	2.91
Std. robustness (against plant extinction)	Intercept	64	17	.56	−4.30	.60	47	−7.16
	Elevation				1.0E–03	3.1E–04	15	3.28	.005	.17	2.24

All models were fitted with elevation as fixed factor and study site as a random term. ΔAICC gives AICC differences of the presented model to a model that includes LUI as single fixed factor. A negative ΔAICC indicates that the LUI model performed better than the model with elevation; |ΔAICC| ≤ 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.

We conducted path analyses to separate the direct and indirect effects of flower abundance, flower richness, pollinator richness, visitation rates, temperature, and habitat area on the community mean of pollinator specialization, d′. We based the a priori structure of the path model on the following hypotheses (Figure 2a). (a) Pollinator specialization is predicted to increase with flower abundance and/or the number of flowering plant species (=resource‐driven specialization), as these factors determine whether specialization is feasible in a habitat, which needs to provide sufficient amount of food for pollinators. MAT, MAP, and LUI will influence flower abundance and flower richness. (b) Pollinator specialization is predicted to increase with the frequency of interactions among pollinators (=pollinator‐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 activity of pollinators and the number of resources. Pollinator richness will depend on MAT, MAP, LUI, area and mean richness and abundance of flowers. Per calculation, visitation rates depend on insect activity and flower abundances, but should further be influenced by ACT, which increases the activity of ectothermic pollinators. (c) Temperature is predicted to have a direct, positive impact on the mean specialization of pollinators (=temperature‐driven specialization). ACT will directly control the costs of an observed foraging flight. In contrast, MAT reflects the temperature that is experienced by communities on a long term and is probably correlated with the number of hours per day suitable for foraging. MAT might further positively correlate with evolutionary rates and with this the likelihood that specialists evolve (Allen et al., 2006; Lin et al., 2019). (d) Pollinator specialization is expected to increase with increasing habitat area (=area‐driven specialization), as specialists typically depend on larger areas than generalists do.

Figure 2

Direct and indirect predictors of mean pollinator specialization on Mt. Kilimanjaro. (a) A priori hypothesized causal structure of the model. Competitive variables within each hypothesis were highlighted with similar colors. Black and colored arrows indicate positive relationship expectations, gray arrows negative relationships. (b) Structure of the full path model after semiautomated preselection of variables. Detailed information on the preselection process are given in the method section. (c) Final path model derived by AICC‐based model selection across all possible paths combinations presented in b. Path coefficients and related p‐Values, as well as both marginal and conditional R 2 values for all response variables are presented. Dashed lines indicate nonsignificant paths. The presented path model is statistically not distinguishable from a model in which flower richness has no impact on pollinator specialization (ΔAICC = 0.05). The global goodness of fit of all path models was estimated with Fisher's C. 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

All variables included in the path model were standardized by z‐transformation using the scale function in R. To avoid correlation between exogenous variables and to improve the ratio between the number of model parameter and the number of observations, we selected variables within and across hypotheses via model selection based on AICC prior to path analysis (Table S2.3). After variable selection, all possible path combinations were tested and compared by the AICC of the respective path model. We used Fisher's C as a goodness of fit parameter. Statistical details on the preselection process and on path analyses are given in Appendix S2.3.

Intraspecific variation of pollinator specialization along the temperature gradient was analyzed for pollinators that occurred on at 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.

RESULTS

In total, we analyzed 4,380 plant–pollinator interactions (3,757 bee, 196 wasp, and 427 hoverfly interactions). We identified 141 plant species (123 species, two species complexes, 16 morphospecies), and 187 pollinator species (84 species, 103 morphospecies). More details on the taxonomic resolution are given in the Appendices S1 (Table S1.4) and S3.

Elevational patterns in plant–pollinator 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 (Table 1). Dependence asymmetry declined with elevation: In the lowlands, positive values of dependence asymmetry indicated higher dependencies of pollinators on plants, while from an elevation of about 2,730 m a.s.l. these dependencies, on average, turned around. Nestedness increased with elevation (Table 1). LUI, which was negatively correlated with elevation (r = −.57), was in case of one network metric a better predictor variable than elevation (Table 1): Matrix size increased with increasing LUI (t‐value: 6.949, p < .001, R 2 = .61).

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.2a). Also, sensitivity analyses confirmed the robustness of the majority of reported patterns: First, all elevational patterns in network metrics were also statistically significant, when analyzing metric means per study site with linear models (N = 18; 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 H2' 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).

Figure 1

Change of plant–pollinator specialization along the elevational gradient on Mt. Kilimanjaro at species and network level. (a) Community mean of pollinator specialization (d'), (b) community mean of plant specialization (d') and (c) plant‐pollinator network specialization (H 2′) decreased with increasing elevation (m a.s.l. = meters above sea level). Dots represent the abundance‐weighted means of species specialization indices (d′) and the H 2′ values per transect walk. Lines represent predicted relationships derived from linear mixed‐effects models with elevation as single predictor variable and site as a random term. Dot colors indicate the strength of land use intensity

While the d′ of pollinators on average decreased with elevation the coefficient of variation (CV) of pollinator d' increased (pollinators: t = 4.76, p < .001, R 2 = .27), indicating more variability in the level of specialization in higher elevations. However, this trend was not detectable for study sites below 2,000 m a.s.l. (p > .05). The CV in the plant d' did not significantly change with elevation (t = 1.59, p = .135).

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 strongest predictor of the community mean of pollinator specialization (cond. R 2 = 79%; Figure 2c). MAT had both a direct positive effect on pollinator specialization, and an indirect positive effect via flower richness, which was weak compared to the direct effect (Figure 2c). In the best path model, (log‐transformed) flower abundance had a positive effect on flower richness. We tested a set of competing alternative path combinations and ranked path models according to their AICC values. A path model that did not include the link from flower richness to pollinator specialization was statistically not distinguishable from the presented model (ΔAICC = 0.05).

Pollinator specialization and species traits

Pollinator specialization (d′) was generally higher in Hymenoptera than in Diptera (lme, t = 2.312, p = .043, Figure 3a), but was neither related to species elevational range sizes (Figure 3b), nor to functional traits (glossa lengths: lme, t = 0.406, p = .687; head width: lme, t = 0.305, p = .762).

Figure 3

Impact of taxonomy and species elevational range on pollinator species specialization (d′). (a) Hymenoptera (bees and wasps) were on average more specialized than Diptera (syrphid flies; t = 2.70, p = .010; bees‐ wasps: t = 0.03, p = .74; bees—syrphid flies: t = −2.52, p = .016; syrphid flies—wasps: t = 3.36, p = .002). While black box plots present common summary statistics (with data medians as black line and means as white asterisks), the surrounding violin plots signal (smoothed) probability density of the data at different values. (b) Pollinator specialization was not related to the elevational range size of pollinator species. Dot color in (b) corresponds to the different taxonomic groups, as introduced in (a). We considered only pollinators that we observed at least three times (87 species)

Robustness

Network robustness against plant extinction was on average lower than network robustness against pollinator extinction (lme, t = −2.70, p = .008), and both increased with elevation (Figure 4, Table 1). Among MAT, LUI, H 2′, and nestedness, network specialization (H 2′) was the best predictor of network robustness, with less specialized networks being more robust (lme, n = 62, df = 44, against plant extinction: t = −3.48, p < .001, R 2 = .23; against pollinator extinction: t = −2.94, p = .004, R 2 = .14).

Figure 4

Network robustness against species extinction along the elevational gradient on Mt. Kilimanjaro. Robustness against pollinator extinctions (black triangles) exceeded robustness against plant extinction (white dots); both metrics increased with elevation. Robustness was standardized via null‐model comparison prior to analysis. Lines represent predicted relationships derived from linear mixed‐effects models with elevation as single predictor variable and site as a random term

DISCUSSION

Plant–pollinator networks showed clear broad‐scale patterns along the elevational gradient of Mt. Kilimanjaro, highlighting the importance of ecological and/or evolutionary drivers in structuring species interactions. Specifically, we found an increase in generalization with elevation on the species (d′ intraspecific), community (community mean of d′), and network level (H 2′). Similar patterns in the specialization of plant–insect pollinator interactions have been reported from biogeographically independent regions, including the Alps (Hoiss et al., 2015, network level), Colorado Rocky Mountains (Miller‐Struttmann & Galen, 2014, intraspecific and network level), the Andes (Ramos‐Jiliberto et al., 2010, network level), and an island volcano in Tenerife, Canary Islands (Lara‐Romero et al., 2019, network and community level). Far less consistent are the patterns reported from plant–hummingbird interactions along elevational gradients, which were either more generalized in the highlands (Maglianesi, Blüthgen, Böhning‐Gaese, & Schleuning, 2015), or in the lowlands (Dalsgaard et al., 2018). Similarly contradicting are the results from meta‐analyses describing specialization patterns along other temperature gradients, that is, latitudinal gradients. 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.

In our study, we concentrated on three important insect pollinator groups: bees, nonbee aculeates (“wasps”), and hoverflies. 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 the power of statistical tools to separate variables that are correlated in nature (e.g., area and MAT; pollinator richness and MAT), are restricted. We also cannot rule out that factors correlating with temperature, which we could not consider in this study, affected specialization. For example, it is conceivable that pollinator generalization in the highlands coincides with the dominance of plant species with open flowers and short nectar‐holding tubes like, for 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.

Network robustness against pollinator extinction was generally higher than robustness against plant extinction. This is in line with simulation studies under climate change scenarios (Schleuning et al., 2016). 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).

CONCLUSIONS

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–2005; IPCC, 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 small‐scale 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).

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

ISD, MKP, and AC designed the study; AC conducted the fieldwork, compiled the data, and conducted statistical analyses with support of MKP. AH selected the study sites and identified all plant species. Insects (bees, syrphid flies, wasps) were identified by CDE, AS and RSP, respectively. AC wrote the first draft of the manuscript, while all authors contributed to the final version.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.