* Genome Atlantic, Department of Mathematics and Statistics, Dalhousie University, Halifax, Nova Scotia, Canada
Genome Atlantic, Canadian Institute for Advanced Research, Program in Evolutionary Biology, Department of Biochemistry and Molecular Biology, Dalhousie University, Halifax, Nova Scotia, Canada
Correspondence: E-mail: susko{at}mathstat.dal.ca.
![]() |
Abstract |
---|
Key Words: inconsistency rates across sites distance methods Neighbor-Joining molecular evolution phylogenetics
![]() |
Introduction |
---|
We consider here whether the topology will be consistently estimated when an incorrect model or no model has been used in constructing the pairwise distances. Sequence evolution in nature is an extremely complex and heterogeneous (across sites and across lineages) process. In most cases, therefore, when inferring phylogeny from sequences, models with simplifying assumptions regarding rate matrices, rates-across-sites, or stationarity of the rates-across-sites process across lineages will always be incorrect to some degree. For this reason, it is of critical importance to investigate the performance of methods when their models are misspecified in various ways. Inconsistency of a method of phylogenetic estimation implies poor performance as the number of sites gets large, and as such it is at least theoretically possible that a method that is inconsistent will perform reasonably with small numbers of sites. This is rarely the case in practice; inconsistency is usually a good indicator of situations in which poor performance can be expected, even with small numbers of sites. Indeed, for the forms of model misspecification similar to those considered through simulation with finite samples by Gaut and Lewis (1995) and Huelsenbeck (1995), the zones of inconsistency match up well with regions of poor performance by the methods considered by those authors.
The results presented indicate that serious difficulties of inference often arise in the presence of inconsistency. In many of the cases considered we show that with probability converging to one, a particular incorrect topology is estimated. As a consequence, with larger numbers of sites, there will be large bootstrap support for this incorrect topology. To illustrate briefly, we simulated data from the tree ((A : 0.8, C : 0.1) : 0.1, (B : 0.8, D : 0.1)) with a Jukes-Cantor (JC; Jukes and Cantor 1969) substitution process. As indicated later, with probability one, NJ with Hamming or p-distance will estimate the long branch topology, characterized by the neighbors (A, B), with a large enough number of sites. Even with 200 sites, the probability is very large that the long branch topology will be estimated. In 1,000 simulations the long branch topology (A, B) was estimated every time. More significantly, the average bootstrap support for the topology (A, B) over the 1,000 simulations was 81.3%, so that there was little information to suggest that there was a problem.
The analytical nature of the results presented lend insight into the reasons for inconsistency. The form of inconsistency most commonly discussed in the literature is long branch attraction in the context of the maximum parsimony method (Felsenstein 1978; Hendy and Penny 1989). Long branch attraction turns out to be the main form of inconsistency that we find occurring as well. The reason for this is that in so many of the cases where an overly simple model or no model is used to construct the distances, the limiting distances turn out to be a concave function of the actual distances. In contrast, when the limiting distances are convex functions of the actual distances, long branch inconsistency will not arise [a continuous concave (convex) function has negative (positive) second derivative everywhere or roughly speaking, has an upside-down (right-side up) bowl shape]. Instead, as we illustrate with an example, the form of inconsistency that can arise in this case is a "long branches repel" type of inconsistency.
Because it is analytically tractable and eases discussion, we primarily consider four taxon topologies. However a form of parallel rates-across-sites variation model misspecification is considered both analytically and, through simulation, with more taxa. Although idealized, the basic notion of parallel rates-across-sites variation is biologically relevant as a model for the evolutionary implications of differing functional constraints for different organisms. For instance, Steel, Huson, and Lockhart (2000) discussed a case in which changes in the numbers of variable sites in 18S ribosomal DNA from groups of Holometabolus insects provided artifactual support for Strepsitera as a sister taxon of Diptera. In this study, we consider a 42-taxon example where a group of intracellular parasitic eukaryotes, the microsporidia, are erroneously placed at the archeabacteria-eukaryote split in an elongation factor 1 phylogeny rather than with their fungal sister taxa (Hirt et al. 1999; Baldauf et al. 2000), possibly because of model misspecification. Whether the simplified model considered actually explains the erroneous placement or not, it illustrates that the analytical results for four taxon topologies have implications for larger numbers of taxa.
The results that we present are shown for the NJ method. In the four-taxon case, however, there is a close connection between the zones for NJ, LS, and ME. An immediate connection comes from the results of Gascuel (1994), which show that ME, with sum of LS branch lengths, and NJ give the same estimated trees. More generally, with concave limiting distances, we show that the zones of inconsistency are the same for NJ and ME, with ME criteria the sums of absolute values of branch lengths, or the sums of positive branch lengths, as well as for constrained LS (branch lengths are constrained to be non-negative). However, the zones of inconsistency for unconstrained LS (negative branch lengths are allowed) are larger than those of NJ. This is consistent with the findings of Kuhner and Felsenstein (1994), who show through simulation that corrections for negative branch lengths improve the accuracy of both topological estimation and branch length estimation.
A number of simulations studies have been conducted that investigate the biases and variability of various phylogenetic methods. These include studies by Gaut and Lewis (1995), Huelsenbeck (1995), Kuhner and Felsenstein (1994), and Bruno and Halpern (1999). On the basis of this and other work, it appears that there is widespread understanding that model misspecification can lead to biased estimation and inconsistency. Theoretical work that explicitly considers consistency in the case of Hamming distances includes reports by DeBry (1992) and Rzhetsky and Sitnikova (1996). Models where rates vary across sites and subtrees have been considered theoretically by Chang (1996) who shows that inconsistency can arise as branch length converges to zero. The work here can be considered an addition to the present work on the topic that (1) lends analytical insight into the reasons for the long branch form attraction that so frequently arises, (2) provides illustrations of alternative forms of model misspecification that can lead to inconsistency (specifically parallel rates-across-sites variation) that are not easily adjusted for, and (3) more closely considers the possible inconsistency of NJ as well as its relationship to LS and ME.
An outline of the article is as follows. The next two sections are technical, introducing assumptions, notation, and the main derivations leading to a zone of inconsistency. The following two subsections give examples of zones of inconsistency arising from a failure to model rates-across-sites and rates-across-subtrees processes. We also show in these sections that failure to account for rates-across-sites variation results in underestimation of evolutionary distances for a variety of models generalizing previous results for the JC model (Golding 1983), GTR model (Kelly and Rice, 1996) and confirming the simulation results of Kuhner and Felsenstein (1994). Other examples are then given of inconsistency corresponding to misspecifications in conventional model assumptions and resulting in long branches repelling. We conclude with a discussion of the relationships between the zones of inconsistency for NJ and those of LS and ME methods.
![]() |
Limiting Distances and Assumptions |
---|
We assume generally that dij(F(i,j), ) is a continuous function of F(i,j) and
. Let
be the probability that character states x and y arise in the two taxa across all sites. Since
converges to
as the number of sites goes to infinity, and by assumption
converges to
we get that
|
It will be valuable to think of incorrectly specified distances in terms of their dependence on the correctly specified distances. To do this, we must make the additional assumption that the distance between two taxa is a monotonically increasing transformation g(t) of the correctly specified distance or branch length. This assumption is equivalent to the statement that larger branch lengths should be associated with larger distances, which should always be the case for any reasonable distance. As an example, for the JC model, the Hamming distance along a branch is related to branch length through the transformation g(t) = 3/4 3 exp(4t/3)/4. In most of the examples that we consider, g(t) is a concave transformation of t, and so, unless stated otherwise, this is assumed as well. We will assume throughout most of the discussion that the true tree is ((A : b, C : a) : a, (B : b, D : a)):
|
|
![]() |
Zones of Inconsistency: Neighbor Joining |
---|
|
|
|
Substituting the limiting distances from table 1 into (1) and (2) we obtain that, with probability 1, as the number of sites gets large, (A, C) is preferred to (A, D) if
|
|
|
Thus h(b;a) is a decreasing function of b and so either h(b;a) < 0 for all a or h(b0(a);a) = 0 for some b0(a) in which case h(b;a) > 0 for all b < b0(a) and h(b;a) < 0 for b > b0(a). The zone of inconsistency can then be found by determining the solution, b0(a) of the equation h(b;a) = 0; we set b0(a) = if no solution exists. For fixed a the topologies estimated by NJ in the limit are as follows:
|
On the boundary, when b = b0(a), regardless of the number of sites, both of the topologies (A, C) or (A, B) will have positive probability of being estimated.
Although in most cases b0(a) must be determined numerically, for some distances an explicit solution for b0(a) can be calculated. For instance, if the true model is a JC model but Hamming distance is used, then b0(a) is determined as
|
|
|
![]() |
Effects of Ignoring Rates-Across-Sites Variation |
---|
|
|
To establish notation, we fix the pair of taxa i and j under consideration and leave as implicit the dependence of various quantities on i and j. Let t denote the true distance between the pair, the diagonal matrix with the stationary probabilities of the character states along its diagonal, and Q the rate matrix for the model under consideration. The number of additional parameters in Q varies depending on the model under consideration. For instance, for the JC model there are no additional parameters and for the F84 model there is a transition-transversion ratio. The matrix of pairwise probabilities, with xyth entry the probability of observing character states x an y for the two taxa, is then
E[eQrt], where expectation is with respect to the random rate r. The estimate
of
can be constructed from the observed frequencies or using ML estimation, but we assume that the estimated distance
and the rest of the parameters in the estimated rate matrix
are obtained through ML estimation without assuming a RAS model.
Now, the estimated distance always satisfies that
|
|
|
![]() |
Parallel Rates-Across-Sites Variation |
---|
We start by considering a four-taxon tree. The true tree is assumed to be of the form ((A : b, C : a) : a, (B : b, D : a). Rates are assumed constant along the middle branch and the branches leading to C and D. The rate variation is along the branches leading to A and B, and the rates-across-sites assigned to these branches are assumed to be positively correlated. Specifically, we assume that rates at any given site are assigned to the branches according to the probabilities of table 3.
|
|
|
|
|
Here
|
|
|
|
|
|
![]() |
The Archaebacteria-Eukaryotes Split |
---|
|
The model for rate change assumes that parallel changes occur along the branch A leading to the archaebacteria and the branch B leading to the microsporidia. Specifically, we assume the following:
The model of rate variation in branches A and B is the same model of parallel rate change that gave rise to long branch attraction in the four-taxon case. The main difference between this model and the model in the four-taxon case is that now a large proportion of 0 rate sites are allowed in all other branches of the tree. Under the model above, the expected proportion of sites with 0 rates is for A and B, and it is
00 +
(1
00) for all other branches. This provides a simple model for parallel rate variation that could arise because of a loss of function in genes. The cellular structure of microsporidia is highly degenerated, probably because of their parasitic lifestyle. It is also likely that the reduction of these organisms could be extended to the molecular level, since only 1,997 open reading frames appear to be encoded in the E. cuniculi genome (Katinka et al. 2001). Suppose for instance, that some sites in elongation factor 1
(EF-1
) are constrained in Eukaryotes but not in archaebacteria because of functions carried out by eukaryotic EF-1
but never present in archaebacterial EF-1
. If those functions were later lost in microsporidia, the corresponding sites would be free to vary, in contrast to the rest of the eukaryotes and in parallel with the archaebacteria. Eukaryotic EF-1
(eEF 1
) proteins are known to have both the core function as a translation elongation factor that is absolutely indispensable for all cellular systems and non-translational (auxiliary) functions that probably are absent in the archaebacterial orthologs (reviewed in Negrutskii and El'skaya 1998). Therefore it is plausible that microsporidian EF-1
proteins secondarily lost some eukaryotic EF-1
specific auxiliary functions during the reductive evolution in this intracellular parasitic lineage (Inagaki et al. 2004). By these parallel absence of EF-1
auxiliary functions in the microsporidian and archaebacterial proteins, the two sub-trees could show correlated patterns of rate variation, different from other eukaryotes (Inagaki et al. 2004).
To obtain a reasonable set of parameters, 00,
, and
for the simulation, we used the proportions of sites with no change among the subgroups of taxa of interest:
|
A rough estimate of , the probability of a 0 rate for the archaebacterial subtree, is then 0.258. Similarly 0.155 is an estimate for
00, and 0.436 is an estimate of
00 +
(1
00), the probability of a 0 rate in all other branches of the tree. Since 0.155 is an estimate for
00, an estimate for
can be obtained by solving
|
The tree used in the simulations is given in figure 4B. The branch with the microsporidia clade was placed between the fungi and the animals, but all other parts of the topology were fixed as in the topology estimated by NJ with uncorrected PAM distances (fig. 4A). If in fact parallel rate variation is present, then uncorrected ML estimates of the branch lengths for the topology in figure 4B will be biased. We used the regression approach described in the Supplementary Material online as a rough bias adjustment.
Because inconsistency implies that wrong topologies will be estimated with high probability when the number of sites is large, we simulated single data sets of 5,000 sites from the model above. Because of the uncertainty in the parameter and branch length estimates, simulations were conducted for a variety of settings for ,
00, and
near the estimated values of 0.333, 0.155, and 0.258. It was clear from these simulations that (1) larger values of the parameters
00 and
are more likely to give rise to microsporidia being placed with the archaebacteria, (2) changing
from its estimated value of 0.33 did not change results, and (3)
00 had to be close to
for microsporidia to be placed with the archaebacteria. The results for some of the parameter settings are indicated in table 4. For the estimated parameters (
,
00, and
set to 0.333, 0.155, and 0.258), the estimated topology places Entamoeba histolytica closest to archaebacteria, as do most of the parameter settings. Parameter settings with
00 large are the only ones that place microsporidia closest to the archaebacteria. The estimated topologies did not vary except in the placement of microsporidia. It was of interest to conduct simulations with 349 sites, because this was the size of the data set used. In this case, however, incorrect placement may be explained by sampling variation and, indeed, the estimated topologies differed more than just in their placement of microsporidia. Nevertheless, the results with respect to the clade closest to the archaebacteria were the same as indicated in table 4.
|
Regardless of the placement of microsporidia, the example illustrates that results for four taxon trees are likely to have implications for larger number of taxa. In the present case a key feature of the tree used in the simulation was the long branches leading to the microsporidia (0.77) and archaebacteria (0.30). The distance between the nodes corresponding to these branches was not small (0.25), and there were a number of intervening branches but, as indicated in figure 5, the shape of the tree was similar to the shape of the long branch four-taxon trees that have been considered throughout most of this article.
|
![]() |
Transition/Transversion Ratios |
---|
|
|
![]() |
An Example Where Long Branches Repel |
---|
|
|
Substituting these limiting distances into equations (1) and (2), we get that, in the limit, (A, C) is preferred to both (A, B) and (A, D), if
|
Suppose now that a JC substitution model is in effect with a gamma rate distribution having = 1. If the RAS model adopted is incorrectly specified as consisting of a mixture of a fixed rate and a set of invariable sites arising with probability 0.2, then
|
|
![]() |
Connections Between LS, ME, and NJ |
---|
|
For four taxa, it is shown in Gascuel (1994) that ME, with sum of LS branch lengths, and NJ give the same estimated trees and hence share the same zones of inconsistency. Other suggestions for the ME criterion are sums of absolute values of branch lengths (Kid and Sgaramella-Zonta 1971) and sums of positive branch lengths (Swofford et al. 1996). Note that these sums will always be larger than or equal to the unconstrained sum of branch lengths with equality if the branch lengths are all non-negative. Thus if the estimated ME topology with sum of unconstrained branch lengths has non-negative branch lengths, it will be the estimated ME topology under the other two criteria as well. It is shown in the Constrained Least-Squares Estimation subsection of the Supplementary Material online that this is the case with concave g. Hence the NJ and ME zones of inconsistency coincide regardless of the criterion used.
![]() |
Concluding Remarks |
---|
Investigations about inconsistency like the one conducted here provide information about the bias of estimation procedures and not the variance. Through simulation Gaut and Lewis (1995) plot regions where ML estimation performs well with finite samples under model misspecification. For small middle branches a, the shapes of their plots are similar to the shapes given here, suggesting that the poor estimation in these finite sample cases is due in large part to bias. However, for larger middle branches a, the zones where ML performs poorly are much less restricted, indicating, not surprisingly since the variance in distance estimates increases with distance, that variance provides additional problems for estimation in these cases.
For almost all of the examples considered here, it has been shown that inconsistency implies that a particular topology, different from the true one, will be converged on with probability 1. Because a particular topology is being converged on, with a large number of sites, measures of uncertainty, such as for instance bootstrap support, will suggest that with a high degree of certainty this wrong topology is the only one supported by the data. The fact that a particular topology is being converged on is important for studies investigating the use of simple distance measures as well. An overly simplistic distance measure that has a convex limiting distance function g(t), for instance, will converge on the wrong topology if the generating topology has long branches together, but it will converge on the right topology if the long branches are apart. In fact, such a distance measure might actually appear to be very efficient in estimating this topology. This is a point that has been considered in some detail in Bruno and Halpern (1999).
The examples involving parallel rate variation were important for several reasons. First, the application to the archaebacteria-eukaryotes EF-1 data indicate that the four-taxon results presented here likely have implications for larger, more complicated topologies. Second, although the model considered is simple, the basic notion of parallel rates-across-sites variation is biologically relevant. Finally, where many of the other forms of model misspecification can be adjusted for with appropriate distances, adjustments here would require a careful consideration of the processes taking place throughout the entire tree. The issues raised are similar to but different from those in Tuffley and Steel (1997). A correction for distances in the covarion model they consider is available as a rates-across-sites adjustment. Similarly, in the parallel rates-across-sites case, corrections could in principle be made through a rates-across-sites adjustment; in contrast, however, the rate distribution will not be constant across the taxa considered. It seems possible that the problem could be diagnosed on the basis of a pilot estimate of the tree, with separate rate distributions being fit to separate subsets of taxa, but we cannot pursue that further here. In any case, some form of iterative distance estimation, iterating between tree estimation and distance estimation with model adjustments, or perhaps ML estimation, would be required.
![]() |
Acknowledgements |
---|
![]() |
Footnotes |
---|
![]() |
Literature Cited |
---|
Baldauf, S. L., A. J. Roger, I. Wenk-Siefert, and W. F. Doolittle. 2000. A kingdom-level phylogeny of eukaryotes based on combined protein data. Science 290:972-977.
Bruno, W. J., and A. L. Halpern. 1999. Topological bias and inconsistency of maximum likelihood using wrong models. Mol. Biol. Evol. 16:564-566.
Cavalli-Sforza, L. L., and A. W. F. Edwards. 1967. Phylogenetic analysis: models and estimation procedures. Evolution 32:550-570.
Chang, J. T. 1996. Inconsistency of evolutionary tree topology reconstruction methods when substitution rates vary across characters. Math. Biosci. 134:189-214.[CrossRef][ISI][Medline]
Dayhoff, M. O., and R. V. Eck. 1968. Atlas of protein sequence and structure 19671968. National Biomedical Research Foundation, Silver Spring, Maryland.
Dayhoff, M. O., R. M. Schwartz, and B. C. Orcutt. 1979. A model of evolutionary change in proteins. Pp 345352 in M. O. Dayhoff, ed. Atlas of protein sequence and structure, volume 5, supplement 3. National Biomedical Research Foundation, Silver Spring, Md.
DeBry, R. W. 1992. The consistency of several phylogeny-inference methods under varying evolutionary rates. Mol. Biol. Evol. 9:537-551.[Abstract]
Fast, N. M., J. M. Logsdon, and W. F. Doolittle. 1999. Phylogenetic analysis of the TATA-binding protein (TBP) gene from Nosema locustae: evidence for a microsporidia-fungi relationship and spliceosomal intron loss. Mol. Biol. Evol. 16:1415-1419.
Felsenstein, J. 1978. Cases in which parsimony and compatibility methods will be positively misleading. Syst. Zool. 27:27-33.[ISI]
Felsenstein, J. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368-376.[ISI][Medline]
Fitch, W. M., and E. Markowitz. 1970. An improved method for determining codon variability in a gene and its application to the rate of fixation of mutations in evolution. Biochem. Genet. 4:579-593.[ISI][Medline]
Gascuel, O. 1994. Concerning the NJ algorithm and its unweighted version, UNJ. Pp 149170 in B. Mirkin, F. R. McMorris, F. S. Roberts, and A. Rzhetsky, eds. Mathematical hierarchies and biology, DIMACS Series in Discrete Mathematics and Theoretical Computer Science, Vol. 37. American Mathematical Society, Providence, R.I.
Gascuel, O., D. Bryant, and F. Denis. 2001. Strengths and limitations of the minimum evolution principle. Syst. Biol. 50:621-627.[CrossRef][ISI][Medline]
Gaut, B. S., and P. O. Lewis. 1995. Success of maximum likelihood phylogeny inference in the four-taxon case. Mol. Biol. Evol. 12:152-162.[Abstract]
Golding, G. B. 1983. Estimates of DNA and protein sequence divergence: an examination of some assumptions. Mol. Biol. Evol. 1:125-142.[Abstract]
Gu, X., and W. H. Li. 1996. A general additive distance with time-reversibility and rate variation among nucleotide sites. Proc. Natl. Acad. Sci. USA 93:4671-4676.
Hasegawa, M., H. Kishino, and T. Yano. 1985. Dating the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22:160-174.[ISI][Medline]
Hendy, M. D., and D. Penny. 1989. A framework for the study of evolutionary trees. Syst. Zool. 38:297-309.[ISI]
Huelsenbeck, J. P. 1995. Performance of phylogenetic methods in simulation. Syst. Biol. 44:17-48.[ISI]
Hirt, R. P., J. M. Logsdon, B. Healy, M. W. Dorey, W. F. Doolittle, and T. M. Embley. 1999. Microsporidia are related to fungi: evidence from the largest subunit of RNA polymerase II and other proteins. Proc. Natl. Acad. Sci. USA 96:580-585.
Inagaki, Y., E. Susko, N. M. Fast, and A. J. Roger. 2003. Covarion shifts cause a long branch attraction artifact that unites microsporidia and archaebacteria in EF-1 phylogenies. Mol. Biol. Evol. 21:1340-1349.[ISI]
Inagaki, Y., and W. F. Doolittle. 2000. Evolution of the eukaryotic translation termination system: origins of release factors. Mol. Biol. Evol. 17:882-889.
Jukes, T. H., and C. R. Cantor. 1969. Pp 21123 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press. New York.
Katinka, M. D., S. Duprat, E. Cornillot, G. Metenier, F. Thomarat, G. Prensier, V. Barbe, E. Peyretaillade, P. Brottier, and P. Wincker, et al. 2001. Genome sequence and gene compaction of the eukaryote parasite Encephalitozoon cuniculi. Nature 414:450-453.[CrossRef][ISI][Medline]
Keeling, P. J., and W. F. Doolittle. 1996. Alpha-tubulin from early-diverging eukaryotes and the evolution of the tubulin family. Mol. Biol. Evol. 13:1297-1305.
Kelly, C., and Rice, J. 1996. Modeling nucleotide evolution: A heterogeneous rate analysis. Math. Biosci. 133:85-109.[CrossRef][ISI][Medline]
Kid, K. K., and L. A. Sgaramella-Zonta. 1971. Phylogenetic analysis: concepts and methods. Am. J. Hum. Genet. 23:235-252.[ISI][Medline]
Kimura, M. 1980. A simple method for estimating evolutionary rate of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16:111-120.[ISI][Medline]
Kuhner, M. K., and J. Felsenstein. 1994. A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates. Mol. Biol. Evol. 11:459-468.[Abstract]
Li, W. H. 1994. Molecular evolution. Sinauer Associates, Sunderland, Mass.
Lockhart, P. J., D. Huson, U. Maier, M. J. Fraunholz, Y. van de Peer, A. C. Barbrook, C. J. Howe, and M. A. Steel. 2000. How molecules evolve in eubacteria. Mol. Biol. Evol. 15:1183-1188.
Negrutskii B. S., and A. V. El'skaya. 1998. Eukaryotic elongation factor-1: structure, expression, functions and possible role in aminoacyl-tRNA channeling. Proc. Nucl. Acid Res. Mol. Biol. 60:47-78.
Nei, M. 1987. Molecular evolutionary genetics. Columbia University Press, New York.
Nei, M., R. Chakraborty, and P. A. Fuerst. 1976. Infinite allele model with varying mutation rate. Proc. Natl. Acad. Sci. USA. 73:4164-4168.[Abstract]
Page, R. D. M., and E. C. Holmes. 1998. Molecular evolution: a phylogenetic approach. Blackwell Science, Oxford.
Penny, D., B. J. McComish, M. A. Charleston, and M. D. Hendy. 2001. Mathematical elegance with biochemical realism: the covarion model of molecular evolution. J. Mol. Evol. 53:711-723.[CrossRef][ISI][Medline]
Rzhetsky, A., and M. Nei. 1993. Theoretical foundation of the minimum-evolution methods of phylogenetic inference. Mol. Biol. Evol. 10:1073-1095.[Abstract]
Rzhetsky, A., and T. Sitnikova. 1996. When is it safe to use an oversimplified substitution model. Mol. Biol. Evol. 13:1255-1265.[Abstract]
Saitou, N., and M. Nei. 1987. The Neighbor-Joining method: a new method for reconstructing evolutionary trees. Mol. Biol. Evol. 4:406-425.[Abstract]
Steel, M., D. Huson, and P. J. Lockhart. 2000. Invariable sites models and their use in phylogeny reconstruction. Syst. Biol. 49:225-232.[CrossRef][ISI][Medline]
Swofford, D. L., G. L. Olsen, P. J. Waddell, and D. M. Hillis. 1996. Phylogenetic inference. Chapter 11 in D. M. Hillis, C. Moritz, and B. K. Mable, eds. Molecular systematics, 2nd edition. Sinauer Associates, Sunderland, Mass.
Tuffley, C., and M. A. Steel. 1997. Modelling the covarion hypothesis of nucleotide substitution. Math. Biosci. 147:63-91.[CrossRef][ISI]
Uzzel, T., and K. W. Corbin. 1971. Fitting discrete probability distributions to evolutionary events. Science 173:1089-1096.
van de Peer, Y., A. Ben Ali, and A. Meyer. 2000. Microsporidia: accumulating molecular evidence that a group of amitochondriate and suspectedly primitive eukaryotes are just curious fungi. Gene 246:1-8.[CrossRef][ISI][Medline]
Waddell, P. J., and M. A. Steel. 1997. General time-reversible distances with unequal rates across sites: Mixing and inverse gaussian distributions with invariant sites. Mol. Phyl. Evol. 8:398-414.[CrossRef][ISI][Medline]
Yang, Z. 1994. Estimating the pattern of nucleotide substitution. J. Mol. Evol. 39:105-111.[ISI][Medline]