Joan Casanelles-Abella

Stéphanie Muller

Alexander Keller

Cristiana Aleixo

Marta Alós Orti

François Chiron

Nicolas Deguines

Tiit Hallikma

Lauri Laanisto

Pedro Pinho

Stefanie Müller

Roeland Samson

Piotr Tryjanowski

Anskje Van Mensel

Loïc Pellissier

Marco Moretti

+ Affiliations

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

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

KEYWORDS: cavity-nesting bees, feeding behaviour, remote sensing, trap-nests, urbanization, urban biodiversity, species distribution models Title page= 238; Abstract = 333; Manuscript+Captions+References = 8421

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. [START_REF] Potts | Safeguarding pollinators and their values to human well-being[END_REF]. Over the last decades, bee populations have dramatically decreased [START_REF] Zattara | Worldwide occurrence records suggest a global decline in bee species richness[END_REF]. Multiple causes have been identified [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF], and the loss of floral resources is among the most important [START_REF] Goulson | Bee declines driven by combined stress from parasites, pesticides, and lack of flowers[END_REF]. 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 [START_REF] Goulson | Foraging strategies of insects for gathering nectar and pollen, and implications for plant ecology and evolution[END_REF][START_REF] Vanderplanck | How does pollen chemistry impact development and feeding behaviour of polylectic bees?[END_REF]. Consequently, diet specialization and diet preference are two key traits determining the sensitivity to land-use changes (e.g. urbanization; [START_REF] Dharmarajan | Anthropogenic disturbance favours generalist over specialist parasites in bird communities: Implications for risk of disease emergence[END_REF] and the distribution patterns of wild bees (e.g. [START_REF] Fournier | The origin of urban communities: From the regional species pool to community assemblages in city[END_REF].

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 [START_REF] Theodorou | Urban areas as hotspots for bees and pollination but not a panacea for all insects[END_REF], 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 [START_REF] Baldock | Where is the UK's pollinator biodiversity ? The importance of urban areas for flower-visiting insects[END_REF][START_REF] Casanelles-Abella | Applying predictive models to study the ecological properties of urban ecosystems: A case study in Zürich, Switzerland[END_REF], although this ultimately depends on each species' traits and its response to urbanization. Urban ecosystems are warmer [START_REF] Roth | Satellite-derived urban heat islands from three coastal cities and the utilization of such data in urban climatology[END_REF], have higher landscape heterogeneity (Turrine & Knopp, 2015) and are generally less polluted by pesticides [START_REF] Scheyer | Pesticides analysed in rainwater in Alsace region (Eastern France): Comparison between urban and rural sites[END_REF]) 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 [START_REF] Somme | Food in a row: urban trees offer valuable floral resources to pollinating insects[END_REF][START_REF] Tew | Quantifying nectar production by flowering plants in urban and rural landscapes[END_REF]. In both public and private 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 [START_REF] Harrison | Urban drivers of plant-pollinator interactions[END_REF]. 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 [START_REF] Faeth | Trophic dynamics in urban communities[END_REF], 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 [START_REF] Garbuzov | Quantifying variation among garden plants in attractiveness to bees and other flower-visiting insects[END_REF][START_REF] Harrison | Urban drivers of plant-pollinator interactions[END_REF].

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, e.g. 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 [START_REF] Pitteloud | The structure of plant-herbivore interaction networks varies along elevational gradients in the European Alps[END_REF]. So far, diet preferences have mostly been assessed indirectly through observations of adult bee plant visitation (e.g. [START_REF] Marquardt | Evaluation of the importance of ornamental plants for pollinators in urban and suburban areas in Stuttgart, Germany[END_REF] or through the morphological identification of pollen grains [START_REF] Haider | Intra-and interpopulational variation in the ability of a solitary bee species to develop on non-host pollen: Implications for host range expansion[END_REF][START_REF] Sedivy | Closely related pollen generalist bees differ in their ability to develop on the same pollen diet: Evidence for physiological adaptations to digest pollen[END_REF]. Nonetheless, specific sampling methodologies, such as trap-nests that allow standardized sampling [START_REF] Staab | Trap nests for bees and wasps to analyse trophic interactions in changing environments-A systematic overview and user guide[END_REF], combined with metabarcoding techniques promise to be a powerful tool to characterize and study larval bee diets [START_REF] Bell | Pollen DNA barcoding: Current applications and future prospects[END_REF][START_REF] Keller | Evaluating multiplexed next-generation sequencing as a method in palynology for mixed pollen samples[END_REF]. 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 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: [START_REF] Requier | Beyond flowers: including non-floral resources in bee conservation schemes[END_REF] What is the taxonomic and trait-based composition of the bee diets in different urban areas? (2) How consistent are the bee diets across urban areas? (3) 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 level) 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.

MATERIALS AND METHODS

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 Casanelles-Abella, Frey et al. (2021). 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 (Fig. 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 (Fig. 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 ca. 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, 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.

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 We followed the method described by [START_REF] Sickel | Increased efficiency in identifying mixed pollen samples by meta-barcoding with a dual-indexing approach[END_REF] and [START_REF] Vierna | PCR cycles above routine numbers do not compromise high-throughput DNA barcoding results[END_REF] to produce a pooled amplicon library on the ITS2 genomic region for the Illumina platform (Illumina, 2017). See Text S3 for details on the laboratory procedure of pollen metabarcoding. Bioinformatics followed mainly the procedure described in [START_REF] Campos | Standard methods for pollen research[END_REF] with minor modifications: We used VSEARCH v2.14.2 [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF] 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. [START_REF] Wilson | Many small rather than few large sources identified in long-term bee pollen diets in agroecosystems[END_REF]. 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 [START_REF] Keller | BCdatabaser: On-the-fly reference database creation for (meta-)barcoding[END_REF], 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 2015;[START_REF] Keller | Evaluating multiplexed next-generation sequencing as a method in palynology for mixed pollen samples[END_REF]. For reads that were still unclassified, we used SINTAX [START_REF] Edgar | SINTAX: a simple non-Bayesian taxonomy classifier for 16S and ITS sequences[END_REF] 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, 83% of species in the global database (direct classification). Further, 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 performed taxonomic and trait-based metrics on the bee diets at the city and site level. Specifically, we computed the proportion of different plant taxa at the family, genus and species level (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: [START_REF] Requier | Beyond flowers: including non-floral resources in bee conservation schemes[END_REF] for each single city and ( 2) 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 level 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 level. We used the plant species sampled within a 100-m buffer by Casanelles-Abella, Frey et al. ( 2021) and complemented with the plant species recorded in GBIF at each city for the period 2000-2018 (Fig. 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 Chi-squared ( 2) test on the family and trait composition between cities' plant species pools (Text S7, 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 (Fig. S4) for each bee separately. We used an ensemble of two common modelling techniques to account for model uncertainty and specificity [START_REF] Buisson | Uncertainty in ensemble forecasting of species distribution[END_REF]. Specifically, we used two regression-based models, i.e. generalized linear models (GLMs) and generalized additive models (GAMs), and two tree-based models, i.e. 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 1000 trees, and GBMs with an interaction depth of 1, a shrinkage of 0.001 and 1000 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; [START_REF] Allouche | Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS)[END_REF]. 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 [START_REF] Thuiller | Uncertainty in ensembles of global biodiversity scenarios[END_REF]. 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.

RESULTS

Species diet analysis

A total of 41 plant families, 93 genera and 135 species were identified from the nests of the four bee species (Fig. 1, Tables S8-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%) (Figs 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. (Fig. 1, Table S9). At the city level, we found dominance patterns in pollen abundance for some bees (Fig. 1). In O.

bicornis, the most abundant species in pollen were Quercus robur (Antwerp, 70%), 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 %; Zurich, 24 %) and Prunus lusitanica (Paris, 33 %). In addition, very few nests were dominated by a single plant species and mostly in C. florisomne (Fig. S6).

We found different levels of diet conservatism across cities at the plant family and plant genus level, according to the bee specialization degree and taxonomic resolution of the plant taxa (Fig. 2). At the family and genus level, diet consistency was high for C. florisomne and declined with broader feeding niches, particularly at the genus level (Fig. 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 (Fig. 3b).

The extent to which bee diet taxonomic and trait-based composition were conserved also varies according to the degree of specialization. Chelostoma florisomne had the most conserved diet, composed exclusively of native Ranunculus spp. (Figs 1, 2a, Tables S8-S9). Osmia cornuta primarily collected the pollen of native tree and shrub species, mainly with dish-bowl or brush type blossoms, from the families Sapindaceae, Salicaceae and Rosaceae (Figs 1, S7-S10, Tables S4,S5, S8-S11).

Nevertheless, in Paris and Zurich we also found a considerable proportion of Ranunculaceae (Fig. 1b). Both C. florisomne and O. cornuta taxonomic and trait-based metrics showed no or little variation along urban intensity gradients (Fig. 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 (Figs 1, S7-S10, Tables S4,S5, S8-S11), but we found some variation across cities concerning the remaining species in the diet (Figs. 1b,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 (Fig. 1b) and a minor part in the remaining cities (Fig. 1b). Furthermore, exotic species were more frequent for H. communis than for the other three bee species (Figs 3a, S7, Tables S5,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 (Fig. S11, Table S12). Finally, in Paris, an important part of the diet was trees with flag type blossoms (Fig. 3b), in part due to the contribution of Styphnolobium japonicum (Fig 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; Fig. 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; Fig. 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 (Fig. 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 (Figs 4a-b, S13) and increased with higher proportions of grasslands at the local scale (Fig. 4c-d). Strikingly, the probability of occurrence peaked at low proportions of broad-leaved trees (Figs 4c-d, 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 (Figs 4a-b, S13). Moreover, the probability of occurrence was positively affected by the amount of deciduous broad-leaved trees at local scales ( Figs 4c-d, S13).

DISCUSSION

With biodiversity loss increasing at an unprecedented speed [START_REF] Leclère | Bending the curve of terrestrial biodiversity needs an integrated strategy[END_REF] and urbanization expanding worldwide, there is a pressing need to better understand the niche requirements of different taxa in order to satisfy them in urban ecosystems through targeted planning and management [START_REF] Aronson | Biodiversity in the city: key challenges for urban green space management[END_REF]. 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 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. [START_REF] Haider | Phylogeny and floral hosts of a predominantly pollen generalist group of mason bees (Megachilidae: Osmiini)[END_REF].

Diet conservatism regarding plant family composition has also been observed in other polylectic bee species such as Bombus spp. [START_REF] Wood | Global patterns in bumble bee pollen collection show phylogenetic conservation of diet[END_REF] and Osmia spp. [START_REF] Haider | Phylogeny and floral hosts of a predominantly pollen generalist group of mason bees (Megachilidae: Osmiini)[END_REF][START_REF] Vaudo | Pollen Protein: Lipid Macronutrient Ratios May Guide Broad Patterns of Bee Species Floral Preferences[END_REF], indicating that diet conservatism might be widespread despite a lack of data. The observed patterns contrast with the low predictability of diet composition based on the available floral resources at all plant levels. 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; [START_REF] Vaudo | Pollen Protein: Lipid Macronutrient Ratios May Guide Broad Patterns of Bee Species Floral Preferences[END_REF] 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 [START_REF] Fournier | The origin of urban communities: From the regional species pool to community assemblages in city[END_REF], with some species even undergoing rapid phenotypic changes to broaden their diet (e.g. [START_REF] Eggenberger | Urban bumblebees are smaller and more phenotypically diverse than their rural counterparts[END_REF]. 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; [START_REF] Ossola | The Global Urban Tree Inventory: A database of the diverse tree flora that inhabits the world's cities[END_REF] 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 [START_REF] Eckhardt | Pollen mixing in pollen generalist solitary bees: A possible strategy to complement or mitigate unfavourable pollen properties[END_REF][START_REF] Praz | Specialized bees fail to develop on non-host pollen: Do plants chemically protect their pollen[END_REF]. 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 e.g. 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; [START_REF] Peters | Diets maintained in a changing world: Does land-use intensification alter wild bee communities by selecting for flexible generalists? Ecology and Evolution[END_REF]. Rarity in bumblebees has been associated with narrower feeding niches [START_REF] Goulson | Foraging strategies of insects for gathering nectar and pollen, and implications for plant ecology and evolution[END_REF], and this is also the case for other wild bees (e.g. [START_REF] Deguines | Functional homogenization of flower visitor communities with urbanization[END_REF]. Some types of urban greenspaces consistently contain high bee diversity, including several specialists [START_REF] Baldock | Opportunities and threats for pollinator conservation in global towns and cities[END_REF][START_REF] Salisbury | Enhancing gardens as habitats for flower-visiting aerial insects (pollinators): Should we plant native or exotic species?[END_REF]. Still, only a handful of usually broad generalists can colonize various types of urban greenspaces. Hence, these species are widespread within urban ecosystems [START_REF] Fournier | The origin of urban communities: From the regional species pool to community assemblages in city[END_REF], 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; [START_REF] Harrison | Urban drivers of plant-pollinator interactions[END_REF]. 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 [START_REF] Ossola | The Global Urban Tree Inventory: A database of the diverse tree flora that inhabits the world's cities[END_REF] or sampling schemes from research projects [START_REF] Baldock | A systems approach reveals urban pollinator hotspots and conservation opportunities[END_REF][START_REF] Casanelles-Abella | Applying predictive models to study the ecological properties of urban ecosystems: A case study in Zürich, Switzerland[END_REF]. 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; [START_REF] Sutter | Enhancing plant diversity in agricultural landscapes promotes both rare bees and dominant crop-pollinating bees through complementary increase in key floral resources[END_REF], although other wild bee requirements (e.g. nesting mode) must also be satisfied to successfully promote them [START_REF] Requier | Beyond flowers: including non-floral resources in bee conservation schemes[END_REF].

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. [START_REF] Bell | Pollen DNA barcoding: Current applications and future prospects[END_REF][START_REF] Keller | Evaluating multiplexed next-generation sequencing as a method in palynology for mixed pollen samples[END_REF] The semi-quantitative 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 [START_REF] Keller | Evaluating multiplexed next-generation sequencing as a method in palynology for mixed pollen samples[END_REF] . 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 [START_REF] Bell | Pollen DNA barcoding: Current applications and future prospects[END_REF]. In addition, the prevalence of certain plant families and species across cities (Fig. 1), and the existing information on C. florisomne, O. cornuta and O. bicornis [START_REF] Haider | Intra-and interpopulational variation in the ability of a solitary bee species to develop on non-host pollen: Implications for host range expansion[END_REF][START_REF] Haider | Phylogeny and floral hosts of a predominantly pollen generalist group of mason bees (Megachilidae: Osmiini)[END_REF][START_REF] Sedivy | Patterns of host-plant choice in bees of the genus Chelostoma: The constraint hypothesis of host-range evolution in bees[END_REF] 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 [START_REF] Sickel | Increased efficiency in identifying mixed pollen samples by meta-barcoding with a dual-indexing approach[END_REF]. Further, if taking dry weight of the collected pollen samples is an option (e.g. [START_REF] Lihoreau | Radar Tracking and Motion-Sensitive Cameras on Flowers Reveal the Development of Pollinator Multi-Destination Routes over Large Spatial Scales[END_REF], 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 [START_REF] Vanderplanck | How does pollen chemistry impact development and feeding behaviour of polylectic bees?[END_REF][START_REF] Vaudo | Pollen Protein: Lipid Macronutrient Ratios May Guide Broad Patterns of Bee Species Floral Preferences[END_REF]. Evidence seems to indicate important differences in the nutritional value of the available floral resources between urban land cover types [START_REF] Tew | Quantifying nectar production by flowering plants in urban and rural landscapes[END_REF] and possibly between ornamental plant species [START_REF] Garbuzov | Quantifying variation among garden plants in attractiveness to bees and other flower-visiting insects[END_REF]. Incorporating nutritional traits in future dietary studies will better elucidate the mechanisms behind diet composition and trophic niche shifts.

Finally, our results open future research directions. First, further studies could consider whether there are signals of feeding behavioural changes due to anthropogenic factors (e.g. plant or management preferences by stakeholders), as seen in other taxa (e.g. [START_REF] Cucherousset | Freshwater killer whales": beaching behavior of an alien fish to hunt land birds[END_REF], and if these are reflected in the variation of other functional traits (e.g. tongue length; [START_REF] Eggenberger | Urban bumblebees are smaller and more phenotypically diverse than their rural counterparts[END_REF]. To achieve this, it would be critical to study and compare bee diets in both urban and non-urban ecosystems. Second, assessing the effects of intraspecific diet variability on fitness (e.g. number and sex-ratio of offspring, adult and larval survival) would provide important mechanistic insights for anticipating how plant management choices might affect urban bee communities.

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 [START_REF] Baldock | A systems approach reveals urban pollinator hotspots and conservation opportunities[END_REF][START_REF] Kanduth | Red and white clover provide food resources for honeybees and wild bees in urban environments[END_REF][START_REF] Larson | Pollinator assemblages on dandelions and white clover in urban and suburban lawns[END_REF] 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; [START_REF] Eckhardt | Pollen mixing in pollen generalist solitary bees: A possible strategy to complement or mitigate unfavourable pollen properties[END_REF][START_REF] Praz | Specialized bees fail to develop on non-host pollen: Do plants chemically protect their pollen[END_REF].

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 [START_REF] Harrison | Urban drivers of plant-pollinator interactions[END_REF]. For example, non-native plants have been observed to retain the blooming time of their original region [START_REF] Godoy | Flowering phenology of invasive alien plant species compared with native species in three Mediterranean-type ecosystems[END_REF], 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 [START_REF] Somme | Food in a row: urban trees offer valuable floral resources to pollinating insects[END_REF]. 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 [START_REF] Wilk | A Guide for Pollinator-friendly Cities: Environment A guide for pollinator-friendly cities: How can spatial planners and land-use managers create favourable urban environments for pollinators? Guidance prepared by ICLEI Europe for the European Commission[END_REF]. 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. S10.

( 2 in

 2 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).

FIGURES 689 690 FIGURE 1 .

 6901 FIGURES

FIGURE 2 .

 2 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) level. 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 O. bicornis in Antwerp and Tartu was obtained using only four and one samples, respectively. Data supporting Fig. 2 can be found in Table

FIGURE 3 .

 3 FIGURE 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

FIGURE 4 .

 4 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 Fig. S13). See also Fig. S12 for more details on the PCA.

ACKNOWLEDGEMENTS 18

This study was supported by the Swiss National Science Foundation (project 31BD30_172467) within the ERA-Net BiodivERsA project "BioVEINS: Connectivity of green and blue infrastructures: living veins for biodiverse and healthy cities" (H2020 BiodivERsA32015104). We thank AllGenetics and especially A. Vizcaíno for conducting the pollen metabarcoding in the middle of a pandemic. We thank the local authorities of Zurichparticularly Grün Stadt Zurich (GSZ), Antwerp, Paris, Poznan and Tartu for supporting this study. Moreover, we acknowledge L. Roquer, J. Bosch and A. Rodrigo for their input on the trap-nest and pollen collection. We thank E. Eggenberg, K. Kilchhofer, P. Bischof, R. Veenstra, Ł. Myczko and L. Villarroya-Villalba for helping collect and process the bee data. We deeply thank B. Peters, C. Bello and C. Waldock for their input on data analysis, A. Müller for his guidance on bee biology and nutrition, and J. Bascompte for providing valuable comments on the manuscript. The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The raw data on diet composition supporting this paper can be find in the repository EnviDat v.0.1. under the DOI: 10.16904/envidat.249

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 Casanelles-Abella, Frey et al. (2021) and in Supplementary Text S4. Furthermore, we collected information on two functional plant traits sensitive to bee-plant interactions, i.e. growth form (Tables S2-S3) and blossom type (Tables S2-S3) using information available in Casanelles-Abella, Frey et al. (2021). See Text S4 and Tables S2-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, 1600 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 indices and Fig. S2 for the distribution of values of each predictor in each city.

Statistical Analysis

We conducted all analyses in R version 4.0.2 (R Core Team, 2019) with RStudio version 1.4.1106(RStudio Team, 2020).

Species diet analysis

AUTHORS' CONTRIBUTIONS

JCA, MM and LP conceived the ideas. MM, FC, LL, RS, PP and PT developed the site selection methodology. JCA collected the bee data. AK performed the bioinformatic analyses on the pollen. JCA, SM, ND, TH and MAO collected the plant data. MAO, TH and LL provided the local land cover metrics. CA and PP provided the remote sensing data. JCA analysed the data. JCA led the manuscript writing, with editing and supervision by MM and LP. All the authors reviewed and revised the manuscript and gave final approval for publication.