MATERIALS AND METHODS

Likelihood-based 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₀, of n_T and n₀ 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₀ 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₀|T, N_e, n₀, 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 L (T, N_e, x₁, x₂, … , x_m) can be obtained by multiplying the probabilities across samples and across loci. Metropolis-Hastings simulation is then used to integrate out the x_i as described below.

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₀ lineages at t = 0, n_c coalescent events between t = 0 and t = T, and n_f = n₀ − n_c lineages at t = T. Since mutations are assumed not to occur, n_f ≥ k₀. 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₀) (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₀|a_f, n₀, 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₀ − n_c, n_c = 0 … n₀ − k₀ it is possible to calculate the probability of obtaining the sample configuration at time 0 as the sum p(a0∣T,Ne,n0,k,x)=Σaf,ncp(a0∣af,n0,nc)p(af∣x,nf=n0−nc)×p(nc∣Ne,T,n0). (1) Although it is feasible to evaluate this sum for small samples, it is impractical for typical cases. We use here an alternative, efficient method of evaluation described by Griffiths and Tavaré (1994) and used by O'Ryan et al. (1998). For a more detailed discussion of the method see Griffiths and Tavaré (1994), Felsenstein et al. (1999), and Stephens and Donnelly (2000).

We are interested in estimating p(a₀|T, N_e, n₀, k, x). Rewriting (1) in a more general way, p(a0∣T,Ne,n0,k,x)=ΣGp(a0∣G)p(G∣af,nc)p(af∣x,nf=n0−nc)×p(nc∣T,Ne,n0). (2) The summation is over all genealogical histories G, each of which can be represented as a sequence of configurations G = (g₀, g₁, … , g_nc), where, looking back in time from t = 0, g₀ is the configuration prior to the first coalescent event and g_nc is the configuration prior to the sampling event at T (i.e., g_nc = a_f). The states (g_nc, … , g₁, g₀) form a Markov chain where the number of lineages increases by 1 at each step. Let b denote a vector of length k with 1 at position b and 0 elsewhere, b being a particular allelic class. In the following, gr(b) denotes the number of genes in the allelic class b within the rth state. The rth state in the sequence can be generated from Pr(g(r)=g(r+1)+b∣g(r+1))=g(r+1)(b)n(r+1), where n(r+1)=Σjg(r+1)(j) . Thus p(G|a_f, n_c) is the product of these probabilities over n_c coalescent events, and p(a₀|G) = 1 if g₀ = a₀, 0 otherwise.

Equation 2 has the general form p(x) = Σ_yp(x|y)p(y), which can be estimated by the classical Monte Carlo method as 1∕s∑i=1sp(x∣yi) , where the y_i are drawn from p(y) (see, e.g., Tanner 1993). For example, we could simulate the number of coalescent events to obtain n_f, simulate a sample from x, and simulate G. Unfortunately p(a₀|G) = 0 for most G. An alternative approach is to use importance sampling (see, e.g., Tanner 1993), where Σ_yp(x|y)p(y) is estimated as 1∕s∑i=1sp(x∣yi)p(yi)∕p∗(yi) , where p*(y) is chosen to be as close as possible to p(y|x); the variance due to importance sampling would be zero if it could be made equal to p(y|x)—this follows because p(x|y) p(y)/p(y|x) = p(y|x)p(x)/p(y|x) = p(x) for any draw from p(y|x). In general it is not possible to draw exactly from p(y|x) but this can be approximated quite closely. Following the method of Griffiths and Tavaré (1994), we can simulate a sample sequence with known probability as follows. For a configuration at event r, g_r, there are m prior configurations with one lineage less in each of the allelic classes, where m is the number of allelic classes in which there is at least one lineage. The configuration g_r−1 with one lineage less in allelic class b is chosen with probability (gr(b)−1)∕(nr−1)∑j=1m((gr(j)−1)∕(nr−1))=gr(b)−1nr−m, (3) where j subscripts are the allelic classes represented in g_r. The ratio p(G|a_f, n_c)/p*(G) is then given by the product of gr(b)−1nr−1/gr(b)−1nr−m=nr−mnr−1 (4) over all coalescent events in G. Thus to estimate (2) we can use the following algorithm, which calculates 1/s × S, with S=∑i=1sp(a0∣Gi)×Pi×p(af∣x,nf)p(nc∣T,Ne,n0) , and P_i = p(G_i)/p*(G_i). This algorithm involves two nested iterations: a first iteration, which builds the summation S over sampled genealogies, and another iteration, nested in the previous one, which calculates each P_i by sampling the number of coalescent events in the genealogy G_i. Let n_r be the number of lineages, m the number of allelic categories in the current configuration, and τ the total time elapsed since the start of the building of the current genealogy. The genealogies are built starting at τ = 0, with the data configuration a₀, so that for all i, p(a₀|G_i) = 1. Then coalescent events are simulated using a waiting time, t, until the total time elapsed becomes >T/N_e. In the following algorithm, t is a random variable, sampled from an exponential distribution with scale ((nr2) ), which changes each time we refer to it:

1. set S to 0;

2. do s times:

   1. set P_i to 1; set n_r to n₀; set τ to t;

   2. while (τ ≤ T/N_e and n_r ≠ m) do:

      1. 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;

      2. multiply P_i by the importance weight (Equation 4);

      3. set τ to τ + t;

   3. 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

      1. let a_fi be the current configuration (i.e., allele counts among founder lineages), containing n_fi genes;

      2. add to S the product of P_i by p(a_fi|x, n_fi), the multinomial probability of choosing the current configuration from x;

3. evaluate S/s.

Thus the method estimates p(a₀|T, N_e, n₀, 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 L (N_e, T, x₁, … , x_m). We integrate out the x_i using Metropolis-Hastings sampling, following the approach described in O'Ryan et al. (1998).

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 X=L(x′)∕p(x′∣x)L(x)∕p(x∣x′), (5) which gives the likelihood of a candidate value, weighted by the probability of choosing it from the current value, relative to the current value, weighted by the probability of choosing it from the candidate value. If X > 1 we accept x′; otherwise we accept it with probability X, retaining x otherwise. The likelihoods include a weighting by priors. Provided certain conditions hold (see Tanner 1993), the resulting sequence of accepted and retained values has a stationary distribution that is L(x)∕∫L(x)dx . In the case of jointly distributed variables marginal distributions can be estimated by looking at only one variable (e.g., T/N_e here). Since T is known, using this approach we can estimate the posterior distribution p(N_e|a₀, a_T, T), which will be proportional to the likelihood for N_e if a uniform prior for N_e is assumed. We need to specify the upper limit, N_eMAX, because convergence of the MCMC simulation is not otherwise guaranteed. Therefore our approach is more accurately described as Bayesian, with an informative rectangular prior on N_e. Here we assume rectangular priors on N_e between zero and some upper limit (N_eMAX—usually 500) and uniform Dirichlet, D (1, … 1), priors on x. However, since the tests of power and precision given in this article are inherently non-Bayesian (since the answer is known), we use the term “likelihood-based” to describe the general approach and restrict the use of the term Bayesian when comparing the effect on inference of different values of N_eMAX. When real data are analyzed, it is clearly sensible to take a fully Bayesian approach, with a prior that reflects background information (which is then unlikely to be rectangular).

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 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 L (N_e, T, x₁, … , x_m) is estimated as described above. The performance of Metropolis-Hastings sampling when the likelihood is estimated with error has not been intensively studied. As in O'Ryan et al. (1998) and Ciofi et al. (1999), we estimated p(a₀|N_e, T, n₀, k, x) for each locus, using 500 importance samples, and reestimated this for current and candidate values at each iteration of the Metropolis-Hastings algorithm. The number of 500, although giving appreciable error in estimation of L (N_e, T, x₁, … , x_m), had been suggested by pilot simulations in O'Ryan et al. (1998) in which it gave indistinguishable results from simulations with 10,000 iterations for subsets of the samples used in that study, and these latter simulations gave indistinguishable results (using smaller data sets) from Metropolis-Hastings sampling using exact likelihoods calculated according to Equation 1. We have not yet determined a lower limit of the number of iterations that will give reasonable performance. It should be noted that X (Equation 5) in the Metropolis-Hastings simulations is the ratio of importance weights and it is straightforward to show that the ratio of estimates using exactly one iteration of the Griffiths and Tavaré procedure is a valid acceptance criterion (in which case we are integrating over G as well). We have investigated this case and find that simulations using one iteration, however, have a very low acceptance rate.

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 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 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 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 σ=k0.5(q2−q1)n0.2, where n is the sample size (the number of updates), q₂ is the value of the 0.625 quantile, and q₁ is the value of the 0.375 quantile. The constant, k, was set to be 1 for all the simulations. This formula was obtained by trial and error in pilot simulations and was set to tend toward under- rather than over-smoothing. A logistic kernel was used because it was easy to define truncated kernels at 0 and N_eMAX, which reduced the degree of underestimation of the density at the boundaries. For appropriate bandwidths (slightly larger values of k) the method gives results very similar to those obtained using the program Locfit (Loader 1996) implemented in R (Ihaka and Gentleman 1996, http://www.r-project.org/), which is based on local-likelihood density estimation methods.

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 likelihood-based estimator .

F-statistic-based estimator: We estimated N_eFk for each population using the equation NeFk^=T2×(F−1∕2S0−1∕2St) (6) (Nei and Tajima 1981), where T is the number of generations between the two samples and S₀ and S_T are sample sizes (of individuals) in the first and second sample, respectively. This equation is appropriate because the individuals sampled to estimate F are independent of the individuals reproducing to make the next generation (sampling scheme II of Waples 1989).

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 Fk=1A−1⋅Σi=1A(xi−yi)2(xi+yi)∕2 (7) (Krimbas and Tsakas 1971; Pollak 1983), where A is the number of alleles at the locus and x_i and y_i are the frequencies of the ith allele in the first and second samples, respectively.

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 χ² 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 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, Sk2 by V(k_i), which we found to be ~1.91 (instead of 2) for N_b = 10. In this case, N_e ~ 10.2. This difference between N_e and N_b is negligible and is even smaller when N_b > 10. Individuals were sampled under sampling scheme II of Nei and Tajima (1981; Pollak 1983; Waples 1989), in which, at given generation, the sample and the 2N_e gametes representing the effective size are both independent binomial samples from the pool of gametes of the preceding generation.

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:

1. 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.

2. 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).

3. 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.

4. 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 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 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 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, MSE ).

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 N_e was below the lower limit and the proportion of times the true N_e was above the upper limit.