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

Aude E. Caizergues, Jeremy Le Luyer, Arnaud Grégoire, Marta Szulkin, Juan-Carlos Senar, Anne Charmantier, Charles Perrier

► To cite this version:

Aude E. Caizergues, Jeremy Le Luyer, Arnaud Grégoire, Marta Szulkin, Juan-Carlos Senar, et al. Epigenetics and the city: Non-parallel DNA methylation modifications across pairs of urban-forest Great tit populations. Evolutionary Applications, 2022, 15 (1), pp.149-165. 10.1111/eva.13334. hal-03515127

HAL Id: hal-03515127

https://hal.science/hal-03515127

Submitted on 6 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Epigenetics and the city: Non-parallel DNA methylation modifications across pairs of urban-forest Great tit populations

Aude E. Caizergues<fnref n="1" /> ID

Jeremy Le Luyer<fnref n="2" /> ID

Arnaud Grégoire<fnref n="1" /> ID

Marta Szulkin<fnref n="3" /> ID

Juan-Carlos Senar<fnref n="4" /> ID

Anne Charmantier<fnref n="1" /> ID

Charles Perrier<fnref n="5" /> ID

## Correspondence

Aude E. Caizergues, CEFE UMR5175, 1919 route de Mende, 34293 Montpellier Cedex 5, France.

Email: audeemiliecaizergues@gmail.com

## Funding information

European Research Council: Starting grant ERC-2013-StG-337365-SHE (AC); OSU-OREME; Polish National Science Centre grants Sonata BIS 2014/14/E/NZ8/00386 and Opus 2016/21/B/NZ8/03082 (MS); Ministry of Economy and Competitivity, Spanish Research Council grant CGL-2016-79568-C3-3-P (JCS).

## 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.

## KEYWORDS

adaptation, DNA methylation, epigenomics, genomics, urbanization

<footnote n="1">¹CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote n="2">²Ifremer, IRD, Institut Louis-Malardé, Univ Polynésie Française, EIO, Taravao, French Polynesia</footnote>

<footnote n="3">³Centre of New Technologies, University of Warsaw, Warsaw, Poland</footnote>

<footnote n="4">⁴Museu de Ciències Naturals de Barcelona, Barcelona, Spain</footnote>

<footnote n="5">⁵CBGP, INRAe, CIRAD, IRD, Montpellier SupAgro, Univ. Montpellier, Montpellier, France</footnote>
# 1 | INTRODUCTION

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 (Rivkin et al., 2018; 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 (Biard et al., 2017; 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 (Salmón et al., 2021). 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 urban-forest 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.
2 | MATERIAL AND METHODS

2.1 | 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 2015) 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.

## 2.2 | DNA extraction, RAD-seq and reduced-representation 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 by the Montpellier GenomiX (MGX) platform (CNRS, Montpellier), using the enzyme Sbfl. Each individual was identified using a unique six nucleotides tag, individuals were randomly multiplexed in equimolar proportions by libraries of 37 individuals. Each library was sequenced on a lane of an Illumina HiSeq 2500. Paired-end sequencing was used to produce 150 bp reads. This generated an average of 4.9M reads per individual. The DNA of the 60 individuals were processed twice to test for reliability of the genotyping process. The RRBS-sequencing started with DNA digestion using MspI restriction enzymes, which cuts CCGG sites and target regions that are CG rich, permitting to have a high proportion of sequences in promotor regions. Individuals were randomly multiplexed in equimolar proportions by libraries of 10 individuals. Bisulphite treatment converted unmethylated cytosines into uracil, then converted to thymine after PCR amplification. Each library was then sequenced on a lane of an Illumina HiSeq 2500. Paired-end sequencing was used to produce 50 bp reads. This generated an average of 19.3M reads per individual.

2.3 | SNP calling and statistical genomic analyses

Fastp v. 0.19.7 (Chen et al., 2018) was used to trim the RAD-seq reads, keeping reads with a minimum quality of 15 before mapping individual sequences against the reference genome of the Great tit (Laine et al., 2016; GenBank assembly accession: GCA_001522545.3) with BWA v0.7.17 (Li & Durbin, 2010). Genotyping was conducted with stacks v2.41 (Rochette et al., 2019) "gstacks" and "population" functions, using "snp" model, filtering for mapping quality >10, and alpha = 0.05.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Great tit blood sample locations in Europe (in urban and forest sites in and near Barcelona, Montpellier & Warsaw)</figcaption></figure>
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 titis 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*<sub>ST</sub> (Weir & Cockerham, 1984) computed using the StAMPP R package (Pembleton et al., 2013). Average *F*<sub>ST</sub> 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*<sub>ST</sub>-outlier based method (using Bayesian 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.

First, we ran Bayesian (with default parameter options) for each pair of populations (Barcelona, Montpellier and Warsaw) separately to detect SNPs with outlier values of *F*<sub>ST</sub>. As recommended by 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).

## 2.4 | 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.

## 2.5 | 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* (*ggallushausensembl*) 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.

# 3 | RESULTS

## 3.1 | Small genetic and epigenetic average differentiation between urban and forest populations

### 3.1.1 | 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 ex-
plained only a small fraction (i.e. less than 2%) of the total variance
(R² = 0.018, Figure 2a, Table S2). All three variables were signifi-
cant (location: *p* = 0.001; habitat: *p* = 0.004; sex: *p* = 0.001) re-
vealing small genetic structuration between groups. Partial RDA
revealed that the net variation explained by habitat (R² = 0.004,
*p* = 0.001) was inferior to the net variation explained by location
(R² = 0.012, *p* = 0.001) but higher than sex (R² = 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 ef-
fects of other variables remained significant and of similar mag-
nitude (Table S3).

Genome-wide differentiation between populations was relatively
low on average (mean F<sub>ST</sub> = 0.019), in the order of 1%-2% between
habitats for each location (F<sub>ST Barcelona(forest-urban)</sub> = 0.018 &plusmn; 0.001;
F<sub>ST Montpellier(forest-urban)</sub> = 0.012 &plusmn; 0.001; F<sub>ST Warsaw(forest-urban)</sub> = 0.01
8 &plusmn; 0.001, Table 1), suggesting relatively high gene flow and lim-
ited genetic drift among populations. Mean F<sub>ST</sub> on autosomes was
lower (mean F<sub>ST</sub> = 0.014, Table S4) than on the Z chromosome (mean
F<sub>ST</sub> = 0.022, Table S5).

3.1.2 | 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 vari-
ance (R² = 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² = 0.003, *p* = 0.001; sex: R² = 0.003, *p* = 0.001; habi-
tat: R² = 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,
*p* = 0.044), with Montpellier individuals showing lower methylation
levels than Warsaw ones (Figure 3, Tukey test: *p* = 0.04) and no other
difference between pairs of cities. Also, no significant effect of sex
(*p* = 0.263) or habitat (*p* = 0.478) was found, suggesting that urban-
ization did not have an important overall effect on global methylation
levels. For the Z chromosome, we found a strong difference between
sexes, with homogametic males showing 2.98% more methylated Z
than heterogametic females (Figure 3; *p* = 1.45 × 10<sup>-15</sup>), while no sig-
nificant difference between location (*p* = 0.577) or habitat (*p* = 0.915)
was found.

3.2 | Nonparallel yet strong genomic footprints of divergent selection between urban and forest populations and evidence for polygenic adaptation

First, Bayesian 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

<figure><img src="image_3.png" /><figcaption>FIGURE 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</figcaption></figure>
<table><caption>TABLE 1 $F_{ST}$ estimation between pairs of subpopulations</caption><thead><tr><th></th><th>Montpellier forest</th><th>Montpellier-urban</th><th>Warsaw-forest</th><th>Warsaw-urban</th><th>Barcelona-forest</th></tr></thead><tbody><tr><td>Montpellier-urban</td><td>0.012 (0.011-0.013)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Warsaw-forest</td><td>0.009 (0.009-0.010)</td><td>0.018 (0.017-0.019)</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Warsaw-urban</td><td>0.018 (0.017-0.018)</td><td>0.026 (0.025-0.027)</td><td>0.018 (0.017-0.019)</td><td>-</td><td>-</td></tr><tr><td>Barcelona-forest</td><td>0.009 (0.009-0.010)</td><td>0.019 (0.017-0.018)</td><td>0.009 (0.009-0.010)</td><td>0.018 (0.017-0.019)</td><td>-</td></tr><tr><td>Barcelona-urban</td><td>0.025 (0.024-0.027)</td><td>0.034 (0.033-0.035)</td><td>0.025 (0.025-0.026)</td><td>0.034 (0.033-0.034)</td><td>0.018 (0.017-0.019)</td></tr></tbody></table>

Note: 95% confidence intervals in brackets were computed using StAMPP package with 1000 bootstrap.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 Distributions of mean per cent of methylation on autosomes and on the Z chromosome, compared between habitats (forest vs. urban), sexes and locations</figcaption></figure>

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.

GO analyses revealed the existence of overrepresented ontologies (*p* < 0.05 and at least 3 genes per GO). Among the most promising GO terms we found functions related to the nervous system (GO:0048846, axon extension involved in axon guidance; GO:0035418 protein localization to synapse; GO:0021987 cerebral cortex development; GO:0007274, neuromuscular synaptic transmission), the blood system (GO:0045777 positive regulation of blood pressure; GO:0042311, vasodilatation), hormonal response (GO:0071277, cellular response to oestrogen stimuli) and stress response (GO:0033555, multicellular response to stress), revealing functions potentially involved in adaptation to urban habitats. Detailed GO results are presented in Table S9 and Figure S3.

## 3.3 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 ($\chi^2$ = 11.25, *p* < 0.001), while no difference was found for Montpellier ($\chi^2$ = 0.941, *p* = 0.332) or for Warsaw ($\chi^2$ = 0.011, *p* = 0.917). DMRs were distributed across all the 32 chromosomes as well as on 37 unplaced scaffolds. 203 of the 224 different DMRs (91%) overlapped genes or 5 kb flanking regions. Among these 203 DMRs, 126 (62%) were directly located in gene bodies (*n* = 53, 26.1%), 48 (23.6%) in promotor (3 kb upstream to 300 bp downstream of genes, excluding TSS) or TSS sequences (300 bp upstream to 50 bp downstream of genes), and 26 (12.8%) in both gene body and promotor/TSS.

Following the procedure previously described, GO analyses on the pooled genes list revealed an overrepresentation of modules associated with the nervous system (GO:2000300, regulation of synaptic vesicle exocytosis; GO:0050804 modulation of synaptic transmission), immunity (GO:005728, negative regulation of inflammatory response; GO:0050852, T cell receptor signalling pathway), metabolic activity (GO:006816 calcium ion transport; GO:0055072, iron ion homeostasis, GO:0043087, regulation of GTPase activity), behaviour (GO:0007626, locomotory behaviour) and endocrine processes (GO:0044060: regulation of endocrine process). All enriched GO are presented in 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 58 for Barcelona, 81 for Montpellier and 99 for Warsaw. Warsaw presented significantly more hyper than hypomethylated (in females compared to males) DMRs ($\chi^2$ = 5.878, *p* = 0.015), but it was not
<figure><img src="image_5.png" /><figcaption>FIGURE 4 Manhattan plot of mean FST between urban and forest populations along the Great tit genome for (a) Barcelona, (b) Montpellier and (c) Warsaw. Dark orange (a), purple (b) and light orange points (c) represent significant outlier SNPs identified by the FST outlier test Bayesian scan for each population pair, given with their associated genes in 5 kb. Green points represent outliers found with the multivariate RDA approach. A few SNPs were identified by both the Bayesian scan and the RDA methods and signalled as a green point circled with the colour used for the considered pair</figcaption></figure>

the case for Barcelona ($\chi^2 = 0, p = 1$) nor Montpellier ($\chi^2 = 1.25, p = 0.264$). DMRs were distributed on 29 chromosomes and 35 unplaced scaffolds. On a total of 206 DMRs, 181 (57.3%) were on genes or in a 5 kb upstream/downstream region around genes. GO analyses revealed enrichment of genes involved in development, growth and morphogenesis, among others (see detailed enriched GO Table 3, 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: $\chi^2$-squared = 3.885, *p* = 0.049). 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).

4 | DISCUSSION

The urban sprawl is a worldwide phenomenon deeply affecting the environment and thus requiring fast adaptive responses in city dwellers. While a large body of literature already describes a myriad of phenotypic shifts in urban populations of numerous species (Chamberlain et al., 2009; Lowry et al., 2013; Sepp et al., 2018), the molecular bases of these shifts and their evolutionary implications remain yet to be documented and understood. This study uses genomic and epigenomic analyses to decipher the potential molecular bases implicated in phenotypic shifts and adaptation in several urban populations of a passerine bird, the Great tit. Note that this species shows largely parallel phenotypic shifts across its range in terms of morphology and life history (Biard et al., 2017; Thompson et al., 2021). In our study, genomic analyses revealed weak yet significant average differentiation between urban and forest populations, suggesting ongoing gene flow and limited drift in urban populations. These analyses also identified a limited number of loci putatively under strong selection, nonrepeated between pairs, and numerous loci supposedly under weaker selection, compatible with a polygenic model of evolution. On the epigenomic side, while we found weak average differentiation of the methylome between urban and forest birds, suggesting an absence of genome-wide epigenetic
<figure><img src="image_6.png" /><figcaption>FIGURE 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</figcaption></figure>

deregulations, we identified several differentially methylated regions between urban and forest birds, mostly nonrepeated between pairs. Genes associated with either genomic footprints of divergent selection or differentially methylated regions had relatively similar functions, related to the nervous system, metabolism, immunity and behaviour, that have been repeatedly convicted in other studies (Riyahi et al., 2017). Hence, by identifying nonrepeated genetic

and epigenetic responses among replicated forest-urban population pairs, our findings support the hypothesis of mostly nonparallel rapid de novo adaptation to similar environments via both genetic and epigenetic mechanisms. Our results are in line with accumulating evidence that polygenic adaptation and epigenetic reprogramming may be involved in quick phenotypic shifts in response to rapidly emerging constraints such as urbanization.
<table><caption>TABLE 2 Significantly enriched GO associated with genes overlapping 5 kb windows around forest and urban habitats</caption><thead><tr><th>GO.ID</th><th>Term</th><th>Annotated</th><th>Significant</th><th>Expected</th><th>Weight</th><th>Barcelona</th><th>Montpellier</th><th>Warsaw</th></tr></thead><tbody><tr><td>GO:0050804</td><td>Modulation of chemical synaptic transmission</td><td>190</td><td>10</td><td>2.92</td><td>0.016</td><td>CACNG3, CPLX2, OPHN1, PPP1R9B, PRKCB</td><td>CDH11</td><td>DRD3, EPHB1, NAPB, SNAPIN</td></tr><tr><td>GO:0042981</td><td>Regulation of apoptotic process</td><td>619</td><td>8</td><td>9.52</td><td>0.004</td><td>DIP2A, HERPUD1, NR4A1, SOD1</td><td></td><td>ANP32E, RSL1D1, UACA</td></tr><tr><td>GO:0007626</td><td>Locomotory behavior</td><td>122</td><td>6</td><td>1.88</td><td>0.005</td><td>ADCY6, ELAVL4, SOD1</td><td></td><td>CDH23, DRD3, SEZ6L</td></tr><tr><td>GO:0006816</td><td>Calcium ion transport</td><td>173</td><td>6</td><td>2.66</td><td>0.024</td><td>CACNB3, CACNG3, PRKCB</td><td></td><td>CDH23, DRD3, STIM1</td></tr><tr><td>GO:0050728</td><td>Negative regulation of inflammatory response</td><td>61</td><td>5</td><td>0.94</td><td>0.011</td><td>PSMB4, SOD1</td><td>CD276, WFDC1</td><td>PBK</td></tr><tr><td>GO:0046677</td><td>Response to antibiotic</td><td>91</td><td>5</td><td>1.4</td><td>0.012</td><td>ADCY6, SOD1</td><td>RSRC1</td><td>AOC1, DRD3</td></tr><tr><td>GO:0043087</td><td>Regulation of GTPase activity</td><td>237</td><td>5</td><td>3.65</td><td>0.025</td><td>GAPVD1, OPHN1, RANBP3L, SOD1</td><td>PROM2</td><td>PROM2</td></tr><tr><td>GO:2000300</td><td>Regulation of synaptic vesicle exocytosis</td><td>40</td><td>4</td><td>0.62</td><td>0.017</td><td>CPLX2</td><td>PRKCB</td><td>NAPB, SNAPIN</td></tr><tr><td>GO:0008217</td><td>Regulation of blood pressure</td><td>74</td><td>4</td><td>1.14</td><td>0.04</td><td>SOD1</td><td>POMC</td><td>ADRA1B, DRD3</td></tr><tr><td>GO:0044060</td><td>Regulation of endocrine process</td><td>22</td><td>3</td><td>0.34</td><td>0.012</td><td>NR3C1</td><td>POMC</td><td>RAB11FIP3</td></tr><tr><td>GO:0055072</td><td>Iron ion homeostasis</td><td>42</td><td>3</td><td>0.65</td><td>0.014</td><td>NEO1, SOD1</td><td></td><td>HMOX2</td></tr><tr><td>GO:0071559</td><td>Response to transforming growth factor beta</td><td>113</td><td>3</td><td>1.74</td><td>0.031</td><td>NR3C1, ZBTB7A, ZFHX3</td><td></td><td></td></tr><tr><td>GO:0050852</td><td>T cell receptor signaling pathway</td><td>52</td><td>3</td><td>0.8</td><td>0.046</td><td>CACNB3, SKAP1</td><td>CD276</td><td></td></tr><tr><td>GO:0007200</td><td>Phospholipase C-activating G protein-coupled receptor signaling pathway</td><td>40</td><td>3</td><td>0.62</td><td>0.048</td><td>ADRA1B</td><td>P2RY1</td><td>DRD3</td></tr></tbody></table>
<table><caption>TABLE 3 Significantly enriched GO associated with genes overlapping 5 kb windows around DMRs between sexes</caption><thead><tr><th>GO.ID</th><th>Term</th><th>Annotated</th><th>Significant</th><th>Expected</th><th>Weight</th><th>Barcelona</th><th>Montpellier</th><th>Warsaw</th></tr></thead><tbody><tr><td>GO:1903706</td><td>Regulation of hemopoiesis</td><td>172</td><td>8</td><td>2.28</td><td>0.050</td><td>KAT7, RARA, RUNX1, TJP2</td><td>CARD11, KAT7, PTN, RASSF2</td><td></td></tr><tr><td>GO:1902930</td><td>Regulation of alcohol biosynthetic process</td><td>20</td><td>3</td><td>0.27</td><td>0.039</td><td>SCAP</td><td>SCAP, SOD1, VDR</td><td>FGFR1</td></tr><tr><td>GO:1902652</td><td>Secondary alcohol metabolic process</td><td>53</td><td>4</td><td>0.7</td><td>0.013</td><td>LSS, SCAP</td><td>LSS, SCAP, SOD1</td><td>FGFR1</td></tr><tr><td>GO:0060173</td><td>Limb development</td><td>123</td><td>4</td><td>1.63</td><td>0.028</td><td>RARA</td><td>PTN, TBX3</td><td>FGFR1, TBX3</td></tr><tr><td>GO:0060117</td><td>Auditory receptor cell development</td><td>13</td><td>3</td><td>0.17</td><td>0.026</td><td>WHRN</td><td>SOD1</td><td>FGFR1</td></tr><tr><td>GO:0046887</td><td>Positive regulation of hormone secretion</td><td>46</td><td>3</td><td>0.61</td><td>0.016</td><td></td><td></td><td>BLK, FGFR1, NMB</td></tr><tr><td>GO:0046849</td><td>Bone remodeling</td><td>46</td><td>3</td><td>0.61</td><td>0.023</td><td>TPH1</td><td>PTN, RASSF2</td><td></td></tr><tr><td>GO:0043087</td><td>Regulation of GTPase activity</td><td>237</td><td>6</td><td>3.15</td><td>0.006</td><td>MLST8</td><td>PROM2, SOD1</td><td>LRRK2, PROM2, RGS14, STARD8</td></tr><tr><td>GO:0042981</td><td>Regulation of apoptotic process</td><td>619</td><td>19</td><td>8.22</td><td>0.002</td><td>DIP2A, RARA, UBE2Z</td><td>BAG5, CARD11, CIDEc, DIP2A, GRK5, RASSF2, SOD1, TBX3</td><td>ANP32E, BLK, FGFR1, HERPUD1, KDM2B, LRRK2, PDIA3, PHLDA1, TBX3</td></tr><tr><td>GO:0035264</td><td>Multicellular organism growth</td><td>93</td><td>4</td><td>1.23</td><td>0.033</td><td>RARA, SELENOM</td><td>ANKRD11, SOD1</td><td>ANKRD11, SELENOM</td></tr><tr><td>GO:0032436</td><td>Positive regulation of proteasomal ubiquitin-dependent protein catabolic process</td><td>46</td><td>4</td><td>0.61</td><td>0.018</td><td></td><td>CSNK1A1</td><td>CSNK1A1, HERPUD1, LRRK2, STUB1</td></tr><tr><td>GO:0032091</td><td>Negative regulation of protein binding</td><td>50</td><td>3</td><td>0.66</td><td>0.028</td><td></td><td>ZFPM1</td><td>LRRK2, STUB1</td></tr><tr><td>GO:0031647</td><td>Regulation of protein stability</td><td>152</td><td>5</td><td>2.02</td><td>0.048</td><td>LSS</td><td>BAG5, LSS, RASSF2</td><td>LRRK2, STUB1</td></tr><tr><td>GO:0009790</td><td>Embryo development</td><td>625</td><td>15</td><td>8.3</td><td>0.049</td><td>MDFI, PCDH8, RARA, WHRN</td><td>ANKRD11, HOXB9, PTN, SOD1, TBX1, TBX3, ZFPM1</td><td>ANKRD11, ASCL2, DNMT3A, FGFR1, KDM2B, PCDH8, TBX3</td></tr><tr><td>GO:0009267</td><td>Cellular response to starvation</td><td>85</td><td>4</td><td>1.13</td><td>0.044</td><td>PIK3C2B</td><td>XPR1</td><td>LRRK2, WDR24</td></tr><tr><td>GO:0007218</td><td>Neuropeptide signaling pathway</td><td>46</td><td>3</td><td>0.61</td><td>0.023</td><td></td><td></td><td>GALR2, NMB, POMC</td></tr><tr><td>GO:0003151</td><td>Outflow tract morphogenesis</td><td>46</td><td>3</td><td>0.61</td><td>0.046</td><td>RARA</td><td>TBX3, ZFPM1</td><td>TBX3</td></tr><tr><td>GO:0001666</td><td>Response to hypoxia</td><td>91</td><td>4</td><td>1.21</td><td>0.023</td><td>NPEP5, SCAP</td><td>NPEP5, SCAP</td><td>ASCL2, LIMD1</td></tr><tr><td>GO:0000122</td><td>Negative regulation of transcription by RNA polymerase II</td><td>454</td><td>12</td><td>6.03</td><td>0.017</td><td>RARA</td><td>CBFA2T2, SCAF4, TBX1, TBX3, VDR, ZFPM1</td><td>CSNK1A1, HERPUD1, LRRK2, STUB1</td></tr></tbody></table>

Note: Bold genes represent when a DMR overlapped a gene boy and underlined genes, when DMR overlapped TSS or promoter regions.
<figure><img src="image_7.png" /><figcaption>FIGURE 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</figcaption></figure>

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
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 (Sprau & Dingemanse, 2017). 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 (Salmón et al., 2021; 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 promoter or TSS regions. This latter result contrasts with Watson and collaborators (Watson et al., 2021) who recently found an underrepresentation of DMS in both gene body and regulatory regions in urban great tits from Malmö (Sweden). Across the three cities, 62% of DMR were directly localized in gene bodies and/or TSS or promoter regions, suggesting a potential role in gene expression modulation. Direction of methylation patterns did not follow any consistent pattern (no over-representation of hypo- or hypermethylated DMR in urban birds, Figure 5), in line with Watson and collaborators' analyses on blood sample. However, birds in Barcelona presented significantly more hypomethylated DMR than hypermethylated ones.

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. Additionally, cities are different from each other because of a wide array of climatic, cultural, historical and socioeconomic factors (Szulkin, Garroway, et al., 2020; Szulkin, Munshi-South, et al., 2020). 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 black-bird (Garroway & Sheldon, 2013), in the black swan (Van Dongen et al., 2015) and in the burrowing owl (Mueller et al., 2020)). In this study, while urban great tits show higher levels of aggressiveness in at least two of the cities (A. E. Caizergues, A. Grégoire, R. Choquet, S. Perret, & A. Charmantier, unpublished data; Riyahi et al., 2017) we found no DMR associated with these two genes in either of the three forest-city pairs. However, we found a significant urban-related 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) blood methylation patterns might not reflect the ones in other tissues, (2) methylation shifts are known to change gene expression in the tissue sampled, but little is known on their effects on other tissues

and (3) DNA methylation patterns can undergo seasonal variation (Viitaniemi et al., 2019) which is especially true for blood because of rapid cellular turnover. Hence, analyses on multiple tissues and life-stages replicated in multiple pairs of urban and forest populations are required, to draw a broader view of the impact of urbanization on global methylation patterns and to understand replicated parallel occurrence across cities. However, tissue-specific and age-specific analyses in multiple individuals across several pairs of urban and forest environments pose major technical, budget and ethical limitations and should be coordinated very carefully. Additionally, specific drivers of shifts in methylation remain to be disentangled to understand which environmental factors are responsible for which change in methylation. To do so, experimental settings manipulating environmental factors such as performed by Mäkinen et al. (2019) would be particularly useful. More integratively, information on how shifts in methylation patterns affect phenotype, fitness and adaptation often remain elusive (Sepers et al., 2019). To our knowledge, a limited number of studies attempted to link methylation and expression levels in natural population contexts (Derks et al., 2016; Laine et al., 2016), and even fewer in urban habitats (but see e.g. McNew et al., 2017; Watson et al., 2021). Hence, future work might need to tackle the question of the origin and adaptive significance of these variations in a controlled framework.

Finally, the consideration of sample size in terms of number of individuals sampled per population as well as number of replicated populations is undoubtedly of a major importance. Although this study is in the upper range of sample sizes found in previous studies on parallelism, we may still have limited power to detect parallel marks of evolution if methylation shifts at these marks are relatively low. To date, no power analysis has been developed specifically for (epi)genomic parallelism and such methodological development could be extremely useful to evaluate the confidence in our studies. Future studies should aim at increasing individual sample size per population and conducting broader geographic scale analyses with more replicated population pairs (see e.g. Salmón et al. (2021) across nine European cities) to increase the power to detect parallelism in epigenomic responses to urbanization.

This study identified both genomic footprints of selection and differentially methylated regions between urban and forest populations, suggesting that both genetic and epigenetic processes could play a role in rapid adaptation to urban habitat. To our knowledge, our study is the first to use replicated pairs of city and forest populations to study modifications of methylation in urban habitats. This study design revealed limited evidence for parallelism between cities both at the genetic and epigenetic levels, suggesting that cities might not present exactly similar environmental conditions or that different genetic and epigenetic pathways are involved in adaptation to urban environmental conditions, although possibly associated with similar biological functions. This study finally highlights the need to unravel both environmental origins and evolutionary implications of methylation shifts, in order to understand to which extent environmentally induced methylation can contribute to adaptation.
## ACKNOWLEDGEMENTS

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.

## CONFLICT OF INTEREST

We declare we have no competing interests.

## DATA AVAILABILITY STATEMENT

DNA sequences: RAD-seq GenBank accessions SRR17125662-SRR17125721, RRBS-seq GenBank accessions SRR17145241-SRR17145300; NCBI project PRJNA786007. Scripts: Github repository https://github.com/AudeCaizergues/Epigenetics_and_the.city

## ORCID

Aude E. Caizergues https://orcid.org/0000-0003-4467-3912

Jeremy Le Luyer https://orcid.org/0000-0001-9409-3196

Arnaud Grégoire https://orcid.org/0000-0001-6103-355X

Marta Szulkin https://orcid.org/0000-0002-7355-5846

Juan-Carlos Senar https://orcid.org/0000-0001-9955-3892

Anne Charmantier https://orcid.org/0000-0002-0691-2647

Charles Perrier https://orcid.org/0000-0001-5820-9374

## REFERENCES

Alberti, M., Correa, C., Marzluff, J. M., Hendry, A. P., Palkovacs, E. P., Gotanda, K. M., Hunt, V. M., Apgar, T. M., & Zhou, Y. (2017). Global urban signatures of phenotypic change in animal and plant populations. *Proceedings of the National Academy of Sciences of the United States of America*, 114(34), 8951–8956. https://doi.org/10.1073/pnas.1606034114

Alexa, A., & Rahnenführer, J. (2009). *Gene set enrichment analysis with topGO*. Bioconductor Improv, 27.

Baird, N. A., Etter, P. D., Atwood, T. S., Currey, M. C., Shiver, A. L., Lewis, Z. A., Selker, E. U., Cresko, W. A., & Johnson, E. A. (2008). Rapid SNP discovery and genetic mapping using sequenced RAD markers. *PLoS One*, 3(10), e3376. https://doi.org/10.1371/journal.pone.0003376

Biard, C., Brischoux, F., Meillère, A., Michaud, B., Nivière, M., Ruault, S., Vaugoyeau, M., & Angelier, F. (2017). Growing in cities: An urban penalty for wild birds? A study of phenotypic differences between urban and rural great tit chicks (*Parus major*). *Frontiers in Ecology and Evolution*, 5, 79. https://doi.org/10.3389/fevo.2017.00079

Boyle, E. A., Li, Y. I., & Pritchard, J. K. (2017). An expanded view of complex traits: From polygenic to omnigenic. *Cell*, 169(7), 1177–1186. https://doi.org/10.1016/j.cell.2017.05.038

Brondizio, E. S., Settele, J., & Díaz, S. (2019). *IPBES. 2019 Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services*. Bonn, Germany.

Buitenhuis, B., Hedegaard, J., Janss, L., & Sørensen, P. (2009). Differentially expressed genes for aggressive pecking behaviour in laying hens. *BMC Genomics*, 10(1), 544. https://doi.org/10.1186/1471-2164-10-544

Caizergues, A. E., Charmantier, A., Lambrechts, M. M., Perret, S., Demeyrier, V., Lucas, A., & Grégoire, A. (2021). An avian urban morphology: How the city environment shapes great tit morphology at different life stages. *Urban Ecosystems*, 24(5), 929–941. https://doi.org/10.1007/s11252-020-01077-0

Caizergues, A. E., Grégoire, A., & Charmantier, A. (2018). Urban versus forest ecotypes are not explained by divergent reproductive selection. *Proceedings of the Royal Society B: Biological Sciences*, 285(1882), 20180261. https://doi.org/10.1098/rspb.2018.0261

Cárdenas Rodríguez, M., Dupont-Courtade, L., & Oueslati, W. (2016). Air pollution and urban structure linkages: Evidence from European cities. *Renewable and Sustainable Energy Reviews*, 53, 1–9. https://doi.org/10.1016/j.rser.2015.07.190

Chamberlain, D., Hatchwell, B., & Gaston, K. J. (2009). Avian productivity in urban landscapes: A review and meta-analysis. *Ibis*, 151, 1–18. https://doi.org/10.1111/j.1474-919X.2008.00899.x

Chen, S., Zhou, Y., Chen, Y., & Gu, J. (2018). fastp: An ultra-fast all-in-one FASTQ preprocessor. *Bioinformatics*, 34(17), i884–i890. https://doi.org/10.1093/bioinformatics/bty560

Corsini, M., Schöll, E. M., Di Lecce, I., Chatelain, M., Dubiec, A., & Szulkin, M. (2020). Growing in the city: Urban evolutionary ecology of avian growth rates. *Evolutionary Applications*, 14(1), 69–84. https://doi.org/10.1111/eva.13081

Crews, D., Gore, A. C., Hsu, T. S., Dangleben, N. L., Spinetta, M., Schallert, T., Anway, M. D., & Skinner, M. K. (2007). Transgenerational epigenetic imprints on mate preference. *Proceedings of the National Academy of Sciences of the United States of America*, 104(14), 5942–5946. https://doi.org/10.1073/pnas.0610410104

Danchin, E., Charmantier, A., Champagne, F. A., Mesoudi, A., Pujol, B., & Blanchet, S. (2011). Beyond DNA: Integrating inclusive inheritance into an extended theory of evolution. *Nature Reviews Genetics*, 12, 475–486. https://doi.org/10.1038/nrg3028 https://doi.org/10.1038/nrg302810.1038/nrg3028

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S. T., McVean, G., & Durbin, R. (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156–2158. https://doi.org/10.1093/bioinformatics/btr330

Derks, M. F. L., Schachtschneider, K. M., Madsen, O., Schijlen, E., Verhoeven, K. J. F., & van Oers, K. (2016). Gene and transposable element methylation in great tit (*Parus major*) brain and blood. *BMC Genomics*, 17(1), 332. https://doi.org/10.1186/s12864-016-2653-y

Ding, R., Jin, Y., Liu, X., Ye, H., Zhu, Z., Zhang, Y., Wang, T., & Xu, Y. (2017). Dose- and time- effect responses of DNA methylation and histone H3K9 acetylation changes induced by traffic-related air pollution. *Scientific Reports*, 7(1), 1–9. https://doi.org/10.1038/srep43737

Donihue, C. M., & Lambert, M. R. (2015). Adaptive evolution in urban ecosystems. *Ambio*, 44(3), 194–203. https://doi.org/10.1007/s13280-0-014-0547-2

Dubin, M. J., Zhang, P., Meng, D., Remigereau, M.-S., Osborne, E. J., Paolo Casale, F., Drewe, P., Kahles, A., Jean, G., Vilhjálmsson, B., Jagoda, J., Irez, S., Voronin, V., Song, Q., Long, Q., Rätsch, G., Stegle, O., Clark, R. M., & Nordborg, M. (2015). DNA methylation in Arabidopsis has a genetic basis and shows evidence of local adaptation. *eLife*, 4(May), e05255. https://doi.org/10.7554/ELIFE.05255

Duncan, E. J., Gluckman, P. D., & Dearden, P. K. (2014). Epigenetics, plasticity, and evolution: How do we link epigenetic change to phenotype? *Journal of Experimental Zoology Part B: Molecular and Developmental Evolution*, 322(4), 208–220. https://doi.org/10.1002/jez.b.22571

Foll, M., & Gaggiotti, O. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. *Genetics*, 180(2), 977–993. https://doi.org/10.1534/genetics.108.092221

Forester, B. R., Jones, M. R., Joost, S., Landguth, E. L., & Lasky, J. R. (2016). Detecting spatial genetic signatures of local adaptation
in heterogeneous landscapes. *Molecular Ecology*, 25(1), 104–120.
https://doi.org/10.1111/mec.13476

Foust, C. M., Preite, V., Schrey, A. W., Alvarez, M., Robertson, M. H., Verhoeven, K. J. F., & Richards, C. L. (2016). Genetic and epigenetic differences associated with environmental gradients in replicate populations of two salt marsh perennials. *Molecular Ecology*, 25(8), 1639–1652. https://doi.org/10.1111/mec.13522

Garroway, C. J., & Sheldon, B. C. (2013). Urban behavioural adaptation. *Molecular Ecology*, 22(13), 3430–3432. https://doi.org/10.1111/mec.12351

Grimm, N. B., Faeth, S. H., Golubiewski, N. E., Redman, C. L., Wu, J., Bai, X., & Briggs, J. M. (2008). Global change and the ecology of cities. *Science*, 319(5864), 756–760. https://doi.org/10.1126/science.1150195

Head, J. A. (2014). Patterns of DNA methylation in animals: An ecotoxicological perspective. *Integrative and Comparative Biology*, 54(1), 77–86. https://doi.org/10.1093/icb/icu025

Husby, A. (2020). On the use of blood samples for measuring DNA methylation in ecological epigenetic studies. *Integrative and Comparative Biology*, 60(6), 1558–1566. https://doi.org/10.1093/icb/icaa123

Jaenisch, R., & Bird, A. (2003). Signals, epigenetic regulation of gene expression: How the genome integrates intrinsic and environmental. *Nature Genetics*, 33, 245–254.

Janowitz Koch, I., Clark, M. M., Thompson, M. J., Deere-Machemer, K. A., Wang, J., Duarte, L., Gnanadesikan, G. E., McCoy, E. L., Rubbi, L., Stahler, D. R., Pellegrini, M., Ostrander, E. A., Wayne, R. K., Sinsheimer, J. S., & von Holdt, B. M. (2016). The concerted impact of domestication and transposon insertions on methylation patterns between dogs and grey wolves. *Molecular Ecology*, 25(8), 1838–1855. https://doi.org/10.1111/mec.13480

Kim, Y. H., Peregrine, J., & Arnold, A. P. (2006). The distribution of expression of doublecortin (DCX) mRNA and protein in the zebra finch brain. *Brain Research*, 1106(1), 189–196. https://doi.org/10.1016/j.brainres.2006.05.080

Krueger, F., & Andrews, S. R. (2011). Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. *Bioinformatics*, 27(11), 1571–1572. https://doi.org/10.1093/bioinformatics/btr167

Kvist, L., Martens, J., Higuchi, H., Nazarenko, A. A., Valchuk, O. P., & Orell, M. (2003). Evolution and genetic structure of the great tit (*Parus major*) complex. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, 270(1523), 1447–1454. https://doi.org/10.1098/rspb.2002.2321

Laine, V. N., Gossmann, T. I., Schachtschneider, K. M., Garroway, C. J., Madsen, O., Verhoeven, K. J. F., de Jager, V., Megens, H.-J., Warren, W. C., Minx, P., Crooijmans, R. P. M. A., Corcoran, P., Sheldon, B. C., Slate, J., Zeng, K., van Oers, K., Visser, M. E., & Groenen, M. A. M. (2016). Evolutionary signals of selection on cognition from the great tit genome and methylome. *Nature Communications*, 7(1), 1–9. https://doi.org/10.1038/ncomms10474

Lande, R. (1979). Quantitative genetic analysis of multivariate evolution, applied to brain: Body size allometry. *Evolution*, 33(1), 402. https://doi.org/10.2307/2407630

Li, H., & Durbin, R. (2010). Fast and accurate long-read alignment with Burrows-Wheeler transform. *Bioinformatics*, 26(5), 589–595. https://doi.org/10.1093/bioinformatics/btp698

Li, Z., Zheng, M., Abdalla, B. A., Zhang, Z., Xu, Z., Ye, Q., Xu, H., Luo, W., Nie, Q., & Zhang, X. (2016). Genome-wide association study of aggressive behaviour in chicken. *Scientific Reports*, 6(1), 1–11. https://doi.org/10.1038/srep30981

Liu, Q. A. (2013). The impact of climate change on plant epigenomes. *Trends in Genetics*, 29(9), 503–505. https://doi.org/10.1016/j.tig.2013.06.004

Losos, J. B. (2011). Convergence, adaptation, and constraint. *Evolution*, 65(7), 1827–1840. https://doi.org/10.1111/j.1558-5646.2011.01289.x

Lowry, H., Lill, A., & Wong, B. B. M. (2013). Behavioural responses of wildlife to urban environments. *Biological Reviews*, 88(3), 537–549.
https://doi.org/10.1111/brv.12012

Mäkinen, H., van Oers, K., Eeva, T., Laine, V. N., & Ruuskanen, S. (2019, November 26). The effect of experimental lead pollution on DNA methylation in a wild bird population. *BioRxiv*, p. 851998. bioRxiv.
https://doi.org/10.1101/851998

Malhi, Y. (2017). The concept of the Anthropocene. *Annual Review of Environment and Resources*, 42(1), 77–104. https://doi.org/10.1146/annurev-environment-102016-060854

McNew, S. M., Beck, D., Sadler-Riggleman, I., Knutie, S. A., Koop, J. A. H., Clayton, D. H., & Skinner, M. K. (2017). Epigenetic variation between urban and rural populations of Darwin's finches. *BMC Evolutionary Biology*, 17(1), 183. https://doi.org/10.1186/s12862-017-1025-9

McNulty, S. E., Barrett, R. M., Vogel-Ciernia, A., Malvaez, M., Hernandez, N., Davatolhagh, M. F., & Wood, M. A. (2012). Differential roles for Nr4a1 and Nr4a2 in object location vs. object recognition long-term memory. *Learning and Memory*, 19(12), 588–592. https://doi.org/10.1101/lm.026385.112

Meaney, M. J., & Szyf, M. (2005). Environmental programming of stress responses through DNA methylation: Life at the interface between a dynamic environment and a fixed genome. *Dialogues in Clinical Neuroscience*, 7(2), 103–123. https://doi.org/10.31887/dcns.2005.7.2/mmeaney

Mueller, J. C., Carrete, M., Boerno, S., Kuhl, H., Tella, J. L., & Kempenaers, B. (2020). Genes acting in synapses and neuron projections are early targets of selection during urban colonization. *Molecular Ecology*, 29(18), 3403–3412. https://doi.org/10.1111/mec.15451

Mueller, J. C., Kuhl, H., Boerno, S., Tella, J. L., Carrete, M., & Kempenaers, B. (2018). Evolution of genomic variation in the burrowing owl in response to recent colonization of urban areas. *Proceedings of the Royal Society B: Biological Sciences*, 285(1878), 20180206. https://doi.org/10.1098/rspb.2018.0206

Mueller, J. C., Partecke, J., Hatchwell, B. J., Gaston, K. J., & Evans, K. L. (2013). Candidate gene polymorphisms for behavioural adaptations during urbanization in blackbirds. *Molecular Ecology*, 22(13), 3629–3637. https://doi.org/10.1111/mec.12288

Pembbleton, L., Cogan, N., & Forster, J. (2013). StAMPP: An R package for calculation of genetic differentiation and structure of mixed-ploidy level populations. *Molecular Ecology Resources*, 13, 946–952.
https://doi.org/10.1111/1755-0998.12129

Perrier, C., Caizergues, A. E., & Charmantier, A. (2020). Adaptation genomics in urban environments. In M. Szulkin, J. Munshi-South, & A. Charmantier (Eds.), *Urban evolutionary biology* (pp. 72–88). Oxford University Press.
https://doi.org/10.1093/oso/9780198836841.003.0005

Perrier, C., Lozano del Campo, A., Szulkin, M., Demeyrier, V., Gregoire, A., & Charmantier, A. (2018). Great tits and the city: Distribution of genomic diversity and gene-environment associations along an urbanization gradient. *Evolutionary Applications*, 11(5), 593–613.
https://doi.org/10.1111/eva.12580

Quinlan, A. R., & Hall, I. M. (2010). BEDTools: A flexible suite of utilities for comparing genomic features. *Bioinformatics*, 26(6), 841–842.
https://doi.org/10.1093/bioinformatics/btq033

Reid, N. M., Proestou, D. A., Clark, B. W., Warren, W. C., Colbourne, J. K., Shaw, J. R., Karchner, S. I., Hahn, M. E., Nacci, D., Oleksiak, M. F., Crawford, D. L., & Whitehead, A. (2016). The genomic landscape of rapid repeated evolutionary adaptation to toxic pollution in wild fish. *Science*, 354(6317), 1305–1308. https://doi.org/10.1126/science.aah4993

Reik, W., Dean, W., & Walter, J. (2001). Epigenetic reprogramming in mammalian development. *Science*, 293(5532), 1089–1093. https://doi.org/10.1126/science.1063443

Richards, C. L., Bossdorf, O., & Pigliucci, M. (2010). What role does heritable epigenetic variation play in phenotypic evolution? *BioScience*, 60(3), 232–237. https://doi.org/10.1525/bio.2010.60.3.9
Rivkin, L. R., Santangelo, J. S., Alberti, M., Aronson, M. F. J., de Keyzer, C. W., Diamond, S. E., Fortin, M.-J., Frazee, L. J., Gorton, A. J., Hendry, A. P., Liu, Y., Losos, J. B., Maclvor, J. S., Martin, R. A., McDonnell, M. J., Miles, L. S., Munshi-South, J., Ness, R. W., Newman, A. E. M., ... Johnson, M. T. J. (2018). A roadmap for urban evolutionary ecology. *Evolutionary Applications*, 12(3), 384–398. https://doi.org/10.1111/eva.12734

Riyahi, S., Björklund, M., Mateos-Gonzalez, F., & Senar, J. C. (2017). Personality and urbanization: Behavioural traits and DRD4 SNP830 polymorphisms in great tits in Barcelona city. *Journal of Ethology*, 35(1), 101–108. https://doi.org/10.1007/s1016 4-016-0496-2

Riyahi, S., Sánchez-Delgado, M., Calafell, F., Monk, D., & Senar, J. C. (2015). Combined epigenetic and intraspecific variation of the DRD4 and SERT genes influence novelty seeking behavior in great tit *Parus major*. *Epigenetics*, 10(6), 516–525. https://doi.org/10.1080/15592 94.2015.1046027

Robinson, M. R., Wray, N. R., & Visscher, P. M. (2014). Explaining additional genetic variation in complex traits. *Trends in Genetics*, 30(4), 124–132. https://doi.org/10.1016/j.tig.2014.02.003

Rochette, N. C., Rivera-Colón, A. G., & Catchen, J. M. (2019). Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics. *Molecular Ecology*, 28(21), 4737–4754. https://doi.org/10.1111/mec.15253

Salmón, P., Jacobs, A., Ahrén, D., Biard, C., Dingemanse, N. J., Dominoni, D. M., & Isaksson, C. (2021). Repeated genomic signatures of adaptation to urbanisation in a songbird across Europe. *Nature Communications*, 12(1), 1–14. https://doi.org/10.1038/s41467-021- 23027-w

Santangelo, J. S., Miles, L. S., Breitbart, S. T., Murray-Stoker, D., Rivkin, L. R., Johnson, M. T. J., & Ness, R. W. (2020). Urban environments as a framework to study parallel evolution. In M. Szulkin, J. Munshi-South, & A. Charmantier (Eds.), *Urban evolutionary biology* (pp. 36– 53). Oxford University Press.

Santangelo, J. S., Rivkin, L. R., & Johnson, M. T. J. (2018). The evolution of city life. *Proceedings of the Royal Society B: Biological Sciences*, 285(1884), 20181529. https://doi.org/10.1098/rspb.2018.1529

Senar, J. C., & Björklund, M. (2020). Recent spread of blue tits into the Barcelona urban environment: Morphological differences and the role of balanced dispersal. *Evolutionary Ecology*, 35(1), 83–99. https://doi.org/10.1007/s10682-020-10087-5

Senar, J. C., Garamszegi, L. Z., Tilgar, V., Biard, C., Moreno-Rueda, G., Salmón, P., Rivas, J. M., Sprau, P., Dingemanse, N. J., Charmantier, A., Demeyrier, V., Navalpotro, H., & Isaksson, C. (2017). Urban great tits (*Parus major*) show higher distress calling and pecking rates than rural birds across Europe. *Frontiers in Ecology and Evolution, 5*, 163. https://doi.org/10.3389/fevo. 2017.00163

Sepers, B., van den Heuvel, K., Lindner, M., Viitaniemi, H., Husby, A., & van Oers, K. (2019). Avian ecological epigenetics: Pitfalls and promises. *Journal of Ornithology*, 160(4), 1183–1203. https://doi. org/10.1007/s10336-019-01684-5

Sepp, T., Mcgraw, K. J., Kaasik, A., & Giraudeau, M. (2018). A review of urban impacts on avian life-history evolution: Does city living lead to slower pace of life? *Global Change Biology*, 24(4), 1452–1469. https://doi.org/10.1111/gcb.13969

Seress, G., Sándor, K., Evans, K. L., & Liker, A. (2020). Food availability limits avian reproduction in the city: An experimental study on great tits *Parus major*. *Journal of Animal Ecology*, 1365–2656, 13211. https://doi.org/10.1111/1365-2656.13211

Shochat, E., Warren, P. S., Faeth, S. H., McIntyre, N. E., & Hope, D. (2006). From patterns to emerging processes in mechanistic urban ecology. *Trends in Ecology & Evolution*, 21(4), 186–191. https://doi. org/10.1016/j.TREE.2005.11.019

Sih, A., & Del Giudice, M. (2012). Linking behavioural syndromes and cognition: A behavioural ecology perspective. *Philosophical

*Transactions of the Royal Society B: Biological Sciences*, 367(1603), 2762–2772. https://doi.org/10.1098/rstb.2012.0216

Sol, D., Lapiedra, O., & González-Lagos, C. (2013). Behavioural adjustments for a life in the city. *Animal Behaviour*, 85(5), 1101–1112. https://doi.org/10.1016/j.anbehav.2013.01.023

Sprau, P., & Dingemanse, N. J. (2017). An approach to distinguish between plasticity and non-random distributions of behavioral types along urban gradients in a wild passerine bird. *Frontiers in Ecology and Evolution, 5*, 92. https://doi.org/10.3389/fevo.2017.00092

Sun, L., Miao, X., Cui, J., Deng, J., Wang, X., Wang, Y., Zhang, Y., Gao, S., & Yang, K. (2018). Genome-wide high-resolution mapping of DNA methylation identifies epigenetic variation across different salt stress in Maize (*Zea mays* L.). *Euphytica*, 214(2), 1–15. https://doi. org/10.1007/s10681-017-2076-0

Szulkin, M., Garroway, C. J., Corsini, M., Kotarba, A. Z., & Dominoni, D. M. (2020). How to quantify urbanization when testing for urban evolution? In M. Szulkin, J. Munshi-South, & A. Charmantier (Eds.), *Urban evolutionary biology* (pp. 13–35). Oxford University Press.

Szulkin, M., Munshi-South, J., & Charmantier, A. (2020). *Urban evolutionary biology*. Oxford University Press.

Taufique, S. K. T., Prabhat, A., & Kumar, V. (2018). Illuminated night alters hippocampal gene expressions and induces depressive-like responses in diurnal corvids. *European Journal of Neuroscience*, 48(9), 3005–3018. https://doi.org/10.1111/ejn.14157

Teranishi, M., Shimada, Y., Hori, T., Nakabayashi, O., Kikuchi, T., Macleod, T., & Mizuno, S. (2001). Transcripts of the HMH region on the chicken Z chromosome accumulate as non-coding RNA in the nucleus of female cells adjacent to the DMRT1 locus. *Chromosome Research*, 9(2), 147–165. https://doi.org/10.1023/A:10092 35120741

Theodorou, P., Radzevičiūtė, R., Kahnt, B., Soro, A., Grosse, I., & Paxton, R. J. (2018). Genome-wide single nucleotide polymorphism scan suggests adaptation to urbanization in an important pollinator, the red-tailed bumblebee (*Bombus lapidarius* L.). *Proceedings of the Royal Society B: Biological Sciences*, 285(1877), 20172806. https:// doi.org/10.1098/rspb.2017.2806

Thompson, C. J., Ahmed, N. I., Veen, T., Peichel, C. L., Hendry, A. P., Bolnick, D. I., & Stuart, Y. E. (2017). Many-to-one form-to-function mapping weakens parallel morphological evolution. *Evolution*, 71(11), 2738–2749. https://doi.org/10.1111/evo.13357

Thompson, K. A., Rieseberg, L. H., & Schlüter, D. (2018). Speciation and the city. *Trends in Ecology & Evolution*, 33(11), 815–826. https://doi. org/10.1016/j.tree.2018.08.007

Thompson, M. J., Capilla-Lasheras, P., Dominoni, D. M., Réale, D., & Charmantier, A. (2021). Phenotypic variation in urban environments: Mechanisms and implications. *Trends in Ecology & Evolution*. https://doi.org/10.1016/J.TREE.2021.09.009

Torné-Noguera, A., Pagani-Núñez, E., & Senar, J. C. (2014). Great Tit (*Parus major*) breath rate in response to handling stress: Urban and forest birds differ. *Journal of Ornithology*, 155(1), 315–318. https:// doi.org/10.1007/s10336-013-1025-5

Van Dongen, W. F. D., Robinson, R. W., Weston, M. A., Mulder, R. A., & Guay, P. J. (2015). Variation at the DRD4 locus is associated with wariness and local site selection in urban black swans. BMC Evolutionary Biology, 15(1), 253. https://doi.org/10.1186/s1286 2-015-0533-8

Viitaniemi, H. M., Verhagen, I., Visser, M. E., Honkela, A., van Oers, K., & Husby, A. (2019). Seasonal variation in genome-wide DNA methylation patterns and the onset of seasonal timing of reproduction in great tits. *Genome Biology and Evolution*, 11(3), 970–983. https:// doi.org/10.1093/gbe/evz044

Waterland, R. A., & Jirtle, R. L. (2003). Transposable elements: Targets for early nutritional effects on epigenetic gene regulation. *Molecular and Cellular Biology*, 23(15), 5293–5300. https://doi.org/10.1128/ mcb.23.15.5293-5300.2003
<figure><img src="image_8.png" /></figure>

Waters, S. A., Capraro, A., McIntyre, K. L., Graves, J. A. M., & Waters, P. D. (2018). The methylome of vertebrate sex chromosomes. *Genes*, 9(5), 230. https://doi.org/10.3390/genes9050230

Watson, H., Powell, D., Salmón, P., Jacobs, A., & Isaksson, C. (2021). Urbanization is associated with modifications in DNA methylation in a small passerine bird. *Evolutionary Applications*, 14(1), 85–98. https://doi.org/10.1111/eva.13160

Weaver, I. C. G., Cervoni, N., Champagne, F. A., D'Alessio, A. C., Sharma, S., Seckl, J. R., Dymov, S., Szyf, M., & Meaney, M. J. (2004). Epigenetic programming by maternal behavior. *Nature Neuroscience*, 7, 847–854. https://doi.org/10.1038/nn1276 https://doi.org/10.1038/nn127610.1038/nn1276

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. *Evolution*, 38(6), 1358. https://doi.org/10.2307/2408641

Zhou, X., & Stephens, M. (2012). Genome-wide efficient mixed-model analysis for association studies. *Nature Genetics*, 44, 821–824. https://doi.org/10.1038/ng.2310

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Caizergues, A. E., Le Luyer, J., Grégoire, A., Szulkin, M., Senar, J.-C., Charmantier, A., & Perrier, C. (2022). Epigenetics and the city: Non-parallel DNA methylation modifications across pairs of urban-forest Great tit populations. *Evolutionary Applications*, 00, 1–17. https://doi.org/10.1111/eva.13334