LikelihoodBased Estimation of the Effective Population Size Using Temporal Changes in Allele Frequencies: A Genealogical Approach
 ^{*}Laboratoire de Biologie des Populations d'Altitude, UMR CNRS 5553, Université Joseph Fourier, F38041 BP53 Cedex 9 Grenoble, France
 ^{†}Computational and Molecular Population Genetics Laboratory, Zoologisches Institut, Universität Bern, 3012 Bern, Switzerland
 ^{‡}School of Animal and Microbial Sciences, University of Reading, Reading RG6 6AJ, United Kingdom
 ^{§}Laboratoire de Modélisation et de BiologieÉvolutive, INRAURLB, 34090 Montpellier, France
 1 Corresponding author: Computational and Molecular Population Genetics Laboratory, Zoologisches Institut, Universität Bern, Baltzerstrasse 6, 3012 Bern, Switzerland. Email: pierre.berthier{at}zoo.unibe.ch
Abstract
A new genetic estimator of the effective population size (N_{e}) is introduced. This likelihoodbased (LB) estimator uses two temporally spaced genetic samples of individuals from a population. We compared its performance to that of the classical Fstatisticbased N_{e} estimator () by using data from simulated populations with known N_{e} and real populations. The new likelihoodbased estimator () showed narrower credible intervals and greater accuracy than () when genetic drift was strong, but performed only slightly better when genetic drift was relatively weak. When drift was strong (e.g., N_{e} = 20 for five generations), as few as ~10 loci (heterozygosity of 0.6; samples of 30 individuals) are sufficient to consistently achieve credible intervals with an upper limit <50 using the LB method. In contrast, ~20 loci are required for the same precision when using the classical Fstatistic approach. The estimator is much improved over the classical method when there are many rare alleles. It will be especially useful in conservation biology because it less often overestimates N_{e} than does and thus is less likely to erroneously suggest that a population is large and has a low extinction risk.
THE effective size of a population can be defined as the size of an ideal population (WrightFisher model) in which genetic drift occurs at the same rate as in the studied population (Wright 1931). The effective population size (N_{e}) is an important parameter in evolution and conservation biology because it influences the amount of genetic drift in populations. Genetic drift influences the rate of loss of genetic diversity, the rate of fixation of deleterious alleles, and the efficiency of natural selection at maintaining beneficial alleles.
Unfortunately, N_{e} is notoriously difficult to estimate by demographic methods in the field (Frankham 1995), because it requires data such as variance in lifetime reproductive success, which is difficult to obtain for many wild populations. Genetic methods of estimating N_{e} are becoming more widely used because of the in creasing availability of polymorphic molecular markers (for reviews, see Waples 1991; Schwartzet al. 1998). However, a serious problem with existing genetic estimators of N_{e} is their poor precision; e.g., their confidence intervals often include infinity (Hill 1981; Waples 1991; Luikart and Cornuet 1999; Luikartet al. 1999).
The most widely used genetic method consists of measuring the variance of allele frequencies between generations to estimate the “variance effective size,” N_{eV} (Roberdset al. 1991; Hedgecocket al. 1992; Husband and Barrett 1992; Tayloret al. 1993; Burczyk 1996; Jorde and Ryman 1996; Miller and Kapuscinski 1997; Saadreva 1997; Sitnikovet al. 1997; Laikreet al. 1998; Planes and Lecaillon 1998; Tarret al. 1998; Fiumeraet al. 1999; Funket al. 1999; Kantanenet al. 1999). Given two genetic samples from one population, spaced by a known number of generations, estimation of N_{eV} can be conducted by momentbased methods (Waples 1989) and some recently published likelihoodbased (LB) methods (Williamson and Slatkin 1999; Andersonet al. 2000).
LB estimators should provide better precision compared to momentbased estimators, because they use more of the information from the data (Edwards 1972). For example, the Williamson and Slatkin (1999) LB method showed better precision and accuracy compared to the Fstatistic estimator of N_{eV} (Waples 1989). This initial study was restricted to diallelic loci because of the numerical complexity of using loci with a large number of alleles. However, Anderson et al. (2000) have developed an “importance sampling” approach that enables multiallelic loci to be analyzed. This approach is based on analysis of the WrightFisher model, in which the gene frequencies in the entire population in each generation are considered. We propose here a coalescentbased approach, which has different properties compared to that of Anderson et al. (2000). In many mating systems and life histories, including that of the WrightFisher model, the gene genealogy tends to the coalescent as the population size becomes large (Donnelly and Tavaré 1995; Mohle 2000). Thus the coalescentbased approach will give very similar answers to that of Anderson et al. (2000) when the population sizes are large. Since the method deals only with the genealogy and not with the population gene frequencies it is potentially more efficient than the method of Anderson et al. for larger population sizes. In addition, in populations with low N_{e} that do not follow the WrightFisher model, the gene genealogy might be better approximated by the coalescent. In this article we demonstrate that, in fact, the method performs satisfactorily with data generated from a WrightFisher model with small N_{e} sampled over a single generation.
The objectives of this article are twofold: (i) to present a new likelihoodbased estimator of N_{e}, on the basis of coalescent theory, which allows the use of multiallelic markers and can be extended to incorporate Bayesian prior information about the N_{e} value (e.g., knowledge that N_{e} < 500) and (ii) to evaluate the accuracy and precision of this method in comparison to the classical estimator based on Fstatistics (Krimbas and Tsakas 1971; Nei and Tajima 1981; Pollak 1983; Waples 1989), which uses the same genetic data from temporally spaced samples. Our evaluation is conducted empirically, using data from simulated populations with known N_{e} and from real populations. The model used to simulate populations is individual based, which provides realistic samples, complementary to the real data sets, for evaluating the usefulness of the estimators in natural populations.
The biggest problem with existing methods is their large confidence intervals (Waples 1991; Luikartet al. 1999). For example, the χ^{2} approximation method (used for ) is known to be slightly too conservative (i.e., confidence intervals are too wide; Waples 1989). Thus it is worth comparing the performance of the new method that we introduce to the widely used, “classical” method. We tested whether the precision of our LB method can be improved by incorporating Bayesian prior knowledge about N_{e}, e.g., by setting N_{eMAX} (the maximum possible effective size) to 500 or 5000. Like numerous LB estimators, the method we present is computationally intensive. On a Pentium II, 400 MHz PC, it requires 2–10 hr to estimate N_{e} for one data set with 10–20 loci and samples of 30 individuals. The program is slower with more loci or individuals. The slow speed is mainly of concern for studies like ours involving thousands of population replicates (typical users can simply run the program overnight to get a single estimation). The speed compares favorably with that reported by Anderson et al. (2000). Because of the slow speed, we could not evaluate many scenarios [i.e., combinations of N_{e}, sample size (S), time between the two samples (T), number of loci, and allele frequency distributions]. We tried to focus on realistic scenarios having the following characteristics: (i) multiallelic loci [e.g., five alleles with heterozygosity (H) of 0.6] because they provide more precise and accurate estimations (Luikartet al. 1999) and are becoming increasingly available for natural populations (e.g., microsatellites); (ii) sets of 5–20 unlinked >loci, typical of most field studies; and (iii) samples of 30 or 60 individuals separated by one or five generations (i.e., one or five episodes of genetic drift). However, we also conducted some evaluations with different allele frequencies (five alleles with equal frequencies, H = 0.8) and diallelic loci (H = 0.2) that more closely correspond to allozyme loci or single nucleotide polymorphisms (SNPs).
In materials and methods, we describe (i) the new likelihoodbased estimator for N_{e} (noted ), (ii) the classical Fstatisticbased estimator (noted ), (iii) the model used to simulate data sets, (iv) the real data sets, and (v) the methods we used to assess the relative performance of both estimators.
MATERIALS AND METHODS
Likelihoodbased estimator: The method used to estimate likelihoods for N_{e} given the data is based on the genealogical approach described in O'Ryan et al. (1998) and Beaumont and Bruford (1999), which is very similar to that described in Nielsen et al. (1998) and Saccheri et al. (1999; see also Ciofiet al. 1999; Chikhiet al. 2001). The model is illustrated in Figure 1. We assume that samples are taken from a closed population at two different times. In principle the method can be easily extended to deal with samples taken at many times, but only a pair of samples is considered here. Following the conventions of genealogical modeling we take time to be increasing into the past, with the most recent sample taken at time t = 0 and the earlier sample taken at time t = T. For a particular locus, the two samples at times t = T and t = 0 consist of two vectors of counts, a_{T} and a_{0}, of n_{T} and n_{0} chromosomes distributed among k distinct alleles. The number of distinct alleles is taken to be that observed in both samples combined; each individual sample may have <k alleles, and this number is denoted k_{0} and k_{T} for the samples taken at times 0 and T. The samples are assumed to be sampled independently with probabilities p(a_{T}n_{T}, k, x) and p(a_{0}T, N_{e}, n_{0}, k, x), where N_{e} is the effective population size and x is the parametric gene frequency at time T. The samples at each of m loci are also assumed to be sampled independently and therefore an overall likelihood
Probability of the first sample: It is assumed that the n_{T} chromosomes are sampled with replacement from a population having unknown gene frequency x. Hence a_{T} follows the multinomial distribution with probability p(a_{T}n_{T}, k, x).
Probability of the second sample: The sample at t = 0 is assumed to have a genealogy, described by the standard coalescent model, with n_{0} lineages at t = 0, n_{c} coalescent events between t = 0 and t = T, and n_{f} = n_{0} − n_{c} lineages at t = T. Since mutations are assumed not to occur, n_{f} ≥ k_{0}. Let a_{f} be the vector of counts of the n_{f} distinct lineages that remain at t = T, distributed among the k allelic classes. Those lineages are sampled with replacement from the population with probability p(a_{f}n_{f}, x). It is possible to calculate the probability of obtaining n_{c} coalescent events, p(n_{c}T, N_{e}, n_{0}) (Tavaré 1984), which depends on the effective population size, N_{e}. In the case of fluctuating populations, N_{e} is the harmonic mean effective population size over the interval T as discussed in O'Ryan et al. (1998). The distribution of coalescent events in an interval, and hence the joint likelihood of sample configurations, is identical between all models of population change, including stable populations, as long as N_{e} is measured as the harmonic mean N_{e} over the interval (Marjoram and Donnelly 1997).
Given the configuration in the founder lineages, a_{f}, it is possible to calculate the probability of obtaining the configuration in the sample, p(a_{0}a_{f}, n_{0}, n_{c}), by n_{c} successive iterations of choosing a lineage uniformly at random and duplicating it (Slatkin 1996). Thus by enumerating all possible a_{f} for each value of n_{f} = n_{0} − n_{c}, n_{c} = 0 … n_{0} − k_{0} it is possible to calculate the probability of obtaining the sample configuration at time 0 as the sum
We are interested in estimating p(a_{0}T, N_{e}, n_{0}, k, x). Rewriting (1) in a more general way,
Equation 2 has the general form p(x) = Σ_{y}p(xy)p(y), which can be estimated by the classical Monte Carlo method as
set S to 0;
do s times:
set P_{i} to 1; set n_{r} to n_{0}; set τ to t;
while (τ ≤ T/N_{e} and n_{r} ≠ m) do:
with probability given by Equation 3, choose an allelic category in the current configuration and decrease the count by 1. Decrease n_{r} by 1;
multiply P_{i} by the importance weight (Equation 4);
set τ to τ + t;
if (n_{r} = m and τ ≤ T/N_{e}) then set P_{i} = 0 (probability of obtaining the data is zero with this number of coalescences) else
let a_{f}_{i} be the current configuration (i.e., allele counts among founder lineages), containing n_{f}_{i} genes;
add to S the product of P_{i} by p(a_{f}_{i}x, n_{f}_{i}), the multinomial probability of choosing the current configuration from x;
evaluate S/s.
Thus the method estimates p(a_{0}T, N_{e}, n_{0}, k, x) with some error depending on the number of iterations, and the number we actually used is discussed below.
Posterior distribution: Assuming independence the probabilities can be multiplied across samples and loci as described above to give a likelihood
In MetropolisHastings sampling we propose candidate values x′ from some distribution conditional on the current value x, p(x′x) and calculate the ratio
Candidate values for N_{e} are proposed, using a lognormal distribution centered around the current value. Candidate values for x for each locus are proposed by randomly partioning the alleles into two groups and using a betadistribution, as described in Ciofi et al. (1999) (although the Dirichlet method described in O'Ryanet al. 1998 also works well), and the likelihood
Convergence of the MetropolisHastings procedure: For the data sets described here, the MetropolisHastings procedure was run for single chains of 20,000 iterations, sampling every 5 iterations. The validity of this number was assessed for four representative data sets, by running five independent chains from different starting points for each set, and using the GelmanRubin criterion to assess convergence, implemented in CODA (Bestet al. 1995). The GelmanRubin statistic takes the values from the last onehalf of the sampled points from the chains and estimates the square root of the ratio of the variance in N_{e} when the chains are combined to the average of the variance within each chain. The idea is that initially, with independent starting points, the variance across chains will be substantially greater than the variance within, reflecting the different starting points. However, when the chains are run long enough the variances should converge and the ratio should tend to 1. Gelman (1996) suggests as a “rule of thumb” that the value of the ratio of variances should be <1.1–1.2 (i.e., value of the statistic <1.05–1.1). The method also calculates the upper 97.5% credible limit for the statistic. The four data sets used and the results obtained for the GelmanRubin criterion are presented in Table 1. These results show that we can be reasonably confident that error due to Markov chain Monte Carlo (MCMC) estimation is a small fraction of the variability in estimation of N_{e} across data sets.
Point estimates and confidence intervals: From the output of the MCMC simulations the following summary statistics are estimated: the mode and 0.05 and 0.95 quantiles (giving a 90% credible interval). The 0.05 and 0.95 quantiles are obtained from the ranked output in the standard way. The mode is obtained by kernel density estimation. We used a logistic function as a kernel with a bandwidth (standard deviation) given by
Many data sets may have little information on the upper limit of N_{e} (i.e., the likelihood function tends to a nonzero constant for large N_{e}), in which case the 90% credible interval depends strongly on the prior used. We use a rectangular prior in the study here, but for real data a more smoothly varying prior may reflect background information more closely. This estimator was implemented in a computer program (TM3) available from http://www.rubic.rdg.ac.uk/~mab/software.html. A program for converting GENEPOP files to input format for the TM3 program is also available. The mode estimated from the MCMC output is referred to as the likelihoodbased estimator .
Fstatisticbased estimator: We estimated N_{eFk} for each population using the equation
We estimated F as did Waples (1989) as the mean of F_{k} over loci. F_{k} was calculated for each locus from the equation
Other estimators of F give generally similar results to those obtained with F_{k} (Waples 1989; Richards and Leberg 1996; Luikartet al. 1999). Confidence intervals (90%) for N_{e} were computed using a χ^{2} approximation, known to be unbiased, but too conservative (Waples 1989; Luikartet al. 1999).
Model for simulating data sets: The simulated populations were generated using an individualbased model with Mendelian inheritance. Populations were initiated by randomly sampling alleles from independent loci (defined by an array of allelic frequencies; see Table 2). Simulations were based on a WrightFisher model, with two modifications: (i) separate sexes (with an equal sex ratio) and (ii) strict allogamy (no selfing). These modifications should lead to an N_{e} slightly larger than the actual number of breeders (N_{b}, individuals), as shown in Caballero (1994, Equation 16). We verified this fact by estimating the variance of reproductive success,
We studied only cases with relatively small N_{e} because it is known that estimators of N_{e} give reasonably small confidence intervals when N_{e} is large and when using realistic sample sizes (e.g., 10–20 loci and 30–60 individuals; Luikartet al. 1999; Waples 2000). Because the N_{e} is usually <N (the census size) in natural populations (see Frankham 1995; Schwartzet al. 1998), it is realistic to model a population where N_{e} is only 10 or 20 and to sample 30 or more individuals. This model was implemented in a computer program available from pierre. berthier{at}zoo.unibe.ch.
Real data: To test further the N_{e} estimators, we applied them to real data, which potentially follow a model very different from the ones used both by the estimators and by our simulation model. For example, real data can show migration (killifishes), selection and linkage among loci (Drosophila; see discussion), overlapping generations (otter), and nonstable population size (mosquito fish). Our purpose here is not to use these data sets to quantify the problems introduced by those violations of the model underlying the N_{e} estimators, but rather to illustrate that the violation of assumptions can lead to changes in the estimators' performance. We used microsatellite data from four animal populations:
Spanish killifish Aphanius iberus: A captive population was founded in 1994 using 83 individuals from a wild population and a subsequent incorporation of 68 new individuals occurred 1 year later from the wild into the captive population (S. Schönhuth, unpublished data). The first sample was taken among the 83 founders (originating from the wild population) and the second sample was taken from the captive population after the incorporation of the 68 new individuals.
Two Drosophila melanogaster populations (1spm24 and 100a), which originated from the wild (at the Tyrrells Winery near Sydney, Australia), were genotyped (first sample); submitted to a bottleneck for a duration of 1 and 57 generations for, respectively, 1spm24 and 100a (see Table 3); and then genotyped a second time 10 generations later, during a rapid expansion period to a population size of 500 and 750 individuals, respectively. The N_{e} value given in Table 3 is the harmonic mean of the N_{e} for the generations between samples. The N_{e} for each generation is estimated from the pedigrees (England 1998).
Eurasian otter Lutra lutra (Dallaset al. 1999), for which samples from 1983 to 1988 were pooled to represent one sample, as were samples from 1991 to 1997.
Two western mosquito fish (Gambusia affinis) populations (A5 and B4), established experimentally by a bottleneck of one pair and eight pairs, respectively, from a wild population “source” (Spenceret al. 2000). For each of these two populations, the first temporal sample was taken in the source population. Experimental populations were then allowed to expand after the founder event and the second temporal sample was taken two to three generations later. Female founders were obtained from a general laboratory stock and may have not been virgin at the release time, so that the N_{e} at establishment could be underestimated by the actual number of founders.
Comparison of estimators: The first goal of this article is to describe a new estimator for N_{e}. The second goal is to provide an evaluation of its performance in comparison with a widely used estimator, so that researchers can use it with some confidence, on the basis of realistic examples. We take a frequentist approach and use simulations to compare the two methods. Although this is inconsistent with the Bayesian paradigm behind the new estimator, it is the most practicable way to compare the two approaches. There are two issues to consider here. Having an upper limit of N_{e} = 500 in the likelihoodbased method will in itself increase the accuracy of the method. However, this effect is generally small, and cases where it has an effect are clear in the results and discussion. In addition, the credible intervals from the likelihood estimator are a priori not expected to have the same coverage properties as the confidence intervals from the Fstatisticbased estimator—the former is a fixed interval, and the “truth” is a random variable, while the latter is a random interval, and the truth is fixed. However, since the assessment is frequency based (where the truth is fixed) it seems reasonable to assess the coverage properties of the two estimators, while accepting that they mean different things. It is also worth noting that a general result in classical theory is that the credible intervals and confidence intervals will converge to have the same coverage properties, providing that the posterior distribution is asymptotically normal (see, e.g., Gelmanet al. 1995, Chap. 4), and this is likely to be the case here with increasing numbers of loci (but not sample size).
In assessing the performance of each estimator we attempted to answer the following questions: (i) Are the point estimates biased?, (ii) how accurate are the two estimators?, (iii) how large can we expect the confidence intervals to be?, and (iv) how often do the confidence intervals contain the true value? We answer these questions by using simulated data sets and compare the performance of the estimators using summary statistics computed from the set of replicates obtained with each tested scenario. The parameters investigated include the true N_{e} value, the sample size (S), the number of loci (L), the allele frequencies in the population from which is taken the first sample (AF; Table 2), and the number of generations between samples (T). Each scenario (i.e., parameter set) was evaluated using 200 populations simulated independently (i.e., replicates).
The bias of the point estimates was investigated by reporting the median of the 200 point estimates from the simulations. The median was used rather than the mean because the distribution of point estimates can be strongly skewed (for example, in the case of the Fstatisticbased estimator, there are occasionally point estimates of infinity). We also report the standard error of the mean. The accuracy of the point estimators was investigated by reporting the square root of the mean of the squared differences between the estimation values and the true value (the square root of the mean square error,
The properties of the estimated confidence intervals (or credible intervals in the case of the likelihoodbased estimator) were assessed by reporting two types of summary statistic. First, to summarize the overall spread of the estimated intervals, we report the 5th percentile of the lower limits of the intervals (that is, 5% of the 200 confidence or credible intervals had lower bounds that were lower than this) and the corresponding 95th percentile of the upper limits. The interval defined by these two percentiles is a summary of the overall spread of confidence or credible intervals in the 200 simulations. Second, to summarize the coverage properties of the estimated intervals, we report the proportion of times the true N_{e} was below the lower limit and the proportion of times the true N_{e} was above the upper limit.
RESULTS AND DISCUSSION
Accuracy of estimates: The likelihoodbased (LB) method appeared generally more accurate than the Fstatisticbased method (Table 4). For example, in 13 out of 16 sets of simulations
On the basis of the medians of the estimates from each set of 200 simulations, patterns in the bias of the estimates are discernible in Table 4. In general, appears to systematically overestimate N_{e}, whereas appears to slightly underestimate it. The overestimation by is due mainly to the loss of alleles in early generations of the time period between the two samples (see Richards and Leberg 1996; Luikartet al. 1999). These earlier observations suggest that the bias is greater with increasing drift and is also greater when there are many rare alleles, which are likely to be lost through drift (Waples 1989; Luikartet al. 1999). Table 4 supports the first observation: With N_{e} = 20 and T = 5, overestimates N_{e} by up to 25%, whereas it generally overestimates by <10% in the other two sets of simulations. Interestingly, the degree of overestimation appears also to be affected by N_{e} or T, because the amount of drift is otherwise the same in these two sets of simulations. The Drosophila data in Table 3 support the sec ond observation on the effect of rare alleles on . Here substantially overestimates N_{e}, whereas does not ( was 191.7 and was 83.3 when true N_{e} was 100.8).
The underestimation by is probably due to the violation of the assumptions of the coalescent in the simulated data sets and appears to be strongest when N_{e} is small. For example, looking at the simulations where
Precision of estimates: Because the accuracy is generally good (see Table 4 and examples above), it is usually less problematic than the poor precision (see, for example, the N_{e} estimates for the otter data in Table 3). Thus, it is very informative to assess the relative performance of the methods used to compute confidence or credible intervals (CIs): the χ^{2} approximation (Waples 1989) for the F_{k}based estimator of N_{e}(N_{eFk}) and the 5–95th percentiles for the likelihoodbased estimator of N_{e}(N_{eLB}).
The credible intervals for the likelihoodbased estimator of N_{e}() gave better precision than the Fstatisticbased estimator () in the cases with high drift, whereas the difference in performance appeared variable (depending on N_{e} and on the number of loci) in the cases with relatively low drift. In the high drift cases, for example, when true N_{e} = 20 (and T = 5, L = 10, AF = A), the summaries of CIs were 9.7–53.9 for N_{eLB} and 11.3–74.7 for (Table 4). With real microsatellite data from a Drosophila population of known N_{e} (
When drift is weak, it is known that the lack of drift signal in allele frequency data leads to bad performance of the N_{e} estimators studied here (e.g., Luikartet al. 1999). In those weak drift cases, gave more reliable confidence intervals than , but both estimators are similar when at least ~20 loci are used (see Table 4). When considering the scenario with N_{e} = 10, T = 1, and 5 loci (first line in Table 4), it is noteworthy that N_{eLB} performed better when doubling the sample size (second line) than when doubling the number of loci used (third line). On the other hand, showed similar benefits from these two possible strategies for increasing the sampling effort. It is also important that the CIs performed slightly worse as N_{e} became smaller (from 50 to 20 to 10). This probably results from violation of the assumption in coalescent methods that N_{e} is large.
Because our simulations used a fixed heterozygosity (H) level for all loci, we investigated if a variance of H among loci changes precision of N_{e} estimates. We used a set of five loci having a mean
Table 5 shows that the N_{eLB} CIs are more often too low than too high (they are biased low), whereas N_{eFk} CIs are more often too high than too low (they are biased high). The same pattern appears with real data (see Table 3); for example, with Drosophila populations, when the actual N_{e} was 18.8, the LB method credible interval was 7.9–17.6, whereas the Fstatisticbased method CI was 21.6–80.7. This is interesting for conservation biology applications because overestimation is likely to delay the detection of a small N_{e}.
Guidelines: To improve precision, one has four possibilities: (i) Use more loci or loci with more alleles, as this provides more independent observations of drift; (ii) sample more individuals, as it gives better allele frequency estimates; (iii) increase the time between samples, as this increases the number of drift events observed (T); and (iv) use prior knowledge on N_{e} in the case of the likelihoodbased method. The first three should increase the signaltonoise ratio in the data (Waples 1991).
Typical users wishing to increase the precision of N_{e} estimators would like to know whether increasing the number of individuals sampled or the number of loci is more likely to help. Indeed, new promising molecular markers (e.g., SNPs) may allow use of dozens of loci in the near future and could double or triple the number of loci typically used today (5–20). However, the cost and time needed for developing them might counterbalance the cost of sampling more individuals while keeping classical genetics markers. This is why performance analyses like ours are needed to help make such decisions about sampling investment. It is difficult to find dozens of unlinked loci in most real populations, and a linkage disequilibrium between loci is possible in small populations. Thus, more study is needed to quantify the potential problems caused by statistical association between loci on the N_{e} estimators we studied.
For improving precision, the total number of independent alleles (Σ_{i}A_{i} − 1, where A_{i} is the number of alleles at the locus i) is more important than the number of loci. For example, using five diallelic loci (thus 5 independent alleles) leads to infinite CIs with both N_{eFk} and N_{eLB}, whereas using five loci with 5 alleles each (20 independent alleles) provided far smaller CIs (Table 4, lines 6 and 10, for example). In this study, adding more loci or more individuals in the sample dramatically increased the precision on N_{eFk} and N_{eLB} estimates (Figure 2), but doubling the number of loci did not lead to obviously better results than doubling the number of individuals (for loci having 5 alleles), except for N_{eLB} in the cases of low drift. Our results provide rough ideas about this issue when using realistic data sets.
It is interesting that increasing the number of alleles in the samples by having more loci with rare alleles is now less problematic with than with . For example, is less biased and more precise than when using loci with some rare alleles (e.g., scheme A in Table 2). This is important because most alleles in natural populations are at low frequency except in severely bottlenecked populations (Luikartet al. 1998).
Bayesian methods give more direct information than classical statistics about the estimated parameter N_{e}, as one gets a posterior distribution for N_{e} (Figure 3). One advantage of a posterior distribution is that it allows the visualization of the information brought by the data. Thus, a flat posterior means that the data contain no information in connection with the model used: See, for example, Figure 3c. The symmetry of the curve allows us to roughly estimate the signaltonoise ratio (large if symmetrical, for example, Figure 3a; small if skewed or flat, for example, Figures 3, b and c).
The Bayesian aspect of using N_{eMAX} = 500 (i.e., setting a narrow prior probability distribution for N_{e}: N_{eMAX} = 500 instead of 5000) provided smaller CIs when drift was relatively weak, but no obvious difference when drift was strong. For example, when N_{e} was 10 (and T = 1, with 10 loci), setting N_{eMAX} to 5000 gave a summary of CIs of 3.6–3732.2 for percentiles, whereas it gave 3.6–130.6 with N_{eMAX} = 500. This is a far better improvement than in the case of high drift (N_{e} = 20, T = 5, with 10 loci), for which case we obtained 6.3–107.8, when N_{eMAX} was 5000 and 6.3–105.6 when N_{eMAX} was 500. The better improvement observed under low drift situations was expected because likelihood curves obtained in these cases are skewed toward high N_{e}, and using the N_{eMAX} prior information is likely to remove a greater amount of the area under the curve (i.e., in the right tail of the curve) than in the high drift cases. Further research is needed to study the influence of using different prior distributions for N_{e}, e.g., smaller N_{eMAX} values.
In the context of management the Bayesian methodology presented here may be further refined to incorporate the methods of decision theory. For example, a loss function can be defined, which enables the posterior risk of management decisions based on N_{e} to be quantified. Nevertheless, for conservation biology purposes, the method should already be useful as it gives narrower credible limits, with less chance to be biased high than low. In conservation biology, it is critical to not overestimate N_{e} because overestimation could lead managers to not detect a small N_{e} (and an associated excessive loss of genetic variation) and to consequently underestimate the population extinction risk.
Acknowledgments
We sincerely thank P. England for the Drosophila data, S. Schönhuth for the killifish data, J. Dallas for the otter data, and C. Spencer for the mosquito fish data. Many thanks to E. Anderson, M. Schwartz, L. Excoffier, N. Valière, S. Manel, and M. Gaudeul for their helpful comments on the manuscript. P.B. was funded by a DEA grant from French Research Ministry. M.A.B. is supported by a University of Reading RETF fellowship. The Bureau des Ressources Génétiques (INRA) provided funding for G.L. and J.M.C. G.L. was also supported by a grant from the European Union (BIO4 CT96 1189).
Footnotes

Communicating editor: G. A. Churchill
 Received October 15, 2000.
 Accepted October 19, 2001.
 Copyright © 2002 by the Genetics Society of America