Bumblebee diversity and pollination networks along the elevation gradient of Mount Olympus, Greece

We studied bumblebee diversity and bumblebee pollination networks along the altitudinal gradient of Mt. Olympus, a legendary mountain in Central Greece, also known for its exceptional flora.

cies, remaining active in cool conditions when other pollinators are inactive (Heinrich, 1979). As a result, bumblebees are often the most important pollinators in cool environments, such as the Arctic (Løken, 1973), montane (Bingham & Orthner, 1998) and temperate (Ricketts et al., 2008) regions. Bumblebees are present at altitudes even above 4,000 m a.s.l., providing invaluable pollination services to the alpine flora (Williams, 1991). A rise of the ambient temperature may act against their adaptive advantage to exploit floral resources at higher altitudes by increasing competition with other pollinators, with unknown consequences for the current species distributions (Martinet et al., 2015). This is particularly likely in the eastern Mediterranean, because although it has long been considered a hotspot for wild bee diversity (Michener, 1979(Michener, , 2000, the diversity of bumblebees is comparatively low (Balzan et al., 2016;Nielsen et al., 2011;Petanidou & Ellis, 1993;Varnava et al., 2020). Research on bumblebee diversity in Mediterranean mountains is therefore paramount in order to assess the role of the genus in structuring mountain communities at different elevations and disclose the conservational importance of these habitats for the genus.
Here, we study the patterns of diversity and pollination networks of bumblebees along an altitudinal gradient of Mt. Olympus  (Strid, 1980;Strid & Tan, 1986, 1991. The mountain retained a certain level of glaciation in Pleistocene interglacial periods and even up to the Holocene (Smith, Damian Nance, & Genes, 1997;Styllas et al., 2018), thus constituting a refuge for glacial relics (Médail & Diadema, 2009), such as the local endemic Jankaea heldreichii (Boiss.) Boiss. (Gesneriaceae) characterized as a "living fossil" (Vokou, Petanidou, & Bellos, 1990). Because of its high altitude, the rich zonation along its elevation gradient, separation from other mountain ranges, high diversity of flora and fauna and glacial history, Mt. Olympus is considered a model system for ecological, evolutionary and biogeographical studies (Strid, 1980). High biodiversity and endemism suggest population differentiation and speciation in past geological times when Mt. Olympus was a species refuge (Hewitt, 2004;Médail & Diadema, 2009). Until today, mainly plant-focused studies have been conducted on the mountain, with at most partial consideration of plant-pollinator interactions, pollinator diversity, and community structure, and only at particular altitudes (Blionis & Vokou, 2005 and references cited therein; Makrodimos, Blionis, Krigas, & Vokou, 2008;Vokou et al., 1990).
In this study, we aim to explore the bumblebee patterns of diversity and function, that is involvement in pollination services within communities occurring along the entire altitudinal gradient of Mt.
Olympus. This approach helps to disclose the present and historical biogeography of the bumblebees in the Mediterranean, and highlights conservation priorities regarding both species and habitat types. It is noteworthy that among all studies investigating both pollinator distribution/diversity and function along a mountain's altitudinal gradient, only two focused on bumblebees: Egawa and Itino (2019), who studied the elevation gradient 700-2,600 m a. s. l. of Mt.
Our study is unique in that it explores bumblebee diversity patterns and involvement in pollination networks along the entire altitudinal gradient of a mountain.
In particular, we address the following questions:

| Study area
The study was conducted on the north-eastern slope of Mt. Olympus, following the path route from Litochoro town to Mousses Plateau ( Figure 1). Along this route and at different altitudes and vegetation zones (see below), we established ten sampling sites with a NNE aspect, including both open areas and forest openings when inside a forest; between successive sites, the average altitudinal difference was ca. 250 m and average distance in straight line ca. 1.5 km; in the two cases the distance between sites was <1 km (in fact between 500 and 600 m), these sites were taken by necessity as the only ones available ( Figure 1, Table S1). Almost half of the mountain, specifically its eastern side, is part of the Olympus National Park, historically the first national park ever designated in Greece (1938), part of which constitutes a European Natura 2000 site (GR 1250001). Our study sites cover all four major vegetation zones of Mt. Olympus, according to Strid (1980): (a) evergreen-sclerophyllous (Mediterranean) scrub (300-600 m), (b) mixed beech and montane coniferous forests (700-1,500 m), (c) cool temperate coniferous forests (1,500-2,500 m) and (d) alpine meadows (2,500-2,917 m). The climate is Mediterranean at lower altitudes and temperate at intermediate to higher ones, with snow covering areas above 2,000 m from late October to late May (Strid, 1980). Afternoon thunderstorms are common in summer, and sudden snowstorms or hail may also occur. There is no permanent water source above 1,100 m; thus, humidity from fog or low clouds is of crucial importance. North-eastern slopes are more humid due to their vicinity to the sea (Strid, 1980).

| Sampling of insects and their flowering plant hosts
Insect surveys were carried out within an area of ca. 0.1 ha at each site; they consisted of a collection of all insects using pan traps (years 2013 and 2014) and additional random transect observations assisted by hand netting collection of unknown bumblebees (years 2013, 2014 and 2016). Pan traps were placed as sets of three plastic bowls (triplets), each pan trap painted with a different UVbright colour: blue, yellow and white (Nielsen et al., 2011;Westphal et al., 2008). Triplets, five in each site and set at a distance of >5 m from each other, were placed on height-adjustable metal stands in order to be levelled with the height of the continuously growing vegetation during the sampling season, and thus being always visible to flying insects. Each pan trap was filled with a mixture of propylene glycol and water (1:1 in volume) and was left open to capture insects for a week. The use of propylene glycol was adopted for two reasons: to minimize evaporation and to preserve the collected specimens, thereby allowing a longer operational time of the pan traps.
Bumblebee surveys were carried out using hand netting in variable transect walks during which both bumblebees and the visited flowering plants were recorded (Nielsen et al., 2011;Westphal et al., 2008). Specimens of unknown bumblebees were collected together with the respective forage plants for later identification. The observation and collection of specimens was conducted in days with fine weather favourable to bee activity and during the active bumblebee foraging hours (10:00-16:00). During the survey, we tried to collect and observe bumblebees as thoroughly as possible and to disturb them as little as possible. If occasional disturbances may have occurred and some interactions missed, we trust this was evenly distributed over all sites along the altitudinal gradient. As some bumblebee species are long-distance flyers, movement between sites could potentially have occurred, although we believe this was most unlikely to happen due to the topology and the existing dense forests between sites. Several factors, like difficult accessibility (all sites above 1,200 m had to be accessed on foot), accommodation restrictions (only one refuge along the path) and the erratic mountain weather, influenced the data collection, affected both the choice of sampling methods and the sampling effort employed within a year (once per month) and in different years (variation between years). To have comparable results, such differences (viz. different sampling methods and duration) were taken into account in the analyses. Finally, because sites were selected not far from the hiking path, carefully prepared bilingual informational signs were set to minimize disturbance by hikers.
All collected insect specimens were processed in the laboratory and identified to species level, except for Bombus lucorum s.l., a complex consisting of three species that are indistinguishable with a microscope: B. lucorum, B. cryptarum and B. magnus; yet, because B. terrestris was also difficult to distinguish under field conditions, visit records of all the above species were pooled together as belonging to the same group (Bossert, Gereben-Krenn, Neumayer, Schneller, & Krenn, 2016;Wolf, Rohde, & Moritz, 2009). Bumblebee specimens are kept at the Melissotheque of the Aegean (Laboratory of Biogeography and Ecology, Department of Geography, University of the Aegean, Greece) . Voucher plant specimens are kept at the Herbarium of the Laboratory of Systematic Botany, Agricultural University of Athens, Greece.

| Plant diversity and flower cover
Number of plants in flower and their flower cover were measured within 1 × 1 m squares (n = 25) randomly selected at each site and round .

| Bumblebee diversity
Using the combined pan trap and hand net data for the years 2013 and 2014 separately, we computed three metrics that describe the diversity of bumblebee species in the study communities: (i) Species richness: the total number of species within a community.
(ii) Shannon-Wiener diversity index (Hʹ): a widely used metric of α-diversity, which accounts both for the diversity and evenness of species within a community.
(iii) Effective Number of Species (ENS): defined as the number of species in a community with equal abundances that would give the same as the observed value of α-diversity, and considered as the true species diversity within a community (Jost, 2007).
Metrics (ii) and (iii) were calculated using the R package vegan 2.5-2.

| Network properties
We calculated the following network properties, using all hand net- (ii) Nestedness: the weighted NODF was used, which describes the degree to which specialist species tend to interact with generalist ones (Almeida-Neto & Ulrich, 2011).
(iv) Shannon's diversity of interactions: the Shannon index calculated for the different links of the network.
(v) Selectiveness: we used the H 2 ʹ index, which describes the selectiveness in the network, that is to which extent observed interactions deviate from random selection of partners. The more selective the species in the network, the higher the H 2 ʹ index (Blüthgen, Menzel, & Blüthgen, 2006).
All the above metrics were calculated using the R package bipartite 2.11.

| Statistical analysis
To test the relationship between altitude and bumblebee community composition, we used multivariate-response generalized linear mixed models (MGLM) (R package mvabund 3.13.1). This approach fitted a separate GLM (binomial family with "cloglog" link) to the distribution matrix of each insect of the network, using altitude as explanatory variable, and a resampling-based hypothesis testing (Wang, Naumann, Wright, & Warton, 2012), and was used to indicate those species whose distribution showed a significant effect to altitudinal change. The multidimensional response variable was the presence/absence matrix of bumblebee species distribution among the ten sites, compiled by using the combined data from pan trapping and hand netting of the years 2013 and 2014. The statistical significance of the fitted model was assessed with ANOVA (likelihood ratio tests) using 999 bootstrap iterations via PIT-trap residual resampling. We used AIC to compare the best-fitted model with the null model (y ~ 1). Univariate tests were subsequently performed to determine which response variables (bumblebee species) showed significant effects in response to altitude.
To test the relationship between altitude and the diversity of bumblebees on Mt. Olympus, we used Generalized Linear Mixedeffects Models (R package lme4 1.1-18-1) considering all pan trapping and hand netting data combined. To avoid pseudoreplication, we used "site ID" nested within "year" as categorical random variables. As dependent variables, we used (i) the bumblebee species richness, (ii) the Shannon-Wiener diversity index and (iii) the Effective Number of Species (ENS), calculated for each of the sampled sites. For Shannon-Wiener index and ENS, we used models with Gaussian distribution, while for the bumblebee species richness, we used Poisson distribution and log link function (in the latter case, the model was not overdispersed; that is, the ratio between the residual deviance and residual degrees of freedom was <1.1 [Zuur, Ieno, Walker, Saveliev, & Smith, 2009]). As fixed variables (predictors), we included both the linear and the quadratic term of the altitude of each site. To avoid collinearity, we standardized the variable "altitude" (x = 0, σ = 1) before calculating its quadratic term.
In addition, we used plant Shannon-Wiener index as a predictor. To select the best combination of predictors (including the interaction terms of plant diversity and altitude), we used a backward selection process based on AICc (R package AICcmodavg 2.1-1). We also used the ΔAICc in order to compare the best-fitted model with the null model (y ~ 1).
The same modelling approach (GLMMs) was used to test the relationship between altitude and species richness of the plants in flower, "site ID" nested within "year" as categorical random variables, and the linear and the quadratic term of the altitude of each site, as predictors. We used Poisson distribution and log link function (in the latter case, the model was not overdispersed; that is, the ratio between the residual deviance to residual degrees of freedom was <1.1 [Zuur et al., 2009]).
To test the relationship between altitude and the b-p visitation network properties, we considered the hand netting data alone and applied GLMMs. To do so, we used the abovementioned network properties (viz. network size, H 2 ʹ index, weighted NODF, linkage density and Shannon's diversity index) as dependent variables. As predictors, we used the linear and the standardized (see above) quadratic term of altitude, and flower abundance. As random factors, we used site ID, nested within year (2013,2014,2016). We performed analyses using (a) the Gaussian distribution for Shannon's diversity and linkage density, (b) binomial distribution (link "logit") for NODF and H 2 ʹ, and (c) Poisson distribution for network size (link "log") (the model was not overdispersed). To determine whether there was a significant altitudinal trend in the network properties, we used the ΔAICc between the best-fitted model and the null model (y ~ 1).

| RE SULTS
During the three study years, we recorded 5,255 bumblebees belonging to 22 species and one species complex (Table S2 in which species identification is of the lowest taxonomic level), which visited the flowers of 94 plant species (Table S3)   Interestingly, bumblebee diversity correlated negatively with plant diversity (Figure 3d), probably as a result of the inverse diversity trends of the two groups along the altitudinal gradient (i.e., diversity of plants peaked at lower to intermediate altitudes, while that of bumblebees at high ones [cf. above and Table S4]).
The structure of b-p networks varied along the altitudinal gradient for each of the three years of sampling (2013,2014,2016) ( Figure 5). Networks were larger (i.e. included more links) between 1,500 and 1,800 m in all sampling years (GLM; χ 2 1,3 = 18.4, p < .0001). (Note that 2013 networks were overall smaller than the other two years, probably due to the lower sampling effort in that year; see Section 2). Linkage density also peaked at ca. 1,800 m (LM; F 1,24 = 20.0, p < .0001); year 2014 did not follow this trend, probably because metrics for the lowest altitude site (S-327) could not be computed (only one bumblebee species was recorded, viz. B. terrestris/lucorum s.l. group) (cf. Figure 5). Shannon's link diversity peaked at around 1,800 m (LM; F 1,24 = 13.5, p = .001). Finally, the b-p networks were more nested at the extreme altitudes of the mountain, that is <600 m and >2,400 m, showing an opposite trend than the other properties (GLM; χ 2 1,3 = 7.1, p = .008). The selectiveness index (H 2 ʹ) was found to correlate neither with altitude nor with flower abundance. In fact, no property was correlated with flower abundance in the study sites: based on AICc, flower abundance was not included in any of the best-fitted models.

| Bumblebee diversity and faunistics
Three years of systematic survey in ten sites along the altitudinal gradient of Mt. Olympus yielded 22 identified species and one species complex of bumblebees, consisting of three cryptic species.
Such a rich bumblebee diversity is comparative with that recorded in a much variable set of temperate habitats in NW Greece surveyed for a wider period of time, which encompassed a total of 29 species out of the 33 currently known for Greece (Anagnostopoulos, 2005(Anagnostopoulos, , 2009, 2017; cf. Table 1). Mediterranean systems are notorious for their paucity in bumblebee diversity (Balzan et al., 2016;Nielsen et al., 2011;Petanidou & Ellis, 1993;Varnava et al., 2020) as compared with temperate regions in the European continent (Table 1) The present study discovered one new species for Greece, viz. Bombus quadricolor which is a rare parasite specialist of B. soroeensis. This parasitic species has a wide distribution range from the Cantabrian Mountains and the Pyrenees through to the Balkan mountains and Fennoscandia, with Mt. Olympus constituting its southernmost distribution area . Furthermore, this study documented the range expansion of another three bumblebee species, for which Mt. Olympus constitutes their southernmost distribution in Europe: B. monticola, whose range extends from the Pyrenees to Fennoscandia, including the Alps and some Balkan mountains (Martinet et al., 2018); B. hypnorum, a species with a wide distribution, from the Pyrenees and Balkan mountains in the south, to the Barents Sea in the north ; and B. haematurus, a species that reaches Slovakia and Romania in the north (Bossert & Schneller, 2014;Rasmont et al., 2015). Given their latitudinal restriction, the above four species can be considered as vulnerable to climate change, especially B. monticola being restricted to the highest sites of Mt. Olympus (>2,400 m a.s.l., cf. Table S2) (Kerr et al., 2015;Rasmont et al., 2015). Olympus as compared to other higher mountains and its low latitude, and (d) climatological factors, for example the influence of the Atlantic to Cantabrian Mountains.

F I G U R E 4 Altitude as predictor of plant species richness on Mt.
Olympus. The three fitted lines correspond to different sampling years, which represent random-effect factors in the fitted models (***p ≤ .001) Six of the above seven species (viz. B. mesomelas, B. pratorum, B. ruderarius, B. rupestris, B. soroeensis, B. vestalis) are associated to high altitudes of Mt. Olympus (cf. Figure 2, Table S2 showing the abundances of these species in different altitudes). Among all these species, however, only B. mesomelas is exclusively associated with very high altitudes. Furthermore, as shown in Figure 6,  There is hard evidence on the decline of bumblebee species worldwide (Cameron et al., 2011;Kosior et al., 2007) which is associated not only with climate change, but also with land use change, pesticide use, pathogens, invasive species and combinations among these factors (Cameron & Sadd, 2020;Goulson, Lye, & Darvill, 2008a). Within the Olympus National Park protected area, And because climate change is considered as the main cause of bumblebee decline by most recent studies (Kerr et al., 2015;Soroye, Newbold, & Kerr, 2020), we justifiably expect the impact of climate change on Mt. Olympus to be significant.
Communities are more vulnerable to loss of generalists, and bumblebees are generalist pollinators offering pollination services to numerous plant species (Memmott, Waser, & Price, 2004). This implies that the consequences in communities will be critical in case of bumblebee species extinction. We assume, however, global warming effects to vary among communities along the altitudinal gradient, not only because of the existence of microrefugia hosted within a highly heterogeneous microclimatic environment (Suggitt et al., 2018), but also because of the communities' diverse resilience related to a variety of different factors. Such a factor, for instance, is nestedness, which was found to increase community's resilience (the speed at which the community returns to the equilibrium following a disturbance) (Dalsgaard et al., 2013;Thébault & Fontaine, 2010), suggesting that in our case communities at low and high altitudes may be more resilient against global warming impacts (see pollination networks below).
Bumblebees respond to global warming by either restricting their southern (Kerr et al., 2015) or shifting uphill their lower (Ploquin et al., 2013)  Our results indicate that the altitudes between 1,900 and 2,200 m a. s. l. are the most important for conserving the highest diversity of bumblebee species on Mt. Olympus. Such a predominance of high-elevational peak in bumblebee diversity agrees with other studies published so far. In some of them, bumblebee diversity was found to peak at relatively lower altitudes than on Mt. Olympus (at 900-1100 m, Gorce and Tatra mountains, maximum study altitude 1,580 m high, Poland [Goulson et al., 2008b]; at 1,300-1,900 m, Mt. Norikura, maximum study altitude 2,600 m high, Japan [Egawa & Itino, 2019]). But Williams (1991) found bumblebee diversity to occur even higher, between 3,000 and 3,800 m (Mt. Apharwat, maximum study altitude 4,143 m, Kashmir Anagnostopoulos (2005,2009,2017) Turkey ( (Adedoja et al., 2018;Hoiss et al., 2015)

| Pollination networks
We found bumblebee-plant interaction networks on Mt. Olympus to be larger, more diverse, and more generalized at intermediate altitudes (1,500-1,800 m), which is also reflected in the lowest nestedness values (NODF) found for the same altitudes ( Figure 5, Table S5). This elevation range is higher than the range of plant diversity maxima (viz. 500-1,500 m) and lower than that of bumblebee diversity maxima (viz. 1,900-2,200 m) (cf. Figures 3a,b and 4; Table S4). Evidently, the resulting larger, more diverse and more gen-

| CON CLUS IONS
We found bumblebee community structure on Mt Olympus to be Certainly, future research will fill this gap in order to advise us on the key species, specific measures and management schemes including assisted migration, on Mt. Olympus (Hällfors et al., 2018).
Overall, this first systematic study regarding the bumblebees of Mt. Olympus highlights the need for future research in order to deepen our understanding of the diversity and particularities of pollinators and their host plants of this divine mountain. This is to enhance our knowledge regarding not only the present and past status of its biota, including the effects of Pleistocene glaciations, but also how to treat habitat conservation on this mountain in view of global warming.

ACK N OWLED G EM ENTS
We thank the Hellenic Rescue Team, especially N. Parmakis and T. Egiptiadis and "Petrostrouga" refuge staff for providing accommodation facilities; Ch. Tsantilis and Zagori Rainbow Waters S.A.
for supplying plastic containers to transport propylene glycol; M. Olympus who facilitated the completion of this study in various ways.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13138.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in the