Reproductive isolation via polygenic local adaptation in sub-divided populations: Effect of linkage disequilibria and drift

Abstract

This paper considers how polygenic local adaptation and reproductive isolation between hybridizing populations is influenced by linkage disequilibria (LD) between loci, in scenarios where both gene flow and genetic drift counteract selection. It shows that the combined effects of multi-locus LD and genetic drift on allele frequencies at selected loci and on heterozygosity at neutral loci are predicted accurately by incorporating (deterministic) effective migration rates into the diffusion approximation (for selected loci) and into the structured coalescent (for neutral loci). Theoretical approximations are tested against individual-based simulations and used to investigate conditions for the maintenance of local adaptation on an island subject to one-way migration from a differently adapted mainland, and in an infinite-island population with two habitats under divergent selection. The analysis clarifies the conditions under which LD between sets of locally deleterious alleles allows these to be collectively eliminated despite drift, causing sharper and (under certain conditions) shifted migration thresholds for loss of adaptation. Local adaptation also has counter-intuitive effects on neutral (relative) divergence: FST is highest for a pair of subpopulations belonging to the same (rare) habitat, despite the lack of reproductive isolation between them.

Fig 1. Local adaptation under mainland-island migration.

A. Expected frequency $1-\mathbb{E}\left[p\right]$ of the locally favoured allele on the island vs. m/s, the migration rate relative to selection per locus, for various L (different colors) for s=0.02, Ns=2, and μ/s=0.005. Inset: Expected load vs. Ls (which is varied by changing L) for various values of m/s. The maximum possible load Ls (dashed line) is also shown for reference. B. Expected frequency of the locally favoured allele vs. m/s for various Ns (different colors) for s=0.02, L=40 and μ/s=0.005. Symbols depict results of individual-based simulations in both Fig 1A and 1B (obtained by averaging over 100–200 simulation replicates for each point). Colored solid lines show theoretical predictions that account for both LD and drift (obtained from Eq (7) together with Eq (4)); colored dashed lines in 1B show LE/single-locus predictions that only account for drift (and are obtained from Eq (1)). Fig 1B also shows deterministic predictions that account for LD (solid black line) as well as the LE/single-locus deterministic prediction (dashed black line). See main text for how these are calculated. Note that there are no simulation results for the deterministic case (as individual-based simulations are always affected by drift). C. Distribution of allele frequencies shown by plotting the fraction of loci with frequency of locally deleterious allele between p and p+Δp, vs. p (for Δp=0.05). The different colors show distributions for N=100, 200, 400, 800 (which correspond to Ns=2, 4, 8, 16, for s=0.02), with m/s chosen in each case such that the expected frequency of the locally deleterious allele is 0.3. Theoretical allele frequency distributions (lines) match well with those from individual-based simulations (symbols), with some (moderate) deviation in larger populations. Theoretical predictions are obtained using Eq (1) with m replaced by m_e, which depends on $\mathbb{E}\left[p\right]$, which is determined numerically, as above. D. Genotype frequencies P_y, which represent the probability that a randomly chosen genotype in the population carries y deleterious alleles, vs. y, for two different values of m/s, for Ns=2 and s=0.02. Symbols depict results of individual-based simulations; dashed lines show deterministic predictions and solid lines predictions under LE (see text for more details about the two kinds of predictions). Other parameters for C. and D. are: L=40 and μ/s=0.005.

Fig 2. Local adaptation in the infinite-island model with two habitats.

A–B. Expected equilibrium frequency of the locally adaptive allele in the rare habitat vs. m/s for ρ=0.1 (Fig 2A) and ρ=0.3 (Fig 2B), for 2 different values of L (10 and 40; squares vs. circles), and two different population sizes (corresponding to Ns=2 and Ns=4; blue vs. black). Symbols depict results of individual-based simulations; solid lines depict theoretical predictions that account for both LD and drift (obtained using Eq (8) together with Eq (4)); dashed lines depict LE (i.e., single-locus) predictions that only account for drift (obtained from Eq (2)). Selective effect per deleterious allele is s=0.02 in both plots. The number of simulated demes is D=500 in all individual-based simulations; the average allele frequency is obtained by averaging over all L loci and all ρD islands in the rare habitat, across 5 simulation replicates. C–D. Theoretical predictions for m_c/s, the critical migration threshold scaled by the per-locus selection coefficient, vs. Ls for ρ=0.1 (Fig 2C.) and ρ=0.3 (Fig 2D.), for s=0.02 and N=50, 100, 200, 400 (corresponding to Ns=1, 2, 4, 8 respectively). Here, m_c is the critical migration threshold above which local adaptation cannot be maintained in the rare habitat. Theoretical predictions are obtained by solving for the polymorphic equilibrium of Eq (8) (using Eq (4)) and determining the value of m above which no such equilibrium exists. The short horizontal colored lines along the vertical axis represent the approximate LE (single-locus) prediction $\frac{m_c}{s}\approx \frac{1}{1-2\rho }(1-\frac{1}{2Ns}\text{log}[\frac{1-\rho }{\rho }\left]\right)$ (see [38]). The exact deterministic predictions for m_c/s (obtained by solving coupled deterministic equations for p_r and p_c; see eq. 12 in Section 5 in S1 Text) are shown using solid black lines. The critical migration rate m_c/s is constant for small Ls, but then starts increasing with Ls beyond a threshold (Ls)_*. The deterministic prediction for (Ls)_* is depicted by vertical dotted lines and depends only on the habitat fraction ρ (see text). In addition, we also show approximate deterministic predictions for m_c/s (triangles and circles)—for Ls>(Ls)_* and in the highly polygenic limit s → 0, L → ∞ with Ls constant, the deterministic m_c/s is given by Eq (10a) and is shown using triangles. Predictions that are more accurate at somewhat larger s are obtained in Section 5 in S1 Text (see eq. 17B in S1 Text). These are shown using circles and agree well with the numerically obtained deterministic m_c/s (solid black lines).

Fig 3. Neutral divergence in the infinite-island model.

A. Average F_ST for a single deme in the rare $\left(F_{ST}^{\left(r\right)}\right)$ and common $\left(F_{ST}^{\left(c\right)}\right)$ habitats vs. Nm, for s=0.02, L=40, Ns=4, ρ=0.1. Here, F_ST is measured relative to the whole metapopulation. Symbols show results of individual-based simulations; dashed lines represent theoretical predictions (obtained from Eq (9a), (9b) and (9f) and using Eq (5)); the solid black line represents F_ST=1/(1+2Nm)– the prediction in the absence of local adaptation. B. Average F_ST for a pair of demes vs. Nm, for the same parameters as in Fig 3A. Here, both demes within the pair may belong to the rare habitat ($F_{ST}^{\left(rr\right)}$), or both to the common habitat ($F_{ST}^{\left(cc\right)}$), or one to the rare and the other to the common habitat ($F_{ST}^{\left(rc\right)}$). The plots show simulation results (symbols) as well as theoretical predictions (lines; obtained from Eq (9c)–(9f) and using Eq (5)). Inset (Fig 3B): The neutral diversity within demes π_W vs. neutral divergence between demes ${\pi }_B$ for rare/rare, common/common and rare/common pairs of demes (shown using red, blue and black points respectively), for Nm=1.0, as measured at a single timepoint in an individual-based simulation. Each point represents a pair of demes (i, j); π_W is computed as $1/40{\sum }_{k=1}^{40}(2p_{i,k}{q}_{i,k}+2p_{j,k}{q}_{j,k})/2$ and π_B as $1/40{\sum }_{k=1}^{40}(2p_{i,k}{q}_{j,k}+2p_{j,k}{q}_{i,k})/2$, where p_i,k represents the allele frequency at the k_th neutral locus in deme i and q_i,k=1−p_i,k. The solid lines represent (π_W, π_B) combinations that would correspond to F_ST values of 0 (orange), 0.213 (blue), 0.320 (black), 0.448 (red): the last three are the predicted $F_{ST}^{\left(cc\right)}$, $F_{ST}^{\left(cr\right)}$ and $F_{ST}^{\left(rr\right)}$ at Nm=1 respectively. C– D. $\beta =\frac{(1/F_{ST}^{\left(r\right)})-1}{(1/F_{ST}^{\left(c\right)})-1}$ vs. Nm for ρ=0.1 (Fig 3C) and ρ=0.3 (Fig 3D), for different values of L and Ns (depicted by the different colors) and s=0.02. The quantity β gives the ratio of the effective number of immigrants per unit time into an island in the rare habitat to the corresponding number for an island in the common habitat. Theoretical predictions (solid lines) match simulation results (symbols) across all parameter combinations. The black dashed line in each plot represents the threshold β_min=ρ/(1−ρ), which is the expected β under complete RI (wherein immigrants from the dissimilar habitat have zero RV). The short horizontal arrows along the vertical axis represent the threshold β=(β_min+ e^−2Ls)/(1+β_min e^−2Ls) for L=10 (upper arrow) and L=40 (lower arrow). This is the expected value of β when allele frequency divergence between habitats is maximum (see text). F_ST values in simulations are computed from 40 unlinked neutral loci.