Robust Estimation of Recent Effective Population Size from Number of Independent Origins in Soft Sweeps

Abstract

Estimating recent effective population size is of great importance in characterizing and predicting the evolution of natural populations. Methods based on nucleotide diversity may underestimate current day effective population sizes due to historical bottlenecks, whereas methods that reconstruct demographic history typically only detect long-term variations. However, soft selective sweeps, which leave a fingerprint of mutational history by recurrent mutations on independent haplotype backgrounds, holds promise of an estimate more representative of recent population history. Here, we present a simple and robust method of estimation based only on knowledge of the number of independent recurrent origins and the current frequency of the beneficial allele in a population sample, independent of the strength of selection and age of the mutation. Using a forward-time theoretical framework, we show the mean number of origins is a function of θ=2Nμ and current allele frequency, through a simple equation, and the distribution is approximately Poisson. This estimate is robust to whether mutants preexisted before selection arose and is equally accurate for diploid populations with incomplete dominance. For fast (e.g., seasonal) demographic changes compared with time scale for fixation of the mutant allele, and for moderate peak-to-trough ratios, we show our constant population size estimate can be used to bound the maximum and minimum population size. Applied to the Vgsc gene of Anopheles gambiae, we estimate an effective population size of roughly 6×107, and including seasonal demographic oscillations, a minimum effective population size >3×107, and a maximum <6×109, suggesting a mean ∼109.

Keywords: Anopheles; demographic oscillations; effective population size; recurrent mutation; soft sweeps.

Fig. 1.

Time series of the frequency of each independent origin of the same recurrent mutant (range of different colors). (A) $N={10}^6,\hspace{0.17em}2N\mu =1$, and s = 0.05, (B) same as (A), but with $2N\mu =10$. Solid black line is the sum of all mutant frequencies ($x(t)=\sum _k{x}_k(t)$), dashed black line the frequency of the wild type ($1-x(t)$), and the solid red line is the deterministic time-course given by equation (3).

Fig. 2.

Average number of origins for population sizes of $N={10}^6,\hspace{0.17em}N={10}^7$, and $N={10}^8$. The filled symbols show the simulation results and standard error bars for the parameter combinations shown in the legend; for $N={10}^6$ and $N={10}^7$, the simulations used multinomial sampling of the Wright–Fisher drift process with 50 and 10 replicates, respectively, for each parameter combination, whereas for $N={10}^8$, the multinomial sampling is replaced by the multivariate Gaussian distribution approximation of the drift process (see the Methods section above), where 100 replicates are used in this plot. The solid thick lines are the predictions for the same parameter combination of the semideterministic theory described in this article (Methods), whereas the thin lines represent the prediction of Pennings and Hermisson (2006a), based on Ewens’ sampling theory (Ewens 2010).

Fig. 3.

Average number of origins for population size of $N={10}^8$ on linear-log scale, for $2N\mu =\{1,10,100\}$ and $s=\{0.05,0.005\}$ showing that plateau number of origins is independent of s. The filled symbols show the simulation results and standard error bars for the parameter combinations shown in the legend. The solid thick lines are the predictions for the same parameter combination of the semideterministic theory described in this article (Methods).

Fig. 4.

Distribution of the number of origins for simulations with various mutation rates for $N={10}^8$ and s = 0.05 (open circles) compared with theory in this article equation (12) and (7) (solid lines) and Ewens’ sampling formula (dotted lines), both with n_s = 1,000. For the mutation rates $2N\mu =\{0.1,1,10\}$, the corresponding typical fixation time (eq. 5) is $t^{*}\approx \{370,320,280\}$ generations.

Fig. 5.

${\hspace{0.17em}\text{log}\hspace{0.17em}}_{10}$-error in estimating the true effective population size, for (A) haploid populations with $N={10}^8$, (B) diploid populations with $N=5\times {10}^7$, for various selection coefficients, mutation rates, and dominance coefficients (diploid only) from Wright–Fisher simulations (100 replicates for each parameter combination). (A) We use equations (12) and (7) to determine the maximum likelihood estimate. (B) For the diploid population, we use the same Poisson likelihood function, but with mean given by equations (13) and (14) in the Supplementary Information, where we assume perfect knowledge of T (squares) and also compare to the case where we have a systematic error in our knowledge of T, where the true time is $T/2$ instead T (circles), and we see the estimates are unchanged. In addition, for the diploid population we use the haploid likelihood function (eqs. 13 and 7) with $\theta =4N\mu$ to estimate N (plus signs) and find again excellent agreement.

Fig. 6.

Mean number of origins for haploid simulations with preexisting mutations (A), where the black hexagram symbols represent simulations without preexisting simulations, and (B) ${\hspace{0.17em}\text{log}\hspace{0.17em}}_{10}$-error in maximum likelihood estimate of the true effective population size $N={10}^8$ from Wright–Fisher simulations with various values of the deleterious selection coefficient s_d (100 replicates for each parameter combination).

Fig. 7.

Likelihood (normalized) of the number of origins as function of effective population size given an observed number η  =  10 and samples size n_s = 1,530 chromosomes, corresponding to that found for the Ag1000 project (Anopheles gambiae 1000 Genomes Consortium 2017) for the Vgsc resistance locus. As shown in the legend, the semideterministic theory in this article, assuming a current day frequency of x = 0.78 (as observed) is compared with assuming x = 1 and the Ewens’ sampling theory equation (14), which only has applicability for x = 1. The 95% confidence intervals (gray dotted lines) and maximum likelihood effective population size (red dotted line) are shown for the semideterministic likelihood function with x = 0.78.

Fig. 8.

The mean number of origins from Wright–Fisher simulations (1,000 replicates) for oscillating population size with period $\Delta T=12$ generations, selection coefficient s = 0.05, $2N\mu =1$, and with the geometric mean (green) and harmonic mean (purple) of N_max and N_min constrained to $\sqrt{N_{\text{max}}{N}_{\text{min}}}=2{(1/N_{\text{max}}+1/N_{\text{min}})}^{-1}=N={10}^8$, for different peak-to-trough ratios. Black squares represent constant population size simulations.