Local selection signals in the genome of blue tits emphasize regulatory and neuronal evolution

Understanding the genomic landscape of adaptation is central to understanding microevolution in wild populations. Genomic targets of selection and the underlying genomic mechanisms of adaptation can be elucidated by genome‐wide scans for past selective sweeps or by scans for direct fitness associations. We sequenced and assembled 150 haplotypes of 75 blue tits (Cyanistes caeruleus) of a single Central European population by a linked‐read technology. We used these genome data in combination with coalescent simulations (i) to estimate an historical effective population size of ~250,000, which recently declined to ~10,000, and (ii) to identify genome‐wide distributed selective sweeps of beneficial variants probably originating from standing genetic variation (soft sweeps). The genes linked to these soft sweeps, but also those linked to hard sweeps based on new beneficial mutants, showed a significant enrichment for functions associated with gene expression and transcription regulation. This emphasizes the importance of regulatory evolution in the population’s adaptive history. Soft sweeps were further enriched for genes related to axon and synapse development, indicating the significance of neuronal connectivity changes in the brain potentially linked to behavioural adaptations. A previous scan of heterozygosity–fitness correlations revealed a consistent negative effect on arrival date at the breeding site for a single microsatellite in the MDGA2 gene. Here, we used the haplotype structure around this microsatellite to explain the effect as a local and direct outbreeding effect of a gene involved in synapse development.

within a population, which includes both global and local adaptation.
Arguably, the direct selection history of a population can better predict the population's future evolutionary trajectory than a comparative history.
Large lists of specific protein-coding genes known to be affected by local selection have been compiled (e.g., for humans, Rees et al., 2020). Since the early 1960s, however, scientists have argued that much of biological diversity may have resulted from changes in gene expression rather than from changes at the protein sequence level (reviewed in Stern, 2016). Among the first evidence in favour of this idea was the observed weak correlation between the rate of protein sequence evolution and the rate of phenotypic evolution among different clades (Wilson et al., 1974). A multitude of studies have subsequently shown the importance of regulatory evolution by linking specific phenotypic trait variation to changes in specific gene expression networks (Stern & Orgogozo, 2008). Evidence has further accumulated showing that regulatory sequences that flank the protein-coding gene of interest (cis-acting), as well as specific regulatory proteins encoded elsewhere in the genome (trans-acting), control when expression of the gene of interest occurs, at what level it is expressed, under what environmental conditions, and in which cells or tissues (Wray et al., 2003). Most of these studies considered variation across species, while only a few investigated genomic mechanisms of adaptation within species (e.g., among ecotypes, Jones et al., 2012).
Here, we analysed historical and current selection signals in a large and stable natural bird population without a history of known specific environmental changes such as those experienced during island or urban colonizations (Johnson & Munshi-South, 2017).
Nevertheless, the population probably experienced unknown adaptive events during its long history, presumably resulting from abiotic or biotic environmental changes (e.g., human-induced habitat change, climate change, disease epidemics). We therefore expect to detect signals of selection at genes involved in regulation, but also at genes related to other known hotspots of selection, such as immune system genes (Shultz & Sackton, 2019) or genes related to behavioural changes . The latter is expected, because changes in behavioural traits may expose individuals to new environments and, vice versa, new environments may promote changes in behaviour (Zuk et al., 2014). Among the probable candidates for adaptive evolution of behaviour are genes that regulate neuronal connectivity in the brain .
We scanned the genome for signatures of selection in a wild population of blue tits (Cyanistes caeruleus) breeding in a natural, oak-dominated forest in Central Europe. The blue tit has long been a model species for microevolutionary studies, and direct measures of fitness in terms of reproduction and survival are available from long-term field studies (Culina et al., 2020). In previous studies we tested for general and local inbreeding and outbreeding effects using a genome-wide set of 83-106 microsatellite loci (Botero-Delgadillo et al., 2020;Olano-Marin et al., 2010;Olano-Marín et al., 2011a, 2011b. In the population considered here, a single microsatellite locus (H095_ TG05_053) showed a consistent negative correlation between individual heterozygosity and arrival date (Botero-Delgadillo et al., 2020). Even though this blue tit population is nonmigratory, individuals differ consistently in the date of yearly arrival at the breeding site . Moreover, arrival date can be considered a key fitness-related trait, because it correlates with yearly reproductive success .
Here, we characterized this potential negative outbreeding effect using sequence data in addition to the de novo scan for selection signatures.
For this study, we analysed 150 haplotypes from 75 blue tits sampled in a single population. The haplotypes were sequenced and assembled using a linked-read technology and the sequence data were mapped against the blue tit reference (Mueller et al., 2016).
We used these data to address the following questions: (i) What is the current and historical effective population size of the Central European blue tit population? The resultant demographic model will be used as the null model for the following selection tests. (ii) How abundant are genomic signatures of selective sweeps in the blue tit population? The number of different types of selective sweeps (hard and soft sweeps of de novo mutations or standing variation, respectively) will inform about the general background of adaptive events in this wild population. (iii) Are the detected selective sweeps located around genes with specific functions and is there evidence for functional hotspots of selection? (iv) Does the haplotype structure around the microsatellite locus H095_ TG05_053 with a known negative heterozygosity-fitness correlation help to explain an outbreeding effect? This will provide information for an alternative mode of selection, in addition to positive selection as evidenced by selective sweeps.

| Study population and sampling
The forest Westerholz (48°08′26″N, 10°53′29″E) represents a highquality breeding habitat for blue tits in southern Germany with an oak-dominated natural forest. In 2007, we installed 277 nestboxes and since then we have monitored the breeding biology of 60-170 breeding pairs each year. All birds were caught for measurement and blood sampling, either during the nestling period, during winter roosting in nestboxes or with mist-nets at feeders . Botero-Delgadillo et al., 2020). We strived for a balanced sex ratio in each sampling category (Table S1). The 35 randomly selected birds were used to infer demography and selective sweeps, whereas all 75 birds were used to characterize the haplotype structure at the microsatellite with known heterozygosity-fitness effect.

| Sequencing, reference mapping, SNP calling and phasing
Blood samples were stored in Queen's lysis buffer and DNA was extracted with the Blood QuickPure kit (Macherey-Nagel). Sequencing libraries for each bird were individually prepared using the linkedread Chromium technology of 10× Genomics at the National Genomics Infrastructure (NGI), Sweden. Samples were then individually sequenced on a NovaSeq6000 device with a 2 × 151-bp setup and controlled for quality at the same service.
We used the LongRanger pipeline for barcode processing, mapping against the blue tit reference genome (Mueller et al., 2016), variant calling and phasing of each sample (Zheng et al., 2016). The pipeline makes use of the Burrows-Wheeler aligner (BWa; Li, 2013) and Lariat aligner (Bishara et al., 2015) and the GATK HapLotypeCaLLer  (Table S1).
We merged all variants that passed the standard quality filter of the Longranger across all samples in a single variant call format file (vcf) using BCftooLs (Danecek et al., 2021). This merging produced missing genotype codes at population-wide variants in samples in which these variants were homozygous for the reference allele and thus not called. We replaced these missing codes with a homozygous reference code (0/0) when read depth at the site was ≥6 using fixVCfMissinggenotypes (Lindenbaum, 2015). After this procedure, we quality-filtered the data selecting autosomal variants with fewer than seven missing genotype codes, a minor allele count >1 and a read depth smaller than the overall mean plus 5SD. In the final step, we phased the consecutive phase blocks produced by Longranger (characterized by the PS tag in vcf files) within each scaffold using sHapeit4 (Delaneau et al., 2019). The final set of biallelic single nucleotide polymorphisms (SNPs) and indels within scaffolds larger than 110 kb (relevant for selective sweep analysis, see below) comprised 20,461,613 loci.

| Demographic inference
We used the software package popsizeaBC with an approximate Bayesian computation approach to estimate historical effective population sizes (N e ) (Boitard et al., 2016). First, we simulated 102,000 data sets of 70 15-Mb-long haplotypes using the coalescent-based simulator sCrM (Staab et al., 2015) assuming a mutation rate of 4 × 10 −9 and a recombination rate of 2.5 × 10 −8 . These values are estimates from a study on a closely related species, the great tit (Parus major) (Laine et al., 2016;van Oers et al., 2014). Effective population size is allowed to change randomly by a factor of 10 between 100 and ~400,000 individuals at 20 time points between three and 500 generations before the present. For each simulated data set, we then extracted the allele frequency spectrum (AFS) and average allelic linkage-disequilibrium (LD) values of 20 distance bins ranging from 41.4 kb to 13.5 Mb. We included long-range LD values to estimate recent population size, because theory predicts that LD between markers of more than 10 Mb distance on chromosomes with a recombination rate of 2.5 × 10 −8 is already affected by population size two generations ago (Boitard et al., 2016;Hayes et al., 2003).
Second, we extracted the observed AFS and LD values for the same distance bins from biallelic variants of the 23 scaffolds that were longer than 15 Mb, representing 51% of the assembled genome.
Third, we compared the simulated summary statistics with the observed ones and identified the most similar simulation with the simple rejection method (acceptance rate 0.0001) using the R package aBC (Csillery et al., 2012;R Core Team, 2020). We plotted the median N e over time assuming an average generation time of 2 years (own observations, see also Laine et al., 2016 for the great tit). We tested the robustness of demographic inference using two alternative methods which are insensitive to potential phase-switch errors. The first one is based on the same procedure as described above (pop-sizeaBC), but used genotypic LD instead of allelic LD as the summary statistic. The second is a sequentially Markovian coalescent method which includes information of sample frequency spectra to estimate more distant histories. Here, we inferred population size history between 100 and 100,000 years before the present using the program sMC++ (Terhorst et al., 2017). Estimates were performed for each of the 23 scaffolds larger than 15 Mb after excluding all stretches of unknown nucleotides (Ns).

| Scanning for selective sweeps
We scanned the genome for selective sweeps in consecutive 10-kb windows using dipLos/HiC in haploid mode (Kern & Schrider, 2018).
Two types of sweeps were distinguished: hard sweeps originating from de novo beneficial mutations, and soft sweeps originating from standing variation which becomes beneficial (Messer & Petrov, 2013). To characterize sweep evolution, the supervised machine-learning algorithm employs 12 major summary statistics (normalized values) of selective sweeps including diversity measures, allele frequency spectra, LD patterns and haplotype structure in the focal and 10 neighbouring 10-kb windows (Kern & Schrider, 2018). For the learning step we produced 46,000 genomic simula-

| Functional enrichment analysis
We tested for enrichment of associated gene functions of all windows according to their selective sweep probabilities. We are particularly interested in functions associated with higher mean ranks than expected by chance. Thus, we used the Wilcoxon rank-sum test implemented in the R package gofunCr (Grote, 2018) on all the windows ranked according to their sweep probabilities. Windows were annotated to gene ontology (GO) terms via overlapping (or lying within a maximum distance of 5 kb to) genes (whole regions including exons and introns) of our blue tit genome assembly (NCBI blue tit annotation release 100). We used the reduced list of 3336 GO-slim terms defined by pantHer (Huaiyu et al., 2021) as being informative of function and evolutionarily conserved. Using this method, we annotated 39,915 of the 79,334 windows (50.3%) to GO terms. We corrected for testing multiple GO terms within each of the three main GO families "biological process," "cellular component" and"molecular function" by evaluating 10,000 simulated data sets with permuted rankings. The family-wise error rate (FWER) was defined as the fraction of random sets where the lowest p-value across all GO terms was lower than or equal to the raw p-value of the tested GO term. To reduce redundancy of significant GO sets, we simplified by grouping similar terms based on their semantic similarity according to the Relevance parameter (Schlicker et al., 2006) using the R package rrVgo (Sayols, 2020).

| Evaluation of a local heterozygosityfitness effect
As mentioned above, individual heterozygosity at a single microsatellite locus (H095_TG05_053) correlated negatively with arrival date at the breeding site (a fitness-related trait) in the study population Because the indels did not pass the quality filter, we checked and corrected the automatic calling by manually scoring the genotypes through direct inspection of the phased read alignments (phased bam-files). Using the same procedure we further phased the indels in relation to the surrounding already-phased variants (see above).
We assessed LD in the genomic region of interest by calculating pairwise allele and heterozygosity correlations between the neighbouring variants and the indels or direct scores of the microsatellite locus using the R package pegas (Paradis, 2010). We further investigated the genetic relationship among haplotypes in the microsatellite region by calculating simple sequence distances between haplotypes using the R packages pegas and HapLotypes (Aktas, 2020). Each different allele of SNPs or indels was scored as one distance unit. Coulson et al. (1999) suggested that within-individual allele distance reflects the degree of outbreeding better than simple heterozygosity. Thus, we repeated the single-locus analyses described in Botero-Delgadillo et al. (2020) for the larger data set (635 observations from 459 PIT-tagged adult individuals), but using the internal distance between the two alleles instead of heterozygosity as a predictor of arrival date. For each of the 83 microsatellite loci we fitted mixed-effects models with arrival date (standardized within each season) as the response variable and allele distance (absolute value), sex and whether the individual bred in the population in previous years (yes or no) as predictors using the R package LMertest (Kuznetsova et al., 2017). Individual identity was entered as a random factor. F I G U R E 1 Estimated effective population size (N e ) from the present to more than 1000 years ago (~500 generations) for our Central European blue tit population. Shown are the median (solid black line) with its 90% credible intervals (dotted lines) estimated by the program popsizeaBC in haploid mode and the simplified trajectory as used for the simulations with selection (orange)

| Demographic inference, selection signals and functional enrichment
The historical effective population size of our Central European blue tit population was estimated to be ~250,000 individuals with a declining trend starting about 180 years ago (Figure 1). The estimated current effective population size is ~10,000 individuals.
These estimates are based on the summary statistics of fully phased haplotypes (popsizeaBC in haploid mode). Similar trajectories of demographic dynamics were estimated by two alternative methods which are independent of phasing information (Figures S1 and S2). This indicates an insignificant rate of phase switch errors in our data for haplotypic analyses.
We estimated sweep probabilities for 79,334 10-kb windows in 119 autosomal scaffolds of sufficient size and variant density.
Of these, 310 windows showed evidence for hard sweeps and 1775 windows for soft sweeps with prediction probabilities >.8.
The remaining windows were predicted to be neutral or linked to sweeps, or the evidence was inconclusive (low probabilities for all scenarios). The hard and soft sweep windows are distributed across the whole genome ( Figure 2). There is also no clustering of hard or soft sweeps as shown by the distributions of intrachromosomal distances ( Figure S3).
The higher ranking of hard sweep windows (higher probabilities) was significantly associated with genes that are functionally enriched for a total of 128 GO terms (Table S2). In other words, the windows associated with these GO terms showed on average higher sweep prediction probabilities than expected by chance. Most of the significant gene-rich GO terms within the major ontology "biological process" relate to the function of "gene expression regulation" (Figure 3). This includes functions of "transcription regulation," "RNA processing and splicing" and post-translational regulation ("ubiquination"). Other, less gene-rich groups of significant GO terms relate to "cell cycle/division," "protein complex assembly," "cell/cilium movement," "response to hypoxia," "calcium ion transport" and "reproduction/gamete generation." Soft sweep windows highly ranked by prediction probabilities were strongly enriched for genes associated with 28 GO terms (Table S3). Similar to hard sweeps, most blue tit genes annotated to GO terms within the ontology "biological process" were linked to "gene expression" and "transcription regulation" (Figure 4). Another prominent group of genes were related to "neuron projection/axonogenesis" and "synapse organization." Less gene-rich groups were associated with "cell/cilium movement" and "calcium ion transport," similar to hard sweeps.

| Evaluation of a local outbreeding effect
We summed the total length of each sequenced haplotype considering the four indels in the originally amplified microsatellite region of H095_TG05_053 (162 bp in the reference). These combined allele lengths were consistent with the previously determined microsatellite allele length scores for 87.3% of the 150 alleles in the 75 sequenced blue tit individuals. In terms of heterozygosity (yes/no), 90.7% of the 75 individuals were consistently scored between the two methods. The observed discrepancies between alleles scored by sequencing and by fragment amplification are probably due to some failed variants and mis-scoring of indels in the poorly mapped region with many sequence repeats and mononucleotide stretches.
Linkage-disequilibrium in the microsatellite region was generally low and the correlations between the microsatellite alleles and the alleles of the neighbouring SNPs and indels were generally weak ( Figure 5a). This was also expressed in the low correlation of heterozygosity between the loci (Figure 5b). Haplotypes of 5 kb length across the microsatellite region (chr5_part1: 1366000-1371000) formed one diffuse cluster in terms of genetic relationship ( Figure S4).
We therefore conclude that the microsatellite locus stands isolated F I G U R E 2 Manhattan plot showing for each tested genomic window the selective sweep probability of hard sweeps (in red) and soft sweeps (in blue)

Scaffolds < 10 Mb
Scaffolds in terms of general genetic variation, but also in terms of potential rare genetic variation (i.e., there are no rare distantly related haplotypes; Figure S4). This finding precludes a simplification of the microsatellite scoring by a specific haplotype grouping. We therefore tested for potential outbreeding effects in the larger sample (635 observations of 459 individuals) by using internal allele distance of each of 83 microsatellites as predictors for arrival date.
H095_TG05_053 showed the strongest and the sole panel-wide significant association between allele distance and arrival date (corrected p = .026, Table S4). Individuals with higher allele distances arrived later in the season.

| DISCUSS ION
Understanding the genomic landscape of signals of past selection is important for the assessment of genomic variation in microevolutionary and population genetic studies in model species such as the blue tit. To analyse the selection history from genomic variance data it is essential to simultaneously infer demographic history, because selection and demographic processes can both produce similar genomic patterns, and interact in the formation of such patterns (Parul et al., 2020). We applied a method with long-range LD measures (up to 13.5 Mb), which allowed us to estimate the recent demographic history of a blue tit population from a few generations ago to more than 1000 years ago. Our ancestral estimate of an effective population size of ~250,000 individuals lies within the range of estimates from a set of models of divergence history between mainland and island populations of blue tits (Perrier et al., 2020). These authors estimated a mean effective population size of 397,685 for the mainland in southern France. Our estimated decline in effective population size, starting about 180 years ago, could be a consequence of the increased replacement of natural mixed, oakrich forests (and other suitable habitats for blue tits) by pure coniferous forests in the area (Freitag, 2012).

F I G U R E 3
Treemap of significantly higher ranked GO terms according to hard sweep probabilities in comparison with randomly ranked terms (based on simulations, see Section 2). Semantically similar terms of the ontology "biological process" are grouped by colour and named by a leading term (in white). The space used by each GO term is proportional to the number of annotated blue tit genes (N total = 3381 genes) By comparing 12 major parameters as indicators of selection across the focal and flanking windows between empirical data and simulations based on the observed demographic history, we identified about 2000 windows with high probabilities of either a hard or a soft selective sweep. These selective sweep windows, both hard and soft, were widely distributed across all chromosomes and showed no clear clustering. Soft sweeps were more common than hard sweeps (85% vs. 15%). This is in line with recent empirical studies showing the predominance of soft selective sweeps in populations of Drosophila and humans (Garud et al., 2015;Schrider & Kern, 2017, but see also Jensen 2014). Abundant recent soft sweeps have also been described in a set of closely related species of southern capuchino seedeaters (Sporophila spp.), albeit in the context of sexual selection on pigmentation genes following speciation (Hejase et al., 2020). It has been argued that selection on standing genetic variation may be more likely than selection on new mutants, because established variants have already been filtered and shaped by selection in their ecological context (Hejase et al., 2020). Here, alleles under balancing selection have great potential to mediate quick genetic adaptations when the population encounters new environments, because they necessarily affect phenotype and fitness, and are already at intermediate frequencies at the onset of selection (De Filippo et al., 2016;Rees et al., 2020).
The genes associated with the highly ranked selective sweeps (ranked by probability) showed strong functional enrichments.
Functions related to gene expression were strongly overrepresented in many genes associated with hard and soft sweeps. The results suggest that the full range of functional modules related to the regulation of gene expression are overrepresented, that is from transcription initiation, elongation, termination, transport, RNA processing and differential alternative splicing to post-translational modification such as protein ubiquitination. Ubiquitination in F I G U R E 4 Treemap of significantly higher ranked GO terms according to soft sweep probabilities in comparison with randomly ranked terms (based on simulations, see Section 2). Semantically similar terms of the ontology "biological process" are grouped by colour and named by a leading term (in white). The space used by each GO term is proportional to the number of annotated blue tit genes (N total = 955 genes) different forms and degrees (the so-called "ubiquitin code") has recently emerged as an omnipresent signalling mechanism that controls diverse physiological processes including gene expression via histone ubiquitination (Kliza & Husnjak, 2020;Swatek & Komander, 2016). A large set of over-represented genes is also specifically related to the transcription of mRNA by the typical RNA polymerase II (Sims III et al., 2004). As expected, functions related to the first steps in gene expression such as transcription initiation were often over represented, because regulation at these steps avoids the costs of follow-up processes such as translation (Wray et al., 2003). Overall, these results emphasize the importance of regulatory evolution in the history of this blue tit population. At the between-species level, a direct relationship between when, where or how strong a gene is expressed and functionally significant phenotypes has been shown for many genes and gene regulatory networks (Wray et al., 2003). Regulatory changes also appear to predominate adaptive changes among populations or ecotypes of threespine sticklebacks (Gasterosteus aculeatus) that have repeatedly transitioned from marine to freshwater environments (Jones et al., 2012). Our method of GO annotations only allows us to detect trans-acting gene expression regulation, such as via transcription factors or other gene products. Interestingly, trans-regulatory variation contributes more strongly to expression differences within species than between species (Hill et al., 2021). Given the larger mutational target size, trans-regulatory factors might be less affected by the depletion of genetic variation through selective sweeps within populations than cis-regulatory regions. Overall, significant regulatory evolution might have been Indeed, there is strong evidence for compensatory evolution of gene expression in many organisms, also among trans-acting regulatory factors (Signor & Nuzhdin, 2019).
High-ranked soft sweep windows also showed a strong association with many genes acting in transcription regulation. This further supports the importance of selection on regulatory functions in the evolutionary history of this population. In addition, functions related to neuron development, in particular axon and synapse development, were also over represented around soft sweeps. In general, both axon and synapse development can modify brain connectivity and signalling and thus can influence behavioural variation. There is abundant evidence in vertebrates that the rich variation in synapse and neurogenesis genes (and proteins) associates with variation in several behavioural and cognitive traits (Grant, 2016;Lam et al., 2017). A previous study using a genome scan showed that genes acting in synapses and neuron projections are early targets of selection during recent urban colonizations by burrowing owls (Athene cunicularia) . Several other studies have reported allele frequency shifts in specific synaptic genes during colonization or urbanization in various bird species (van Dongen et al., 2015;Mueller et al., 2013Mueller et al., , 2017. Selection tests across different habitats in blue tits and great tits also revealed an enrichment of neuronal functions among the significant associated genes (Laine et al., 2016;Perrier et al., 2020).
In addition to signals of past positive selection, expressed by selective sweeps, we also identified a different type of selection, namely a negative correlation between heterozygosity at a single microsatellite locus and the fitness-related trait "arrival date at the breeding ground" (Botero-Delgadillo et al., 2020). Here we show through haplotype analysis that this negative outbreeding effect of the specific microsatellite is local and potentially a direct cis-acting effect in the 3′UTR of the gene MDGA2 (MAM domain containing glycosylphosphatidylinositol anchor 2). Linkage with multiple rare partial recessive deleterious alleles-another potential cause for negative heterozygosity-fitness-correlations-appears unlikely in this case . Interestingly, MDGA2 interacts with neurolignins and neurexins, which promote synapse development and validation by forming trans-synaptic bridges spanning the synaptic cleft (Gangwar et al., 2017). MDGA2 is overexpressed in the brain and testis, and variants have been associated with harm avoidance and personality disorders in humans (Heck et al., 2011). Of note, Heck et al. (2011) found a negative heterosis effect for one of the intronic SNPs in this gene. Thus, MDGA2 seems a good candidate gene for the study of local adaptation in behaviour.
In summary, a genome-wide scan for signals of selection revealed evidence for selection on gene regulation in the evolution of a blue tit population. Our data further support the idea that genes modulating neuronal connectivity in the brain are important targets of recent soft selective sweeps and outbreeding effects.

ACK N OWLED G EM ENTS
We thank A Türk, A Wittenzellner, P Loës, P. Skripsky and the many assistants for fieldwork, S. Kuhn, A. Girg, C. Baumgartner, F.
Taborsak-Lines for laboratory work, R-A. Olsen for bioinformatic work, Y. Pei for the drawing of the blue tit in Figure 1 and two anonymous reviewers for constructive comments. Open access funding enabled and organized by ProjektDEAL.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The blue tit genome browser can be found at http://publi c-genom