Abstract
A new genetic estimator of the effective population size (Ne) is introduced. This likelihood-based (LB) estimator uses two temporally spaced genetic samples of individuals from a population. We compared its performance to that of the classical F-statistic-based Ne estimator () by using data from simulated populations with known Ne and real populations. The new likelihood-based 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., Ne = 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 F-statistic 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 Ne 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 (Wright-Fisher model) in which genetic drift occurs at the same rate as in the studied population (Wright 1931). The effective population size (Ne) 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, Ne 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 Ne 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 Ne 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,” NeV (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 NeV can be conducted by moment-based methods (Waples 1989) and some recently published likelihood-based (LB) methods (Williamson and Slatkin 1999; Andersonet al. 2000).
LB estimators should provide better precision compared to moment-based 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 F-statistic estimator of NeV (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 Wright-Fisher model, in which the gene frequencies in the entire population in each generation are considered. We propose here a coalescent-based approach, which has different properties compared to that of Anderson et al. (2000). In many mating systems and life histories, including that of the Wright-Fisher model, the gene genealogy tends to the coalescent as the population size becomes large (Donnelly and Tavaré 1995; Mohle 2000). Thus the coalescent-based 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 Ne that do not follow the Wright-Fisher 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 Wright-Fisher model with small Ne sampled over a single generation.
The objectives of this article are twofold: (i) to present a new likelihood-based estimator of Ne, on the basis of coalescent theory, which allows the use of multiallelic markers and can be extended to incorporate Bayesian prior information about the Ne value (e.g., knowledge that Ne < 500) and (ii) to evaluate the accuracy and precision of this method in comparison to the classical estimator based on F-statistics (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 Ne 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 Ne, e.g., by setting NeMAX (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 Ne 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 Ne, 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).
Model used by the new likelihood-based estimator. We assume that samples are taken from a closed population at two different times. The population is represented by the dotted arrow. Following the conventions of genealogical modeling we take time to be increasing into the past, with the most recent sample (a0) taken at time t = 0 and the earlier sample (aT) taken at time t = T. a0 is assumed to have a genealogy, described by the standard coalescent model, and af represents its founders.
In materials and methods, we describe (i) the new likelihood-based estimator for Ne (noted ), (ii) the classical F-statistic-based 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
Likelihood-based estimator: The method used to estimate likelihoods for Ne 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, aT and a0, of nT and n0 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 k0 and kT for the samples taken at times 0 and T. The samples are assumed to be sampled independently with probabilities p(aT|nT, k, x) and p(a0|T, Ne, n0, k, x), where Ne 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 nT chromosomes are sampled with replacement from a population having unknown gene frequency x. Hence aT follows the multinomial distribution with probability p(aT|nT, 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 n0 lineages at t = 0, nc coalescent events between t = 0 and t = T, and nf = n0 − nc lineages at t = T. Since mutations are assumed not to occur, nf ≥ k0. Let af be the vector of counts of the nf 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(af|nf, x). It is possible to calculate the probability of obtaining nc coalescent events, p(nc|T, Ne, n0) (Tavaré 1984), which depends on the effective population size, Ne. In the case of fluctuating populations, Ne 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 Ne is measured as the harmonic mean Ne over the interval (Marjoram and Donnelly 1997).
Given the configuration in the founder lineages, af, it is possible to calculate the probability of obtaining the configuration in the sample, p(a0|af, n0, nc), by nc successive iterations of choosing a lineage uniformly at random and duplicating it (Slatkin 1996). Thus by enumerating all possible af for each value of nf = n0 − nc, nc = 0 … n0 − k0 it is possible to calculate the probability of obtaining the sample configuration at time 0 as the sum
We are interested in estimating p(a0|T, Ne, n0, k, x). Rewriting (1) in a more general way,
Equation 2 has the general form p(x) = Σyp(x|y)p(y), which can be estimated by the classical Monte Carlo method as
set S to 0;
do s times:
set Pi to 1; set nr to n0; set τ to t;
while (τ ≤ T/Ne and nr ≠ m) do:
with probability given by Equation 3, choose an allelic category in the current configuration and decrease the count by 1. Decrease nr by 1;
multiply Pi by the importance weight (Equation 4);
set τ to τ + t;
if (nr = m and τ ≤ T/Ne) then set Pi = 0 (probability of obtaining the data is zero with this number of coalescences) else
let afi be the current configuration (i.e., allele counts among founder lineages), containing nfi genes;
add to S the product of Pi by p(afi|x, nfi), the multinomial probability of choosing the current configuration from x;
evaluate S/s.
Thus the method estimates p(a0|T, Ne, n0, 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 Metropolis-Hastings 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 Ne 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 beta-distribution, as described in Ciofi et al. (1999) (although the Dirichlet method described in O'Ryanet al. 1998 also works well), and the likelihood
Data sets used to test for the convergence of the Metropolis-Hastings procedure embedded in NeLB
Convergence of the Metropolis-Hastings procedure: For the data sets described here, the Metropolis-Hastings 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 Gelman-Rubin criterion to assess convergence, implemented in CODA (Bestet al. 1995). The Gelman-Rubin statistic takes the values from the last one-half of the sampled points from the chains and estimates the square root of the ratio of the variance in Ne 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 Gelman-Rubin 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 Ne 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 Ne (i.e., the likelihood function tends to a nonzero constant for large Ne), 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 likelihood-based estimator .
F-statistic-based estimator: We estimated NeFk for each population using the equation
We estimated F as did Waples (1989) as the mean of Fk over loci. Fk was calculated for each locus from the equation
Other estimators of F give generally similar results to those obtained with Fk (Waples 1989; Richards and Leberg 1996; Luikartet al. 1999). Confidence intervals (90%) for Ne 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 individual-based 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 Wright-Fisher model, with two modifications: (i) separate sexes (with an equal sex ratio) and (ii) strict allogamy (no selfing). These modifications should lead to an Ne slightly larger than the actual number of breeders (Nb, 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 Ne because it is known that estimators of Ne give reasonably small confidence intervals when Ne is large and when using realistic sample sizes (e.g., 10–20 loci and 30–60 individuals; Luikartet al. 1999; Waples 2000). Because the Ne is usually <N (the census size) in natural populations (see Frankham 1995; Schwartzet al. 1998), it is realistic to model a population where Ne 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 Ne 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 Ne 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:
Allele frequency arrays used for the initiation of each locus for simulated data sets
Effective population size estimates and confidence intervals for empirical data sets
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 Ne value given in Table 3 is the harmonic mean of the Ne for the generations between samples. The Ne 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 Ne 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 Ne. 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 Ne = 500 in the likelihood-based 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 F-statistic-based 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 Ne 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 F-statistic-based 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 likelihood-based 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 Ne was below the lower limit and the proportion of times the true Ne was above the upper limit.
Ne estimates and confidence intervals for data simulated with true Ne = 10, 20, or 50 for one or five generations (T)
RESULTS AND DISCUSSION
Accuracy of estimates: The likelihood-based (LB) method appeared generally more accurate than the F-statistic-based 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 Ne, 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 Ne = 20 and T = 5,
overestimates Ne 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 Ne 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 Ne, whereas
does not (
was 191.7 and
was 83.3 when true Ne 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 Ne is small. For example, looking at the simulations where
over
.
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 Ne 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 Fk-based estimator of Ne(NeFk) and the 5–95th percentiles for the likelihood-based estimator of Ne(NeLB).
The credible intervals for the likelihood-based estimator of Ne() gave better precision than the F-statistic-based estimator (
) in the cases with high drift, whereas the difference in performance appeared variable (depending on Ne and on the number of loci) in the cases with relatively low drift. In the high drift cases, for example, when true Ne = 20 (and T = 5, L = 10, AF = A), the summaries of CIs were 9.7–53.9 for NeLB and 11.3–74.7 for
(Table 4). With real microsatellite data from a Drosophila population of known Ne (
(7.9–17.6) was also much smaller than that for
(21.57–80.75; Table 3). We note that the
When drift is weak, it is known that the lack of drift signal in allele frequency data leads to bad performance of the Ne 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 Ne = 10, T = 1,
and 5 loci (first line in Table 4), it is noteworthy that NeLB 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 Ne became smaller (from 50 to 20 to 10). This probably results from violation of the assumption in coalescent methods that Ne is large.
Percentage of confidence intervals lower (L) and higher (H) than the true Ne for 200 simulation estimates
Because our simulations used a fixed heterozygosity (H) level for all loci, we investigated if a variance of H among loci changes precision of Ne estimates. We used a set of five loci having a mean
Table 5 shows that the NeLB CIs are more often too low than too high (they are biased low), whereas NeFk 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 Ne was 18.8, the LB method credible interval was 7.9–17.6, whereas the F-statistic-based 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 Ne.
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 Ne in the case of the likelihood-based method. The first three should increase the signal-to-noise ratio in the data (Waples 1991).
Typical users wishing to increase the precision of Ne 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 Ne estimators we studied.
For improving precision, the total number of independent alleles (ΣiAi − 1, where Ai 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 NeFk and NeLB, 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 NeFk and NeLB 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 NeLB 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 Ne, as one gets a posterior distribution for Ne (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 signal-to-noise ratio (large if symmetrical, for example, Figure 3a; small if skewed or flat, for example, Figures 3, b and c).
Confidence intervals for NeLB estimates from independently simulated populations with the same effective population size (Ne = 20). All estimates used loci with five alleles, an initial heterozygosity of 0.6, and samples of 30 individuals separated by five generations. a was generated using only 5 loci, whereas b uses 10 loci. Dashed lines show the true Ne. (a and b) The single vertical bars represent the confidence intervals (the 5–95th percentiles of the posterior distribution) for 40 independent simulation estimates. The summary box chart at the right is made from 200 independent simulations. Ninety-five percent of the 200 upper support limits obtained are below the upper horizontal bar (see “95th percentile” arrow), and 95% of the 200 lower support limits obtained are above the lower horizontal bar. Similarly, 80% of the upper support limits are below the upper edge of the box, and 80% of the lower support limits are above the lower box edge.
The Bayesian aspect of using NeMAX = 500 (i.e., setting
a narrow prior probability distribution for Ne: NeMAX = 500 instead of 5000) provided smaller CIs when drift was relatively weak, but no obvious difference when drift was strong. For example, when Ne was 10 (and T = 1, with 10 loci), setting NeMAX to 5000 gave a summary of CIs of 3.6–3732.2 for percentiles, whereas it gave 3.6–130.6 with NeMAX = 500. This is a far better improvement than in the case of high drift (Ne = 20, T = 5, with 10 loci), for which case we obtained 6.3–107.8, when NeMAX was 5000 and 6.3–105.6 when NeMAX was 500. The better improvement observed under low drift situations was expected because likelihood curves obtained in these cases are skewed toward high Ne, and using the NeMAX 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 Ne, e.g., smaller NeMAX values.
Posterior distributions (likelihood curves) for samples from three independent simulated populations with (a and b) Ne = 20 and 5 loci with five alleles (H = 0.6) and with (c) Ne = 10 and 10 loci with two alleles (H = 0.2). All estimates use samples of 30 individuals separated by (a and b) five generations or by (c) one generation. Confidence intervals are often narrow (a and b), but the confidence intervals are occasionally very large (c). Please note that c uses different scales on the axis than a and b.
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 Ne 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 Ne because overestimation could lead managers to not detect a small Ne (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