How city traits affect taxonomic and functional diversity of urban wild bee communities: insights from a worldwide analysis

Land-use change, including urbanization, is known to affect wild bee (Hymenoptera: Apoidea) diversity. However, while previous studies have focused on differences across local urbanization gradients, to the best of our knowledge, none focused on differences among cities at a wide geographical scale. We here used published data for wild bee communities in 55 cities across the globe, in order to explore how city traits (population density, city size, climate and land-use parameters) affect both taxonomic (diversity, distinctness, dominance) and functional (body size, nesting strategy, sociality, plant host specialization) profile of urban bee communities. By controlling for sample size and sampling effort, we found that bigger cities host few parasitic and oligolectic species, along with more above-ground-nesting bees. Cities with highly fragmented green areas present a lower proportion of oligolectic species and a higher proportion of both social species and large-bodied bees. Cities with more impervious surfaces seem to host a lower proportion of below-ground-nesting bees. Hotter cities present both a lower richness and diversity, with functional diversity highest at intermediate precipitation values. Overall, it seems that high levels of urbanization—through habitat modification and the “heat island” effect—lead to a strong simplification of the functional diversity of wild bee communities in cities. Our results may help explain the previously observed variable response of some bee community traits across local urbanization gradients.


INTRODUCTION
Bees (Hymenoptera: Apoidea) are considered the most effective pollinator insects, thus guaranteeing one of the most valuable ecosystem services (Porto et al. 2020) by allowing plant reproduction (Ollerton et al. 2011), maintaining food security (Montoya et al. 2021) and ultimately having a high relevance for global Sustainable Development Goals (Patel et al. 2021;Díaz et al. 2015). It is thus especially worrying the recently documented decline of wild bees (all bee species except the domesticated honeybee) in different countries (Biesmeijer et al. 2006), though the magnitude of such decline is still debated given that the conservation status of most wild bee species remains unknown to date (see Nieto et al. 2014 for Europe). Land-use change, including the expansion of cities (urbanization), seems however to be particularly detrimental (Goulson et al. 2015). Insects are among the organisms with the largest diversity in urbanized environments (Corcos et al. 2019) and provide several ecosystem services (Hall et al. 2017). However, bees (Biesmeijer et al. 2006) and other non-hymenopteran flower visitors (Deguines et al. 2012) have been reported to decline as a consequence of increasing urbanization (but see Hall et al. 2017 andBaldock et al. 2019), with potential impacts on pollination service.
The expansion of cities causes natural habitats to be reduced, fragmented and substituted mostly by impervious surfaces (i.e. covered by concrete) (Ayers and Rehan 2021;Geslin et al. 2016), hence leading to variable loss of natural cover. On one hand, the reduction in natural land cover has long been known to impact biodiversity (Brooks et al. 2002), either by directly leading to species loss when habitats are altered or causing indirect effects such as resources depletion (Kehoe et al. 2021). Another effect of urbanization is the habitat fragmentation of green areas, which leads urban landscapes to be configured as a mosaic of small and isolated natural patches (Semper-Pascual et al. 2021). As a result, following the island biogeography theory (Wilson and MacArthur 1967), the remaining natural patches should no longer sustain biodiversity as bigger natural areas, primarily because they can host fewer resources and of lower quality (Ayers and Rehan 2021). Moreover, in a highly fragmented habitat, there are few viable links between green patches, making it difficult for wildlife to move between patches (Zanette et al. 2005) and hence possibly reducing their interactions, fitness (Banks et al. 2007) or chances to find a food source (Ingala et al. 2019). With the projected increased level of urbanization, a reduction in the potential supply of ecosystem services might be expected (Eigenbrod et al. 2011), though a moderate level of disturbance seems to enhance biodiversity, as suggested by the intermediate disturbance hypothesis (Connell 1978).
In the light of their crucial ecological importance and the current expansion of urbanized landscapes, the scientific community has focused its attention on urban wild bee assemblages (Dicks et al. 2021). Not only species richness and diversity (Banaszak-Cibicka et al. 2018), but more recently also functional traits (Wilson and Jamieson 2019) have begun to be investigated in urbanized landscapes. The need to study both taxonomic and functional diversity is testified by the fact that wild bee abundance and functional trait diversity have been shown to synergically promote pollination (Woodcock et al. 2019). To highlight the urgency of the topic, more than a thousand studies have been published on the ecology of wild bees in urban landscapes (Buchholz and Egerer 2020). More recently, Prendergast et al. (2022) reviewed how abundance and species richness vary across different landscape types. The bulk of this body of literature has focused on comparing natural, agricultural (or rural) and urban landscapes.
Often, idiosyncratic conclusions have been drawn on whether cities can filter for some functional traits but not for others. For instance, whether oligolecty (host-plant specialist pollen collectors) or polylecty (host-plant generalist pollen collectors) is more favoured in urban environments is strongly debated. Some studies found a reduction in oligolectic species richness in cities (Twerd and Banaszak-Cibicka 2019), while others showed no significant differences (Wray et al. 2014). Conversely, increasing fragmentation appears to favour above-groundnesting bees (i.e. those nesting in holes in wood or plant stems) over below-ground nesters (i.e. those nesting in soil) (Fortel et al. 2014). Contrasting results have emerged also when investigating the role of different urban elements in shaping wild bee assemblages, as no single factor consistently influenced the abundance or richness of urban wild bees (Prendergast et al. 2020). Indeed, despite the detrimental effects of urbanization on natural habitats, cities can still constitute refugia for bees (Hall et al. 2017). For example, cities can provide novel food resources or nesting substrates (Turrini and Knop 2015). Thus, depending on the resource quality or the properties of the surrounding matrix, urban areas have the potential to support rich pollinator assemblages (Baldock et al. 2019). Things are, consequently, far from being clear on the effects of urbanization on wild bee communities in cities.
Variation in city traits can be potentially one of the reasons for the contrasting results that emerged so far for urban bee communities since these traits can affect how cities filter for taxonomic composition and functional traits relative to city-surrounding natural or rural areas.
Hence, an analysis across cities over a large geographical scale may help in understanding better bee community responses to urbanization. Here, we make a first attempt to fill this gap by performing an analysis based on previously published data for bee communities of 55 cities across the globe, in order to explore how city traits (population density, city size, climate and land-use parameters) affect both taxonomic (diversity, distinctness, dominance) and functional (body size, nesting strategy, sociality, plant host specialization) profile of urban bees ( Figure 1).

Bibliographical review and wild bee community composition
To retrieve published papers on wild bee communities in cities, we searched Scopus (https:// www. scopus. com/) using the keywords "Bees" + "City" as well as "Bees" + "Urban". The search was stopped on the 30th of November 2021. Over one thousand published papers were initially filtered through this search, and each of them was inspected in order to exclude those not useful for our study. First, we excluded works whose sampling area covered two or more distant cities or a gradient of urbanization from agricultural to urbanized landscapes and abundance or occurrence data were not separated across sites. Then, we discarded those studies with a sampling period of over a decade; this is because, over long periods, cities might have heavily transformed. We then removed works focused on a single taxon of bees. Finally, we removed works with only trap-nesting as a sampling method since these are strongly biased towards solitary species and catch exclusively cavity-nesting species (Prendergast et al. 2020). The final selection included 71 studies. Most included abundances of bee species, while some included presence/ absence data only. Few studies included more than one dataset (more than one city analysed) so the total of analysed cases was 74.
From each selected work, we characterized the composition of the wild bee community. For each sampled bee species, we associated abundance (or presence in case of binary information) and four different functional traits: sociality, nesting strategy, host-plant specialization and female body size (see Appendix 1 for a full description). For study cases with available abundance data, we calculated two diversity indices: Shannon-Wiener diversity (H') and Gini-Simpson index (GS, i.e. 1-Simpson index), which are commonly used to characterize the species diversity in a community (Shannon and Weaver 1949;Simpson 1949). H' measures uncertainty about the identity of species in the sample, and its units quantify information, while GS measures the probability that two individuals, drawn randomly from the sample, will be of different species (i.e. the lower the index, the higher the dominance). Additionally, to account for taxonomic distance among the species in each community, we calculated the taxonomic diversity (∆) (i.e. the expected path length between any two randomly picked individuals from the sample) and taxonomic distinctness (δ) (i.e. the average path length between two randomly chosen but taxonomically different individuals) (Clarke and Warwick 1998). Both ∆ and δ could be calculated also for presence/absence data (in this case, diversity and distinctness have the same value). We entered the taxonomic information on three levels: species, genus and family. When possible, the information on functional traits was taken directly from the analysed literature (Appendix 2). Where information was missing in the literature, information on functional traits was searched on a variety of additional literature sources or websites (Appendix 2). For some species, the information for traits is not known. Since we aimed to evaluate the effects of city traits on wild bees, we did not include domestic honeybee (Apis mellifera L.) in the analysis.

Demographic, geographic and climatic characterization of cities
We retrieved city size (surface), population size, population density and altitude for each city from www. wilki pedia. en. In addition, each city was associated with the 19 bioclimatic variables from WorldClim version 2.1 (released in January 2020). These variables describe different climatic measures and are the average for the years 1970-2000 (see Appendix 3 for descriptions) (Hijmans et al. 2005) and were retrieved through R version 4.0.3 using raster (Hijmans and van Etten 2012) and sp libraries (Pebesma and Bivand 2005) with a resolution of 2.5 arcsec (~ 5km 2 ) to cover as much of the city surface as possible. Each bioclimatic variable was retrieved for the central longitude and latitude of each analysed city.

Landscape characterization
For the landscape characterization, we downloaded land cover maps representing spatial information from Copernicus Global Land Service (Buchhorn et al. 2020). Within these maps, different classes of land-use are colour coded following the Land Cover Classification System (LCCS) developed by the United Nations (UN), for a total of 22 different classes (Buchhorn et al. 2017). After, we clipped each land cover map with a shapefile with the borders of the cities considered. The city borders were downloaded from different sources: Europe from "Eurostat" https:// ec. europa. eu/ euros tat/ web/ main/ home; Argentina, Brazil, India, Mexico and Turkey from "GADM version 3.6" https:// gadm. org/ index. html; Australia from "Australian Bureau of Statistics" https:// www. abs. gov. au/; Canada from "Statistics Canada" https:// www. statc an. gc. ca/ en/ start; Colombia and Sri Lanka from "The Humanitarian data exchange" https:// data. humda ta. org/; and USA from "United States Census Bureau" https:// www. census. gov/ en. html. With this operation, we obtained the landuse data for every city.
Since LCCS present different types of vegetation, we decided to aggregate all natural vegetation classes (i.e. ID 20,30,111,112,113,114,115,116,121,122,123,124,125,126) into a bigger class called "Green area". Using the GIS package LecoS (Jung 2016), we extracted different metrics for each LCCS class: land cover, landscape proportion, edge length and edge density (see Appendix 3 for descriptions). In addition to these land-use parameters, we calculated the ratio between green and impervious surface and the normalized edge density of green patches (i.e. the normalization of edge density with respect to the total green area (McGarigal 1995;Ma et al. 2013)). All land-use data mining has been performed in QGIS 3.16.

Statistical analysis
The 74 selected studies for the analysis had data on bee communities obtained through a variety of sampling methods, most often pan traps and netting on flowers, or a combination of these two techniques. To exclude a possible bias in the wild bee species richness and abundance due to different sampling methods (as a categorical variable), we performed an ANOVA test. This test was performed to control if there were any statistically significant differences between the number of sampled species (richness) and individuals (abundance) across different sampling techniques. The results showed no biases (Appendix 4).
We performed correlation analysis (using Pearson's r) to reduce the number of bioclimatic and land-use variables, in PAST 3.04 (Paleontological Statistics Software Package) (Hammer et al. 2001). We decided to keep variables the least correlated between each other (p > 0.05) and possibly impactful on bee ecology and appropriate to discriminate different cities (Appendix 5 for full correlation table). The reason to exclude highly correlated variables is due to the possible emergence of concurvity between variables, possibly causing the narrowing of confidence intervals (Jiang et al. 2018). We selected the surface of the city (km 2 ), its population density (people/km 2 ) and mean altitude (m) as descriptors of the spatial and demographic composition of each city. To describe land-use cover, we selected normalized green edge density and the ratio between the green and impervious surfaces. The first has been shown to play a key role in shaping wild bee assemblages in urban environments (Theodorou et al. 2020a). The second was correlated with both green and impervious surfaces, and thus summed up well the overall land cover of the city. As bioclimatic variables, we chose annual mean temperature (°C*10) as in WorldClim (variable BIO1) and annual precipitation (mm) (variable BIO12). They overall distinguish cities from a climatic point of view. In addition, the temperature is a key factor in shaping wild bee communities since the highest diversity of wild bees is recorded in the Mediterranean biome and warm temperate environments (Ascher and Pickering 2020). The two chosen climatic variables were highly correlated with the rest of the climatic parameters, and the annual mean temperature was also linked with the latitude of the city. Taken together, these seven variables (city surface, population density, altitude, annual mean temperature, annual precipitation, normalized green edge density and green/impervious surface ratio) had a maximum r value of 0.4; thus, we assumed the absence of a strong correlation that would distort the following analyses. Finally, we included the number of sampling months as a control for different sampling efforts and the total number of sampled wild bees (i.e. total abundance, N), in an attempt to correct for possible biases induced by different durations of collecting activities and the amount of collected individuals. We graphically checked for the skewness of the distributions of the selected variables, and all were log-transformed to reduce their asymmetry.
To analyse how different functional traits were explained by the seven selected city characteristics, we performed generalized additive models (GAM), a semi-parametric expansion of generalized linear models (Rigby and Stasinopoulos 2005;Hastie 2017). This method has clear advantages over linear regression models. GAM does not assume a priori any form of relationships between dependent and explanatory variables, and thus can be used to reveal non-linear effects between them. This is particularly useful when dealing with land cover and climatic variables in predictive models (Virkkala et al. 2005;Feng et al. 2018;Falťan et al. 2020). The explanatory variables included the final selection of the seven land-use and climatic variables, plus total abundance (sample size) and sampling period length. The dependent variables, each tested with an individual GAM, were 19: species richness, the four diversity indices, the proportion of parasitic species (or individuals) over the total number of species (or individuals), the ratios between proportions of species (or individuals) with binary functional traits (oligolectic/polylectic, above-ground nester/below-ground nester, social/ solitary) and the proportion of species (or individuals) in each of the three body size ranks. For each dependent variable, the sub-optimal model was selected with a backward stepwise regression, removing one at a time the least significant explanatory variable (i.e. with the highest p value), until reaching the final model with the lowest value for the Akaike's information criterion (AIC) (Marra and Wood 2011). The goodness of the model was also tested by the increasing value of R 2 . For each predictor, a value of k = 5 was applied to be large enough to be reasonably sure of having enough degrees of freedom to represent the underlying smooth relationships. GAM analyses were performed in R version 4.0.3 using the mgcv package (Wood 2012). The complete dataset used in the analyses can be found in the supplementary file DATAset.xlsx.

Overview on cities and their wild bee communities
The 74 studied cases covered 55 different cities from 19 different nations and five continents (supplementary file DATAset.xlsx). The bulk of the studies has been performed in the northern hemisphere: 35 in North America and 29 in Europe (overall 86%). Only six studies were conducted in Central and South America, three in Asia and one in Oceania. The cities were very different from each other in terms of surface and population density. Almost all cities are below 500 m a.s.l. (Table I). The proportion between green and impervious surfaces varied considerably (Table I). Since most of the cities are in the temperate region, annual mean temperature and precipitation were more homogeneous than demographic and landuse variables (Table I).
Six out of the seven known bee families were represented in the analysed sample: Andrenidae, Apidae, Colletidae, Halictidae, Megachilidae and Melittidae (only the Australian endemic family Stenotritidae did not occur in the samples). Apidae resulted to be the richest family with 458 species, while Melittidae only had 11 species in the samples. Within all the families together, 158 different genera had been reported, the most abundant being Andrena with 212 species. Overall, a total of 1460 bee species have been reported in the analysed works (Appendix 2), of which most were small-to medium-size solitary, polylectic and below-ground-nesting species (Table II).

City size, demography and altitude effects
We did not find any significant effects of city surface on three of the four diversity indices considered (Shannon-Wiener diversity, taxonomic diversity and taxonomic distinctness) (Table III). Only Gini-Simpson dominance was significantly affected by city surface (Table III). This index was greater (and thus dominance lower) in cities with an intermediate surface of around 200-300 km 2 , while it seemed to decrease rapidly in very large cities (> 800 km 2 ) and, to less extent, in small ones (< 100 km 2 ) ( Figure 2A). The proportion of parasitic species ( Figure 2B) seemed to be negatively correlated with the city surface (Table III). The ratio between the proportions of oligolectic over polylectic species also seemed to be negatively correlated with the city surface (Table III, Figure 2C). In particular, small cities (around 10-50 km 2 ) seemed to host higher proportions of oligolectic species, with a decreasing trend in larger cities. However, from cities of ≥ 90 km 2 , such proportion seemed to remain constant with a slight increase in the biggest cities. Additionally, the proportion of mediumsized species seemed to be negatively affected by the city surface (Table III, Figure 2F). On the other hand, city surface has a positive effect on other functional traits. The proportion between above-and below-ground-nesting individuals seems to increase with city size, with larger cities (> 300 km 2 ) hosting fewer below-groundnesting bees (Table III, Figure 2D), even though this trend seems to be driven only by a few observations. Also, the proportion of large-sized species (Table III) and the proportion of individuals from medium-sized species (Table III) appeared to increase in larger cities ( Figure 2E, G).
Population density had a statistically significant effect on two tested variables. The number of large-sized bees seems to be strongly negatively affected by population density (Table III, Figure 3A), while a strong opposite trend emerged for the number of small-sized bees (Table III, Figure 3B).
Finally, altitude had a negative effect on the Gini-Simpson index (Table III, Figure 3C), the proportion of parasitic species (Table III, Figure 3D) and the proportion of medium-sized species (Table III, Figure 3F). In addition, the proportion of individuals from large-sized species seemed to be greater for intermediate values of altitude (Table III, Figure 3E), while the proportion of individuals from small-sized species had a wiggling trend, with a peak of small bees' abundances in high-elevation cities (Table III, Figure 5G).

Land-use effects
The ratio between the proportions of social and solitary species' abundances was positively affected by the edge density of green patches (Table III, Figure 4A); a similar effect was found on the proportion of individuals from large-sized species (Table III, Figure 4D). To a lesser extent, the proportion of medium-sized species (Figure 4E) increased in cities with medium to high values of edge density (Table III) and remained constant to its lowest values in cities with less fragmented green patches. Conversely, the ratio between the proportions of oligolectic and polylectic species was negatively affected (Table III) by the increasing level of fragmentation ( Figure 4B); the same trend is shown by the ratio between the proportions of above-ground and below-groundnesting species (Table III, Figure 4C).
The green over impervious surface ratio had more mixed effects on the tested variables. Though significant, its effect on the proportion of individuals from parasitic species seemed swinging (Table III, Figure 4F). Such proportion decreased in cities with much more to little more impervious than green surface, while a balanced ratio between both types of surfaces seemed to slightly positively affect the proportion of individuals from parasitic species. However, in cities with large green areas, the number of individuals from parasitic species seemed to remain quite Table II Descriptive statistics of the variables used to describe the taxonomic and functional diversity of the analysed wild bee communities, the sample size (total abundance) and the variable chosen to represent sampling effort (sampling period length). All wild bee functional traits are expressed as a proportion out of the total number of wild bees sampled in each study included in the analysis. S, number of species; N, number of individuals; Soc, social; Sol, solitary; Oli, oligolectic; Pol, polylectic; Ab, above-ground nester; Bel, belowground nester; large-sized, bees with body length ≥ 15 mm; medium-sized, bees with body length of 9-14 mm; small-sized, bees with body length ≤ 8 mm.   constant. Below-ground-nesting species are favoured in cities with more green than impervious surfaces (Table III). As a result, cities with smaller green patches are strongly associated with more above-than below-ground-nesting species ( Figure 4G). Finally, the proportion of individuals from small-sized species seemed to increase (Table III) in cities with more green than impervious surfaces ( Figure 4H).

Climatic effects
Temperature influenced particularly the diversity and the body size of wild bees. Both the species richness (Table III) and Shannon-Wiener diversity (Table III) decreased linearly as the annual mean temperature increases ( Figure 4A, B). The ratio between the proportions of oligolectic and polylectic species seemed to be negatively affected (Table III) by temperatures > 20 °C, but for low to medium mean annual temperatures, this functional trait followed an oscillatory pattern ( Figure 5C). The ratio between the proportions of above-and below-groundnesting species remained constant at temperatures ≤ 16 °C, and then it steeply increased (Table III) ( Figure 4D). The proportion of large-sized species was higher in warmer cities (Table III, Figure 5E), but the proportions of individuals from large-sized species (Table III) and from medium-sized species (Table III) were lower in warmer cities (Figure 5F and G respectively). Similarly, the proportion of individuals from medium-sized bees (Table III) decreased in warmer cities ( Figure 5H). Oligolectic species seemed to be favoured in terms of species richness at intermediate precipitation values (Table III, Figure 5I), while the ratio between the proportions of above-and below-ground-nesting species was lowest in both very dry and very rainy cities (Table III and shown in Figure 5J).

Sample size and sampling effort effects
The sample size (i.e. total abundance of sampled individuals, N) positively affected species richness (Table III, Figure S1A),  A, B, C, D, E) and the ratio between green and impervious surfaces (F, G, H) on different functional traits. Y-axis is the partial effect of the variable with 95% confidence intervals (grey shading). In brackets is reported the estimated degrees of freedom (e.d.f.) for each tested variable. All the explanatory variables were log-transformed. S: number of species, N: number of individuals, Soc: social, Sol: solitary, Oli: oligolectic, Pol: polylectic, Ab: above-ground nester, Bel: below-ground nester, parasitic: % of cuckoo bees, large-sized: bees with body length ≥ 15 mm, medium-sized: bees with body length of 9-14 mm, small-sized: bees with body length ≤ 8 mm. Shannon-Wiener diversity ( Figure S1B) and the proportion of parasitic species (Table III, Figure S1C). In contrast, sample size negatively affected the ratio between both the proportions of social and solitary species (Table III, Figure S1D) and the proportions of individuals from social and solitary species (Table III, Figure S1E). Finally, the proportion of individuals from both large-sized and small-sized species (Table III, Figure S1F) (Table III, Figure S1G) showed a wiggling pattern with increasing sampling size. The ratio between the proportions of above-and below-ground-nesting individuals seemed to be positively affected by the number of sampling months (Table III, Figure S1H). However, the number of sampling months seemed to be especially significant to explain part of the variation of body size. Indeed, the proportion of individuals from medium-and large-sized species, as well as the proportion of medium-sized species, all increased with increasing number of sampling months (Table III, Figure S1I-K). On the other hand, the proportion of individuals from small-sized species was negatively affected by the number of sampling months (Table III, Figure S1L).

DISCUSSION
Here, studying both taxonomic metrics (indices based on occurrences or abundances of species) and the diverse ecological traits related with the bee functional roles (Blondel 2003) allowed us to better evaluate how characteristics of cities influence the bee communities they host. Overall, our results suggest that variations in demographic, topological, land-use and climatic city traits have a variable role in shaping urban wild bee communities, as discussed in detail below.
Demographic traits of cities seemed to poorly affect their hosted bee communities, with only Gini-Simpson dominance being influenced by city surface. Instead, either variation in species richness and Shannon-Wiener diversity was affected by sample size or variation in both taxonomic diversity and taxonomic distinctness was not explain by any city trait, sample size or sampling effort. This important effect of sample size on species richness and Shannon-Wiener diversity is not surprising, since both are known to be strongly dependent on this factor (Konopiński 2020;Gotelli and Colwell 2001). On the other hand, Gini-Simpson index was greater in cities with an intermediate surface. Particularly in the larger cities of our dataset, communities tended to have a greater degree of dominance (lower Gini-Simpson index). This could be explained by the increasing presence of urban-adapting species (McKinney 2006) which may deplete most of the resources, possibly underpinning an increase in the homogenization of communities in larger cities (Ferenc et al. 2018). To less extent, the Gini-Simpson dominance decreased also in smaller cities. Cities with a reduced area could host a limited amount of resources (i.e. smaller green patches), hence causing again a possible impoverishment and homogenization of wild bee communities (Matteson et al. 2013).
Larger cities seemed to host only a few parasitic species when compared to smaller ones. In general, cuckoo bees seem to be a minor group in urban environments (Prendergast et al. 2022). Because cleptoparasitic species may play a stabilizing role in wild bee communities (Combes 1996), the lack of a rich and diverse cleptoparasitic component may underline a possible future disturbance of the whole community, since this guild is proposed as an indicator guild for the health of the entire wild bee assemblage (Sheffield et al. 2013). In addition, the largest cities seemed to host fewer oligolectic species than smaller ones, in accordance with several previous studies (e.g. Makinson et al. 2017or Lanner et al. 2020. Oligolectic bees heavily rely on the plant taxa they forage on (Praz et al. 2008) and bigger cities often host more exotic and ornamental plants than native ones (Threlfall et al. 2016), for example in botanical gardens or private properties. These plants could promote the abundances of generalist species rather than specialist ones. Hence, oligolectic bees could be less resilient to urbanization (Prendergast et al. 2022), since the narrow pollen preference may make these species more vulnerable to possible local extinctions (Buchholz and Egerer 2020). In turn, urban areas would favour generalist species (Garbuzov et al. 2017) that may outcompete oligolectic species for food resources. A loss of specialists could be deleterious since oligolectic bees are sometimes considered more effective pollinators than polylectic ones (Parker 1981) and could lead to a simplification of plant-pollinator networks in urban landscapes (Martins et al. 2017). Availability and distribution of nesting resources play a crucial role in shaping wild bee communities (Potts et al. 2005;Fortel et al. 2016), and we have found that largest cities host more above-ground-nesting species, perhaps because in large urban environments, man-made structures provide numerous cavities that could serve as nesting sites for above-ground nesters ) such as cavity renters and mud nesters (Häusler 2014). Bee size distribution within urban communities seemed also influenced by city surface, with richness of large species greater-and richness of medium-sized species lower-in larger cities. Banaszak-Cibicka and Zmihorski (2012) conducted a study in Poznań, a rather small city compared with bigger metropolis we included in our study, and indeed found the opposite pattern. Thus, the filter towards small-bodied bees could be caused by the property of the city rather than the increasing level of urbanization.
We have found urban-driven fragmentation to have a greater impact than the ratio between green and impervious surfaces on city bee communities. Bees are central place foragers: they depart from and return to their nests many times during their life (Charnov 1976). Social bee species tend to have a broad pollen diet spectrum (Danforth et al. 2019), limiting the effects of fragmentation since more workers can disperse in the environment to find food resources (Kaluza et al. 2017). Moreover, the maximum foraging range from their nests (Westrich 1996) is often larger for social than for solitary species (Gathmann and Tscharntke 2002), thus limiting the impact of fragmentation of foraging. Accordingly, we have found communities composed of many social species especially in cities with greater fragmentation. This hypothesis finds support in Banaszak-Cibicka et al. (2018) that found social species to be more frequent in highly impervious habitats compared to suburban sites.
In our dataset, highly fragmented green areas seemed to host fewer oligolectic bees probably due to a reduction in the quantity and quality of floral resources (Theodorou et al. 2020b). These highly fragmented habitats would probably be insufficient in hosting the specific plant species on which oligolectic species feed. Additionally, oligolectic species do not often have high dispersal abilities (López-Uribe et al. 2019). Consequently, small green fragments can hardly sustain specialist wild bees unless they host the required floral resources (Török et al. 2021). This, once again, could lead to a functional simplification of the communities and of the plant-pollinator networks (Martins et al. 2017). However, we cannot exclude that the reduction in oligolectic species could just be the consequence of the increase in social species (which are generally polylectic) in more fragmented cities. Cities with highly fragmented green areas seemed to host also larger bees. The ability of a bee to fly for longer time and farther from the nest is directly dependent on its body size (Gathmann and Tscharntke 2002;Steffan-Dewenter and Tscharntke 1997). Thus, fragmentation could favour larger bees due to their better dispersal ability, as they can more easily fly across patches to feed (Greenleaf et al. 2007). However, larger bees require more energy (Müller et al. 2006), and thus, there should be a tradeoff between within-patch food sources and flight range. Possibly, highly fragmented areas could disfavour small bees with a low dispersal ability (Wright et al. 2015) even though they have lower energy demands. In addition, habitat fragmentation reduces the possibility to find a mate within the green patch, especially for those species, often solitary, which mate on flowers (Paxton 2005). This may also help explain why social species are more abundant than solitary ones in urban highly fragmented habitats (Exeler et al. 2010).
More "green" cities seemed to host more below-ground-nesting bees compared with more impervious cities, in accordance with what was highlighted by Buchholz and Egerer (2020). An appropriate nesting substrate is a key factor for bees (Potts et al. 2005). Ground-nesting is a feature shared by more than 60% of non-parasitic bees (Cane and Neff 2011) and is probably the ancestral state for Apoidea. Cities with low levels of urbanization and more green areas are expected to host more ground-nesting bees (Tonietto et al. 2011) as these bees can find bare soil to dig the nest. However, bare soil was not a consistent class of land-use across the analysed cities, and thus, we cannot draw a specific conclusion about the co-occurrence of ground-nesting bees with barer soils. Results showed how an increase in impervious surface negatively correlated with the proportion of ground-nesting bees. Many ground-nesting bee species have been shown to prefer sandy soils (Antoine and Forrest 2021) and such a texture is typical of rural landscapes, becoming increasingly rarer in more urbanized habitats. However, increasing vegetation diversity, as it can be found in "greener cities", did not show an increase in nesting sites for ground-nesting bees in an agricultural landscape (Sardiñas et al. 2016). These mixed results could come from the fact that many aspects of the behavioural ecology of ground-nesting bees are still widely understudied and further investigations are certainly needed (Antoine and Forrest 2021). Additionally, cities with more green than impervious surface harbour more small bees. A reason could be that large green areas may provide bare soils for many small-sized bees as well as higher abundance of other nesting substrates, such as plant stems, typically used by small-sized bee species (Danforth et al. 2019).
Temperature, which largely affects several life-traits aspects such as the adult activity (Woods et al. 2005) and larval development rates (Forrest 2017), significantly affected diversity of the studied urban bee communities. Warmer cities showed a lesser richness and diversity than cooler cities, possibly reflecting the limits of thermal tolerance of bees (Maia-Silva et al. 2021). However, the upper thermal tolerance limit (CT MAX ) largely covaries with life-history traits in bees (Hamblin et al. 2017). Hotter cities seemed to host more above-ground-nesting species. In cities, cemented surfaces are hotter than vegetated lands (Herb et al. 2008), and thus, the combination of an increasing proportion of impervious surface, as we previously discussed, with higher temperature might overall reduce the abundance of ground-nesting species (as well as extreme dry or wet conditions) (Burdine and McCluney 2019). Hotter cities seemed also to favour large-bodied bees. An increase in body mass was shown to be positively correlated with higher upper critical temperatures (Oyen et al. 2016;Heinrich and Heinrich 1983). However, whether larger individuals are favoured in hotter climates is still under debate (Eggenberger et al. 2019). On the other hand, sociality in wild bees has been shown to be correlated with higher CT MAX (Hamblin et al. 2017). This could be the reason why hotter cities (as well as drier cities) seemed to host less oligolectic species as they are all solitary species (Danforth et al. 2019). However, we did not find temperature to affect the proportion between social and solitary species. Overall, hotter cities seem to host less diverse wild bee communities composed of large, social, polylectic and/or above-groundnesting species. These patterns are relevant from a conservation point of view since the mean temperature of cities around the world is greater than the surrounding areas (the "heat island" effect) and could even increase by over 2 °C in the near future (Bastin et al. 2019), jeopardizing the survival of wild bees in urban environments (Ayers and Rehan 2021).
In conclusion, our analysis overall suggests that, to some extent, some patterns (e.g. variation in lecty, nesting strategy and sociality) previously found in urban-natural and urban-rural clines can be observed in large-scale clines along gradients of city features. Hence, some city traits may amplify some patterns of variation in bee community found at local scale (from rural/natural to urban). On the other hand, the response of some bee functional traits (e.g. body size) to urbanization may be more variable at local scale depending on characteristics of the considered city. While patterns remain, thus, still difficult to generalize, actions to improve the colonization of cities by diverse bee functional groups are certainly needed, particularly in the form of nature-based solutions (Xie and Bulkeley 2020), such as increasing richer flower strips (Hofmann and Renner 2020) and green spaces (Antoine and Forrest 2021).