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

can


| INTRODUC TI ON
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 land-use 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;Casanelles-Abella, 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 flower-rich 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 trap-nests 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). Trap-nests 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 bee-friendly cities. 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 trait-based 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.

| 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 . Overall, we selected sites from the urban green areas mapped and defined in the pan-European 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).

| Bee sampling
At each site, we installed trap-nests 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 trap-nests with reeds and cardboard tubes ( Figure S1). Our sampling trap-nests 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 cavity-nesting 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 trap-nests 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 trap-nests and stored them at c. 5°C until February 2019, and then transferred them to a new room at ambient temperature to recreate spring-like conditions.
Bees hatched and were identified to the species level from February to June 2019.

| Study organisms
We collected pollen from the nests of four solitary bee species, 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.

| Pollen identification
We extracted a total of 464 pollen samples (Table S1)  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 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 de-novo 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 cut-off 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 Keller et al., 2015). For reads that were still unclassified, we used SINTAX (Edgar, 2016a(Edgar, , 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.

| 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  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 Casanelles- .
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 remote-sensing-based 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 broad-leaved trees and then calculating their proportions at different spatial scales (i.e. 8, 16 and 32 m) from the focal trap-nest (see Text S5 for additional details). Finally, we used remote-sensing-based 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 remote-sensing-based indices and Figure S2 for the distribution of values of each predictor in each city.

| Species diet analysis
We performed taxonomic and trait-based 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 trait-based 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 trait-based 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  Table S6).

| Species distribution of urban gradients
We studied bee distribution patterns with species distribution models (SDMs). We assembled occurrence matrices indicating the 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 (D 2 ; 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 regression-based models, that is, generalized linear models (GLMs) and generalized additive models (GAMs), and two tree-based 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 first-order 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.8-30, RandomForest version 4.6-14 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.

| 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 dish-bowl 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) Table S9). At the city level, we found dominance patterns in pollen abundance for some bees (Figure 1).
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 par-  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).

| 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 (Figure 4c,d). Strikingly, the probability F I G U R E 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 F I G U R E 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 of occurrence peaked at low proportions of broad-leaved trees  Figure S13).

F I G U R E 3
Trait-based 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) first-order 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 trait-based composition and change along urban gradients of the other three bee species F I G U R E 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

| D ISCUSS I ON
With biodiversity loss increasing at an unprecedented speed (Leclère et al., 2020) and urbanization expanding worldwide, there is a pressing need to better understand the niche requirements of different taxa to satisfy them in urban ecosystems through targeted planning and management (Aronson et al., 2017). Here, using trap-nests and pollen metabarcoding techniques, we studied the larval diet and distribution patterns of four wild bee species in five European cities.
The larval diet of the studied oligolectic and intermediate polylectic bees was conserved in terms of plant family composition, representing successful strategies in urban ecosystems as an alternative to generalism. As hypothesized, diet consistency decreased with decreasing wild bee degree of specialization and was relatively high at the plant family level and to lower extend at the plant genus level. Moreover, for C. florisomne, O. cornuta and O. bicornis, bee diet at the plant family level in the studied cities was consistent with prior studies in non-urban ecosystems (e.g. Haider et al., 2014). Diet conservatism regarding plant family composition has also been observed in other polylectic bee species such as Bombus spp. (Wood et al., 2021) and Osmia spp. Vaudo et al., 2020), indicating that diet conservatism might be widespread despite a lack of data. Finally, diet consistency at the plant species level, with the exception of highly specialized C. florisomne, was very low, possibly because the specific trait values (e.g. pollen properties imposing cognitive or physiological restrictions; Vaudo et al., 2020) driving plant-pollinator interactions in less specialized wild bees might be relatively consistent within plant families or genera.
Across different animal taxa, broad feeding niches (i.e. low diet conservatism) have been identified as a key characteristics to pass the urban ecological filtering, and hence to thrive in cities (Fournier et al., 2020), with some species even undergoing rapid phenotypic changes to broaden their diet (e.g. Eggenberger et al., 2019). While widely distributed species inside cities are typically generalists, especially when highly mobile and with specific nesting and sociality modes (i.e. cavity nesting or social), intermediate specialization might be an alternative advantageous strategy in urban ecosystems.
Focusing on wild bees, species whose preferred plant families are selected and facilitated by stakeholders (in Central European cities, e.g. Rosaceae, Fagaceae, Salicaceae; Ossola et al., 2020) could have very diversified resources, due to the co-occurrence of spontaneous and cultivated native and exotic plant species. Here, we show that this is the case for O. cornuta and O. bicornis. Consequently, wild bees could still switch pollen hosts at the species level and better exploit existing resources while avoiding the costs associated with broad niches, as some pollen types are indigestible or toxic (Eckhardt et al., 2014;Praz et al., 2008). On the other hand, strict specialization, as in C. florisomne, can still be a successful strategy in urban ecosystems when pollen hosts are highly facilitated and widespread in urban ecosystems. Nonetheless, strict specialists are vulnerable to the partial or total loss of pollen hosts, due to, for example, urban sprawl across the habitat of the bee or pollen host species, changes in social investment when gardening or managing greenspaces, or the arrival of pests.
The degree of specialization in bees is associated with distribution patterns along urban gradients. Our results suggest that increasing specialization leads to a higher sensitivity towards increased urban intensity. Diet specialization determines the possibility of occupying new patches, and greater specialization represents a strong limitation when the nutritional requirements are not met (e.g. with agricultural intensification; Peters et al., 2021). Rarity in bumblebees has been associated with narrower feeding niches (Goulson, 1999), and this is also the case for other wild bees (e.g. Deguines et al., 2016).
Some types of urban greenspaces consistently contain high bee diversity, including several specialists (Baldock, 2020;Salisbury et al., 2015). Still, only a handful of usually broad generalists can colonize various types of urban greenspaces. Hence, these species are widespread within urban ecosystems (Casanelles-Abella, Chauvier, et al., 2021;Fournier et al., 2020), even though not only the degree of specialization but also other correlated functional traits are involved in the response to urbanization (e.g. stress tolerance, dispersal; Harrison & Winfree, 2015). To promote a larger number of bees in previously unoccupied areas of the city, wild bee habitat must be strategically increased, including enhancing the availability of high floral resource diversity. A first step to achieve this is to map the existing floral resources within the different urban land covers, making use of ongoing inventories (Ossola et al., 2020) or sampling schemes from research projects (Baldock et al., 2019;. These products could be combined and compared with wild bee diet preferences to detect where and what kind of planning actions can be taken. For instance, strategically increasing plant diversity rather than overall quantity with targeted taxonomic and trait groups is a measure that has been successful in other ecosystem types (e.g. agroecosystems; Sutter et al., 2017), although other wild bee requirements (e.g. nesting mode) must also be satisfied to successfully promote them (Requier & Leonhardt, 2020).

| Limitations and prospects
Pollen metabarcoding only yields relative abundance data and can be subject to PCR or taxon-dependent biases as discussed in several studies. (e.g. Bell et al., 2016;Keller et al., 2015) The semiquantitative abundances obtained in metabarcoding analyses are however still useful to identify the relative proportions of taxa and differentiate dominant, common and rare contributions to a mixed pollen sample by showing correlations (even though not perfect) between read and grain numbers . It is further agreed that pollen metabarcoding can identify to deep taxonomic levels, is able to detect rare taxa and is well comparable between studies (Bell et al., 2016). In addition, the prevalence of certain plant families and species across cities (Figure 1), and the existing information on C. florisomne, O. cornuta and O. bicornis (Haider et al., 2013Sedivy et al., 2008) showing patterns similar to those observed here support our findings. Combining molecular and morphological approaches would provide more robust relative abundance estimates as suggested by Sickel et al. (2015). Furthermore, if taking dry weight of the collected pollen samples is an option (e.g. Lihoreau et al., 2012), these relative data could be converted to absolute abundances per sample. Furthermore, we could only use morphological plant traits to investigate diet composition and shifts along urban gradients. However, bee choices when selecting larval food are likely to be influenced by the nutritional properties of the pollen, such as the content of protein, sugar or other essential nutrients (Vanderplanck et al., 2014;Vaudo et al., 2020). Evidence seems to indicate important differences in the nutritional value of the available floral resources between urban land cover types (Tew et al., 2021) and possibly between ornamental plant species (Garbuzov & Ratnieks, 2014). Incorporating nutritional traits in future dietary studies will better elucidate the mechanisms behind diet composition and trophic niche shifts.

| Policy implications
The uncovered taxonomic and functional diet preferences can support the planning and management of urban greenspaces to promote wild bees, particularly in wild bee species where the diet patterns are consistent across cities. Typical common urban weeds (e.g. Taraxacum officinale, Bellis perennis, Trifolium pratensis, T. repens) that are important floral resources in urban areas (Baldock et al., 2019;Kanduth et al., 2021;Larson et al., 2014) contributed little or nothing to the larval diet of the studied bees. Conversely, both native and exotic woody species proved to be a widely used floral resource.
For example, O. cornuta and O. bicornis collected large amounts of tree pollen from different plant families, which might enable them to exploit secondary pollen hosts by mixing high-quality pollens with less digestible ones (e.g. Ranunculus spp., Asteraceae;Eckhardt et al., 2014;Praz et al., 2008). Thus, the maintenance of different vegetation and urban habitat types (e.g. meadows, street trees, shrublands) is of major importance for preserving bees. In addition, our results indicate that the occurrence of specific plant taxa (e.g. at the family or genus level) or trait values is more important than the origin of the plants (Harrison & Winfree, 2015). For example, non-native plants have been observed to retain the blooming time of their original region (Godoy et al., 2009), and thus plant species from the Northern Hemisphere (e.g. planted urban trees or shrubs such Acer spp., Salix spp., Crataegus spp and Quercus spp.) may bloom in synchrony with natives species, providing additional resources for wild bees. Furthermore, urban trees have been shown to be an important resource for several wild bee species (Somme et al., 2016).
For instance, in our study, H. communis increasingly collected more pollen from ornamental trees (mainly Styphnolobium japonicum) with decreasing amounts of greenspace. Because the tree distribution in cities is mostly driven by anthropogenic factors, it represents an important point of action for greening strategies, specifically in densely urbanized city areas with limited herbaceous vegetation.
Overall, floral preferences obtained from pollen metabarcoding, particularly when combined with existing information on available floral resources, could help to improve current strategies for developing bee-friendly cities (e.g. in the EU, see Wilk et al., 2020). In particular, characterizing bee diets could inform planning, management and decision-making (e.g. what species, genera or families to plant), involving stakeholders, policymakers, nurseries and plant centres, for urban greenspaces to preserve and further promote wild bees in urban ecosystems.