Estimating the Effective Population Size from Temporal Allele Frequency Changes in Experimental Evolution Free

Abstract

DURING experimental evolution studies, populations are maintained under specific laboratory conditions (Kawecki et al. 2012; Long et al. 2015; Schlötterer et al. 2015). In sexually reproducing organisms, the census population size is typically kept fixed at fairly low numbers, rarely exceeding 2000 individuals. With such small population sizes, genetic drift causes stochastic fluctuations in allele frequencies. Under neutrality, the level of random frequency changes is determined by the effective population size (⁠$Ne$⁠) (Wright 1931). Furthermore, the efficacy of selection is influenced by $Ne.$ For weakly selected alleles, the probability of fixation is directly proportional to the product of $Ne$ and the intensity of selection (Fisher 1930; Kimura 1964). As changes in allele frequency are greatly affected by the population size, it is fundamental to estimate $Ne$ accurately to understand molecular variation in experimental evolution studies.

Krimbas and Tsakas (1971) estimated $Ne$ using the standardized variance of allele frequency (F, see also Falconer and Mackay 1996) from longitudinal samples in natural populations of olive flies. As F was calculated from these samples, they accounted for the sampling variance that also contributed to the true allele frequency variance. This approach was further improved and used by several authors (Nei and Tajima 1981; Pollak 1983; Waples 1989; Jorde and Ryman 2007).

With the widespread availability of powerful computers, also maximum-likelihood-based methods became popular (Williamson and Slatkin 1999; Anderson et al. 2000; Wang 2001; Hui and Burt 2015) in addition to the moment-based approaches discussed above. Although these methods show less bias than the moment-based approaches (Wang 2001), they are still computationally demanding, in particular for the large numbers of markers typically obtained with novel sequencing technologies (Foll et al. 2015).

Estimating $Ne$ with temporal methods requires samples collected at least at two time points. Alternative methods that use only a single time point are based on linkage disequilibrium (LD) (Hill 1981; Waples and Do 2008, 2010; Waples and England 2011), heterozygote excess (Pudovkin et al. 1996), molecular coancestry (Nomura 2008), sibship frequencies (Wang 2009, 2013), or combinations of summary statistics using approximate Bayesian computation (Tallmon et al. 2008). LD-based methods are widespread but require haplotype or unphased diploid genotype information, which limits their applicability.

Although the cost for sequencing has dropped considerably, the separate sequencing of thousands of individuals in replicate populations in experimental evolution studies is still out of reach. Sequencing samples in pools (Pool-seq) can provide a cost-effective alternative (Schlötterer et al. 2014). Pool-seq has also been shown to outperform individual-based sequencing in estimating allele frequencies and inferring population genetic parameters under several conditions (Futschik and Schlötterer 2010; Zhu et al. 2012; Gautier et al. 2013). For these reasons, Pool-seq has become the basis of many experimental evolution “evolve and resequence” (E&R) studies (Turner et al. 2011; Schlötterer et al. 2015). Following the emergence of E&R, many population genetic estimators have been adjusted to handle the properties of Pool-seq data (Futschik and Schlötterer 2010; Kofler et al. 2011a,b; Kolaczkowski et al. 2011; Boitard et al. 2013; Ferretti et al. 2013). To the best of our knowledge, no $Ne$ estimators have been developed so far that properly deal with Pool-seq data.

In this article, we present a novel temporal method to estimate $Ne$ from pooled samples. We show that previously proposed estimators can lead to substantial bias, as they neglect the variance component due to pooled sequencing. We introduce a new model accounting for the two-step sampling process associated with Pool-seq data. In the first sampling step individuals are drawn from the population to create pooled DNA samples. In the second step, pooled sequencing is modeled as binomial sampling of reads out of the DNA pool. We show on simulated data that our method outperforms standard temporal $Ne$ estimators. For real data, we also suggest to use a segmentation algorithm, to partition the genome-wide sequence data into stretches of DNA with significantly different $Ne$ estimates. Finally, we present an application to a genome-wide experimental evolution data set from Drosophila melanogaster (Franssen et al. 2015).

Materials and Methods

Sampling schemes

Nei and Tajima (1981) investigated the sampling properties of temporal $Ne$ estimators and proposed two different sampling schemes. Under the first scheme (plan I), individuals are either sampled after reproduction or returned to the population after their genotypes have been examined. In contrast, under the second scheme (plan II) sampling takes place before reproduction and the sampled individuals are permanently removed from the population and their genotypes do not contribute to the next generation. By assuming different sampling distributions, they derived separate $Ne$ estimators under sampling plans I and II.

Waples (1989) unified the calculations under the two plans by assuming binomial sampling out of an infinitely large parental gamete pool for both sampling schemes. He concluded that the measure of variance under the two sampling plans differs only in a covariance term. For plan I, there is a positive correlation between allele frequencies sampled t generations apart because they are both derived from the same population at generation 0. In contrast, for plan II, the initial sample and individuals contributing to the next generation can be considered as independent binomial samples; thus sample frequencies at generations 0 and t are uncorrelated.

For a typical E&R study, outbred experimental populations are created by mixing a large number of inbred lines (e.g., Turner and Miller 2012; Huang et al. 2014; Franssen et al. 2015). The populations are then propagated under the desired experimental conditions while keeping the census size of the population controlled through time (Figure 1). However, the experimenter has no direct influence on the effective population size, which is in general lower than the census size. In E&R studies with Drosophila, the census size rarely exceeds some hundreds of individuals, and sampling usually takes place after reproduction according to plan I. For organisms maintained at larger sizes, such as yeast, the sample for genetic analysis is not returned to the population (Burke et al. 2014). Plan II applies to such cases.

Figure 1

Two-step sampling in experimental evolution with Drosophila. In E&R studies, populations are propagated at a census size N defined by the experimenter, which is in general larger than the effective population size $Ne.$ Using temporal methods, $Ne$ can be estimated from the variance in allele frequency between samples taken t generations apart. To get an accurate representation of allele frequencies in population genetic studies, a large number of individuals $Sj$ (⁠$j∈{0,t}$⁠) are sampled and pooled. Sampling can take place according to sampling plan I or II based on the mode of reproduction. Pooled samples are then subjected to high-throughput sequencing. Sequenced reads are subsequently aligned to the reference genome (shown at the bottom). We represent pool sequencing by an additional sampling step (called sampling step 2). We correct for both sampling steps when estimating $Ne$ in pooled samples. Additionally, we take into account variable coverage levels across the genome (coverage $Rij$ for site i at $T=j,$  $j∈{0,t}$⁠) when correcting for the variance coming from sequencing.

Open in new tabDownload slide

In E&R studies, sampled individuals are often pooled together for DNA extraction (Schlötterer et al. 2014). The size of the pool can be as large as the whole population. Depending on the experimental design, it is also possible that only a fraction of the population is sequenced, for instance, only females (Tobler et al. 2014; Franssen et al. 2015). Pooled individuals are used to create DNA libraries, which are, in turn, subjected to high-throughput sequencing.

We consider two separate sampling steps when estimating $Ne$ from Pool-seq samples (Figure 1). In the first step, we model the sampling of individuals out of the population. This can take place according to either plan I or plan II. In the second step, we model the sequencing of a DNA pool by drawing reads at random with replacement from the first-step sample. The allele frequency variance inferred from the sample is corrected for the additional variance coming from the two-step sampling and used for estimating $Ne.$

Notation

We assume that the experimental population is propagated at a constant census size N and that $N≥Ne.$ We use genome-wide single-nucleotide polymorphism (SNP) data sampled t generations apart to estimate $Ne$ (Figure 1) and denote the estimated effective population size by $N^e.$ Multiallelic sites in populations with low mutation rates, such as Drosophila, exist but are rare and likely to be sequencing errors (Burke et al. 2010). Therefore we consider only biallelic SNPs at n polymorphic sites. At each site i (⁠$i=1…n$⁠) the true population allele frequency is denoted by $pij$ at time $T=j,$ where $j∈{0,t}.$ To obtain allele frequency estimates for an unknown $pij,$ the population is subjected to sampling. We consider two sampling steps (Figure 1). At $T=j,$ we first sample $Sj$ individuals out of the population to create a pooled DNA library for sequencing. Note that the number of sampled individuals is constant over the n sites, and therefore the index i is omitted here. Sampling individuals can take place according to either plan I or plan II, as described above (also shown in Supplemental Material, Figure S1). As the second sampling step, we model Pool-seq by drawing $Rij$ reads out of the pooled DNA sample at each site i (⁠$i=1…n$⁠). This allows for variation in sequence coverage. Below we derive the variance in allele frequency for a given site. To keep notation simple, we omit again the index i and denote the unknown sample allele frequency among the $S0$ individuals at the first sampling time point (⁠$T=0$⁠) by x and the subsequent allele frequency estimate obtained via pool sequencing from $R0$ reads by $x^.$ Similarly, at some $T=t,$ the respective frequencies are denoted by y and $y^.$ Note that under pool sequencing only $x^$ and $y^$ are observed.

Estimating $Ne$ from temporal allele frequency changes

Possible sources of bias in $Ne$ estimators were later investigated by Jorde and Ryman (2007). The authors pointed out that the expectation over F is typically approximated by taking the expected values separately for the numerator and the denominator (Turner et al. 2001). They suggested a different weighting scheme of alleles leading to an alternative less-biased estimator to measure temporal frequency change.

Correction for two-step sampling

We consider a random-mating population with discrete generations. Neutral evolution is assumed with no selection, migration, and mutation. Samples are drawn from the population at generations $T=0$ and t. Throughout the derivation we consider diploid populations, and therefore a sample of $Sj$ individuals leads to $2Sj$ sequences at times $T=j∈{0,t}.$ Sampling is assumed to be binomial with parameters $2Sj$ and $pj$ (Waples 1989). In the second sampling step at time $T=j,$ sequencing a random pool $Rj$ of reads is also modeled as binomial sampling.

Long time series have recently become available for some E&R experiments (Barrick et al. 2009; Burke et al. 2010, 2014). Standard $Ne$ estimators (Krimbas and Tsakas 1971; Nei and Tajima 1981; Waples 1989) assume a small number of generations $(t)$ and approximate $Ne$ using $2Ne≈t/F.$ If, however, $t/Ne$ is larger, using this approximation can lead to severe bias (Figure S2). To avoid such a bias, we use Equation 2 to estimate $Ne.$

Simulations

We evaluate the performance of our estimator on data simulated under the neutral Wright–Fisher model. With a given population size of $2Ne,$ we simulate the frequency trajectory of n independent SNPs. As we focus on biallelic SNPs, we assume two possible nucleotides (alleles) to be present in the population with given frequencies at the start. To create a new generation, nucleotides are drawn independently at random with a probability given by their respective allele frequencies in the old generation. The population is propagated at a constant size of $2Ne$ for t nonoverlapping generations. The effective population size is then estimated from allele frequencies inferred from Pool-seq samples taken from the population at the start and after t generations. The sampling of individuals to create the pooled DNA library is simulated by using sampling without replacement. To model the uneven coverage of genome-wide sequence data, we simulate a random coverage for each site, using a Poisson distribution with parameter equal to the given mean coverage. For every position, reads are then generated by binomial sampling from the library with sample size equal to the local coverage.

We assess the performance of our estimator for various combinations of $Ne,$ pool size, coverage, number of SNPs, and distribution of starting allele frequencies. Additional to these parameters, we also test how the ratio between census and effective population size (⁠$r=N/Ne$⁠) affects the accuracy of the proposed estimator. For this purpose, we increment the population size to a desired level of N in each generation while keeping the allele frequencies unchanged to avoid introducing additional sampling variance. We simulated each scenario 100 times.

Linkage disequilibrium between loci can reduce the number of independent SNPs, thereby increasing the variance of the estimate. The impact of dependence between SNPs is investigated based on 10 replicated whole-genome forward simulations with recombination, using the software tool MimicrEE (Kofler and Schlötterer 2014). As a starting population for the forward simulations, we sampled 2000 haploid genomes out of 8000 genomes simulated with fastsimcoal v.1.1.2 (Excoffier and Foll 2011; Bastide et al. 2013). The parameters used to generate the genomes mimic a wild population of D. melanogaster from Vienna (Fiston-Lavier et al. 2010; Bastide et al. 2013; Kofler and Schlötterer 2014). Allele counts are subjected to binomial sampling to mimic Pool-seq with a given sequence coverage. $Ne$ is estimated in nonoverlapping windows, each containing a fixed number of SNPs.

Estimating $Ne$ on simulated data

We denote our estimator corrected for the additional sampling step, i.e., pooling, by $Ne(P).$ We compare $Ne(P)$ to the standard estimators $Ne(W)$ and $Ne(JR)$ proposed by Waples (1989) and Jorde and Ryman (2007) that correct only for a single sampling step.

We illustrate experimental sampling procedures by considering two major scenarios: (i) The full population is sequenced as one large pool and (ii) only a subset of the population is used to create pooled samples. Under scenario (i) we simulate only a single binomial sampling step to represent sampling reads out of the DNA pool. The pool size is set to be equal to the census size of the population (⁠$Sj=N$⁠), and the number of sampled reads (⁠$Rij$⁠) represents the per-site coverage. For estimators that correct only for a single sampling step, we use the coverage (⁠$Rij$⁠) as the sample size. For scenario (ii), the sampled individuals (⁠$Sj$⁠) and the read number (⁠$Rij$⁠) represent the pool size and coverage for $Ne(P).$ The coverage (⁠$Rij$⁠) is taken as the sample size for the $Ne(W)$ and $Ne(JR)$ estimators, as these methods consider only one sampling step.

Change point inference for genome-wide estimates

The effect of genetic drift on the variance in allele frequency specified in Equation 1 holds only under the assumptions of Wright–Fisher neutral evolution. Deviations from the Wright–Fisher model, such as the presence of selection or demography, may cause systematically different changes in allele frequency, affecting the variance and causing locally variable patterns in genetic diversity. Furthermore, the effect of selection on one site of the genome may cause changes in the behavior of variants at nearby sites (Maynard Smith and Haigh 1974; Barton 2000; Comeron et al. 2008). As a result, the estimates of $Ne$ at different locations of the genome may deviate from the true number of breeding individuals in the population (Kimura and Crow 1963; Charlesworth 2009). For example, regions under background selection are associated with reduced $N^e$ values that extend to linked sites due to the Hill–Robertson effect (Charlesworth 1996, 2012a; Comeron et al. 2008). Similarly, selectively favorable alleles can also drag nearby neutral sites to high frequency (Maynard Smith and Haigh 1974), causing a local reduction in the estimated $Ne$ (Liu and Mittler 2008). Such an event is also known as a selective sweep (Berry et al. 1991). On the other hand, we expect the opposite pattern, i.e., a local elevation of $N^e$ for types of selection such as balancing selection that maintain variation in the genome (Baysal et al. 2007; Charlesworth 2009).

To detect such patterns in $N^e,$ we apply a segmentation algorithm to partition the genome into locally homogeneous $N^e$ stretches. We use a method related to an approach suggested by Futschik et al. (2014) for partitioning DNA sequences with respect to GC content. It is based on a statistical multiscale criterion and provides statistical error control, in the sense that the estimated number of windows will not exceed the true one except for a small error probability α to be specified by the user. With our $Ne$ estimates, we use a criterion proposed by Frick et al. (2014) for normally distributed responses. It is implemented as part of the R package stepR (Frick et al. 2014). By using simulations with selection we also illustrate that this method is able to capture the signal of locally variable $N^e$ along the chromosome.

Data availability

We estimated N_e in an E&R study with D. melaongaster, published in Orozco-terWengel et al. (2012) and Franssen et al. (2015). Pool-seq read libraries from these studies are available at the European Sequence Read Archive at http://www.ebi.ac.uk/ena/ under accession nos. ERP001290 and ERS460611–ERS460613.

Results and Discussion

Two-step correction is vital to avoid large bias in $N^e$ with Pool-seq data

Methods that do not correct for the additional sampling step caused by pooling can lead to substantial bias in $N^e$ as illustrated in Figure 2. Using simulated data, we compare our proposed estimator $Ne(P)$ to two commonly used estimators $Ne(W)$ (Waples 1989) and $Ne(JR)$ (Jorde and Ryman 2007) that provide highly accurate estimates when only a single sampling event is simulated (Figure S3). Figure 2 shows that the additional correction substantially decreases the bias for almost all scenarios (see also Figure S4, Figure S5, and Figure S6). Under plan I, $Ne(P)$ is nearly unbiased. The plan II version of the estimator has a slight upward bias when applied on data simulated under plan I, if the samples are taken at very close time points.

Figure 2

Effective population size estimated with different methods. Sixty generations of Wright–Fisher neutral evolution with $Ne=100$ diploid individuals were simulated for n = 2000 unlinked loci (SNPs). Prior to sampling, the population was increased to a census size of $N=500$ individuals at each generation. At the starting population and at each indicated time point a sample was taken to create a pool of $S=100$ individuals. The pool was sequenced to an average coverage of $R=50$ and $Ne$ was estimated on the resulting data set by separately contrasting allele frequencies at generation 0 to each of the evolved generations denoted on the x-axis, using $Ne(P),$  $Ne(W)$ (Waples 1989), and $Ne(JR)$ (Jorde and Ryman 2007). Each box represents results from 100 simulations with identical parameters. The dashed gray line shows the true value of $Ne.$ Data are simulated under plan I assumptions and the results of plan I and II estimators are shown in the left and right panels, respectively.

Open in new tabDownload slide

As an alternative approach, we also estimated $Ne$ separately for each locus, using $F′c$ in Equations 9 and 11. We then calculated the effective population size across the n loci as the harmonic mean over the single-locus $N^e$ estimates (⁠$N^e*(P)$⁠) (Peel et al. 2013). In our simulations, the harmonic mean estimator shows an accuracy similar to that of the original $N^e(P)$ (Figure S7). However, for t lying in the midrange of the simulated interval (⁠$t=$ 15–40), $N^e*(P)$ is slightly more biased under plan I.

Because of the additional sampling variance, both $Ne(W)$ and $Ne(JR)$ have a downward bias in particular for small t. Furthermore, $Ne(W)$ is upwardly biased for larger values of t, probably reflecting that alleles closer to fixation or loss are contributing less to the variance (Waples 1989). The drift variance accumulates with an increasing number of generations, while the sampling variance stays constant, making the initial bias of $Ne(JR)$ less pronounced for larger t. When samples are collected only a few generations apart, the variance of $Ne(P)$ estimators tends to be larger than that of $Ne(W)$ and $Ne(JR)$ under both plans.

Plan I and II estimators differ by a factor resulting from the covariance between the sample frequencies at generations 0 and t (Equation 10), which is inversely proportional to the census population size. Consequently, the difference between plans I and II becomes smaller for increasing N. Waples (1989) investigated how the ratio between census and effective population size (⁠$r=N/Ne$⁠) affects the accuracy of the estimators and concluded that the ratio of $r≥2$ is sufficient to reach similar estimates for both sampling schemes. We tested the performance of $Ne(P)$ on simulated data with $Ne=100$ and $N:Ne$ ratios of $r=1,2,5$ with different coverages and pool sizes (Figure S4, Figure S5, and Figure S6). When $N=Ne,$ the $Ne(P)$ plan I method achieves highly accurate estimates for all time points in contrast to the other methods (Figure S4). If, however, the $Ne(P)$ plan II estimator is applied to data simulated under plan I, we observe an upward bias for small t, which improves with an increasing number of generations. This pattern is not unexpected since the missing covariance term becomes less influential in view of the increasing drift variance after several generations. When the entire population is sequenced as a single pool (⁠$S=100$⁠), the plan II estimators of Waples (1989) and Jorde and Ryman (2007) perform similarly to the $Ne(P)$ plan I estimator because the correction for pooling in $Ne(P)$ cancels out the additional covariance term when $S=N,$ making the term used as F approximately identical to that of $Ne(JR).$ This is a general pattern irrespective of r.

For $r≥2,$  $Ne(P)$ plan I remains highly accurate (Figure S5 and Figure S6). Furthermore, when increasing the census size under a constant $Ne$ (equivalent to increasing r), the covariance between sample allele frequencies decreases, making the difference between plans I and II almost negligible (Waples 1989). The sampling variance becomes proportionally smaller compared to the drift variance with an increasing number of generations between the samples. This improves our ability to accurately estimate $Ne.$

Correcting for the additional variance inherent to Pool-seq leads to an improved performance of $Ne(P)$ compared to the standard methods for both plans. In general, with Pool-seq data the extent of the bias of the $Ne(W)$ and $Ne(JR)$ estimates depends on the ratio between N and S, smaller sample sizes (S) leading to a larger bias. As we accounted for the sequencing step with these estimators (Estimating $Ne$  on simulated data), decreasing the coverage at a given pool size does not change the bias much but rather increases the variance of the estimators.

In most of the experimental studies the investigator has control over the census population size; thus requiring the knowledge of N for $Ne(P)$ plan I does not restrict the analysis. We illustrate the performance of $Ne(P)$ plan I only when $Ne=N$ in the main text, but according to Figure S5 and Figure S6, $Ne(P)$ plan I is also highly accurate when $r≥2.$

We show the coefficient of variation (CV) of the $Ne(P)$ plan I estimator in Figure 3. The CV is defined as the ratio between the standard deviation and the mean (⁠$CV=σ^/μ^,$ where both $σ^$ and $μ^$ are estimated from the sample). It measures the relative dispersion of the distribution of the estimated values. $Ne(P)$ estimators are highly precise in nearly all cases, except when the drift variance is negligible compared to the sampling variance (Figure 3; see also Figure S9 and Figure S11 where $Ne=1000, t<30, S≤100,$ and $R=50$⁠). The bias is coming from a few outlier estimates, but the median shows more robust results (Figure S13). For plan II estimators, the behavior of the method is similar (Figure S8, Figure S10, Figure S12, and Figure S14). Note that the simulations underlying Figure S8, Figure S10, Figure S12, and Figure S14 have been done under plan I.

Figure 3

Coefficient of variation of $Ne(P)$ under plan I for various parameter values. Neutral Wright–Fisher simulations were performed with various combinations of the parameters: effective population size (⁠$Ne=100, 500, 1000$ diploid individuals), pool size (⁠$S=100,50$⁠), and coverage (⁠$R=150, 100, 50$⁠). $Ne$ was estimated with $Ne(P)$ under plan I, using $n=2000$ SNPs. $S=N$ indicates scenarios when the whole population is sequenced as a single pool. For all simulations, we assumed $N=Ne.$ Each value is calculated over 100 simulations. When the coefficient of variation exceeds one, the inset shows the actual value.

Open in new tabDownload slide

Increasing the number of SNPs reduces the variance of $Ne(P)$

We test how the number of loci (n) used to infer $Ne$ affects the accuracy and the precision of the estimates by gradually increasing the number of independent SNPs from 100 to 10,000 (Figure 4). We observe a larger variance and a slight downward bias for a small number of SNPs (100 SNPs). Both the bias and the variance become smaller with a larger number of SNPs. Some further improvement is obtained when >10,000 SNPs are used (not shown), but the benefit of additional independent SNPs levels off. We conclude that $n=2000$ SNPs usually provide sufficient precision and accuracy. However, when linkage disequilibrium is present in a genome-wide data set, the number of truly independent SNPs per window is reduced and a larger number of loci is recommended.

Figure 4

Effect of the number of SNPs used for estimating $Ne.$ The effective population size is estimated using $Ne(P)$ plan I on simulated data with $Ne=N=100.$ A total number of $S=100$ individuals are pooled and sequenced at a mean coverage of $R=50.$ Based on 100 simulation runs, $Ne$ is estimated using different numbers of SNPs at multiple time points.

Open in new tabDownload slide

A skewed starting allele frequency distribution only moderately increases the variance of $Ne(P)$

In natural populations, the neutral site frequency spectrum is skewed toward allele frequencies close to the boundaries. $Ne(P)$ uses a weighting scheme that is not very sensitive to this skew (see also Jorde and Ryman 2007). This makes it robust with respect to the shape of the starting allele frequency distribution. We illustrate this with a simulated data set having a starting allele frequency distribution that is skewed toward low- and high-frequency variants (Beta(0.2, 0.2)) as expected under neutrality. The estimates of $Ne$ from such data sets are compared to simulated data with matching parameters but uniform starting allele frequency distribution (Figure 5). We observe a very slight upward bias with neutral starting allele frequencies compared to uniform and a moderate increase in the variance given $t≥15.$ With an increasing number of generations, the difference becomes negligible.

Figure 5

Influence of the starting allele frequency distribution on $Ne(P)$ under plan I. A comparison between uniform and Beta(0.2, 0.2)-distributed (neutral) starting allele frequencies is shown. The simulation parameters match those of the genome-wide simulations in Figure 6.

Open in new tabDownload slide

The presence of linkage disequilibrium does not have a large effect on the precision of $Ne(P)$

We investigated the sensitivity of our estimator to linkage disequilibrium between loci, using genome-wide neutral simulations with recombination (Kofler and Schlötterer 2014). We simulated data with three different rates of recombination: high, normal, and no recombination. For the first case, the recombination rate is set to mimic the behavior of almost independent SNPs. In the normal recombination rate scenario, we use D. melanogaster recombination rates (Fiston-Lavier et al. 2010). The effective population size was estimated in nonoverlapping windows with a fixed number of $n=10,000$ SNPs (Figure 6). Different levels of linkage disequilibrium affect the number of independent loci per window. Nevertheless, we observe only a slight increase in the precision of the $Ne$ estimates with increasing recombination rate (Figure 6).

Figure 6

Effect of linkage disequilibrium on $N^e.$ The effect of linkage disequilibrium on our estimator was evaluated based on a whole-genome forward simulation with recombination using the software MimicrEE (Kofler and Schlötterer 2014). Three sets of simulations were performed with different rates of recombination: high, normal, and no recombination. For each parameter setup, a genome-wide simulation is replicated 10 times. The effective population size was estimated with $Ne(P)$ (plan I) in nonoverlapping windows of n = 10,000 SNPs for each replicate. The box plots show the distribution of $Ne$ estimates across replicates and windows.

Open in new tabDownload slide

Heterogeneous $N^e$ along the genome in an E&R study with D. melanogaster

We estimated $Ne$ in a recent E&R study with D. melanogaster (Orozco-terWengel et al. 2012; Franssen et al. 2015). In this experiment replicate populations of 1000 individuals were subjected to a fluctuating hot environment for 59 generations. Allele frequency estimates were obtained for founder and evolved populations, using Pool-seq. $Ne$ was estimated based on the allele frequency changes between founder and latest evolved populations in nonoverlapping windows containing 10,000 SNPs, using $Ne(P)$ under plan I. To determine the number of DNA stretches with different $N^e$ along the genome, we use a segmentation algorithm provided in the software tool by Frick et al. (2014). This method requires homogeneity of variances. Since the variance of estimates increases with $Ne,$ the estimates were log-transformed before applying the partitioning procedure. The obtained step function was back-transformed to the original scale and is shown for three biological replicates (Figure 7).

Figure 7

Genome-wide $N^e$ from an E&R study with D. melanogaster. $Ne$ is estimated based on the allele frequency changes between founder and evolved populations at generation 59 (Franssen et al. 2015). In the top panel, we show genome-wide estimates calculated with $Ne(P)$ (plan I), using $N=1000$ as census size and $S=500$ as pool size (Orozco-terWengel et al. 2012) and nonoverlapping windows of 10,000 SNPs. Chromosome-wide mean estimates across replicates are shown by the dashed lines and also calculated separately for each replicate in Table 1. DNA stretches with significantly different $N^e$ are determined using the stepR software package (Frick et al. 2014) (bottom panel). Lower and upper $1−α$ confidence bands are shown as shaded areas. α controls the error, i.e., the probability for overestimating the number of change points, and is calculated automatically as described in Frick et al. (2014). The colors indicate different biological replicates.

Open in new tabDownload slide

The mean and the median estimates for each chromosome arm as well as across the genome are stable across replicates (see Table 1 and Table S1). As experimental evolution studies often aim to find signals that are consistent across replicates, this can be an important check of the experimental setup. On the other hand, we see differences between chromosome arms. For example, the mean is clearly lower for 3R, emphasizing the added value of spatial analysis compared to genome-wide estimates.

Genome-wide mean $N^e$ from an E&R study with D. melanogaster

In D. melanogaster  $N^e$ ranges between ∼100 and 400. Around the centromere of chromosome 2, the estimated $Ne$ decreases by two-thirds in replicates 1 and 3, which is in agreement with the expectation of low diversity and, as a consequence, low $Ne$ in regions with reduced recombination (Begun and Aquadro 1992; Presgraves 2005; Haddrill et al. 2007; Campos et al. 2012). Furthermore, $N^e$ is low on the entire chromosome arm 3R and also on parts of 3L. Overall, these patterns can be attributed to strong LD, caused either by low recombination rates around the centromeres (Chan et al. 2012) or by segregating inversions (Kapun et al. 2014) in combination with selection potentially on rare variants. The reduction in $N^e$ is also well captured by the segmentation algorithm (Figure 7), which shows a similar pattern when applied on simulated data with selection (Figure S15). These results are consistent with those of Tobler et al. (2014), who observed a massive amount of outlier SNPs around the centromere of chromosome 2 and on 3R. Interestingly, certain regions of the genome show extensive differences in $N^e$ between the replicates, which might be reflecting different selection histories or differences in demography, such as replicate-specific bottlenecks.

$N^e$ may also vary as a result of differences in the modes of transmission of different components in the genome. For example, on the X chromosome, $Ne$ is equal to three-quarters of the autosomal population size (Vicoso and Charlesworth 2006, 2009). Interestingly, our estimates in the E&R experiment do not reflect this expectation of reduced effective population size. Instead, we estimate $Ne$ to be as high as $N^e$ on the autosomes. Unequal sex ratio between males and females can be a source of such a pattern (Charlesworth 2009); however, unbalanced sex ratio has not been reported in this experiment. Another possible explanation for increased $N^e$ on the X can be the presence of background selection as suggested by Charlesworth (2012b). He argues that because of the lack of recombination in male Drosophila, the effect of background selection is more effective on the autosomes than on the X chromosome. Orozco-terWengel et al. (2012) reported differences in the number of putatively selected sites between the X and autosomes. They found that candidate SNPs were underrepresented on the X. Their selection scan identifies signatures of deviation from neutral expectation, which is also reflected in the reduction in $N^e$ on the autosomes, indicating higher selection pressure.

Recommendations for genome-wide data sets

Most of the methods proposed previously are not designed for genome-wide high-density SNP data sets. However, the method of Jorde and Ryman (2007) was successfully used for genome-wide data by Foll et al. (2014). Reed et al. (2014) also used a similar approach to estimate $Ne$ for whole-genome data, using sliding windows. We estimated $Ne$ in windows with a fixed number of SNPs. Using windows of fixed lengths in base pairs would affect the variance of the estimator (Figure 4) but does not distort the mean. All these approaches, however, do not account for the ruggedness of the recombination landscape and can lead to windows with different levels of linkage disequilibrium in them. To overcome this problem it would be possible to define windows based on recombination distance. Unfortunately, the lack of haplotype information in the Pool-seq data makes it difficult to infer linkage disequilibrium. One way to infer linkage information from pooled sequence data is provided by the software LDx (Feder et al. 2012). For model organisms such as Drosophila, readily available recombination maps can also be used as a proxy (Przeworski et al. 2001; Kulathinal et al. 2008; Fiston-Lavier et al. 2010). If only a single genome-wide $Ne$ estimate is required, one can alternatively use a set of randomly distributed SNPs over the genome to obtain $N^e.$

Temporal methods make a number of assumptions, which, if violated, can lead to biased $Ne$ estimates. For example, in our simulations, we considered only effective population sizes that are constant over time. Fluctuating $Ne$ is a frequent phenomenon in natural populations and can be an important component of an experimental design. For example, in repeatedly bottlenecked populations, the smallest population size dominates $N^e$ (Luikart et al. 1999; Charlesworth 2009). But even in strictly controlled populations the experimental regime can induce changes in $Ne.$ When the population changes in size, the estimated $Ne$ is generally interpreted as the harmonic mean of the effective population sizes over the generations (Wright 1938; Nei and Tajima 1981; Waples 1989). However, if time series allele frequency data are available, such changes can be detected by estimating $Ne$ from pairwise comparisons between consecutive time points.

All evolutionary forces (selection, demography, etc.) that lead to deviations from the neutral expectation will also affect our estimate. Nevertheless, systematic forces that result in locally different values of $N^e$ can be detected with a sliding-window approach, as illustrated with simulations under selection (Figure S15). The D. melanogaster data set also illustrates this point; i.e., the hypothesized regions under selection coincide with regions of reduced $Ne$ (Orozco-terWengel et al. 2012; Tobler et al. 2014; Franssen et al. 2015). For this to be detected, however, most of the allele frequency change has to occur over the sampled time span.

In the E&R study with D. melanogaster, shown above, the criterion of nonoverlapping generations, assumed by temporal methods, is met (see Orozco-terWengel et al. 2012 for details on experimental design). However, for samples from an age-structured population, the resulting $N^e$ can be biased (Waples and Yokota 2007). In these cases, as suggested by Waples and Yokota (2007), larger spacing between samples maximizes the drift signal compared to sampling biases associated with age structure.

Using a small number of generations can lead to outlier estimates

In general, $Ne(P)$ has a lower bias but larger variance, especially when t is small. As pointed out by Jorde and Ryman (2007) our weighting scheme leads to an increased variance but a smaller bias compared to other schemes. We observe outlier estimates among replicates at early generations (generation 5, Figure 2, Figure S4, Figure S5, and Figure S6) for $Ne(P).$ When the sampling variance is large compared to the drift variance (⁠$Ne=1000,$  $S≤100,$ and $R=50,$  Figure S11 and Figure S12), the deviation of the outlier estimates from the true $Ne$ is particularly large. For a few cases, we even observe large negative estimates. Negative estimates, in general, can be interpreted as $Ne$ being infinity, that is, no evidence of genetic drift (Peel et al. 2013). In our simulations this is plausible when $Ne$ is large and t is small, such that drift has not had a large effect on the population allele frequencies yet. Note that the harmonic mean estimator (⁠$N^e*(P)$⁠) has smaller variance for large $Ne$ (Figure S16). This estimator, however, is less accurate than $Ne(P)$ for small $Ne$ as shown in Figure S17.

To eliminate potential outliers and an inflated variance we recommend increasing the signal-to-noise ratio by pooling a sufficient number of individuals. Using later generations or increasing the number of SNPs in the analysis also helps to avoid outlier estimates. When none of these strategies can be applied, we suggest using the genome-wide median $Ne(P)$ estimates or the harmonic mean estimator, as these are more robust to extreme outliers.

Conclusions

Effective population size is an important parameter for describing evolutionary dynamics, making its accurate estimation essential for population genetic studies. Several methods have been designed to estimate $Ne,$ and their performance was comprehensively evaluated on simulated as well as real data (Barker 2011; Serbezov et al. 2012; Baalsrud et al. 2014; Holleley et al. 2014; Gilbert and Whitlock 2015). These studies mainly focused on genetic data collected from natural populations, which usually differ from experimental studies in terms of the census population size and sampling scheme. We designed a method that accurately infers the effective population size in genome-wide data from experimental populations sequenced in pools. Our approach improves temporal methods by explicitly correcting for two stages of sampling introduced by pooling and sequencing. Our results on simulated data confirm that methods that fail to properly account for the two stages of sampling inherent to Pool-seq can lead to severely biased $Ne$ estimates.

Pool-seq data are often considered to be overdispersed, i.e., displaying more variability than is predicted by the binomial sampling model (Yang et al. 2012). However, Zhu et al. (2012) and Futschik and Schlötterer (2010) validated that the error in allele frequency estimates is reasonably well approximated by binomial sampling given that a large enough number of individuals are pooled. Nevertheless, if overdispersion is present in the data, that will lead to additional variance, which is not modeled in our framework and will result in a downward bias of the estimated $Ne.$ If the level of overdispersion can be inferred for the data (see, e.g., Gautier et al. 2013; Illingworth 2015), it is possible to introduce a parameter that accounts for the additional between-pool variation (see File S1, Equation S8).

We also illustrate the applicability of our method for estimating $Ne$ from experimental data of D. melanogaster and show that in combination with a recursive partitioning method we can infer patterns of local variation in $Ne$ along the genome. Additionally, it is possible to calculate confidence intervals based on the $χ2$ distribution (Waples 1989) or alternatively apply a nonparametric bootstrap approach.

Software availability

Our proposed estimators along with standard methods from the literature are implemented within the R package Nest. The package is currently available at https://github.com/ThomasTaus/Nest.

Acknowledgments

We thank Mads Fristrup Schou, Susanne U. Franssen, and Neda Barghi for helpful comments on the Nest software package and Robin S. Waples and an anonymous reviewer for their constructive comments that greatly improved the manuscript. A.J. and T.T. are members of the Vienna Graduate School of Population Genetics, which is funded by the Austrian Science Fund (FWF, W1225). C.S. is also supported by the European Research Council grant ”ArchAdapt,” and T.T. is a recipient of a Doctoral Fellowship (DOC) of the Austrian Academy of Sciences.

Footnotes

Literature Cited