Population-genomic inference of the strength and timing of selection against gene flow

Here, we exploit the genome-wide variation in recombination rate and its effect on genetic divergence. Divergent selection reduces effective gene flow at neutral sites, and this effect decreases with the recombinational distance from the loci under selection. We conceptualize this relationship in terms of the effective migration rate and the expected pairwise between-population coalescence time (Fig. 1A). The latter directly connects to the absolute genetic diversity between populations, a quantity that is readily estimated from DNA sequence data. Our model considers two populations of diploids with effective sizes $N_1$ and $N_2$ and nonoverlapping generations. In population 1, a balance between one-way gene flow from population 2 at rate $m$ per generation and local directional selection is maintained for $\tau$ generations before the present. In this “migration–selection” (MS) phase (Fig. 1A), selection against maladaptive immigrant alleles acts at an arbitrary number of biallelic loci that we call “MS polymorphisms” (MSPs). At each MSP, one allele is favored in population 1 over the other by an average selection coefficient $s$, whereas the deleterious allele is introduced by gene flow. We assume additive fitness and no dominance.

Before the MS phase, we assume a “panmictic” (P) phase in an ancestral population of effective size $N_0$ that starts $\tau$ generations ago and extends into the past (Fig. 1A). We call this the (MS)P demographic scenario. The P phase can be exchanged for an ancestral “migration” (M) phase with gene flow at rate $m_0$. Here, we use the (MS)P and (MS)M scenarios to describe our approach. We provide extensions to more general scenarios with an intermediate “isolation” (I) phase in SI Appendix, part 1.

We denote the per-base-pair recombination rate by $r_{\mathrm{bp}}$ and assume that the MSPs occur at a constant rate, $\nu$, per base pair, such that the distance between consecutive MSPs is exponentially distributed with mean $1/\nu$ base pairs.

Selection against maladapted immigrant alleles acts as a barrier to gene flow in the MS phase. At a focal neutral site, the baseline migration rate $m$ is reduced to an effective migration rate $m_e$ (19, 20). This reduction in effective gene flow increases with the strength of selection at the MSPs and decreases with their recombinational distance from the focal neutral site (Fig. 1A and SI Appendix, Eq. S1.1). To extrapolate from a given neutral site to the entire genome, we need to average over the possible genomic locations and selection coefficients of the MSPs. For simplicity, we assume an infinite chromosome with a linear relationship between physical and genetic map distance. Given an exponential distribution of selection coefficients, the expected effective migration rate depends on $s$, $\nu$, and $r_{\mathrm{bp}}$ exclusively through $\sigma /r_{\mathrm{bp}}$, where $\sigma =s\nu$ is the product of the mean selection coefficient times the density of MSPs (SI Appendix, part 1). Note that $\sigma /r_{\mathrm{bp}}$ has the meaning of a selection density per genetic map unit. For instance, conditioning on two MSPs on each side of an average neutral site, we findThis equation is a good approximation if $\sigma /r_{\mathrm{bp}}\lesssim 0.1$—that is, if recombination is at least 10 times stronger than selection, at which point effective gene flow is reduced by $\sim$50% (SI Appendix, Fig. S2). Eq. 1 shows that the mean effective gene flow decreases with selection density and increases with recombination rate. Adding increasing numbers of MSPs has a diminishing effect on $𝔼[m_e]$ (SI Appendix, Fig. S2A), so that Eq. 1 captures the essential pattern if $\sigma /r_{\mathrm{bp}}\lesssim 0.1$. The exclusive dependence of $𝔼[m_e]$ on selection and recombination through the compound parameter $\sigma /r_{\mathrm{bp}}$ holds for any number of MSPs and so applies to the genome-wide average of $m_e$ (SI Appendix, Eq. S1.8). Our results imply that doubling the number of MSPs has the same effect on average effective gene flow as doubling the mean selection coefficient. We therefore anticipate that, in practice, $s$ and $\nu$ can be inferred only jointly as $\sigma$ from population-genomic data in our framework.

To facilitate parameter inference from population-genomic data, we phrase our theory in terms of the expected coalescence time of two lineages, one from each population. The expectation of this coalescence time under neutrality, $T_B$, depends on the baseline migration rate $m$ (SI Appendix, Table S2). We incorporate the effect of selection by substituting the effective migration rate for $m$. Averaging over all possible numbers and genomic locations of MSPs, we obtain $𝔼[T_B]$ and can predict the between-population diversity ${\pi }_B$ as $2u𝔼[T_B]$, where $u$ is the mutation rate per base pair and generation.

To better reflect real genomes, we now assume a finite genome size and define $r_f=0.5$ as the recombination rate that corresponds to free recombination, such that MSPs located more than $k_f=1/(2r_{\mathrm{bp}})$ base pairs from a neutral site are unlinked. We start by assuming that $\nu$ is so small that, at most, a single, nearest-neighboring MSP is linked to any focal neutral site. In the simplest case of the (MS)M scenario with $m_0=m$ (SI Appendix, Fig. S1C), the expected pairwise between-population coalescence time is approximatelywhere $D$ and $F$ depend on $m$, $\tau$, and $\nu$ (Materials and Methods). The first two terms in Eq. 2 are the expectation without selection (21) (SI Appendix, Table S2). The third and fourth terms reflect the increase in coalescence time if the MSP is linked and unlinked to the neutral site, respectively. Importantly, the term accounting for a linked MSP shows that $\sigma /r_{\mathrm{bp}}$ strongly determines $𝔼[T_B]$, although $s$, $\nu$, and $r_{\mathrm{bp}}$ also enter Eq. 2 independently. Indeed, given $r_{\mathrm{bp}}$ and in the parameter range where Eq. 2 is a good approximation (i.e., for $\nu \ll r_{\mathrm{bp}}/s,m,\tau$), the effect of selection on $𝔼[T_B]$ is entirely captured by $\sigma$ (SI Appendix, Fig. S4). For details and other demographic scenarios, see SI Appendix, part 1. In this simplified model, the effect of all other MSPs is absorbed by $m$ as a genome-wide reduction in gene flow that is independent of recombination.

In practice, we want to explicitly account for all MSPs possibly present in the genome, as well as for the average physical chromosome length. Finding $𝔼[T_B]$ in this more realistic setting amounts to averaging over all possible numbers and genomic locations of the MSPs. We wrote a C++ program to do this integration numerically (SI Appendix, part 1). The result agrees well with individual-based forward simulations (Fig. 1B and SI Appendix, Fig. S5). As with a single MSP, if $r_{\mathrm{bp}}$ is given, $𝔼[T_B]$ depends on $s$ and $\nu$ effectively only through the selection density $\sigma$ (Fig. 1C). In fact, returning to the idealizing assumption of a global linear relationship between physical and genetic map distance (i.e., $r_f\to \mathrm{\infty }$), we show that this dependence holds exactly (SI Appendix, Eq. S1.68). This finding corroborates $\sigma$ as a key parameter and natural metric to quantify genome-wide divergent selection in the face of gene flow.

We developed an inference procedure based on our theory and applied it to two datasets from the predominantly outcrossing yellow monkeyflower (Mimulus guttatus), an important model system for speciation and local adaptation (22). For both datasets, we fit the model with multiple MSPs to the empirical relationship between recombination rate ($r_{\mathrm{bp}}$, estimated from a linkage map) and putatively neutral between-population diversity (${\pi }_B$, estimated from fourfold degenerate coding sites), after correcting the latter for genomic correlates and divergence to the outgroup Mimulus dentilobus (SI Appendix, part 3). Our procedure computes the sum of squared deviations (SSD) across genomic windows between these observed values of ${\pi }_B$ and those predicted by our model, given the estimate of $r_{\mathrm{bp}}$ for each window and a set of parameter values. Minimizing the SSD over a large grid of parameter values, we obtained point estimates for the selection density ($\sigma$), baseline migration rate ($m$), and duration of the MS phase ($\tau$). We estimated 95% nonparametric confidence intervals (CIs) for the parameters by doing a block-bootstrap over genomic windows (SI Appendix, part 3). For both datasets, we explored two alternative demographic scenarios, but focus here on the one that provided more plausible parameter estimates and tighter 95% CIs. Unless otherwise stated, we only report results obtained with genomic windows of 500 kb because results for windows of 100 and 1,000 kb were very similar (SI Appendix, part 4).

The genomes of incompletely isolated species and locally adapted populations have long been thought of as mosaics of regions with high and low divergence (31, 32). This pattern is due in part to variation in effective gene flow along the genome, created by an interaction of divergent selection and recombination (33, 34). The recent explosion of genome-wide DNA sequencing data allows us to directly observe this mosaic and has spurred theoretical and empirical studies aiming to better understand the mechanisms underlying local adaptation and speciation (e.g., refs. 35–38). However, an explicit, model-based framework linking observed genome-wide patterns of divergence with the underlying mechanism has hitherto been missing.

Here, we developed such a framework by merging the concept of effective migration rate with coalescence theory. We showed that a genome-wide negative correlation of between-population diversity with recombination rate (14, 17) can be described by the compound parameter “selection density,” such that very different genomic mosaic patterns are compatible with the same aggregate effect of divergent selection and gene flow: A large number of weak genetic barriers to gene flow (MSPs) are equivalent to a much smaller number of strong barriers. Our approach is most sensitive to polygenic selection and therefore complements existing genome scans for empirical outliers of population divergence (39–42), which tend to identify only strong barriers to gene flow. It also provides a better null model for such genome scans, because outliers could be judged against the appropriate background level of divergence.

Our approach is inspired by earlier work exploiting the genome-wide relationship between recombination rate and genetic diversity within a population for quantitative inference about genetic hitchhiking (43, 44) and BGS (45, 46). In fact, we have used an established model of BGS (47, 48) to filter out any confounding effect of BGS and gene density in a first step, before fitting our model of divergent selection against gene flow to the relationship between recombination rate and diversity between populations.

We have assumed that MSPs occur at a constant rate $\nu$ along the genome. This assumption could be relaxed by making $\nu$ depend on the functional annotation of genomes (e.g., exon coordinates), which might allow $\nu$ and $s$ to be estimated separately (49). We explored this solution heuristically for the serpentine dataset by setting $\nu$ proportional to the local density of exonic sites. Point estimates of $\sigma$ were on the same order as before (${10}^{-10}$ to ${10}^{-8.25}$). However, the 95% CIs became much wider, and the variance explained ($R^2$) dropped from $\sim$15 to 5%. This reduced goodness of fit might suggest that selection against gene flow is not acting exclusively on coding (exonic) variation.

Our model currently does not account for the clustering of locally adaptive mutations arising in tight linkage to previously established MSPs and the synergistic sheltering effect among MSPs that protects them from being swamped by gene flow (20, 50). If accounted for, this clustering would lead to an even more pronounced uptick of between-population diversity in regions of low recombination. Therefore, one might be able to use deviations from our current model in regions of low recombination as a way of detecting the presence of clustering in empirical data. At the very least, our parameter estimates would indicate whether and in what genomic regions one should expect clustering of MSPs to have evolved.

An inherent limitation of our approach is that enough time must have passed for between-population divergence to accumulate. Otherwise, there is no power to detect variation in divergence among genomic regions. This limitation constrains the temporal resolution of our method, in particular if the duration of the MS phase is short, or if strong reproductive isolation evolved so quickly that gene flow was completely and rapidly reduced across the entire genome. Another potential limitation is a relatively low resolution to infer the duration of the MS phase. A genome-wide negative correlation of recombination rate with between-population diversity will persist for a long time, even after gene flow has come to a complete halt, because subsequent neutral divergence will just add uniformly to the existing pattern. Our inference approach should therefore still provide good estimates of the strength of selection and gene flow, even after speciation has completed, as long as these estimates are interpreted as averages over the inferred time $\tau$. In this sense, our approach is likely robust to the specifics of the most recent demographic history of the populations or species of interest. To better resolve the timing of events, we suggest using the additional information contained in the entire distribution of pairwise coalescence times (and pairwise sequence differences), rather than relying on their mean, as we currently do.

The opposing roles of gene flow and selection in speciation and local adaptation have a long and contentious history in evolutionary biology and population-genetics theory (3, 51). We anticipate that the type of genome-wide quantitative inference developed here, applied to the growing amount of whole-genome polymorphism and recombination data, will help to resolve how gene flow is constraining divergent selection.

In Eq. 2, $D=\mathrm{Ei}[(1-g_f)m\tau ]-\mathrm{Ei}[(1-g_{\circ })m\tau ]$ and $F=\mathrm{Ei}[-g_{\circ }m\tau ]-\mathrm{Ei}[-g_{f}m\tau ]+\mathrm{Ei}[-\nu k_f]-\mathrm{Ei}[-\nu k_{\circ }]$, where $\mathrm{Ei}[z]=-{\int }_{-z}^{\mathrm{\infty }}{e}^{-t}/t\mathrm{d}t$ is the exponential integral. Here, $g_f={[1+s/r_f]}^{-1}$ and $g_{\circ }={[1+s/(k_{\circ }{r}_{\mathrm{bp}})]}^{-1}$ are the contributions to the reduction in gene flow if the MSP is unlinked ($k_1=r_f/r_{\mathrm{bp}}$) or fully linked ($k_1=k_{\circ }$, $0<k_{\circ }\lesssim 1/[r_{\mathrm{bp}}\tau ]$); $k_{\circ }$ is a small positive lower limit for the physical distance to the MSP. For details of our model, theory, and simulations, see SI Appendix, part 1. Statistical data analyses, bias corrections, and the inference procedure are described in SI Appendix, part 3. The Mimulus datasets (sampling design, DNA sequencing, and quality filtering) and the linkage map are discussed in SI Appendix, part 2. For complementary results, including tests of partial correlation between diversity and recombination rate, as well as the inference of BGS, see SI Appendix, part 4.

We thank Ben Blackman for help with collections; and Yaniv Brandvain, Lex Flagel, Amanda Kenney, Andrea Sweigart, Kevin Wright, Chenling Xu, and members of the G.C., Ross-Ibarra, and Schmitt laboratories at University of California, Davis for discussions and comments that improved our study. This work was supported by National Institute of General Medical Sciences of the National Institutes of Health Grants R01-GM083098 and R01-GM108779 (to G.C.); National Science Foundation (NSF) Grant 1354688 (to J.H.W.); NSF Grant 1353380 (to G.C.); NSF DDIG Grant 1110753 (to J.P.S.); and Swiss National Science Foundation Advanced Postdoc Mobility Fellowship P300P3_154613 (to S.A.). The computational results presented have been achieved in part by using the Vienna Scientific Cluster.