Measuring the degree of starshape in genealogies – summary statistics and demographic inference

Summary

The degree of starshape of a genealogy is readily detectable using summary statistics and can be taken as a surrogate for the effect of past demography and other non-neutral forces. Summary statistics such as Tajima's D and related measures are commonly used for this. However, it is well known that because of their neglect of the genealogy underlying a sample such neutrality tests are far from ideal. Here, we investigate the properties of two types of summary statistics that are derived by considering the genealogy: (i) genealogical ratios based on the number of mutations on the rootward branches, which can be inferred from sequence data using a simple algorithm and (ii) summary statistics that use properties of a perfectly star-shaped genealogy. The power of these measures to detect a history of exponential growth is compared with that of standard summary statistics and a likelihood method for the single and multi-locus case. Statistics that depend on pairwise measures such as Tajima's D have comparatively low power, being sensitive to the random topology of the underlying genealogy. When analysing multi-locus data, we find that the genealogical measures are most powerful. Provided reliable outgroup information is available they may constitute a useful alternative to full likelihood estimation and standard tests of neutrality.

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. (2003), Hamblin et al. (2004), Haddrill et al. (2005), Ometto et al. (2005) and Thornton & Andolfatto (2006). 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, 1992). 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, 1991; Harpending, 1994; Schneider & Excoffier, 1999). 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, 1989; Galtier et al., 2000).

Approaches to demographic inference fall into three broad categories; for a review see Emerson et al. (2001). 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é, 1994; Kuhner et al., 1995). 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., 1995; Pybus et al., 2002). 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, 1989) and Fu & Li's D (hereafter referred to as D ₂) (Fu & Li, 1993), the haplotype distribution, e.g. Fu's F _S (Fu, 1996; Innan et al., 2005), or the mismatch distribution, e.g. the raggedness statistic (Slatkin & Hudson, 1991). 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., 1995; Fu, 1996; Ramos-Onsins & Rozas, 2002; Depaulis et al., 2003; Ramirez-Soriano et al., 2008). However, joint analysis of multiple loci is not only necessary to distinguish between selective and demographic events (Galtier et al., 2000) 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., 1998; Reich et al., 1999) but are now routinely used to analyse sequence data from multiple loci (Pluzhnikov et al., 2002; Haddrill et al., 2005; Heuertz et al., 2006) or even species (Hickerson et al., 2006).

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, 2002; Innan et al., 2005; Ramirez-Soriano et al., 2008). Earlier, Felsenstein made a theoretical argument for the inferiority of pairwise measures (Felsenstein, 1992). Their large variance under neutrality arises both from their sensitivity to the last coalescence event and the random genealogical topology (Tajima, 1983). 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, 1992), 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, 1992), 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, 1993; Tajima, 1989; Fay & Wu, 2000) θ=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 ₂ 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 ₂ is to distinguish between two classes of mutations: those found on terminal branches, η_e and those on internal branches, η_i (Fig. 1) (Fu & Li, 1993). 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 ₂) (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.

Fig. 1. 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 ₂, 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 ₂, 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, 1997, Appendix).

One possibility is to only consider the shorter of the two rootward branches that has exactly length T ₂. Thus the number of mutations found on this branch, η_Rmin, over θw constitutes a very simple measure of starshape.

(1)

Such genealogical ratios have first been employed to study the effect of balancing selection on plant incompatibility loci (Uyenoyama, 1997). 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:

(2)

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:

(3)

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 (2ⁿ⁻¹−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:

(4)

(Ramos-Onsins & Rozas, 2002, eqn (2)). R _2E has proven superior to a wide range of summary statistics in detecting histories of bottlenecks (Ramos-Onsins & Rozas, 2002). 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:

(5)

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, 2002) 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 ₂, R _2E and the new statistics defined above, F _S, (Fu, 1996) and H (Fay & Wu, 2000) 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, 1996; Ramos-Onsins & Rozas, 2002). H was conceived as a test for the effect of selection on linked neutral sites (Fay & Wu, 2000) 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 ₀:

(6)

Following standard practice, this exponential growth is incorporated through a re-scaling of time (Slatkin & Hudson, 1991). We define a rescaled time T _coal relative to N ₀ and α:

(7)

This represents the total amount of genetic drift that has occurred. It is convenient to define a growth rate relative to N ₀ as A=2N ₀α, which gives:

(8)

(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, 1993; Braverman et al., 1995; Ramos-Onsins et al., 2007). 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, 2002). 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., 2002) 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., 2000). 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é, 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).

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 ₂ (Ramos-Onsins & Rozas, 2002; Ramirez-Soriano et al., 2008). 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 ₁ has higher power than D ₂ and the power of X ₂ 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.

Fig. 3. 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 ₁ 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., 1984). 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 ₁ 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).

Fig. 4. 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 ₂ 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 ₁ 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).

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, 1983, 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 ₂ has higher power than all other summary statistics (Fig. 6). The power of X is slightly lower than that of X ₁ (not shown). Analogously to the results for a single locus, power increases both with more extreme growth scenarios and larger n (not shown).

Fig. 6. 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 ₁ 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.

Fig. 7. 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 ₁ 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 ₁ under growth. Genealogies with more than one possible root-partition generally have a very low value of X ₁, since we take an average over all possible partitions most of which will be associated with X ₁=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, 1992; Kliman et al., 2000; Wakeley, 2004).

(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, 1996; Ramos-Onsins & Rozas, 2002; Ramirez-Soriano et al., 2008). 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, 1992). 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., 2003; Innan et al., 2005). However, it has also been noted that F _S sometimes behaves erratically (Fu & Li, 1993; Ramos-Onsins & Rozas, 2002). 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, 2000; Ramirez-Soriano et al., 2008). 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. (2000) or Jennings & Edwards (2005). 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, 2005). 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, 2003). 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, 2006). 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.