1. Introduction
The motivation for studying the impact of past demography on sequence data is two-fold. Firstly, changes in population size are interesting in their own right, being intimately linked to processes such as speciation or geographic range shifts. Secondly, the standard neutral model (SNM) of a randomly mating Wright–Fisher population of constant size and discrete generations, hardly ever describes the patterns of diversity found in natural populations. Thus, studies aiming to detect loci under selection are faced with the considerable challenge of fitting realistic demographic models against which selection can be tested e.g. Glinka et al. (Reference Glinka, Ometto, Mousset, Stephan and De Lorenzo2003), Hamblin et al. (Reference Hamblin, Mitchell, White, Gallego, Kukatla, Wing, Paterson and Kresovich2004), Haddrill et al. (Reference Haddrill, Thornton, Charlesworth and Andolfatto2005), Ometto et al. (Reference Ometto, Glinka, De Lorenzo and Stephan2005) and Thornton & Andolfatto (Reference Thornton and Andolfatto2006). Since the rate of coalescence is inversely proportional to the effective population size, it is clear that demographic changes must leave a detectable signature in genealogies (Felsenstein, Reference Felsenstein1992). In general, positive population growth distorts genealogies towards a starshape with shorter internal branches, resulting in more low frequency variants and a unimodal rather than multi-peaked mismatch distribution (Slatkin & Hudson, Reference Slatkin and Hudson1991; Harpending, Reference Harpending1994; Schneider & Excoffier, Reference Schneider and Excoffier1999). In contrast to selective processes that act on single genetic variants, demography affects the whole genome, so one expects to find a concordant signature across loci (Tajima, Reference Tajima1989; Galtier et al., Reference Galtier, Depaulis and Barton2000).
Approaches to demographic inference fall into three broad categories; for a review see Emerson et al. (Reference Emerson, Paradis and Thebaud2001). Firstly, likelihood methods, which are available for bottleneck and exponential growth models, make use of all the information in a sample by integrating over a large set of likely genealogies (Griffiths & Tavaré, Reference Griffiths and Tavaré1994; Kuhner et al., Reference Kuhner, Yamato and Felsenstein1995). Although optimal in terms of statistical power and accuracy, likelihood estimation is computationally intensive and requires a fully specified alternative model. Therefore realistic growth histories often remain analytically intractable. Secondly, there are tree-based methods, which take the branch length information of a reconstructed tree as their starting point. Assuming that sequence evolution is clock-like, the number of lineages can be plotted against time and the shape of this trajectory compared with its neutral expectation (Nee et al., Reference Nee, Holmes, Rambaut and Harvey1995; Pybus et al., Reference Pybus, Rambaut, Holmes and Harvey2002). Despite their conceptual appeal, these methods neglect any uncertainty in tree topology and are thus only as good as the reconstructed tree they are based on. Furthermore they cannot deal with recombination by definition. Finally, there are classical neutrality tests, most of which do not explicitly consider the genealogy but instead use more immediate aspects of the data such as the frequency spectrum of mutations, e.g. Tajima's D (Tajima, Reference Tajima1989) and Fu & Li's D (hereafter referred to as D 2) (Fu & Li, Reference Fu and Li1993), the haplotype distribution, e.g. Fu's F S (Fu, Reference Fu1996; Innan et al., Reference Innan, Zhang, Marjoram, Tavaré and Rosenberg2005), or the mismatch distribution, e.g. the raggedness statistic (Slatkin & Hudson, Reference Slatkin and Hudson1991). Compared with likelihood estimation, summary statistics are straightforward to calculate and their distribution can be simulated under almost any growth model.
Considering the zoo of statistics available and their wide use, there are surprisingly few studies that systematically compare their power, and those that do mainly consider bottlenecks and single locus data (Simonsen et al., Reference Simonsen, Churchill and Aquadro1995; Fu, Reference Fu1996; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Depaulis et al., Reference Depaulis, Mousset and Veuille2003; Ramirez-Soriano et al., Reference Ramirez-Soriano, Ramos-Onsins, Rozas, Calafell and Navarro2008). However, joint analysis of multiple loci is not only necessary to distinguish between selective and demographic events (Galtier et al., Reference Galtier, Depaulis and Barton2000) but also potentially far more powerful than inferences based on a single locus. An added advantage of multi-locus analysis is that both means and variances of summary statistics can be used for testing. Variance based tests were first developed for microsatellite data (Di Rienzo et al., Reference Di Rienzo, Donnelly, Toomajian, Sisk, Hill, Petzl-Erler, Haines and Barch1998; Reich et al., Reference Reich, Feldman and Goldstein1999) but are now routinely used to analyse sequence data from multiple loci (Pluzhnikov et al., Reference Pluzhnikov, Di Rienzo and Hudson2002; Haddrill et al., Reference Haddrill, Thornton, Charlesworth and Andolfatto2005; Heuertz et al., Reference Heuertz, De Paoli, Kallman, Larsson, Jurman, Morgante, Lascoux and Gyllenstrand2006) or even species (Hickerson et al., Reference Hickerson, Dolman and Moritz2006).
A general conclusion that has emerged from simulation studies is that tests based on the number and distribution of haplotypes have more power to detect bottlenecks than statistics based on π, in particular Tajima's D (Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Innan et al., Reference Innan, Zhang, Marjoram, Tavaré and Rosenberg2005; Ramirez-Soriano et al., Reference Ramirez-Soriano, Ramos-Onsins, Rozas, Calafell and Navarro2008). Earlier, Felsenstein made a theoretical argument for the inferiority of pairwise measures (Felsenstein, Reference Felsenstein1992). Their large variance under neutrality arises both from their sensitivity to the last coalescence event and the random genealogical topology (Tajima, Reference Tajima1983). Under the SNM more symmetric genealogies are on average associated with higher π and more ragged mismatch distributions than asymmetric genealogies. It is important to realize that this topological variance is independent of the already large variance in coalescence times inherent in the genealogical process. In other words ‘despite their aura of robustness’ (Felsenstein, Reference Felsenstein1992), statistics based on π suffer from an unnecessarily large variance under neutrality, and hence have comparatively low power. Despite these results, D and mismatch distributions continue to be the methods of choice for demographic inferences in population genetics and phylogeography, respectively.
Following Felsenstein's recommendation that ‘there is much to gain from explicitly taking the genealogical relationship of a sample into account’ (Felsenstein, Reference Felsenstein1992), the aim of this study is to consider how genealogical information can be used for demographic inference in a summary statistics framework. Our premise here is that the mutation rate is sufficiently high relative to the per site recombination rate such that non-recombining blocks of sequences can be easily identified and treated as independent loci.
Given that there is usually not enough information in within-species sequences data to infer the full topology unambiguously it seems important to ask which part of the topology yields most information. The first part of the paper introduces some simple measures of starshape, which are based on the properties of a rooted genealogy. Using simulations, their power to detect a history of exponential growth is compared with standard neutrality tests for both the single and multi-locus cases. We focus on the exponential growth model for two reasons. Firstly, although it is a frequently used demographic model, the power of summary statistics to detect exponential growth has been little investigated. Secondly, likelihood methods are available, which can be taken as an absolute ‘upper bound’ of power for comparison. Such a direct comparison between summary statistics and the optimal likelihood methods is lacking so far.
2. Summary statistics
Several neutrality tests compare two different estimators of the scaled mutation rate (Fu & Li, Reference Fu and Li1993; Tajima, Reference Tajima1989; Fay & Wu, Reference Fay and Wu2000) θ=4N eμ, where μ is the mutation rate and N e the effective population size, which capture different aspects of the data. Most prominently, Tajima's D is defined as the difference between θ estimated as π, and θw=S/a n (Watterson's θ, where 
, n is the sample size and S the total number of polymorphic sites in the sample), normalized by the standard deviation of this difference. Genealogies from growing populations typically have relatively more low frequency variants and hence tend to have a negative D.
While neutrality tests are commonly based on the frequency spectrum and π, it is instructive to consider departures from the SNM in terms of their effect on the genealogy. Such tree-thinking necessarily underlies summaries that make use of outgroup information, e.g. D 2 has a straightforward genealogical interpretation. Below two different ways of employing genealogical information in the construction of summary statistics are considered.
(i) Genealogical ratios
The rationale behind D 2 is to distinguish between two classes of mutations: those found on terminal branches, ηe and those on internal branches, ηi (Fig. 1) (Fu & Li, Reference Fu and Li1993). Suppose that some limited topological information can be inferred from the data. In particular, we will for now assume that the placement of the root is known. It is then possible to distinguish mutations found on the two rootward branches, which we shall denote ηR. Under the infinite sites assumption, these are all derived mutations that are shared by all individuals in either of the two sub-clades defined by the root. The advantage of considering the proximity of mutations to the root rather than the tips is twofold: firstly, rootward branches cover a greater proportion of the time to the most recent common ancestor of the sample (T MRCA) and should, in general, be more informative about past changes in population size. Under the SNM, on average half of the T MRCA is taken up by the coalescence of the last two lineages (T 2) (Fig. 1), whereas in a growing population, the smaller population size in the past forces the last two lineages to coalescence much more rapidly. Secondly, the average length of a branch connected to the root is less dependent on the sample size than the average length of a terminal branch.
Random genealogy of a sample of 20 sequences. The root partitions the sample into two subclades of size 3 and 7. Rootward branches are shown as bold, terminal branches as dotted lines, mutations are represented as crosses. The time interval until the last coalescence event, T 2, is shorter than average under the SNM. In this example S=30, ηR=7, ηRmin=2 and ηe=14.
Ideally, one wants to know the total number of mutations that have occurred during T 2, rather than the number of mutations on both rootward branches, ηR which is larger and depends on the topology, i.e. the order of the first node on the longer of the two branches (Uyenoyama, Reference Uyenoyama1997, Appendix).
One possibility is to only consider the shorter of the two rootward branches that has exactly length T 2. Thus the number of mutations found on this branch, ηRmin, over θw constitutes a very simple measure of starshape.
Such genealogical ratios have first been employed to study the effect of balancing selection on plant incompatibility loci (Uyenoyama, Reference Uyenoyama1997). Being based on a single random event, X clearly neglects much of the information contained in the genealogy. Its power is limited by the probability of observing ηRmin=0 under neutrality. In other words, X is unlikely to be of much use in the case of a single locus.
Alternatively, one can ignore the uncertainty in node order and take the number of mutations found on both rootward branches relative to θw:
It is possible of course to construct various composite measures from the number of mutations found on different parts of the genealogy. Here, we only consider one additional statistic, the relative difference between rootward and terminal mutations:
The X statistics assume some knowledge of the tree topology that is usually unknown. Of course one could use some standard method of tree reconstruction and infer ηR and ηRmin from the most likely topology. However, not only is it inefficient to reconstruct the full topology when all that is required is the placement of the root, conditioning on a single tree also ignores any topological uncertainty. We have therefore developed a simple scheme of inferring the root in a sample of polarized sequences that circumvents these problems.
Under the infinite sites assumption, a necessary criterion for the root-node is that no mutations are shared between the two subsets on either side. One can show that if both branches connected to the root carry mutations, i.e. ηRmin>0 there exists exactly one bipartition of the sample with no mutational overlap. If however one or both of the rootward branches of the genealogy carry no mutations there may be multiple bipartitions that meet this criterion. In this case ηRmin=0 and the tree reconstructed from such a sample would have an unresolved polytomy at its base. To incorporate the topological uncertainty about the placement of the root, we compute the average value of ηR over all partitions that are compatible with the criterion of no mutational overlap. Note that in contrast to most tree reconstruction algorithms that join similar sequences (i.e. start from the tips down the tree), our scheme is divisive (i.e. it starts from the root). To avoid having to consider all possible bipartitions of the sample (2n−1−1), we make use of the fact that any sequences that share mutations have to be on the same side of the root. By first binning sequences that share at least one mutation, we can directly calculate ηR and the number of possible partitions.
(ii) Starting from the limiting case
A different approach is to construct summaries that measure departures from the limiting case of a perfectly star-shaped genealogy. Star-shaped genealogies have some convenient properties that can be used for this. Assuming that outgroup information is available, one can record the number of terminal mutations in each sequence i (because lineages are exchangeable, the labelling is arbitrary), V i. In a perfectly star-shaped genealogy, all mutations must fall onto terminal branches by definition. Thus one expects the number of derived mutations in a sequence to be half the average pairwise diversity, i.e. E[V i]=π/2. The statistic R 2E proposed by Ramos-Onsins and Rozas measures the average departure from this expectation:
(Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002, eqn (2)). R 2E has proven superior to a wide range of summary statistics in detecting histories of bottlenecks (Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). However, because of its dependence on π, one may suspect it to suffer from a large variance under neutrality. We therefore consider a similar statistic that uses the observed S rather than π to assess the degree of starshape. Consider the total number of derived mutations in each sequence, D i. Note that 
, in terms of the unfolded frequency spectrum, where ξi denotes derived mutations that occur i times in the sample. Using the fact that E[D i]=S/n in a star-shaped genealogy we can define a new statistic:
Since under neutrality a large proportion of mutations will be found on inner branches, i.e. be shared by many sequences, E[D i]=S/n. In other words, R S is such that smaller values are expected under a history of growth.
3. Methods
(i) Summary statistics and demographic model
We carried out coalescent simulations in ms (Hudson, Reference Hudson2002) to compare the power of a range of summary statistics to distinguish between the SNM and a history of exponential growth. In addition to D, D 2, R 2E and the new statistics defined above, F S, (Fu, Reference Fu1996) and H (Fay & Wu, Reference Fay and Wu2000) were considered. F S is based on the number of haplotypes in the sample and has previously been found to be more powerful than statistics based on the frequency distribution (Fu, Reference Fu1996; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). H was conceived as a test for the effect of selection on linked neutral sites (Fay & Wu, Reference Fay and Wu2000) and is not expected to have power to detect continuous growth. However, other demographic scenarios such as moderate bottlenecks may perturb genealogies in ways similar to genetic hitchhiking resulting in significant values of H. We assume that the population size has grown exponentially with rate α to its current size N 0:
Following standard practice, this exponential growth is incorporated through a re-scaling of time (Slatkin & Hudson, Reference Slatkin and Hudson1991). We define a rescaled time T coal relative to N 0 and α:
This represents the total amount of genetic drift that has occurred. It is convenient to define a growth rate relative to N 0 as A=2N 0α, which gives:
(ii) Power test
Critical values of 5% confidence for each statistic were determined from 10 000 replicate genealogies simulated under the SNM for each of a wide range of S values (1–250) (Hudson, Reference Hudson, Takahata and Clark1993; Braverman et al., Reference Braverman, Hudson, Kaplan, Langley and Stephan1995; Ramos-Onsins et al., Reference Ramos-Onsins, Mousset, Mitchell-Olds and Stephan2007). Genealogies from growing populations were simulated conditional on θ. For each replicate the alternative hypothesis of positive growth was tested by comparing the observed value of a statistic to the critical value given the observed S. Power was estimated as the proportion of 10 000 replicate genealogies for which a statistic was below its critical value in a one-tailed test. Power to reject the SNM was recorded for a large range of parameter combinations. We compared the performance of statistics for different growth rates, (0<A<50), sample sizes (n=10, 50) and values of θ (5–50). When varying θ, we chose a fixed value of A=8. This seems compatible with growth rates estimated from empirical data. For example, variation at silent sites in the Adhr region and X-linked genes in Drosophila pseudoobscura is consistent with A=7 (Schaeffer, Reference Schaeffer2002). While θ can be arbitrarily high for mitochondrial data, θ=20 may be unrealistic for nuclear loci in out-crossing species. Therefore, power was evaluated for a range of θ values (5–50) again keeping the growth rate fixed at A=8.
When using means and variances of summary statistics across loci, power was determined analogously to the single locus case. Critical values of 5% confidence of means and variances of statistics were determined from 10 000 replicate sets of loci with the exact same combination of S values. Although computationally expensive, this avoids making any assumptions about the distribution of mutation rates between loci. However, given that mutation rates vary along the genome assuming the same θ for all loci to simulate the alternative history of growth seems unrealistic and may lead to overestimation of power. We checked for the influence of heterogeneity in mutation rates on power by repeating the multilocus power tests with θ values drawn from a gamma distribution with α=2 (Pluzhnikov et al., Reference Pluzhnikov, Di Rienzo and Hudson2002) and a scale parameter equivalent to a mean of θ=20. This combination of growth and mutation rates is roughly comparable to mutation rate estimates for nuclear loci in Drosophila melanogaster (Galtier et al., Reference Galtier, Depaulis and Barton2000). As before we assumed no recombination within loci as well as absence of linkage between loci, i.e. replicate genealogies were simply treated as multiple loci.
(iii) Likelihood method
In practice, both θ and A are unknown, and their likelihood should, in principle, be estimated jointly. However, because of the non-independence of these two parameters, this is not a practical option. Following standard practice we alternated between maximum likelihood estimation of A and θ (Griffiths & Tavaré, Reference Griffiths and Tavaré1994). First a maximum likelihood estimate (MLE) for θ under the SNM was estimated using the program GENETREE (http://www.stats.ox.ac.uk/griff/software.htm). In a second step, this MLE for θ was fixed to run a likelihood surface for A. Finally, the MLE value for A was used to re-evaluate θ. This scheme yields two MLEs for θ for each replicate, one under the assumption of no growth and one given the most likely growth rate, which were compared in a likelihood ratio test (LRT). We did not find that the MLE estimates for A and θ improved upon repeated re-evaluation suggesting that a single round of estimation is sufficient for this moderate growth scenario. 100,000 runs were performed for each likelihood surface evaluation. Again, the proportion of replicate genealogies for which the null hypothesis could be rejected was taken as a measure of statistical power. Due to the long computing time, 100 replicates per parameter combination were used.
4. Results
(i) Single locus
In general, both the likelihood method and summary statistics have low power to detect a history of moderate (A<8) exponential growth for n=10 (Fig. 2). As expected, the likelihood method is most powerful overall, although its superiority is surprisingly small. For example, based on the LRT the SNM is rejected for 30% of genealogies simulated under exponential growth of A=4. In comparison, R S and R 2E detect this history of growth in 23% of cases (Fig. 2).
Power of summary statistics and likelihood method against exponential growth rate A=0–50. n=10, θ=20. Each point is based on 10 000 replicate simulations. The power of the likelihood method was estimated from 100 replicates (see large filled circles and error bars).
Consistent with previous results, F S, R 2E, and the new measure R S, are considerably more powerful than both D and D 2 (Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Ramirez-Soriano et al., Reference Ramirez-Soriano, Ramos-Onsins, Rozas, Calafell and Navarro2008). For θ=20, F S is the most powerful statistic. The new measure R S has consistently higher power than R 2E. As expected, H and X have no power to distinguish between the SNM and the growth case (not shown). However, the other two genealogical ratios perform surprisingly well. X 1 has higher power than D 2 and the power of X 2 is between that of R 2E and R S (Fig. 2). The complete lack of power of D for n=10 is somewhat surprising. Comparison with the result for n=50 (Fig. 3) reveals that its performance is strongly dependent on sample size. We ran additional simulations (not shown) and found that for n<15 extremely negative values of D are more likely under neutrality than under growth resulting in a rejection rate of the SNM of less than 5%. In other words, when n is small, the variance of D under neutrality is too large to detect exponential growth.
Power of summary statistics against exponential growth rate A=0–50. n=50, θ=20. Note the different range (0–1) on the y-axis compared with Fig. 2.
In general, all statistics have considerably higher power for n=50 (Fig. 3). Interestingly, it never reaches 100% even when growth is extreme (A=50). However, the relative effect of the sample size on power differs between statistics. For instance, X 1 improves relatively little in comparison to other measures. This is to be expected given that even small samples are likely to include the deepest split in the genealogy of the whole population (Saunders et al., Reference Saunders, Tavaré and Watterson1984). For n=10, the power of all statistics decreases for histories of extreme growth (A>25) (Fig. 2). This is due to the overall shortening of genealogies under rapid growth.
The mutation rate has a relatively small influence on power. In general, the power of all measures increases with θ (Fig. 4). However, the trajectories X 1 and F S level off while the power of the other statistics continues to improve with increasing values of θ. The power of F S is limited by the number of haplotypes (which cannot exceed n).
Power of summary statistics to detect a history of exponential growth (A=8) against θ. n=10.
To check how statistics are affected by the topological variance, genealogies simulated under the alternative history of growth were sorted according to the bipartition by the root and the proportion of significant values determined for each topology class. Figure 5 clearly shows that the two statistics based on π, D and R 2E as well as D 2 are sensitive to asymmetric topologies. The chance of observing a significant value increases markedly with topological asymmetry. This effect is most pronounced for D, which has no ‘power’ to reject the SNM unless genealogies are very asymmetric and growth is weak. In contrast, the dependency of X 1 on the rootward partition is relatively slight and in the opposite direction, i.e. the chance of rejecting the SNM is smaller for asymmetric genealogies (Fig. 5).
The effect of topological asymmetry on statistical power (simulation parameters as in 2). Genealogies of Fig. 2 were sorted according to the partition by the root (shown above plot). Only the most asymmetrical partition (9, 1) (a) and one other case (7, 3) (b) are shown. Results for the other three partitions were very similar to (b). Note that since lineages are exchangeable all asymmetrical partitions have the same probability P a=2/(n−1) (Tajima, Reference Tajima1983, eqn (2)).
(ii) Multiple loci
Compared with the relatively subtle effect both θ and n have on statistical power, increasing the number of loci improves power dramatically. In the mean-based test, all statistics apart from D have a power of close to 100% to detect a history of moderate exponential growth (A=8) for 10 loci. However, the relative performance of statistics changes slightly compared with the single locus case. Notably, X 2 has higher power than all other summary statistics (Fig. 6). The power of X is slightly lower than that of X 1 (not shown). Analogously to the results for a single locus, power increases both with more extreme growth scenarios and larger n (not shown).
Power of summary statistics to detect a history of growth A=8 using the mean across multiple loci against the number of loci, n=10 (A) and θ=20 (B). Assuming mutational rate heterogeneity (θ gamma distributed with α=2 and E[θ]=20).
As one may suspect, the increase in power with the number of loci is slower for the variance test. More importantly, the relative performance of statistics is very different. By far the most powerful statistic in the variance test is X 1 followed by D and X (Fig. 7). This indicates a general trade-off. Statistics with a high variance under the SNM have comparatively low power in the single-locus case and the mean test, but high power in the variance test and vice versa.
Power of summary statistics in the variance-based tests across multiple loci for three different growth rates (from left to right A=2, 4, 8). (A) θ=20. (B) Assuming mutational rate heterogeneity (θ gamma distributed with α=2 and E[θ]=20).
Allowing for heterogeneity in mutation rates between loci affects both the relative performance of summary statistics and their overall power. As one may expect, heterogeneity in θ generally results in a decrease in power. In the mean-based test, the three X statistics are most affected. However, in the variance test the performance of X 1 is little affected. This statistic even has slightly higher power when mutation rates vary between loci. This appears to be due to the non-normal distribution of X 1 under growth. Genealogies with more than one possible root-partition generally have a very low value of X 1, since we take an average over all possible partitions most of which will be associated with X 1=0.
5. Discussion
It is important to distinguish between the general limitations that genealogical and mutational stochasticity impose on demographic inference from genetic data and problems associated with particular methods. Two main conclusions emerge from comparing the performance of the new ‘genealogical statistics’ to classical neutrality tests and the LRT.
(i) General limits to demographic inference
The signatures that changes in population size leave in genealogies are typically subtle compared with the randomness of the ancestral process. Thus all methods have low power to distinguish between the SNM and histories of moderate growth in the single locus case. A surprising finding of this study was that the full likelihood method only works marginally better than the most powerful summary statistics. Changes in N e disproportionally affect the length of the basal branches of a genealogy. However, because these rootward branches also contribute most to the variance in total tree length, inferences based on a single locus will be weak at best. It is telling that the X statistics which only considers the last coalescence events in the history, outperform standard neutrality tests in the variance test when multiple realizations of this event, i.e. loci, are available. As has been argued before, most statistical power can be gained by increasing the number loci, which represent independent realizations of the ancestral process, rather than the sample size or the length of sequence (Felsenstein, Reference Felsenstein1992; Kliman et al., Reference Kliman, Andolfatto, Coyne, Depaulis, Kreitman, Berry, McCarter, Wakeley and Hey2000; Wakeley, Reference Wakeley2004).
(ii) Pairwise measures
Independent of the general limits to demographic inference, pairwise measures such as D have particularly low power to infer demography. This has been found in previous simulation studies, which consider other demographic scenarios such as strong bottlenecks and rapid logistic growth (Fu, Reference Fu1996; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Ramirez-Soriano et al., Reference Ramirez-Soriano, Ramos-Onsins, Rozas, Calafell and Navarro2008). The fundamental flaw of pairwise measures can be best understood in terms of the underlying genealogy. In contrast to selection and population structure, changes in N e on their own only alter the distribution of branch lengths without affecting the topology, which can be regarded as a random nuisance parameter. While the full topology can rarely be reconstructed, there is potentially a lot of topological information in sequence data. Thus, the challenge that any efficient inference method has to meet is to separate this topological information from the relevant branch length information while taking topological uncertainty into account. Tree-based methods such as lineage-through time plots clearly fall short of the latter because they rely on a fully resolved topology. Pairwise measures on the other hand simply ignore the confounding effect of the topology (Felsenstein, Reference Felsenstein1992). It is thus easy to see why D has power only when sample sizes are large. While increasing sample size adds increasingly shorter external branches and therefore little additional information, it does reduce the chance of extremely asymmetric bipartitions by the root which are responsible for much of the variance in π and hence D.
Perhaps worryingly, this sensitivity to the topology not only translates into a loss of statistical power but also means that negative D values may in fact be more informative about the topological asymmetry of the genealogy (which may be caused by other non-neutral forces, e.g. selection) underlying the sample than about past growth. In order to distinguish between the effects of selection and demography, topology needs to be separated from branch length information. One approach is to explicitly account for the topology information if possible. For instance, one could determine confidence intervals of statistics conditional on the bipartition by the root if this is known. Not surprisingly, this improves the power of D, but has little effect on statistics that are not based on π (not shown). The alternative is to use measures that are less sensitive to the topology. F S and other haplotype statistics have previously been shown to be more powerful than frequency spectrum statistics for this very reason (Depaulis et al., Reference Depaulis, Mousset and Veuille2003; Innan et al., Reference Innan, Zhang, Marjoram, Tavaré and Rosenberg2005). However, it has also been noted that F S sometimes behaves erratically (Fu & Li, Reference Fu and Li1993; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). As mentioned earlier, its power levels off with increasing θ (Fig. 4), because the sample size sets an upper bound to the number of haplotypes.
(iii) Recombination and topological uncertainty
The X statistics presented here fall somewhere in between tree-based methods and classical summary statistics. They exploit the fact that changes in population size disproportionally affect the relative length of the deepest branches in the genealogy and make use of topological information, without sacrificing the simplicity of the summary statistics framework. Given their high power in the multilocus case, how useful are such genealogical ratios in practice?
Recombination presents a fundamental problem to tree-based methods like the X statistics, which are defined only for non-recombining sequences. Similarly, likelihood methods that can deal with recombination are currently not available. To wrongly reconstruct trees from recombining data can potentially be severely misleading especially in the context of demographic inference. In fact, genealogical ratios similar to the ones presented here have been used to show that recombination can mimic the effect population growth has on the shape of inferred genealogies. Internal branches will appear relatively shorter and the tree overall more star-shaped (Schierup & Hein, Reference Schierup and Hein2000; Ramirez-Soriano et al., Reference Ramirez-Soriano, Ramos-Onsins, Rozas, Calafell and Navarro2008). Ideally one would like to model recombination explicitly when making demographic inferences. However estimates of recombination rates are usually associated with a large uncertainty. Furthermore, it is notoriously difficult to distinguish between recombination and back-mutations.
One approach to circumvent these problems is to test for recombination beforehand (e.g. using the four gamete test) and exclude recombinant regions from the analysis if necessary. One can then both condition on there being no within-locus recombination and afford to use more powerful statistics such as the ones presented here. This strategy of identifying non-recombining stretches of sequence is increasingly used to analyse multilocus data, e.g. Galtier et al. (Reference Galtier, Depaulis and Barton2000) or Jennings & Edwards (Reference Jennings and Edwards2005). Fortunately, many organisms appear to have lower recombination rates than model species such as Drosophila. For instance in a recent study on Australian birds only 6 out of 30 loci of intergenic sequence showed evidence for recombination (Jennings & Edwards, Reference Jennings and Edwards2005). How profitable this scheme is ultimately depends on the relative magnitude and distribution of recombination and mutation rates. Before the genealogical ratios can be used on multiple loci, which have been pruned to exclude recombinant stretches, both the potential bias of such pruning and the effect of undetected recombination events on the genealogical ratios need to be properly evaluated. Interestingly, our method of inferring the root does in itself constitute a test for recombination and may help to focus on those recombination events that matter to the statistical test.
A related problem concerns the infinite sites assumption. Although the algorithm we have developed to compute the X statistics takes topological uncertainty into account, ignoring the possibility of back-mutations may underestimate the length of basal branches (Baudry & Depaulis, Reference Baudry and Depaulis2003). Although this source of error has been ignored here it should in principle be possible to account for back-mutations considering that they are independent of the assumptions of the genealogical process. In fact, any mutational model can be used to define statistics analogous to the genealogical ratios presented here. The problem with more complicated mutation models is in estimating the basal topology needed to calculate these measures.
(iv) Conclusions
In summary, the results confirm that only the most extreme demographic events leave a sufficient signature to be detectable in single locus data. Still, instead of the excessive and often non-quantitative employment of mismatch distributions, phylogeographic studies could benefit from using more powerful statistics such as R S and R 2E to test demographic hypotheses. Conversely, population genetics studies of sequence data from multiple, unlinked loci could benefit from using summary statistics that incorporate genealogical information explicitly. When outgroup information is available and the assumptions of no within-locus recombination and infinite sites mutations can be justified, simple genealogical ratios are potentially more powerful than standard statistics. In taking the relative number of mutations found on specific parts of the genealogy as a measure of the degree of starshape, the demographic signal can be separated from irrelevant and confounding topological information. Extensions of this approach are feasible. For instance, one could consider the covariance between the number of basal and terminal mutations. Such simple statistics may be profitable for approximate likelihood or Bayesian approaches (Thornton & Andolfatto, Reference Thornton and Andolfatto2006). There remains a need to understand the effect of pruning and undetected recombination events on tree reconstruction in general and tree-based measures such as the X statistics presented here in particular.
Many thanks to N. Barton, P. Haddrill, J. Polechova, D. Charlesworth, Kai Zeng and D. Obbard for helpful advice and discussion. Thanks also to Kelly Dyer for help with Genetree. Detailed comments and valuable suggestions from two anonymous reviewers on an earlier version of this manuscript greatly improved this work. K.L. is funded by a studentship from the Biotechnology and Biological Sciences Research Council. J.K. is funded by EPSRC.
