Narrow habitat breadth and late-summer emergence increases extinction vulnerability in Central European bees

1. Introduction

Bee diversity in the USA and Europe appears to be declining due to intensive agriculture, habitat loss, and the increased use of herbicides and insecticides, all of which negatively affect food sources and nesting opportunities [1–4]. Climate warming per se should benefit, not harm, bee diversity because most bees do well in Mediterranean-type climates [5], with the exception of a few cool-adapted groups, such as bumblebees [6]. A recent study of bee diversity in a protected urban garden in Germany indeed supported an increase in warm-adapted species with climate warming, based on inventories taken 20 years apart [7]. In non-protected areas, however, reduced floral and nesting resources due to human land use negatively impact bee diversity and abundance [1,3,8,9], potentially masking the expected positive effect of shorter winters and warmer springs and summers.

The effect of reduced and temporally shifted floral resources (because of earlier flowering under climate change) on the decline of wild bees in Europe is evident from a study of population trends in 57 species (from 10 genera) in the Netherlands, which found that late-flowering pollen and nectar sources have declined more strongly than early-flowering ones, leading to stronger declines in late-flying bees [10]. Earlier studies on bumblebee decline also implicate changes in agricultural practices. Thus, in Ireland, late-nesting bumblebee species associated with grassy habitats are declining, probably due to the agricultural trend of replacing hay with silage, leading to earlier mowing [11]. Similarly, bumblebee species in Britain, Canada, and China with late-active queens are especially susceptible to decline [12]. All this differentially affects the relative extinction risks of early- and late-flying bees. In addition, climate warming is uneven among seasons [13,14]. For instance, climate warming from 1960 to 2016 across 53 German weather stations shows strong March–May and July–August warming, while June, September, and October have hardly changed over those 56 years (inset figure 2a). If spring-warming and the early mass-flowering of herbaceous crops [10] benefit spring-flying bees, while agricultural activities during summer and autumn are detrimental to late-flying bees, this might lead to a shift in the relative vulnerability of spring- versus summer-active species.

Here, we use the German bee fauna as a system in which to study extinction risk predictors. Germany harbours 561 bee species in 35 genera [15], most of which have been barcoded with multiple specimens per species [16]. In terms of both taxonomy and ecology, Germany's bees are comparatively well known, with red-listing having been carried out for at least 40 years [15,17]. We, therefore, decided to rely on Red List status and whether status assessments have changed since 1977 as a proxy for species' extinction risks. Of the 445 bees included in this study, about 23% are pollen specialists (oligolectic), 227 (51%) build their nests in the ground and flight times fall between late February (first species of Andrena and Osmia, early-emerging overwintering queens of Bombus) and mid-October (e.g. species of Colletes and freshly emerged sexuals of Lasioglossum and Halictus).

Specifically, we set out to investigate whether early-season and late-season species differ in their Red List status and how their abundances and long-term population trends have changed over the past 40 years. Our expectation was that bees active in the fast-warming spring might do better than those later in the year. In addition, we were interested in how bees of different habitats may have fared over the past half-century. About 51.6% of Germany's land surface consists of arable land [18], 14% of urban areas [18], and 32% of forests, percentages that have changed little over the past 20 years, except for a slight increase in the percentage of urban area [19] (see electronic supplementary material, figure S1: land use 2000–2017). If urbanization and climate warming are less of a threat to bees than is intensive agriculture, this might lead to bees in urban or forest habitats doing better than those in open, non-urban habitats. We used hierarchical Bayesian (HB) models, which allowed us to control for phylogenetic structure in the data, to analyse the simultaneous effects of the biotic and abiotic traits relevant to our questions, namely habitat breadth (using a fine-grained categorization that included all types of vegetation mentioned in the autecological literature on the bees in question, namely meadows, forests, hedgerows, urban areas, raw-soil sites, and wasteland as well as altitudinal distribution), pollen specialization, body size, nesting preference, mean flight time, and duration of flight activity.

2. Material and methods

(a) Phylogeny

Cytochrome oxidase subunit I sequences of most bees that occur in Germany are available from the barcoding study of Schmidt et al. [16]. Where necessary, we downloaded additional sequences from GenBank in April 2018 (electronic supplementary material, table S3), resulting in a matrix of 553 species (TreeBase accession 23291). Maximum-likelihood tree inference relied on RAxML v. 8.0 [20]. Phylogenetic signal in species-level extinction risk and habitat breadth was measured with Pagel's λ [21], using the ‘phylosig’ function in the R package ‘phytools’ v. 0.2-1 [22] (figure 1).

• Download figure
• Open in new tab
• Download PowerPoint

(b) Trait scoring

We obtained information on life-history traits and extinction vulnerability status for 445 out of 561 German bee species, thus covering 79% of the German bee fauna [15] (electronic supplementary material, table S1). For each species, we scored its current Red List status either as an ordinal vector (1, not threatened; 2, near-threatened; 3, threatened; 4, highly threatened; 5, almost extinct; 6, extinct; figure 2a) or as a binary character (0, not threatened; 1, threatened; electronic supplementary material, figure S2) as well as abundance and long-term population trends from Westrich et al. [15]. The Red List status was determined based on the data of seven bee specialists, each with a different regional focus: P. Westrich for Baden-Wuerttemberg, U. Frommer for Hesse, K. Mandery for Northern Bavaria, H. Riemann for Bremen and Lower Saxony, H. Ruhnke for Saxonx-Anhalt, C. Saure for Berlin and Brandenburg, and J. Voith for Southern Bavaria and the German Alps. In terms of Red List status, 16 species are considered extinct, 26 almost extinct, 61 highly threatened, 100 threatened, 38 near-threatened, and 187 not threatened. For 17 species, the data were insufficient to assess their Red List status, leaving 428 red-listed species in our statistical models. Information on current commonness status from the Red Lists was available for all 445 species, long-term population trends, based on data of the past 50–150 years, were known for 404 species (table 1).

• Download figure
• Open in new tab
• Download PowerPoint

We scored species’ pollen preferences (pollen specialization) (polylectic, oligolectic, or cuckoo) based on Westrich [23] and Scheuchl & Willner [24]. In a few cases, data were complemented with information from the International Union for the Conservation of Nature (IUCN) Red List [25]. Mesolectic species, which use pollen from a few genera of plants, were included in the category polylectic; kleptoparasitic and social parasitic species were placed in the category cuckoo. In total, 218 species were assigned the category ‘polylectic’, 102 ‘oligolectic’, and 125 ‘cuckoo’. Species' nest location (above or below ground) was scored based on the literature [23,24], with 213 species nesting solely below ground and 232 above or below ground or in host nests. Sociality was coded as a binary trait, and all species not explicitly stated to be social in Westrich [23] or Scheuchl & Willner [24] were considered solitary, regardless of whether they are parasitic, form nest aggregations, or nest alone (268 species were social, 178 solitary).

The mean body length of females (as a proxy for body size) was recorded in millimetres based on Amiet et al. [26–30], von Hagen & Aichhorn [31], Schmid-Egger & Scheuchl [32], Dathe & Scheuchl [33], and Scheuchl [34,35]. For bumblebees, the size of female workers, not queens, was used for the analysis.

The average duration of flight activity per year (duration of flight activity; see electronic supplementary material, table S2) and the average month of flight activity (mean flight month) were scored based on Westrich [23] and Scheuchl & Willner [24] who have for many years monitored relative bee abundances and flight times in southwestern Germany. To calculate the duration of activity for bivoltine species, we used the cumulative span of occurrence, summing the duration of activity of the first and the second generation. To obtain the mean flight month, we calculated the mean of the first and last month of activity.

Habitat preferences were taken from Westrich [23] and Scheuchl & Willner [24] with six categories: forests and heaths, including alluvial forests (forests); raw-soil sites with little vegetation, such as sand dunes, heathland, steppes, and sand or gravel pits; urban areas, including gardens and parks; meadows (mown at least once a year); hedgerows; wasteland and nutrient-poor sites. Species occurring everywhere or with no specific habitat preferences were scored as ubiquists. In a second step, we scored habitat breadth as the number of habitats a species is found in (thus ranging from 1 to 6). Ubiquists were assigned a habitat breadth value of 6. Species' altitudinal distribution (altitude) was scored as lowlands from 0 to 300 m, colline from 300 to 800 m, montane from 800 to 1600 m, subalpine from 1600 to 2100 m, and alpine from 2100 to 3000 m [23,24].

(c) Hierarchical Bayesian modelling

We applied an HB approach (following [36]) for testing effects of species-level traits on extinction risk. This allowed us to control for phylogenetic structure in the data and to simultaneously fit slope parameters of biotic and abiotic traits relevant to our questions without concerns of p-value correction or multiple testing. To test for multicollinearity of our 14 predictor variables, we used a variance inflation factor (VIF) analysis, implemented in the R function ‘vif’, from the package ‘HH’ [37]. This analysis showed high variance inflation (greater than 10) for habitat breadth, because the six habitat types (urban areas, hedgerows, forests, raw-soil sites, meadows, and wasteland) together explain a significant proportion of habitat breadth. Thus, we decided to analyse specific habitat preferences in a separate analysis. After having separated habitat breadth and habitat preferences, all VIF were smaller than 2, indicating sufficient independence among predictor variables. We also calculated a correlation dendrogram, showing that Pearson's correlation coefficients among predictor variables are less than 0.5 (electronic supplementary material, figure S3).

We applied both binary and ordinal models, treating Red List status as a binary character (see above; figure 2a). To account for phylogenetic structure in the data, genus- and family-level random effects were incorporated in the models (note that the phylogeny could not be incorporated as a distance matrix in the analysis because this is only possible for continuous, not binary or ordinal, dependent variables). Six continuous and eight binary variables were used as predictor variables (see Trait scoring); continuous variables: body size, mean flight month, duration of activity, altitude, habitat breadth; binary variables: pollen specialization (oligolecty versus polylecty/cuckoo), sociality (social versus solitary), nest location (belowground versus aboveground ± belowground), urban areas (yes/no), hedgerows (yes/no), forests (yes/no), raw-soil sites (yes/no), meadows (yes/no), and wasteland (yes/no).

Regression components of the ordinal models are of the form:

Life-history traits and habitat breadth model:

The duration of activity, mean flight month, body size, pollen specialization, sociality, nest location, habitat breadth, urban areas, hedgerows, altitude, wasteland, meadows, raw-soil, and forest refer to species values (i) in electronic supplementary material, table S1, β refers to the estimated slopes of the respective variable (figure 2) and genus and family refer to the genus- and family-level random intercept effects inserted in the model.

Regression components of the binary models are of the form:

Life-history trait and habitat breadth model:

We additionally applied the same models, using current abundance status or long-term population trends instead of Red List status as response variables (electronic supplementary material, figures S1 and S2). Again, we applied both binary and ordinal models, treating abundance and population trend either as ordinal vectors or binary characters (abundance ordinal vector: 1, abundant; 2, common; 3, frequent; 4, occasional; 5, rare; 6, extremely rare; 7, extinct; abundance binary: 0, abundant to frequent; 1, rare to extinct; population trend ordinal vector: 1, increase; 2, no change; 3, small decline; 4, decline; 5, large decline; 6, extinct; population trend binary: 0, no decline; 1, decline).

To examine relative effect sizes of predictor variables, we standardized all variables by subtracting their mean and dividing by 2 s.d. before analysis [38]. The resulting posterior distributions are a direct statement of the relative effect of each tested trait on species-level extinction vulnerability. The effective posterior means for the respective relationships are shown in figure 1a; electronic supplementary material figures S2 and S3. To parametrize our models, we used the JAGS implementation [39] of Markov chain Monte Carlo (MCMC) methods in the R package R2JAGS [40]. The R code for our HB models is provided in the electronic supporting material as Note S1. We ran three parallel MCMC chains for 20 000 iterations with a 5000-iteration burn-in and evaluated model convergence with the Gelman and Rubin [41] statistic. Non-informative priors were specified for all parameter distributions, including normal priors for α and β coefficients (fixed effects; mean = 0; variance = 1000), and uniform priors between 0 and 100 for the variance of the genus random intercept effect, based on De Villemereuil et al. [42].

In table 1, we summarize the statistical results.

(d) Classification rates and recursive partitioning analysis

Multivariate generalized linear models were applied to test how well the combination of scored traits explains species-level extinction risk, i.e. testing for type 1 and type 2 errors (species erroneously considered as not endangered by the model or species erroneously considered endangered by the model). The models included all 14 predictor variables and an interaction for each pair of them to account for possible interactions among traits and Red List status, abundance, and population trend as binary dependent variables (table 1). Overall, Red List status was predicted with an accuracy of 85%; false positives (type 1 error) occurred for 9% of species and false negatives (type 2 error) for 6%. For recursive partitioning analyses (figure 3) of the four best predictor variables (habitat breadth, urban areas, duration of activity, and mean flight month), we used the ‘rpart’ option in the R library [43], setting the minimum node size to 15 (minimum number of species contained in each terminal node).

• Download figure
• Open in new tab
• Download PowerPoint

All statistical analyses were conducted in R 3.2.2 [44].

3. Results

Early-emerging bees had the lowest extinction probabilities, late-emerging bees the highest (figures 2 and 3), and extinction risk was also significantly increased in species with narrow habitat breadths (figures 1 and 2) or short flight periods, irrespective of whether extinction risk was coded as an ordinal (figure 2) or a binary character (electronic supplementary material, figure S2). Large-bodied species were slightly more vulnerable to extinction (figure 2a). Pollen preference, sociality, and nest location had no statistical effect on vulnerability. When we ran the analysis again, replacing the response variable ‘extinction risk’ with either ‘bee abundance’ or ‘long-term population trends’, the logistic regression (table 1) and HB models (compare figure 2a with electronic supplementary material, figure S4a,b) showed similar results, with similar effect sizes of predictor variables.

The fine-grained habitat scoring possible because of the relatively well-known autecology of Germany's bees revealed that occurrence in urban habitats, for example gardens and parks, as well as hedge habitats, strongly lowers a species' extinction risk compared with occurrence in any of the other five habitat types (figure 2a). High-altitude bees (coded as a continuous trait; Material and methods) also have a reduced extinction risk (figure 2a). Occurrence in forests and nutrient-poor sites, such as raw-soils, meadows, and wastelands had no statistical effect on vulnerability.

A logistic regression model, using the predictor variables (figure 2a) and pairwise interaction terms among them, correctly predicted the extinction risk of 85% of the 428 species for which we had information on Red List status (table 1). Nine per cent (38 species) were erroneously predicted as endangered (type 1 error) and 6% (26 species) were erroneously predicted as not endangered (type 2 error). Whether a species was misclassified or not did not show a phylogenetic signal (electronic supplementary material, figure S5).

4. Discussion

The three strongest extinction-predicting factors for the bees with data on both life-history traits and Red List status over the past 40 years (445 species or 79% of the German bee fauna) were narrow habitat preferences, a short flight time, and late-summer emergence (figures 1b and 3). By contrast, spring emergence reduced extinction risk. Shorter winters, earlier springs, and increased average monthly temperatures in Central Europe should benefit many species of bees because the superfamily Apoidea worldwide is most species-rich in Mediterranean-type climates [5]. For Germany, March–May and July–August have warmed more strongly than other months (inset figure 2a), a pattern that per se cannot explain why spring-active species should have a lower extinction risk than species active later in the year. The phenological mismatch between pollinators and flowers is also an unlikely explanation because current data indicate strongly buffered synchronization mechanisms between bees’ flight times and the flowering time of their pollen and nectar sources, probably because mutualistic interactions depend on such buffering for their persistence [45,46]. This leaves lack of pollen and nectar sources as the most probable explanation for why late-flying bees are declining more strongly than bees that reproduce early in the year when mass-flowering crops (e.g. Brassica napus), flowering trees and shrubs (e.g. Salix spp., Prunus spp.), and flowering herbs in improved grasslands and field margins (e.g. Taraxacum spp., Cardamine pratensis, and Ranunculus spp.) still provide early-season floral resources [10]. Intensive land use (with monocultures, fertilizers, insecticides, and herbicides) and habitat fragmentation, rather than changes in the relative area of arable or built-up land and forest, thus emerge as the most probable factors driving the decline of wild bee species in Germany.

A caveat applying to this study is that local extinction may be common in rare species, including species for which the area under study (Germany) represents the edge of their distribution ranges. Rare species might be less reliably monitored than more abundant species, so that their changing Red List status over the past 40 years may be less well understood than that of common species. However, models that included abundance and population trends as separate parameters yielded the same results (see electronic supplementary material, figure S4). Occurrence in urban areas lowered a species' extinction risk, under both the ordinal and the binary model (figure 2; electronic supplementary material, S2). Between 2000 and 2017, Germany's urban areas have been increasing in surface (electronic supplementary material, figure S1), and they support numerous bees that find nest sites, pollen, and nectar in city gardens, probably benefitting from lower pesticide use there compared to arable land [47,48]. However, the monitoring of changes in land use in Germany is at a coarse-scale, and it is therefore problematic to link them to bee decline.

That high-altitude bees have low extinction risks may be due to Germany's montane and alpine regions being less impacted by modern human land use than are its lowlands. The higher extinction risk of large-bodied species (figure 2a) agrees with the studies of Bartomeus et al. [46] and Scheper et al. [10] who both found that body size has a negative effect on wild bee population persistence, regardless of whether bumblebees were included or excluded in the analyses. A plausible explanation is that the larger pollen requirements of larger species [49] result in stronger food limitations compared to smaller bees. Larger bees also have larger foraging ranges [50–52], which may increase their exposure to pesticides.

An earlier study of traits affecting bees' extinction vulnerability that also incorporated phylogenetic relationships did not rely on Red Lists, but instead used relative abundances inferred from museum specimens (of 438 species in 47 genera) collected over 140 years in the northeastern USA [10]. The results revealed greatest declines (of collected specimens) in species with a narrow dietary breadth (pollen specialization), a short flight time (days of adult activity per year), and a large body size (intertegular distance in millimetres). Flight times were not included as a risk factor, so results cannot be compared with the present study. Meta-analyses of datasets collected either throughout the world ([8]: 19 studies) or throughout Europe ([3]: 30 studies) also did not include time of flight activity. The one study that did address a possible effect of time of flight activity on population trends—in 57 species from 10 genera in the Netherlands—found stronger declines in late-flying bees than early-flying ones [10] in agreement with our findings.

Overall, these results from a large sample of species from Central Europe that have been monitored for a comparatively long time highlight that intensive agriculture, with its negative repercussions for summer-active bees, is currently shifting faunistic diversity (at least for bees) towards warm-adapted, early-flying, city-dwelling species.

Data accessibility

All DNA sequences have been submitted to NCBI (https://www.ncbi.nlm.nih.gov) under the accession numbers listed in electronic supplementary material, table S3. All trait scoring is shown in electronic supplementary material, tables S1 and S2.

Authors' contributions

M.M.H. gathered data and observations; C.M.Z. conducted statistical analyses; S.S.R. and M.M.H. designed the study and wrote the first draft; all authors worked on the final manuscript.

Competing interests

The authors declare no competing financial interests.

Funding

Not applicable.

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4407731.