Abstract
Background
We compared two methods of rooting a phylogenetic tree: the stationary and the nonstationary substitution processes. These methods do not require an outgroup.
Methods
Given a multiple alignment and an unrooted tree, the maximum likelihood estimates of branch lengths and substitution parameters for each associated rooted tree are found; rooted trees are compared using their likelihood values. Site variation in substitution rates is handled by assigning sites into several classes before the analysis.
Results
In three test datasets where the trees are small and the roots are assumed known, the nonstationary process gets the correct estimate significantly more often, and fits data much better, than the stationary process. Both processes give biologically plausible root placements in a set of nine primate mitochondrial DNA sequences.
Conclusions
The nonstationary process is simple to use and is much better than the stationary process at inferring the root. It could be useful for situations where an outgroup is unavailable.
Background
Several approaches for inferring a phylogenetic tree from the substitution patterns in multiply aligned sequences are available; they include maximum parsimony, distancebased, maximum likelihood and Bayesian methods [1]. Typically, the inferred tree is unrooted, because the explicit or implicit substitution process used is usually timereversible. An effective way to put the root on the unrooted tree is to perform a phylogenetic analysis on the sequences of interest together with an outgroup, which is a set of distantly related sequences [2,3]. If the ingroup is monophyletic in the combined phylogenetic tree, then the point where the outgroup touches the ingroup tree is the estimated root. The practical challenge is to find suitable outgroups, and if no such outgroup is available, then one is forced to root the tree using just the ingroup. Several such methods include the molecular clock and nonreversible substitution processes. It seems clear that compared to the outgroup method, the success of these methods is more dependent on the extent to which the accompanying assumptions about the substitution process are satisfied in the data. For example, the molecular clock method should work well if the lineages indeed evolved more or less at the same rate. Likewise, as shown by Huelsenbeck et al. [4], a nonreversible process is more likely to succeed the less reversible the real substitution process is.
The nonreversible substitution process, introduced by Yang [5], is stationary, i.e., the sequence composition is unchanged in time, and is equal to the equilibrium distribution of the rate matrix Q. The consensus is that it does not have enough power to discriminate among the candidate rooted trees. In this paper, we investigate a slightly more general, nonstationary process: in which the initial distribution π may not be the equilibrium distribution of the rate matrix Q. A priori, giving up stationarity is expected to produce a much better fit to data, since sequence composition is known to evolve, and should be accounted for. Indeed, substitution models where each branch has its own rate matrices had been used to resolve deep splittings in certain phylogenetic trees; see Yang and Roberts, and Galtier and Gouy [6,7]. Our process, which to our knowledge has not been investigated in this context, may be viewed as the simplest case of such nonstationary processes, with many fewer parameters. Thus, it can be used to decide whether the substitution processes on certain branches should be modeled differently. The input to our procedure is a multiple alignment and the topology of an unrooted binary tree. For each rooted tree associated with the given unrooted tree, we seek the maximum likelihood (ML) estimates of the branch lengths, π and Q. The rooted trees are then ranked in descending order of likelihoods. We model systematic variation in substitution rates among sites by assigning sites into several classes, and the relative rate for each class is estimated by ML; this is equivalent to the combined analysis framework of Yang [8].
We compared the ability of the stationary and nonstationary processes to place the root in three groups of species where the answer is considered wellknown: (1) human, chimpanzee and gorilla, (2) human, chimpanzee, gorilla and orangutan, (3) human, mouse, chicken and frog (xenopus laevis). The analyses were based on all available mitochondrial proteincoding genes, as well as two nuclear proteincoding genes. Next, we applied the methods to a set of primate mitochondrial DNA sequences.
Results
Verification studies
We fitted the nonstationary (NONSTA), stationary (STA) and reversible (REV) substitution models to all available mitochondrial proteincoding genes, as well as the nuclear genes albumin and cmyc, for three groups of organisms: (1) human, chimpanzee and gorilla, (2) human, chimpanzee, gorilla and orangutan, and (3) human, mouse, chicken and frog (xenopus laevis). The sequences were downloaded from Genbank and aligned using the CLUSTALW alignment of the amino acid sequences. Most alignments looked quite solid [see Additional files]. The beginning of the alignments for the genes COX1, CYTB, ND1 and ND6 were slightly adjusted. The root positions are assumed to be on the (1) gorilla, (2) orangutan, and (3) frog branch, respectively. The branches on a tree are referred to by the organism names, except for the case of four taxa, where there is an internal branch (Figure 1). For groups (2) and (3), it was assumed that human was most closely related to chimpanzee and mouse respectively; thus the unrooted tree is determined.
In group 1, the NONSTA and STA processes correctly placed the root in 8 and 6 genes respectively, out of 13 genes (Table 1). In group 2, NONSTA correctly placed the root in 9 genes out of 13 genes, compared to 2 genes for STA (Table 2). In group 3, NONSTA correctly placed the root in 11 genes out of 15 genes, compared to 7 genes for STA (Table 3). Furthermore, NONSTA gives stronger signal, or has better discriminative power: the highestscoring rooted tree often has noticeably higher log likelihoods than competing rooted trees; this is not so with STA. Thus, NONSTA is much better than STA in placing the root at the individual gene level. Combining the log likelihoods across genes yields overall evidence for the root placements. Table 4 shows that NONSTA is unambiguously correct in all three analyses, while STA only gets the root correctly in group 3, and the signal is weak.
Additional File 1. A text file containing the amino acid sequence alignments for group 1.
Format: ALN Size: 20KB Download file
Additional File 2. A text file containing the amino acid sequence alignments for group 2.
Format: ALN Size: 25KB Download file
Additional File 3. A text file containing the amino acid sequence alignments for group 3.
Format: ALN Size: 32KB Download file
Figure 1. Unrooted tree with four taxa The four branches adjacent to leaf nodes will be referred to by the corresponding taxon names.
Table 1. Human, chimpanzee and gorilla Loglikelihoods (rounded to closest integer) of the MLEs for three rooted trees under the nonstationary (NONSTA), stationary (STA) and reversible (REV) models. If NONSTA or STA places the root correctly, the corresponding log likelihood appears in bold.
Table 2. Human, chimpanzee, gorilla and orangutan Loglikelihoods (rounded to closest integer) of the MLEs for five rooted trees under the nonstationary (NONSTA), stationary (STA) and reversible (REV) models. If NONSTA or STA places the root correctly, the corresponding log likelihood appears in bold.
Table 3. Human, mouse, chicken and frog Loglikelihoods (rounded to closest integer) of the MLEs for five rooted trees under the nonstationary (NONSTA), stationary (STA) and reversible (REV) models. If NONSTA or STA places the root correctly, the corresponding log likelihood appears in bold.
Table 4. Combined analysis Combined log likelihoods over all genes under the nonstationary (NONSTA), stationary (STA), and reversible (REV) models. If NONSTA or STA places the root correctly, the corresponding log likelihood appears in bold.
The nuclear genes albumin and cmyc and three mitochondrial genes, COX1, COX2 and ATP6 from group 3 (with some mouse genes replaced with rat genes) were studied by Huelsenbeck et al. [4]. For these five genes, NONSTA and STA performed equally, getting all the correct root placements, except for ATP6, with NONSTA again noticeably more discriminative.
Primate mitochondrial DNA
Brown et al. and Yang [5,9] studied a set of mitochondrial DNA (mtDNA) sequences from human, chimpanzee, gorilla, orangutan, gibbon, crabeating monkey, squirrel monkey, tarsier and lemur. The topology of Yang's unrooted tree and the branch labels are shown in Figure 2. The mtDNA sequences consist of two proteincoding fragments, separated by three RNA genes. Thus, four site classes are required. Analysis with NONSTA shows that the root is most likely on the tarsier branch, followed closely by the lemur and "f" branches, and the corresponding log likelihoods are quite different from the others (see Table 5). Under STA, the most likely root placements are on the squirrel monkey and lemur branches. Thus, both processes give predictions that are consistent (NONSTA more than STA) with the idea that the root should be somewhere near tarsier and lemur. However, as observed before, NONSTA has much greater discriminative power, and fits the data much better, than STA.
Figure 2. Unrooted tree for nine primate mtDNA sequences The assumed unrooted tree is that presented in Yang [5]. The branches adjacent to leaf nodes are referred to by the corresponding organisms, while the interior branches are labelled a through f as indicated.
Table 5. Nine primates Loglikelihoods (rounded to closest integer) of the MLEs for 15 rooted trees under the nonstationary (NONSTA), stationary (STA) and reversible (REV) models.
Discussion
Our results confirmed earlier findings that the stationary process (STA) is not very good at discriminating among rooted trees corresponding to the same unrooted tree. In contrast, the nonstationary (NONSTA) process seems much more effective, with individual genes, and with combined genes. It is quite clear that the difference in log likelihoods between fitting STA and the reversible process (REV) is often small, and statistically insignificant, based on the likelihood ratio test, while those between NONSTA and STA, and between NONSTA and REV, are often large, and statistically very significant. Though the chisquare distribution may be inappropriate [10], it seems to be satisfatory in practice [11]. This indicates that NONSTA fits the data much better than STA and REV. Thus it appears that allowing an initial distribution that is uncoupled with the rate matrix gives a better description of the data, and that the greater capacity of NONSTA over STA at estimating the root placement may stem from the ability of NONSTA to allow for some amount of evolution in base composition.
Although Huelsenbeck et al.'s analysis using STA failed to place the root correctly in any of the genes albumin, cmyc, COX1, COX2 and ATP6, there are some differences between the analyses. The raw data were different: the rat albumin and cmyc genes were used by Huelsenbeck et al.; since mouse and rat are very similar, this is not likely to matter much. Secondly, the alignments were probably different, though since the sequences are quite similar, this should not be too important. It is plausible that most of the discrepancies between the results is due to the difference in the estimation procedure (maximum likelihood vs. Bayesian) and to the fact that in Huelsenbeck et al., site variation was modeled by the gamma distribution [12], whereas here we only accounted for the codon position effect.
Estimates of the relative rates are quite independent of the model used, and their relative magnitudes are largely within expectations. In particular, for group 3, the relative rates for codon positions 1, 2, and 3 fall between .2 and 1.1, .1 and .6, and 1.5 and 2.7 respectively. For all genes, the third codon position evolved the fastest, followed by the first and second positions. To gauge the contribution from the third codon position, we left out the corresponding bases in group 3 and reran the analysis with NONSTA. This gave the correct root placement in only three genes: albumin, cmyc and ND2, showing the usefulness of the third codon position in this dataset, despite its markedly higher substitution rates. We also found that the pairwise identity at the third codon positions for all genes in groups 3 ranges from 34% to 61%. Base composition being generally nonuniform, the expected pairwise identity at saturation (i.e., infinite evolutionary distance) is lower than 25%. This seems to indicate that the third codon position is not saturated, and hence the phylogenetic information from this position is not just the base composition at each taxon. In addition, the base composition at the third codon position for some genes is quite different from the other positions. Our model does not fit these genes as well as a model where separate processes are associated with the codon positions. Such a model will be investigated in future.
The NONSTA process is only slightly more complicated to apply, compared to the STA and REV processes. The fact that it works quite well in the verification studies and predicts biologically plausible roots for the nineprimate data demonstrates its utility and perhaps argues for its use in routine phylogenetic analysis. In any case, if no suitable outgroup is available, it could be worthwhile to try it. Though the NONSTA process is the most general timehomogeneous Markov process, it is still simplistic and imposes a severe constraint on the evolution of base composition: if two leaf nodes are at the same distance from the root, then the process stipulates that the corresponding sequences must have the same composition. This is patently unrealistic: once lineages split, they should evolve quite independently, and may explain the failure of the process at estimating the root placement for some genes. However, it is still valuable even if it does not always work, in that it can serve as a base from which exploration of richer models can be launched. For instance, one could identify lineages where the evolution significantly deviates from expectations, and then allow these lineages to have different rate matrices, which brings us closer to the very rich models of [6,7,13,14].
Conclusions
The nonstationary substitution process is simple to use, has much greater power at estimating the root compared to the stationary process, and also fits data much better than the stationary and reversible processes. It seems feasible to use this process in analyses where a suitable outgroup is not easily available. It is also a good starting point for conducting more sophisticated phylogenetic analysis with richer models.
Methods
Substitutions in DNA sequences are assumed to occur independently at each site according to a Markov process, i.e., given the present base, future substitutions are independent of past substitutions.
Furthermore, it is assumed that the process is timehomogeneous, i.e., substitution rates stay constant in time. As usual, the substitution rate from base a to b is the (a, b)entry in a 4 × 4 rate matrix Q; the diagonal entries are such that each row sums to 0. For any t > 0, the transition probability P(t) is given by P(t) = exp(Qt). Let π be a probability distribution on the DNA bases. The pair (π, Q) defines a substitution process on a rooted tree, as follows: pick a base at the root according to π, then run the substitution process according to Q down the tree, splitting into independent copies whenever a branching is encountered. The joint probability of the observed bases at the leaf nodes can be computed using almost exactly the same algorithm by [15].
There are two important special cases of the timehomogeneous process (π, Q). Associated with the rate matrix Q is a unique distribution π_{Q}, called the equilibrium distribution of Q, such that the matrix product π_{Q }× Q is the zero vector. The process (π_{Q}, Q) is stationary, i.e., the sequence composition remains unchanged through time, and is described by π_{Q}. Q is said to be reversible if it satisfies the detailed balance condition:
Π_{Q}Q = Q'Π_{Q}
where Π_{Q }is the diagonal form of π_{Q }and Q' is the transpose of Q. The process (π_{Q}, Q) is then reversible, i.e., statistically the process looks the same in forward and backward time. In particular, as shown in [15], the joint distribution of the leaf bases is the same regardless of where the root is placed on the tree. The reversible process is known as the REV or timereversible process in the molecular evolution literature [5,16,17]. Special cases of the REV process include those by Jukes and Cantor, Kimura, Felsenstein (two processes), Hasegawa, Kishino and Yano, and Tamura and Nei [15,1822]. The nonreversible stationary process was first explored by Yang [5], and subsequently by Huelsenbeck et al. [4]. Yang referred to this process as "unrestricted", but we use the abbreviation STA here. We shall refer to the nonstationary process as NONSTA. The numbers of free parameters in the NONSTA, STA and REV processes are respectively 15 (3 in π and 12 offdiagonal entries in Q), 12 (offdiagonal entries in Q) and 9 (3 in π_{Q }and 12 offdiagonal entries in Q, minus 6 detailed balance constraints). Since the models are nested, the likelihood ratio test can be used to assess the relative goodnessoffit of the MLEs. It is standard practice to allow only calibrated rate matrices, i.e., Q satisfies
so that a branch length is the average number of substitution events per site. We adopt this practice, and remark that for the nonstationary process (π, Q), with calibrated Q, since in general π ≠ π_{Q}, it is not true that the expected number of substitutions in 1 time unit is 1, but the difference gets arbitrarily small as time goes to infinity.
The sites in a DNA sequences can have very different substitution rates, the most wellknown example being coding sequences, where the third codon positions evolved much faster than the others because of the degeneracy of the genetic code. In cases where the assignment of sites into several classes is known in advance, such as a coding sequence, the easiest way to deal with it is to associate to class i an unknown positive number r_{i}, with the constraint that
where n_{i }is the number of sites in class i. The relative rate r_{i }either expands or shrinks the tree depending on whether it is more or less than 1. The constraint gives a new interpretation of a branch length: it is now the average over all sites of their expected number of substitutions. Thus, this approach is similar to [8]: effectively, the classes are treated as separate datasets. In this study, coding sequences are divided into three classes by codon position. In the last dataset consisting of nine primate mitochondrial sequences, an additional class is created to account for the RNAcoding bases. Another source of site variation is related to the threedimensional structure of the protein. For example, hydrophilic residues are usually exposed, hence tend to evolve faster than hydrophobic residues which are deeply buried. Our present approach does not model this and other less obvious sources of site variation. Possible remedies include using the gamma distribution [12] or the hidden Markov model [23].
Given a rooted tree relating aligned coding sequences, we seek the ML estimates of the branch lengths, the substitution parameters, and the relative rates. For other sequences, the relative rates are not estimated. Gradientbased methods are perhaps the most efficient at finding the maximum. The EM algorithm [24] is another possibility. We implemented the simplex method [25], which is slower but is less likely to be misled to local maxima than gradientbased methods. To further reduce the chance of being fooled by local maxima, different initial estimates were used, and the final estimates with the highest likelihood was picked. The initial estimates were obtained by first deriving a reversible rate matrix from a pairwise comparison of two sequences, then using the associated REV process to find the most likely branch lengths and relative rates; all pairwise comparisons were used in this study, so that, for example, four taxa give six initial estimates.
The estimation procedure was implemented in C, and the source code can be requested from the first author.
Authors' contributions
The idea was conceived by the first author and was inspired and refined by the second author. The first author composed the code and performed the data analysis.
Acknowledgements
Ziheng Yang kindly provided the alignment of the primate mtDNA sequences and entertained numerous enquiries. Several programs in the PAML package were used in the data analysis. We thank the referees for many comments and suggestions, and agree that more work is needed to investigate the notion of evolutionary time in nonstationary substitution processes, as well as the utility of our method in rooting larger trees.
References

Felsenstein J: Inferring Phylogenies. Sunderland, Massachusetts: Sinauer Associates, Inc; 2004.

Maddison WP, Donoghue MJ, Maddison DR: Outgroup analysis and parsimony.

Wheeler WC: Nucleic acid sequence phylogeny and random outgroups.

Huelsenbeck JP, Bollback JP, M LA: Inferring the root of a phylogenetic tree.
Sys Biol 2002, 51:3243. Publisher Full Text

Yang Z: Estimating the pattern of nucleotide substitution.
J Mol Evol 1994, 39:105111. PubMed Abstract

Yang Z, Roberts D: On the use of nucleic acid sequences to infer early branchings in the tree of life.
Mol Biol Evol 1995, 12:451458. PubMed Abstract  Publisher Full Text

Galtier N, Gouy M: Inferring pattern and process: Maximumlikelihood implementation of a nonhomogeneous model of DNA sequence evolution for phylogenetic analysis.
Mol Biol Evol 1998, 15:871879. PubMed Abstract  Publisher Full Text

Yang Z: Maximumlikelihood models for combined analyses of multiple sequence data.
J Mol Evol 1996, 42:587596. PubMed Abstract  Publisher Full Text

Brown WM, Prager EM, Wang A, Wilson AC: Mitochondrial DNA sequences of primates, tempo and mode of evolution.
J Mol Evol 1982, 18:225239. PubMed Abstract

Goldman N: Statistical tests of models of DNA substitution.
J Mol Evol 1993, 36:182198. PubMed Abstract

Yang Z, Goldman N, Friday A: Maximum likelihood trees from DNA sequences: A peculiar statistical estimation problem.

Yang Z: Maximumlikelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites.
Mol Biol Evol 1993, 10:13961401. PubMed Abstract  Publisher Full Text

Barry D, Hartigan JA: Statitstical analysis of hominoid molecular evolution.

Chang JT: Full reconstruction of Markov models on evolutionary trees: Identiflability and consistency.
Mathematical Biosciences 1996, 137:5173. PubMed Abstract  Publisher Full Text

Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach.
J Mol Evol 1981, 17:368376. PubMed Abstract

Lanave C, Preparata G, Saccone C, Serio G: A new method for calculating evolutionary substitution rates.
J Mol Evol 1984, 20:8693. PubMed Abstract

Tavaré S: Some probabilistic and statistical probles in the analysis of DNA sequences.
Lectures on Mathematics in the Life Sciences 1986, 17:5786.

Jukes TH, Cantor C: Evolution of protein molecules. In In Mammalian Protein Metabolism. Edited by Munro HN. Academic Press; 1969:21132.

Kimura M: A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences.
J Mol Evol 1980, 16:111120. PubMed Abstract

Kishino H, Hasegawa M: Converting distance to time: Application to human evolution.
Methods in Enzymology 1990, 183:550570. PubMed Abstract

Hasegawa M, Kishino H, Yano T: Dating the humanape splitting by a molecular clock of mitochondrial DNA.
J Mol Evol 1985, 22:160174. PubMed Abstract

Tamura K, Nei M: Estimation of the number of nucleotide substitutions in the control region of nitochondrial DNA in humans and chimpanzees.
Mol Biol Evol 1993, 10:512526. PubMed Abstract  Publisher Full Text

Felsenstein J, Churchill GA: A hidden Markov model approach to variation among sites in rate of evolution.
Mol Biol Evol 1996, 13:93104. PubMed Abstract  Publisher Full Text

Holmes IP, Rubin GM: An expectation maximization algorithm for training hidden substitution models.
J Mol Biol 2002, 317:753764. PubMed Abstract  Publisher Full Text

Nelder JA, Mead R: A simplex method for function minimization.