Continent-wide genomic signatures of adaptation to urbanisation in a songbird across Europe

Genetic diversity and population structure across European urban and rural populations

A total of 192 great tits from the nine paired urban–rural populations were genotyped at 517,603 filtered SNPs, with 10–16 individuals per sampling site (Supplementary Table 1). We quantified the relative degree of urbanisation for each site (urbanisation score: PC_urb, from principal component analysis, PCA; see “Methods”, Fig. 1b, Supplementary Fig. 1 and Supplementary Table 1) to inform our genetic downstream analyses. Population structuring based on 314,351 LD (linkage-disequilibrium)-pruned SNPs (excluding small linkage groups and the Z-chromosome) was overall low across the 18 studied sites (Supplementary Fig. 2), with each of the first two principal components explaining <3% of the overall variation across populations (Fig. 2a and Supplementary Fig. 3a). Thus, we used a UMAP (Uniform Manifold Approximation and Projection) approach to summarise the genetic variation along the first 20 PC axes. The UMAP analysis revealed the presence of distinct genetic clusters for some of the localities, although there was still a strong clustering of individuals from Gothenburg, Munich, Milan and Paris (Fig. 2a and Supplementary Fig. 3b; results were comparable when including the Z-chromosome, Supplementary Fig. 3c). This analysis also suggested that the levels of divergence differ strongly between each urban and adjacent rural population (Fig. 2a and Supplementary Fig. 3b). Furthermore, the population pairs from Glasgow and Lisbon showed the highest levels of divergence, with Lisbon and Glasgow separating along PC1 and PC2, respectively (Supplementary Fig. 3a), and also separating strongly in the UMAP plot (Fig. 2a). This increased divergence of Glasgow and Lisbon, both located in the range edge of the great tit distribution, could be explained by their slightly reduced heterozygosity, particularly in the urban populations (Supplementary Table 1). Indeed, population tree analyses using TreeMix supported the presence of increased drift in both urban populations (Fig. 2b). However, overall, we did not find lower heterozygosity levels in any of the nine-urban compared to the rural populations, which suggests that urban colonisation was in general not associated with significant bottlenecks (see Supplementary Table 1 for details; Wilcoxon test: W = 30, P = 0.377).

Interestingly, the population structure analysis indicated that multiple urban populations, namely Glasgow, Lisbon, Madrid, Malmö, Milan, and to weaker extent Barcelona and Gothenburg, formed distinct genetic clusters and independent drift units together, in some cases, with their adjacent rural counterparts (Fig. 2a, b). Nonetheless, the TreeMix analysis also suggests the presence of migration events across some populations, including long distance migration events, e.g., from Lisbon to Glasgow (Supplementary Fig. 4b). Although, there was no clear pattern in relation to predominant gene flow between urban populations (Supplementary Fig. 4). In contrast to the mentioned population clusters, we did not detect any particular signs of genetic divergence between populations from Munich and Paris (neither PC or UMAP axis see, e.g., Fig. 2a and Supplementary Fig. 3). Likewise, we detected significant genetic differentiation between all pairwise comparisons (P < 0.05; F_ST permutation test; see “Methods”) except for the urban populations in Munich and Paris (Fig. 2e and Supplementary Table 2). Moreover, the genetic population differentiation (F_ST) within population pairs (urban vs adjacent rural population), as well as between urban populations (urban vs urban), was on average higher than the differentiation between rural populations (rural vs rural) across Europe (Fig. 2c). These results support the idea that gene flow is generally stronger between rural compared to urban habitats in the nine population pairs, and indicates reduced gene flow between urban and adjacent rural populations. This is further supported by the observed isolation-by-distance patterns, which showed slightly stronger IBD for urban (Mantel’s r = 0.46, P = 0.008) compared to rural populations (r = 0.43, P = 0.007; Fig. 2d). Lastly, analysis of the effective migration surfaces confirmed reductions in gene flow in many parts of Europe, compared to neutral expectations, including those in close proximity, e.g., Gothenburg and Malmö, and simultaneously highlighted strong gene flow in Central/Western Europe, i.e., Munich and Paris (Fig. 2e).

Overall, our results reveal weak population structuring across Europe, with slightly increased genetic differentiation between urban populations compared to the more admixed rural populations. Previous results in another European songbird, the blackbird (Turdus merula) showed that urban birds likely independently colonised multiple times European cities from forest sources⁴². In our study, the finding of independent clusters for some urban populations also point towards a similar scenario, i.e., an independent and repeated colonisation of urban habitats from largely admixed rural populations. However, the overall weak genetic divergence and the detection of significant migration events across distant populations suggests a possible role of gene flow in the facilitation of urban adaptation, via the spread of adaptive alleles. Still, in our subsequent analyses we treated all our urban–rural population pairs as independent, including Munich and Paris, as selection pressures might differ across cities and result in idiosyncratic genomic responses to selection despite any ongoing gene flow.

Identification of urbanisation-associated allele frequency shifts

To identify genomic regions with consistent allele frequency shifts associated with urbanisation, we used two complementary genotype–environment association (GEA) approaches, LFMM (latent-factor mixed models) and an additional Bayesian approach (using BayPass). Testing for GEAs with LFMM (using PC_urb as a continuous habitat descriptor, Fig. 1b) revealed 2758 SNPs associated with urbanisation (0.52% of the full SNP dataset, false-discovery rate (FDR) < 1%; Fig. 3a). Urbanisation-associated SNPs were widely distributed across the genome and did not cluster in specific regions. Larger chromosomes (R² = 0.97) and those with more genes (R² = 0.90; Fig. 3b) contained more urbanisation-associated SNPs, highlighting the polygenic nature of urban adaptation. A PCA (PC_GEA), based on all urban-associated SNPs detected by LFMM, clearly separated urban and rural populations along PC1_GEA (proportion of variance explained—PVE—by PC1 = 1.98%; Supplementary Fig. 5c), suggesting consistent allele frequencies in those loci across European cities.

In contrast, the results obtained by BayPass only identified 70 SNPs strongly associated with urbanisation (Bayes factor ≥ 20; Supplementary Fig. 6a). The lower number of identified SNPs might respond to the stronger population structure correction applied by BayPass, which correlates to some extent with the direction of local adaptation⁴³. Of these significant SNPs, 34 were shared with the LFMM analysis (more than expected by chance: χ² = 31.20, P = 2.35 × 10⁻⁰⁸; Fig. 3a and Supplementary Fig. 6b). These shared SNPs, which we term “core urbanisation SNPs” (Supplementary Table 3), are likely involved in the local adaptation of great tits to urban habitats, and indeed, they more strongly discriminated urban and rural individuals across Europe (PVE by PC1_GEA-shared = 11.7%; Fig. 3c and Supplementary Fig. 7a).

In order to gain a better understanding of the importance of habitat (i.e., urban vs rural) on allele frequency shifts in the detected SNPs, we assessed using univariate linear models, the explanatory power of habitat, locality (city), and their interaction on the direction and strength of allele frequency shifts per-SNP across all paired populations. A significant “habitat” term indicates consistent allele frequency shifts across cities, particularly when the effect size (partial η²) is higher compared to the effect size of a significant “habitat × locality” interaction term, which describes differences in the direction and magnitude of the allele frequency change between cities⁴⁴. Applying these models to each urbanisation-associated SNP from the LFMM model (2758 SNPs), we detected large variation in allele frequency shifts across localities (Fig. 3d and Supplementary Fig. 7b), with a slightly larger proportion of SNPs showing differences in allele frequency shifts across localities (“habitat × locality” η² > “habitat” η²). In contrast to this, most of the “core urbanisation SNPs” showed similar allele frequency differences across localities, with 76% of them showing a main effect of the urban habitat (Fig. 3d), and with the same allele increased in frequency in seven or more urban populations (Supplementary Fig. 7c; see Supplementary Fig. 8 for the minor allele frequency trajectory between habitats). We refined this analysis by accounting for the effect of allele frequency correlations between SNPs through the use of the principal component axis from all the LFMM urbanisation-associated SNPs (PC1_GEA) and the “core urbanisation SNPs” (PC1_GEA-shared, Supplementary Fig. 7d). In both cases, the “habitat” term explained the majority of the total variation in allele frequency divergence (η²_{PC1 GEA, habitat} = 0. 87, P < 0.001; η²_{PC1 GEA-shared, habitat} = 0. 73, P < 0.001). In comparison, both the effect of “locality”, which corresponds to the distinct evolutionary history of the local populations (η²_{PC1 GEA, locality} = 0.59, P < 0.001; η²_{PC1 GEA-shared, locality} = 0.20, P < 0.001), and the interaction of “habitat × locality” (η²_{PC1 GEA, habitat × locality} = 0.13, P = 0.002; η²_{PC1 GEA-shared, habitat × locality} = 0.13, P = 0.001), explained much smaller proportions of the total variation (smaller effect size) than the habitat term by itself. In contrast, genetic variation along PC2 and PC3 was mostly inconsistent across localities and likely specific for each city (Supplementary Fig. 7d). Overall, the PCA-based results are in line with those obtained by the per-SNP results, reinforcing the idea of consistent continent-wide responses to urbanisation, in particular regarding the “core urbanisation SNPs”.

Together, we show that there is local adaptation to urban habitats in great tits, and that this occurred across Europe through shifts in allele frequency of the same loci. This might have occurred via the standing genetic variation observed for the species, as putatively adaptive alleles were shared across large parts of the species’ European distribution⁴¹, or the exchange of adaptive alleles through gene flow across urban populations. Nonetheless, in line with a polygenic basis of adaptation via subtle allele frequency shifts in several loci²⁶, the detected urban–rural differences were generally small (Supplementary Fig. 9) with only a few SNPs showing relatively strong allele frequency shifts with |∆AF| (allele frequencies difference between urban and rural) >0.5.

Signatures of selection in urban populations

While GEA approaches are powerful tools for identifying even subtle allele frequency shifts associated with local adaptation⁴³, they have difficulty detecting selective sweeps unique to one or few populations or sweeps of different haplotypes associated with the same gene. Selective sweeps have been associated with urban adaptation in other species, for example, in the case of New York City brown rats (Rattus norvegicus) and white-footed mice (Peromyscus leucopus)^4,45. Therefore, we further performed genome-wide scans of differentiation and selective sweep analyses in each of our studied populations to test if these contributed to the urban adaptation in great tits and, importantly, to explore if selective sweeps were population specific or widely shared across Europe.

Since the exchange of genetic material between urban and rural populations, and the selection pressures are likely recent and ongoing in great tits, we opted to use a cross-population statistic that can identify ongoing and recently completed hard and soft selective sweeps (XP-nSL)⁴⁶. Using this approach, we found that between 436 and 700 genomic windows (200 kb sliding windows with 50 kb steps, see “Methods”), which clustered into 127–173 wider genomic regions, showed signatures of selection in urban populations (Fig. 4a). Genomic outlier regions had an average size of 355.8 kb ± 251.6 SD, with the largest region (3.05 Mb) detected on chromosome 1 in Glasgow and were widely distributed across the genome with a few distinct population-specific peaks (Fig. 4a). The large number of genomic outlier regions are in line with the high number of significant GEAs across the genome and confirms the inferred polygenic nature of urban adaptation in great tits (Fig. 3a). However, in contrast to the highly concordant allele frequency shifts of GEA loci, XP-nSL values were in general weakly correlated across populations (Fig. 4b), and outlier windows were shared across a maximum of five populations (Supplementary Fig. 10). The top three shared windows clustered in a 300 kb region on chromosome 11 (16.05–16.35 Mb) and were detected in Gothenburg, Munich and the Iberian Peninsula (Barcelona, Madrid and Lisbon), with one partially overlapping outlier window detected in Glasgow. In general, the most pronounced sharing of outlier windows did not seem to occur between the geographically closest populations, as expected for shared genetic variation or adaptive introgression (see Supplementary Fig. 10). For example, Gothenburg and Paris, and Lisbon and Glasgow, showed the highest number of shared outlier windows (53 and 47, respectively). Indeed, the number of shared outlier windows between urban populations was not correlated with their genetic (Mantel’s r = 0.28, P = 0.150) or geographic distance (Mantel’s r = 0.05, P = 0.400), suggesting that shared adaptive genetic variation through introgression and gene flow between neighbouring urban populations likely did not facilitate urban adaptation. However, we cannot fully exclude the role of gene flow as migration was generally strong, even across large distances (Supplementary Fig. 4).

The observed landscape of variation in XP-nSL is in agreement with the genomic landscape of genetic differentiation (standardised Z-transformed F_ST [ZF_ST]) between each adjacent urban and rural pair across the same windows (Fig. 4b, Supplementary Figs. 11 and 12). ZF_ST outlier windows (ZF_ST > 4, equivalent to four SD) were largely population specific (Supplementary Fig. 12), suggesting that population-specific selective and demographic processes have partly shaped the genomic landscape⁴⁷. To strengthen our inference of the selective sweeps landscape associated with urban adaptation in our species, we additionally searched for older and completed selective sweeps using the haplotype-based Rsb statistic⁴⁸. Similarly to the results obtained using the XP-nSL statistic and ZF_ST, window-based Rsb scores were weakly and inconsistently correlated across populations (Fig. 4b and Supplementary Fig. 13). We detected fewer Rsb outlier windows (209–538) and clustered genomic outlier regions (39–98; average size: 396.9 kb ± 329.2 SD; max = 3.45 Mb on chromosome 3 in Gothenburg) compared to XP-nSL. This result suggests that ongoing/recent sweeps outnumber older and completed selective sweeps in urban great tits. Similar to XP-nSL, the Rsb outlier windows were only shared across a maximum of five populations and in general widely distributed across the genome (Supplementary Fig. 14).

Overall, the selection analyses showed that selective sweeps, and in particular those that are recent and ongoing (XP-nSL), were pervasive across the genomes of urban great tits, supporting our inference of a highly redundant polygenic adaptation to urban environments. The underpinning genomic changes associated with urban adaptation were largely population specific, likely the result of differential selective pressures, adaptation through complex multivariate phenotypic changes and/or selection on standing genetic variation^49,50. Nonetheless, we still detected shared selective sweeps across urban populations, suggesting the presence of common genomic responses to the urban selective pressures, either through repeated selective sweeps or sharing of adaptive genetic variants.

Evolutionary drivers of genetic differentiation and signatures of selection

The largely population-specific landscapes of differentiation are in line with observations in other young divergences, and have been previously explained by the weaker effect of linked selection at early stages of the divergence process⁴⁷. However, to exclude a driving role of non-adaptive processes, we corroborated that the detected signatures of selection were indeed caused by divergent selection and not by genetic drift or linked selection in low-recombination regions⁴⁷. To achieve this, we evaluated the correlation between signatures of selection (XP-nSL and Rsb) and differentiation (F_ST) with two proxies of genomic recombination rate, i.e., LD and intronic GC-content (both in 200 kb non-sliding windows). LD is correlated with population-specific recombination rates and the strength of selection, while intronic GC-content is a surrogate for the long-term recombination rate variation across the great tit genome, i.e., a higher GC-content will indicate a higher recombination rate^51,52. To confirm that both proxies provide similar pictures of the recombination landscape, we tested their correlation across all population pairs. Indeed, LD and GC-content, were (negatively) correlated in all populations (Fig. 4c), suggesting that the general effect of recombination on the genomic diversity landscape is broadly consistent.

Under a scenario of linked selection in low-recombination regions, we would expect strong and consistent positive correlations between GC-content or LD with our estimates of genetic differentiation and selection, i.e., F_ST, XP-nSL and Rsb^53,54,55. Yet, this was not the case, (Fig. 4c), which indicates that processes other than linked selection are the main driver for population differentiation in our study. Indeed, the overall patterns are more in line with divergent selection, as e.g., the population-specific selection peak on chromosome 1 in Glasgow is not associated with low-recombination regions (low GC-content) or a pronounced peak in LD (Fig. 4d). However, it is important to note that these correlations are positive in some populations, see e.g., XP-nSL and GC-content (Fig. 4c), and that some population-specific peaks are located on low-recombination regions, and likely caused by linked selection or other non-adaptive processes (Fig. 4d). While the contribution of non-adaptive processes likely differed across populations, we did not detect a general signature of linked selection as a major driver of divergence between urban and rural populations.

Divergent selection has been suggested to be the main driver in the early stages of differentiation in other avian systems^47,56, and it is not surprising that a similar incipient pattern is observed between urban–rural populations. Although other non-adaptive processes, such as stochasticity via population-specific drift, could also generate variable genomic landscapes in an urban context⁵⁷, that scenario would be characterised by a decreased gene flow and increased genetic differentiation with geographical distance between adjacent urban and rural populations. However, genetic differentiation between paired urban and rural populations did not significantly increase with geographic distance (R² = −0.11, P = 0.645; Fig. 4e), further excluding drift as the main process underneath the variability in the genomic and the genetic differentiation landscape between urban and rural populations. Nonetheless, genetic drift potentially led to increased genetic differentiation in some urban populations, i.e., Glasgow and Lisbon, which both showed reduced heterozygosity.

Hence, we conclude that the mosaic of genomic signatures of differentiation between urban and rural populations and selective sweeps in the studied populations are most likely driven by genomic responses to selection on standing genetic variation, with potentially minor contributions of linked selection in low-recombination regions or other idiosyncratic non-adaptive processes.

Genes and pathways associated with adaptation to urbanisation

While responses to selection at the haplotype level were highly variable across urban populations, selective sweeps might affect the same genes or different genes with similar functional impacts (e.g., same pathway). Such consistent genomic responses on the gene or functional level could explain the similarity in phenotypic responses to environmental variation despite a large variation on the haplotype level^14,15,58. Following this rationale, we identified genes associated with signatures of selective sweeps in each urban population, but without requiring outlier windows to overlap with each other. Between 976 and 1366 genes were associated with signatures of ongoing or recent selective sweeps (XP-nSL), and 362–923 with completed selective sweeps (Rsb). Of these genes, 64–207, and 5–76, were detected in at least two urban populations (based on XP-nSL or Rsb, respectively; Supplementary Figs. 15 and 16). The higher number of genes associated with recent and ongoing selective sweeps compared to complete selective sweeps is in line with the colonisation pattern of urban habitats, under ongoing or recent gene flow.

Although it is statistically unlikely to detect selective sweeps in the same gene in more than three urban populations by chance alone (permutation test; χ²₁ = 77.95, P < 2.2 × 10⁻¹⁶), we conservatively focused our functional interpretation only on those genes that were associated with signatures of selection (XP-nSL or Rsb) in more than half of all urban populations (≥5 urban populations), as these likely play more important roles in urban adaptation. In total, we detected 42 genes associated with recent or ongoing selective sweeps in at least five populations (XP-nSL: max. six populations), and 15 genes associated with older completed selective sweeps (Rsb: max. seven populations; Fig. 5a and Supplementary Table 4), with none of the genes detected by both selection statistics. It is noteworthy that we did not detect an increase in the number of shared genes under selection in relation to geographical (Mantel’s r = 0.27, P = 0.080) or genetic distance (Mantel’s r = 0.05, P = 0.480), which is in accordance with the results regarding shared outlier windows and further supports that many selective sweeps likely occurred independently in multiple urban populations. In addition, 12 of the 57 shared genes under selection were also associated with urbanisation in the LFMM analysis (Supplementary Table 4). Interestingly, two of these, GMDS and SLC6A15, were not only associated with complete selective sweeps (Rsb scores) in seven and six populations, respectively, but also previously shown to be differentially expressed in blood and/or liver between urban and rural birds from Malmö⁵⁹. Moreover, following Malmö as an example, one gene found to be differentially expressed in Watson et al.⁵⁹, VPS13A, was highly differentiated (ZF_ST) and associated with selective sweeps in that particular population, but importantly, also under selection in Madrid (Rsb) and Lisbon (XP-nSL) and in general, associated with urbanisation (LFMM analysis). This exemplifies that some genes show signatures of urban adaptation on a continental scale, making them strong candidates for adaptation in these urban centres.

It is noteworthy that despite finding signatures of selection overlapping the same gene in five or more urban populations (e.g., GMDS, HTR7 or CDH18), the exact locations of these signatures were not consistent (Fig. 5b, c). Under a shared selective sweep or adaptive introgression scenario, we would have expected the sharing of urban haplotypes between populations. However, haplotype plots do not show any noticeable clustering across urban individuals (e.g., GMDS or HTR7, Fig. 5b), suggesting that different adaptive haplotypes swept to high frequency in urban populations across Europe, which is in line with soft selective sweeps from standing genetic variation⁶⁰. Furthermore, we find that the strongest selective sweep signatures are, in some cases, not within the genes but up- or downstream and without a consistent pattern among populations. For example, signatures of selection around HTR7 were downstream in Glasgow and Munich, but upstream in Malmö, Paris, Madrid and Lisbon (Fig. 5c). This might suggest that regulatory changes rather than protein-coding ones are associated with urban adaptation and the nature of regulatory variation likely differs across populations. However, high-resolution genomic studies based on whole-genome resequencing coupled with multi-tissue functional genomic analyses are needed to resolve this question, which is beyond the scope of this study.

Many of the genes associated with urbanisation have previously been linked to behavioural divergence, and cognitive and learning functions, suggesting adaptive phenotypic shifts related to behaviour. In particular, HTR7 (Chr 7) codes for a receptor of the neurotransmitter serotonin, a pathway consistently identified as target of changes linked to urbanisation. Indeed, two previous studies in birds, one on European blackbirds using a candidate gene approach and another on burrowing owls (Athene cunicularia), using full-genome resequencing, also described changes in this pathway as the main difference between urban and rural populations^3,22. Both species also show consistent urban–rural differences in two behaviours associated with serotonin: down-regulation of the stress axis and altered nocturnal behaviour under artificial light at night, respectively^9,61,62. The same gene, HRT7, has also been directly linked to behavioural differences between migratory and resident rainbow trout (Onchorynchus mykiss)⁶³, supporting its potential role in behavioural changes in urban great tits. The CDH18 gene (Chr 2), is part of a superfamily of membrane proteins involved in synaptic adhesion and revealed as a candidate gene in phonological alterations in humans⁶⁴. The PTPRD and VPS13A genes (both in Chr Z), are suggested to be involved in hippocampal neural development⁶⁵ and were previously linked to bird navigation, flight performance and migratory behaviour^66,67,68. Moreover, two other genes under selection and associated with urbanisation are involved in the regulation of cognitive processes and behavioural disorders in vertebrates (DLG2 (ref. ⁶⁹) and NRXN3 (ref. ⁷⁰), Chr 1 and 5, respectively), and were also previously linked to urban divergence in burrowing owls²².

Thus, our selection analyses suggest that natural selection repeatedly acts on behavioural traits and sensory and cognitive performance, all previously shown to be among the most widespread differences between urban and rural wildlife populations^71,72. These findings were furthermore supported by gene ontology (GO) analysis of the 2758 urbanisation-associated SNPs (LFMM analysis), linked to 984 genes (1501 SNPs in genic regions). Accordingly, most of the GO terms were related to neural functioning and development (e.g., GO:0016358, FDR = 4.30 × 10⁻⁶), cell adhesion (e.g., GO:0098742, FDR = 3.43 × 10⁻⁶) and sensory perception (e.g., GO:0050954, FDR = 3.42 × 10⁻³; Fig. 5d and Supplementary Table 5). These GO terms were mainly clustered into two interacting networks, one related to sensory recognition and the other to neural development and cell adhesion (Supplementary Fig. 17). These findings reinforce the previous idea on cognitive and behavioural changes as key responses to urbanisation, as song structure and escape or distress behaviour have been previously shown to differ between urban and rural great tit populations across Europe^7,35,73. Nonetheless, whether this is the result of a genetic response to selection or phenotypic plasticity was to a large extent still unknown. Our present study suggests a strong genetic component for these putatively urban adaptative phenotypes in great tits across Europe. Detailed functional genomic and phenotypic analyses are now needed to understand the role of these genes and pathways in the adaptive divergence of urban and rural great tits and other songbirds.

Our study demonstrates clear genomic signals of local adaptation to urban habitats in a common songbird on a continent-wide scale across Europe. We found that a combination of polygenic allele frequency shifts and spatially varying selective sweeps are associated with adaptation to urban environments. Our results strongly suggest that a few genes, which have known neural developmental and behavioural functions, experienced selective sweeps in urban populations. This suggests a strong consistency in the functional processes associated with urbanisation, despite the fact that the underlying haplotypes are often not shared. Thus, our study exemplifies evolutionary adaptation to urban environments on a European scale, and highlights behavioural and neurosensory adjustments as important phenotypic adaptations in urban habitats.

Source: Ecology - nature.com