Urban fragmentation leads to lower floral diversity, with knock-on impacts on bee biodiversity

Bees and flowering plants are two closely interacting groups of organisms. Habitat loss and fragmentation associated with urbanisation are major threats to both partners. Yet how and why bee and floral richness and diversity co-vary within the urban landscape remain unclear. Here, we sampled bees and flowering plants in urban green spaces to investigate how bee and flowering plant species richness, their phylogenetic diversity and pollination-relevant functional trait diversity influence each other in response to urban fragmentation. As expected, bee abundance and richness were positively related to flowering plant richness, with bee body size (but not bee richness and diversity) increasing with nectar-holder depth of flowering plants. Causal modelling indicated that bottom-up effects dictated patterns of bee-flower relationships, with urban fragmentation diminishing flowering plants richness and thereby indirectly reducing bee species richness and abundance. The close relationship between bees and flowering plants highlights the risks of their parallel declines in response to land-use change within the urban landscape.

monthly between June and August 2017 using transect sampling in eight urban semi-natural sites. From a total of 1440 min of transect sampling, we collected a total of 845 wild bee specimens representing 63 species from 23 genera within six families (Supplementary Table S1). The majority of these specimens (326, 39%) were bumblebees (Bombus spp.), 140 (17%) were from the genus Halictus, 108 (13%) were from the genus Lasioglossum, 106 (13%) were from the genus Hylaeus, 36 (4%) were from the genus Megachile, 26 (3%) were from the genus Dasypoda, 25 (3%) were from the genus Panurgus, 22 (3%) were from the genus Andrena, and 20 (2%) were from the genus Colletes. The number of bee species sampled per sampling site and per transect walk ranged from 6 to 14 (mean ± SD = 9.91 ± 2.60). Additionally, from our transect sampling we collected 291 honey bee individuals. Due to the potential effect of honey bees on wild bee communities 44 , the abundance of honey bees per site was used as a predictor in all downstream analyses. Honey bee abundance was not, however, an important predictor in any of our statistical models.
The majority of bee species sampled were ground nesters ( We sampled a total of 58 flowering plant species in flower (Supplementary Table S2). The most widespread were Hypericum perforatum, Daucus carota, Rubus fruticosus, Hieracium lachenalii and Tanacetum vulgare. The number of plant species in flower per sampling site and per transect walk ranged from 5 to 17 (mean ± SD = 10.92 ± 3.19). Most had an open flower shape (36,62.06%), followed by plants with papilionaceous (19,32.75%) and tubular flower shape (3,5.17%). The majority of the plants sampled had a yellow flower colour (22,37.93%), followed by violet (17,29.31%), white (16,27.58%) and red (3,5.17%). Furthermore, most of the plants were perennials (38,65.51%) and allogamous (38,65.51%) with protandrous (30,51.72%) or homogamous (22,37.93%) flower sex timing. Local (patch) and landscape scale effects on bee communities. We calculated flowering plant and bee diversity metrics (i.e. abundance, richness, functional diversity and phylogenetic diversity) and quantified local and landscape environmental variables (i.e. bare soil, landscape composition and fragmentation). We first Scientific Reports | (2020) 10:21756 | https://doi.org/10.1038/s41598-020-78736-x www.nature.com/scientificreports/ explored the main flowering plant community diversity metrics and environmental drivers of bee diversity. We found a variety of relationships between bee communities and flowering plant diversity, bare soil cover and landscape heterogeneity. Firstly, bee species richness was positively related to flowering plant richness (Table 1; Fig. 1a). Overall bee abundance was also positively associated with flowering plant richness and the percentage of bare soil cover (Table 1; Fig. 1b,c). Bee functional diversity decreased with increasing community weighted  Fig. 1d). The CWM of bee intertegular distance (ITD) was, though, positively related to the CWM of plant nectar holder depth (Table 1; Fig. 2) as well as local flower richness (Table 1). Therefore diverse flowers equated with diverse bees. Yet as mean corolla length increased, bee communities became less diverse and with a larger average body size. We did not find any important predictor for bee phylogenetic diversity (Table 1). Landscape composition, edge density, habitat fragmentation, overall plant phylogenetic and functional trait diversity were not included in the most parsimonious models of overall bee abundance and richness, bee functional and phylogenetic diversity. As bee species responses to the availability and distribution of resources could be trait-specific, we used fourth-corner analysis and explored the relationships between the abundance of each functional group and local flowering plant diversity, bare soil cover and landscape scale factors. Due to the strong correlation among bee traits (i.e. polylectic bees were mostly social, long-tongued bees were usually large, and above-ground nesters were mostly solitary and with a long tongue (Supplementary Table S3), we ran fourth-corner analysis only for lecty, nesting and body size. Fourth-corner analysis revealed numerous associations between the abundance of different bee functional groups, floral and environmental variables (model deviance = 48.03, P = 0.004; Fig. 3). Oligolectic bee species abundance was positively related to the proportion of residential cover, including domestic gardens (Fig. 3). In contrast, oligolectic bee species abundance was negatively related to the CWM of nectar holder depth (Fig. 3). Polylectic bee species abundance was positively associated with local flowering plant richness, the proportion of allotment gardens, and the proportion of parks (Fig. 3). Above ground nesting bee abundance was positively associated with the proportion of residential cover (domestic housing with gardens) (Fig. 3). Ground nesting bee abundance was positively associated with local flowering plant richness and the percentage of bare soil cover (Fig. 3). Bee body size (ITD) was also positively associated with local flowering plant richness and the CWM of nectar holder depth (Fig. 3).  www.nature.com/scientificreports/ Bee and landscape scale effects on flowering plant communities. We also explored the main bee community diversity metrics and environmental drivers of plant communities. Plant species richness was positively related to bee abundance and negatively related to landscape fragmentation (Table 1; Fig. 4a,b). Flowering plant functional diversity was also lower in more fragmented landscapes (Table 1; Fig. 4c). Plant phylogenetic diversity was also negatively related to urban landscape fragmentation (Table 1; Fig. 4d). As expected, the community level weighted mean of nectar holder depth was positively associated with the CWM of bee intertegular distance (Table 1). Landscape composition, edge density, overall bee phylogenetic and functional trait diversity were not included in the most parsimonious models of flowering plant richness, functional and phylogenetic diversity.

Structural equation modeling.
To reveal those putative causal factors, or interactions among them, that influenced plant and bee communities and to infer causality in plant-bee relationships within the urban ecosystem, we used structural equation modeling (SEM). Our piecewise SEM method confirmed the strong relationships between bee biodiversity, local ground nesting resources, floral richness and urban fragmentation. In the best version of the SEM, with a stable fit to our data (AIC plants→bees = 33.387, Fisher's C = 9.387, d.f. = 8, P = 0.311; Fig. 5), flowering plant richness had a positive impact on overall bee richness (P = 0.02; Fig. 5). In addition,  www.nature.com/scientificreports/ bee abundance was driven by both flowering plant richness and bare soil cover (P < 0.01; Fig. 5). Flowering plant richness itself was negatively impacted by urban habitat fragmentation (P = 0.02, Fig. 5). The relationship between bee richness and abundance and flowering plant richness resulted in an indirect negative effect of habitat fragmentation on bees (Sobel test; − 0.06, P = 0.04 for bee richness and − 0.09, P = 0.02 for bee abundance) mediated by its impact on flowering plant richness. The alternative SEM with the direction of relationship from bees to flowering plants had a higher AIC and considerably reduced support (AIC bees→plants = 38.328 versus AIC plants→bees = 33.387; AIC bees→plants − AIC plants→bees = 4.941).

Discussion
Our sampling in urban semi-natural sites revealed strong relationships between flowering plant and bee diversity. We furthermore found urban fragmentation to be negatively related to flowering plant species richness, functional and phylogenetic diversity. In accordance to our expectations, our causal modeling identified bottom-up effects (from flowering plants to bees) to dominate the relationship between flowering plant and bee diversity. Most importantly, our causal modeling revealed an indirect negative effect of urban fragmentation on bee diversity that was mediated via its negative effects on flowering plant richness. We expand on these results and discuss their ecological implications. We found marked relationships between overall bee richness and local flowering plant richness and between overall bee abundance, flowering plant richness and bare soil cover. This reflects the strong dependence of bees on the local availability not only of food (the richness of flowering plants) but also of nesting resources, namely the availability of bare soil cover for below-ground nesting bees, as has been previously reported for semi-natural and urban ecosystems 10,11,13,14,32,[45][46][47] . Additionally, our results are in line with previous studies suggesting that bees might respond positively to local floral and nesting resources, irrespective of surrounding land-use change 14 .
Resource availability is a strong determinant of community assembly, and bee species may respond differently depending on how resources are distributed at different scales 48 . Investigating patterns of bee biodiversity across broad scales could potentially mask trends in bee functional guilds at smaller scales, knowledge of which could be important for bee conservation actions within the urban environment. Bees are "central place foragers" (i.e. females carry resources back to their nest); as such, their resource needs must be met within their flight range limits. Large-bodied species with large foraging ranges might respond to the resource distribution not only at small but also at larger scales. In contrast, small-bodied species with short home ranges might only be able to respond to resource availability locally. Furthermore, species that depend on two or more resources (e.g. nesting and food resources) that vary in their distribution patterns might respond to each resource at a different scale. In addition to the availability of local floral and nesting resources, our data also showed that many of our sampled species also responded to the availability of resources at larger scales in a trait-dependent manner. Several of the surrounding urban green land-uses had a positive effect on the abundance of many of the bee functional groups sampled. Our results suggest that parks, allotment gardens and residential areas can all be important habitats for bees 49 and point to the importance of the availability of resources (i.e. food and nesting) at both the local and the landscape scales for bee communities within the urban ecosystem.
We found that bee responses to urbanisation were dependent on species life history traits, especially those related to foraging and nesting preferences. Previous studies have documented that urban areas benefit cavity-(above ground) nesting and generalist bee species 26,27,[32][33][34] . Yet, in our study we found a greater abundance of ground nesting bees in urban areas. This suggests that not only urban residential, public and private gardens 33,49,50 but that low management fragments of semi-natural vegetation could also be important habitats within the urban ecosystem which benefit ground nesting bees. Furthermore, we found more polylectic compared to oligolectic bee species. Although flowering plant richness could be high in urban areas 37 , it seems that it may only benefit generalist foragers 14 .
Urban landscape structure did not appear to relate to bee functional or phylogenetic diversity. Indeed, the CWM of flowering plant nectar holder depth was the only significant predictor of bee functional diversity, which decreased as the CWM of flowering plant nectar holder depth increased. Functional traits related to size are important predictors of mutualistic networks and can act as barriers preventing interactions 22,51 . The increase in the CWM nectar holder depth seems to act as a major filter on the wild bee community, leading to trait convergence favouring long-tongued, large and social bee species, traits that we also found to be strongly correlated among bee species. The lack of a relationship between overall flowering plant functional diversity and bee functional diversity and the strong relationship between the CWM of bee intertegular distance and the CWM of flowering plant nectar holder depth suggest that the majority of the flowering plant traits analysed were unimportant for structuring bee communities. They suggest that size constraints have the strongest effect on bee and floral community assembly.
During urbanisation, semi-natural areas are converted into those with impervious surfaces, resulting in habitat loss 3 . Furthermore, the development of roads, railways and impervious surfaces results in habitat fragmentation 52 . Urban habitat loss and fragmentation reduce the size and increase the isolation of plant populations, which makes species more vulnerable to extinction 53 . In contrast to the bees, we found a strong negative and direct impact of urban fragmentation on flowering plant richness. Moreover, fragmentation resulted in flowering plant functional convergence and phylogenetic clustering.
Although we did not find direct effects of fragmentation on bee biodiversity, our findings suggest plant-mediated effects of urban fragmentation on bees. Bees are well adapted to shift between patchy foraging and nesting resources and often benefit from habitat edges between patch types 28,54 . Nonetheless, depletion and extensive fragmentation of floral food resources at local and landscape scales could lead to local bee population extinctions and limit recolonizations, with negative effects on bee abundance and richness 55  www.nature.com/scientificreports/ reduced dependency on pollinators than pollinators do on plants, potentially due to the alternative reproduction strategies of many plants such as clonal propagation, self-fertilisation and seed banking. Our results point to the importance of considering indirect effects mediated by species interactions when assessing species' responses to land-use changes, including that of urbanisation, and further highlight the risks of parallel declines between wild bees and flowering plants. From an applied perspective, our results have important implications for pollinator conservation in cities. The strong effects of flower resource richness on bee communities suggest that simple local management practices to increase the diversity of native plants could benefit wild bees in urbanised areas.

Materials and methods
Study sites. Our  Using land-cover maps obtained from Geofabrik GmbH (http://www.geofa brik.de/) and Quantum GIS 56 , we pre-selected eight independent sites with semi-natural vegetation that differed in their surrounding proportion of urban cover (i.e. residential cover and amount of impervious surfaces) (Supplementary Fig. 1; Supplementary  Table S4). All sites were of low management, without obvious differences in shading or sources of pollution. No mowing was performed prior to or during our sampling in any of the sites. We ensured a minimum distance of 1.5 km between sites, sufficient for them to be considered independent 57 .
Sampling flowering plant and bee diversity. Local communities of flowering plants and bees were sampled between 10th June and 18th August 2017 using 30 min transect walks at each of the eight sites. In every month (June, July and August), we conducted two rounds of sampling in each site, one round in the morning between 9:00 and 12:00 and one round in the afternoon between 13:00 and 16:00. Hence, every site was visited a total of six times. For the transect sampling we used the variable transect walk method, as described by Westphal et al. 58 . In each of the urban semi-natural sites, the collector was not restricted to a fixed transect line; instead they walked for 30 min at a slow steady pace among any potential floral resource patch and collected bees visiting flowers as well as the flower that was visited. In comparison to fixed transect walks, this method overcomes potential undersampling due to the spatiotemporal variation in bees and floral resources 58 . Furthermore, in a comparative study across geographic regions in Europe, variable transects walks were found to be more efficient with respect to sample coverage, bee species richness and abundance in comparison to fixed transect walks. We performed our transect walks only on sunny days with clear skies, wind speed less than 7 m/s and a temperature above 17 °C. Sampling was performed by the same collector. Bees were collected in tubes with 96% ethanol and stored at − 20 °C for later identification. We collected whole botanical specimens, including leaves and reproductive structures of the flowering plants. We identified bees to species using identification keys [59][60][61] and DNA barcoding of the COI gene 62 (Supplementary Table S1; Supplementary Methods 1). All flowering plants were identified to species using identification keys for the local flora 63,64 (Supplementary Table S2). Bee specimens are kept at the General Zoology research group, Martin-Luther University Halle-Wittenberg. Sequences obtained from COI barcoding were submitted to the NCBI GenBank database and are available under Accession Numbers MW065567-MW065776.
Flowering plant and bee diversity metrics. We calculated four metrics of community diversity: abundance, species richness, functional diversity and phylogenetic diversity. Abundance was calculated as the total number of every bee individual sampled in each site. Bee and flowering plant species richness were calculated as the number of species of each taxon in each site. To estimate functional diversity, for every individual bee sampled we quantified body size as intertegular distance (ITD) to the nearest 10 μm using an OLYMPUS SZX7 stereo-microscope and the software cellSens standard v.1.6 (Tokyo, Japan). In addition, using the available literature 65 , we obtained five bee functional traits for each species (sociality, tongue length, nesting behaviour, voltinism and lecty; Table 2).
For every flower collected, we measured the nectar holder depth (the maximal depth between the base of flower, where nectar is stored, and where the flower opens and is accessible to all insects) using a Vernier calliper to the nearest 10 μm. Additionally, we used two qualitative measures to characterise every flower: (1) flower restrictiveness (the accessibility of a flower for a plant-flower visitor due to its shape) and (2) flower colour. Flower restrictiveness was classified into three categories: (i) open or bowl-shaped flowers, (ii) bell shaped, short tubular or long narrow tubular flowers, and (iii) papilionaceous flowers. Using the TRY database v. 4.1 66 (see Supplementary Table S2 for detailed references) we extracted three additional pollination-related plant functional traits: breeding system, flower sex timing and plant longevity ( Table 2).
As a measure of functional diversity we used Rao's quadratic entropy (Rao's Q) 67,68 , calculated using the R package FD 68 . Community weighted mean (CWM) trait values for intertegular distance and nectar holder depth were also calculated using the R package FD. As functional diversity is partly confounded by species richness, we used a null model approach to obtain an unbiased standardised functional diversity estimate per site 69 . For a given number of species, this approach creates a random (null) distribution of functional diversity values. Holding species richness constant for each site, it uses randomly selected species from the entire species pool to calculate a null functional diversity for each richness level. We repeated the procedure 1000 times to produce a distribution of null values. The standardised Rao's Q was calculated for each site as; (N − N̅ r )/σ Nr , where N is the observed value of Rao's Q and N̅ r and σ Nr are the mean and standard deviation, respectively, for the 1000 replicated randomised (null) Rao's Q values. Standardised Rao's Q values > 0 suggest that traits are more divergent To estimate phylogenetic diversity we build bee and plant phylogenies based on two genes for each taxon (bee: COI and elongation factor 1-alpha [ef1-alpha]; plant: ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit [rbcL] and maturase K [matK]), which were downloaded from NCBI GenBank (see Supplementary Methods 2). To estimate the total evolutionary history shared across all species within a community (i.e. phylogenetic diversity), we used Faith's phylogenetic diversity index (PD; 70 ). PD is the phylogenetic analogue of taxon richness and is defined as the minimum total length of all the phylogenetic branches required to span a given set of taxa on the phylogenetic tree 70 . Similarly to functional diversity, PD and species richness are strongly correlated and we thus used the same null model approach used for functional diversity to obtain standardize PD measurements. The standardized PD was calculated by dividing the difference between the observed and expected PD by the standard deviation of the null distribution. Standardised PD values > 0 indicate phylogenetic dispersion and values < 0 indicate phylogenetic clustering.

Local (patch) and landscape variables.
To quantify local (patch) environmental variables, we randomly placed ten 1 m 2 quadrats in every site and on every sampling visit. Within each quadrat, we measured the percentage of exposed soil, quantified as a surrogate for the availability of nesting opportunities for ground nesting bees like Andrenidae, Colletidae and Halictidae 45 . Additionally, we quantified patch size (i.e. the size of the continuous area of green or semi-natural vegetation at our sampling site) for every site using Google Earth Pro v.7.3.0 (Supplementary Table S4).
At each sampling site we estimated habitat heterogeneity at six spatial scales (250, 500, 750, 1000, 1250, and 1500 m) using Quantum GIS 56 and land cover data obtained from Geofabrik GmbH (http://www.geofa brik.de/). The accuracy of the data was confirmed using onsite surveys and by matching the shapefile's features with Google earth images. We calculated landscape diversity (H s ) for each site at each radius as: H s = − ∑p i × ln p i , where p i is the proportion of each land cover of type i.
We quantified landscape heterogeneity with a number of metrics known to impact flowering plant-bee interactions 14,49 . These metrics were: (i) the proportion of parks (ii) the proportion of forest cover, (iii) the proportion of allotment gardens (iv) the proportion of semi-natural areas (meadows, vacant lots, remnants), (v) the extent of agricultural cover, (vi) the proportion of residential cover (domestic housing with gardens), (vii) edge density, as the total length of green patch edges divided by the total area, which quantifies landscape configuration, and (viii) number of disconnected green patches divided by their total surface area as an indicator of the degree of fragmentation.
To detect the spatial scale at which land cover had the most power to explain flowering plant and bee occurrence, we correlated flowering plant and bee diversity with landscape Shannon diversity in each of our eight sites at all six spatial scales and compared the resulting correlation coefficients. Correlation coefficients peaked at 1000 m radius for both interacting partners (Supplementary Table S5), which was then chosen as the spatial scale for subsequent statistical modelling. This scale is within the foraging range of many bee species 57 .
Statistical analyses. We tested for local (patch) and landscape scale effects on each bee and plant diversity metrics using linear mixed effects models (LMMs) and generalised linear mixed models (GLMMs). We used LMMs for functional diversity and phylogenetic diversity and GLMMs with a Poisson error structure for abundance and richness. Plant biodiversity metrics (i.e. richness, functional and phylogenetic diversity) were included as predictors to the models fitted to bee metrics and vice versa. Bee abundance was included as a covariate in the models for bee and plant species richness to control for sample size effects. In all models, month and Table 2. Description of (a) bee functional traits and (b) flowering plant functional traits used in this study. www.nature.com/scientificreports/ site were included as random effect factors. We used an automated model selection approach (all subsets) based on the second-order Akaike Information Criterion corrected for small sample size (AICc), and allowing only up to three variables to avoid over-fitting, to identify relevant predictors for each biodiversity metric using the 'dredge' function in the R package MuMIn 71 . We used a cut-off ΔΑΙCc value of 2 72 and, if more than one model was retained, we used model averaging (function 'model.avg'; 71 ). All mixed model analyses were performed using the package lme4 73 .
To explore relationships between bee functional traits and environmental variables (flowering plant diversity metrics, percentage of exposed soil and landscape heterogeneity) we used fourth-corner analysis 74 using the 'traitglm' function within the R package mvabund 75 . Multivariate generalised linear fourth-corner models were fitted with a negative binomial distribution and a least absolute shrinkage and selection operator's (LASSO) penalty (i.e. method = glm1path). LASSO is an approach which penalises coefficients that do not reduce Bayesian Information Criteria (BIC) to zero. Analysis of model deviance was estimated using a Monte-Carlo resampling procedure (1000 resamples) to evaluate the global significance of trait-environment relationships.
In addition to our multiple regressions models, we performed piecewise structural equation modelling (SEM) in an attempt to infer causality in plant-bee biodiversity relationships. We hypothesised that local and landscape composition and configuration might have affected bee and flowering plant biodiversity directly and also indirectly through affecting each trophic level. We performed the piecewise SEM analysis using the R package piecewiseSEM 76 . Based on a priori knowledge of interactions and our multi-level modelling, we used the d-separation (d-sep) test to evaluate whether the non-hypothesised independent paths were significant and whether the models could be improved with the inclusion of any of the missing path(s). We used Fisher's C statistic for evaluating the fit of piecewise SEM 77 . We constructed alternative structural equation models, changing the direction of causality from bees to flowering plants to that of flowering plants to bees. Models were compared using the Akaike information criterion (AIC) and Δ i . Δ i is the difference between the AIC value of a given model against the value of the model with the lowest AIC min (Δ i = AIC i − AIC min ). Models with a Δ i > 4 have little support 72 . Path coefficients and deviance explained were then calculated for the best structural equation model. We used Sobel's method to test for significant indirect effects 78 .
Prior to all analyses we standardized all quantitative predictors. This was done to minimise potential effects of collinearity and to derive comparable estimates. Additionally we used variance inflation factors (VIFs) to check for collinearity among our explanatory variables. VIF was lower than 3 for all predictors, indicating no major effects of collinearity 79 . The residuals of all models were tested for spatial autocorrelation using Moran's I, implemented in the R package ape 80 . Residuals were not found to be autocorrelated (P > 0.05 for all models). All statistical analyses were performed in R v. 3.5.2 81 . All model (GLMMs and LMMs) assumptions were checked visually and were found to conform to expectations (e.g. normality of the distribution of residuals, homogeneity of variances, linearity, non-overdispersion).

Data availability
All data are included as supplementary material.