1. Introduction
Monte Carlo simulation based on the coalescent is widely used in population genetics. This approach enables researchers to generate data for a given sample size under a panmictic neutral model or other evolutionary models. Using tests of neutrality or other summary statistics, the observed and simulated data can then easily be compared.
In order to simulate a sample under a strict neutral panmictic model, it is necessary to know the population mutation parameter θ (where θ=4Nμ, N being the effective population size and μ the mutation rate). The population mutation parameter, however, is generally unknown and must be estimated (Watterson, Reference Watterson1975; Tajima, Reference Tajima1983; Fu, Reference Fu1994; Griffiths & Tavaré, Reference Griffiths and Tavaré1994; Kuhner et al., Reference Kuhner, Yamato and Felsentein1995). To avoid the uncertainty in using an estimate of θ (usual estimation methods are mostly inefficient, e.g. mean pairwise difference is inconsistent, Watterson's estimate converges only asymptotically), several approaches have been proposed (e.g. see Hudson, Reference Hudson, Takahata and Clark1993; Markovtsova et al., Reference Markovtsova, Marjoram and Tavaré2000; Simonsen et al., Reference Simonsen, Churchill and Aquadro1995). The most popular approximate procedure simulates samples by fixing the number of segregating sites (S) instead of using the mutational parameter θ (Hudson, Reference Hudson, Takahata and Clark1993). However, using this approach, population samples with short genealogical trees tend to exhibit high mutation rates, and conversely, long genealogical trees have low mutational rates.
This latter procedure should not be confused with approaches that condition on S taking into account the parameter θ. A procedure that conditions on S (instead of fixing S) considers for a given θ all possible trees weighted by their probabilities of giving rise to the observed number of segregating sites (Tavaré et al., Reference Tavaré, Balding, Griffiths and Donnelly1997; Markovtsova et al., Reference Markovtsova, Marjoram and Tavaré2000, Reference Markovtsova, Marjoram and Tavaré2001; Jakobsson et al., Reference Jakobsson, Hagenblad, Tavaré, Säll, Halldén, Lind-Halldén and Nordborg2006), and thus takes into account the uncertainty of θ (see Kelly, Reference Kelly1997; Depaulis et al., Reference Depaulis, Mousset and Veuille2001, equation 2). In the case of conditioning on S, the internode times are not independent. Therefore, the procedure proposed by Hudson (Reference Hudson, Takahata and Clark1993) might be considered from two different points of view: (i) as an approximation of the standard procedure that is based on a (known) θ value, and (ii) as an approximation of a rigorous procedure that conditions on S and takes into account the uncertainty of θ.
Although fixing the number of segregating sites does not strictly follow the assumptions of coalescent theory (e.g. independence of the genealogical and mutational phases: see Kingman, Reference Kingman1982a, Reference Kingmanb; Hudson, Reference Hudson, Futuyma and Antonovics1990; Donnelly & Tavaré, Reference Donnelly and Tavaré1995; Nordborg, Reference Nordborg, Balding, Bishop and Cannings2001), this approximate procedure is reasonably accurate for obtaining critical values of statistics for the standard neutral model, provided that the true θ value is well supported by the data (Depaulis et al., Reference Depaulis, Mousset and Veuille2001; Markovtsova et al., Reference Markovtsova, Marjoram and Tavaré2001; Wall & Hudson, Reference Wall and Hudson2001). Thus, the probability of obtaining a θ value that is not supported by the data is expected to be low (Depaulis et al., Reference Depaulis, Mousset and Veuille2001), and the critical values of tests obtained by fixing the number of segregating sites appear to be quite accurate. However, for alternative evolutionary models the accuracy of this approximation is unknown, although many publications have used a fixed number of segregating sites for testing purposes (e.g. Braverman et al., Reference Braverman, Hudson, Kaplan, Langley and Stephan1995; Depaulis & Veuille, Reference Depaulis and Veuille1998; Wall, Reference Wall1999; Fay & Wu, Reference Fay and Wu2000; Przeworski, Reference Przeworski2002; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Glinka et al., Reference Glinka, Ometto, Mousset, Stephan and De Lorenzo2003). It is therefore necessary to study the accuracy of this procedure for a wide range of models.
2. Simulation methods
Simulations have been performed with the program mlcoalsim (Ramos-Onsins & Mitchell-Olds, Reference Ramos-Onsins and Mitchell-Olds2007), which is available at http://www.ub.edu/softevol/mlcoalsim.
(i) Evolutionary models
The following models are considered:
(a) The neutral panmictic model with constant population size.
(b) The symmetric finite-island model (Wright, Reference Wright1943) with a constant number of populations (d), equal population size for each population, and a symmetric and constant migration parameter M among islands (with M=4Nm, where N is the effective population size of a deme and m the migration rate).
(c) The logistic growth model (Fu, Reference Fu1997). In this the population size changes with time as follows:
(1)
Here N(t) is the population size at time t (expressed in N 0 generations), N 0 is the population size at time t 0, and N 1 is the population size at time t 1. We used γ=10/(t 1−t 0). Coalescence times were calculated by integrating (1) over t. Expansion and reduction of population size was studied using this model. We used 10- and 100-fold differences between N 1 and N 0. We also used recent expansion processes, where t 1−t 0=0·1 N 0 generations, and t 0=0.
(d) A bottleneck model using a constant population size initially for a period of time t d (expressed in 4N generations), then a sharp reduction for a time t b followed by a size increase. Reduction and expansion of population size also follow a logistic model. The conditions used here are bottlenecks of 10- or 100-fold differences in effective population size.
(e) Hitchhiking model: We followed essentially the algorithm described in Braverman et al. (Reference Braverman, Hudson, Kaplan, Langley and Stephan1995) to generate hitchhiking genealogies, and allowed recombination within the locus of interest during the selective and neutral phases (see also Fay & Wu, Reference Fay and Wu2000; Kim & Stephan, Reference Kim and Stephan2002; Przeworski, Reference Przeworski2002). As in Fay & Wu (Reference Fay and Wu2000), we used as parameters the selection coefficient (s), the recombination rate between the selected locus and the studied locus (c), the intragenic recombination rate (r), and the time at which an advantageous mutation is fixed (t f). For the selective phase, we calculated the time to the most recent common ancestor (instead of checking at small time increments) for both the selected and unselected ‘subpopulations’ using the reasoning of Nordborg (Reference Nordborg, Balding, Bishop and Cannings2001, equation 7). The selective phase starts at time t f with a frequency of the selected allele of 1−1/2N, and ends when the frequency of x(t)<1/N. The value x(t) was calculated deterministically using equation 1 in Kim & Stephan (Reference Kim and Stephan2002) (see also equation 3a in Stephan et al., Reference Stephan, Wiehe and Lenz1992). The computer code was tested by comparing the results with those of table 2 in Stephan et al. (Reference Stephan, Wiehe and Lenz1992) and also by comparison with the ssw program (Kim & Stephan, Reference Kim and Stephan2002). The parameter values used in this work were t f=0, 4Ns=2×104, N=106 and ε=1/4N.
(ii) Monte Carlo methods
In all procedures, we used sample data (n lines) obtained from a diploid species. The procedures are as follows:
(a) The Fθ procedure (Fixed θ procedure): This is the original standard coalescent procedure. We used a fixed value of the population mutation parameter θ. We placed a number of mutations on the tree according to the Poisson distribution with the mean value θ times the total length of the tree.
(b) The FS procedure (Fixed S procedure): We placed a fixed number of segregating sites uniformly on each tree generated under the models mentioned above (Hudson, Reference Hudson, Takahata and Clark1993). The procedure for generating trees is identical to the Fθ procedure.
(c) Procedure based on fixing the number of segregating sites and the value of θ (FSθ): This procedure employs the rejection algorithm #2 justified and described by Tavaré et al. (Reference Tavaré, Balding, Griffiths and Donnelly1997).
(d) Procedure based on fixing the number of segregating sites but taking into account the uncertainty of the value of θ (FSθprior procedure): We use the RAU procedure (Rejection Algorithm using a Uniform prior). This procedure employs the rejection algorithm #2 described by Tavaré et al. (Reference Tavaré, Balding, Griffiths and Donnelly1997) but sampling of θ is done from a given prior distribution instead of having a fixed θ value.
(iii) Mutational parameter θ
We assume a uniform distribution over some arbitrarily chosen interval [θmin; θmax] (Depaulis et al., Reference Depaulis, Mousset and Veuille2001). Thus the prior density of θ is
When we do not have information about the value of θ, the assumption of a uniform density for any value of θ is reasonable. This is a commonly used strategy to estimate θ given observed data (e.g. Watterson, Reference Watterson1975; Wright & Charlesworth, Reference Wright and Charlesworth2004; Wright et al., Reference Wright, Vroh Bi, Schroeder, Yamasaki, Doebley, McMullen and Gaut2005; Haddrill et al., Reference Haddrill, Thornton, Charlesworth and Andolfatto2005). If the researcher has available information about the distribution of θ, then this information should be used instead of assuming a uniform distribution (e.g. Pritchard et al., Reference Pritchard, Seielstad, Perez-Lezaun and Feldman1999; Przeworski, Reference Przeworski2003), but we do not consider other prior distributions in this paper.
In order to avoid biologically unrealistic values of θ, we used as a minimum bound of θ per nucleotide a value of 0·0005, and as a maximum bound a value of 0·05. Numerous publications show that these numbers are realistic. Different bounds might modify some of our results, but do not change the main conclusions.
(iv) Statistical methods
From the simulated trees we calculated Tajima's D (Tajima, Reference Tajima1989), here named TD, Fu & Li's D and F (Fu & Li, Reference Fu and Li1993), here named FD and FF, respectively, and Fay & Wu's H (Fay & Wu, Reference Fay and Wu2000), abbreviated as H, the statistic B (Wall, Reference Wall1999) (considered sensitive to structured populations), F S (Fu, Reference Fu1997) and R2 (Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002), which are sensitive to population size expansion, as well as the number of haplotypes (here divided by the sample size) K w (Strobeck, Reference Strobeck1987; Fu, Reference Fu1996; Depaulis et al., Reference Depaulis, Mousset and Veuille2001; Wall, Reference Wall1999) and the haplotype diversity H w (Depaulis & Veuille, Reference Depaulis and Veuille1998). In total, nine summary statistics were computed. The calculation of these statistics was examined with the software package DnaSP 3.51 (Rozas & Rozas, Reference Rozas and Rozas1999; Rozas et al., Reference Rozas, Sanchez-DelBarrio, Messeguer and Rozas2003).
For multilocus analyses, we treated each locus independently because we assumed that the studied loci are unlinked. When we used the RAU procedure, we chose independent θ values for each locus in order to perform simulations. For a given statistic, we recorded the P value independently for each locus and combined them as in Voight et al. (Reference Voight, Adams, Frisse, Qian, Hudson and Di Rienzo2005) but calculated each tail separately.
To avoid excessively liberal critical values (95% interval of the null sampling distribution) for discrete distributions, we have used the following procedure: for the 2·5% interval of the upper tail, we took the first value that is larger than the observed value at 2·5% of the distribution (see Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). The same logic was applied to the 2·5% interval of the lower tail and also for comparing discrete distributions. Because of this issue, the realized type I error rates are slightly lower than the expected 5% and have to be assessed for every statistic and method.
(v) Determining the accuracy of the FS procedure in relation to Fθ
We have used the approach described in Wall & Hudson (Reference Wall and Hudson2001) to determine the level of accuracy of the FS procedure compared with the Fθ procedure. In this approach, a number of simulations using a large range of values of S (e.g. from S=1 to S=120) is obtained by the FS procedure, and its critical values and the false positive rate (called the ‘nominal’ size) for each statistic are calculated for a given evolutionary model. Then, the ‘true’ θ value for the same evolutionary model is obtained, which is the value that gives an average of S=20 in 10 000 iterations under the Fθ procedure. Next, a simulation with the Fθ procedure (i.e. the null hypothesis) is performed using the ‘true’ fixed value of θ. For each iteration, the value of S is calculated, and acceptance or rejection of the null hypothesis is determined for every statistic by comparing its value with the critical values of the corresponding FS distribution. The rejection rate (i.e. the false positive rate given the FS critical values, here called the ‘true’ size) for each statistic is stored. The choice of the ‘true’ θ value, although somewhat arbitrary, is suggested by the design of the study, which fixes the number of segregating sites instead of fixing the value of θ (which is unknown to the researcher). Using the following method, we try to be conservative in the sense of minimizing the differences between the two different procedures (see Wall & Hudson, Reference Wall and Hudson2001; Depaulis et al., Reference Depaulis, Mousset and Veuille2003, Reference Depaulis, Mousset, Veuille and Nurminsky2005).
(vi) Determining the accuracy of the FS procedure in relation to FSθU
Here we assume that the θ value is unknown. The FSθU procedure is considered the ‘true’ procedure because it takes into account the uncertainty assumed by the researcher. In this approach, a simulation for a given evolutionary model is performed for the FS procedure (for a fixed S=20) and its critical values and the false positive rate (called the ‘nominal size’) of each statistic are stored. Then, a simulation using the FSθU procedure (for S=20) is performed for the same evolutionary model (the null hypothesis). The ‘true’ false positive rate is calculated in the following way: for each iteration calculated with the FSθU procedure, rejection of the null hypothesis is determined for every statistic by comparing its value with the critical values of the FS distribution. Finally, the rejection rate for each statistic is stored.
(vii) Parameter values
We used n=20 for a locus of 1000 nucleotides. For procedures that fix S, we preferentially used S=20. The reason for choosing n=20 and S=20 is largely historical, as the value of θ (given the neutral equilibrium model) is then approximately 0·005 (which has been used in a number of the theoretical population genetics studies). All θ values under any evolutionary model studied here were in the range we considered biologically realistic. For n=20 and S=50 we found similar results.
3. Analytical approaches
(i) Effects of recombination
Recombination is difficult to take into account in analytical models. However, one can consider the extreme case of free recombination (e.g. a sequence consisting of m freely recombining fragments of equal size such as m nucleotides). This is the limiting case when the recombination parameter R=4Nr becomes very large. Let the parameter T be the vector of the coalescent times T k, then 
is the total length of the coalescent tree (by summing the length of all branches) measured in units of 4N generations. Assuming m is large, the central limit theorem shows that the prior distribution of the average length 
of the m independent coalescent trees will converge towards the normal distribution with mean a 1(n) and variance a 2(n)/m (defined in equations (A3) and (A4) of Appendix 1, respectively). This is,
where 
is the density of the average length 
. This last equation can be used in equations (A7)–(A8) (Appendix 1) to address the effects of the simulation procedures Fθ, FS, FSθ, and FSθprior on the posterior density of the length of the coalescent tree in the limiting case of freely recombining sequences.
(ii) Shape of the coalescent trees
In the present work, we focus on procedures where simulations use a fixed number of segregating sites S. In these procedures, unlike in the Fθ procedure, the total length of a tree does not play a role as long as the shape of the tree remains the same (trees are only scaled). In order to assess the impact of the rejection algorithm on the shape of the trees, we study the ratio of the branch lengths in the upper and lower parts of the tree X/Y n (Fig. 1). The mean value of this ratio is given by
where f X and 
denote the prior densities of X and Y n (see equations (A5) and (A6)) and
and 
is a normalizing constant.
Fig. 1. Shape index of a coalescent tree. The shape index is the ratio of the branch lengths in the upper and lower parts of the tree, X/Y n (see text).
4. Simulation results
We determined the accuracy of the approximate procedure proposed by Hudson (Reference Hudson, Takahata and Clark1993) (which fixes the number of segregating sites while ignoring the value of the population mutation value θ; the FS procedure) in comparison with the two rigorous procedures: (i) the standard procedure fixing θ (Fθ) and (ii) procedures that while conditioning on S take into account the uncertainty of θ (FSθprior, here using a uniform prior, named FSθU). We consider several different evolutionary models including the neutral panmictic model, population subdivision, population size bottleneck, expansion and genetic hitchhiking. As a measure of accuracy, we calculated the difference in the false positive rate between procedures.
(i) Comparison of the FS procedure with the standard procedure Fθ
We examined the accuracy of the FS procedure for different alternative models. To do this we also re-examined the type I error for the standard neutral model (neutral panmictic population with zero recombination) for the nine statistics studied here, since the type I error for single-locus data was also studied in Depaulis et al. (Reference Depaulis, Mousset and Veuille2001), Markovtsova et al. (Reference Markovtsova, Marjoram and Tavaré2001) and Wall & Hudson (Reference Wall and Hudson2001).
Table 1 shows the difference between the nominal size (the false positive rate of a test given the FS procedure) and the ‘true’ size (the false positive rate of a test given the Fθ procedure) for a 5% critical region. Our analysis confirmed the small difference in the type I error rate observed by Wall & Hudson (Reference Wall and Hudson2001) and Depaulis et al. (Reference Depaulis, Mousset and Veuille2001) under a neutral model for a single locus, although we used some other statistics and slightly different conditions (n=20 and θ=0·0057 for the Fθ procedure; i.e. the estimate of θ for S=20). Fig. 2 A shows the differences in size for each S value among the procedures for the nine studied statistics under the neutral model. Negative values indicate that the critical values of the FS procedure are more liberal than expected (increased type I error). The large fluctuations observed in Fig. 2 are also a consequence of the discrete distribution of values obtained when the number of segregating sites is fixed (see Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). Important differences are observed in the proportion of acceptance/rejection when the observed number of mutations is far from the average (S=20), as indicated in Wall & Hudson (Reference Wall and Hudson2001) (see also Depaulis et al., Reference Depaulis, Mousset and Veuille2003, Reference Depaulis, Mousset, Veuille and Nurminsky2005). Nonetheless, the more extreme values contribute very little such that when summed the differences cancel each other out, thus leading to a good accuracy of the FS procedure for the standard neutral model (Table 1). Therefore, we consider the FS procedure as sufficiently accurate under the neutral panmictic model for single-locus data.
Fig. 2. Differences between the sizes of FS and Fθ procedures for each S separately. (A) Neutral model: In the upper panel the distribution of S values obtained with the Fθ procedure using θ=0·0057 is shown. In the middle and lower panels the size differences for statistics for the upper (97·5%) and the lower (2·5%) tails, respectively, are presented. Note that large S values, although causing large differences, contribute very little to the total. (B) Hitchhiking model: In the upper panel the distribution of S values obtained with the Fθ procedure using θ=0·0127 is shown. In the middle and lower panels the size differences for statistics for the upper (90%) and the lower (10%) tails, respectively, are presented. Abbreviations are explained in Section 2. Critical values are not calculated for B (lower tail), and for K w, and H w (upper tail), because these statistics are not conservative with recombination.
Table 1. Difference (in percentage) between the FS and Fθ procedures for several critical valuesFootnote a
a Abbreviations are explained in Simulation methods. The critical values studied are 2·5% and 97·5% for the neutral model and 10%, 50% and 90% for alternative models (see Simulation methods). The maximum precision error detected is around ±1%.
b Island model with d=10 subpopulations, 4Nm=0·5 and 4Nr=10. The samples are all obtained from a single population. In the FSθU procedure, the bounds of θ were reduced d times for a better comparison with the other models. In the Fθ procedure, θ=0·00062.
c Logistic expansion model with a 10-fold growth of population size. The expansion process started 0·4N 0 generations ago and finished at present. 4N 0r=10. In the Fθ procedure, θ=0·03384.
d Logistic contraction model: In contrast to the logistic expansion model, the population is reduced 10 fold. In the Fθ procedure, θ=0·00059.
e Bottleneck model: Population size is maintained constant during 0·5N 0 generations since present. It follows a logistic, 100-fold reduction of population size for 0·01N 0 generations, maintained for 0·01N 0 generations, then an instantaneous recovery to the present population size. 4N 0r=10. In the Fθ procedure, θ=0·00860.
f Hitchhiking model: The selective phase is completed at present time (tf=0). N=106, 4Ns=104, the population recombination parameter between the selected and observed locus is 4Nc=400, c/s=0·02, and intragenic recombination is indicated in the table. In the Fθ procedure, θ=0·00127.
Next we consider the accuracy for alternative models. Table 1 shows the difference between nominal (FS procedure) and true (Fθ procedure) distributions of statistical tests for certain critical values. The reason to use different critical values for alternative models (10%, 50% and 90% instead 2·5% and 97·5%) is because we are interested in whether the probability distribution for a given statistic can be accurately represented using the FS procedure. The analysis is performed for five different alternative models (subdivision, expansion, contraction, bottleneck and hitchhiking) and for some arbitrary parameters. The strong differences observed between the two procedures indicate that the FS procedure is inaccurate when alternative models are used. Fig. 2 B shows considerable differences between nominal and true sizes for several of the studied statistics for each S separately, leading to a large difference in total (Table 1).
(ii) Comparison of the FS procedure with a procedure that conditions on S and takes into account the uncertainty of θ
We have examined the accuracy of FS (the nominal procedure) by comparison with the procedure FSθU (the ‘true’ procedure, see Section 2). With regard to the type I error for the standard neutral model, we have studied the 95% interval of the nine neutrality test statistics for the FS procedure and the procedure FSθU by simulation for n=20 and a range of S values from S=2 to 80. Fig. 3 shows the difference between the nominal and ‘true’ sizes. Liberal tests are Fay & Wu's H with regard to the upper tail, and Fu's F S, Ramos-Onsins & Rozas' R2, and Fu & Li's D and F for the lower tail. In the case of Fu & Li's F, the test is, in fact, less liberal than shown, because the critical value for the FS procedure is sometimes much lower than 2·5% (given the discrete distribution of values; not shown). We observed that all statistics at intermediate S values show a difference lower than 1·5%. Therefore, in this comparison the FS procedure is sufficiently accurate for statistical inferences under the neutral panmictic model for single-locus data.
Fig. 3. Differences between the sizes of the FS and FSθU procedures. Critical values are obtained for each of the nine neutrality tests for the 2·5% upper and lower tails. Abbreviations are indicated in Section 2.
In contrast, for alternative models there are strong differences between the nominal (FS procedure) and true (FSθU procedure) distributions of statistical tests for 10%, 50% or 90% critical values (Table 2). The important differences observed between the two procedures indicate that the FS procedure is also inaccurate when it is compared with the FSθU procedure and when alternative models are used.
Table 2. Difference (in percentage) between the FS and the FSθU procedures for several critical valuesFootnote a
a As in Table 1.
We performed a multilocus analysis for the nine test statistics using the combined P values for different numbers of loci (see Section 2) for the standard neutral model. The difference between the nominal and true size is examined in Fig. 4. The results indicate that for a larger number of loci the FS procedure becomes inaccurate (through the accumulation of small departures in the single-locus tests). Thus, the small differences observed for one locus between FS and FSθU become very important when multilocus studies are performed. Some statistics, like K w, H w and Fay & Wu's H, are extremely conservative for the lower tail, as are most statistics (except for Fay & Wu's H) for the upper tail (see Fig. 4 A). On the other hand, statistics such as Tajima's D, Fu's F S, Fu & Li's D and F, and Ramos-Onsins & Rozas' R2 are too liberal for the lower tail of the distribution, and Fay & Wu's H is also too liberal for the upper tail of the distribution.
Fig. 4. Effect of the number of loci on the probability of rejecting the neutral panmictic model for nine neutrality tests. (A) Differences between the sizes of FS and FSθU procedures with no recombination in the upper 2·5% tail, given different numbers of loci. (B) Differences between the sizes of FS and FSθU procedures with no recombination in the lower 2·5% tail, given different numbers of loci. Abbreviations are explained in Section 2. Critical values are not calculated for B (lower tail), and for K w, and H w (upper tail), because these statistics are not conservative with recombination. S was fixed at 20 for each locus. Plots obtained using S values from a distribution compatible with θ=0·0057 gave equivalent results (not shown) although for a small number of loci we observed a large variance.
5. Discussion
We have compared the FS procedure with the standard procedure using a fixed θ value (Fθ) and with a procedure that conditions on S taking into account the uncertainty of θ (FSθU). Our results show that the FS procedure is inaccurate in both comparisons when alternative models are used.
(i) Causes of the discrepancy between the FS and the Fθ or FSθU procedures
In the comparison of FS with Fθ, differences in the length (L n) and topology (X/Y n) of the trees for the FS and FSθ procedures can explain the inaccuracy of FS, given that the critical values of FS are compared with the FSθ distribution for each S. The distribution of the total length of the coalescent trees of a non-recombining sequence is shown in Fig. 5. The shapes of the FS distribution and the distribution obtained with the FSθ procedure (taking S=20) are clearly different. The distribution of L n obtained with the FSθ procedure has lower variance than the FS procedure, which can be explained by the lack of uncertainty in the value of θ (see also Table 3). The mean and standard deviation of the shape index X/Y n for the different simulation procedures are shown in Table 4. Higher values are obtained for trees with long internal branches whereas smaller ratios denote shorter internal branches. We observed that for S=20 the trees obtained with the FSθ procedure have smaller ratios X/Y n than the trees obtained with the FS procedure. These observations can be explained as follows. For a given θ, when S is lower than the average for that θ, this will often have been caused by a tree that is shorter than average. Although the internode times are a priori independent (because the total tree length is dominated by the last few internode times, given a short tree), it will often be that the final internode time is particularly short. Thus X/Y n is smaller than the unconditional expectation. The situation is reversed for S greater than the average for a given θ. Thus, these results show how the simulation procedure may affect the shape of the sampled coalescent trees and may lead to different outcomes of neutrality tests.
Fig. 5. Posterior densities of the total length L n of a coalescent tree with n=20 are shown for the different procedures and no recombination. Results based on the FS, FSθ and FSθU procedures are displayed. For FSθ the θ parameter was arbitrarily set to 20/a 2(20).
Table 3. Mean and standard deviation (in parentheses) of the posterior distribution of L 20
a Assuming m=50 recombining fragments (equation 3).
b Assuming θ=20/a1(20).
Table 4. Mean and standard deviation (in parentheses) of the tree shape
a See Figure 1. Assuming no recombination and n=20.
b Assuming θ=20/a1(20).
To understand the causes of the inaccuracy of the FS procedure in comparison with the procedure FSθU, we examined differences in the trees by considering the total length of the tree (L n) analytically and by simulations. The distribution of the total length of the coalescent tree of a non-recombining sequence is shown in Fig. 5. Some inaccuracies of the FS procedure in comparison with the procedure FSθU are observed under a panmictic neutral model with no recombination. Differences in the distributions of L n are small but apparent, with longer and slightly broader distributions for the procedure FSθU than for FS. It is noteworthy that this procedure samples shorter trees with a lower variance of L n than the FS procedure (Table 3). When we analysed the shape of the coalescent trees (Table 4), the coalescent trees obtained by the FSθU procedure are skewed towards shorter trees with smaller internal branches than trees sampled with the FS procedure.
Figure 6 shows the L n distribution for four of the alternative models analysed with the FS and FSθU procedures. There are strong differences in the distribution of L n for the FS procedure in comparison with FSθU for most of the alternative models. In particular, in the case of subdivided populations, large differences are observed for the parameter values studied. Bottleneck and hitchhiking models show also apparent differences in the trees that are used in FS in relation to FSθU. On the other hand, models of population size expansion exhibit smaller differences among procedures. In summary, differences in the length (and in the topology) of the trees explain the inaccuracy of the FS procedure and these differences are large in most of the alternative models studied.
Fig. 6. Probability distribution of L n for four alternative models using n=20 and S=20. The FS and FSθU procedures are displayed in each case. Parameter values are the same as in Table 1 but intragenic recombination was set to 4Nr=10 (except for panel D where it was zero). (A) Island model. (B) Logistic expansion model. (C) Bottleneck model. (D) Hitchhiking model.
(ii) The effect of recombination
Analysing the effect of recombination is also important, as the true value of the recombination parameter is generally unknown. Recombination has an important effect on the distribution of segregating sites because it breaks up the correlation among contiguous positions. Differences in L n are observed when recombination is added to the model (Fig. 7). Average values of the parameter L n are constant under the FS procedure, as expected. Indeed, the FS procedure uses all output trees, like the standard coalescent procedure uses a fixed θ value. The expected average value of L n for the FS procedure value is given by equation (A3) in Appendix 1. Therefore, for n=20, E(L n)=3·548. On the other hand, the procedure FSθU has lower average L n values for zero recombination (L n is 3·16 for n=20), and this value increases for larger recombination values. The procedure FSθ leads to a similar pattern as that shown for FSθU, but with fewer differences relative to FS.
Fig. 7. Effect of recombination on average L n for the FS, FSθ (θ=0·0057) and FSθU procedures in a neutral panmictic population using n=20 and S=20.
The means and standard deviations of L n in the limiting case of free recombination between fragments are given in Table 3. Recombination tends to bring the posterior distributions closer to the prior distribution of L n. Whereas the difference between the FS and FSθU simulation procedures remains clear, the outcomes of the FS and FSθU simulation procedures can hardly be distinguished (independent of the assumed prior for θ), confirming the observation of Fig. 7.
The consequences of having differences in the average L n given no recombination might be important, because it indicates that the zero-recombination neutral panmictic model may have an average of a given statistic that is different from that of a model with recombination, and therefore may not be conservative. Zero recombination in a neutral panmictic model leads to the largest deviation in average L n in relation to the FS procedure. We have observed that for single loci this difference can be tolerated, but becomes too large for multilocus analyses. That is, for multilocus analyses, the FS procedure should not be used unless recombination is quite high for each locus.
The prior distribution of the mutational parameter θ is shown to have an important effect on the posterior distribution of coalescent trees (Fig. 5), although, as previously pointed out by Wall & Hudson (Reference Wall and Hudson2001), the recombination parameter also has a strong effect on this distribution. The improvement obtained from using sophisticated techniques such as the FSθU procedure may, however, be negligible in comparison with the errors caused by the uncertainty in the recombination parameter R (Wall & Hudson, Reference Wall and Hudson2001). A procedure considering also the uncertainty in the recombination parameter would give more appropriate distributions, but is beyond the scope of this paper.
6. Appendix 1
(i) The branch lengths in a coalescent tree
In the standard neutral coalescent process without recombination, the waiting time T n (in units of 4N generations) until two lineages among n have a common ancestor is exponentially distributed with parameter n (n−1) (Kingman, Reference Kingman1982b). The probability density of T n is thus
The density of the sum L n of the lengths of the branches of a coalescent tree with n tips is given by Tavaré (Reference Tavaré1984, unlabelled equation at the top of p. 153) as
The mean and variance of L n are
In order to characterize the shape of a coalescent tree, we introduce a shape index X/Y n where X is the length of the two branches from the second-last coalescent to the last (in the upper part of the coalescent tree) and Y n is the length of the rest (the branches in the lower part) of the coalescent tree (see Fig. 1). The density of X follows from equation (A1):
As shown in Appendix 2, the density of Y n is
(ii) Branch lengths in the procedures Fθ and FS
In the standard procedure Fθ, but also in the FS simulation procedures, all simulated coalescent trees are used and thus the posterior density of tree length is the same as the prior density 
(equation A2).
(iii) Branch lengths in the FSθ and FSθprior procedures
In the FSθ simulation procedure, trees are sampled from the prior distribution according to the probability of observing exactly S mutations given a known mutational parameter θ. The number of mutations in the coalescent tree follows a Poisson distribution with parameter θL n. Thus the probability of observing S mutations in a sample of n lines given θ is (Tavaré, Reference Tavaré1984, unlabelled equation on p. 153)
and the posterior density of L n given S and θ follows from the definition of the posterior density:
In the FSθ procedures that consider the uncertainty of the value of θ (FSθprior), a prior distribution of θ with density g is assumed. Thus the probability of observing S mutations in a sample of n lines is
and the posterior density of L n given S and the prior density g is
If a uniform density g u is assumed for θ, g u is constant (equation 2) and the posterior density in equation (A8) becomes
Using θmin=0, we see that for n⩾3 a limit of the posterior distribution of L n exists when θmax is very large (this condition is necessary for the integral in the denominator of equation A9 to converge):
It is noteworthy that this distribution is independent of the observed number of mutations S and may be denoted as 
.
7. Appendix 2. Proof of equation (A6)
(i) Notation and preliminary results
In this section the density function of a random variable X will be denoted as f X. We introduce a new random variable P k=kT k (see equation A1). It is straightforward to show that P k is exponentially distributed with parameter k−1, thus
and Y n is defined by 
(see Fig. 1). Y n is the sum of n−2 independent exponentially distributed random variables with parameters 2, …, (n−1). As we show below, the density function of Y n can be obtained using simple order statistics.
The following properties of the exponential distribution will be used:
• Minimum of independent exponentially distributed random variables: Consider k independent random variables
exponentially distributed with parameters λi (i=1, …, k). Then is exponentially distributed with parameter .• Memoryless property of the exponential distribution: Consider a random variable X λ exponentially distributed with parameter λ:
exponentially distributed with parameters λi (i=1, …, k). Then 
is exponentially distributed with parameter 
.
(ii) Density function of Y n
We consider n−1 independent random variables E i(i=1, …, n−1), exponentially distributed with parameter 1, and denote the smallest one as E (1)⋅ E (1) is exponentially distributed with parameter n−1, and the n−2 remaining random variables are independent and such that E i−E (1) is exponentially distributed with parameter 1. Then the difference between the second and the first smallest random variables E (2)−E (1) is exponentially distributed with parameter n−2 and the n−3 remaining random variables are independent and such that E i−E (2) is exponentially distributed with parameter 1. If we define E (0)=0, it is straightforward to show that the difference between two successive sorted random variables (E k−E (k−1)) is exponentially distributed with parameter n−k. E (n−2) can be rewritten as 
. This is the sum of n−2 independent exponentially distributed random variables with parameters 2, …, (n−1), thus E (n−2) and Y n have the same distribution. Because E (n−2) is the (n−2)th random variable among the n−1 random variables E i sorted in increasing order, the density function of E (n−2) and Y n is (Pitman, Reference Pitman1992, p. 326):
We thank J. M. Comeron, J. Rozas, K. Schmid, J. Zavala and A. Navarro-Quezada for helpful comments on the manuscript. S. R. O. thanks S. Gebauer-Jung and S. Wright for help with the calculation of Hudson's recursive equation, J. Fay for supplying information to check the selective model and Y. Kim for his assistance in the analysis of the selective model. Special thanks go to two anonymous reviewers for their input into this manuscript. This research was supported by the Max Planck Society, a Marie Curie fellowship from the EU (MCFI 2002 01461) to S. M., and grant Ste 325/5 from the DFG to W. S.
