How wild bees find a way in European cities: Pollen metabarcoding unravels multiple feeding strategies and their effects on distribution patterns in four wild bee species

J Appl Ecol. 2021;00:1–14. | 1wileyonlinelibrary.com/journal/jpe Received: 23 July 2021 | Accepted: 12 October 2021 DOI: 10.1111/1365-2664.14063 RESEARCH ARTICLE How wild bees find a way in European cities: Pollen metabarcoding unravels multiple feeding strategies and their effects on distribution patterns in four wild bee species Joan CasanellesAbella1,2 | Stefanie Müller1,3 | Alexander Keller4 | Cristiana Aleixo5 | Marta Alós Orti6 | François Chiron7 | Nicolas Deguines7,8 | Tiit Hallikma6 | Lauri Laanisto6 | Pedro Pinho5 | Roeland Samson9 | Piotr Tryjanowski10 | Anskje Van Mensel9 | Loïc Pellissier2,11 | Marco Moretti1 This is an open access article under the terms of the Creat ive Commo ns Attri bution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Journal of Applied Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society Loïc Pellissier and Marco Moretti— Joint senior position. 1Biodiversity and Conservation Biology, Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, Switzerland 2Institute of Terrestrial Ecosystems, ETH Zurich, Zurich, Switzerland 3Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zurich, Switzerland 4Organismic and Cellular Interactions, Biocenter, Faculty of Biology, LudwigMaximiliansUniversität München, Martinsried, Germany 5Centre for Ecology, Evolution and Environmental Changes (cE3c), Faculdade de Ciências, Universidade de Lisboa, Lisboa, Portugal 6Institute of Agricultural and Environmental Sciences, Estonian University of Life Sciences, Tartu, Estonia 7Université ParisSaclay, CNRS, AgroParisTech, Ecologie Systématique Evolution, Orsay, France 8Laboratoire Ecologie et Biologie des Interactions, Equipe Ecologie Evolution Symbiose, Université de Poitiers, UMR CNRS, NouvelleAquitaine, France 9Laboratory of Environmental and Urban Ecology, Department of Bioscience Engineering, University of Antwerp, Antwerp, Belgium 10Department of Zoology, Poznan University of Life Sciences, Poznań, Poland Abstract 1. Urban ecosystems can sustain populations of wild bees, partly because of their rich native and exotic floral resources. A better understanding of the urban bee diet, particularly at the larval stage, is necessary to understand biotic interactions and feeding behaviour in urban ecosystems, and to promote bees by improving the management of urban floral resources. 2. We investigated the larval diet and distribution patterns of four solitary wild bee species with different diet specialization (i.e. Chelostoma florisomne, Osmia bicornis, Osmia cornuta and Hylaeus communis) along urban intensity gradients in five European cities (Antwerp, Paris, Poznan, Tartu and Zurich) using two complementary analyses. Specifically, using trapnests and pollen metabarcoding techniques, we characterized the species' larval diet, assessed diet consistency across cities and modelled the distribution of wild bees using species distribution models (SDMs). 3. Our results demonstrate that urban wild bees display different successful strategies to exploit existing urban floral resources: not only broad generalism (i.e. H. communis) but also intermediate generalism, with some degree of diet conservatism at the plant family or genus level (i.e. O. cornuta and O. bicornis), or even strict specialization on widely available urban pollen hosts (i.e. C. florisomne). Furthermore, we detected important diet variation in H. communis, with a switch from an herbaceous pollen diet to a tree pollen diet with increasing urban intensity. 4. Species distribution modelling indicated that wild bee distribution ranges inside urban ecosystems ultimately depend on their degree of specialization, and that broader diets result in less sensitivity to urban intensity. 2 | Journal of Applied Ecology CASANELLESABELLA Et AL. 1 | INTRODUCTION Wild bees are responsible for major ecosystem functions and make many contributions relevant to people, including pollination and maintenance of ecosystem stability, and they represent social and cultural values (e.g. Potts et al., 2016). Over the last decades, bee populations have dramatically decreased (Zattara & Aizen, 2021). Multiple causes have been identified (Goulson et al., 2015), and the loss of floral resources is among the most important (Goulson et al., 2015). Wild bees critically depend on large amounts of nectar and pollen to survive and reproduce during their life cycle, and most species display some degree of fidelity to specific plant taxa (Goulson, 1999; Vanderplanck et al., 2014). Consequently, diet specialization and diet preference are two key traits determining the sensitivity to landuse changes (e.g. urbanization; Dharmarajan et al., 2021) and the distribution patterns of wild bees (e.g. Fournier et al., 2020). Urban ecosystems can harbour large and diverse bee communities, helping to preserve and promote wild bee diversity. Although urbanization has major negative impacts on biodiversity (Theodorou et al., 2020), a significant number of wild bee species can thrive in cities. Documentation of wild bees in urban ecosystems has frequently indicated diverse wild bee communities (Baldock et al., 2015; CasanellesAbella, Chauvier, et al., 2021), although this ultimately depends on each species' traits and its response to urbanization. Urban ecosystems are warmer (Roth et al., 1989), have higher landscape heterogeneity (Turrini & Knop, 2015) and are generally less polluted by pesticides (Scheyer et al., 2007) than intensive agricultural areas. Moreover, while intensified agricultural systems have impoverished floral resources, in cities these resources might be maintained, thanks to social investment, high availability of woody species (e.g. street trees) in highly urbanized areas, and the presence of flowerrich habitats (Somme et al., 2016; Tew et al., 2021). In both public and private urban greenspaces, there are important efforts to establish and maintain flowering plant assemblages, with each phase reflecting the preferences and needs of the specific owners and managers (Harrison & Winfree, 2015). Therefore, there is a major opportunity to promote wild bee fitness and reproduction by increasing and improving wild bee habitats. Urban ecosystems can induce dietary changes in species, due to their distinct availability of various food resources. In urban ecosystems, natural food resources are complemented with anthropogenic food resources (Faeth et al., 2005), whose accessibility is modulated by each species' diet specialization. Urban floral resources are especially diverse in cities, due to gardening and horticultural activities, with many native and exotic species planted for different purposes. Some of these species provide additional sources of food for pollinators within their range of foraging preference, phenology and trait matching (Garbuzov & Ratnieks, 2014; Harrison & Winfree, 2015). Therefore, generalist wild bee species with a broad dietary range might be better able to exploit the existing urban resources, access and forage on a greater variety of patches, and consequently be more widely distributed. Knowledge on bee diet preferences could reveal which plant species are important for their survival and reproduction and could be translated into important decisions concerning the planning and management of floral resources, for example, what species to plant. Plant identification with DNA metabarcoding techniques has increased in diet studies and provides new knowledge about the feeding preferences of animals, which can help us to understand their distribution along environmental gradients (Pitteloud et al., 2021). So far, diet preferences have mostly been assessed indirectly through observations of adult bee plant visitation (e.g. Marquardt et al., 2021) or through the morphological identification of pollen grains (Haider et al., 2013; Sedivy et al., 2011). Nonetheless, specific sampling methodologies, such as trapnests that allow standardized sampling (Staab et al., 2018), combined with metabarcoding techniques promise to be a powerful tool to characterize and study larval bee diets (Bell et al., 2016; Keller et al., 2015). Trapnests target larval pollen and thus can better describe bee diet preferences than measurements of adult visitation, while pollen metabarcoding techniques can identify a larger number of taxa with a higher taxonomic resolution than pollen morphological identification; this is particularly useful in urban ecosystems with unique and rich plant pools. Metabarcoding techniques reduce the need for taxonomic expertise 5. Policy implications. Satisfying larval dietary requirements is critical to preserving and enhancing wild bee distributions within urban gradients. For high to intermediate levels of feeding specialization, we found considerable consistency in the preferred plant families or genera across the studied cities, which could be generalized to other cities where these bees occur. Identifying larval floral preferences (e.g. using pollen metabarcoding) could be helpful for identifying key plant taxa and traits for bee survival and for improving strategies to develop beefriendly cities. KEYWORDS cavitynesting bees, feeding behaviour, remote sensing, species distribution models, trapnests, urban biodiversity, urbanization 11Land Change Science, Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, Switzerland Correspondence Joan CasanellesAbella Email: [email protected] Funding information Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung, Grant/Award Number: 31BD30_172467; Narodowe Centrum Nauki, Grant/Award Number: NCN/2016/22/Z/ NZ8/00004; ERANet BiodivERsA, Grant/ Award Number: BiodivERsA32015104 Handling Editor: Margaret Stanley | 3 Journal of Applied Ecology CASANELLESABELLA Et AL. associated with pollen morphological identification, thus broadening its application across multiple sites. Here, we investigated the larval diet and distribution of four widespread wild bee species of urban ecosystems, representing a gradient of decreasing diet specialization (Chelostoma florisomne, Osmia bicornis, Osmia cornuta and Hylaeus communis), along urban intensity gradients in five European cities (Antwerp, Paris, Poznan, Tartu and Zurich). In particular, we asked the following questions: (a) What is the taxonomic and traitbased composition of the bee diets in different urban areas? (b) How consistent are the bee diets across urban areas? (c) How does diet specialization influence the bee distribution in urban ecosystems? We expected that specialized bee species (i.e. C. florisomne) would have strong preferences for specific plant taxa and thus a highly consistent diet across urban areas and within urban gradients. Conversely, we predicted that more generalist bees (i.e. O. cornuta, O. bicornis and H. communis) would have a more flexible diet and be capable of switching to alternative floral resources, including exotic taxa, and thus have a higher turnover in the diet composition (at the plant family, genus and species levels) and a less consistent diet. Finally, we hypothesized that bee species with greater diet specialization would have low flexibility in terms of switching their diet to other plant taxa and thus would be more sensitive to urban intensity. 2 | MATERIALS AND METHODS 2.1 | Cities and study sites We investigated wild bee diets in five cities in Europe: Antwerp (Belgium), Greater Paris (France, hereinafter referred as Paris), Poznan (Poland), Tartu (Estonia) and Zurich (Switzerland), covering a large part of the climatic variability in mainland Europe. Site selection followed CasanellesAbella, Frey, et al. (2021). Overall, we selected sites from the urban green areas mapped and defined in the panEuropean Urban Atlas (EEA, 2012). We used an orthogonal gradient of patch area and connectivity. In particular, we calculated connectivity using the proximity index (PI), which considers the area and the distance to all nearby patches with a favourable habitat, within a given search radius. We considered as favourable habitat the land cover classes urban green areas, forest and low urban density with <30% impervious surface. We set the search radius to 5 km from each focal patch, the maximum value possible with the available cartography. This resulted in the final selection of 80 sites: 32 in Zurich and 12 in each of the remaining four cities (Figure S1; Table S1). We maintained a minimum distance of 500 m between selected sites (except for two sites in Zurich selected based on their position in the patch and connectivity gradient, which were separated by 260 m). 2.2 | Bee sampling At each site, we installed trapnests in trees, and in three cases (one in Paris, one in Tartu and one in Zurich) in other vertical structures (e.g. lamp post). We constructed trapnests with reeds and cardboard tubes (Figure S1). Our sampling trapnests consisted of a standardized wood box with three plastic pipes 15 cm in diameter and 20 cm long. We assembled the first two pipes using 200– 300 reeds from Phragmites australis (Cav.) Trin. and 5– 10 bamboo reeds with diameters of 1– 10 mm and a length of 20 cm to cover all requirements of the cavitynesting bee community. We assembled the third pipe only using cardboard tubes of 7.5 mm diameter specific for Osmia spp. (WAB Mauerbienenzucht; Konstanz, Germany). We installed trapnests at 2.5– 3.5 m height with direct sunlight and SE or SW exposition, and kept them in the field from January until October 2018. In October, we collected the trapnests and stored them at c. 5°C until February 2019, and then transferred them to a new room at ambient temperature to recreate springlike conditions. Bees hatched and were identified to the species level from February to June 2019. 2.3 | Study organisms We collected pollen from the nests of four solitary bee species, Chelostoma florisomne (Linnaeus, 1758), Hylaeus communis (Nylander, 1852), Osmia bicornis (Linnaeus, 1758) and Osmia cornuta (Latreille, 1805). These species encompass a gradient of diet specialization (i.e. the number of different plant families exploited as resources, from the oligolectic C. florisomne to the highly polylectic H. communis), differ in phenology, and are common species in urban areas in Europe. In our study, each species was present in at least three of the studied cities. For more details about the ecology of these four wild bee species, see Text S1. 2.4 | Pollen identification We extracted a total of 464 pollen samples (Table S1) from undeveloped cells (i.e. where the larva had died) in nests where at least one adult had emerged and thus taxonomic identification of the bees was possible. Specifically, for C. florisomne we used 121 samples distributed in 3 cities and 18 sites (2 in Antwerp, 1 in Tartu and 15 in Zurich), for O. cornuta we used 66 samples distributed in 3 cities and 20 sites (6 in Antwerp, 5 in Paris and 9 in Zurich), for O. bicornis we used 176 samples distributed in 5 cities and 37 sites (3 in Antwerp, 10 in Paris, 8 in Poznan, 1 in Tartu and 15 in Zurich), and for H. communis we used 101 samples distributed in 5 cities and 33 sites (4 in Antwerp, 6 in Paris, 6 in Poznan, 9 in Tartu and 8 in Zurich). DNA metabarcoding (isolation, amplification and sequencing) was performed by AllGenetics laboratories (AllGenetics & Biology SL; A Coruña, Spain). We followed the method described by Sickel et al. (2015) and Vierna et al. (2017) to produce a pooled amplicon library on the ITS2 genomic region for the Illumina platform (Illumina). See Text S3 for details on the laboratory procedure of pollen metabarcoding. Bioinformatics followed mainly the procedure described in Campos et al. (2021) with minor modifications: We used VSEARCH 4 | Journal of Applied Ecology CASANELLESABELLA Et AL. v2.14.2 (Rognes et al., 2016) to join paired ends of forward and reverse reads. We also used VSEARCH to remove reads shorter than 200 bp, complete quality filtering (EE < 1; Edgar & Flyvbjerg, 2015) and denovo chimera filtering, and define amplicon sequence variants (ASVs), as previously done successfully for pollens (e.g. Wilson et al., 2021). The ITS2 rDNA reads were first directly mapped with VSEARCH global alignments and an identity cutoff threshold of 97% against a floral ITS2 reference database generated with the BCdatabaser (Keller et al., 2020), which consisted of plants recorded within the study regions. For the remaining unclassified reads, we first used global alignments against a global reference database (Ankenbrand et al., 2015; Keller et al., 2015). For reads that were still unclassified, we used SINTAX (Edgar, 2016a, 2016b) to assign taxonomic levels as deep as possible but a maximum of genus level with the same global reference database. In total, 82% of species recorded at the sites were present in the local database and 83% of species in the global database (direct classification). Furthermore, 92% of genera were covered by the global database for hierarchical classification. Please note that the global database contains 112,115 unique species, 11,321 genera and 710 families in total, with a very high likelihood of coverage for species and genera of any interest for anthropogenic use, including exotic garden species. 2.5 | Environmental variables We assembled variables that were potential drivers of bee diets and distributions and that represented different aspects of the urban environmental gradients. Specifically, we focused on proxies of stress (particularly thermal stress), amount of habitat and resource availability at different spatial scales. We inferred resource availability at the local scale by performing floristic inventories on standardized plots, as explained in CasanellesAbella, Frey, et al. (2021) and in Supplementary Text S4. Furthermore, we collected information on two functional plant traits sensitive to bee– plant interactions, that is, growth form (Tables S2 and S3) and blossom type (Tables S2 and S3) using information available in CasanellesAbella, Frey, et al. (2021). See Text S4 and Tables S2 and S3 for additional information on the definition of the traits. We used local and landscape connectivity metrics, local land cover metrics and landscape remotesensingbased indices to infer thermal stress and the amount of available habitat, particularly regarding resource availability. As connectivity metrics, we used patch size and the proximity index. We obtained the local land cover metrics by mapping grasslands, artificial surfaces, bare land, coniferous trees and deciduous broadleaved trees and then calculating their proportions at different spatial scales (i.e. 8, 16 and 32 m) from the focal trapnest (see Text S5 for additional details). Finally, we used remotesensingbased indices on land surface temperature, impervious surfaces, soil, water and vegetation at different spatial scales (50, 100, 200, 400, 800, 1,600 m). Specifically, we used land surface temperature (LST), the urban index (UI), the colour index (CI), the normalized difference water index (NDWI) and the normalized difference vegetation index (NDVI), which can be used to characterize existing vegetation and urban infrastructure. In addition, we performed a principal component analysis (PCA) on the explanatory variables to define new meaningful underlying variables while reducing the dimensionality of the data set (see Section 2.6). See Text S6 for details on the calculation of the remotesensingbased indices and Figure S2 for the distribution of values of each predictor in each city. 2.6 | Statistical analysis We conducted all analyses in R version 4.0.2 (R Core Team, 2021) with RStudio version 1.4.1106 (RStudio Team, 2020). 2.6.1 | Species diet analysis We performed taxonomic and traitbased metrics on the bee diets at the city and site levels. Specifically, we computed the proportion of different plant taxa at the family, genus and species levels (Table S4). Furthermore, we calculated the species, genus and family richness and the Shannon diversity index. Concerning traitbased responses, we calculated the proportion of the different categories of the three studied traits (Table S5). For each of the four studied bees, we performed Pearson correlations to investigate the relationships between the taxonomic and the traitbased diet metrics with the proxies of urban intensity, habitat amount and resource availability. We assessed these relationships (a) for each single city and (b) for all the cities combined. We calculated the pairwise correlations between cities for each bee species to study diet consistency. Specifically, we first assembled binary trophic interaction matrices between the four bee species and the plant taxa at the family, genus and species levels and then calculated the Pearson correlations of the binary trophic interaction matrices between pairs of cities for each bee and plant level. However, the trophic interaction matrix for a given city, and thus the pairwise correlations between cities, is influenced by the available plant pool. To avoid effects of plant composition, we first created a list with the plant pool occurring in each city at the family, genus and species levels. We used the plant species sampled within a 100m buffer by CasanellesAbella, Frey, et al. (2021) and complemented with the plant species recorded in GBIF (2021) at each city for the period 2000– 2018 (Figure S3). If a plant family, genus or species was missing in one of the plant pools of a pair of cities, we removed the interaction when performing the correlations. Moreover, we computed a Chisquared (χ2) test on the family and trait composition between cities' plant species pools (Text S7, Table S6). 2.6.2 | Species distribution of urban gradients We studied bee distribution patterns with species distribution models (SDMs). We assembled occurrence matrices indicating the | 5 Journal of Applied Ecology CASANELLESABELLA Et AL. occurrence of each bee species in the different sites. From the candidate environmental predictors, we evaluated the statistical relevance of each predictor using the predictive power (D2; Table S7) and then manually picked three predictors that had correlations <0.7 to avoid collinearity (Figure S4) for each bee separately. We used an ensemble of two common modelling techniques to account for model uncertainty and specificity (Buisson et al., 2010). Specifically, we used two regressionbased models, that is, generalized linear models (GLMs) and generalized additive models (GAMs), and two treebased models, that is, gradient boosting machines (GBMs) and random forests (RFs), that show a higher complexity in their fitting procedures than GLMs and GAMs. We used city as a fixed factor to account for the nested structure of the data with a binomial probability distribution. We parameterized each modelling technique in the following way: we calibrated GLMs with firstorder polynomials, GAMs with a spline smoothing term of intermediate complexity (k = 4), RFs with a node size of 5 (nodesize = 5) and 1,000 trees, and GBMs with an interaction depth of 1, a shrinkage of 0.001 and 1,000 trees. We ran the models using the r packages mgcv version 1.830, RandomForest version 4.614 and gbm version 2.1.5. We randomly split the species records of the four bees into two sets containing 80% of the data for model calibration and 20% of the data for model evaluation. We repeated the procedure five times. We assessed model performance with the True Skill Statistic (TSS; Allouche et al., 2006). TSS evaluates model skill in distinguishing absences from presences. The predictive performance of the different models was deemed acceptable when TSS > 0.4, following a commonly used minimum threshold (Thuiller et al., 2019). Thus, we discarded models with TSS values lower than 0.4. We used the selected models of each studied bee species to predict the probability of occurrence over the environmental space of the studied cities. 3 | RESULTS 3.1 | Species diet analysis A total of 41 plant families, 93 genera and 135 species were identified from the nests of the four bee species (Figure 1; Tables S8 and S9). Over half of the species were native (55%), there were more herbs (42%) than trees (34%), and dishbowl blossoms were more common (56%) (Figures S6– S8; Table S9). The number of plant species per bee nest was similar among bee species (Table S9). The total number of collected plant taxa varied greatly among bee species, reflecting their differences in diet specialization: 1 family and 4 species in C. florisomne, 12 families and 33 species in O. cornuta, 18 families and 51 species in O. bicornis, and 32 families and 81 species in H. communis (Figure 1; Table S9). At the city level, we found dominance patterns in pollen abundance for some bees (Figure 1). In O. bicornis, the most abundant species in pollen were Quercus robur (Antwerp, 70%) and Acer pseudoplatanus (Paris, 64%; Poznan, 44%; and Zurich, 33%). In H. communis, Styphnolobium japonicum was the most abundant species in pollen but only in Paris (52%) and Poznan (32%), with the vast majority of species representing 1%– 14% of the pollen abundance. Interestingly, in C. florisomne, the most abundant Ranunculus spp. in pollen changed between cities (R. acris in Antwerp, R. repens in Tartu and R. bulbosus in Zurich). Finally, in O. cornuta, no species made up more than 37%, being the most abundant ones A. pseudoplatanus (Antwerp, 37%; Paris, 21%; and Zurich, 24%) and Prunus lusitanica (Paris, 33%). In addition, very few nests were dominated by a single plant species and mostly in C. florisomne (Figure S6). We found different levels of diet conservatism across cities at the plant family and plant genus levels, according to the bee specialization degree and taxonomic resolution of the plant taxa (Figure 2). At the family and genus levels, diet consistency was high for C. florisomne and declined with broader feeding niches, particularly at the genus level (Figure 2a; Table S10). Conversely, we found major variation at the plant species level, which was particularly prominent for the broad generalist H. communis, which switched from herbaceous pollen to tree pollen with increasing urbanization (Figure 3b). The extent to which bee diet taxonomic and traitbased composition were conserved also varies according to the degree of specialization. Chelostoma florisomne had the most conserved diet, composed exclusively of native Ranunculus spp. (Figures 1 and 2a; Tables S8 and S9). Osmia cornuta primarily collected the pollen of native tree and shrub species, mainly with dishbowl or brush type blossoms, from the families Sapindaceae, Salicaceae and Rosaceae (Figure 1; Figures S7– S10; Tables S4, S5, S8– S11). Nevertheless, in Paris and Zurich, we also found a considerable proportion of Ranunculaceae (Figure 1b). Both C. florisomne and O. cornuta taxonomic and traitbased metrics showed no or little variation along urban intensity gradients (Figure S11; Table S12). In O. bicornis, native tree species with dishbowl or brush blossoms from the families Sapindaceae and Fagaceae represented a large part of the diet (Figure 1; Figures S7– S10; Tables S4, S5, S8– S11), but we found some variation across cities concerning the remaining species in the diet (Figure 1b; Figure S10). Hylaeus communis had the most diverse and variable diet. The Fabaceae family represented 34% in Paris and 42% in Poznan of the species found in the larval diet (Figure 1b) and a minor part in the remaining cities (Figure 1b). Furthermore, exotic species were more frequent for H. communis than for the other three bee species (Figure 3a; Figure S7; Tables S5 and S9). In addition, we found family richness, species richness and pollen diversity to be positively correlated with NDVI for H. communis for all cities except Poznan (Figure S11; Table S12). Finally, in Paris, an important part of the diet was trees with flag type blossoms (Figure 3b), in part due to the contribution of Styphnolobium japonicum (Figure 1a; Table S8), which became more dominant in the diet with increasing urban intensity (e.g. decreasing NDVI and increasing UI and CI at different scales; Figure 3; Table S12). 6 | Journal of Applied Ecology CASANELLESABELLA Et AL. | 7 Journal of Applied Ecology CASANELLESABELLA Et AL. 3.2 | Species distributions along urban gradients The PCA conducted on the explanatory variables returned two main axes that explained 38% and 12% of the variation, respectively. The first axis was composed of remotely sensed variables, with larger values of the PC axis indicating less vegetation (i.e. higher UI, CI and LST, and lower NDVI; Figure S12) independently of the landscape scale considered. The second axis was mostly composed of local land cover variables and metrics representing the available floral resources. Specifically, larger values on the PC axis indicated larger proportions of grasslands and lower proportions of deciduous trees and artificial surfaces, independently of the local scale considered (Figure S12). We found two distinct distribution patterns of the four bee species along urban intensity gradients. The first type of response was composed of C. florisomne, O. cornuta and O. bicornis. The probability of occurrence of these three bee species decreased rapidly with increasing urban intensity at the landscape scale (Figure 4a,b; Figure S13) and increased with higher proportions of grasslands at the local scale (Figure 4c,d). Strikingly, the probability FIGURE 1 Bee larval diet composition in the studied cities. (a) For each bee species, the collected plant species in each city where the bee species was recorded (three cities for Chelostoma florisomne and Osmia cornuta, five cities for Osmia bicornis and Hylaeus communis) are shown. The size of the circle represents the mean relative abundance of plant species contributing to pollen samples per city and bee species. (b) For each bee species, the proportion in the pollen of the different collected plant families in the studied cities is shown (mean relative abundance of plant species contributing to pollen samples per city and bee species). Only families with a proportion in pollen ≥ 0.01 are plotted, whereas the remaining ones are represented in the category ‘Other families’. Note that the proportion in pollen for O. bicornis in Antwerp and Tartu was obtained using only four and one samples, respectively. For each bee and city, we provide the number of sites where pollen samples were taken, and the total number of samples. Information on the computation of the phylogenetic tree can be found in Text S2 and Figure S5. Data supporting (a) can be found in Table S8. Ast, Asteraceae; Api, Apiaceae; Fab, Fabaceae; Fag, Fagaceae; Sal, Salicaceae; Br, Brassicaceae; Ma, Malvaceae; Sap, Sapindaceae; Ra, Ranunculaceae; Pa, Papaveraceae FIGURE 2 Pairwise correlations of the larval diet composition among cities. For each of the four studied bee species, the city pairwise correlations of the collected plant taxa are shown at the family (a), genus (b) and species (c) levels. The colour of the dots indicates the value of the correlation, with lower and higher values in orange and blue, respectively. Note that the correlation values are expressed as absolute values. Note also that the pollen for Osmia bicornis in Antwerp and Tartu was obtained using only four and one samples, respectively. Data supporting Figure 2 can be found in Table S10 8 | Journal of Applied Ecology CASANELLESABELLA Et AL. | 9 Journal of Applied Ecology CASANELLESABELLA Et AL. of occurrence peaked at low proportions of broadleaved trees (Figure 4c,d; Figure S13), even though the diets of O. cornuta and O. bicornis were composed largely of tree pollen. By contrast, the probability of occurrence of H. communis remained constant for larger values of urban intensity (Figure 4a,b; Figure S13). Moreover, the probability of occurrence was positively affected by the amount of deciduous broadleaved trees at local scales (Figure 4c,d; Figure S13). FIGURE 3 Traitbased larval diet composition in Hylaeus communis. (a– c) Composition of the diet according to the origin status (a), growth form (b) and blossom class (c) of the plant species in the larval pollen. (d– g) firstorder GLMs of the proportion of the different plant trait levels in relation to the first (d– f) and second (g– i) PC axes for the origin status (d and f), growth form (e and g), and (f and i) blossom class. Grey shaded bands indicate 95% confidence intervals. Higher PC1 values indicate less vegetation (lower normalized difference vegetation index) and more artificial surfaces (urban index, land surface temperature, colour index). Higher PC2 values indicate a larger proportion of grasslands and lower proportion of deciduous trees at local scales. PC1 explained 38% of the variation and PC2 explained 12% of the variation. See Figures S7– S9 for the traitbased composition and change along urban gradients of the other three bee species FIGURE 4 Bee distribution along urban gradients. (a and c) Loess smoothing of the mean predicted probability of occurrence of the four bee species in relation to the (a) first PCA axis (PC1) and (c) second PCA axis, performed on the explanatory variables, representing 38% and 12% of the variation, respectively. The mean predicted probability of occurrence results from the predicted probabilities of occurrence of the models with TSS > 0.4. Bands represent 95% confidence intervals. (b and d) Variation in the explanatory variables contributing the most to PC1 (b) and PC2 (d). (b) Larger values of PC1 correspond to higher values of impervious surfaces (urban index, UI), bare land (colour index, CI) and land surface temperature (LST) and lower vegetation cover (normalized difference vegetation index, NDVI) at different landscape scales (i.e. 100 and 400 m). (d) Lower values of PC2 correspond to higher proportions of deciduous trees and lower proportions of grasslands and floral resources (plant species richness) at local scales. Other scales and variables have been omitted here for simplicity (see Figure S13). See also Figure S12 for more details on the PCA