Applying predictive models to study the ecological properties of urban ecosystems: A case study in ZÃ1⁄4rich, Switzerland

Cities are human dominated ecosystems providing novel conditions for organisms. Research on urban biodiversity is rapidly increasing, yet it is still hampered by the partial spatial coverage of cities and because of existing taxonomic biases. Predictive models have proved to be a key tool to solve this shortfall. However, predictive models have rarely been used in urban ecosystems due to either the lack of sufficient species records or high-quality predictors (e.g. meaningful ecological maps). Here, we assemble a large cross-taxa inventory of 1446 species from 12 taxonomic groups, including several understudied invertebrate groups, sampled in 251 sites in Zürich, Switzerland. We investigate the species diversity distributions and the structure of species assemblages along artificial urban ecological gradients by applying predictive models. We find that the general species diversity distribution law, where assemblages are dominated by a few very abundant and frequent species, applied consistently across all taxonomic groups (3% of the species accounting for approximately 50% of abundance). Furthermore, only species of intermediate abundance and frequency are spatially structured along urban intensity gradients, with rare species numbers keeping constant even in the most urbanised parts of the city. In addition, we show that green areas with low mowing regimes are associated with higher species diversity in the majority of taxonomic groups. Hence, this suggests management relaxation as a low-cost solution to promote species richness. Our study demonstrates the potential of predictive modelling for addressing ecological questions in urban environments and to inform management and planning. * Corresponding author at: Conservation Biology, Biodiversity and Conservation Biology, Swiss Federal Institute for Forest, Snow and Landscape Research WSL, 8903 Birmensdorf, Switzerland. E-mail address: joan.casanelles@wsl.ch (J. Casanelles-Abella). 1 These two authors share senior authorship.

Cities are human dominated ecosystems providing novel conditions for organisms. Research on urban biodiversity is rapidly increasing, yet it is still hampered by the partial spatial coverage of cities and because of existing taxonomic biases. Predictive models have proved to be a key tool to solve this shortfall. However, predictive models have rarely been used in urban ecosystems due to either the lack of sufficient species records or high-quality predictors (e.g. meaningful ecological maps). Here, we assemble a large cross-taxa inventory of 1446 species from 12 taxonomic groups, including several understudied invertebrate groups, sampled in 251 sites in Zürich, Switzerland. We investigate the species diversity distributions and the structure of species assemblages along artificial urban ecological gradients by applying predictive models. We find that the general species diversity distribution law, where assemblages are dominated by a few very abundant and frequent species, applied consistently across all taxonomic groups (3% of the species accounting for approximately 50% of abundance). Furthermore, only species of intermediate abundance and frequency are spatially structured along urban intensity gradients, with rare species numbers keeping constant even in the most urbanised parts of the city. In addition, we show that green areas with low mowing regimes are associated with higher species diversity in the majority of taxonomic groups. Hence, this suggests management relaxation as a low-cost solution to promote species richness. Our study demonstrates the potential of predictive modelling for addressing ecological questions in urban environments and to inform management and planning.

Introduction
Urbanization is increasing worldwide, creating novel ecological conditions for biodiversity. Anthropogenic activities have influenced all ecosystems, culminating in the emergence of unprecedented systems around the globe (Boivin et al., 2016), of which urban ecosystems, such as cities, are hallmarks (Kowarik, 2011). Cities are social-ecological systems, where biophysical processes are integrated with human ones, and hence are composed of a fine-grained mosaic of grey, green and blue land covers  . Consequently, ecological patterns, processes and functions, such as species diversity, biotic interactions and biogeochemistry, are substantially altered compared with non-urban ecosystems, yet the same ecological laws operate (Fournier, Frey, & Moretti, 2020;Niemelä, 1999). Ecosystems are characterized by specific ecological properties, which include: structured patterns of diversity across sites, particularly dominance, frequency and rarity (MacArthur, 1965) as well as specific biotic and abiotic conditions that select for certain traits and species (Tansley, 1935) and form ecological gradients (Hawkins, 2001). Understanding how biodiversity is distributed along ecological gradients is essential to perform effective management (e.g. see Aronson et al., 2017). Thus, unravelling the ecological properties of urban ecosystems is a key step in order to make ecologically-based recommendations for urban biodiversity management.
Species diversity patterns are a main property of any ecosystem (MacArthur, 1965). Species abundance distributions across all types of ecosystems have been consistently found to be composed by a majority of scarce species (McGill et al., 2007). In addition, species abundance and site occurrence are known to be positively correlated, and thus rare species tend to be locally distributed. However, the specific diversity pattern might vary across ecosystem types, according to their specific properties. Changes in the amount of available habitat, disturbance and biotic interactions, together with population dynamics, translate to dominance and rarity shifts. For instance, intermediate levels of disturbance can boost the number of species that are rare and locally occurring, resulting in distributions with longer right tails (Magurran, 2004), that is, with a larger proportion of rare species (see Fig. 1). Conversely, increased disturbances dramatically raise the dominance of some species, leading to distributions with shorter right tails (Magurran & Phillip, 2008), that is, with a lower proportion of rare species (see Fig. 1). Finally, loss of habitat amount or increased stress also lead to a decline in rare species and also to distributions with shorter right tails (Magurran, 2004, see also Fig. 1).
Many urban ecosystems (i.e. cities, towns) are characterized by frequent disturbance, high levels of environmental stress and high degrees of habitat isolation. These urban ecosystem features are expected to represent strong environmental filters (Williams et al., 2009), through which only a subset of species from the regional species pool is able to passoften due to specific traits (Fournier et al., 2020), e.g. a broad feeding niche. Thus, a skewed abundance and site occurrence is expected, with species assemblages mainly composed of dominant and widespread species, with the total number of species largely dependent on the harshness of the environmental conditions. Rare species should be constrained to moderately urbanized areas where the effects of habitat isolation, stress and disturbance are lower. However, the effects of urban land cover are neither linear nor constant, varying both within and among cities (e.g. because of differences in habitat amount, management intensity or age, see Beninde, Veith, & Hochkirch, 2015;Sattler, Obrist, Duelli, & Moretti, 2011).
In fact, urban ecosystems are spatially heterogeneous, and often include strong environmental gradients (Fournier et al., 2020) along which biodiversity is distributed (McDonnell & Pickett, 1990). Patterns of species richness in non-urban ecosystems typically follow abiotic gradients of temperature and precipitation, which are particularly marked at large biogeographical scales (e.g. biomes). At smaller scales, however, a subtle combination of abiotic and biotic factors, including disturbances, play a role in shaping diversity (e.g. fire activity (He, Lamont, & Pausas, 2019), elevation (Rahbek, 1995), resource availability (Waldrop, Zak, Blackwood, Curtis, & Tilman, 2006)). Urban ecosystems are also characterized by such small-scale patterns, which are all imbued with human influence and can be expected to impact species diversity. Influential urban ecological gradients include the amount of available habitat, which is determined by the composition and configuration of land cover (e.g. Young & Jarvis, 2001); the disturbance regime (e.g. mowing, Smith, Broyles, Larzleer, & Fellowes 2014); the degree of environmental stress (e.g. heat island effect, Zumwald, Knüsel, Bresch, & Knutti, 2021); resource availability (e.g. Tew et al. 2021); and species interactions (e.g. altered competition, Ropars et al. 2019, mutualism, Harrison & Winfree, 2015, or herbivory, Just, Dale, Long, and Frank 2019. These urban gradients are likely associated with gradients in species diversity, which typically decreases with increasing urban intensity (see Blair, 1996Blair, , 1999Luck & Smallbone, 2013;Sol, Bartomeus, González-Lagos, & Pavoine, 2017; but see also Guenat, Kunin, Dougill, & Dallimer, 2019). These underlying ecological drivers are usually linked to warmer, drier and more polluted conditions, and to fewer and more isolated available habitats. Nonetheless, few studies have formally tested the existence of continuous diversity gradients across taxa and urban land cover types (McDonnell & Hahs, 2013).
As in other fields in ecology, urban research suffers from the incompleteness of species richness data, both taxonomically and spatially. First, the taxonomic coverage is fragmented and biased towards some groups, especially birds, with many taxa remaining understudied (McIntyre, 2006). Second, biodiversity has mostly been assessed at the sampling point scale within a single subset of existing land covers, and analyses of the full panel of urban land cover types are lacking. An improved understanding of which factors drive urban biodiversity and its distribution is therefore needed. Multi-taxa datasets (Hortal et al., 2015) with high-resolution spatial predictor variables describing abiotic and biotic urban gradients are beginning to be available and will likely prove valuable in this regard. Integrating such resources with algorithms commonly used for Species Distribution Models (SDMs) might open interesting prospects in predicting and highlighting biodiversity patterns associated with fine-grained ecological gradients of urban ecosystems impossible to obtain using data from only explicit sampling points. Additionally, such analyses across multiple taxonomic groups are expected to provide essential information to better characterise the ecological properties of urban biodiversity at the Three examples of potential ecological processes are given: disturbance regime, habitat amount and stress. For each example, the specific color legend is depicted. Note that the three exemplary processes are independent and do not represent all possible processes affecting species abundance distributions. Note also that the shapes of the curves are a simplification based on Bazzaz (1975), Maguran (2005) and Kempton (1979). citywide level, which can in turn enhance biodiversity management and urban planning.
Here, using predictive modelling we investigate the ecological properties of an urban ecosystem and explore whether such properties could inform urban biodiversity management. We use a unique species dataset collected during the last decade across 251 locations within the city of Zürich, Switzerland. The dataset is composed of 1446 species belonging to birds and invertebrate groups. We examine biodiversity patterns based on modelled and predicted citywide distributions and ask the following questions: (i) Do urban ecosystems exhibit ecological organisation equivalent to those found in other ecosystems? In particular, are urban communities composed of a large number of rare and local species and few very abundant and widespread ones? (ii) Does biodiversity change along ecological gradients and, if so, what factors generate the gradients?; and (iii) what management and planning recommendations can be derived from the observed biodiversity patterns and predicted citywide distribution of urban biodiversity? We also discuss the potential of predictive models for ecology, urban planning and society.

Study area
The study area is the city of Zürich, Switzerland (47 • 22 ′ 0 ′′ N, 8 • 33 ′ 0 ′′ E). The municipality of Zürich covers approximately 92 km 2 and has a population of over 400,000 inhabitants, corresponding to the third quartile of city population size in Europe (Fournier et al., 2020). The city's climate is temperate, having an annual mean precipitation of 1,134 mm and mean temperature of 9.3 • C.

Bird and invertebrate sampling
We used species occurrence records and abundances collected by  Table S1. Combined, these datasets contain records of urban fauna from 251 sampling sites, corresponding to widely distributed types of urban green spaces in the world , that is, parks, allotments, gardens, brownfields and green roofs (see Section "2.2.2 Urban green space types" in Fournier et al., 2020). Sampling sites had a minimum area of around 20 m 2 and were separated on average by ca. 4000 m, with a minimum distance of 100 m and a maximum distance of 11000 m (see Fournier et al., 2020, for details). Moreover, these datasets contain records of over 1500 species of both birds and invertebrate groups (Table S1), representing one of the largest datasets sampled in an urban ecosystem.
Bird species were sampled using the point count method in the early morning during the breeding season that spanned between April and June 2007 (Fontana et al., 2010) and April and June 2019 (Moretti et al. unpublished data). Each bird sampling site was visited six times in Fontana et al., 2011 and two times in Moretti et al. (unpublished data) alternating the visiting hour, which ranged between 1 h before sunrise and 5 h after sunrise (Fontana et al., 2011).
Invertebrate species were sampled using activity traps, including pan traps and pitfall traps (see details in Frey et al., 2016;Sattler, Duelli, et al., 2010;Tresch et al., 2019) and additionally using trap-nests for cavity-nesting bee and wasp species. Pan traps were used to sample air-dispersing arthropods (e.g. bees, wasps, hoverflies, beetles). Two types of pan traps were used. First, in parks, brownfields and green roofs, pan traps consisted of a window interception trap made of Plexiglas with a yellow pan trap (diameter 44 cm). In both allotments and gardens, pan traps consisted in three 1-liter bowl traps (diameter 20 cm without an interception window), each coloured in UV-bright blue, white or yellow paint. All pan traps were placed on wooden poles at around 1.5 m height. Pitfall traps were used to sample ground invertebrates (e.g. ground beetles, spiders, snails). Pitfall traps were placed in parks, brownfields, green roofs, allotments and gardens. Particularly, each sampling location had a set of three pitfall traps of 72 mm diameter, covered with transparent roofs at 10 cm. Both pan and pitfall traps were emptied on a weekly basis (see details in Table S1). For cavity-nesting bees and wasps, we installed a trap-nest at each sampling location following Staab, Pufal, Tscharntke, and Klein (2018) in allotments, gardens and parks. Particularly, trap-nests consisted in three 12 cm diameter plastic tubes. The first two were assembled using 200-300 reeds of Phragmites australis (Cav.) Trin. with diameters ranging from 1 to 10 mm and a length of 22 cm and 5-10 bamboo reeds to cover the whole requirements of the cavity-nesting bee and wasp community. The third trap-nest was assembled only using cardboard tubes of 7.5 mm specific for Osmia spp. In allotments and gardens, trap-nests were placed in wooden pools at around 1.5 m height. In parks, trap-nests were installed in vertical structures (mainly trees and occasionally light poles) at a height of 2.5-3.5 m. In all cases, trap-nests were exposed S-E and with direct sun light. The trap-nests were placed from January until October 2016 (Frey et al., 2016) and 2018 (Casanelles-Abella et al. unpublished data).

Environmental data
The distribution of individual species, as well as the richness of taxonomic groups, is typically modelled using environmental variables representing biotic and abiotic conditions relevant for the specific taxa (Chauvier et al., 2020;Descombes, Pradervand, Golay, Guisan, & Pellissier, 2016). Urban ecosystems do have abiotic and biotic conditions, but in contrast with other ecosystem types, these are contained within the human dimension, and hence urban ecosystems are considered social-ecological systems (Alberti, 2015;Des Roches et al., 2020). Following established procedures in macroecology and conservation (Guisan, Thuiller, & Zimmermann, 2017;Pearson, 2010), but also aiming to integrate the unique social-ecological properties of cities, we gathered available high-resolution predictors as proxies of the urban intensity. The predictors represented (i) climate, (ii) pollution, (iii) vegetation structure, and (iv) urban land cover to be tested as candidate variables. These are four main proxies of the social-ecological conditions of cities and specifically of urban intensity, measuring environmental stress, disturbance, habitat heterogeneity and amount of available habitat Frey et al., 2018; for instance, see Sattler, Borcard, et al., 2010), and thus expected to be important determinants of urban species distributions. Specifically: • Climate. Cities are known to have higher air temperatures than the surrounding landscape because of the high proportion of impervious surfaces and artificial materials that retain heat during the day and release it at night. This phenomenon is known as the urban heat island effect, and its eco-evolutionary effects are starting to become apparent (see Fournier et al., 2020 andPiano et al., 2017 for examples on the effects of overwarming on community assembly of bees and ground beetles). While increased temperatures in cities might benefit species with a relatively high heat tolerance, phenotypic plasticity or pre-adaptation, they might also have detrimental effects on species not adapted to hot and dry conditions (Fournier et al., 2020). In addition, altered temperatures together with the building landscape might have an effect on wind currents. This might affect the mobility of certain animals, such as birds or flying insects. Consequently, we included two climatic variables: (i) a model of local overwarming from 2010 (Parlow, Fehrenbach, & Scherer, 2010). Specifically, overwarming (i.e. overheating) measures the increased heat load of air temperature (above the ambient temperature) due to heat exchange with urban surfaces, and (ii) a model of daily average wind speed above a height of 50 m from the ground surface from 2014 (Kanton Zürich, 2019). • Pollution. Vehicles produce a large amount of urban pollution in the form of exhaust (e.g. NO 2 , CO, SO 2 ), which can negatively impact several urban species and have cascading effects. We used two models of NO 2 , one from 2010 (Parlow et al., 2010) and one from 2015 (Kanton Zürich, 2019), and one of particulate matter (PM) from 2015 (Parlow et al., 2010) that modelled the concentration of the exhaust particles. • Vegetation structure. We used a publicly available light detection and ranging (LiDAR) remote sensing dataset to describe two main dimensions of vegetation structure: height and cover. These attributes affect the quality and amount of habitat available to our target species (Frey et al., 2018), e.g. via the provisioning of suitable microclimatic conditions (Zellweger, De Frenne, Lenoir, Rocchini, & Coomes, 2019). Vegetation height was defined as the 95th percentile of all raw LiDAR points classified as vegetation, and vegetation cover was defined as vegetation return heights of greater than 1 m above ground. Layers were acquired between March and April 2014, with an average density of 8 points per m 2 , a footprint size of 0.2 m, and a vertical accuracy of 0.1 m. The raw data have the form of a classified point cloud including the categories buildings, ground points and vegetation points. • Land cover. The available land cover cartography consisted of the urban habitats of Zürich developed by , included in a very detailed raster layer. It included: (i) the habitat map of the city of Zürich at 1 m resolution, and (ii) the detailed layers of buildings and streets from the Swiss Federal Office of Topography (see  for details). In total, 20 land cover classes were described (i.e. 1 for water bodies, 8 for grey structures and 11 for urban green areas) at high resolution, accounting for the precise composition of habitats within the city that have proved to be important predictors of urban biodiversity . The 20 classes were merged into 7 (for modelling) and 6 (for plotting) categories summarizing major grey and green land cover types. Specifically, all grey classes that were not buildings were combined together (Grey and Other grey surfaces, see Table S6). Green land cover aggregation was done following Aronson et al. (2017, particularly as described in WebTable 1; see also Table S6). Particularly, green land covers were aggregated based on their (i) stakeholders, (ii) native species input and (iii) management regimes of the herbaceous vegetation (see Table S6), which in most temperate cities such as Zürich is mainly done via mowing (Smith et al., 2014). We computed two types of landscape metrics on the different categories. First, each new category was then converted and aggregated by sum to five distinct raster layers, from 1 m to 10, 20, 50, 100 and 200 m resolution to account for local and landscape effects. Second, Euclidean distances within each layer were calculated using GDAL (https://gdal.org/). All water bodies, annual crops and forest cover located in the hills surrounding the city were excluded to improve model projections.

Species diversity distributions
2.5.1. Classification of species according to rank-abundance and rankoccurrence distributions We studied the frequency and dominance patterns of the different taxonomic groups separately (Fig. S3-S4). For dominance, we constructed rank-abundance diagrams and classified species according to the proportion of abundance they contained. Specifically, a species was considered abundant for a given taxonomic group if it was contained between the upper 90th (Gaston, 1994) and the lower 50th (Ter Steege et al., 2013) percentile of the rank-abundance diagram. In addition, species that appeared above the 50th percentile of the rank-abundance diagram were classified as hyperabundant. The remaining species were classified as rare. Similarly, for frequency we created rankoccurrence diagrams and classified species according to their proportion of occurrence (i.e. the proportion of sites where the species was recorded). In this case, a species of a given group was considered widespread if it was contained between the 75th (Gaston, 1994) and 25th percentile of the rank-occurrence diagram, and hyperwidespread if it belonged to the upper 25th percentile of the rank-occurrence diagram. The remaining species were classified as local. As site occurrences and abundances were positively correlated (Table S2), we used three general categories: hypercommon, common and rare.

Selection of models to study species abundance distribution
Currently, several models are available to describe species abundance distribution using rank-abundance diagrams, and related fits depend mainly on the evenness of the data and the sampling intensity (see Slik et al., 2015, for a detailed discussion). To incorporate this range of possibilities, we fitted our species abundance and occurrence distributions to four main model types, that is, log series, log normal, broken stick and Pareto (power law) distributions. The Pareto (power law) model consistently provided the best fit of the data for all taxonomic groups. We tested the goodness of the fit using two common procedures (Slik et al., 2015): (i) using maximum likelihood tools to rank the models according to Akaike's Information Criterion (AIC), and (ii) graphical exploration using Whittaker (rank-abundance diagrams) plots. The Pareto model always had the lowest AIC (Table S3) and visually fitted the observed data best (Fig. S1). All analyses were conducted in the R environment (version 3.6.1 (R Core Team, 2019)) using the package sads version 0.4.2 (Prado, Miranda, & Chalom, 2013).

Species records
We obtained presence and absence data for each species of the 12 selected taxonomic groups. The total number of observations per species depended on the number of datasets it was sampled in (i.e. if a given taxonomic group was studied in two datasets, but a species of that taxonomic group was only sampled in one dataset, the number of observations of that species was calculated only in the dataset it was sampled in). For each species, presences and absences were each reduced to one record per sampling site.

Response data for species distribution models and species richness models
To ensure better model distribution performances, we performed two types of predictive modelling: (i) Species Distribution Models (SDMs), and (ii) Species Richness Models (SRMs). (i) SDMs allowed each species to be modelled individually. We filtered the species list to include only those that occurred in a minimum of 30 sampling sites and that had at least 30% of the absences. The filtered dataset included 272 species. These species mostly belonged to the dominant and frequent groups within rank-abundance and rank-occurrence distributions (86% of coincidence with both rank-abundance and rank-occurrence classifications, 93% coincidence with rank-occurrence alone; Table S4 & Data file S1).
(ii) SRMs allowed overall species diversity to be modelled i.e. accounting also for rare species not included in SDMs. For that, we computed raw species richness at each sampling location and for each of the 12 taxonomic groups including three distinct species richness sets: (i) including all species; (ii) including only species with a low proportion of presences, i.e. mostly rare species according to both rank-abundance and rank-occurrence distribution; (iii) including only species with a low proportion of absences, i.e. mostly hypercommon species. In total, 36 species richness sets (12 taxonomic groups * 3 species sets) were modelled using SRMs.

Variable selection
Variable selection was done separately four times (i.e. for the common species SDMs, and the hypercommon, rare and all species SRMs) and followed a 2-step procedure following existing literature (Descombes et al., 2016;Guisan & Zimmermann, 2000). First, we preselected a large set of predictors based on their statistical relevance measured using the predictive power (D 2 ). Second, we manually picked six key variables of ecological relevance (i.e. to infer environmental stress, disturbance, habitat heterogeneity and amount of available habitat), as well as inter-correlations < 0.7 (Dormann et al., 2013). Each candidate list of predictors was balanced across the four proxies, were kept at or converted to a 10 m spatial resolution and projected to the standard Lambert projection (EPSG:2056). More details on the final selection of candidate variables from the available climate, pollution, vegetation structure and urban land cover variables can be found in Table S7.

Model calibration, evaluation and ensemble
Depending on the type of model (i.e. SDMs or SRMs), two distinct probability distributions were used: Binomial for SDMs and Poisson for SRMs. Each SDM and SRM was calibrated using an ensemble of four common modelling techniques to account for model uncertainty and specificity (Buisson, Thuiller, Casajus, Lek, & Grenouillet, 2010): Generalized Linear Models (GLM), Generalized Additive Models (GAM), Random Forests (RF) and Gradient Boosting Machines (GBM). GLMs and GAMs are models based on linear regression. While GLMs assume a parametric relationship between response and predictors, GAMs focus on flexible nonparametric smoothing functions. RFs and GBMs are defined as tree-based models and show a higher complexity in their response than GLMs and GAMs. Each modelling technique was parameterized in the following way: GLMs were calibrated with second-order polynomials, GAMs with a spline smoothing term of intermediate complexity (k = 3), RF with a node size of 5 (nodesize = 5) and 1000 trees, and GBM with an interaction depth of 1, a shrinkage of 0.001 and 1000 trees. In addition, we set GBM shrinkage at 0.001. The models were computed using the R packages mgcv version 1.8-30 (Wood, 2001), RandomForest version 4.6-14 (Liaw & Wiener, 2018), and gbm version 2.1.5 (Ridgeway, 2007).
For each SDM and SRM, species records were split randomly into two sets containing 80% and 20% of the data. The former was used to calibrate the model and the latter for evaluation. This procedure was replicated five times. Model performance was assessed with True Skill Statistic (TSS, Allouche, Tsoar, & Kadmon, 2006) for SDMs, and with Pseudo-R 2 (see Thomas et al., 2018) for SRMs. TSS evaluates model skill in distinguishing absences from presences. Pseudo-R 2 provides a measure of predictive performance by determining the ratio between model error and variance of the response variable (Thomas et al., 2018). Both performance metrics range from 0 to 1, with 1 indicating perfect models. Models were then filtered according to their predictive performances. For SDMs, model predictive performance was considered reliable when TSS > 0.4, a commonly used minimum threshold (Thuiller, Guéguen, Renaud, Karger, & Zimmermann, 2019). For SRMs, the quartile distribution of performance metrics was calculated and models with the 25% worst performance were removed.
For each of the 272 species retained SDMs among the initial 20 (5 repetitions * 4 algorithms) were projected over the study area. Each model prediction was then converted into binary output using the value that yielded the maximum TSS as a threshold. Binary layers of presences/absences (PA) were stacked and the final species PA layer was formed by applying a threshold of 50%, above which cells were assigned species presence (Araújo & New, 2007). Finally, we combined the distribution maps of the 272 species into the respective taxonomic groups to form group richness maps of common species.

Citywide predictions of the species richness models
Citywide maps of predicted species richness were then created for the 12 taxonomic groups, considering all species, as well as hypercommon species only and rare species only (Fig. S4). We stacked the species richness maps of the 12 groups to generate a single citywide map for overall species richness (Fig. 4 C).
Additional methodological information on model calibration, validation, ensemble and predictions can be found in Table S7. Detailed information on model performance is given in Data files S1-S2.

Predicted species richness per land cover type and city district
We obtained the predicted species richness for the six main categories of land cover (i.e. buildings, human-made surfaces, urban woody patches, gardens, green areas with high management regimes and green areas with low management regimes; Fig. 4, see also Table S6 for the definitions of land covers and the specification on management regimes), for the overall biodiversity (Table S5) and for each taxonomic group (Data file S4). To analyse the diversity distribution within the city and to discuss potential management and planning strategies, we classified the city quartiers (i.e. the minimum administrative unit of the city of Zürich, see Stadt Zürich, 2020a) into four main regions. The four regions were defined based on the urban history and on the available metrics of recent urban development, that is, the proportion of new build lodgements and population density between 2011 and 2020 (see Stadt Zürich, 2020a). They included (i) the old townthe oldest region of the city and currently not undergoing major changes (district 1) -(ii) the former industrial quartier (districts 4-5) -slightly less old than the old town, where industries used to be located, and undergoing rapid urban development -(iii) the peripheral districts undergoing strong urban development (i.e. where the proportion of new build apartments between 2011 and 2020 was larger than the city average of 10.9%, comprising districts 9, 11 and 12; Stadt Zürich, 2020a), and (iv) the remaining peripheral districts (i.e. where the proportion of new build apartments between 2011 and 2020 was smaller than the city average of 10.9%, comprising districts 2, 3, 6, 7, 8, 10; Stadt Zürich, 2020a). In each region, we calculated the area and proportion of the six land cover categories.

Species diversity distributions
Our results demonstrate that the universal diversity distribution patterns of non-urban ecosystems also occur in urban ecosystems across different taxonomic groups. We show that urban ecosystems follow the abundance-distribution relationship, where rank-abundance and rank-occurrence diagrams rendered a similar classification of species (Pearson correlation = 0.75 ± 0.13 and 89% identical classification; Fig. 2 and Table S2) and fitted a Pareto distribution (Fig. S1 and Table S3). A small percentage of the species (206 spp.) were highly abundant, contributing up to 90% of the total abundance (Fig. 2 A). Among these, 27 species were disproportionately abundant (i.e. hyperabundant species; Fig. 2 A), representing half of the sampled individuals (mean ± sd = 2921 ± 2121, min = 1307, max = 10,200; Table S4).
A small percentage of species (328 spp., 22.4%) accumulated 75% of the recorded occurrences (Fig. 2 B). Among these, 55 species were recorded on most of their study sites (mean ± sd = 82 ± 10% of the sites; Table S4), which indicates that they are distributed over almost the entire city (i.e. hyperwidespread species; Fig. 2 B). In contrast, 89% of the species (1257 spp.) had low abundances and represented only 10% of the sampled individuals (mean ± sd = 13 ± 22, min = 1, max = 117; Table S4), indicating that urban ecosystems are composed of an extremely large number of rare and scarce species (Fig. 2 A). Likewise, the majority of species recorded (1136 spp., 77.6%) were locally occurring and therefore were recorded solely in a small fraction of the study sites where they were sampled (mean ± sd = 4 ± 3%, min = 1%, max = 16%; Table S4). Moreover, the proportion of rare and locallyoccurring species in each sampling site was generally low. Because the vast majority of species were not widespread, our results suggest an important taxonomic turnover among urban patches. We found some groups (bees, spiders, true bugs) to have a particularly large number of both hyperabundant and hyperwidespread species, which supports the idea that urban ecosystems might especially promote certain taxa (Baldock et al., 2019;Concepción, Moretti, Altermatt, Nobis, & Obrist, 2015;Theodorou et al., 2020). Hereafter we use three general categories (see Methods and Data file S1): 'hypercommon' for hyperabundant and hyperwidespread species (46 spp.), 'common' for widespread and abundant species (272 spp.), and 'rare' for local and rare species (1128 spp.).

Biodiversity gradients in cities
We found a prominent biodiversity gradient for common species and urban intensity, which was consistent for the vast majority of taxonomic groups and well represented by the variable local overwarming (Fig. 2), one of the urban intensity proxies used (see section 2.4 in Methods). Local overwarming measures the degrees of overheating within the city, mainly due to impervious surface cover, and is highly correlated with other urban intensity proxies (e.g. pollution, land cover composition, vegetation heterogeneity, see Fig. 4). Species richness followed a negative nonlinear relationship with local overwarming that was consistent and widespread across most of the taxonomic groups (Fig. 3). In particular, there was a marked decline in the richness of common species at around 1 • C of overwarming, with pronounced declines in the richness of bees, beetles, spiders and true bugs (Fig. 3). Most of the city has values higher than 1 • C of overwarming (Fig. 4 C). In contrast, species richness of rare groups did not follow a detectable ecological gradient and their richness kept relatively constant along the local overwarming gradient, suggesting a high taxonomic turnover given that rare species represent the vast majority of the species pool (Fig. 2). Conversely, the richness of hypercommon species showed a weak (true bugs, other beetles) or no relationship with local overwarming (Fig. 3).

Fig. 2.
Dominance and frequency patterns of the studied taxa. Rankabundance (A) and rank-occurrence diagrams (B), and the abundance distribution relationship (C) of the 1446 species from 12 taxonomic groups sampled at 251 study sites (D). For the rank-abundance diagram (A), species are classified into the categories hyperabundant (dark blue), abundant (green) and rare (light grey). For the rankoccurrence diagram (B), species are classified into the categories hyperwidespread (dark blue), widespread (green) and local (light grey) for frequency. Species classification was done separately for each taxonomic group (for details, see Methods and Data file S1). Each bar in (A) and (B), and each dot in (C) corresponds to one species. Each dot in (D) corresponds to one sampling location, coloured according to taxonomic groups sampled. Additional information can be found in Table S2-S4 and Data file S1. Rankabundance and rank-occurrence diagrams for each individual taxonomic group are shown in Figs. S3-S4. This lack of relationship was also observed for several of the other proxies of urban intensity, such as air pollution and land cover metrics.

Citywide distribution of urban biodiversity
We found major differences in species diversity across the investigated land cover types offering management opportunities. Our estimates indicate that the highest species richness appears in green areas with low management regime (mean ± sd = 190 ± 23, min = 104, max = 229; Fig. 5 and Table S4), e.g. meadows, pastures and ruderal patches, all of which are characterized by a lower vegetation management frequency and intensity than other green land covers, particularly concerning mowing regimes (see also Table S6). The other types of green areas had similar numbers of predicted species, including urban woody patches (mean ± sd = 167 ± 24, min = 104, max = 227; Fig. 4 and Table S5), green areas with high management regime (mean ± sd = 167 ± 18, min = 105, max = 225; Fig. 5 and Table S5) and gardens (mean ± sd = 165 ± 17, min = 103, max = 230; Fig. 5 and Table S5). Green areas with low and high management regime covered a similar percentage of the city (ca. 18% and ca. 20%, respectively).
Although green areas with low management regimes consistently had a higher richness for most taxonomic groups, our results show that differences may occur among taxonomic groups. Specifically, we found millipede and snail richness to peak in gardens and wasp richness to peak both in gardens and on green areas with high management regimes (Fig. S4). Furthermore, our results demonstrate an unequal distribution of green land covers and consequently of biodiversity within the city, with higher values in the peripheral districts (ca. 50%; Fig. 5) than in the core ones (ca. 14% and 25%; Fig. 5).

Ecological properties of urban ecosystems
Urban ecosystems display ecological properties, in terms of diversity distributions patterns and spatial structure along ecological gradients, equivalent to those of other ecosystems and can harbour diverse species assemblages. We tested for our 12 studied taxonomic groups two main ecological rules and found consistently that (1) the vast majority of individuals sampled within and across sampling sites belong to a few species, and (2) that highly abundant species are also the most widespread. Thus, our results demonstrate that one of the few universal ecological laws (McGill et al., 2007) also applies in cities (Smith, Warren, Thompson, & Gaston, 2006). Particularly, the recovered urban diversity patterns differ from highly managed (Bazzaz, 1975) or stressed (Kempton, 1979) ecosystems composed mostly by few dominant species. While the subset of abundant species has been frequently assumed to account for almost all urban biodiversity (McKinney, 2006), our results

Fig. 3.
Predicted species richness along a gradient of urbanization. Response curves of the predicted species richness along a gradient of urbanization, inferred using local overwarming as a proxy (urban heat island effect) for (A) hypercommon, (B) common and (C) rare species of the 12 taxonomic groups. Bands represent the 95% confidence interval.
show that the majority of urban biodiversity is composed of scarce, locally occurring species distributed all over the urban intensity gradients. Thus, a large fraction of the urban species pool is distributed locally and might be more sensitive to patch-scale habitat features (Beninde et al., 2015) or stochasticity (Sattler, Borcard, et al., 2010) rather than citywide social-ecological gradients. Urban ecosystems are composed of a mosaic of small-sized land cover types subjected to large numbers of small-scale management decisions. This can satisfy locally the habitat requirements of several species through, for instance, habitat supplementation and complementation (Colding, 2007), but also lead to a considerable contribution of stochasticity, ecological drift and extinction and colonisation events. Several studies have investigated the diversity distributions and patterns of specific groups in cities, finding that cities can promote some taxonomic groups (Alvey, 2006;Baldock et al., 2015;Sattler et al., 2011;Smith et al., 2006), but not others (for instance, see Aronson et al., 2014;Theodorou et al., 2020). Here, we show that species diversity distribution rules apply to our studied groups regardless of their taxonomic identity and their functional features and roles.
By applying predictive models and obtaining citywide species richness distributions, we showed that diversity was spatially structured along urban intensity gradients. In any ecosystem, biotic and abiotic factors structure the spatial distribution of species, which is particularly noticeable with pronounced gradients such as tidal, alluvial or fireprone and also urban ecosystems. However, environmental gradients in urban ecosystems are distinct as they are of social-ecological nature (Avolio, Pataki, Trammell, & Endter-Wada, 2018;Des Roches et al., 2020;Rivkin et al., 2019), being spatially heterogeneous and dynamic (Ramalho & Hobbs, 2012) with potential consequences on the spatial distribution of diversity. We found that only intermediately abundant and common species are spatially distributed along urban intensity gradients and that their richness follows a nonlinear negative relationship with urban intensity, represented in Fig. 3 as local overwarming. The declines of several groups (e.g. bees, spiders, true bugs, beetles) at relatively low urban intensity values (i.e. with relatively low values of streets and high values of available habitat) are concerning. Common species might encompass species displaying a broad niche and a mobility degree sufficient to access and use several habitat types (Concepción et al., 2015;Fournier et al., 2020), yet still having important constraints. Whether the strong declines are caused by actual overwarming and thus thermal stress or by one or more different proxies of urban intensity related to other urban biodiversity drivers (e.g. amount of available habitat, disturbance) is unclear and deserves future attention to better inform biodiversity management.
On the other hand, the predicted number of rare species is relatively low and constant in the whole city even though they represent most of our species. This challenges the idea that urban intensity gradients uniformly filter overall diversity across taxonomic groups (Knop, 2016;McKinney, 2006). Although this filter clearly exists, its outcomes are modulated by city-specific factors with functional traits mediating the sensitivity of each species to urban intensity (Fournier et al., 2020). Species with a small distribution may be mostly sensitive to processes occurring at the patch scale or to stochastic population dynamics . Finally, as expected, hypercommon species are distributed all over the city, suggesting that they represent an extreme case of species able to exploit the vast majority of the urban environmental conditions. Although our results apply to the city of Zürich, we expect our patterns to be generalisable in other European cities, particularly those with similar social-ecological features (e.g. size, composition of urban land covers, management types). Prior evidence has already demonstrated the existence of rich urban species assemblages (e.g. bees, plants, Baldock et al., 2019;Kühn, Brandl, & Klotz, 2004;Sattler et al., 2011), which might apply to other understudied groups.

Opportunities for biodiversity management and urban planning
Predicted patterns of hypercommon, common and rare species Fig. 4. Local overwarming. Flat violin plots (Allen et al., 2019) and boxplots, each calculated from 2000 randomly selected points, showing the relationship between the local overwarming and the amount of grey surface in a radius of 100 m (A). Correlation between the local overwarming values and 10 other environmental predictors used as proxies to measure the urbanization gradient, including pollution (PM = particulated matter, NO 2 = nitrogen dioxide), remote sensing variables (LiDAR p95, NDVI) and land cover metrics (amount of grey cover and green cover in different buffers) (B). Histogram of the local overwarming in the city of Zurich, from 0 (no overwarming) until 5 • C (C). Map of the local overwarming (number of degrees above the calibrated baseline temperature) distribution in the city of Zürich. The baseline temperature (0 • C) corresponds to the ambient air temperature. Negative values indicate temperatures below the baseline, and correspond to the forested hills (D).
richness have implications for management. Regarding common species, many groups display a higher species richness at low values of local overwarming. Provided that a large part of Zürich has higher values of overwarming than the threshold observed (see Fig. 3 B), this poses important questions on how urban planning should promote common species, specifically regarding the "urban compacting" versus "urban sprawling" debate (i.e. the land sparing and land sharing debate in cities, see Geschke, James, Bennett, & Nimmo, 2018). Amelioration efforts in areas of high urban intensity are a possible strategy and of particular benefit for people. However, preventing densification of less urbanized areas might benefit better common species by maintaining the amount of available habitat (see for instance Geschke et al., 2018) (Allen et al., 2019) and boxplots of the predicted number of species in the six main types of land cover (i.e. buildings, other grey surfaces, gardens, green areas with high management regime and green areas with low management regime). (B) Land-cover composition in four areas of the city: the core districts represented by the former industrial quartier currently containing the main railway network ('Former Industrial') and the old town ('Old Town'); and the peripheral districts divided into those undergoing and important urban development ('Peripheral develop.') and the remaining peripheral districts ('Periphery'). (C) Combined predicted species richness map of the 12 taxonomic groups included in the study. Predicted species richness was modelled with the Species Richness Model (SRM) for each group using all available species per group. White areas correspond to water bodies, annual crops and forest that were not included in the study. See Fig. S2 for detailed predicted species richness maps of each of the 12 taxonomic groups and Fig. S4 for the flat violin plots and boxplots. species likely depend more on the individual local-scale decisions of stakeholders. They might be less suitable target of large-scale conservation management initiatives, and therefore substantially more sensitive to extinctions. Areas of high urban intensity, that negatively affect common species, might be still able to sustain at least temporarily, in specific locations with particular local conditions, a diverse and spatially dynamic community of rare species.
Green areas with low management regimes have the highest predicted species richness, while green areas with high management regimes generally have a lower number of species. Mowing is the main type of management of both public and private non-woody urban green areas in many cities (e.g. Ignatieva & Ahrné, 2013;Ignatieva et al., 2015), resulting in lawns being a dominant element of the cityscape, with important effects on biodiversity Smith et al., 2014). Because green areas with low and high management regimes represent the majority of the green areas in our study area (Fig. 4), a management relaxation could substantially increase the available habitat for many species. A growing body of evidence has shown mowing reduction as a management prescription to boost certain taxa in a cost-effective way in temperate cities (e.g. Ignatieva & Ahrné, 2013;Ignatieva et al., 2015;Ignatieva & Hedblom, 2018;Smith et al., 2014). Our results extend these findings by adding evidence on 12 distinct taxonomic groups (see Fig. S4) over a whole city. Interestingly, instead of traditionally reported pollinator groups such as bees (Baldock et al., 2019;Salisbury et al., 2015), we find gardens to have a remarkable large predicted richness of functionally diverse groups including wasps, millipedes, snails and ground beetles (see Fig. S4). Hence, as noted by Theodorou et al (2020), there is a need of developing taxon-specific conservation management to promote overall urban diversity. This has shown to have social benefits because of urban green space preferences varying among urban dwellers.

Limitations
Although our study attempts a first urban assessment, the available data timespan only allowed our study to showcase patterns of biodiversity in Zürich for the period 2006-2018. Indeed, the compiled observational dataset comprised species records within these years, and environmental predictors spanned the same range but representing individually specific years within. Provided that Zürich has expanded and some environmental conditions such as microclimate might have changed in parts of the city, using recent predictors and species observations within the same year may have been more relevant. However, macroecological and landscape studies recurrently employ compiled environmental and species datawhich often summarize larger timeline to uncover meaningful ecological patterns (Karger et al., 2017;Phillips et al., 2019;Wüest et al., 2020). Doing otherwise here, would have led to a clear lack of data, missing the large scope of our analysis and associated found ecological patterns. Furthermore, cities undergo rapid changes either by expanding into non-urban areas or by compacting available urban land-covers (e.g. vacant lots, abandoned buildings, see Wolff, Haase, & Haase, 2018). Therefore, minor changes may be expected in already urbanised areas where a large part of our data was sampled.

Using predictive models in urban ecosystems: Opportunities for ecology, planning and society
Predictive models may further provide important insights in urban ecosystems with unique social-ecological features. Spatial information on urban ecosystems are increasingly available related to both their social-ecological properties (Egerer, Fouch, Anderson, & Clarke, 2020;Schell et al., 2020) and biodiversity data (Kendal, Dobbs, & Lohr, 2014;Müller & Kamada, 2011;Ossola et al., 2020), allowing the application of predictive modelling. As shown in this study, predictive models can inform species ecology, management and urban planning. Predictive models can both address traditional questions in urban ecology, but also reveal meaningful spatial diversity patterns and gradients difficult to uncover through other more traditional approaches (Thomas et al., 2018;Zhang, Chen, Xu, Xue, & Ren, 2019). Also, predictive models may help to understand quantitatively and spatially the contribution of different types of predictors on species distribution in urban ecosystems. For instance, several existing methods (Chauvier et al., 2020;Meier et al., 2010) might allow the relative importance of abiotic, biotic and social drivers to be disentangled across the cityscape. Finally, citywide species richness maps could help detect both cold-and hotspots of species in the city, a current gap in cities (Kantsa, Tscheulin, Junker, Petanidou, & Kokkini, 2013; but see also Planillo et al., 2021, andStas et al., 2020). As many cities are developing plans to adapt to climate change that involve increasing green cover (Butt et al., 2018;Seto, Golden, Alberti, & Turner, 2017;Stadt Zürich, 2020b;Yao et al., 2019), outputs of predictive modelling provide an opportunity to also secure urban biodiversity (Butt et al., 2018;Nilon et al., 2017) and improve the accessibility to urban nature and ecosystem services in future urban planning.

Conclusions
While urban ecosystems display features that make them substantially distinct from non-urban ecosystems, they still share main ecological properties. Here, we showed that biodiversity in the city of Zürich follows a similar ecological structure as other non-urban ecosystems with a large proportion of rare and locally occurring species. Furthermore, we found marked biodiversity variation along urban intensity gradients, but even highly urbanized areas can contain rare species. To effectively maintain biodiversity in cities, the unique biophysical and human dimensions of urban ecosystems should be studied in combination. Historically, cities have been treated as non-ecosystems dominated by a handful of synanthropic species, of which many are regarded as pests. Nonetheless, urban ecosystems can have a relatively large amount of biodiversity that shows structured ecological organization. Cities are densifying to palliate urban sprawl particularly in low intensity urban areas. Predictive model outputs may inform urban planning by including biodiversity targets in the existing sustainability and climate adaptation goals, usually more orientated to increase solely human well-being. In conclusion, developing a solid mechanistic understanding of urban ecosystems is key to ensuring that cities remain functional and resilient for both nature and humans in the face of global change.

Funding
This research was supported by the Swiss National Science Foundation (project 31BD30_172467) within the programme ERA-Net Bio-divERsA project "BioVEINS: Connectivity of green and blue infrastructures: living veins for biodiverse and healthy cities" (H2020 BiodivERsA32015104). F.Z. received funding from the Swiss National Science Foundation (project 172198). Y.C. acknowledges the ANR-SNF bilateral project OriginAlps (grant number 310030L_170059).

CRediT statement
J.C.A., M.M. and L.P. conceived the study; Y.C. helped develop the methodology; J.C.A., D.F. and M.M. provided the biodiversity data; F.Z. and C.G. processed and provided the ALS data; J.C.A., Y.C. and P.V. analysed the data; J.C.A. wrote the first version of the manuscript; M.M. and L.P. helped draft the manuscript. All authors commented on the manuscript and gave final approval for publication.
Dawes for the language editing. We thank E. Eggenberg, A. Zannetta, S. Tresch, L. Villarroya-Villalba and K. Kilchhofer for their assistance in collecting biodiversity data. We thank Prof. Dr. E. Parlow for providing the local overwarming data and some of the pollution data. We especially thank S. Fontana, C. Bello, R, Wüest, P. Descombes, C. Waldock, L. Force Seguí and S. Müller for their useful input in data analysis and their support. We would also like to thank two anonymous reviewers for their helpful comments.