Epigenetics and the city: Non‐parallel DNA methylation modifications across pairs of urban‐forest Great tit populations

Abstract Identifying the molecular mechanisms involved in rapid adaptation to novel environments and determining their predictability are central questions in evolutionary biology and pressing issues due to rapid global changes. Complementary to genetic responses to selection, faster epigenetic variations such as modifications of DNA methylation may play a substantial role in rapid adaptation. In the context of rampant urbanization, joint examinations of genomic and epigenomic mechanisms are still lacking. Here, we investigated genomic (SNP) and epigenomic (CpG methylation) responses to urban life in a passerine bird, the Great tit (Parus major). To test whether urban evolution is predictable (i.e. parallel) or involves mostly nonparallel molecular processes among cities, we analysed both SNP and CpG methylation variations across three distinct pairs of city and forest Great tit populations in Europe. Our analyses reveal a polygenic response to urban life, with both many genes putatively under weak divergent selection and multiple differentially methylated regions (DMRs) between forest and city great tits. DMRs mainly overlapped transcription start sites and promotor regions, suggesting their importance in modulating gene expression. Both genomic and epigenomic outliers were found in genomic regions enriched for genes with biological functions related to the nervous system, immunity, or behavioural, hormonal and stress responses. Interestingly, comparisons across the three pairs of city‐forest populations suggested little parallelism in both genetic and epigenetic responses. Our results confirm, at both the genetic and epigenetic levels, hypotheses of polygenic and largely nonparallel mechanisms of rapid adaptation in novel environments such as urbanized areas.


| INTRODUC TI ON
Identifying mechanisms involved in rapid adaptation to novel environmental conditions is a central theme in evolutionary biology and a pressing concern in the context of global changes characterizing the Anthropocene (Malhi, 2017). The vast majority of studies investigating mechanisms involved in rapid adaptation to new environments have focused on phenotypic plasticity on the one hand and on genetic responses to selection on the other hand. At their crossroad, recent work underlines the potential role of epigenetics in rapid adaptation to new environments (Liu, 2013). In particular, environmental variations can induce differences in DNA methylation patterns and hence modulate gene expression and upper-level phenotypes (Duncan et al., 2014;Jaenisch & Bird, 2003). Such methylation-linked phenotypic variation can occur during an individual's lifetime, especially early on during the organism's development (Waterland & Jirtle, 2003;Weaver et al., 2004). Although methylation changes acquired across an individual's lifetime may often be nonheritable (Reik et al., 2001) but see (Crews et al., 2007;Janowitz Koch et al., 2016), epigenetically induced phenotypic shifts may nevertheless enhance individual fitness in new environments. Moreover, during the course of evolution, divergent genetic variants regulating epigenetic modifications may also be under selection, hence promoting the evolution of divergent epigenotypes and epigenetically linked phenotypic variation (Richards et al., 2010). While epigenetic studies focused on human diseases and medical topics are now abundant, studies in an ecological context are still rare (Derks et al., 2016). A few epigenetic studies in natural plant populations revealed that DNA methylation shifts might play a determinant role in local adaptation to environmental variation (Dubin et al., 2015;Foust et al., 2016), however regulation and effects of DNA methylation are quite different between plants and vertebrates and methylation studies in vertebrates are rare (McNew et al., 2017).
There is hence an urgent need for further empirical investigations of simultaneously rapid genetic and epigenetic evolution in response to environmental change (Danchin et al., 2011).
Urbanization rapidly and irreversibly changes natural habitats into human-made environments and is considered a major threat to biodiversity (Brondizio et al., 2019). For species who appear to cope with urbanization, urban habitats present a myriad of novel environmental conditions compared to the habitat where they evolved, including high levels of chemical, light and sound pollution, high proportion of impervious surfaces, high habitat fragmentation, low vegetation cover and high human densities (Grimm et al., 2008;Szulkin, Garroway, et al., 2020;Szulkin, Munshi-South, et al., 2020). Such extreme environmental changes compared to natural areas are expected to result in numerous novel selection pressures for city-dwelling species (Szulkin, Garroway, et al., 2020;Szulkin, Munshi-South, et al., 2020).
Accordingly, rates of recent phenotypic change, concerning multiple types of traits related to behaviour, morphology, phenology and physiology, were found greater in urban areas than in any other habitat types, including nonurban anthropogenic contexts (Alberti et al., 2017;Thompson et al., 2018). The exploration of the molecular mechanisms implicated in urban-driven phenotypic changes has only begun, with both genetic (Mueller et al., 2013;Perrier et al., 2018;Salmón et al., 2021), and epigenetic investigations (McNew et al., 2017;Riyahi et al., 2015;Watson et al., 2021). For instance, DNA methylation variations have been associated in vertebrates with high levels of traffic-related air pollution (Ding et al., 2017). Yet, epigenetic studies have been performed at relatively small genomic resolution. In addition, very little is known about the level of parallelism and hence of the predictability of genetic and epigenetic evolution in response to urbanization in distinct cities Santangelo et al., 2018). So far, a small number of studies have provided evidence for a range of situations: from local adaptation despite strong gene flow (e.g. in the Great tit: Salmón et al., 2021; in the red-tailed bumblebee Bombus lapidaries: Theodorou et al., 2018) to restricted gene flow and independent colonization in different cities by a few founders, followed by adaptation (e.g. in the burrowing owl Athene cunicularia, Mueller et al., 2018).
Providentially, recent genomic tools of high resolution and the multitude of cities around the globe offer unique opportunities to compare simultaneously individuals' genomic and epigenomic responses in several cities and thereby study the parallelism and predictability in molecular mechanisms implicated in rapid adaptation to urbanization (Perrier et al., 2020;Santangelo et al., 2020).
In this study, we used both genome-wide and epigenome-wide sequencing approaches to compare genetic and epigenetic responses among three pairs of great tit Parus major populations in urban and forest habitats. At the European level, population monitoring of Great tits revealed parallel phenotypic shifts in city birds compared to their forest conspecifics, with in particular smaller and lighter urban birds laying earlier and smaller clutches Caizergues et al., 2021;Chamberlain et al., 2009;Corsini et al., 2020). In addition, genomic analyses showed patterns of genome-wide differentiation between urban and forest birds (Perrier et al., 2018) while a large-scale analysis revealed some parallel footprints of adaptation to urbanization across nine European cities . At the epigenetic level, a preliminary Great tit study recently described methylation shifts associated with urbanization (Watson et al., 2021). However, this analysis focused on a single location (n = 6 urban and 6 forest males) and hence could not test for a potential parallelism in the urban-related epigenomic response. In order to advance our understanding of the genome-wide and epigenome-wide responses to urbanization, and its putative spatial parallelism, we here searched for genomic footprints of divergent selection and for DMRs between three pairs of urbanforest populations across Europe, thereby complementing the work of Salmón et al. (2021) with an epigenomic insight into great tit adaptation to urbanization. Our results show that despite limited genetic differentiation and few genomic footprints of divergent selection between forest and urban populations, urban life was associated with numerous differentially methylated regions notably associated with neural development, behaviour and immunity. Hence, this study suggests that shifts in DNA methylation patterns could play a role in adaptation to urbanization. Importantly, we found little parallelism between cities in both the genomic and the epigenomic responses to urbanization, possibly confirming the hypothesis that multiple evolutionary ways exist to independently cope with similar novel environmental conditions.

| Study sites and sampling
Three pairs of great tit populations in urban and forest environments were sampled in the three European cities of Barcelona (Spain), Montpellier (France) and Warsaw (Poland), Figure 1. For each location, 10 individuals were sampled within the city and 10 individuals were sampled from nearby forest. Blood samples were collected from breeding individuals during spring between 2016 and 2018 (except 2 individuals for Barcelona city, collected in 2014 and and kept in 96% Ethanol or Queen's Lysis Buffer. Samples had balanced sex ratio (5 males and 5 females for each population) except for the forest population of Barcelona where 6 females and 4 males were sampled. The sample design did not include any closely related individuals (i.e. no full-or half-sibling nor parents/offsprings): (1) we chose to sequence individuals that were not full-sibs or half-sibs, based on field pedigrees and (2) we subsequently confirmed that neither full-sibs nor half-sibs were present in the dataset, by inspecting identity-by-state measured between all pairs of individuals using the RADseq data.

| DNA extraction, RAD-seq and reducedrepresentation bisulphite sequencing
We used QIAGEN DNeasy blood and tissue kits to extract DNA from blood samples, following the provided instructions for nucleated blood samples. DNA was quantified using a NanoDrop ND8000 spectrophotometer and a Qubit 2.0 fluorometer with the DNA HS assay kit (Life Technologies). DNA quality was examined on agarose gels.
We then performed RAD-sequencing and RRBS-sequencing using standard protocols. For RAD-sequencing (restriction-site-associated DNA sequencing, Baird et al., 2008), the library preparation was done

F I G U R E 1
Great tit blood sample locations in Europe (in urban and forest sites in and near Barcelona, Montpellier & Warsaw) Still using the population module, we kept SNPs with minor allele frequency >0.01, observed heterozygosity <0.65, and genotyped in at least 90% of individuals in each population. We obtained 181,041 SNPs. Then, using vcftools v0.1.15 (Danecek et al., 2011), loci with minor allele frequency <0.05 were removed, shortening the dataset to 75,246 SNPs. Still using vcftools, SNPs with extremely low or high coverage (5%-95% of the distribution) were removed, resulting in 74,137 SNPs retained for subsequent population genomic analyses.
To document genomic variation among urban and forest great tits from the three locations we used a redundancy analysis (RDA), with location (Barcelona, Montpellier and Warsaw), environment (urban or forest) and sex as explanatory variables. Partial RDA was also produced to test for each variable effect (environment, location or sex) alone after controlling for all other variables. The effect of a given factor was considered significant with a p-value < 0.05.
To estimate genome-wide differentiation between populations we used Weir and Cockerham's F ST (Weir & Cockerham, 1984) computed using the StAMPP R package (Pembleton et al., 2013).
Average F ST was estimated using all SNPs, and confidence intervals were assessed using 1000 bootstrap replicates.
We used two methods to investigate outlier SNPs potentially under divergent selection between forest and urban populations: an F ST -outlier based method (using Bayescan v2.1, Foll & Gaggiotti, 2008), and a multivariate method (using a RDA, Forester et al., 2016) aiming respectively at identifying strong outliers indicating footprints of differentiation for each population pair and weaker outliers that represent footprints of divergent selection typically expected in polygenic adaptations in response to complex environmental heterogeneity across several population pairs studied at once.  Foll and Gaggiotti (2008) we considered SNPs as outliers when they displayed a q-value above the 0.1 threshold. Second, following a similar procedure as described above for the RDA analysis, we used a constrained RDA to investigate the effect of habitat (forest vs. city) and to identify outliers SNPs that displayed more than 3 times SD from the mean score on the constrained axis (Forester et al., 2016).

| Methylation calling and epigenomic statistical analyses
The RRBS reads were first trimmed using fastp software v0.19.7 (Chen et al., 2018), and quality filtered to keep only reads with a quality >15.
The BISMARK software v0.20.0 (Krueger & Andrews, 2011) was used for mapping reads on the masked reference genome with default parameters and a maximum of one allowed mismatch (see Table S1 for mapping and filtration details). Note that BISMARK implements methylation calls for the entire R1 and for the R2 until reads overlap, to prevent double counts. Methylation information for cytosines in a CpG context with sufficient coverage (≥10×) was extracted.
Similarly to genetic differentiation, an RDA was performed to describe epigenomic variation (using level of methylation measured at each CpG position of the genome) across location, habitat and sex. A partial RDA was also conducted to test for the habitat effect alone.
Additionally, we investigated more finely whether individual methylation on CpG cytosines varied across location, habitat (urban vs. forest), sex, and a location × habitat interaction, using an ANOVA, run on autosomes and Z chromosome separately.
We used the MethylKit R package to identify differentially methylated regions (DMRs) between groups of individuals (habitat or sex) for each location. With a logistic regression including sex as covariate ("calculateDiffMeth" function) on normalized coverage data ("normalizeCoverage" function), we looked for 1000 bp regions (produced by cutting the genome into 1000 bp units; mean ± SD = 33,682 ± 664 tiles per location; mean ± SD = 28.940 ± 17.906 CpG/tiles) differing by at least 10% (as usually reported in the literature) of methylation between urban and forest individuals (regions present in at least 9 over 10 individuals). Methylkit was used with the option destrand = TRUE, collapsing both strands because (1) it is recommended for studies that focus on CpG only and (2) do not explore hemimethylation. We kept only DMRs with a q-value < 0.001, as suggested in Methylkit's vignette, that were considered as significantly different between habitats. Similarly, we ran a logistic regression to identify DMRs between sexes, including habitat as covariable.

| Genes associated with SNP outliers and DMRs, and gene ontology analyses
We investigated whether genomic outliers and DMRs overlapped genes in 5 kb upstream and downstream regions using BEDTools v2.28.0 (Quinlan & Hall, 2010). Gene ontology analyses were performed with the R package topGO (Alexa & Rahnenführer, 2009), with GO referenced for the chicken Gallus gallus (ggallus_gene_ensembl) and using a background list of genes covered in our dataset, to identify potential statistically enriched genes ontologies among the lists of genes extracted. In order to take into account hierarchy and nonindependence in GO terms we used the "weight01" algorithm with a "fisher" statistic (options algorithm="weight01" and statistics="fisher" of the runTest() function). GO terms were considered significantly enriched when having a weightFisher <0.05.

| Genetic differentiation
The redundancy analysis (RDA) performed on 74,137 SNPs obtained by RAD-sequencing, including location (Barcelona, Montpellier or Warsaw), habitat (urban vs. forest), and sex as explanatory variables was highly significant (p < 0.001) but explained only a small fraction (i.e. less than 2%) of the total variance (R 2 = 0.018, Figure 2a, Table S2). All three variables were significant (location: p = 0.001; habitat: p = 0.004; sex: p = 0.001) revealing small genetic structuration between groups. Partial RDA revealed that the net variation explained by habitat (R 2 = 0.004, p = 0.001) was inferior to the net variation explained by location (R 2 = 0.012, p = 0.001) but higher than sex (R 2 = 0.002, p = 0.004, Table S2). As expected, when removing the Z chromosome from the data, sex became nonsignificant (p = 0.260), whereas the effects of other variables remained significant and of similar magnitude (Table S3).
Genome-wide differentiation between populations was relatively low on average (mean F ST = 0.019), in the order of 1%-2% between habitats for each location (F ST Barcelona(forest-urban) = 0.018 ± 0.001;

| Methylation variation
Similarly to genetic data, we performed an RDA on methylation level of 157,741 CpG sites to describe epigenetic variation among individuals in relation to location, habitat and sex (Figure 2b). The model was significant but explained less than 1% of the total variance (R 2 = 0.007, p = 0.001). All variables contributed significantly (location: p = 0.001, habitat: p = 0.03, sex: p = 0.001, Table S6).
Partial RDA revealed that location and sex explained a similar proportion of the total variance which was higher than habitat (location: R 2 = 0.003, p = 0.001; sex: R 2 = 0.003, p = 0.001; habitat: R 2 = 0.001, p = 0.3). When removing the sex chromosome from analyses, results remained similar (Table S7), showing that the difference in methylation was not entirely driven by sexual chromosomes.
We then investigated more finely whether individual methylation on CpG cytosines varied across location, habitat (urban vs. forest), sex, and location × habitat interaction, using an ANOVA, run on autosomes and Z chromosome separately (Figure 3, Table S8). For autosomes, we detected a significant effect of location (F = 3.319,

| Nonparallel yet strong genomic footprints of divergent selection between urban and forest populations and evidence for polygenic adaptation
First, Bayescan identified 15 outliers for Barcelona, 11 for Montpellier and 10 for Warsaw, distributed across 15 chromosomes and associated with 13 genes in 5kb upstream or downstream regions ( Figure 4, q-value < 0.1, see Figure S1). None of these outliers was shared between the three population pairs, revealing no convergence between cities. Second, the multivariate approach based on F I G U R E 2 Rdancy analyses (RDA) on (a) genomic data (74,137 filtered SNPs) and (b) methylation levels (based on methylation levels observed at 157,741 positions). Triangles represent forest habitats, circles represent urban habitats, empty and solid symbols represent females and males respectively. ***p-value < 0.001, **p-value < 0.01 and *p-value < 0.05, related to the explanatory factors an RDA revealed a list of 1163 loci with outlier loading score ( Figure   S2) suggesting a high number of loci putatively undergoing weaker selection. These 1163 loci were associated with 561 genes. Detailed GO results are presented in Table S9 and Figure S3.

| Evidence for mostly nonparallel differentially methylated regions between urban and forest environments
We identified a total of 224 distinct DMRs between urban and forest great tits: 80 for Barcelona, 68 for Montpellier and 93 for Warsaw. Only 14 DMRs (6.25%) were found repeatedly in at least two comparisons, and only 3 were common to the three cities. 7 of these 14 parallel DMRs were in the same direction of methylation in urban compared to forest areas ( Figure 5, Figure S4a).
Barcelona urban birds presented significantly more hypomethylated than hypermethylated DMRs compared to Barcelona forest birds (χ 2 = 11.25, p < 0.001), while no difference was found for  Table 2 and Figure S5.
We also searched for DMRs between sexes, following the same procedure. We identified 206 DMRs associated with sex, of which   Figure 6 and Figure S6).
Almost twice more sex DMRs were shared between locations (11.7%) than between habitats (6.25%, see Figure S4a,b; z-test: When taking into account the direction of methylation difference, 7 sex DMRs were shared between at least two cities (9.7%) which was three times more than for habitat DMRs (3.1%; z-test: X-squared = 7.904, p = 0.005).
Finally, we also explored patterns of methylation associated with candidate genes linked to behavioural, morphological and phenological traits that are known to differ between forest and city great tits, and which were covered by our sequencing. While none of these genes shows strong methylation difference, we provide informative plots that could be used by other studies that wish to focus on these genes and provide informative plots that could be used by other studies that wish to focus on these genes (see Appendix 1 and Figures S7-S11 for CLOCK, COL4A5, DRD4, EGR1 and FOXP2).

| DISCUSS ION
The urban sprawl is a worldwide phenomenon deeply affect-

F I G U R E 5 Circos plot of differentially methylated regions (DMRs) identified between populations of forest and urban great tits in and near
Barcelona, Montpellier and Warsaw (from inner to outer circles). Red points show hypermethylated regions in urban great tits relatively to forest birds, and blue points show hypomethylated regions. For graphical clarity, only a subset of genes are represented: genes associated with the 10% most extreme DMR (triangles) and genes found associated with DMR in at least two cities (stars). Names of the genes found within 5 kb of the represented DMRs are given Overall genetic differentiation between populations was relatively low (F ST ranging from 0.009 to 0.034), although higher than what has been found at a much larger scale across the species distribution (e.g. F ST around 0.01 between the UK and Spanish or between French and Spanish populations, Laine et al., 2016). Low but significant differentiation levels are in line with previously documented genetic divergences between city and forest great tit populations (Perrier et al., 2018;Senar & Björklund, 2020), and altogether suggests important gene flow, large effective population sizes and limited genetic drift at multiple spatial scales (Kvist et al., 2003). This overall genetic context is particularly suitable to search for genomic footprints of divergent selection between urban and forest populations, which would easily be identifiable above the neutral level of genetic differentiation.
We found a limited number of strong footprints of divergent selection, which is in line with previous results in Montpellier (Perrier et al., 2018), and with the broader work of Salmón et al. (2021) across nine European cities. Similarly to low levels of parallelism in allele F I G U R E 6 Circos plot of differentially methylated regions (DMRs) identified between females and males great tits in and near Barcelona, Montpellier and Warsaw (from inner to outer circles). Red points show hypermethylated regions in female great tits relatively to males, and blue points show hypomethylated regions. For graphical clarity, only a subset of genes are represented: genes associated with the 10% most extreme DMR. Names of a subset of the genes found within 5 kb of the represented DMRs are given frequency changes between cities observed by previous studies (Reid et al., 2016;Salmón et al., 2021), and despite similarities in phenotypic shifts, none of these outliers were shared between cities. This result suggests limited parallel evolution, supporting a scenario of independent de novo evolution between cities and/or different selection pressures between cities. Indeed, there may be multiple evolutionary solutions to the same environmental challenges (Losos, 2011) and multiple traits are linked to the same functional outcome (Thompson et al., 2017). Besides, the identification of numerous outliers by the multivariate framework applied at the scale of all six sampling sites supports a model of polygenic urban adaptation implicating multiple genes, biological pathways and phenotypic traits (Boyle et al., 2017).
Polygenic adaptation is a reasonable expectation in urban habitats since the multiple new environmental conditions in cities most probably result in many novel selective pressures acting on a multitude of functional traits (Shochat et al., 2006), and because many of these traits may be quantitative, and genetically correlated (Lande, 1979).
Further polygenic analyses on more samples and more markers (i.e. whole genome data) are however required in order to estimate the potential effect of each genetic variant implicated in the adaptation to life in the city (Robinson et al., 2014;Zhou & Stephens, 2012).
Several genomic footprints of divergent selection between urban and forest environments were in, or in the vicinity of, genes that have already been described as playing a role in neuronal development, behaviour or cognitive abilities. In particular, the NR4A2 gene plays an important role in recognition of novel objects and memory in mice (McNulty et al., 2012). Reaction to novel objects and novel food is known as one of the main factors determining the capacity of a species to thrive in an urban environment (Lowry et al., 2013). The DCX gene is related to neuronal plasticity (Kim et al., 2006) and experimental approaches revealed that artificial light at night induces an overexpression of this gene linked to a change in behaviour and expression of depressive-like behaviour in crows (Taufique et al., 2018). Finally, the CHRNA1 gene is associated with aggressive behaviour in chicken (Buitenhuis et al., 2009), and higher aggressiveness is commonly observed and hypothesized as adaptive in urban bird populations .
Besides, the gene ontology enrichment analysis, performed on the entire set of genes identified via the outlier genome scan, reinforced these findings since multiple enriched GO terms were associated with the nervous system and stress response as well as hormonal response (Table 2). These results are informative on the type of traits involved in avian urban adaptation in cities and corroborate previous results from Sih & Del Giudice, 2012;Sol et al., 2013) suggesting that natural selection repeatedly acted on neuronal, behavioural and cognitive traits that could contribute to the phenotypic shifts described in urban great tits (i.e. more aggressive and exploratory birds, with higher breath rate; A. E. Caizergues, A. Grégoire, R. Choquet, S. Perret, & A. Charmantier, unpublished data;Senar et al., 2017;Torné-Noguera et al., 2014).
Contrary to the common prediction that living in cities is likely to influence epigenomes (McNew et al., 2017;Watson et al., 2021), no genome-wide pattern of differentiation in methylation between urban and forest great tits was detected. However, we observed a difference in mean methylation level between birds from Warsaw and Barcelona on their autosomes, as well as between males and females on the Z chromosome, showing that methylation differences were identifiable. In addition, we found no mean difference in methylation level between habitats, revealing that urbanization did not strongly affect overall methylation levels in a specific direction. This overall low differentiated methylation context is perfectly suitable to investigate more localized zones that could differ in their levels of methylation. Note that the strong methylation contrast between males and females on sex chromosomes (Figure 2) is in line with previous reports in vertebrates (Teranishi et al., 2001;Waters et al., 2018) showing that methylation plays a major role in sex differentiation via regulation of gene expression and genetic imprinting.
Despite the nonsignificant effect of habitat on overall methylation levels, we found a large number of DMRs between pairs of forest and urban populations, suggesting that urbanization did affect particular regions of the genome. DMR were less likely to occur within a gene body than by chance, but it was not the case for promotor or TSS regions. This latter result contrasts with Watson and collaborators (Watson et al., 2021)  Only a limited number of urbanization-linked DMR were shared between two or more locations (note that the particular sampling pattern in Barcelona (see methods) could potentially affect the results found for this city). In contrast, three times more sex-linked DMRs were found in two locations or more. This comparison suggests that urbanization-linked epigenetic modifications most probably do not occur in a parallel way across cities, but rather that each city might have its particular epigenetic response. Indeed, in the emerging field of urban evolutionary biology, cities are often regarded as valuable replicates of human-altered habitats (Donihue & Lambert, 2015;Santangelo et al., 2020), and it is often expected that parallel adaptive responses to similar selective pressures will occur.
This expectation is particularly strong when phenotypes show parallel changes, as is the case for the Great tit, which is consistently smaller and lays earlier and smaller broods in the various cities where it has been studied, compared to forest habitats (Caizergues et al., 2018;Seress et al., 2020). However, as discussed above, parallel adaptation to similar environmental conditions should not be expected in the case of independent evolution, especially for multilocus traits.
In fact, besides the obvious differences in cities' climatic conditions depending on their position on the globe, land use, fragmentation and pollutants levels can also largely vary across cities (Cárdenas Rodríguez et al., 2016). In a general way, pollutants are known to affect DNA methylation and result in both hypo-and hypermethylation, but the patterns of change observed largely rely on the pollutant involved (Head, 2014). Hence differences in cohorts of pollutants present in cities could be responsible for differences in patterns of methylations. Taken together, the results of the present study highlight the importance of questioning the assumption that cities are replicated environments that can be considered similar, and parallel evolution across cities may be the exception rather than the norm.
As mentioned earlier, increasing evidence suggests that DNA methylation can be associated with environmental and stress factors (env: Foust et al., 2016, stress: Sun et al., 2018 especially during early development (Meaney & Szyf, 2005). Here, we found four genes (POMC, ADAMTS3, PAPD4 et GCC1), associated with DMR that were previously described in great tits as undergoing major changes in methylation levels in case of exposure to higher levels of pollutants (Mäkinen et al., 2019). Notably, the functions of these genes remain to be determined, and they could thus be interesting to target in future studies. In addition, the past literature has repeatedly found SERT and DRD4 as two major genes involved in urban-specific avian human avoidance (or wariness) behaviours (see for example, in the Great tit (Riyahi et al., 2015), in the blackbird (Garroway & Sheldon, 2013), in the black swan (Van Dongen et al., 2015) and in the burrowing owl (Mueller et al., 2020) we found no DMR associated with these two genes in either of the three forest-city pairs. However, we found a significant urbanrelated change in methylation linked to the DRD3 gene, belonging to the same gene family as DRD4 and known to be similarly involved in chicken aggressive behaviour (Li et al., 2016). In line with these results, GO analyses revealed enrichment in genes associated with neuronal functions, behaviour, but also blood, immune and endocrine systems (Table 2, Figure S5), revealing the potential need of physiological adjustments in urban habitats. Surprisingly, a recent study on great tit differences of methylation between city and forest habitats in another European city found no GO enrichment in blood, while some in liver tissue (Watson et al., 2021; note that they investigated DMSs, Differentially Methylated Sites, which differs from DMRs identified here and do not provide the same information). These contrasted results highlight the fact that methylation patterns highly depend on the analysed tissues (Derks et al., 2016), and show, once more, that urban linked methylation might not be similar from one city to another. As reviewed by Husby (2020), the use of blood sample in methylation studies comes with several potential limitations. For instance, as previously mentioned (1)

ACK N OWLED G EM ENTS
We thank all the people that contributed to fieldwork and blood sample collection, in particular the members of the PLT platform of the CEFE, as well as the fieldwork teams of Barcelona and Warsaw.
We also thank Patricia Sourouille and the CEFE GEMEX platform for laboratory support, and Enrique Ortega-Aboud as well as Rémi Allio for bioinformatics support. Finally, we thank the editor and two reviewers for their constructive comments.

CO N FLI C T O F I NTE R E S T
We declare we have no competing interests.