Replicated, Urban-Driven Exposure to Metal Pollutants in Two Passerines

While there are increasing examples of phenotypic and genotypic differences between urban and non-urban populations of plants and animals, few studies identied the mechanisms explaining those dissimilarities. The characterization of the urban landscape, which can only be achieved by measuring variability in relevant environmental factors within and between cities, is a keystone prerequisite to understand the effects of urbanization on wildlife. Here, we measured variation in bird exposure to metal pollution within 8 replicated urbanization gradients and within 2 agship bird species in urban evolutionary ecology: the blue tit (Cyanistes caeruleus) and the great tit (Parus major). We report on a highly signicant, positive linear relationship between the magnitude of urbanization – inferred as either tree cover, impervious surface cover, or an urbanization score computed from several environmental variables, and copper, zinc and lead concentrations in bird feathers. The reverse relationship was measured in the case of Hg while Cd and As did not vary in response to the urbanization level. This result, replicated across multiple cities and two passerine species, strongly suggests that urban metal pollution is likely to trigger the emergence of parallel responses at the phenotypic and/or genotypic level between urban environments worldwide.


Introduction
There is increasing evidence that urbanization is associated with modi cations to plant and animal communities and populations [1][2][3][4][5][6][7][8][9] . For instance, birds inhabiting urbanized environments tend to suffer from lower reproductive outputs 4 and higher physiological costs (e.g., higher oxidative stress 10 and shorter telomeres 11 ). While chemical, light and noise pollution, as well as human presence and altered food availability and quality are often suggested as the main potential drivers of those phenotypic changes in urban areas [12][13][14] , few studies actually identi ed the environmental factors responsible for such modi cations, e.g., 15,16 . One of the reasons is that the majority of studies focussing on wildlife ecology and biology in the urban space uses an extremely simpli ed urban ecology framework that often lacks (i) adequate replication stemming from a comparison of multiple cities and (ii) knowledge stemming from multiple and contrasted urban habitats contributing to the urban mosaic. While awareness of these limitations is recently growing 17,18 , all too often urban ecology inference focuses on the phenotypes of individuals caught within one location in a city (often urban parks) and outside of this same city (usually in forests adjacent to the city) 19 . Thus, such experimental design ignores the diversity of environments within and between urbanized areas 20 , which prevents from i) establishing rm conclusions about the effects of urbanization per se on wildlife, ii) disentangling the effect of the different abiotic and biotic environmental factors that vary in response to urbanization 21 , and iii) drawing universal conclusions about the impact of urbanization on the biology of wild organisms at a continental or global scale. For this reason, recent reviews in urban ecology urge future research to focus on replicated and continuous gradients of disturbance 2,3,14,19 . Indeed, measuring how potential environmental stress factors vary within the urban mosaics across multiple cities is without a doubt a prerequisite to further understand the effects of urbanization on wildlife.
Metallic/metalloid Trace Elements (MTEs; e.g., lead, cadmium, copper) are a major class of pollutants that may have lethal and sublethal effects on organisms 22 . In birds, individuals nesting close to metallurgic smelters show reduced reproductive outputs [23][24][25] . While less documented, similar trends were measured in bird exposed to urban MTE pollution 16,26 . MTEs are mainly emitted by anthropogenic activities 27 . In urban environments, road tra c, residential heating, coal burning, and industry are the main sources of MTE pollution 28 . Literature on bird exposure to MTEs is abundant, although only a minority of studies focussed on urban MTE pollution 29 , and all of those studies but three 15,16,30 compared MTE concentrations in individuals using the coarse dichotomy of urban versus non-urban areas [31][32][33][34][35][36][37][38][39][40] . Those studies measured MTEs in feathers, blood, liver, kidney, bones or eggshell. While most studies reported higher levels of lead at urban sites than at non-urban sites, the impact of urbanization on other MTEs varied 15,16,30−40 . All in all, while MTE pollution might be a signi cant driver of phenotypic and/or genotypic changes triggered by urbanization, we currently lack knowledge on MTE pollution levels within complex urban-rural gradients.
To ll in this gap, we measured 6 of the most common MTE pollutants (i.e. copper, zinc, lead, cadmium, arsenic and mercury) in the feathers of 179 individual blue tits (Cyanistes caeruleus) and great tits (Parus major). Importantly, feathers were proven to be one of the most reliable non-invasive material for the biomonitoring of MTE exposure in birds 16,29,41 . Birds have been sampled in a continuous gradient of urbanization replicated across 8 cities (i.e. densely populated areas with more than 50,000 inhabitants 42 ) in Poland. Here, we de ne urbanization using high-resolution environmental data. This quantitative approach was compared with a qualitative approach where sampling sites were sorted into 5 habitat categories. Thanks to this unique experimental design, we address whether bird exposure to MTEs, captured for two passerine species, varies consistently, linearly and in a replicated fashion along multiple urbanization gradients.  1). Within each city, individuals were sampled in 3 distinct urban environments, namely the city centre, a residential area, and an urban park. For the city of Warsaw, those 3 habitats and the suburban forest have been additionally replicated three times. Birds were attracted to a mist-net using a loud-speaker playing calls from the two focal species as well as a dummy of a great tit. The birds were aged (one-year-old or older) and sexed based on their plumage features 43,44 . Additionally, from each bird, the second tail feather from the left side was plucked and stored in individual paper bags until MTE analyses. The protocol was performed in accordance with the Directive 2010/63/EU of the European Parliament and of the Council of 22 September 2010 on the protection of animals used for scienti c purposes. Moreover, this study was approved by the Local Ethical Committee nr I for Animal Experimentation in Warsaw (I Lokalna Komisja Etyczna ds. Doświadczeń na Zwierzętach w Warszawie).

MTE quantitative analysis
The feathers of 179 males (97 blue tits and 82 great tits) were analysed for their MTE content, out of a larger dataset of 350 individuals (140 male and 56 female blue tits, and 111 male and 43 female great tits) sampled in 2017 within the project. Feather selection was based on two criteria: rst, due to contrasted dispersal strategies (males dispersing shorter distances from their natal or previous breeding site than females [45][46][47][48] we standardized our data set by selecting only males; those represented over 70% of the individuals that were sampled. Second, out of the 251 feathers from males, a sub-sample of 179 feathers were selected in a way to maximise dataset balance in order to generate comparable sample sizes between species per habitat category (X² = 1.44, df = 4, P = 0.838; Table D1) and age class per habitat category (X² = 6.51, df = 4, P = 0.164). Feathers were prepared for metallic/metalloid trace element (MTE) analyses using the protocol from 16 . The following MTEs were quanti ed: lead, zinc, copper, cadmium, arsenic, mercury. Brie y, feathers were pooled and washed alternatively with 0.25 M NaOH solution and ultrapure water (Milli-Q puri ed, Merck KGaA, Darmstadt, Germany) to remove external contamination, then dried 12 h at 50°C to dry mass. Feathers were digested in 1 mL of HNO 3 30% for 24h at 80°C. The product of digestion was transferred into plastic tubes and ultrapure water was added to reach a nal 1% acid concentration. Spectrometer, Perkin Elmer SCIEX, USA). A conventional Mainhardt nebulizer and a quartz cyclonic spray chamber were used for sample introduction. Each isotope was measured three times and each sample was analysed two times. For concentrations above quanti cation limits, measurements with relative standard deviation above 10% were excluded. The ICP-MS was calibrated before performing measurements with the use of multi standard solutions (ICP Calibration Standard from Merck). During the measurements, the parameters of calibration were checked using the standard containing mercury at the concentration of 1 µg/L in 1% nitric acid. The blanks and Certi ed Reference Materials (CRMs; trace elements in water 1643f from LGC Standards and SPS-SW1 batch 112 from SpectraPure Standards) were prepared and analysed using the same methods as the samples. The recovery of the CRMs ranged from 90 to 110%. Concentrations measured in the blank were extremely low: they were 1.401 ppb for zinc, 0.508 ppb for arsenic and below quanti cation limits for copper, lead, cadmium and mercury. All measurements were performed at the Biological and Chemical Research Centre (NCBCh, University of Warsaw, Poland).

Quantifying urbanization
Impervious surface cover and tree cover in a 100 m radius around each mist-net sampling point was extrapolated via satellite imagery following the method described in 19 . Brie y, tree cover (i.e. the percentage of trees) and impervious surface cover (i.e. the percentage of soil sealing and built-up areas) were downloaded from Copernicus Land Monitoring Services; the basic maps referred to 2015 and are of 20 m pixel resolution. Distance to the closest road and to the city centre were calculated using GIS 2.8.2 and Google maps, respectively. For each city, the coordinates of the city centres were extracted from Wikipedia. Urbanization was quanti ed as (i) tree cover, (ii) impervious surface cover and (iii) an urbanization score computed from a Principal Component Analysis on the four environmental variables that were measured in this study: tree cover, impervious surface cover, distance to the closest road and distance to the city centre. Those three quantitative indexes of urbanization are commonly used in urban ecology 19,49,50 . Tree cover, distance to the closest road and distance to the city centre were all positively correlated, while impervious surface cover was negatively correlated with the other environmental variables ( Figure A1). Based on the Kaiser-Guttman criterion, one component, hereafter named "Urbanization score", was retained in the PCA. It accounts for 61.4% of the variance in the data set; it is positively correlated to the impervious surface cover (r = 0.85) and negatively correlated to tree cover (r = -0.80), the distance to the closest road (r = -0.69) and the distance to the city centre (r = -0.68). While a composite multivariate metric (here urbanization score) is the most accurate index of urbanization, univariate metrics that are highly correlated to such multivariate metric (here tree cover and imperviousness surface cover) are preferred as they are unambiguous and readily comparable between studies 19 . In addition to the three quantitative indexes of urbanization listed above, we also categorized the environment where the birds have been caught into 5 habitats with contrasted environmental features and land use; when arranged from the lowest to the highest urbanization level, the habitats are ordered as follows: protected forest, suburban forest, urban park, residential area and city centre (Fig. 2).

Statistical analyses
Statistical analyses were performed using R software (version 4.0.3) 51 . The percentage of data above the MTE quanti cation limit were similar between blue tits and great tits. Overall, they were 100% for Zn, 99% for Cu, 92% for Pb, 91% for As, 59% for Cd and 21% for Hg. Data below the quanti cation limits were given a value using a regression on order statistic ('ROS' function of the 'NADA' package); this estimation was done separately in blue tits and great tits. To minimize the in uence of possible spurious outliers on the distribution of MTE concentrations, values more than 1.5 times the interquartile range from the quartiles (i.e. either below Q1 − 1.5IQR, or above Q3 + 1.5IQR) were removed; the procedure was done separately for blue tits and great tits. In total, 16 values (ca. 1.5% of the values) were identi ed as outliers.
MTE concentrations in blue tit and great tit feathers were compared using linear mixed-effects models with MTE concentrations (i.e. Cu, Zn, Pb, As, Cd or Hg) as the response variable and urbanization level (computed as either impervious surface cover, tree cover, urbanization score or habitat category) as the explanatory variable. MTE concentrations may vary differently in response to bird moulting pattern, foraging behaviour and/or seasonal movements. Those are known to differ between species and age 43,44,52−54 . Therefore, species, age, and their interactions with the urbanization index were added as explanatory variables and the single terms "species" and "age" were considered as covariates. The location (i.e. either the city or the protected forest) was added as random intercept. Lmer were tted with restricted maximum likelihood (REML) method using the 'lme4' package. Normality of model residuals was validated using quantile-quantile plots. For each model, we performed a backward stepwise selection using the AIC 55 . A Type III Wald chisquare test Anova was used to determine the signi cance of retained variables in the nal models. When discrete explanatory variables were retained in the models, contrasts among groups were tested using least-square mean pairwise comparisons (contrast function of the 'lsmeans' package) 56 .
The proportion of variance in MTE concentrations that is explained by urbanization level (computed either from tree cover, impervious surface cover or urbanization score) was calculated using the 'r2_nakagawa' function of the 'performance' package. It was calculated using two metrics of relative importance: i) the difference between the conditional r-squared of the full model and the conditional rsquared of the model without the urbanization index as explanatory variable (de ned as the "last" metric in the 'relaimpo' package) and ii) the marginal r-squared of the model including the urbanization index only as explanatory variable (de ned as the " rst" metric in the 'relaimpo' package) 57 . The "last" and the " rst" metrics tend to underestimate and overestimate, respectively, the variance explained by the variable of interest (here the urbanization index), meaning that the exact variance explained by the urbanization index falls between the two values computed from those two metrics 57 . To further investigate what environmental variable(s) better explain(s) the variation in MTE concentrations, the two metrics were also computed for the distance to the closest road and the distance to the city centre using similar models but replacing the urbanization index by either the distance to the closest road or the distance to the city centre (Table C2).

Variation in MTEs concentrations along continuous urbanization gradients
We recorded an unequivocally clear impact of replicated urbanization gradients on the concentration of multiple MTEs: whatever the species and age of the individuals, Cu, Zn and Pb increased while Hg decreased with increasing the urbanization level (i.e. with decreasing tree cover but with increasing impervious surface cover or urbanization score); Cd and As did not vary in response to the urbanization level (see Table 1 and Fig. 3 for results on tree cover; the results of the models using impervious surface cover and urbanization score as a proxy of urbanization level are presented in Table C1). Urbanization level explained a maximum of 16%, 30%, 12% and 8% of Cu, Zn, Pb and Hg variation, respectively (Table 1  and Table C1). Models including "distance to centre" or "distance to road" systematically tted worse than the models including one of the urbanization indexes (Table C1). Moreover, MTE concentration in feathers are species-and age-speci c: Cu, Pb and As were higher in the feathers of blue tits than of great tits (Table 1). Pb and Zn were higher in the feathers of one-year-old birds than of older birds, although the relationship was only marginally non-signi cant for Zn (Table 1). Table 1 Results of the best tting statistical models testing the link between MTE concentrations (Cu, Zn, Pb, Cd, As and Hg) and urbanization level (here tree cover) while taking into account the species, the age and the location, tested in 8 cities and 4 protected forests. The proportion of the variance in MTE concentrations that is explained by the urbanization level -r² -is comprised between two values computed from the metrics " rst" and "last" [57]. For Hg, the random effect "city" had a zero variance, preventing to calculate the coe cient of determination -R² -of the model; for this reason, we report for Hg the adjusted coe cient of determination from the linear model. Signi cant effects (P < 0.05) are highlighted in bold. Results are strikingly similar when using impervious surface cover and urbanization score as metric for urbanization quanti cation; these are reported in Table C1.

Variation in MTE concentrations between habitat categories
Cu, Zn and Pb exhibited considerable variation from one habitat to another, and between urban and rural sites: overall, they were higher in city centres and residential areas than in urban parks, suburban forests, and protected forests ( Table 2 and Fig. 4). For instance, Cu, Zn and Pb were ca. 74%, 82% and 135% higher in urban centres than in adjacent suburban forests, respectively (Fig. 4). Similar to the previous models, Cu, Pb and As were higher in the feathers of blue tits than of great tits ( Table 2) and Pb and Zn were higher in the feathers of one-year-old birds than of older birds ( Table 2). Raw MTE concentrations per species and per habitat are detailed in Table D1.

Discussion
For the rst time, we report on a highly signi cant, positive linear relationship between the magnitude of urbanization -inferred as either impervious surface cover, the non-tree cover or an urbanization score computed from several environmental variables, and copper, zinc and lead concentrations in bird feathers; moreover, the reverse relationship was measured in the case of mercury. Importantly, those relationships were similar when quanti ed in two passerine species, and for both one-year-old and older birds. Strikingly, the signi cant covariation between urbanization and speci c metallic trace elements (MTEs) in avian feathers were measured across 8 replicated urbanization gradients, which included 8 cities, 8 suburban forests and 4 protected forests. This result con rms the pervasive impact on urbanization on the presence of MTEs in wild organisms, irrespective of the fact that the 8 cities differed in terms of size and population density (Fig. 1). Interestingly, the city of Warsaw, where 80% of the buildings have been destroyed during the second world war, is undergoing a recent densi cation and expansion; this suggests that bird exposure to urban MTE pollution is likely to vary in similar ways in every urbanized landscape, including relatively recent ones.
When arranged from the least to the most urbanized site, habitat categories were ordered as follows: protected forest, suburban forest, urban park, residential area, and city centre (Fig. 2). Therefore, it is not surprising that we measured higher concentrations of copper, zinc and lead in the feathers of the birds sampled at the most urbanized habitats than in the birds caught in urban parks were not signi cantly different from the concentrations measured in their conspeci cs sampled in suburban or protected forests. This result, as well as the fact that urban parks are the least urbanized urban habitats, are a striking demonstration that comparing bird populations between urban parks and forests, a still common experimental design in urban ecology 19 , is not appropriate to study i) how MTE exposure varies in response to urbanization and ii) the ecology and evolution of animal populations.
Importantly, the fact that we measured a signi cant and consistent association between urbanization level and the concentration of copper, zinc, lead and mercury in the feathers of two bird species suggests that blue tits and great tits have limited movements between their birth and their rst breeding attempt (in the case of one-year-old individuals) or between two consecutive breeding seasons (in the case of older individuals). Indeed, in this study, the birds were caught at the beginning of their reproductive season (i.e. in March and April). Yet, MTE concentrations in their feathers re ect the concentrations of those same MTEs in the environment where the birds moulted (i.e. between June and September), which is the nest where they hatched (if rst-year breeders), or their previous breeding ground (for 2nd -year breeders or older) in most bird species with seasonal moulting, including the blue tit and the great tit 43 . This, in case the birds had moved during their natal or breeding dispersal to an environment that is characterized by a different urbanization level, we would not have been able to measure such consistent correlation between MTE concentrations in the feathers and the urbanization level. This is a crucial result as it implies that 1) urban areas generate a systematic "signature" of MTE concentrations that can be detected in avian tissues such as feathers, 2) dispersal patterns are likely to be considerably reduced in urban birds relative to those that hatched in a natural environment [45][46][47][48] , with non-trivial consequences to our understanding of gene ow and local adaptation in the urban space, and nally 3) this signature is a promising tool to infer natal and breeding dispersal within urbanized landscapes, e.g., 58 .
Altogether, those results suggest that the effects of urban copper, zinc and lead pollution on birds (e.g., on their survival, reproductive success, immunity, corticosterone levels and plumage colouration 15,16,26,59−62 ) are expected to increase from less to more urbanized environments. In line with this prediction, rates of phenotypic change are greater in urban environments compared with non-urban landscapes 2 and numerous studies measured phenotypic divergence between urban and rural populations 3 . Consequently, we argue that some of those phenotypic changes (e.g., oxidative balance, xenobiotic metabolism, immunity) are likely to be triggered by urban MTE pollution. For instance, in the great tit, genes involved in antioxidant defences are expressed at a higher rate in individuals sampled within the city of Malmö (Sweden) than in their counterparts sampled in an adjacent forest; those individuals also exhibited a higher expression of the genes involved in xenobiotic metabolism 63 . Similar patterns were measured in the white-footed mouse (Peromyscus leucopus) 64 . The toxicity of several MTEs, including lead, resulting partly from the fact that they induce oxidative damages 65 , urban MTE exposure may trigger such modi cation in gene expression rate 66 .
Phenotypic differences can result from adaptive or non-adaptive adjustments 1,67,68 . The exposure to urban MTEs being often associated with detrimental effects on bird tness 16,26,59 , and bird movements being limited within the urban mosaic, urban MTE pollution is likely to exert selective pressures on populations inhabiting more urbanized environments. Importantly, this study reveals that urban MTEtriggered selective pressure is likely to be directional and consistent across all urban sites, and thereby conducive to the emergence of parallel responses at the genetic and / or phenotypic level 17 . In invertebrates, MTE exposure, although not in a context of urbanization, has been con rmed to drive considerable genetic variations 69 : for instance, wild populations of common fruit y (Drosophila melanogaster) that have been exposed to MTE pollution are endowed with a duplicated gene coding for metallothionein, a protein involved in the binding, transport and detoxi cation of MTEs 70 . In the same way, we could expect urban MTE exposure to select for individuals assimilating lower amounts and/or excreting higher amounts of toxic MTEs. As a matter of fact, there is some evidence that more melanistic individuals transfer higher amounts of MTEs in their teguments (i.e. feathers and skin) 59,71 and are positively selected in environments polluted with MTEs 59,71−73 . Stronger antioxidant defences 65 or faster pace of life 69 may also be selected in more urbanized environments where MTE pollution level is high. All in all, the association between the urbanization level and bird exposure to copper, zinc, lead and mercury exposure measured in this study across multiple cities suggests that urban MTE pollution is likely to trigger the emergence of similar genotypes across multiple cities. Therefore, comparing phenotypic traits (e.g., tegument melanin-based colouration, MTE gastro-intestinal absorption rates, antioxidative response, metallothionein expression) and allele frequency at locus associated with the metabolism of MTEs (e.g., GCLC, GCLM and MT2A genes) 74,75 along replicated urbanization gradients appears as a unique opportunity to study parallel evolution mediated by urbanization 17 .
Declarations Figure 2 Mean ± se urbanization level, either percent tree cover, percent impervious surface cover or urbanization score (i.e. computed from imperviousness, tree cover, distance to the closest road and distance to the city centre) per habitat category. Signi cant differences of urbanization level between habitats are indicated by different letters.

Figure 3
Relationship between MTE concentrations (i.e. Cu, Zn, Pb, Cd, As or Hg after log-transformation; in ppm) and the urbanization level (here tree cover). For the x axis to positively correlates with the urbanization level, we highlight the percent of non-tree cover (100 -percent of tree cover). We highlight the concentration for each single individual (grey dots), the mean ± se concentration per percent of non-tree cover (in black) and the regression line (in blue) and its con dence interval (in grey). The statistical signi cance of the relationship is highlighted with asterisks. Species-speci c relationships between MTE concentrations and the urbanization level are displayed in Figure B1. Mean ± se MTE concentrations (Cu, Zn, Pb, Cd, As or Hg after log-transformation; in ppm) per habitat category (protected forest, suburban forest, urban park, residential area or city centre). Signi cant differences of MTE concentrations between habitat categories are indicated by different letters.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.