Body-size shifts in aquatic and terrestrial urban communities

Body size is intrinsically linked to metabolic rate and life-history traits, and is a crucial determinant of food webs and community dynamics1,2. The increased temperatures associated with the urban-heat-island effect result in increased metabolic costs and are expected to drive shifts to smaller body sizes3. Urban environments are, however, also characterized by substantial habitat fragmentation4, which favours mobile species. Here, using a replicated, spatially nested sampling design across ten animal taxonomic groups, we show that urban communities generally consist of smaller species. In addition, although we show urban warming for three habitat types and associated reduced community-weighted mean body sizes for four taxa, three taxa display a shift to larger species along the urbanization gradients. Our results show that the general trend towards smaller-sized species is overruled by filtering for larger species when there is positive covariation between size and dispersal, a process that can mitigate the low connectivity of ecological resources in urban settings5. We thus demonstrate that the urban-heat-island effect and urban habitat fragmentation are associated with contrasting community-level shifts in body size that critically depend on the association between body size and dispersal. Because body size determines the structure and dynamics of ecological networks1, such shifts may affect urban ecosystem function. The urban-heat-island effect drives community-level shifts towards smaller body sizes; however, habitat fragmentation caused by urbanization favours larger body sizes in species with positive size–dispersal links.

Publisher: NPG; Journal: Nature: Nature; Article Type: Biology letter ms: 2017-07-09177B Body size is intrinsically linked to metabolic rate and life-history traits, and is a crucial determinant of food webs and community dynamics 1,2 . The increased temperatures associated with the urban-heat-island effect result in increased metabolic costs and are expected to drive shifts to smaller body sizes 3 . Urban environments are, however, also characterized by substantial habitat fragmentation 4 , which favours mobile species. Here, using a replicated, spatially nested sampling design across ten animal taxonomic groups, we show that urban communities generally consist of smaller species. In addition, although we show urban warming for three habitat types and associated reduced community-weighted mean body sizes for four taxa, three taxa display a shift to larger species along the urbanization gradients. Our results show that the general trend towards smaller-sized species is overruled by filtering for larger species when there is positive covariation between size and dispersal, a process that can mitigate the low connectivity of ecological resources in urban settings 5 . We thus demonstrate that the urban-heat-island effect and urban habitat fragmentation are associated with contrasting community-level shifts in body size that critically depend on the association between body size and dispersal. Because body size determines the structure and dynamics of ecological networks 1 , such shifts may affect urban ecosystem function.
Body size is a fundamental species trait relating to space use and key life-history features such as longevity and fecundity 6 . It also drives interspecific relationships, thus affecting ecological network dynamics 1 . Size-biased species loss has profound effects on ecosystem function 7,8 . Ectotherms rely on ambient conditions to achieve operational body temperatures 9 . Because higher ambient temperature increases metabolic rates and the associated costs for a given body size 2 , global climatic warming is expected to drive shifts to communities consisting of smaller species 3 .
Our planet is urbanizing quickly 10 , which is a primary example of human-induced rapid environmental change. Cities are urban heat islands characterized by increased temperatures that are decades ahead of global averages 11 . Not only are cities warmer than surrounding areas, but they also experience extensive fragmentation of (semi-)natural habitats, and both of these effects increase with percentage built-up cover (BUC; a proxy for urbanization) 12,13 . This provides an opportunity to study the opposing effects of sizedependent thermal tolerance and dispersal capacity, as larger body size favours dispersal in some, but not all, taxa.
Publisher: NPG; Journal: Nature: Nature; Article Type: Biology letter ms: 2017-07-09177B Here we test the hypothesis that urbanization causes shifts in community-level body size, and that these shifts are dictated by the community-specific association between body size and dispersal. We generally expect the urban-heat-island effect to drive shifts to species with smaller body sizes in communities of ectothermic species, in line with Atkinson's temperature-size rule 14 . For taxa characterized by a positive association between body size and dispersal, however, we also expect a filtering in favour of larger-bodied species associated with habitat fragmentation 5,15 . Filtering for increased mobility has been demonstrated for urban ground beetle and plant communities 16,17 . Hence, for taxa characterized by a positive body-size-dispersal link, we predict that the general communitylevel pattern of smaller species with increasing urbanization may be neutralized or even reversed.
To test our hypothesis, we engaged in an analysis of community-level shifts in body size across a broad range of both terrestrial and aquatic taxa along the same systematically sampled urbanization gradients. We studied the direction of change of community-level body size in ten taxa using a replicated, highly standardized and nested sampling design that covers urbanization gradients at seven spatial scales (50-3,200 m radii; Fig. 1). We sampled each taxon at up to 81 sites, sampling 95,001 individuals from 702 species, with species-specific body size varying by a factor of 400 (0.2-80 mm; Extended Data Table 1). Three of the ten groups are characterized by a positive association between body size and dispersal capacity (see Extended Data Table 1).
We show that the local temperature of pond, grassland and woodland habitats significantly increases with urbanization (linear mixed regression models, P < 0.002; Extended Data Table 2). The intensity of these urban-heat-island effects is consistently larger during night and summer, in accordance with slower night-time city cooling and higher irradiation levels in summer 18 (Fig. 2, Extended Data Fig. 1, Extended Data Table 2). We also show that increased urbanization is linked to significant declines in habitat amount and the patch size of terrestrial habitats, and significant increases in distances among patches for both terrestrial and aquatic habitats (Pearson's r correlations, P ≤ 0.020; Extended Data Fig. 2).
Confirming our metabolism-based prediction that interspecific mean body size decreases with increasing temperature, urban communities for four out of the seven taxa (ground spiders, ground beetles, weevils and cladocerans) that did not have a positive sizedispersal link display reduced community-weighted mean body size (CWMBS). For ostracods, bdelloid rotifers and web spiders, no relationship with urbanization is found. By  Table 3). The positive shifts in size observed for these taxa are in line with our prediction that increased urbanization-mediated habitat fragmentation selects for larger species in taxa with positive size-dispersal links.
The Benjamini-Hochberg procedure 19 , which controls for false positives, confirms that all seven responses are significant at the study-wide level. Comparing the percentage changes in body size over a percentage BUC gradient of 0-25% shows a marked difference between taxa with a positive size-dispersal link (13.6% ± 8.3% (mean ± s.e.m.) body size increase) versus the other taxa (15.6% ± 5.3% body size decrease) (weighted two-sided analysis of variance (ANOVA): F 1,8 = 12.38; P = 0.0079). These community-level shifts in body size occur independently of shifts in species abundance and diversity along the urbanization gradients. For example, reduced diversity is apparent for taxa that display positive and negative size shifts, as well as for web spiders that lack a size shift. By contrast, cladocerans show size reduction without diversity change (Extended Data Table 4). For butterflies, macro-moths and orthopterans (that is, taxa with a positive size-dispersal link), the increase in the CWMBS ranges from 7% to 21% depending on the taxon, whereas size reductions of ground beetles, weevils and ground spiders (that is, terrestrial taxa with nonpositive size-dispersal links) range from −18% to −21% over an urbanization gradient of 0-25% BUC (Fig. 3). The cladocerans display the largest size reduction (−44%), in accordance with the temperature-size response generally being stronger in aquatic species than in terrestrial species as a result of the greater oxygen limitation in water 20 . However, the size reduction for the ostracods is much smaller (-13%) and non-significant (linear mixed regression model, P = 0.10), and for the rotifers no size shift is found. The absence of a size shift for the microscopic rotifers might indicate that their small size allows for sufficient oxygen exchange between warm, low-oxygen environments and body tissues, so that no community shift to smaller body sizes is induced by increased temperature. The absence of a size shift for web spiders may be explained by behavioural flexibility in their extended phenotype, as modified web designs help web-spider communities to adapt to urbanizationinduced lower average body size of aerial dipteran prey 21 .
Our multi-scale approach allows the pinpointing of the spatial scales at which urbanization best explains the observed effects. During winter, the urban-heat-island effect fades with increasing spatial scale during the day but not at night, whereas during summer both diurnal and nocturnal urban-heat-island effects are more pronounced at small scales (  Table 2). The spatial scale at which most of the variation in CWMBS is explained varied considerably among taxa, with effects for smallersized taxa prevailing at small spatial scales (Figs. 3, 4, Extended Data Table 3).
Urbanization induces biodiversity loss and biotic homogenization 10 (see also Extended Data Table 4). Here, we demonstrate that urbanization also leads to community-wide shifts in body size for the majority of studied species groups. The size reductions within aquatic and terrestrial taxa follow metabolic rules in line with the urban-heat-island effect, especially as our data on various pollutants suggest no correlation with percentage BUC (data not shown).
By contrast, the increased fragmentation that is a result of urbanization appears to cause size increases for taxa with positive size-dispersal links. Hence, our multi-taxa study provides evidence of bi-directional shifts in community body size. In addition to the interspecific patterns reported here, shifts in body size can also occur at the intraspecific level, through phenotypic plasticity and genotypic change 22-24 . Our results should enable mechanistic studies that elucidate the cause of the variation in the observed shifts in body size along urban gradients and quantify their functional effects in urban ecosystems. A better insight into the mechanisms behind shifts in body size will allow prediction of the intertwined effects of climate change and urbanization on the body-size distribution of communities.
The size-biased species loss reported here is expected to strongly affect ecosystem function 7,8 . If taxa in urban areas are represented by smaller or larger species, ecosystem structure and function will be affected in several ways. Metabolic theory and a recent artificial-selection experiment predict that shifted size distributions affect whole-ecosystem properties such as primary productivity, carbon cycling and decomposition 25,26 . Shifts in body size also translate into altered life histories, demographic rates and interspecific relationships 1,2 . For example, consumer-resource dynamics have recently been modelled for warming-related intraspecific size shifts mediated by phenotypic plasticity 27 . A clear-cut effect of shifts in body size on ecosystem function can be predicted for cladoceran zooplankton. Smaller-sized cladoceran communities are typified by reduced densities of large Daphnia species (highly efficient filter feeders that consume phytoplankton), and are thus less able to maintain top-down control on algal blooms than larger-sized communities 28 . Also, the observed shifts in macro-moth body-size distributions may be functionally linked to flowering plant diversity through pollination 29,30 .
The shifts in body size that we observe across a range of animal taxa will be directly relevant to future efforts to understand, predict and mediate population resilience, trophic Publisher: NPG; Journal: Nature: Nature; Article Type: Biology letter ms: 2017-07-09177B interactions, and ecosystem function in urban ecosystems 31,32 . Such insights will be essential to design the biodiverse towns and cities of the future. For example, urban planners could mitigate the micro-climatic effects and habitat fragmentation that result from urbanization with measures implemented at multiple spatial scales. Such interventions could involve the creation and/or modification of urban ponds and urban green infrastructure to increase the amount and quality of habitats 33 . Doing so would reduce the urban-heat-island effect and favour dispersal, and hence gene flow, in urban animal populations. Our results indicate that such impacts would maintain variation in the body-size distributions of urban communities and potentially mitigate the effect that shifts in body size may have on ecosystem function.  19310-19314 (2012). Basemap, copyright Esri.
Temperature averages were analysed in relation to site-specific percentage BUC values using linear mixed regression models. We calculated mean changes in temperature over a 0-25% BUC gradient, on the basis of model-estimated intercepts and slopes. Error bars represent the range of temperature change on the basis of these slopes with 95% confidence intervals.

Sampling design
Sampling was performed according to a nested design in which a local urbanization gradient (three classes: non-urban, semi-urban and urban) was repeatedly sampled within landscapes distributed along a landscape-scale urbanization gradient (three classes: non-urban, semiurban and urban). For each of ten taxa a total of up to 81 local-scale subplots (200 × 200 m) were sampled within 27 landscape-scale plots (3 × 3 km) situated in an 8,140 km 2 study area in northern Belgium (Fig. 1, Extended Data Table 1). The average human population density of this highly urbanized area amounts to 693 individuals per km 2 , with cities and urban sprawl embedded within an agricultural and semi-natural matrix 34 . As a proxy for urbanization we used percentage BUC, which was assessed in a geographic information system (GIS) using an object-oriented reference map of Flanders with the precise contours of all buildings, excluding roads and parking infrastructures, as a vectorial layer. Given that only buildings are considered, 15% BUC can be considered highly urbanized. Within each of the nine urban (BUC > 15%), nine semi-urban (5% < BUC < 10%) and nine non-urban (BUC < 3%) plots, one urban, one semi-urban and one non-urban subplot were chosen using identical BUC cutoff values, for a total of 81 subplots. Within each subplot, and for each of the ten taxa, a single grassland, woodland or pond habitat patch was targeted for sampling during the most appropriate season for each taxon (Extended Data Table 1). As each taxon was sampled in only one of three habitat types (that is, grassland, woodland or ponds), it was often impossible to sample all taxa within the same 200 × 200 m subplots. As such, independent subplots containing the corresponding habitats were sometimes selected among taxa, but these subplots were always of the same urbanization level and located within the same 3 × 3 km plot.
The classification of urban, semi-urban and non-urban (sub)plots on the basis of BUC cut-off values was used to establish the nested sampling design, which allowed samples to display a wide range of urbanization values at both local (subplot) and landscape (plot) scales.
To increase precision in the data analysis, however, we moved away from having BUC as a class variable with three levels, and instead quantified BUC as a continuous variable, at seven spatial scales around the sampling site (50, 100, 200, 400, 800, 1,600 and 3,200 m radii).
Owing to our nested design, BUC values at small scales were not correlated with values at large scales, hence allowing the pinpointing of the scales at which the effects of urbanization are most pronounced.
Publisher: NPG; Journal: Nature: Nature; Article Type: Biology letter ms: 2017-07-09177B Using this highly replicated, nested sampling design, our sampling effort involved counting and assigning 95,001 individuals to 702 species in ten taxa: (i) aquatic: cladocerans and ostracods sampled in pond habitats; (ii) limno-terrestrial: aquatic bdelloid rotifers sampled within the water layers of terrestrial Xanthoria lichens; and (iii) terrestrial: butterflies, orthopterans (that is, grasshoppers and bush crickets), macro-moths, ground beetles, weevils, web spiders and ground spiders sampled in grassland and woodland habitats (Extended Data Table 1).

Urban-heat-island effect
The urban-heat-island effect was quantified using hourly temperature readings that were collected automatically across 104 sampling sites for the three habitat types in which the ten taxa were sampled: ponds, grasslands and woodlands. March.

Habitat fragmentation
Correlations between urbanization (BUC) and three habitat fragmentation variables (that is, habitat coverage, mean size of habitat patches, and mean nearest-neighbour distance among habitat patches) were quantified using Pearson's r coefficients (Extended Data Fig. 2). This was done at a 3 × 3 km plot scale, on the basis of detailed land-use data from all 27 sampling plots (Fig. 1), and separately for terrestrial (that is, all types of (semi-)natural habitat) and aquatic habitat (that is, all pond types) 35,36 . Eutrophied, mono-specific intensive grasslands as well as orchards, plantations and conifer woodlands were not retained for analyses.

Statistical analyses
Temperature averages were analysed in relation to site-specific urbanization (BUC) values and habitat type (grassland, woodland and pond) using linear mixed regression models (R Publisher: NPG; Journal: Nature: Nature; Article Type: Biology letter ms: 2017-07-09177B package lme4). We ran separate models for both seasons (summer and winter) and for both day and night conditions (diurnal and nocturnal). Site ID and date (nested within year) were included as random factors. We used a multi-scale approach, running separate models with BUC values quantified at seven spatial scales (50-3,200 m radii). P values for the fixed effects were obtained using likelihood-ratio tests of nested models that were fitted with maximum-likelihood and parameter estimates from restricted maximum-likelihood models.
Residual plots were always visually inspected to evaluate the fit of models, and we compared maximum-likelihood-based AICc values (R package AICcmodavg) to select a confidence set of models whose AICc values did not differ substantially from the value of the best-fitting model, using ∆AICc ≤ 2 as a criterion 37 .
CWMBS was calculated for a given site as the average of the species-specific body sizes (mm) for all locally sampled species, weighted by species abundance. The raw data for calculating this metric are species-level count data for all taxa in all sites (based on taxonspecific sampling and identification protocols) and mean species-specific body-size values extracted from the literature or, in the case of web spiders and cladocerans, from our own measurements (Extended Data Table 1). An increase in CWMBS with increasing urbanization implies that the species assemblage of the site is increasingly composed of individuals belonging to larger species along the gradient from communities in more rural sites to communities in more urban sites. Our CWMBS index hence reflects the relative composition of large versus small species in local communities, and we use it here to quantify community response to urbanization. Although every sampling method introduces some bias in relative species abundances, the extent of the bias should be similar for non-urban and urban sampling sites. Therefore, using the relative species abundances that we obtained via sampling to calculate the CWMBS is appropriate to look into the relative effects of urbanization.
CWMBS was analysed for each taxon in relation to site-specific urbanization (BUC) values using linear mixed regression models with restricted maximum-likelihood estimation (R package lme4). Plot ID was used as a random variable to account for potential spatial autocorrelation of variables among sites belonging to the same landscape-scale plot. CWMBS values were log 10 -transformed for cladocerans and ostracods. For ostracods, we also transformed BUC values by taking the arcsine of square-rooted BUC values, which resulted in residual plots with a more homogeneous distribution. Analyses for the other taxa were run with untransformed data as residual plots proved to be homogeneous. The residual plots for orthopterans, ostracods and ground beetles each displayed one outlying data point, and the Publisher: NPG; Journal: Nature: Nature; Article Type: Biology letter ms: 2017-07-09177B residual plot for weevils displayed two such points. Because these five data points are legitimate (that is, they are not due to measurement, data or sampling errors) we assessed their effect on the consistency of the regressions in the model output. Filtering these data points out of the regressions showed (i) that the best-fitting models remained linked to the identical spatial scales; (ii) that the positive slope for orthopterans remained positive and the negative slopes for the other taxa remained negative, and (iii) that the significance levels stayed equal for ground beetles and ostracods, got stronger for weevils, and decreased but remained significant for orthopterans. Because those five data points are legitimate and do not have a qualitative effect on the output, we opted to retain them in the analyses. We used a multi-scale approach, running separate models with BUC values quantified at seven spatial scales (50-3,200 m radii). For each taxon, we then selected the model (and hence the spatial scale) that fitted the data best using maximum-likelihood-based AICc values (R package AICcmodavg).
We also retained a confidence set of the models whose AICc values did not differ substantially from the value of the best model using ∆AICc ≤ 2 as a criterion 37 .
For each taxon, and at the spatial scale of the best-fitting model, we calculated the percentage change (with 95% confidence interval) in CWMBS over a 0-25% BUC gradient, on the basis of the modelled intercept and slope, or of back-transformed values for ostracods and cladocerans (Fig. 3). These values were then contrasted for taxa with a positive sizedispersal link against all other taxa using two-sided ANOVA, with the inverse of the taxonspecific error bars as weights to account for the differences among taxa in the variance of the estimated percentage change. This weighted ANOVA allows testing of the percentage change values for taxa with a positive size-dispersal link to determine whether they are significantly different from those of all other taxa.

Reporting summary
Further information on experimental design is available in the Nature Research Reporting Summary linked to this paper.

Data availability
The datasets generated and analysed during the current study are available from the corresponding author upon reasonable request.  Table 2 | Model output of average temperature in relation to urbanization and habitat type Output of linear mixed regression models that test the relationship between average local ambient temperatures and the interaction between percentage BUC (%BU) and habitat type (pond, grassland and woodland) (n = 104 sites). Only the output for the confidence set of models (∆AICc ≤ 2) is given, with scale referring to the associated radius scale (in metres) of percentage BUC. Model estimates (±s.e.m.) for percentage BUC regression coefficients are provided. Model output consistently shows clear temperature differences among habitats and a clear positive effect of urbanization on temperature, irrespective of habitat type (as shown by the non-significant interactions).

Extended Data Table 3 | Model output of CWMBS in relation to urbanization
Output of linear mixed regression models testing the relationship between CWMBS and percentage BUC at multiple scales (radii in metres), for ten taxa (n = 76, 12, 75, 80, 62, 60, 81, 81, 68 and 80 biologically independent communities, top to bottom). Confidence sets of models (∆AICc ≤ 2) have grey shading and the best-fitting model output is given in bold. Modelled intercepts and slopes (±s.e.m.) are given. Credits for animal outlines are listed in the legend of Fig. 3.