Main

Cities represent novel anthropogenic environments, leading to significant mismatches with the evolutionary history of most organisms¹. As urban areas are set to expand two- to sixfold over the 21^st century², urbanization poses a substantial and accelerating threat to global biodiversity³. Yet species vary widely in their responses: some are completely extirpated (‘urban avoiders’), others persist through behavioral or ecological acclimation (‘urban adapters’) and/or rapid adaptive evolution to cities^4,5, and some even thrive in urban environments⁶ (‘urban exploiters’). Quantifying these divergent responses is key for understanding how urbanization will continue to shape global biodiversity and urban ecosystem functioning (e.g., urbanization reshapes a food web^7,8). However, the young field of urban ecology has typically relied on small-scale studies, and/or has often focused on relatively few species constrained to particular taxonomic groups such as birds⁹ or mammals¹⁰. Thus, whether general patterns of urban preference or avoidance emerge across species within and across regional pools remains untested.

Species likely do not respond randomly to urbanization as traits may modulate species’ responses. Trait-based approaches are powerful for identifying generalizable patterns across diverse taxa, particularly at broad spatial and phylogenetic scales¹¹. Although many traits (e.g., diet, reproductive strategy, dispersal ability) may influence species-specific responses to urbanization, body size stands out as both the most widely measured and most taxonomically integrative. It is a key trait of any organism, relating to space use, life-history, and metabolic rate across the tree of life^12,13,14, making it a plausible direct or indirect predictor of species’ responses to key stressor gradients in urban landscapes, including resource fragmentation, heat stress from urban heat islands (UHI), and anthropogenic disturbance (Fig. 1).

Fig. 1:

Body size correlates with urban mobility, but in complex ways¹⁵. For instance, while higher mobility allows some animal species to better cope with the fragmented nature of urban habitat resources, reduced mobility might be advantageous for exploiting localized resources and avoiding risks associated with the urban matrix^16,17. In plants, body size—often indexed by height—correlates with competitive ability, light acquisition, and reproductive strategy, with larger species generally having greater resource needs but potentially higher resilience to urban stressors¹⁸. Notably, the UHI may drive shifts to smaller-bodied animal species due to elevated metabolic costs^14,15, echoing similar patterns observed under global warming^19,20. However, these shifts toward smaller size can be overruled by requirements for increased urban mobility in taxa where mobility increases with body size¹⁵. All in all, despite body size’s fundamental role in modulating species responses to their environment, a systematic, cross-taxon assessment of how urbanization filters body size distributions is lacking.

Here, we present a global synthesis of urban tolerance across more than 30,000 animal and plant species, using occurrence records from the Global Biodiversity Information Facility (GBIF) and remotely sensed night-time light intensity as a proxy for urbanization. We define a continuous metric of urban tolerance as a species’ relative affinity to urban areas, with negative values indicating avoidance and positive values indicating preference²¹. To account for regional variation in both species pools and in species’ urban tolerance across their range^22,23, we calculated urban tolerance within 52 geographic subrealms (Fig. S1). Our dataset spans 47 taxonomic classes—from flowering plants (8,980 species) to insects (8,674 species) to mammals (648 species)—and represents a uniquely broad assessment across the tree of life. We first introduce the concept of “Species Urbanness Distribution” (SUD) to characterize the composition of a regional community in terms of its species’ urban affinities. We then use Bayesian hierarchical modeling to test whether body size predicts urban tolerance across 17,722 species from 371 families, across multiple subrealms. These analyses allow us to quantify both consistency and heterogeneity of body size effects on urban tolerance across multiple taxonomic levels, and across regions. This synthesis addresses two core questions: (1) What is the shape of species’ urban tolerance distributions within regional communities? (2) Does body size consistently correlate with species’ urban tolerance across taxonomic groups and biogeographic contexts?

Results and Discussion

Species Urbanness Distributions (SUDs)

We identified a consistent pattern in the extent to which regional communities—defined as all species within a given biogeographic subrealm—consist of urban avoiders and exploiters. We term this the Species Urbanness Distribution (SUD), which captures the full distribution of urban tolerance values across species in a community. SUDs exhibited a characteristic shape: a pronounced peak at negative values and a long right-skewed tail toward high positive values, indicating that most species are ‘urban avoiders’, while only a few are ‘urban exploiters’. This pattern was consistent across subrealms and several taxonomic levels (Fig. 2). Conceptually, SUDs parallel the well-known Species Abundance Distributions²⁴ (SADs), which describe the general law of communities being composed of many rare and few common species^24,25.

Fig. 2:

Similarly, much like the skewed distributions observed in SADs^24,26, the skewed shape of SUDs indicates that only a select few species seem to possess traits that enable them to thrive in urban environments. While previous studies have reported a similar pattern at smaller spatial scales²⁷, our synthesis demonstrates this applies across the globe.

The skewed shape of SUDs suggests that traits enabling species to tolerate urban environments are both rare and specific—an indication that adaptation (sensu lato) to urban settings is uncommon, yet occurs geographically widespread across taxa. These traits may be partly found in human-commensal species, which are found globally and often include invasive non-native species²⁸. Such traits “pre-adapted” to urban conditions allow for some species to not only persist but thrive in urban environments where most species cannot. Indeed, the median urban tolerance is most often below one, indicating widespread avoidance among species. This is expected, given that urbanization typically involves severe loss and fragmentation of habitat resources for most organisms¹, which negatively impacts numerous taxa, including birds²⁹, beetles³⁰, and moths³¹.

Additionally, the UHI—typical of most cities³²—imposes thermal stress, especially during extreme heat events³³, which may drive community shifts toward heat-tolerant species, as observed in ants³⁴, bees³⁵, and plants³⁶. Yet, urban environments also frequently contain green infrastructure and remnants of natural habitats that can act as sanctuaries for rare and endangered species³⁷, though these spaces are often also dominated by widespread generalist species. The multitude of environmental stressors, combined with competition dynamics between urban-tolerant and intolerant species, likely influence the shape of SUDs (Fig. 2). Future ecological research should aim to disentangle these contrasting drivers across taxa and regions to better understand the ecological filters at play.

The relationship between body size and urban tolerance

To test whether there is a consistent pattern of body size filtering associated with urban tolerance across the tree of life, we integrated our species-level urban tolerance values with body size measures collated from the literature. We performed this analysis for 17,722 species spanning both plants (i.e., plant height) and animals (e.g., body length, mass, wingspan). We first evaluated the overall relationship between body size and urban tolerance across all taxa, fitting two models (see Methods for details). We found that the effect of body size on urban tolerance is stronger in plants than in animals (0.64 vs. 0.21), although uncertainty was substantial in both groups, with zero-overlapping 95% credible intervals (plants: -0.68 to 1.84; animals: -0.25 to 0.66). Variation in urban tolerance was most pronounced at the family level (sd = 2.43), followed by class and order level (sd = 2.34 and 1.75, respectively). Similarly, the effect of body size varied most among families (sd = 1.24), compared to orders (sd = 0.73) and classes (sd = 0.27), suggesting that body size predicts urban tolerance most consistently within the finer family-level taxonomic scale, rather than at higher scales where other factors may obscure relationships. Based on these patterns, we focused subsequent analyses at the family level, fitting Bayesian models for families with data on at least ten species to independently explore the relationships between urban tolerance and body size.

Across 371 families (93 plant and 278 animal families), we found no universal relationship between body size and urban tolerance; instead, patterns were heterogeneous (Fig. 3). Using 80%, 90%, and 95% credible intervals to infer weak, moderate, and strong effects, respectively, we found evidence of body size filtering in 83 families (29 strong, 19 moderate, and 35 weak). These included 48 animal families (17% of animal families) and 35 plant families (38% of plant families), indicating that for most families, body size does not consistently predict urban tolerance. This suggests that other traits or ecological factors likely underpin the emergence of urban exploiters and the skewed shape of SUDs. Nonetheless, among plant families, the effect of body size on urban tolerance was predominantly (N=77; 83%) positive, with 23 families showing strong or moderate effects. In contrast, animal families (N=278) exhibited more variable responses: 58% showed positive effects and 42% negative, with only a subset showing strong or moderate effects (positive: 18 families; negative: 6 families). Notable families with the strongest positive associations included Salicaceae (willows) and Columbidae (doves and pigeons), where larger species are more urban-tolerant. In contrast, Dipsadidae (snakes) and Accipitridae (birds of prey) showed the strongest negative associations. Further discussion of notable order-level effects is provided in the Supporting Information. While the overall patterns for family-level effects are visible in Fig. 3, we also developed an interactive online visualization of family-level effect sizes, available here.

Fig. 3:

The variability in both the direction and strength of urban body size filtering across taxa likely reflects the diversity of urban stressors and the contrasting taxon-specific pathways through which body size mediates sensitivity to them. In raptors, for example, the negative association between body size and urban tolerance—consistent with previous findings^38,39—likely stems from large species’ need for extensive hunting territories, whose availability and access are constrained by urbanization-driven habitat loss and fragmentation. Similarly, the strong negative effect observed in Coleoptera supports the hypothesis that the UHI drives shifts to smaller body sizes in urban-tolerant ectotherms, consistent with Atkinson’s temperature-size rule⁴⁰. Larger beetles may also be more vulnerable to urban disturbance due to lower reproductive output and reduced dispersal capacity⁴¹. In contrast, urban communities of Lepidoptera (moths and butterflies) typically shift toward larger—and hence more mobile—species than rural communities¹⁵, likely because greater mobility facilitates persistence in typically fragmented urban landscapes⁴².

Order-level analyses (see Supplementary Text) offer further insight into the mechanisms shaping body size responses to urbanization. In plants, excluding grasses and ferns, we consistently observed a positive association between body size and urban tolerance. Taller species may gain a competitive edge in urban settings⁴³, by exploiting higher temperatures, increased CO₂ and nitrogen availability⁴⁴, and typically reduced O₃ concentrations⁴⁵. In contrast, grasses and ferns favored smaller species, likely stemming from high disturbance regimes in urban grasslands and ruderal environments. Among ectotherm animals—including Diptera, Coleoptera, Hemiptera, and Squamata—the general expectation is a shift toward smaller species, consistent with metabolic constraints imposed by the UHI (Fig. 1, Fig. 3). However, exceptions occur when urban fragmentation requires greater mobility, favoring larger and more mobile species¹⁵. Urban spiders also exhibit body size reductions, driven by both elevated temperatures and reduced prey size⁴⁶. For endotherms such as rodents and birds, we typically anticipate a shift toward larger species because of increased habitat fragmentation, especially when body size correlates positively with mobility. While this pattern held for most bird families, notable exceptions—like Piciformes and Accipitridae—demonstrate negative associations, possibly due to constraints related to large home range requirements and prey size^47,48. Future research should further explore such trophic consequences of urban size filtering, as these effects can either align across trophic levels or become mismatched. For instance, urban environments may favor larger predators—such as insectivorous birds and mammals—while simultaneously selecting for smaller insect prey species. Such mismatches can lead to nutritional stress, reduced reproductive success, and ultimately a more severe homogenization of urban predator guilds^49,50. Another ecological consequence of size mismatches may occur when urban pollinator assemblages become disproportionally large or small relative to the floral traits of native plants, potentially reducing pollination efficiency⁵¹—although rapid evolutionary shifts in floral morphology could help mitigate these effects⁵².

Current limitations and future directions

Our analysis, incorporating over 17,000 species, represents the most extensive assessment to date of how body size relates to urban tolerance across the tree of life. Such scale was enabled by the growing mobilization of biodiversity data through GBIF, particularly from citizen science platforms like iNaturalist⁵³. Still, our analysis largely focused on common species, as we applied a cutoff of at least 100 observations per species per subrealm. As a result, rare or less frequently observed species are underrepresented, and their responses to urbanization remain an important avenue for future study. This is especially important given the growing recognition of cities as sanctuaries for various species³⁷. Nevertheless, our findings highlight the relevance of common species to understanding macro-ecological patterns (sensu⁵⁴). We focused here on broad-scale patterns, but future work should further explore the context-dependency of body size filtering— for instance, how it varies with city size, urban green infrastructure configuration, regional climate, and human cultural preferences. For example, one possible explanation for the observed trend of larger, urban-tolerant plants may lie in the widespread use of non-native woody ornamentals, such as alien trees and tall shrubs, which reflect human aesthetic preferences^55,56— suggesting human preferences are an additional filter on plant size. In parallel, while our study examined interspecific variation, growing evidence also points to important intraspecific responses, where body size shifts within species can modulate urban tolerance^31,57,58.

We anticipate that the urban tolerance scores presented here will form a valuable foundation for future research and local-scale planning efforts, in particular those that use citizen science to track restoration progress. For example, these tolerance scores could inform urban greenspace surveys aimed at calculating integrity indices—providing a repeatable, quantitative tool to assess and monitor the long-term success of urban ecological restoration by measuring the ‘urbanness’ of local species communities in urban greenspaces⁵⁹. In this context, incorporating taxon-specific body size-urban affinity relationships can enhance conservation outcomes by tailoring strategies to the size dynamics of assemblages. For instance, in assemblages where larger species are disproportionally filtered out (i.e., a negative size-tolerance link), efforts should focus on increasing the size and quality of habitat patches and mitigating the UHI. Conversely, in assemblages where smaller, less-mobile species are more vulnerable, improving functional connectivity between habitat patches should be prioritized⁶⁰. For butterflies in particular, Pla-Narbona et al.⁶⁰ advocate for such context-specific management strategies—balancing habitat patch quality and connectivity—to promote more diverse urban butterfly communities. Future research should also test the role of other key functional traits, like seed size (Fig. S2), dispersal capacity (Fig. S3), or thermoregulatory ability⁶¹, as analyses involving multiple functional traits will likely offer a better understanding of the mechanisms underlying species’ urban tolerance responses.

Conclusions

Our global analysis of over 17,000 common species provides new insight into how diverse taxa respond to urbanization, revealing a consistent skew in urban tolerance—characterized by many urban avoiders and few exploiters—across biogeographic regions and taxonomic groups. We introduce the concept of Species Urbanness Distributions (SUDs) as a novel framework to quantify and compare the impact of urbanization on species communities. Much like Species Abundance Distributions (SADs) are essential for ecology and biodiversity research, SUDs offer a generalizable lens through which the community structure and ecological resilience of urban biotas can be assessed. Understanding the processes that generate SUDs, and understanding the ecological impacts of differently shaped SUDs, represent key avenues for future research. Although body size emerged as a predictor of urban tolerance, we found not only substantial heterogeneity across families and orders, but also that body size filtering alone is unlikely to explain the consistently skewed SUD shape. Nevertheless, trait-based approaches—especially those integrating multiple functional traits—hold strong potential for uncovering the processes driving the diversity of species’ urban responses and for interpreting the shape and skew of SUDs. Moreover, trait-based predictions of species’ vulnerability could be used to formulate effective strategies to promote biodiversity in urban environments⁶². Looking ahead, the continued growth of citizen science data platforms, such as iNaturalist, will play an increasingly critical role in tracking and understanding the mechanisms driving biodiversity responses to urbanization. These data streams will enable broader quantification of species’ urban tolerance, including rare and currently underrepresented taxa. Combined with trait-based modeling, this growing dataset could be leveraged to identify urban-vulnerable species and hence guide urban habitat restoration initiatives. In this way, SUDs and trait-informed predictions offer powerful tools for more effective urban strategies to conserve biodiversity on an increasingly urban planet.

Methods

Our methodological approach can be broken down into three key steps: (1) quantifying species-specific urban tolerance scores, stratified by subrealm (i.e., biogeographic region); (2) collating measures of body size for as many species as possible for which we were able to quantify urban tolerance; and (3) quantifying the relationship between urban tolerance and body size at various taxonomic levels.

Quantifying urban tolerance

Our goal was to derive a continuous metric of the extent to which species vary in their responses to urban areas (hereafter urban tolerance), ranking species from those that are largely intolerant (i.e., urban avoiders) to those that thrive in urban areas (i.e., urban exploiters). To quantify species-specific urban tolerances, we followed the approach described in detail in Callaghan et al.^21,63,64.

Following Callaghan et al.⁶³, we downloaded data from the Global Biodiversity Information Facility (hereafter GBIF) (https://doi.org/10.15468/dl.4dcbgt) on February 4th, 2021, including ∼1.4 billion biodiversity records from >24,000 datasets. GBIF data were filtered to only include observations that were recorded as species, and an additional step was taken to ensure that the listed genus in GBIF matched the genus of the species name, ensuring validity of the taxonomic nomenclature of GBIF. Due to uncertainty in matching observations with remotely-sensed products, any GBIF observation with a coordinate uncertainty > 1 km was removed. We only included species from a taxonomic Class (i.e., the rank between Phylum and Order) that had at least 10 species reported to GBIF, focusing on the most common classes for downstream analyses. We then overlaid GBIF species occurrence observations with a remotely-sensed layer representing a continuous proxy of urbanization—Visible Infrared Imaging Radiometer Suite (VIIRS) night-time lights⁶⁵. We used this proxy as it is a continuous measure of urbanization commonly used to represent urban extent^66,67,68. Previous work has shown that VIIRS night-time lights is negatively correlated with greenness measured through the Enhanced Vegetation Index (EVI) and positively correlated with population density^63,65. Further, previous work has also shown that measures of urban tolerance calculated from VIIRS night-time lights, population density, and global human modification index are all strongly correlated⁵⁹. We used Google Earth Engine⁶⁹ for our geospatial processing and data were obtained for the VIIRS Stray Light Corrected Nighttime Day/Night Band Composites product, representing monthly composites, (i.e., this dataset in Google Earth Engine: NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG).

We took the median of all monthly composites for each pixel to calculate a pixel-level urbanization value, measured in average radiance, and used imagery from January 2015 to January 2021 to calculate this median. GBIF data were filtered from January 1st, 2010 to February 4th, 2021. We acknowledge that these two temporal scales do not precisely match, but because urban conversion at the landscape/regional scale happens relatively slowly, we assume that the level of urbanization between 2015 and 2021 accurately captures the relative differences between regions with low and high urban land cover.

To avoid computational overload, we used geohash encoding to assign every record in GBIF to a “pixel” representing the VIIRS average radiance. We used geohash7 encoding to divide the geographic area into grid cells (referred to as blocks), each representing an approximate spatial area of 150 m². Each GBIF record was assigned to a corresponding geohash7 block based on its geographic coordinates. The VIIRS night-time lights data, with a native resolution of ∼500 m², was then matched to these blocks by assigning each geohash7 block the average VIIRS radiance value that intersects it. While the geohash resolution is finer than the VIIRS data, this approach allows us to efficiently associate highly precise GBIF occurrence points with nearby VIIRS pixels using spatial indexing. This avoids unnecessary redundancy in the data while maintaining the original VIIRS resolution.

After each record from GBIF was assigned a measure of VIIRS, described above, for each species, we took the mean of the distribution of VIIRS values as a relative species-specific measure of a response to urban land cover (hereafter referred to as urban score). This urban score represents the relative ranking of that species to other species in response to urban land cover. To account for potential intraspecific differences in species’ response to urbanization throughout their global range, we stratified our analysis to a total of 52 subrealms⁷⁰ (Fig. S1). Therefore, a species could be represented across multiple subrealms with differing measures of urban tolerance (Fig. S4). Importantly, this continuous metric of urban tolerance is a relative measure of a species’ preference, or affinity, to urban areas: it should be interpreted only within each subrealm (i.e., urban tolerance scores are not interpretable across subrealms), providing a comparative or relative ranking of how a species responds to urbanization compared with other species within that subrealm. Within each subrealm, we further accounted for the potential of different levels of urbanization by scaling each species’ urban score by subtracting the mean VIIRS of all observations in the subrealm (this value is hereafter referred to as urban tolerance). This ‘urban tolerance’ (Fig. S5) value can be negative—when species under-occupy urban areas suggesting they actively avoid them—or positive—when species over-occupy urban areas suggesting they prefer them (i.e., ranging from urban avoiders to urban exploiters, respectively).

We applied the above procedure to any species that had at least 100 observations in a subrealm. Previous work has shown that this cutoff approximates the variability of a species’ response to urbanization and ensures that enough of a species urban habitat use has been captured^21,71. After this filtering we were left with a total of 56,181 data points—i.e., unique urban tolerance measures of a species in a subrealm—for potential inclusion in our analysis. These records represented 30,373 species that had at least one urban tolerance measure, from 47 classes and across 51 subrealms (Fig. S6; Data S1). The number of subrealms a species was found in was mostly (64%) 1 but ranged from 1 to 34 (Fig. S7). All urban tolerance measures (i.e., urban tolerance by subrealm measures) which we produced as part of our workflow are made available, as we hope to advance further testing of urban tolerance among different taxonomic groups and regions. But correct interpretation of these urban tolerance values is essential before use.

Sampling biases are common in citizen science data. In GBIF, urban areas are better sampled compared to remote regions, as contributors typically concentrate around areas with high human density. However, sampling bias across the urban modification gradient was accounted for because each species’ urban tolerance value was calculated relative to the observations across its geographic range (i.e., within the subrealm), compared with other species within the same geographic range. Importantly, while disproportionate sampling in human occupied areas may inflate the urban tolerance value, each species in the dataset was subject to similar biases, enabling comparisons across species. And hence, the focus here should be on comparisons among species. Previous work has found that this approach of assigning urban tolerance is strongly correlated with occupancy-detection models²¹. As illustrative examples, these data can be used to better understand the differences in urban tolerance among Hymenoptera in Northeast American Forests (Fig. S8), Lepidoptera in Southeast Asian Forests (Fig. S9), or Asterales in Scandinavia & West Boreal Forests (Fig. S10).

Species urbanness distributions

Our analysis of species urbanness distributions was descriptive in nature, relying on data visualization using ggplot2 to make histograms and density histograms of urban tolerance values for each subrealm. In these analyses at the subrealm level, species are treated equally without any differentiation taxonomically, compared with the taxonomically-stratified analyses where species within each taxonomic rank are treated as similar.

Aggregating measures of body size

Our objective was to compile a dataset of species-specific measures of ‘body size’, or ‘organism size’ more broadly to integrate with our measures of species-specific urban tolerance described above. We aimed to find measures for as many species as possible, maximizing our sample size for modeling. For plants, we used ‘plant height’ as a proxy for body size, and data were downloaded from the TRY database—a database compilation of plant traits from many different datasets⁷². For animals, body size can be measured in different ways, depending on the taxonomic group of interest (e.g., body length, wingspan, biomass, shell size, radius, forewing length, body mass, or biovolume). Therefore, our dataset was purposefully broad in its composition of body size estimates. Because ‘body size’ is the predominant form referred to in the literature, we use this throughout, but take it to broadly mean ‘organism size’ as discussed above.

To develop this dataset, we first performed a semi-systematic broad literature search, using Web of Science, Google Scholar, Zenodo, Dryad, and Figshare with variations of the key words: “body size”, “organism size”, “body length”, “body mass” and “interspecific”, “global patterns”, “animals”, and higher-level groups such as classes (e.g., “Aves”, “Birds”, “Mammals”, “Mammalia”). We started by downloading known compilations for various taxa, for example mammals⁷³ or Odonata⁷⁴. In these instances, we spent minimal further time searching for specific body size data for these larger taxonomic groups, assuming that these compilations had already aggregated the majority of data available. Following this, we conducted more detailed literature searches, predominantly using Google Scholar and Web of Science, which focused on our a priori list of species for which we had a measure of urban tolerance (Data S1). These were conducted at class, order, or family taxonomic levels, but not at taxonomic levels below family. Multiple keywords were used for each group while searching when it was obvious to do so. For example, when searching for body size measures in Coccinellidae both “ladybirds”, “ladybugs”, and “Coccinellidae” were used in search strings, as well as “wing length”, “body length”, or “body size”. Search effort was qualitatively proportional to the number of species we had in each taxonomic level. For example, because we only had 43 species with urban scores from the class Diplopoda but 3,515 species from the order Lepidoptera, we spent proportionately more time and effort searching for body size data for Lepidoptera. Similarly, we spent more time on hard-to-find taxa, namely groups from Insecta. In our searches we generally ignored studies that were focused on intraspecific variation in body size for one or a few species, unless we found that the species was included in our urban tolerance list. In addition to literature searches, we emailed networks of colleagues, looked at reference lists in papers that had data, and emailed corresponding authors of papers who had potentially relevant data but for which data were not accessible.

After this semi-systematic searching, we re-assessed the remaining species for which we had an urban tolerance measure but not a measure of body size. We then performed a second-level targeted search (i.e., more systematic) for any family with at least 10 species remaining for which we had a measure of urban tolerance. First, we searched Web of Science and Google Scholar one family at a time and looked at the first 100 hits, with the family name ‘and’ “body size” used as search terms. Second, we used Google search engine to search for these remaining species individually (N∼3,100). We focused on including species for which Google immediately returned an estimate of body size (e.g., through Wikipedia, or an aggregating website), and did not go to original species descriptions. This detailed searching, and the aggregation of body size measures, concluded in June 2022.

Measurements of body size

In our searching for body size, we were agnostic to the type of body size, given that many different measures are used for quantifying body size, or proxies for body size. For each dataset we downloaded or processed, we kept track of the following key metadata properties: ‘type’, ‘units’, and ‘measurement detail’. Type refers to the different types of body size measurement, ranging from body length and body mass to Weber’s length in ants. Table S1 shows a list of the different body size types encompassed in our analysis (see statistical details below for how different body size types were accounted for in our analyses). For each dataset, we also kept track of the units (e.g., mm, g, mg cm^-2, or mm^3). And lastly, wherever possible, we noted the measurement detail of the body size measures, whether the mean, median, or maximum was used, or whether it was for males or females only. Many times, however, this information was not readily available, and thus this metadata field was not required in our dataset composition.

Taxonomic harmonization and dataset integration

Because body size data were aggregated across many potentially different datasets, we performed a taxonomic harmonization step to regain some species that did not match due to typos, spelling mistakes, or synonyms. For this, we treated Plantae and Animalia separately and used the taxize package⁷⁵ in R to search through the list of species and find any synonyms or fuzzy matches of other taxonomic entities, always matching back to the GBIF taxonomic backbone, given that our urban tolerance measures follow this taxonomy (Data S1). With the possible synonyms and matches, we then integrated the urban tolerance measures with the body size dataset using either direct matches or matches of synonyms. The taxonomic harmonization step regained a total of 831 (see Fig. S11) species that were integrated in the dataset for analysis. Data S3 provides details on the 223 different ‘datasets’ (i.e., downloaded data, or manually added data) that were processed as measures of body size. The most prevalent class in our final dataset (Fig. S12) was Insecta (N=7,526), followed by Magnoliopsida (N=4,872) and Aves (N=3,313).

Our final dataset for potential analysis and modeling (Data S2) included a total of 94,087 observations (i.e. unique combination of species urban tolerance values, subrealm, and body size measure) for 20,957 species that had at least one measure of urban tolerance and at least one measure of body size.

Quantifying the relationship between urban tolerance and body size

Our objective here was to understand the relationship between urban tolerance (the response variable) and body size (the predictor variable). Our dataset had potentially multiple measures of body size for an individual species, and differing levels of metadata. For example, some datasets may have specified whether the mean was used, or others might have specified for females only (see Table S1). Therefore, before modeling, our goal was to minimize the undue effect of having too many random effect levels. First, because an individual species could have multiple measures of body size, we only selected the body size measure that was most common among the species in each family. For example, if a species had both a measure of body length and of body mass, we selected the measurement type that was most commonly available among the species in each family and the largest source of data. Second, we further reduced the potential number of levels and variation in body size measurement by creating a synthetic “metadata” variable that was the combination of ‘type’, ‘units’, and ‘measurement detail’. By doing this, we decreased the potential number of levels of a random effect from 223 (i.e., the number of data sources) to 71 (i.e., the combination of type, unit, and measurement detail). Third, several papers/data sources could have all reported the same metadata (e.g., maximum body length for females in mm), and instead of being treated as individual random effects because they were from different data sources, we assumed that these represented similar data and collapsed them into one random effect level (here, for maximum body length for females in mm).

Overall kingdom-level model construction

We used Bayesian hierarchical models to test for the effect of body size on urban tolerance. In our models, urban tolerance was the response variable and body size was the predictor variable of interest. To evaluate the overall relationship between body size and urban tolerance, we fitted two comprehensive models that incorporated hierarchical taxonomic levels and the relationship with body size. The first model included fixed effects for kingdom and log-transformed body size due to its skewed distribution, along with their interaction, to investigate differences in the body size effect on urban tolerance between plants and animals. The random effects accounted for hierarchical taxonomic levels, incorporating random intercepts and slopes for class, order, and family to capture variation at each level. Additional random effects for metadata and subrealm were included as random intercepts to control for study-specific differences in data and regional differences in urban tolerance. The model formula was: urban tolerance ∼ body_size_scaled_log10 * kingdom + (1 + body_size_scaled_log10 | class/order/family) + (1 + body_size_scaled_log10 | metadata) + (1 | subrealm). The second model built upon this structure but only included an interaction term between log-transformed body size and kingdom as well as kingdom as a fixed effect to determine if the effect of body size on urban tolerance varied consistently across the plant and animal kingdom. It maintained the same hierarchical random effects structure as the first model, with class, order, and family nested to account for hierarchical variation, and metadata and subrealm included as additional grouping factors. The formula was: urban tolerance ∼ kingdom + body_size_scaled_log10:kingdom + (1 + body_size_scaled_log10 | class/order/family) + (1 + body_size_scaled_log10 | metadata) + (1 | subrealm).

Family-specific model construction

However, for our main analyses we chose to employ a ‘many models’ approach, where we fit individual models for each family that met a set of minimum criteria because many families do not have a known association or published information on the extent to which urban tolerance associates with body size. A family was only included if there were at least 10 unique species with an urban tolerance measure and body size measure, and a subrealm was only included if it had at least two species included to specify a random slope. For example, if there was only one species of a particular family included in our dataset, this family would not have received an independent model. We decided to focus our main analysis at the family level for three key reasons. First, this was done because for most of the taxonomic families included in our analyses there is little a priori expectation about the effect of body size on urban tolerance (i.e., little to no literature on the subject), and therefore we wanted each model to represent an advance in knowledge of that relationship for that family. Second, because there were large imbalances in the number of data points available for each taxonomic group, for instance some families having hundreds of species with available data and others with only a handful of species, we wanted to avoid shrinkage where the model results were primarily driven by the taxonomic groups with a larger number of data points represented. Third, we wanted to avoid an instance where a family (or other taxonomic grouping) could have little or no relationship between body size and urban tolerance, but if included together in one model could result in a positive effect (see Fig. S13 for a simulated example). This was a smaller concern at the family level, which we present as the main analysis, because of the smaller imbalance in data availability across the taxa within each family.

When we fit the family-specific models, there were different requirements for the random effect structure depending on how many subrealms and metadata sources (i.e., types of body size measures). In the simplest case, all species in a family were just from one subrealm and had just one data source for body size. But this was not always the case, and therefore we developed four different ‘types of models’ which varied in their treatment of random effect structure. Each ‘type’ of model, however, had the same overall objective—to quantify the relationship between body size and urban tolerance. The four model types are as follows. Type 1: had at least two subrealms and at least two types of synthetic metadata (see above) sources included. Type 2: had data from only one subrealm and only one type of synthetic metadata source included (but could have data from more than one initial data source). Type 3: had at least two subrealms and only one type of synthetic metadata source included. Type 4: had data from only one subrealm and at least two types of synthetic metadata sources included. For all types of model fits, the response variable was the urban tolerance measure, described above, and the fixed effect of interest was body size which was log10-transformed and scaled and centered by the synthetic metadata source. Type 1 had a random intercept and slope for subrealm, allowing the intercept and slope to vary across different levels of subrealm as well as a random intercept for the synthetic metadata source. Type 2 had no random effects. Type 3 had only a random intercept and slope for subrealm, allowing the intercept and slope to vary across different levels of subrealm. Type 4 had only a random intercept for the synthetic metadata source. Modeling was done using the ‘brm()’ function in R from the brms package^76,77,78. For all model types we specified a standard normal prior, with a mean of 0 and a standard deviation of 1 applied to the fixed effects of the model. Models were fit with 1000 warmup iterations, 6000 iterations for the Markov Chain Monte Carlo sampling, 4 chains, and the adaptation target acceptance probability was set to 0.99.

For models of type 1 and type 4, where we combined metadata for body size measurements among different species, we also explored more complex model structures that included a random slope for body size by metadata source, in addition to subrealm. Such a model structure would allow the relationship between body size and urban tolerance to vary across measurement types. However, these models often substantially reduced interpretability and inflated uncertainty around slope estimates, especially when metadata groups were sparsely represented, which was common. Additionally, many of the families only had 2 types of metadata, limiting our statistical power to fit random slopes. Comparative analyses showed that such models produced weaker or more uncertain estimates than our main model structure (i.e., random intercept for metadata of body size). We also note that to account for heterogeneity in body size measurement types, all body size values were log10-transformed, scaled, and centered by their metadata source prior to analysis, following standard practices in allometric studies. We use two illustrative examples (Figures S14 and S15) to illustrate the influence of including a random slope for metadata source.

Additional analyses beyond family level

Although our main analysis was focused on the family level, one potential drawback is that species could be excluded from any of the analyses for a family that has a small number of species. Plus, because some analyses have previously been conducted at higher taxonomic levels (e.g., Lepidoptera), we repeated the above procedure described for family at the order (see Supplementary Text) and kingdom level.

Statement of authorship

CTC, DEB, TM conceptualized the project and developed project methodology. CTC, TM, VS, JHW, BMM conducted analyses. CTC conducted project administration and supervision. CTC, DEB, TM wrote the original draft, and CTC, DEB, VS, BMM, LA, IS, JHW, TM were involved in review and editing.

Supplementary Text

In addition to family-level effect sizes, we performed the analysis at the order level, and found similar heterogeneous results. Among the 128 orders, a total of 45 had some evidence of an effect of body size on urban tolerance (14 weak, 6 moderate, and 25 strong). Some of the orders with the strongest positive effect size included Sapindales (dicots, soapberry plants) and Columbiformes (birds, doves and pigeons), and orders with the strongest negative effect sizes included Coleoptera (insects, beetles) and Piciformes (birds, woodpeckers). Of note are the differential responses down the taxonomic guilds, where some families can show contrasting results, but there may be a strong effect when looking at the order level, or higher (e.g., Fig. S2).

Fig. S1:

Fig. S2:

Fig. S3:

Fig. S4:

Fig. S5:

Fig. S6:

Fig. S7:

Fig. S8:

Fig. S9:

Fig. S10:

Fig. S11:

Fig. S12:

Fig. S13:

Fig. S14:

Fig. S15:

Table S1:

Data accessibility

All data analysis was conducted in R statistical software and relied heavily on tidyverse (Wickham et al. 2019). The doi representing our download is: https://doi.org/10.15468/dl.4dcbgt. Our processing of the GBIF data to quantify urban scores is available here: https://github.com/coreytcallaghan/body_size_vs_urban_tolerance. Not all body size data can be openly shared due to copyright restrictions, however the freely available dataset, cleaned and processed for our analysis is available here: https://github.com/coreytcallaghan/body_size_vs_urban_tolerance. In this repository we did provide reproducible data and code, albeit not for the entire dataset of body sizes. We will archive a cleaned version of this repository in Zenodo upon acceptance of this manuscript.

Acknowledgements

We thank Tim Robertson and the GBIF support team in supplying avro files of the GBIF download (https://doi.org/10.15468/dl.9ufars) for upload to our personal Google Big Query. Henrique M. Pereira provided helpful discussion in the early stages of this project. We thank the following people for help in qualitatively assessing our proxy for urban tolerance: Luis Aguirre, Timothy Bonebrake, Jason Bried, Alex Córdoba-Aguilar, Michael Crossley, Wouter Deconinck, Geert De Knijf, Tim Gardiner, Michael Greenfield, Ralf Gyselings, Frederik Hendrickx, Axel Hochkirch, Roy Kleukers, Anton Kristin, Paulo Lemos, André Victor Lucci Freitas, Onildo João Marini Filho, Daniele Matenaar, Elisa Montes-Fuentemayor, Tara Murray, Howon Rhee, Leslie Ries, Szabolcs Sáfián, Ann Swengel, Scott Swengel, Pieter Vantieghem, Wade Worthen. We thank the following people for supplying unpublished data of body sizes: Roel van Klink, Michelle Tseng, Erik Öckinger, Philip Barton, Werner Ulrich, Axel Höchkirch, Thomas Sherratt, Michiel Wallis de Vries, Lea Heidrich, Jörg Müller, Robert Tropek, Maldwyn Evans, Martin Gossner, Bernhard Hausdorf, Stano Pekar, Carlo Seifert, Onildo João Marini Filho, Sei-Woong Choi, Dominik Rabl, Markus Franzén, Elena Piano, Derek Hennen. CTC, DEB, IS acknowledge funding of iDiv via the German Research Foundation grant DFG FZT 118. CTC was supported in part by a Marie Skłodowska-Curie Individual Fellowship grant 891052 and was supported in part by the intramural research program of the U.S. Department of Agriculture, National Institute of Food and Agriculture, Hatch FLA-FTL-006297. VS was supported by Georgetown University, a U.S. National Science Foundation Graduate Research Fellowship grant 1937959, and a David H. Smith Conservation Research Fellowship. LHA was supported by the Academy of Finland grant 340280.