# The role of recombination dynamics in shaping signatures of direct and indirect selection across the Ficedula flycatcher genome

Madeline Chase, Maurine Vilcot, Carina Mugal

► To cite this version:

Madeline Chase, Maurine Vilcot, Carina Mugal. The role of recombination dynamics in shaping signatures of direct and indirect selection across the Ficedula flycatcher genome. Proceedings of the Royal Society B: Biological Sciences, 2024, 291 (2015), pp.20232382. 10.1098/rspb.2023.2382. hal-04497308

HAL Id: hal-04497308

https://cnrs.hal.science/hal-04497308v1

Submitted on 17 Oct 2024

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>

Distributed under a Creative Commons Attribution 4.0 International License
PROCEEDINGS B

royalsocietypublishing.org/journal/rspb

<figure><img src="image_2.png" /></figure>

Research

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for
updates

**Cite this article:** Chase MA, Vilcot M, Mugal CF. 2024 The role of recombination dynamics in shaping signatures of direct and indirect selection across the *Ficedula* flycatcher genome†. *Proc. R. Soc. B* **291**: 20232382.

https://doi.org/10.1098/rspb.2023.2382

Received: 20 October 2022

Accepted: 14 December 2023

**Subject Category:**

Evolution

**Subject Areas:**

evolution, genomics

**Keywords:**

meiotic recombination, linked selection, direct selection, speciation, Hill-Robertson interference

**Author for correspondence:**

Carina F. Mugal

e-mail: carina.mugal@ebc.uu.se

†An earlier version of this article is available as a preprint on bioRxiv: https://doi.org/10.1101/2022.08.11.503468 [1].

Electronic supplementary material is available online at https://doi.org/10.6084/m9 figshare. c.7005761.

# The role of recombination dynamics in shaping signatures of direct and indirect selection across the *Ficedula* flycatcher genome†

Madeline A. Chase<fnref n="1" />,<fnref n="2" />, Maurine Vilcot<fnref n="1" />,<fnref n="3" /> and Carina F. Mugal<fnref n="1" />,<fnref n="4" />

<figure><img src="image_5.png" /></figure>

MAC, 0000-0002-7916-3560; MV, 0000-0001-8765-1364; CFM, 0000-0003-4220-4928

Recombination is a central evolutionary process that reshuffles combinations of alleles along chromosomes, and consequently is expected to influence the efficacy of direct selection via Hill-Robertson interference. Additionally, the indirect effects of selection on neutral genetic diversity are expected to show a negative relationship with recombination rate, as background selection and genetic hitchhiking are stronger when recombination rate is low. However, owing to the limited availability of recombination rate estimates across divergent species, the impact of evolutionary changes in recombination rate on genomic signatures of selection remains largely unexplored. To address this question, we estimate recombination rate in two *Ficedula* flycatcher species, the taiga flycatcher (*Ficedula albicilla*) and collared flycatcher (*Ficedula albicollis*). We show that recombination rate is strongly correlated with signatures of indirect selection, and that evolutionary changes in recombination rate between species have observable impacts on this relationship. Conversely, signatures of direct selection on coding sequences show little to no relationship with recombination rate, even when restricted to genes where recombination rate is conserved between species. Thus, using measures of indirect and direct selection that bridge micro- and macro-evolutionary timescales, we demonstrate that the role of recombination rate and its dynamics varies for different signatures of selection.

## 1. Introduction

Meiotic recombination is a central evolutionary process with many effects on genome evolution. Recombination can create novel combinations of alleles that aid adaptation, which is hypothesized to contribute to the origin of sexual reproduction [2,3]. However, recombination may also be a disadvantage by breaking apart existing adaptive associations. The suppression of recombination can therefore occasionally promote local adaptation and speciation [4,5], for instance when inversions capture multiple loci with beneficial variation. Besides breaking up linkage among sites, recombination is also associated with the process of GC-biased gene conversion (gBGC), which leads to the preferential fixation of G:C over A:T alleles and can interfere with fitness [6,7]. Thus, recombination can have a multi-faceted impact on fitness.

Variation in recombination rate across the genome can also play a major role in shaping genomic signatures of natural selection. Recombination rate can impact the efficacy of selection via Hill-Robertson interference (HRI), where

© 2024 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

<footnote n="1">¹Department of Ecology and Genetics, Uppsala University, 75236 Uppsala, Sweden</footnote>

<footnote n="2">²Swiss Ornithological Institute, 6204 Sempach, Switzerland</footnote>

<footnote n="3">³CEFE, University of Montpellier, CNRS, EPHE, IRD, 34293 Montpellier 5, France</footnote>

<footnote n="4">⁴Laboratory of Biometry and Evolutionary Biology, University of Lyon 1, CNRS UMR 5558, 69622 Villeurbanne cedex, France</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Outline of the study design and working hypothesis. (a) Topology of the species included in the study, with approximate divergence times based on demographic modelling. Species that were used as outgroups or as reference species for estimating pairwise $F_{ST}$ are represented by dashed black lines. To the right of the topology we indicate the timescale in colour, from a macroevolutionary (blue) to a microevolutionary (orange) timescale. The corresponding colour code is used for the measures of natural selection that we estimate (b). (c) Summary of the predicted relationships between different measures of selection and recombination rate and density of selected sites in the presence of linked selection and Hill-Robertson interference (HRI). Predictions are presented for a stable recombination landscape versus a dynamic recombination landscape, where correlations are shown in parentheses if we expect them to be weakened by a dynamic recombination rate.</figcaption></figure>

linkage between multiple non-neutral mutations causes selective interference [8]. HRI predicts that the efficacy of natural selection should increase with increasing recombination rate, reflected by an increased fixation of deleterious mutations where recombination rate is low and increased fixation of beneficial mutations where recombination rate is high. Typical measures used to assess the presence of HRI are the nonsynonymous over synonymous ratios of diversity $\pi_N/\pi_S$ and divergence $d_N/d_S$, which are predicted to show a negative relationship with recombination rate, and the adaptive substitution rate $\omega_a$, predicted to show a positive relationship with recombination rate (figure 1c). Physical linkage among sites also manifests in a reduction of genetic diversity at neutral sites that are linked to targets of selection. The process by which neutral diversity is reduced through physical linkage to a beneficial mutation was originally referred to as genetic hitchhiking (HH) [9], leaving the signature of a selective sweep [9–11]. The size of a selective sweep depends on both the strength of selection and the recombination rate, with the reduction in diversity rapidly decreasing with increasing recombination [12,13]. This means selective sweep signatures are predicted to correlate negatively with recombination rate (figure 1c). Neutral diversity can also be reduced at sites linked to novel deleterious mutations by background selection (BGS), with a greater impact where recombination rate is lower [14]. We will jointly refer to the impacts of selective sweeps and BGS as indirect selection, to emphasize the feature that they refer to the indirect effects of selection on linked, neutral sites. Although selective sweeps and BGS leave several similar genomic signatures, there are key features that allow the distinction of, in particular, hard selective sweeps from BGS, such as the derived allele frequency spectrum [10]. Nevertheless, the diversity reducing effects of both types of indirect selection can lead to increased genetic differentiation ($F_{ST}$) between species [15,16], resulting in a negative relationship between $F_{ST}$ and recombination rate (figure 1c).

Many studies have investigated the relationship between recombination rate and genomic signatures of indirect selection. Across a wide range of taxa, variation in genetic

diversity is strongly correlated with recombination rate [17–19]. Differentiation islands, genomic regions showing significantly elevated $F_{ST}$ relative to the genomic background, frequently correspond to regions of reduced recombination [20–23]. Reproductive barrier loci and reductions in introgression have also been observed in regions of low recombination [24–27]. Altogether, these observations point to a pervasive role of recombination rate shaping signatures of indirect selection and correspond well with the expectations from population genetic theory.

On the other hand, evidence for recombination rate shaping signatures of direct selection via HRI appears less consistent across different taxa. HRI has been suggested to play a role in shaping patterns of molecular evolution in many systems, including, among others, the invertebrates *Drosophila melanogaster* [28,29], *Drosophila pseudoobscura* [30], *Heliconius melpomene* [31], the vertebrate great tit *Parus major* [32] and the plant sorrel (*Rumex hastatulus*) [33]. However, in other systems, including collared flycatcher (*Ficedula albicollis*), there has been no evidence of HRI shaping the efficacy of selection across the genome [34–37]. The intensity of HRI relies on multiple interacting variables, which can differ greatly among divergent taxa and could explain the mixed evidence for a relationship between recombination rate and signatures of direct selection. Both the density of functional sites and mutation rate contribute to the strength of HRI because interference requires multiple selected mutations to be linked and segregating at the same time. Neutrally evolving sites occurring in between selected sites (i.e. introns in between exons) can lessen the impact of interference by increasing the probability of recombination between interfering sites [38,39]. Thus, the impact of HRI on measures of direct selection may be more apparent in more compact genomes, such as invertebrate genomes, where gene density is higher compared with vertebrates and plants. Additionally, a relationship between recombination rate and genomic signatures of direct selection may not have the chance to build up if the recombination landscape evolves rapidly [30].

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
indirect and direct selection, we collated whole genome re-sequencing data for four species of *Ficedula* flycatchers: 65 taiga flycatchers (*Ficedula albicilla*) [51], 15 red-breasted flycatchers (*Ficedula parva*) [51], 95 collared flycatchers (*Ficedula albicollis*) [52] and 11 pied flycatchers (*Ficedula hypoleuca*) [20]. Additionally, we included one individual of snowy-browed flycatcher (*Ficedula hyperythra*) [20] as an out-group for variant polarization. We performed variant calling on all five flycatcher species and identified in total 51 424 863 single nucleotide variants (SNVs) within a set of 566 724 393 callable sites. We then focused our study on the more distant comparison of taiga flycatcher and collared flycatcher, while red-breasted flycatcher and pied flycatcher are only included for SNV polarization and to estimate pairwise $F_{ST}$ (figure 1).

To assess how the demographic history of our study system might influence our results, we performed demographic modelling of the divergence between the two closely related pairs of species using the software GADMA [53] and reconstructed species-specific demographic history using PSMC [54]. Recombination rate estimates based on linkage disequilibrium (LD) have been shown to be influenced by recent demography, where fine-scale variation and recombination hotspot detection appear to be sensitive to population bottlenecks [55], while broad-scale estimates of recombination rate may be more robust to demography [56]. Recent work also suggests that gene flow between species could confound genome-wide estimates of recombination rate from LD-based methods [57]. Although the best fitting model identified by GADMA included gene flow for both species pairs (electronic supplementary material, table S1 and figure S1), the effective migration rates estimated here appear sufficiently low as not to show an impact on LD-based recombination estimates [57]. However, owing to the sensitivity of fine-scale recombination rate on demography, in this study we focus our analysis on broad-scale variation in recombination rate across the genome, at a resolution of 200 kb windows, which also enables us to validate our findings based on LD-based recombination rate estimates by pedigree-based recombination rate estimates for collared flycatcher [58].

## (b) The evolutionary dynamics of recombination rate among taiga and collared flycatchers

To assess the evolutionary dynamics of recombination rate in *Ficedula* flycatchers, we estimated recombination rate in taiga flycatcher and collared flycatcher based on patterns of LD. LD-based estimates for both taiga flycatcher and collared flycatcher showed significant correlations with pedigree-based estimates for collared flycatcher at 200 kb resolution [58], with Pearson's correlation coefficients *R* = 0.46 and *R* = 0.58, *p*-value < 10<sup>-16</sup>, respectively (electronic supplementary material, figure S2). The strength of the correlation increased with increasing window size and was highest in 5 Mb windows (*R* = 0.77 and *R* = 0.76, *p*-value < 10<sup>-16</sup> for taiga flycatcher and collared flycatcher; electronic supplementary material, figure S3). These observations and the strength of the correlations are in good agreement with observations based on pedigree maps in other organisms [59,60], where it has been shown that the resolution of pedigree-based estimates can limit the strength of the correlation with LD-based estimates of recombination rate [60,61]. We next converted LD-based recombination rate from population-scaled

The dynamics of recombination rate evolution are highly variable across different taxa. In many mammals, a rapid turnover of recombination hotspots is observed. The fast evolution of hotspots is mediated by the zinc-finger protein PRDM9 [40–42], which encodes the location of hotspots and is rapidly evolving in response to the erosion of hotspot motifs [43–45]. In organisms that lack PRDM9, such as birds, recombination hotspots tend to be more stable over evolutionary time and the recombination landscape evolves at a slower pace [46,47]. Broad-scale variation in recombination rate also changes over time, for example, owing to chromosomal rearrangements [48–50]. However, despite a broad acknowledgement that recombination rate is dynamic, only few studies have directly investigated the impact of recombination rate evolution on genomic signatures of natural selection. Given the clear importance of recombination rate in shaping genomic signatures of selection, more investigation is warranted into how the evolutionary dynamics of recombination rate among species impacts these genomic signatures.

Here we investigate the role of recombination rate dynamics in shaping patterns of indirect and direct selection in two species of *Ficedula* flycatchers (figure 1a), taiga flycatcher and collared flycatcher. These two species diverged at an estimated 2.5 Ma, share a negligible amount of polymorphisms (approx. 2%), and show no recent history of gene flow [51], and both have a more closely related sister species, red-breasted flycatcher and pied flycatcher, respectively, for which population re-sequencing data are available to compute $F_{ST}$ landscapes (figure 1a). Therefore, these two species constitute an excellent study system to infer evolutionary changes in local recombination rate in birds and to assess their impact on genomic signatures of selection. Specifically, we estimate lineage-specific recombination rates in taiga flycatcher and collared flycatcher, as well as four measures of selection, with two measures each representing genomic signatures of indirect and direct selection, using both polymorphism-based and divergence-based measures, which reflect different timescales (figure 1). These measures include the composite likelihood ratio (CLR) test for selective sweeps, genetic differentiation $F_{ST}$, which in *Ficedula* flycatchers is suggested to be governed by linked selection [20,51], $\pi_N/\pi_S$, and $d_N/d_S$ (figure 1b). Additionally, we estimate the adaptive substitution rate ($\omega_a$), which relies on a combination of polymorphism and divergence data. We then investigate the relationship between evolutionary changes in recombination rate and these measures of selection. We also study the association between selective sweeps and $\omega_a$, i.e. two distinct genomic signatures of positive selection, which respectively assess indirect and direct selection. In brief, we address the following questions: (1) How does recombination rate impact genomic signatures of indirect and direct selection in *Ficedula* flycatchers? (2) How do the dynamics of recombination rate evolution impact both signatures of selection? and (3) How do signatures of indirect and direct selection compare with one another?

## 2. Results

### (a) Study system and whole genome re-sequencing data

To understand the relationship between the evolutionary dynamics of recombination rate and genomic signatures of
<figure><img src="image_9.png" /><figcaption>**Figure 2.** Recombination and differentiation landscapes. (a,b) LD-based recombination landscapes across the genome for the taiga flycatcher (red) and collared flycatcher (blue), and genome-wide variation in $F_{ST}$ between taiga and red-breasted flycatcher (red) and collared and pied flycatcher (blue). Grey rectangles display regions where recombination rate is divergent between taiga flycatcher and collared flycatcher.</figcaption></figure>

recombination ($\rho = 4N_c r$) into cM Mb$^{-1}$ (figure 2a,b) using the linkage map.

We observed a significant correlation for recombination rate between the two species estimated in 200 kb windows (electronic supplementary material, figure S4; Pearson's R = 0.45, p-value < 10$^{-16}$), which was weaker than correlations observed previously between the more closely related collared and pied flycatchers (Pearson's R = 0.79) [47]. We identified that 24% of 200 kb windows showed differences in the recombination rate between the two species, and classified windows into windows with conserved and divergent recombination rate (figure 2a,b). To investigate if the inferred changes in recombination rate represent an evolutionary signal or mainly reflect uncertainty in LD-based recombination rate estimates, we compared the correlation between LD-based recombination rate of both species and pedigree-based recombination rate of collared flycatcher separately for the conserved and divergent windows. We expect an evolutionary signal to manifest in a lower correlation for divergent compared with conserved windows for LD-based recombination rate of taiga flycatcher but not collared flycatcher. By contrast, if a low signal-to-noise ratio explains divergent windows, then a lower correlation with the pedigree-based estimates is expected for both species in these regions. We observed that collared flycatcher showed similar patterns for both window types (Pearson's R = 0.59 and R = 0.53 in conserved and divergent windows respectively; electronic supplementary material, figure S5), but for taiga flycatcher the correlation coefficient was reduced by half in divergent windows compared with conserved windows (Pearson's R = 0.57 and R = 0.24 in conserved and divergent windows respectively; electronic supplementary material, figure S5). Thus, the correlation analysis provides support that inferred evolutionary changes in recombination rate represent an evolutionary signal.

As an additional assessment, we compared signatures of gBGC between the two species separately for conserved and divergent windows. The strength of gBGC should positively correlate with recombination rate [62]; therefore we would expect divergent windows to show a greater impact

of gBGC for the species with the higher recombination rate. Our analysis revealed that, relative to conserved regions, divergent regions with higher recombination rate in collared flycatcher showed a greater impact of gBGC in collared flycatcher than in taiga flycatcher, and *vice versa* for divergent regions with higher recombination rate in taiga flycatcher (electronic supplementary material, figure S6). Because inferred signatures of gBGC rely on patterns in the site frequency spectrum (SFS) rather than patterns of LD, these observations provide an additional and independent validation that the inferred changes in recombination rate reflect an evolutionary signal rather than noise in LD-based measurements.

### (c) Recombination rate shapes genomic signatures of indirect selection

We examined the relationship between recombination rate and measures of indirect selection based on two statistics, the CLR test for selective sweeps and pairwise $F_{ST}$ between the species pairs taiga and red-breasted flycatcher and collared and pied flycatcher (figure 1a). Logistic regression analysis between lineage-specific recombination rate and CLR estimates in both taiga flycatcher and collared flycatcher supported the predicted relationship between recombination and selective sweep signatures (figure 3; electronic supplementary material, table S2 and figure S7), with signatures of selective sweeps tending to occur in windows with lower recombination rate. To investigate how evolutionary changes in recombination rate impact genomic signatures of indirect selection, we also examined correlations across species, which revealed a weaker relationship between selective sweeps estimated in one species and recombination rate estimated in the other (figure 3; electronic supplementary material, table S2 and figure S7).

Next, we compared how lineage-specific estimates of recombination rate corresponded with the differentiation landscape for the two species comparisons. We observed a negative relationship between recombination rate and $F_{ST}$ for both species pairs, taiga and red-breasted flycatchers
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Summary of relationships between genomic signatures of selection and recombination rate. Shown are the *R*² values for four different selection measures and recombination rate. The diameter of the circles indicates the strength of the relationship and the colour represents the direction of the relationship, positive in blue and negative in red. For both species, three comparisons between recombination rate and signatures of selection are shown: lineage-specific recombination rate estimates across all windows, lineage-specific recombination rate estimates across windows with conserved recombination rate, and an interspecies comparison. Asterisks denote significance of the relationship, with \*p-value < 0.05, \*\*p-value < 0.01, and \*\*\*p-value < 0.001. HRI, Hill–Robertson interference.</figcaption></figure>

(figure 3; electronic supplementary material, figure S8; Spearman's R = -0.71; *p*-value < 2.2 × 10⁻¹⁶) and collared and pied flycatchers (figure 3; electronic supplementary material, figure S8; Spearman's R = -0.59; *p*-value < 2.2 × 10⁻¹⁶), which is in good agreement with the linked-selection prediction. As with the selective sweep analysis, the relationship between recombination rate and *F*<sub>ST</sub> was lower for interspecies comparisons (figure 3; electronic supplementary material, figure S8), and slightly higher when restricted to conserved recombination rate regions (figure 3). We also identified *F*<sub>ST</sub> peaks between both species pairs, where local *F*<sub>ST</sub> was significantly higher than the chromosomal background, which we then categorized as lineage-specific or shared between the species comparisons (electronic supplementary material, table S3). Both taiga flycatcher and collared flycatcher showed significantly reduced recombination rate in shared *F*<sub>ST</sub> peaks (electronic supplementary material, figure S9) compared with windows without *F*<sub>ST</sub> peaks. Collared flycatcher showed significantly lower recombination rate in collared/pied-specific (CP unique) *F*<sub>ST</sub> peaks, but not in taiga/red-breasted-specific (TR unique) peaks (electronic supplementary material, figure S9). Although taiga flycatcher showed a significant reduction in recombination rate for both TR unique peaks as well as CP unique peaks compared with the genomic background, the reduction was greater for shared and TR unique peaks compared with CP unique peaks (electronic supplementary material, figure S9). Overall, these results support that evolutionary changes in recombination rate shape lineage-specific signals of *F*<sub>ST</sub> peaks. Note also that for collared flycatcher, variation in LD-based recombination rate among *F*<sub>ST</sub> peaks recaptures earlier observations using the pedigree-based recombination rate ([51]; electronic supplementary material, figure S6b). Thus, our results suggest that neglecting evolutionary changes in recombination rate can provide a biased picture of the action of natural selection.

Previous studies have demonstrated that recombination rate is often reduced towards the centres of chromosomes [47,63], an observation replicated here for the macrochromosomes (figure 2; electronic supplementary material, figure S10), while the opposite pattern was observed on microchromosomes (figure 2; electronic supplementary material, figure S10). We investigated whether there was also an impact of the relative location along the chromosome on the presence of *F*<sub>ST</sub> peaks. For both species, *F*<sub>ST</sub> peaks were enriched in the ends of microchromosomes, where recombination rate was reduced, suggesting that chromosomal location does impact *F*<sub>ST</sub> peak prevalence (electronic supplementary material, figure S10). In macrochromosomes there was little variation among chromosomal regions for collared flycatcher (electronic supplementary material, figure S10); however, *F*<sub>ST</sub> peaks were more prevalent in chromosome ends for taiga flycatcher (electronic supplementary material, figure S10). It is therefore possible that chromosomal rearrangements that change the relative location of regions along the chromosome may contribute to generating some of the observed changes in recombination rate and associated lineage-specific *F*<sub>ST</sub> peaks.

### (d) Genomic signatures of direct selection are not consistent with HRI

To assess evidence for HRI in taiga flycatcher and collared flycatcher, we compared the relationship between recombination and two lineage-specific signatures of direct selection, π<sub>N</sub>/π<sub>S</sub> and *d*<sub>N</sub>/*d*<sub>S</sub>, on a gene-by-gene basis. To account for the noisiness of gene-based estimates of π<sub>N</sub>/π<sub>S</sub> and *d*<sub>N</sub>/*d*<sub>S</sub>, we also assessed evidence for HRI using a binning approach. For this purpose, we assigned genes to one of three bins with varying recombination rate, after estimating the per gene average recombination rate (electronic supplementary material, table S4). The binning approach permits estimation of the distribution of fitness effects (DFE), particularly the mean strength
of selection against deleterious mutations $N_{es}$, and $\omega_a$, for which there is insufficient power on a gene-by-gene basis. Since we might expect that any relationship between recombination rate and measures of direct selection would be weakened by evolutionary changes in recombination rate, we additionally estimated $\pi_N/\pi_S$, $d_N/d_S$, and $\omega_a$ for the subset of genes with conserved gene-based recombination rate between taiga flycatcher and collared flycatcher. We also examined whether variation in the density of functional sites might obscure patterns of HRI, by dividing recombination rate bins into low and high functional density (i.e. density of exons and conserved non-coding elements (CNEs); electronic supplementary material, table S5). Finally, to account for any potential bias due to gBGC, estimates of $\pi_N/\pi_S$, $N_{es}$, $d_N/d_S$ and $\omega_a$ are based on GC-conservative changes only.

Based on the gene-by-gene analysis as well as the binning approach we did not observe the relationships between recombination rate and signatures of direct selection predicted by HRI in either species, regardless of whether we restricted the analysis to genes with conserved recombination rate or not (figure 3; electronic supplementary material, figures S11 and S12). Accounting for variation in functional density, we again observed no consistent pattern between species that matched the predictions based on HRI for recombination rate and signatures of direct selection (electronic supplementary material, figures S13 and S14). Taken together, our results therefore suggest a negligible role of recombination in shaping genomic signatures of direct selection in *Ficedula* flycatchers.

### (e) Weak association between signatures of direct and indirect selection

Finally, we compared genomic signatures of direct and indirect selection with each other to address the hypothesis that recurrent selective sweeps may contribute to $F_{ST}$ peaks, reflected by an increased rate of adaptive substitutions in $F_{ST}$ peaks compared with the genomic background. For this purpose, we compared estimates of $\pi_N/\pi_S$, $N_{es}$, $d_N/d_S$ and $\omega_a$ for genes overlapping with shared $F_{ST}$ peaks, lineage-specific $F_{ST}$ peaks, and $F_{ST}$ peaks with or without a selective sweep signature against genes not overlapping with $F_{ST}$ peaks as a background reference. Lineage-specific $F_{ST}$ peaks and $F_{ST}$ peaks without a selective sweep signature showed no significant differences in signatures of direct selection compared with the genomic background. Both taiga flycatcher and collared flycatcher showed significantly higher $\pi_N/\pi_S$ in shared $F_{ST}$ peaks relative to the genomic background (table 1), and taiga flycatcher showed significantly larger $d_N/d_S$ and $\omega_a$. For collared flycatcher, estimates of $d_N/d_S$ were elevated in shared $F_{ST}$ peaks but differences were not statistically significant. Similarly, $\pi_N/\pi_S$ was significantly higher for both species within $F_{ST}$ peaks overlapping with selective sweeps, while differences in $d_N/d_S$ were only significantly larger for taiga flycatcher, and differences in $\omega_a$ were not statistically significant for either species (table 1).

## 3. Discussion

### (a) Signatures of indirect selection are shaped by recombination rate dynamics

Speciation genomic studies across a wide range of species have revealed that differentiation islands tend to coincide

<table><caption>**Table 1.** Comparison of measures of indirect and direct selection. Shown are estimates of the efficacy of direct selection for genes grouped by $F_{ST}$ peak category. Categories include genes outside $F_{ST}$ peaks, genes in $F_{ST}$ peaks shared between species comparisons, genes within peaks specific to the focal species, genes within peaks that show a signature of a selective sweep, and genes within peaks that do not show a signature of a selective sweep. *p*-values comparing whether statistics within $F_{ST}$ peaks are significantly different from outside $F_{ST}$ peaks are presented in parentheses. Significant differences are indicated in bold, *p*-value < 0.05. See electronic supplementary material, table S6 for estimates based on all changes.</caption><thead><tr><th rowspan="2"></th><th colspan="3">taiga flycatcher</th><th colspan="3">collared flycatcher</th></tr><tr><th>π<sub>N</sub>/π<sub>S</sub></th><th>N<sub>es</sub></th><th>d<sub>N</sub>/d<sub>S</sub></th><th>ω<sub>a</sub></th><th>genes</th><th>N<sub>es</sub></th><th>d<sub>N</sub>/d<sub>S</sub></th><th>ω<sub>a</sub></th><th>genes</th></tr></thead><tbody><tr><td>outside peaks</td><td>0.17</td><td>1485</td><td>0.19</td><td>0.087</td><td>5862</td><td>0.17</td><td>34852</td><td>0.18</td><td>0.029</td><td>5842</td></tr><tr><td>shared peaks</td><td><strong>0.34</strong> (0.001)</td><td>146 (0.16)</td><td><strong>0.32</strong> (0.021)</td><td><strong>0.21</strong> (0.034)</td><td>245</td><td><strong>0.32</strong> (0.002)</td><td>4789 (0.70)</td><td>0.23 (0.25)</td><td>-0.00084 (0.65)</td><td>241</td></tr><tr><td>lineage-specific peaks</td><td>0.12 (0.15)</td><td>149 (0.24)</td><td>0.17 (0.79)</td><td>0.14 (0.37)</td><td>172</td><td>0.086 (0.16)</td><td>208 (0.53)</td><td>0.091 (0.29)</td><td>0.063 (0.71)</td><td>82</td></tr><tr><td>peaks with sweep overlap</td><td><strong>0.46</strong> (0.001)</td><td>111 (0.21)</td><td><strong>0.31</strong> (0.038)</td><td>0.15 (0.26)</td><td>201</td><td><strong>0.53</strong> (0.006)</td><td>8.3 (0.67)</td><td>0.093 (0.43)</td><td><strong>-0.37</strong> (0.041)</td><td>52</td></tr><tr><td>peaks with no sweep overlap</td><td>0.20 (0.39)</td><td>190 (0.21)</td><td>0.21 (0.72)</td><td>0.13 (0.39)</td><td>206</td><td>0.23 (0.10)</td><td>13818 (0.84)</td><td>0.24 (0.22)</td><td>0.073 (0.45)</td><td>265</td></tr></tbody></table>
selection [9,14–16]. Simulation studies suggest a relationship between the rate of selective sweeps and the adaptive substitution rate, at least in the absence of recombination [66]. Since signatures of selective sweeps rapidly break down with increasing recombination rate [12,13], the predicted relationship is less obvious in recombinating genomes. The *Ficedula* flycatcher system allows us to explore this relationship, since the comparison of adaptive substitution rates in lineage-specific and shared $F_{ST}$ peaks permits assessing the role of recombination dynamics. With variable recombination rate, we did not observe a relationship between lineage-specific $F_{ST}$ peaks and any genomic signatures of direct selection. This may be because there has been no increase in adaptive substitutions in these regions or that any adaptive substitutions occurred too recently to be detected with the methods employed. Shared $F_{ST}$ peaks and $F_{ST}$ peaks coinciding with selective sweeps showed significantly higher $\pi_N/\pi_S$ for both species, which was driven by a greater reduction in $\pi_S$ compared with $\pi_N$ as both values were reduced compared with the background. Hitchhiking due to a selective sweep is known to cause a reduction in linked diversity [9], and theory suggests that the reduction in diversity is greater for synonymous sites than nonsynonymous sites [67]. Although estimates of the adaptive substitution rate were significantly elevated in the shared $F_{ST}$ peaks compared with the genomic background for taiga flycatcher but not for collared flycatcher, the adaptive substitution rate was not significantly larger in $F_{ST}$ peaks coinciding with selective sweeps for either of the two species. These results therefore indicate that differentiation islands and selective sweep signatures may not require significantly elevated rates of adaptation in order to manifest in the genome, and highlight again the central role of recombination dynamics in shaping indirect selection signatures.

## (b) Recombination rate is not a major force in shaping signatures of direct selection

In addition to creating more pronounced signatures of indirect selection, tighter linkage between sites owing to low recombination rate is predicted to increase the impact of HRI. The consequence of this relationship is an increase in the efficacy of natural selection with increasing recombination rate, which has been observed in several systems. However, within both taiga flycatcher and collared flycatcher signatures of direct selection generally showed little variation with recombination rate, and did not reflect the relationships predicted by HRI.

There are several explanations why HRI may be less prevalent in *Ficedula* flycatchers compared with other systems. Selective interference occurs when multiple linked, selected mutations segregate within the population at the same time, and it has been demonstrated that the addition of neutral sites between selected sites may help to alleviate selective interference [39]. We are therefore more likely to observe HRI where there is a greater density of functional sequences relative to the recombination rate. For example, gene density per centimorgan in *D. melanogaster*, where evidence for HRI has been observed, is roughly an order of magnitude larger than in flycatcher [68] (FlyBase release FB2022_04; Ensembl v. 107). Additionally, genetic diversity is much higher in *D. melanogaster* compared with flycatchers [69], meaning multiple selected mutations will be more likely to segregate at the same time. These circumstances could explain

with regions of low recombination as a consequence of a local reduction of genetic diversity owing to indirect selection [20–22,26,51]. Despite the central role of recombination in shaping the differentiation landscape, the role of the evolutionary dynamics of recombination rate in these patterns has largely remained unexplored, often because recombination rate estimates have only been available for one of the studied species, or not available for the focal species ([20,22,26,51] but see [47]). Recent work has, however, demonstrated that not accounting for recombination rate evolution can lead to an underestimation of the impact of recombination rate on signatures of indirect selection [64], particularly in species where *PRDM9* entails highly dynamic recombination rates, such as mammals. Birds, in contrast, lack *PRDM9* and their genomes are characterized by a more stable fine-scale recombination landscape, which may lead to the conclusion that estimates of recombination rate from one species are sufficient to explain patterns of selection also at the broad-scale. Nevertheless, by estimating recombination rate in two *Ficedula* flycatcher species we provide evidence that evolutionary changes in recombination rate are central in driving differences in signatures of indirect selection between species, while evolutionarily stable reductions in recombination rate manifest in differentiation islands that are shared among species. Therefore, we show that even in species that lack *PRDM9*, accounting for variation in recombination rate between species is crucial to interpreting genomic signatures of natural selection. We thereby expand on our understanding from previous studies where recombination rate estimates were only available for a single species.

Characteristic for the *Ficedula* flycatchers is that both lineage-specific and shared signatures of indirect selection generally stretch several hundred kilobases, suggesting that the mechanism behind low-recombining regions acts at the broad-scale rather than the fine-scale. One possible mechanism for broad-scale reductions in recombination rate that could also generate changes in recombination rate between species involves chromosomal rearrangements such as inversions or translocations [48–50]. Inversions can lead to a reduction in recombination when polymorphic owing to the absence of crossing over in heterozygotes, while translocations can relocate chromosomal segments to a location where recombination rate differs. For example, reduced recombination is often observed in chromosome centres [63], and we observed that chromosomal location was associated with $F_{ST}$ peak prevalence, with differences observed between the two species. Aside from chromosomal rearrangements, centromeres are also suggested to contribute to broad-scale $F_{ST}$ peaks in flycatchers [20]. Previous work has shown that centromere shifts occur in birds [65], which could also explain some of the changes in recombination rate between taiga flycatcher and collared flycatcher. Nevertheless, in several instances we find multiple $F_{ST}$ peaks per chromosome, a signature not expected if $F_{ST}$ peaks solely coincide with centromeres but pointing towards another mechanism, such as chromosomal rearrangements. While we can currently only speculate as to the molecular mechanisms behind the recombination landscape dynamics between these species, future research into these hypotheses can benefit from recent advancements in long-read sequencing technology.

Yet, regardless of the molecular mechanism that drives local reductions in recombination rate, differentiation islands or selective sweeps would not arise without the action of
why we find no effect of recombination rate on selection effi-
cacy in the *Ficedula* flycatchers, even in the most functionally
dense regions of the genome. It is important to note, however,
that evidence of HRI in birds is inconclusive [32,34,37], and
generalization among birds asks for further investigations.

An absence of HRI in flycatchers is nevertheless in good
agreement with earlier observations, which suggest that the
intensity of selection on coding sequences is impacted more
by functional characteristics such as gene expression patterns
than by the genomic background [70]. For instance, recent
work found that estimates of *d*<sub>N</sub>/<i>d</i><sub>S</sub> were highly correlated
among distantly related species with divergent genomic
backgrounds [71], suggesting that conserved gene functions
rather than genomic background shape the variation in *d*<sub>N</sub>/
<i>d</i><sub>S</sub> among genes. To investigate if evolution of the recombina-
tion landscape could obscure the expected pattern of HRI in
flycatchers, we limited our analysis to genes with evolutiona-
rily conserved recombination rate between taiga flycatcher
and collared flycatcher. However, we still observed no clear
relationship between recombination rate and signatures of
direct selection consistent with HRI. This suggests that
recombination rate and its dynamics appear to not play a sig-
nificant role in shaping any differences in signatures of direct
selection between the two birds studied here.

(c) Conclusions

By comparing recombination dynamics with genomic signa-
tures of indirect and direct selection, we show that these
different signatures are not shaped by the same evolutionary
processes. While signatures of indirect selection appear to be
strongly shaped by the recombination landscape and its
dynamics, signatures of direct selection are largely unaffected
by the genomic background. Rather, gene function and
expression patterns may play a more central role in shaping
the efficacy of direct selection on protein-coding sequence [70].

4. Material and methods

(a) Samples and genotyping

Whole genome re-sequencing data for 187 flycatchers were col-
lated from previously published work, including 65 taiga
flycatchers [51], 15 red-breasted flycatchers [51], 95 collared fly-
catchers [52], 11 pied flycatchers [20], and one sample of
snowy-browed flycatcher as an outgroup [20]. Base quality
score recalibrated BAM files with reads mapped to the collared
flycatcher reference genome FicAlb1.5 [58] were obtained for
all samples. Genotyping was performed individually with Hap-
lotypeCaller in GATK v.4.1, followed by GenotypeGVCFs with
all samples combined specifying the *-all-sites* flag to genotype
both variant and invariant sites. We removed indels and variable
sites with more than one alternative allele using VCFtools
v.0.1.16 [72], resulting in a dataset containing 974 084 046 sites
in total, including 116 929 990 SNVs.

Here, we focus only on sites located on the autosomes, exclud-
ing unassigned scaffolds, the Z-chromosome, and mitochondria.
We followed the filtering methods applied in [50], resulting in a
final dataset of 566 724 393 callable sites including 51 424 863
SNVs. We polarized variable sites using snowy-browed flycatcher
as an outgroup, and combining taiga and red-breasted flycatcher
into one group and collared and pied flycatcher into another
group. Whenever any two of the three groups were fixed for the
same allele, this allele was considered the ancestral state (electronic
supplementary material, methods). This resulted in 564 393 274
callable sites and 49 121 805 polarized SNVs.

(b) Demographic modelling

We investigated historical fluctuations in population size for the four ingroup flycatcher species using the pairwise sequentially Markovian coalescent (PSMC; [54]) with one individual per species. Sites with a read depth less than 10, sites with more than twice the genome-wide average read depth, and 100 bp blocks with more than 20% missing data were excluded. We used parameters following [73] and performed 100 bootstrap replicates. In addition, we examined the demographic history of the two species pairs with the method of ordinary differential equations using moments software [74] implemented in GADMA [53] with the unfolded joint site frequency spectrum (jSFS) computed from *daði* [75]. We used the following parameters: 100 replicates, *μ* = 4.6 × 10⁻⁹, *g* = 2, effective length *L* = 564 393 274, and *θ*₀ = 4 µl = 10.38. The best model was selected according to both log-likelihood and composite likelihood Akaike information criterion (CLAIC). Confidence intervals were calculated for parameter estimates of the best fitting model by dividing each dataset into 1 Mb segments and resampling from these segments for bootstrapping.

(c) Estimation of recombination rate

Population-scaled recombination rate (ρ = 4N<sub>e</sub>r) was estimated for
taiga flycatcher and collared flycatcher. Haplotypes were inferred
with SHAPEIT2 [76]. We set genome-wide ρ to 0.037 [47], the
effective population size to 300 000 [77], the *-states* parameter to
200 and *-window* to 0.5. Recombination rate was estimated with
LDhelmet [78] using 25 samples with the least missing data (elec-
tronic supplementary material, table S7) and following parameters
in [44] (electronic supplementary material, methods). We calcu-
lated the average LD recombination rate in 200 kb, 1 Mb, and
5 Mb genomic windows, weighted by the distance between SNP
pairs accounting for SNP pairs spanning adjacent windows.
Additionally, we summarized the weighted average of recombina-
tion rate within each protein-coding gene.

Population-scaled recombination rates were converted into
cM Mb⁻¹ by dividing by the genome-wide N<sub>e</sub> following [48],
using the pedigree recombination map of collared flycatcher, in
200 kb, 1-Mb, and 5 Mb windows [58]. This gave an N<sub>e</sub> estimate
of 318 000–458 000 for taiga flycatcher, and 30 000–42 500 for the
Gotland population of collared flycatcher, where the latter is in
line with LD-based N<sub>e</sub> estimates from [52]. We used the value of
N<sub>e</sub> estimated from 5 Mb windows, which showed the strongest
correlation with the linkage map (electronic supplementary
material, figure S3), to convert estimates of recombination rate
into cM Mb⁻¹. For both the window-based and gene-based recom-
bination rates, we divided regions into 'conserved' or 'divergent'
recombination rate between the two species. We performed a prin-
cipal component analysis (PCA) using the recombination rate
estimates for taiga flycatcher and collared flycatcher and identified
regions with conserved recombination rate if they were within one
standard deviation from the mean of the second principal com-
ponent. We validated our inference of changes in recombination
rate by comparing the strength of gBGC in regions with conserved
and divergent recombination rates between species (electronic
supplementary material, methods; [47]).

(d) Estimation of signatures of indirect selection within
and between species

To assess signatures of indirect selection, we performed a selec-
tive sweep scan within taiga flycatcher and collared flycatcher.
A composite likelihood ratio (CLR) test [79] was used to identify
selective sweeps, using polarized SNP data (see §4a) for the two
species. We ran SweepFinder2 [80] with a pre-computed back-
ground SFS from each species (taiga flycatcher *n* = 65 diploid
individuals; collared flycatcher *n* = 95), and a grid of locations
for every variant. A significance threshold of 46.25 was taken from previous work [51], based on simulations using the collared flycatcher annotation [81] and recombination map [58] to obtain a threshold that accounts for BGS. We merged adjacent sites with significant CLR values, removed sweeps with only one significant site or with a site density less than $1 \text{ kb}^{-1}$, and called presence/absence of selective sweeps in 200 kb windows.

We calculated Weir & Cockerham $F_{ST}$ [82] between the species pairs taiga and red-breasted flycatcher and collared and pied flycatcher in 200 kb nonoverlapping windows using VCFtools v.0.1.16 [72]. To identify $F_{ST}$ peaks, we Z-transformed $F_{ST}$ separately for each chromosome to account for differences in mean $F_{ST}$ among chromosomes, and applied a Savitzky-Golay smoothing filter to the $ZF_{ST}$ values using a polynomial of 3 and a filter length of 7. We identified $F_{ST}$ peaks where $ZF_{ST}$ was two standard deviations above the chromosome mean, and then classified these as peaks into shared and lineage-specific, i.e. taiga-red-breasted unique (TR unique) or as collared-pied unique (CP unique), based on any overlap between the two species comparisons. To compare recombination rate in different $F_{ST}$ peak categories, we performed a permutation test, randomizing the peak category and re-estimating the difference in mean recombination rate among 200 kb windows for 1000 permutations.

### (e) Estimation of signatures of direct selection within and between species

To estimate signatures of direct selection on protein-coding sequences we constructed sequence alignments for one-to-one orthologues between taiga flycatcher, collared flycatcher and zebra finch (*Taeniopygia guttata*). Zebra finch sequences were retrieved from Ensembl v.104 for one-to-one orthologues between collared flycatcher and zebra finch. We reconstructed gene sequences for multiple sample sizes of taiga flycatcher and collared flycatcher, to account for the recent divergence time of the two species (approx. 3.76 coalescent units) [51]; see supplementary material, methods for details. For each gene we combined the sequences for all samples into a single fasta file and aligned the sequences using prank v.170427 [83]. Prior to alignment with prank, we ran clustalw v.2.1 [84] to estimate a guide tree for each gene, which we then provided as input to prank, running prank with the *-once* flag to perform one iteration. We estimated signatures of direct selection in taiga flycatcher and collared flycatcher for one-to-one orthologues between flycatcher and zebra finch. Using fourfold and zerofold degenerate sites extracted from the reconstructed ancestral genome, we calculated $\pi_N/\pi_S$ for taiga flycatcher ($n=65$) and collared flycatcher ($n=95$). We subset the variable sites for only GC-conservative variants, i.e. strong-to-strong (C and G) variants or weak-to-weak (A and T) variants, to correct for gBGC.

Using Bio++ v.3.0 [85] we estimated $d_N/d_S$ for taiga flycatcher and collared flycatcher based on a strand-symmetric L95 model [86] implemented in bppml, and mapped substitutions as implemented in mapnh [87]. This allowed us to separate substitutions into different categories and correct for gBGC by restricting the analysis to GC-conservative changes. We excluded genes with an exon length below 200 bp and summarized $d_N/d_S$ across genes by dividing the average value of $d_N$

by the average value of $d_S$. There was little evidence of sample size dependence in the estimates of $d_N/d_S$ (electronic supplementary material, tables S8 and S9; electronic supplementary material, methods), suggesting that divergence between taiga flycatcher and collared flycatcher is deep enough for phylogenetic analysis [88]. All further analyses with $d_N/d_S$ were therefore only based on the 32 haplotype sample size for both taiga flycatcher and collared flycatcher.

Using the fourfold and zerofold site frequency spectra and the results of $d_N/d_S$, we estimated the adaptive substitution rate ($\omega_a$) with the software DFE-alpha [89,90] for bins of genes (electronic supplementary material, tables S4 and S5). The binning approach is outlined in supplementary material, methods. We used the 2-epoch model of demographic history for all analyses, after confirming this was the best fitting model based on all genes (electronic supplementary material, tables S10 and S11). To test for statistically significant differences in measures of selection between gene bins we performed 1000 permutations, reshuffling the genes among bins. We used the same approach to test for differences in measures estimated for different $F_{ST}$ peaks compared with the genomic background.

**Ethics.** We use animal genetic data from earlier publications that specify according ethical approvals.

**Data accessibility.** Sequencing data for all samples are available at the EMBL-EBI European Nucleotide Archive (ENA; http://www.ebi.ac.uk/ena) with the following accession numbers: PRJEB43825 (taiga and red-breasted flycatchers), PRJEB22864 (collared flycatcher), and PRJEB7359 (pied and snowy-browed flycatchers). VCF files used for analyses and recombination rate data are available on Dryad (https://doi.org/10.5061/dryad.q2bvq83nw) [91]. Scripts used for analysis are available at https://github.com/madeline-chase/flycatcher_recom.

Supplementary material is available online [92].

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** M.A.C.: conceptualization, data curation, formal analysis, investigation, validation, visualization, writing—original draft, writing—review and editing; M.V.: formal analysis, investigation, visualization; C.F.M.: conceptualization, funding acquisition, investigation, project administration, supervision, validation, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed herein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** This study was funded by grants from the Swedish Research Council (2013-8271 to Hans Ellegren), and the Knut and Alice Wallenberg Foundation (2014/0044 to Hans Ellegren). The computational resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) were partially funded by the Swedish Research Council through grant agreements no. 2022-06725 and no. 2018-05973.

**Acknowledgements.** We thank Laurent Guéguen for generous advice on running Bio++, and Laurent Duret for helpful feedback on an earlier version of this manuscript. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX).

## References

1. Chase MA, Mugal CF. 2022 The role of recombination dynamics in shaping signatures of direct and indirect selection

across the *Ficedula* flycatcher genome.
*bioRxiv*, 503468. (doi:10.1101/2022.08.11.503468)

2. Barton NH, Charlesworth B. 1998 Why sex and recombination? *Science* **281**, 1986–1990. (doi:10.1126/science.281.5385.1986)
3. Felsenstein J. 1974 The evolutionary advantage of recombination. *Genetics* **78**, 737–756. (doi:10.1093/ genetics/78.2.737)

4. Kirubakaran TG *et al*. 2016 Two adjacent inversions maintain genomic differentiation between migratory and stationary ecotypes of Atlantic cod. *Mol. Ecol*. **25**, 2130–2143. (doi:10.1111/ mec.13592)

5. Lowry DB, Willis JH. 2010 A widespread chromosomal inversion polymorphism contributes to a major life-history transition, local adaptation, and reproductive isolation. *PLoS Biol*. **8**, e1000500. (doi:10.1371/journal.pbio.1000500)

6. Glémin S. 2010 Surprising fitness consequences of GC-biased gene conversion: I. Mutation load and inbreeding depression. *Genetics* **185**, 939–959. (doi:10.1534/genetics.110.116368)

7. Xue C, Chen H, Yu F. 2016 Base-biased evolution of disease-associated mutations in the human genome. *Hum. Mutat*. **37**, 1209–1214. (doi:10. 1002/humu.23065)

8. Hill WG, Robertson A. 1966 The effect of linkage on limits to artificial selection. *Genet. Res*. **8**, 269–294. (doi:10.1017/S0016672300010156)

9. Smith JM, Haigh J. 1974 The hitch-hiking effect of a favourable gene. *Genet. Res*. **23**, 23–35. (doi:10. 1017/S0016672300014634)

10. Stephan W. 2019 Selective sweeps. *Genetics* **211**, 5–13. (doi:10.1534/genetics.118.301319)

11. Murphy DA, Elyashiv E, Amster G, Sella G. 2022 Broad-scale variation in human genetic diversity levels is predicted by purifying selection on coding and non-coding elements. *eLife* **11**, e76065. (doi:10.7554/eLife.76065)

12. Kaplan NL, Hudson RR, Langley CH. 1989 The 'hitchhiking effect' revisited. *Genetics* **123**, 887–899. (doi:10.1093/genetics/123.4.887)

13. Ohta T, Kimura M. 1975 The effect of selected linked locus on heterozygosity of neutral alleles (the hitch-hiking effect). *Genet. Res*. **25**, 313–325. (doi:10.1017/ S0016672300015731)

14. Charlesworth B, Morgan MT, Charlesworth D. 1993 The effect of deleterious mutations on neutral molecular variation. *Genetics* **134**, 1289–1303. (doi:10.1093/genetics/134.4.1289)

15. Cruickshank TE, Hahn MW. 2014 Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. *Mol. Ecol*. **23**, 3133–3157. (doi:10.1111/mec.12796)

16. Wolf JBW, Ellegren H. 2017 Making sense of genomic islands of differentiation in light of speciation. *Nat. Rev. Genet*. **18**, 87–100. (doi:10. 1038/nrg.2016.133)

17. Begun DJ, Aquadro CF. 1992 Levels of naturally occurring DNA polymorphism correlate with recombination rates in *D. melanogaster*. *Nature* **356**, 519–520. (doi:10.1038/356519a0)

18. Corbett-Detig RB, Hartl DL, Sackton TB. 2015 Natural selection constrains neutral diversity across a wide range of species. *PLoS Biol*. **13**, e1002112. (doi:10.1371/journal.pbio. 1002112)

19. Cutter AD, Payseur BA. 2013 Genomic signatures of selection at linked sites: unifying the disparity among species. *Nat. Rev. Genet*. **14**, 262–274. (doi:10.1038/nrg3425)

20. Burri R *et al*. 2015 Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of *Ficedula* flycatchers. *Genome Res*. **25**, 1656–1665. (doi:10.1101/gr. 196485.115)

21. Renaut S, Grassa CJ, Yeaman S, Moyers BT, Lai Z, Kane NC, Bowers JF, Burke JM, Rieseberg LH. 2013 Genomic islands of divergence are not affected by geography of speciation in sunflowers. *Nat. Commun*. **4**, 1827. (doi:10.1038/ncomms2833)

22. Doren BMV, Campagna L, Helm B, Illera JC, Lovette IJ, Liedvogel M. 2017 Correlated patterns of genetic diversity and differentiation across an avian family. *Mol. Ecol*. **26**, 3982–3997. (doi:10.1111/mec.14083)

23. Delmore KE, Lugo Ramos JS, Van Doren BM, Lundberg M, Bensch S, Irwin DE, Liedvogel M. 2018 Comparative analysis examining patterns of genomic differentiation across multiple episodes of population divergence in birds. *Evol. Lett*. **2**, 76–87. (doi:10.1002/evl3.46)

24. Martin SH, Davey JW, Salazar C, Jiggins CD. 2019 Recombination rate variation shapes barriers to introgression across butterfly genomes. *PLoS Biol*. **17**, e2006288. (doi:10.1371/journal.pbio.2006288)

25. Schumer M *et al*. 2018 Natural selection interacts with recombination to shape the evolution of hybrid genomes. *Science* **360**, 656–660. (doi:10. 1126/science.aar3684)

26. Stankowski S, Chase MA, Fuiten AM, Rodrigues MF, Ralph PL, Streisfeld MA. 2019 Widespread selection and gene flow shape the genomic landscape during a radiation of monkeyflowers. *PLoS Biol*. **17**, e3000391. (doi:10.1371/journal.pbio.3000391)

27. Roesti M, Hendry AP, Salzburger W, Berner D. 2012 Genome divergence during evolutionary diversification as revealed in replicate lake–stream stickleback population pairs. *Mol. Ecol*. **21**, 2852–2862. (doi:10.1111/j.1365-294X.2012.05509.x)

28. Comeron JM, Kreitman M, Aguadé M. 1999 Natural selection on synonymous sites is correlated with gene length and recombination in *Drosophila*. *Genetics* **151**, 239–249. (doi:10.1093/genetics/151. 1.239)

29. Comeron JM, Kreitman M. 2000 The correlation between intron length and recombination in *Drosophila*: dynamic equilibrium between mutational and selective forces. *Genetics* **156**, 1175–1190. (doi:10.1093/genetics/156.3.1175)

30. McGaugh SE, Heil CSS, Manzano-Winkler B, Loewe L, Goldstein S, Himmel TL, Noor MAF. 2012 Recombination modulates how selection affects linked sites in *Drosophila*. *PLoS Biol* **10**, e1001422. (doi:10.1371/journal.pbio.1001422)

31. Martin SH, Möst M, Palmer WJ, Salazar C, McMillan WO, Jiggins FM, Jiggins CD. 2016 Natural selection and genetic diversity in the butterfly *Heliconius melpomene*. *Genetics* **203**, 525–541. (doi:10.1534/ genetics.115.183285)

32. Gossmann Tl, Santure AW, Sheldon BC, Slate J, Zeng K. 2014 Highly variable recombinational landscape modulates efficacy of natural selection in birds. *Genome Biol. Evol*. **6**, 2061–2075. (doi:10.1093/gbe/ evu157)

33. Hough J, Wang W, Barrett SCH, Wright SI. 2017 Hill-Robertson interference reduces genetic diversity on a young plant Y-chromosome. *Genetics* **207**, 685–695. (doi:10.1534/genetics.117.300142)

34. Bolívar P, Mugal CF, Nater A, Ellegren H. 2016 Recombination rate variation modulates gene sequence evolution mainly via GC-biased gene conversion, not Hill-Robertson interference, in an avian system. *Mol. Biol. Evol*. **33**, 216–227. (doi:10. 1093/molbev/msv214)

35. Bullaughey K, Przeworski M, Coop G. 2008 No effect of recombination on the efficacy of natural selection in primates. *Genome Res*. **18**, 544–554. (doi:10. 1101/gr.071548.107)

36. Guillén Y, Casillas S, Ruiz A. 2019 Genome-wide patterns of sequence divergence of protein-coding genes between *Drosophila buzzatii* and *D. mojavensis*. *J. Hered*. **110**, 92–101. (doi:10.1093/ jhered/esy041)

37. Ponnikas S, Sigeman H, Lundberg M, Hansson B. 2022 Extreme variation in recombination rate and genetic diversity along the Sylvioidea neo-sex chromosome. *Mol. Ecol*. **31**, 3566–3583. (doi:10. 1111/mec.16532)

38. Comeron JM, Kreitman M. 2002 Population, evolutionary and genomic consequences of interference selection. *Genetics* **161**, 389–410. (doi:10.1093/genetics/161.1.389)

39. Comeron JM, Williford A, Kliman RM. 2008 The Hill-Robertson effect: evolutionary consequences of weak selection and linkage in finite populations. *Heredity* **100**, 19–31. (doi:10.1038/sj.hdy.6801059)

40. Parvanov ED, Petkov PM, Paigen K. 2010 *Prdm9* controls activation of mammalian recombination hotspots. *Science* **327**, 835. (doi:10.1126/science. 1181495)

41. Baudat F, Buard J, Grey C, Fledel-Alon A, Ober C, Przeworski M, Coop G, De Massy B. 2010 PRDM9 is a major determinant of meiotic recombination hotspots in humans and mice. *Science* **327**, 836–840. (doi:10.1126/science.1183439)

42. Myers S, Bowden R, Tumian A, Bontrop RE, Freeman C, Macfie TS, Mcvean G, Donnelly P. 2010 Drive against hotspot motifs in primates implicates the *PRDM9* gene in meiotic recombination. *Science* **327**, 876–879. (doi:10.1126/science.1182363)

43. Ponting CP. 2011 What are the genomic drivers of the rapid evolution of *PRDM9*? *Trends Genet*. **27**, 165–171. (doi:10.1016/j.tig.2011.02.001)

44. Ptak SE, Hinds DA, Koehler K, Nickel B, Patil N, Ballinger DG, Przeworski M, Frazer KA, Pääbo S. 2005 Fine-scale recombination patterns differ between chimpanzees and humans. *Nat. Genet*. **37**, 429–434. (doi:10.1038/ng1529)

45. Winckler W *et al*. 2005 Comparison of fine-scale recombination rates in humans and chimpanzees. *Science* **308**, 107–111. (doi:10.1126/science. 1105322)
46. Singhal S *et al.* 2015 Stable recombination hotspots in birds. *Science* **350**, 928–932. (doi:10.1126/science.aad0843)

47. Kawakami T, Mugal CF, Suh A, Nater A, Burri R, Smeds L, Ellegren H. 2017 Whole-genome patterns of linkage disequilibrium across flycatcher populations clarify the causes and consequences of fine-scale recombination rate variation in birds. *Mol. Ecol.* **26**, 4158–4172. (doi:10.1111/mec.14197)

48. Auton *A et al.* 2012 A fine-scale chimpanzee genetic map from population sequencing. *Science* **336**, 193–198. (doi:10.1126/science.1216872)

49. Mugal CF, Arndt PF, Ellegren H. 2013 Twisted signatures of GC-biased gene conversion embedded in an evolutionary stable karyotype. *Mol. Biol. Evol.* **30**, 1700–1712. (doi:10.1093/molbev/mst067)

50. Rieseberg LH. 2001 Chromosomal rearrangements and speciation. *Trends Ecol. Evol.* **16**, 351–358. (doi:10.1016/S0169-5347(01)02187-5)

51. Chase MA, Ellegren H, Mugal CF. 2021 Positive selection plays a major role in shaping signatures of differentiation across the genomic landscape of two independent *Ficedula* flycatcher species pairs. *Evolution* **75**, 2179–2196. (doi:10.1111/evo.14234)

52. Nadachowska-Brzyska K, Dutoit L, Smeds L, Kardos M, Gustafsson L, Ellegren H. 2021 Genomic inference of contemporary effective population size in a large island population of collared flycatchers (*Ficedula albicollis*). *Mol. Ecol.* **30**, 3965–3973. (doi:10.1111/mec.16025)

53. Noskova E, Ulyantsev V, Koepfli KP, O'Brien SJ, Dobrynin P. 2020 GADMA: genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data. *GigaScience* **9**, giaa005. (doi:10.1093/gigascience/giaa005)

54. Li H, Durbin R. 2011 Inference of human population history from individual whole-genome sequences. *Nature* **475**, 493–496. (doi:10.1038/nature10231)

55. Dapper AL, Payseur BA. 2018 Effects of demographic history on the detection of recombination hotspots from linkage disequilibrium. *Mol. Biol. Evol.* **35**, 335–353. (doi:10.1093/molbev/msx272)

56. Raynaud M, Gagnaire PA, Galtier N. 2023 Performance and limitations of linkage-disequilibrium-based methods for inferring the genomic landscape of recombination and detecting hotspots: a simulation study. *Peer Community J.* **3**, e27. (doi:10.24072/pcjournal.254)

57. Samuk K, Noor MAF. 2022 Gene flow biases population genetic inference of recombination rate. *G3* **12**, jkac236. (doi:10.1093/g3journal/jkac236)

58. Kawakami T, Smeds L, Backström N, Husby A, Qvarnström A, Mugal CF, Olason P, Ellegren H. 2014 A high-density linkage map enables a second-generation collared flycatcher genome assembly and reveals the patterns of avian recombination rate

variation and chromosomal evolution. *Mol. Ecol.* **23**, 4035–4058. (doi:10.1111/mec.12810)

59. Auton *A et al.* 2013 Genetic recombination is targeted towards gene promoter regions in dogs. *PLoS Genet.* **9**, e1003984. (doi:10.1371/journal.pgen.1003984)

60. Myers S, Bottolo L, Freeman C, Mcvean G, Donnelly P. 2005 A fine-scale map of recombination rates and hotspots across the human genome. *Science* **310**, 321–324. (doi:10.1126/science.1117196)

61. Bhérer C, Campbell CL, Auton A. 2017 Refined genetic maps reveal sexual dimorphism in human meiotic recombination at multiple scales. *Nat. Commun.* **8**, 14994. (doi:10.1038/ncomms14994)

62. Duret L, Galtier N. 2009 Biased gene conversion and the evolution of mammalian genomic landscapes. *Annu. Rev. Genom. Hum. Genet.* **10**, 285–311. (doi:10.1146/annurev-genom-082908-150001)

63. Haenel Q, Laurentino TG, Roesti M, Berner D. 2018 Meta-analysis of chromosome-scale crossover rate variation in eukaryotes and its significance to evolutionary genomics. *Mol. Ecol.* **27**, 2477–2497. (doi:10.1111/mec.14699)

64. Booker TR, Payseur BA, Tigano A. 2022 Background selection under evolving recombination rates. *Proc. R. Soc. B* **289**, 20220782. (doi:10.1098/rspb.2022.0782)

65. Kiazim LG, O'Connor RE, Larkin DM, Romanov MN, Narushin VG, Brazhnik EA, Griffin DK. 2021 Comparative mapping of the macrochromosomes of eight avian species provides further insight into their phylogenetic relationships and avian karyotype evolution. *Cells* **10**, 362. (doi:10.3390/cells10020362)

66. Lourenço JM, Glémin S, Galtier N. 2013 The rate of molecular adaptation in a changing environment. *Mol. Biol. Evol.* **30**, 1292–1301. (doi:10.1093/molbev/mst026)

67. Otto SP, Whitlock MC. 1997 The probability of fixation in populations of changing size. *Genetics* **146**, 723–733. (doi:10.1093/genetics/146.2.723)

68. Gramates LS *et al.* 2022 FlyBase: a guided tour of highlighted features. *Genetics* **220**, iyac035. (doi:10.1093/genetics/iyac035)

69. Andolfatto P. 2001 Contrasting patterns of X-linked and autosomal nucleotide variation in *Drosophila melanogaster* and *Drosophila simulans*. *Mol. Biol. Evol.* **18**, 279–290. (doi:10.1093/oxfordjournals.molbev.a003804)

70. BoláVar P, Mugal CF, Rossi M, Nater A, Wang M, Dutoit L, Parsch J. 2018 Biased inference of selection due to GC-biased gene conversion and the rate of protein evolution in flycatchers when accounting for it. *Mol. Biol. Evol.* **35**, 2475–2486. (doi:10.1093/molbev/msy149)

71. Sigeman *H et al.* 2021 Avian neo-sex chromosomes reveal dynamics of recombination suppression and W degeneration. *Mol. Biol. Evol.* **38**, 5275–5291. (doi:10.1093/molbev/msab277)

72. Danecek P *et al.* 2011 The variant call format and VCFtools. *Bioinformatics* **27**, 2156–2158. (doi:10.1093/bioinformatics/btr330)

73. Nadachowska-Brzyska K, Burri R, Smeds L, Ellegren H. 2016 PSMC analysis of effective population sizes in molecular ecology and its application to black-and-white *Ficedula* flycatchers. *Mol. Ecol.* **25**, 1058–1072. (doi:10.1111/mec.13540)

74. Jouganous J, Long W, Ragsdale AP, Gravel S. 2017 Inferring the joint demographic history of multiple populations: beyond the diffusion approximation. *Genetics* **206**, 1549–1567. (doi:10.1534/genetics.117.200493)

75. Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD. 2009 Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. *PLoS Genet.* **5**, e1000695. (doi:10.1371/journal.pgen.1000695)

76. Delaneau O, Marchini J, Zagury JF. 2012 A linear complexity phasing method for thousands of genomes. *Nat. Methods* **9**, 179–181. (doi:10.1038/nmeth.1785)

77. Nater A, Burri R, Kawakami T, Smeds L, Ellegren H. 2015 Resolving evolutionary relationships in closely related species with whole-genome sequencing data. *Syst. Biol.* **64**, 1000–1017. (doi:10.1093/sysbio/syv045)

78. Chan AH, Jenkins PA, Song YS. 2012 Genome-wide fine-scale recombination rate variation in *Drosophila melanogaster*. *PLoS Genet.* **8**, e1003090. (doi:10.1371/journal.pgen.1003090)

79. Nielsen R, Williamson S, Kim Y, Hubisz MJ, Clark AG, Bustamante C. 2005 Genomic scans for selective sweeps using SNP data. *Genome Res.* **15**, 1566–1575. (doi:10.1101/gr.4252305)

80. Degiorgio M, Huber CD, Hubisz MJ, Hellmann I, Nielsen R. 2016 SweepFinder2: increased sensitivity, robustness and flexibility. *Bioinformatics* **32**, 1895–1897. (doi:10.1093/bioinformatics/btw051)

81. Uebbing *S et al.* 2016 Divergence in gene expression within and between two closely related flycatcher species. *Mol. Ecol.* **25**, 2015–2028. (doi:10.1111/mec.13596)

82. Weir BS, Cockerham CC. 1984 Estimating F-statistics for the analysis of population structure. *Evolution* **38**, 1358–1370. (doi:10.2307/2408641)

83. Löytynoja A. 2014 Phylogeny-aware alignment with PRANK. In *Multiple sequence alignment methods* (ed. DJ Russell), pp. 155–170. Totowa, NJ: Humana Press.

84. Larkin MA *et al.* 2007 ClustalW and ClustalX version 2.0. *Bioinformatics* **23**, 2947–2948. (doi:10.1093/bioinformatics/btm404)

85. Dutheil J, Boussau B. 2008 Non-homogeneous models of sequence evolution in the Bio++ suite of libraries and programs. *BMC Evol. Biol.* **8**, 255. (doi:10.1186/1471-2148-8-255)

86. Lobry JR. 1995 Properties of a general model of DNA evolution under no-strand-bias conditions. *J. Mol. Evol.* **40**, 326–330. (doi:10.1007/BF00163237)

87. Bolívar P, Guéguen L, Duret L, Ellegren H, Mugal CF. 2019 GC-biased gene conversion conceals the prediction of the nearly neutral theory in avian genomes. *Genome Biol.* **20**, 5. (doi:10.1186/s13059-018-1613-z)
88. Mugal CF, Kutschera VE, Botero-Castro F, Wolf JBW, Kaj I. 2020 Polymorphism data assist estimation of the nonsynonymous over synonymous fixation rate ratio $\omega$ for closely related species. *Mol. Biol. Evol.* **37**, 260–279. (doi:10.1093/molbev/msz203)

89. Keightley PD, Eyre-Walker A. 2007 Joint inference of the distribution of fitness effects of deleterious mutations and population demography based on

nucleotide polymorphism frequencies. *Genetics* **177**, 2251–2261. (doi:10.1534/genetics.107.080663)

90. Eyre-Walker A, Keightley PD. 2009 Estimating the rate of adaptive molecular evolution in the presence of slightly deleterious mutations and population size change. *Mol. Biol. Evol.* **26**, 2097–2108. (doi:10.1093/molbev/msp119)

91. Chase MA, Vilcot M, Mugal CF. 2024 Data from: The role of recombination dynamics in

shaping signatures of direct and indirect selection across the *Ficedula* flycatcher genome. Dryad Digital Repository. (doi:10.5061/dryad.q2bvq83nw)

92. Chase MA, Vilcot M, Mugal CF. 2024 Data from: The role of recombination dynamics in shaping signatures of direct and indirect selection across the *Ficedula* flycatcher genome. Figshare. (doi:10.6084/m9.figshare.c.7005761)