The effective population size (Ne) is frequently estimated using temporal changes in allele frequencies at neutral markers. Such temporal changes in allele frequencies are usually estimated from the standardized variance in allele frequencies (Fc). We simulate Wright-Fisher populations to generate expected distributions of Fc and of F̅c (Fc averaged over several loci). We explore the adjustment of these simulated F̅c distributions to a chi-square distribution and evaluate the resulting precision on the estimation of Ne for various scenarios. Next, we outline a procedure to test for the homogeneity of the individual Fc across loci and identify markers exhibiting extreme Fc-values compared to the rest of the genome. Such loci are likely to be in genomic areas undergoing selection, driving Fc to values greater (or smaller) than expected under drift alone. Our procedure assigns a P-value to each locus under the null hypothesis (drift is homogeneous throughout the genome) and simultaneously controls the rate of false positive among loci declared as departing significantly from the null. The procedure is illustrated using two published data sets: (i) an experimental wheat population subject to natural selection and (ii) a maize population undergoing recurrent selection.
THE effective population size (Ne), defined as the size of an ideal Wright-Fisher population undergoing the same rate of genetic change as the population under study, is an essential parameter to predict the evolution of a population due to genetic drift in terms of rates of loss of genetic variation, fixation of deleterious alleles, or inbreeding (Wright 1969). However, obtaining direct estimates of Ne from demographic data has often proved difficult. An alternative is to use indirect methods, for instance, those based on the measurement of temporal changes in allele frequencies at neutral markers (Krimbas and Tsakas 1971; Waples 1989a). The foundation of these methods is that the variance of allele frequency due to drift from parents to offspring, V(P1), depends on Ne as follows: V(P1) = P0(1 − P0)/2Ne, where P0 is the frequency in the parental population. After t generations of drift, the expected frequency of the allele is E(Pt) = P0 and the variance of the allele frequency, V(Pt) = E(Pt − P0)2, can be written as a function of Ne: V(Pt) = P0(1 − P0)[1 − (1 − 1/2Ne)t] (Crow and Kimura 1970). If t is not too large (t ≪ Ne), Ne can be approximated by Ne ≅ (P0(1 − P0)t)/(2V(Pt)) and therefore an estimator for Ne based on the standardized variance in allele frequency is (V(Pt))/(P0(1 − P0)). Nei and Tajima (1981) proposed estimating the standardized variance in allele frequency between generation tx and ty for each locus l with Kl alleles as 1where px(i,l) [respectively py(i,l)] represents the frequency of allele i at locus l in the sample of Sx individuals drawn at generation tx (respectively Sy individuals at ty). A weighted mean of F̂c,l-values across several loci, 2is then typically used to estimate Ne via 3
(Waples 1989a). Note that Equation 1 assumes that alleles frequencies are estimated from samples taken prior to reproduction (so-called plan II sampling scheme). We use that sampling in the remainder of this article; allowing for an alternative sampling scheme is straightforward.
Recently, renewed interest in estimating Ne has led to the development of numerous methods using allele frequencies observed in a series of temporally spaced samples of a population. Williamson and Slatkin (1999) and Anderson et al. (2000) introduced a maximum-likelihood approach to estimate Ne. Berthier et al. (2002) proposed a coalescent-based likelihood approach to estimate Ne and Wang (2001) devised a faster approximate version using a pseudo-maximum-likelihood method. These methods, tested by the authors for some values of population parameters (Ne, l, Kl, and px(i,l)), have proved to be slightly more accurate than the Fc method, but are much more computationally intensive. Hence Fc-based estimators of Ne remain frequently used in practice (Fujiio et al. 1999; Turner et al. 1999; Goldringer et al. 2001; Shikano et al. 2001). Properties of Fc-based estimators of Ne and the quality of the confidence intervals around such estimates depend critically on the distribution of F̅c-values. Confidence intervals around Ne have been based on the fact that , with n = ∑l(Kl − 1), is distributed approximately as a chi-square with n d.f. (Lewontin and Krakauer 1973). Hence, assessing the adjustment of the actual F̅c distribution to a chi-square distribution is important. The chi-square approximation has been studied for some special cases, but the effects of initial allele frequencies, of the number of alleles, of the number of loci, and of the number of generations as well as the “true” effective population size on the distribution of F̅c-values and on N̂e are still poorly known (Waples 1989a).
Before averaging estimates of Fc obtained at individual loci to obtain F̅c and an estimate of Ne, it is desirable to test whether all loci used for that study have experienced the same effective population size. This implicit assumption, which underlies all methods of estimation mentioned above, is rarely tested. Several factors can modify the local effective size at a given locus. The recurrent elimination of deleterious variants linked to a marker locus, known as “background selection,” will reduce the effective size locally; this effect depends on the local recombination rates and on genome-wide parameters describing spontaneous mutation and their effect on fitness (Charlesworth et al. 1993). Hitchhiking will also drive higher than expected the temporal variance in allele frequency of markers linked to a positively selected variant (Wiehe and Stephan 1993).
The remainder of this note is organized as follows. In the first part, we study the actual F̅c distribution and the quality of the chi-square approximation. The actual distribution of F̅c, its divergence from a chi-square distribution, and the quality of the Ne estimation based on Fc are studied under various scenarios by varying the initial allele frequencies, the number of loci, the number of generations, the sample size at both generations, and the true effective population size. In the second part, we outline a procedure to identify loci with “extreme” individual F̂c,l-values. We illustrate our approach by reanalyzing two experimental data sets: temporal variation in allele frequencies at 29 markers in an experimental wheat population under natural selection and frequencies at 82 markers in a maize population under recurrent selection.
The actual distribution of F̅c and its consequences for estimating Ne:
To investigate the actual distribution of F̅c, we simulated Wright-Fisher populations using an exact multinomial sampling scheme. We generated expected distributions of temporal variations in allele frequencies conditional on initial allele frequencies. Distributions were based on 3000 independent replicates. In each replication, several loci with the same initial allele frequencies were simulated and F̅c was computed. All simulations were carried out using Mathematica (Wolfram 1996). Our simulations show that substantial departure of the actual distribution of F̅c from a chi-square distribution, as measured through Kullback's (1968) symmetric measure of divergence between both distributions (see Table 1), can be observed under a variety of conditions depending on the parameter values chosen (Table 1). The actual F̅c distribution is closest to a chi-square and thus Ne is best estimated when biallelic marker loci with equal frequencies are used. Conversely, the discrepancy between the actual F̅c distribution and the chi-square approximation is large when allele frequencies are strongly unbalanced (P0 < 0.1 for at least one allele), when the number of alleles per locus is large (K ≥ 5) such as for microsatellite markers, or when the number of generations increases (δT = ty − tx > 15 when Ne = 100 is assumed). Increasing the sample sizes up to 200 or 500 individuals does not diminish the discrepancy, especially for a high number of alleles (data not shown). In most cases, the distribution of simulated values is shrunk compared to the chi-square approximation. In addition, the actual distribution is a bit skewed toward higher values of Fc. As a consequence, the distribution of N̂e in the simulations is often more narrow than the one based on the chi-square approximation. Confidence intervals at the 95% level based on either the chi-square approximation or the actual F̅c distribution are given in Table 1. These are exactly the intervals that would be computed in experimental studies. Chi-square confidence intervals are often wider than the confidence intervals based on the simulated F̅c (Table 1). The width of this interval, which is connected to the precision of the estimation, depends mainly on the number of independent alleles used [L(K − 1)]. Note that the product L(K − 1) is also the number of degrees of freedom of the chi-square used in previous approximations. Confidence intervals derived from simulated data are reduced by 10–25% relative to chi-square-based confidence intervals for scenarios involving unbalanced initial allele frequencies with 5 ≤ L(K − 1) ≤ 10, or very large number of alleles (K = 10), sample sizes <100, δT > 10, or Ne < 75. Otherwise, they are quite close to the chi-square-based intervals (reduction is <10%). Similar results are found with larger sample sizes (200 and 500 individuals).
Except for δT = 5 generations, N̂*e, computed from F̅*c (an average of F̅c over 3000 independent replicates), is always higher than the population size Ne used in simulations. This indicates that F̅c-based estimation tends to return overestimated values of Ne. Richards and Leberg (1996) and Luikart et al. (1999) argued that the overestimation of Ne using F̅c [or Pollak's (1983) estimator, F̂k] is due mainly to the loss of alleles in early generations, suggesting that the bias would be greater with increasing drift and when there are rare alleles. Luikart et al. (1999) focused on the estimation of Ne after very strong bottlenecks (Ne = 4–40), which were not considered here. Whereas our results confirm the existence of a greater bias for rare alleles, for the range of Ne we considered, population size does not appear as critical; however, increasing the number of generations between samples leads to overestimation of Ne. Hence, to improve precision on the Fc-based estimation of Ne, we recommend generating the actual F̅c distribution using simulations based on the estimated value N̂e and obtaining a confidence interval directly on N̂e. With such a gain in accuracy, the performance of the Fc-based estimator of Ne becomes close to those of likelihood-based estimators.
Distribution of Fc at individual loci and the detection of loci departing from pure drift:
Once Ne has been estimated on the basis of F̅c by averaging over several marker loci, it is desirable to test whether all loci considered have undergone the same rate of change in allele frequency. Heterogeneity in individual Fc-values might be used as evidence for selection since “while natural selection will operate differently for each locus and each allele at a locus, the effect of breeding structure (migration, genetic drift, inbreeding) is uniform over all loci” (Lewontin and Krakauer 1973, pp. 176–177). Loci with significantly high F̂c,l-values should be discarded before (re)computing F̅c to yield a more reliable estimate of Ne.
We propose a way to identify, in a series of experimental F̂c,l measurements, markers exhibiting F̂c,l-values significantly higher than expected under pure drift based on N̂e. To do so, each F̂c,l-value should be compared to an expected distribution based on a genome-wide effective size estimated from the remaining loci and the trajectory of allele frequencies at this locus. We exemplify below our method with two published data sets. First we consider individual F̂c,l-values estimated from temporal variations of allele frequencies in an experimental composite wheat population undergoing natural selection. A total of 250 and 213 individuals were sampled at generations 1 and 10, respectively, and genotyped at 29 RFLP loci (Goldringer et al. 2001). For each locus l, we pooled the remaining 28 loci to obtain a global F̅c estimate and an F̅c-based estimate of Ne (described hereafter as the genome-wide average estimate). The expected distribution of Fc at locus l was then obtained (using typically 3000–5000 independent simulations) conditional on N̂e (excluding locus l) and the observed initial allele frequencies (at locus l). We then tested whether the observed temporal variance of allele frequencies at locus l, F̂c,l, was significantly larger than the genome-wide average variations by computing p, the probability for F̂c,l to be greater than or equal to the observed value at this locus on the basis of the simulated distribution described above. Note that one could also test for the presence of loci exhibiting smaller than expected variations in allele frequency. Some loci exhibited some “excess drift” relative to the rest of the loci and accordingly fairly small P-values: Fba242-C (P = 0.021), Fba280-C (P = 0.042), Fba65-D (P = 0.085), and Fba204-A (P = 0.09). However, the distribution of P-values was fairly uniform (data not shown) and to take into account the fact that multiple loci were examined we computed the expected false discovery rates, also known as q-values, using the distribution of P-values (see Storey and Tibshirani 2003 for details). The q-values were calculated using the package QVALUE (http://faculty.washington.edu/jstorey/qvalue/index.html). This analysis suggests that declaring only Fba242-C and Fba280-C as “significant” for excess of drift would still yield an estimated rate of false positives of ∼40% among these two loci.
Next we consider the study published by Labate et al. (1999), where temporal variations in allele frequencies were surveyed at 82 RFLP loci after 12 generations in maize populations undergoing recurrent selection. A P-value was calculated at each locus (Figure 1), using the method described above. The distribution of P-values was then used to calculate corresponding q-values (Figure 2). In contrast with the previous case, the distribution of P-values is clearly L-shaped (Figure 1) and choosing a q-value cutoff of 0.05 yields 10–11 loci exhibiting significant departures from the genome-wide level of drift. This method proved to be somewhat more conservative than the one used by the authors, who declared 14 loci as outliers (Labate et al. 1999). Discarding those outlier loci, one can compute a new genome-wide effective size and check that no more loci exhibit Fc-values departing significantly from the null hypothesis of homogenous drift (data not shown).
Waples (1989b) proposed a method based on the chi-square test of homogeneity to test the hypothesis that observed changes in allele frequencies can be satisfactorily explained by drift alone. This allows one to examine the variation of a particular allele according to the range of possible Ne-values for the population under study. Yet, this test has not been widely used in experimental studies. Indeed, it is rather complicated to implement, particularly in cases of multiple alleles since it is necessary to consider covariances of frequencies for different alleles sampled at different times. Lewontin and Krakauer (1973) provided the theoretical grounds for homogeneity tests of variation in allele frequencies, but they emphasized more spatial variation, and the test proposed for temporal variation, which again relies on the assumption of a chi-square distribution of individual Fc,l-values, is much too restrictive [see Beaumont and Nichols (1996) and Vitalis et al. (2001) for the case of spatial variation in allele frequencies]. The use of simulation-based distributions provides a robust method to test for homogeneity of F̂c,l-values across loci before pooling estimates. Our procedure yields a more reliable genome-wide estimate of the realized Ne and can be used to detect markers exhibiting Fc-values significantly higher than expected on the basis of mean Ne, thereby providing a (formal) way of assessing if selection is operating on any given genomic segment (see also Luikart et al. 2003 for a review of available methods for population structure). One potential caveat of our method is that the distribution of Fc under the null hypothesis is generated using information from the data (to estimate the genome-wide Ne). We verified through simulations (see online supplementary material at http://www.genetics.org/supplemental/) that our procedure is actually fairly robust to uncertainty in the estimation of the genome-wide Ne. A program generating the expected individual or mean Fc distributions used in this note is available upon request as a Mathematica notebook from the authors.
We thank F. Hospital, I. Bonnin, and C. Dillmann for helpful discussions and A. Tsitrone, R. Waples, and an anonymous reviewer for their comments on earlier versions of this article. We thank O. Martin for correcting and improving the English of this article.
Communicating editor: O. Savolainen
- Received December 17, 2003.
- Accepted June 7, 2004.
- Genetics Society of America