Mitochondrial Genes and Mammalian Phylogenies: Increasing the Reliability of Branch Length Estimation

Patrice Showers CorneliGo,* and Ryk H. Ward*{dagger}

*Department of Biology, University of Utah; and
{dagger}Institute of Biological Anthropology, University of Oxford


    Abstract
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 
Since branch lengths provide important information about the timing and the extent of evolutionary divergence among taxa, accurate resolution of evolutionary history depends as much on branch length estimates as on recovery of the correct topology. However, the empirical relationship between the choice of genes to sequence and the quality of branch length estimation remains ill defined. To address this issue, we evaluated the accuracy of branch lengths estimated from subsets of the mitochondrial genome for a mammalian phylogeny with known subordinal relationships. Using maximum-likelihood methods, we estimated branch lengths from an 11-kb sequence of all 13 protein-coding genes and compared them with estimates from single genes (0.2–1.8 kb) and from 7 different combinations of genes (2–3.5 kb). For each sequence, we separated the component of the log-likelihood deviation due to branch length differences associated with alternative topologies from that due to those that are independent of the topology. Even among the sequences that recovered the same tree topology, some produced significantly better branch length estimates than others did. The combination of correct topology and significantly better branch length estimation suggests that these gene combinations may prove useful in estimating phylogenetic relationships for mammalian divergences below the ordinal level. Thus, the proper choice of genes to sequence is a critical factor for reliable estimation of evolutionary history from molecular data.


    Introduction
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 
During the past decade, mitochondrial sequence data have been used to resolve long-standing phylogenetic problems. Mitochondrial DNA provides the ideal test case for these kinds of analyses, because all genes are inherited as a single unit and thus have a single evolutionary history. It has become clear, however, that phylogenies based on single genes (<2 kb) lack reliability and resolution and that different mtDNA genes can generate different trees (Cao et al. 1994Citation ; Waits 1996Citation ; Zardoya and Meyer 1996Citation ). Conflicts among gene trees arise because short or highly conserved sequences may lack the phylogenetic information to detect short, deep branches that distinguish the correct phylogeny from an incorrect one. Single genes also may be subject to selective constraints that differentially bias the phylogenetic signal away from the correct tree. Although it has been well established that these sources of error lead to distorted tree topologies (Huelsenbeck 1997Citation ; Sullivan and Swofford 1997Citation ), excessive error can also yield unreliable estimates of branch lengths, and hence incorrect estimates of divergence times, even when the topology is correct.

The increased efficiency of sequencing technology holds the promise that in the future, much longer sequences will routinely be generated. The potential gain in precision is already apparent for mitochondrial DNA (mtDNA) for which the existence of the entire genomic sequences for a substantial number of taxa (e.g., Xu, Janke, and Arnason 1996Citation ; Arnason, Gullberg, and Janke 1997Citation ; Penny and Hasegawa 1997Citation ) has begun to yield more informative phylogenies, superseding those based on single mtDNA genes. The reliability to be gained by using the entire mtDNA genome has yet to be comprehensively evaluated, since there have been few statistically robust analyses comparing the accuracy of branch lengths estimated by the whole genome with that of those estimated from smaller mtDNA segments (but see Hillis et al. 1992Citation ; Graybeal 1994Citation ).

It is often not appreciated that maximum-likelihood (ML) criteria include both branch length and topology estimates, even though poor estimation of either can interfere with efforts to build a reliable tree. We show that trees that have identical topologies but are estimated from different sequences can differ significantly from one another simply because branch length estimates differ. Hence, the estimates of the timing and extent of evolutionary divergence may be significantly better for some mtDNA sequences than for others. Furthermore, some sequences that recover incorrect topologies provide significantly better branch length estimates than do sequences that recover the correct topology.

In addition, to detect statistically significant differences among trees, it is important that the test of the hypothesis is set up properly. So far, most analyses of the phylogenetic utility of shorter sequences (Cao et al. 1994, 1998Citation ; Zardoya and Meyer 1996Citation ) have not used statistically informative comparisons. These earlier studies fit the expected tree to the individual gene sequence and generally failed to show that ML gene trees are significantly different from the expected tree. Single gene sequences lack sufficient phylogenetic information to provide a powerful contrast, and hence will fail to discriminate among the alternative hypotheses. Ranking genes using tests that lack the power to detect the significance of differences among ML scores is meaningless, since the differences may be attributable to random error. Instead, the more appropriate procedure is to fit trees defined by single-gene (and longer) sequences to the full sequence of all 13 genes (11 kb). When this is done, we can demonstrate that trees built from single genes are significantly poorer than the expected tree generated from the entire mtDNA protein-coding sequence.

Nineteen mammalian sequences, including six hominoid primates (gibbon, orangutan, gorilla, human, and two chimpanzee species), six ungulates, and three carnivores, were chosen because there is general consensus about the phylogenetic relationships within (although not among) these orders. Phylogenies for modern eutherians can be problematic, because most families within modern orders appeared rapidly during the early Paleocene (Flynn and Galiano 1982Citation ; Alroy 1999Citation ) or during a subsequent mid-Cenozoic radiation (Simpson 1945Citation ; Radinsky 1982Citation ; Novacek 1990Citation ). The difficulty arises because sequences must contain enough change to reveal closely spaced branching events without accumulating homoplasies that obscure the true branching pattern; thus, sequences must be relatively conserved and they also must be relatively long. Even the entire mtDNA sequence may be inadequate for resolving phylogenetic trees among very divergent taxa (Sullivan and Swofford 1997Citation ; Naylor and Brown 1998Citation ; this study). On the other hand, the subordinal relationships of taxa included in this study are known. Consequently, these modern eutherians that have an evolutionary history of not much more than 70 Myr represent an ideal set of taxa to evaluate the performance of mtDNA genes.


    Materials and Methods
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 
Sequence Data
Thirteen protein-coding gene sequences were extracted from each of 19 complete mammalian mitochondrial DNA genomes (fig. 1 ). Taxa were chosen to include a range of distances: the 19 species were distributed among 8 orders, 11 families, and 17 genera. Sequences include those of the gibbon, Hylobates lar (X99256, Arnason et al. 1996Citation ); the orangutan, Pongo pygmaeus (D38115; Horai et al. 1992Citation ); the gorilla, Gorilla gorilla (X93347; Xu and Arnason 1996Citation ); the human, Homo sapiens (X93334; Arnason, Xu, and Gullberg 1996Citation ); the common chimpanzee, Pan troglodytes (X93334; Arnason, Xu, and Gullberg 1996Citation ); the pygmy chimpanzee, Pan paniscus (D38116; Horai et al. 1995Citation ); the cow, Bos taurus (V00654; Anderson et al. 1982Citation ); the blue whale, Balaenoptera physalus (X61145; Arnason, Gullberg, and Widegren 1991Citation ); the fin whale, Balaenoptera musculus (X72204; Arnason and Gullberg 1993Citation ); the horse, Equus caballus (X79547; Xu and Arnason 1994Citation ); the Indian rhinoceros, Rhinoceros unicornis (X97336; Xu, Janke, and Arnason 1996Citation ); the white rhinoceros, Ceratotherium simum (Y07726; Xu and Arnason 1997Citation ); the harbor seal, Phoca vitulina (X63726; Arnason and Johnsson 1992Citation ); the gray seal, Halichoerus grypus (X72004; Arnason et al. 1993Citation ); the cat, Felis catus (U20753; Lopez, Cevario, and O’Brien 1996Citation ); the mouse, Mus musculus (J01420; Martens and Clayton 1979Citation ); the rat, Rattus norvegicus (X14848; Gadeleta et al. 1989Citation ); the armadillo, Dasypus novemcinctus (Y11832; Arnason, Gullberg, and Janke 1997Citation ); and the hedgehog, Erinaceus europaeus (X88898; Kretteck, Gullberg, and Arnason 1995Citation ).



View larger version (28K):
[in this window]
[in a new window]
 
Fig. 1.—The "expected" mammalian tree: maximum-likelihood tree (Felsenstein, F84, model lnL = -118,892.07) inferred from the concatenated sequence of 13 mitochondrial DNA protein-coding genes (11,262 bp). Branches are labeled to correspond with the X-axis in figure 3 . All bootstrap values are 100% except for branches 18 (99%), 31 (89%), and 28 (58%). The tree is not significantly different from an alternative that places the armadillo just outside of the ferungulate clade (Arnason, Gullberg, and Janke 1997Citation ; Janke 1997Citation ) or from ones with monophyletic ungulates or with cetartiodactyls sister to the carnivores. The addition of parameters for rate heterogeneity and the proportion of invariant sites (Hasegawa-Kishino-Yano+invariant+gamma, HKY+I+G, model lnL = -107,458.14) improves the model but does not improve discrimination among the alternative trees.

 
Alignment (Pileup; GCG 1997Citation ) of each gene across 19 sequences was straightforward and produced only a few gaps. These were removed along with neighboring nucleotides so that only full codons were retained in the edited sequences. For construction of the expected tree, the 13 individual gene sequences were concatenated into a sequence comprising 11,262 bp. Some redundancy resulted, because ND4 and ND4L overlap, as do the ATPase6 and ATPase8 and the ND5 and ND6 genes.

Data Analysis
ML trees were constructed from each of the 13 protein-coding gene sequences, from the full concatenated sequence, and from seven different combinations of genes. Several criteria were used to combine the seven multiple-gene sequences, which ranged in length from about 2 to 3.5 kb. For example, cytochrome b, widely used for mammalian phylogenies, was combined with the nearby ND5 gene, for which primers are readily available (CYTB/ND5). We used sequence similarity plots (GCG 1997Citation ) to empirically determine the relative rates of sequence evolution and then assembled two combinations of sequences with similar levels of conservation (COI/COII, conserved; ND3/ND4/ND4L, intermediate) and one that mixed fast, intermediate, and conserved genes (ND2/ND1/COI, respectively). Finally, some sequences comprising genes that recovered topologies similar to that of the expected tree (ND1/ND2, ND5 [first and second positions only]/ND4) were assembled and compared with sequences of similar lengths with poorly performing genes (ND3/ND4L/COII/COIII).

While earlier studies focus on the utility of genes for resolving very deep divergences (Cao et al. 1994, 1998Citation ; Graybeal 1994Citation ; Zardoya and Meyer 1996Citation ; Naylor and Brown 1997Citation ), we increased the density of subordinal lineages dating from the mid-Cenozoic and assumed that our results would be more applicable to recent mammalian divergences. The time depth for which each gene is appropriate differs, because the rate of sequence evolution differs among genes (Simon et al. 1994Citation ); less conserved genes (e.g., ND2) may perform better for more recent splits than the highly conserved genes (COI, COII, and COIII). To determine whether a different set of genes is phylogenetically informative for subordinal splits, we used the same set of parameters (DNAML default) (Felsenstein 1993Citation ) used by Zardoya and Meyer (1996)Citation for higher-level relationships. Subsequently, we checked our results using a more complex model (see below).

The default DNAML model, hereinafter referred to as the Felsenstein and Hasegawa-Kishino-Yano model (F84/HKY; Swofford et al. 1996Citation ), accommodates unequal base frequencies (empirical) and a 2:1 transition/transversion ratio (the ML estimate for these data and for this model, which assumes rate homogeneity; details available from P.S.C.). Each tree was constructed using a global search strategy for 10 runs, each with a randomly different order of species input. We estimated the degree to which the data support each clade with majority-rule consensus trees built from ML trees generated from 100 bootstrap replicates (Felsenstein 1985Citation ) of the sequence data. Sequences were analyzed with and without the third codon positions.

To ensure that conclusions from this study were not artifacts of a poor model (e.g., Saccone, Reyes, and Pesole 1998Citation ), we expanded the model to include across-site rate heterogeneity parameters and invariant sites (Hasegawa-Kishino-Yano+invariant+gamma [HKY+I+G]; Swofford et al. 1996Citation ) using ML methods (Swofford 1999Citation ). Since rate heterogeneity and invariant-sites parameters have been shown to significantly increase the likelihood for mammalian trees (Huelsenbeck 1997Citation ; Sullivan and Swofford 1997Citation ; unpublished data), we estimated ML parameters (transition/transversion ratio, base composition bias, the proportion of invariant sites, and gamma rate parameters) separately for each sequence. Subsequently, we used the more complex models to verify the results from the simpler (F84/HKY) model. For each sequence, model choice was based on nested ML tests (Posada and Crandall 1998Citation ) of 12 models ranging in complexity from the 4-parameter F81 to the 10-parameter general-time-reversible (GTR+I+G) model (Swofford et al. 1996Citation ).

To verify that results held across methods, we constructed neighbor-joining (NJ) trees for the combined sequence data and for individual genes from Tajima-Nei distance matrices (GCG 1997Citation ), and using parsimony criteria (Maddison and Maddison 1993Citation ), we compared the fit of the 13 gene tree topologies with the expected tree topology and reasonable alternatives.

The fit of each gene or multiple-gene tree to the full sequence was compared with that of the expected tree using the KH test (Kishino and Hasegawa 1989Citation ), in which the mean difference across sites in log-likelihoods (lnL) was compared with the standard error (SE) of the differences. If the log-likelihood statistic (-lnL difference/SE) was greater than 1.96, then the alternative trees were considered significantly different from one another. For ML methods, the KH test does not (and cannot) explicitly compare alternative topologies (branching pattern), but, rather, compares the branch length estimates conditioned on each topology. The topology that fits the optimal branch length estimate is the ML tree topology (Felsenstein 1981Citation ). Thus, alternative topologies have different branch length estimates, and these differences may or may not be statistically significant. Similarly, two trees that have exactly the same topology but are estimated from different genes can have significantly different branch length estimates. For this reason, for each gene sequence, we used a nested set of likelihood scores to evaluate deviations in branch length estimates not attributable to topology differences. In a manner analogous to the analysis of variance (ANOVA), which partitions the total log-likelihood statistic into individual components of variation, we separated the component of the log-likelihood difference due to changes in topology alone from the additional component due to branch length differences unrelated to topology.

Partitioning the likelihood score is straightforward using DNAML (Felsenstein 1993Citation ) and PAUP (Swofford 1999Citation ), which allow the user to specify a set of trees for the KH test and to either specify (constrain) or re-estimate (unconstrain) the branch lengths for each tree. For the first set of tests, we fit the topology and branch lengths estimated from the single genes to the full sequence and compared the resulting log-likelihood score with that of the expected tree. The full log-likelihood statistic is the difference between log-likelihood scores for the ML tree derived from the full sequence (the expected tree) and the ML tree estimated from the single-gene sequence. Therefore, the size of the full log-likelihood statistic reflects deviations of the gene tree from the expected tree that are attributable to topological differences (if any) and to other branch length differences.

Second, we fit each gene tree topology to the full sequence with the branch lengths unconstrained. Thus, the branch lengths were re-estimated conditioned on the gene tree topology, and the resulting change in log-likelihood deviation reflects only the change in topology. By constraining the topology but allowing the branch lengths to be reoptimized to those that best fit the full- sequence data, we were able to separate the component of the log-likelihood deviations attributable to the fit of the gene tree topology from the total deviation.

To explore the performance of genes in estimating specific branches of the expected tree, we used PHYLIP-generated estimates to compare the accuracy, across genes, of individual branch lengths by determining which estimates fell within the relatively narrow 95% confidence intervals (CIs) estimated for the expected tree. These tests of individual branch length fit cannot provide reliable measures of statistical significance and are meant to be exploratory rather than inferential. Work is now in progress to statistically evaluate hypotheses suggested by this analysis.


    Results
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 
The Expected Tree
The ML tree (fig. 1 ; F84 -lnL = 118,892.1, HKY+I+G = 107,456.8) recovered by the concatenated sequences of all 13 protein-coding genes is consistent with those of other molecular phylogenies using complete mtDNA sequences (Zardoya and Meyer 1996Citation ; Penny and Hasegawa 1997Citation ). It is not significantly different from the tree based only on the first and second codon positions, in which Xenarthra (armadillo; see also Arnason, Gullberg, and Janke 1997Citation ; Janke 1997Citation ; Penny et al. 1999Citation ) is the sister group to the Ferungulata (including carnivores, ungulates, and others; Simpson 1945Citation ), rather than appearing in the more basal position (Novacek 1990Citation ), as it does in the ML tree. Nor is it significantly better than six other alternatives, including one with the armadillo as the outgroup taxon, one with the rodents as sister to the primates (Cao et al. 1998Citation ), and one with monophyletic ungulates. These sequence data, regardless of the model, do not resolve these interordinal relationships, but resolve, with strong statistical confidence, relationships within the orders (ML tests and bootstrap herein; Hendy and Penny 1993Citation ; unpublished spectral analyses). Because the eight trees differ very little in likelihood ({Delta}lnL/SE < 1.96 in each case) and not at all in subordinal relationships, the choice of trees with which to compare the performance of the shorter sequences of protein-coding genes is somewhat arbitrary. The ML tree is in most features noncontroversial, so we refer to it as the "expected" tree.

ML Topologies
No single gene sequence recovers the expected topology (fig. 1 ) regardless of method (ML and NJ, with or without the third- position sites). Among single genes, only the ND1 and ND2 topologies provide a reasonable fit (F84, HKY+I+G, {Delta}lnL/SE = 1.27 and 1.78, respectively) to the full 11-kb sequence (table 1 ). All other single-gene trees, regardless of model, are significantly worse than the expected tree ({Delta}lnL/SE > 1.96 in each case). In contrast, three combinations of genes (Cytb/ND5, ND5 [third positions excluded]/ND4, and ND3/ND4/ND4L) recover the same topology as the full sequence does (table 2 ). The sequence comprising ND2, ND1, and COI recovers a slightly different tree topology, but it is not significantly different from the expected tree (F84, {Delta}lnL/SE = 1.55; HKY+I+G, {Delta}lnL/SE = 0.03; the armadillo joins just outside the ferungulate clade, with a monophyletic rather than paraphyletic ungulate clade). Three combinations of gene sequences (COI/COII, ND3/ND4L/COII/COIII, and ND1/ND2) fail to recover the expected tree, yielding significantly worse topologies than the expected tree (table 2 ; {Delta}lnL/SE > 1.96 in each case). However, by accommodating rate heterogeneity, it can be seen that the ND3/ND4L/COII/COIII and ND1/ND2 differences are spurious ({Delta}lnL/SE = 1.37 and 1.64, respectively). Among the combinations of genes studied, only the COI/COII topology differs significantly from the expected topology.


View this table:
[in this window]
[in a new window]
 
Table 1 Kishino-Hasegawa (KH) Tests for the Fit of Trees Generated by Single-Gene Sequences to the Concatenated Sequence

 

View this table:
[in this window]
[in a new window]
 
Table 2 Kishino-Hasegawa (KH) Tests for the Fit of Trees Generated by Multiple-Gene Sequences to the Concatenated Sequence

 
Branch Length Analysis
In contrast to the well-studied effect of sequence length (Felsenstein 1988Citation ; Graybeal 1998Citation ) on the quality of phylogenetic estimation, we wished to address the effect of gene choice on the quality of branch length estimation. Among trees generated from sequences of similar lengths, some that were topologically identical to the expected tree fit the full sequence data less well than some with the wrong topology (table 2 ). For example, even though the ND5/ND4 sequence (2,565 bp) recovers the correct topology, the full log-likelihood deviation from the expected tree is large. Thus, the deviation in fit from the expected tree (F84, {Delta}lnL/SE = 8.67; HYK+I+G, {Delta}lnL/SE = 13.75) reflects the poor quality of the branch length estimates compared with the expected tree. In contrast, although the topology of the tree generated by the combined ND3/ND4L/COII/COIII sequence (2,150 bp) differs from that of the expected tree, the branch length estimates are significantly better (HYK+I+G, {Delta}lnL/SE = 2.43) than those of the topologically correct ND5/ND4 tree.

Figure 2 plots branch lengths estimated from each of the seven gene combinations against the expected branch lengths estimated from the full sequence. In principle, with perfect agreement between the two sequences, all points will lie on the diagonal line, while large deviations indicate large branch length differences. For example, both the ND1/ND2/COI and the Cytb/ND5 estimates are significantly closer to the expected branch lengths than are those of ND4/ND5 (table 3 ; {Delta}lnL/SE = 5.70 and 5.18, respectively). Although the ND1/ND2/COI tree topology differs slightly from the expected tree, the branch length deviations are so small that the overall fit (F84, {Delta}lnL/SE = 3.01; HKY+I+G, {Delta}lnL/SE = 5.70) is better than that for any other gene tree. Similarly, the ND3/ND4L/ND4 tree is significantly better than the topologically identical ND4/ND5 tree ({Delta}lnL/SE = 4.40), even though ND3/ND4L/ND4 is shorter by 555 bp. Although very important, sequence length alone is not responsible for poor branch length estimation or for large likelihood deviations from the expected branch lengths. Regardless of length, some sequences provide a better fit than others do.



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 2.—Branch lengths estimated from multiple-gene sequences plotted against expected branch lengths. Points along the line of equality represent branch lengths in perfect agreement with the branch lengths estimated from the entire sequence of all genes. Points above (or below) the line indicate branches that are too long (or too short) relative to branch lengths estimated from the full sequence. Despite recovering the correct tree topology, the ND4/ND5 sequence appears to systematically underestimate branch lengths, generating branch length estimates that are significantly poorer ({Delta}lnL/SE = 5.70) than those of the ND1/ND2/COI sequence. Even larger deviations in branch length estimates characterize the topologically incorrect COI/COII tree.

 

View this table:
[in this window]
[in a new window]
 
Table 3 Pairwise Comparisons of Fit Among Multiple-Gene Trees

 
Although the log-likelihood statistic measures the overall fit of the tree to the full sequence, it does not provide measures of fit for the lengths of individual branches within the tree. We compared individual branch length estimates from shorter sequences with the CI for the full sequence (table 4 ). Regardless of time depth, most single genes failed to recover reasonable estimates. In contrast, certain combinations of gene sequences recover high percentages of branch lengths within the CI for the expected (ML) tree. Among the 18 branches describing divergences occurring since the mid-Cenozoic, 83% estimated from the ND1/ND2/COI sequence are within the expected CI. On the other hand, the ND4/ND5 tree had the correct tree topology, but nearly all branch length estimates fell well outside the CI for the expected tree, resulting in an overall fit that is significantly worse than that of the ND1/ND2/COI tree ({Delta}lnL/SE = 5.70).


View this table:
[in this window]
[in a new window]
 
Table 4 Comparison of Branch Lengths Estimated from Single-Gene Sequences with Branch Lengths Estimated from the Full Sequence of All Genes

 

    Discussion
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 
Our approach to evaluating the quality of the phylogenetic information within individual genes was to fit each of the 13 gene trees to the full 11.2-kb sequence (table 1 ) and to compare log-likelihood scores between each gene tree and the expected tree. Conversely, Zardoya and Meyer (1996)Citation fit the expected tree to individual genes. Statistically, this is a much less powerful use of the KH test, because single-gene sequences are shorter than the full sequence and thus are less likely to distinguish the many erroneous trees from the expected tree. Regardless of model complexity (with or without rate heterogeneity assumptions), 13 different gene trees are recovered by the 13 mtDNA protein-coding genes, and none are statistically distinguishable from the expected tree (fig. 3 ; with a multiple-test correction, even ATPase8 and COII differences are not significant). In reality, all except ND1 and ND2 are significantly worse than the expected tree, and these differences are easily detectable using the full sequence (table 1 ), implying that single genes lack the power to distinguish good trees from bad ones. For phylogenies with mid-Cenozoic or earlier divergence times, we believe that ranking single genes by their phylogenetic utility is virtually meaningless and that great caution should be used in interpreting single-gene phylogenies.



View larger version (38K):
[in this window]
[in a new window]
 
Fig. 3.—A, Log-likelihood deviations (F84/HKY model) comparing the fits of the expected tree to the individual gene sequences. Dark bars indicate log-likelihood differences for gene trees calculated for this study; light bars were calculated from Zardoya and Meyer (1996)Citation (table 1 ). The dashed line indicates the critical value (1.96). Alternative trees with values larger than 1.96 are considered significantly different from the expected tree. B, Log-likelihood deviations (HKY+I+G model) comparing the fits of the expected tree to the individual gene sequences. The critical values (dashed lines) show that all single genes, except perhaps COII and ATPase8, lack the phylogenetic information to detect the differences between trees that are easily detectable using the full sequence. Corrected for multiple tests, neither the ATPase8 nor the COII tree provides a significantly better fit to its respective sequence than does the expected tree.

 
We have shown that some trees, estimated from longer sequences, are better estimates of the expected tree than are single-gene trees. However, even among those trees with the correct topology, some are significantly better than others. In particular, the large log-likelihood deviation of the ND4/ND5 tree, resulting from inaccurate branch length estimates, indicates that this sequence may not adequately represent the phylogenetic information of the full sequence of protein-coding genes and thus may be particularly unreliable for inferring unknown phylogenies. In general, the largest branch length deviations are associated with the longest branches (table 4 ; detailed branch length data available from P.S.C.). On the other hand, even small differences on very short branches can prevent convergence on the correct tree. For example, the branch that connects the perrisodactyls to the carnivores (fig. 1 , #28) is so short that the two alternative arrangements of the ferungulates cannot be distinguished from the ML tree ({Delta}lnL/SE < 1.96) (see also Waddell et al. 1999Citation ). Only one of the arrangements is correct, but the branch length error, although small, is evidently large enough to prevent resolution of this part of the tree. In general, the ability to retrieve the expected topology necessarily depends on branch length fit, so gene sequences that fail to give the right branch length may converge on a biologically unlikely tree.

Regardless of the model (HKY or HKY+I+G), the Cytb/ND5 and ND1/ND2/COI trees fit the phylogenetic information in the full 11.2-kb sequence significantly better than do those of any other sequence tested. In particular, ND1/ND2/COI recovered a tree with the smallest overall log-likelihood deviation from the expected tree (table 2 ), suggesting that these three genes might be especially useful for resolving other mammalian phylogenies with post-Eocene or later divergence times. This is when the many modern eutherian groups radiated (Radinsky 1982Citation ; Prothero 1994Citation ), and the depth and rapidity of the divergences renders them difficult to detect with fossil data (Simpson 1945Citation ) or with short sequences (Waits 1996Citation ; Dragoo and Honeycutt 1997Citation ). It is for this time frame, roughly the past 50 Myr, that our results are probably most useful.

We have less confidence that any of these sequences should be used to resolve ancient divergences (see also Graybeal 1994Citation ). Indeed, it is doubtful that any combination of mtDNA genes will be adequate, since even the entire 11.2-kb sequence failed to fully resolve the relationships among ferungulate orders. Given the continuing controversy regarding relationships among the mammalian orders (Arnason, Gullberg, and Janke 1997Citation ; Sullivan and Swofford 1997Citation ; Saccone, Reyes, and Pesole 1998Citation ; Waddell, Okada, and Hasegawa 1999Citation ) and the relative lack of discrimination of the entire genome for these evolutionary depths, it is likely that resolution will require addition of nuclear sequences (Graybeal 1994Citation ), as well as denser taxa sampling at the interordinal level (Graybeal 1998Citation ; Poe 1998Citation ).

On the other hand, for subordinal divergences, ML trees from the full sequence and from some multiple-gene sequences contain all of the expected topological relationships regardless of the model. This may be because ML techniques are relatively robust to violations of model assumptions (Felsenstein 1993Citation ; Swofford et al. 1996Citation ) or because the effects of adding rate parameters are smaller for shorter branches (Yang 1996Citation ). Our preliminary results (unpublished data) suggest that with the inclusion of rate parameters, branch length estimates for divergences of less than 20 MYA are nearly identical to those estimated with the rate homogeneity model and are about 20% longer for divergences occurring about 50 MYA. Only for divergences occurring during the late Cretaceous to early Tertiary (just 2 of 34 branches) are branch lengths dramatically underestimated by the rate homogeneity model. With rate parameters included, these two branches are twice as long as those estimated without the extra parameters.

The cost of the additional parameters is reduced power to distinguish alternative trees from the ML tree, as well as increased CI size on the branch lengths (unpublished data). However, it is important to insure that the simple model of sequence evolution used in analysis is not biased in favor of an incorrect tree. Our preliminary site likelihood comparisons and skew tests (unpublished data) suggest that characters that distinguish one topology from another using the F84/HKY model are simply spurious differences corrected by inclusion of across-sites rate heterogeneity parameters. Similarly, Sullivan and Swofford (1997)Citation identified a larger set of alternative trees than did others using the same data (e.g., Saccone, Reyes, and Pesole 1998Citation ), possibly because of reduced power, but more likely because their more realistic model of sequence evolution discarded spurious differences. Regardless, the lack of power certainly does not warrant using the wrong model, as some imply (Saccone, Reyes, and Pesole 1998Citation ). Such a practice merely provides increased confidence in the wrong tree. Even if a larger set of alternative trees is obtained, the extra parameters render the branch length estimates less biased and therefore more useful for correctly inferring evolutionary history. The next logical step, now underway, is to identify the set of parameters that balance optimal branch length estimation against loss of power for phylogenies dating from the mid-Cenozoic to Recent.


View this table:
[in this window]
[in a new window]
 
Table 4 Extended

 

    Acknowledgements
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 
We are grateful to K. Johnson, K. Crandall, R. Honeycutt, and an anonymous reviewer for comments on the manuscript. D. L. Swofford and D. Posada kindly provided beta versions of PAUP (version 4b2) and MODELTEST, respectively. We also thank D. Wolstenholme, F. Adler, S. George, E. Rickart, and especially J. Seger for helpful discussions during the course of the research and L. Okun for computing advice. The Department of Biology, University of Utah, and the Institute of Biological Anthropology, University of Oxford, provided computer support.


    Footnotes
 
Rodney Honeycutt, Reviewing Editor

1 Keywords: maximum likelihood, mammalian phylogeny, complete mitochondrial genome, branch length estimates, mtDNA protein-oding genes. Back

2 Address for correspondence and reprints: Patrice Showers Corneli, Department of Biology, University of Utah, Salt Lake City, Utah 84112. E-mail: corneli{at}biology.utah.edu Back


    literature cited
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results
 Discussion
 Acknowledgements
 literature cited
 

    Alroy, J. 1999. The fossil record of North American mammals: evidence for a Paleocene evolutionary radiation. Syst. Biol. 48:107–118.[ISI][Medline]

    Anderson, S., M. H. de Bruijn, A. R. Coulson, I. C. Eperon, F. Sanger, and I. G. Young. 1982. Complete sequence of the mammalian mitochondrial genome. J. Mol. Biol. 156:683–716.[ISI][Medline]

    Arnason, U., and A. Gullberg. 1993. Comparison between the complete mtDNA sequence of the blue and fin whale, two species that can hybridize in nature. J. Mol. Evol. 37:312–322.[ISI][Medline]

    Arnason, U., A. Gullberg, and A. Janke. 1997. Phylogenetic analyses of mitochondrial DNA suggest a sister group relationship between Xenarthra (Edentata) and ferungulates. Mol. Biol. Evol. 14:762–768.[Abstract]

    Arnason, U., A. Gullberg, A. Janke, and X. Xiufeng. 1996. Pattern and timing of evolutionary divergences among Hominoids based on analysis of complete mtDNAs. J. Mol. Evol. 43:650–661.[ISI][Medline]

    Arnason, U., A. Gullberg, E. Johnsson, and C. Ledje. 1993. The nucleotide sequence of the mitochondrial DNA molecule of the grey seal, Halichoerus grypus, and a comparison with mitochondrial sequences of other true seals. J. Mol. Evol. 37:323–330.[ISI][Medline]

    Arnason, U., A. Gullberg, and B. Widegren. 1991. The complete nucleotide sequence of the mitochondrial DNA of the fin whale, Balaenoptera physalus. J. Mol. Evol. 33:556–568.[ISI][Medline]

    Arnason, U., and E. Johnsson. 1992. The complete mitochondrial sequence of the harbor seal, Phoca vitulina. J. Mol. Evol. 34:493–505.[ISI][Medline]

    Arnason, U., X. Xu, and A. Gullberg. 1996. Comparison between the complete mitochondrial DNA sequence of Homo and the common chimpanzee based on nonchimeric sequence. J. Mol. Evol. 42:145–152.[ISI][Medline]

    Cao, Y., J. Adachi, A. Janke, S. Pääbo, and M. Hasegawa. 1994. Phylogenetic relationships among Eutherian orders estimated from inferred sequences of mitochondrial proteins: instability of a tree based on a single gene. J. Mol. Evol. 39:519–527.[ISI][Medline]

    Cao, Y., A. Janke, P. J. Waddell, M. Westerman, O. Takenaka, S. Murata, N. Okada, S. Pääbo, and M. Hasegawa. 1998. Conflict Among individual mitochondrial proteins in resolving the phylogeny of eutherian orders. J. Mol. Evol. 47:307–322.[ISI][Medline]

    Dragoo, J. W., and R. L. Honeycutt. 1997. Systematics of mustelid-Like carnivores. J. Mammal. 78:426–443.[ISI]

    Felsenstein, J. 1981. Evolutionary trees from DNA Sequences: a maximum likelihood approach. J. Mol. Evol. 17:368–376.[ISI][Medline]

    ———. 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39:783–791.

    ———. 1988. Phylogenies from molecular sequences: inferences and reliability. Annu. Rev. Genet. 22:521–565.[ISI][Medline]

    ———. 1993. PHYLIP (phylogenetic inference package). Version 3.5c. Distributed by the author, Department of Genetics, University of Washington, Seattle.

    Flynn, J. M., and H. Galiano. 1982. Phylogeny of early Tertiary Carnivora with a description of a new species of Protictis from the Middle Eocene of northwestern Wyoming. Am. Mus. Novit. 2725:1–64.

    Gadeleta, G., G. Pepe, G. DeCandia, C. Quadgliariello, E. Sbisa, and C. Saccone. 1989. The complete mitochondrial genome: cryptic signals revealed by comparative analysis between vertebrates. J. Mol. Evol. 28:497–516.[ISI][Medline]

    GCG. 1997. Wisconsin package. Version 9.1. Genetics Computer Group, Madison, Wis.

    Graybeal, A. 1998. Is it better to add taxa or characters to a difficult phylogenetic problem? Syst. Biol. 47:9–17.

    ———. 1994. Evaluating the phylogenetic utility of genes: a search for genes informative about deep divergences among vertebrates. Syst. Biol. 43:174–193.[ISI]

    Hendy, M. D., and D. Penny. 1993. Spectral analysis of phylogenetic data. J. Classif. 10:5–24.[ISI]

    Hillis, D. M., J. J. Bull, M. E. White, M. R. Badgett, and I. J. Molineux. 1992. Experimental phylogenetics: generation of a known phylogeny. Science 255(5044):589–591.

    Horai, S., K. Hayasaka, R. Kondo, K. Tsigane, and N. Takhata. 1995. Recent African origin of the modern humans revealed by complete sequences of hominoid mitochondrial DNAs. Proc. Natl. Acad. Sci. USA 92:532–536.

    Horai, S., Y. Satta, K. Hayasaka, R. Kondo, T. Inoue, T. Ishida, S. Hayashi, and N. Takahata. 1992. Man’s place in Hominoidea revealed by mitochondrial DNA genealogy. J. Mol. Evol. 35:32–43.[ISI][Medline]

    Huelsenbeck, J. P. 1997. Is the Felsenstein Zone a flytrap? Syst. Biol. 46:69–74.

    Janke, A. 1997. The complete mitochondrial genome of the wallaroo (Macropus robutus) and the phylogenetic relationship among Montremata, Marsupialia, and Eutheria. Proc. Natl. Acad. Sci. USA 94:1276–1281.

    Kishino, H., and M. Hasegawa. 1989. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea. J. Mol. Evol. 29:170–179.[ISI][Medline]

    Kretteck, A., A. Gullberg, and U. Arnason. 1995. Sequence analysis of the complete mitochondrial DNA molecule of the hedgehog, Erinaceous europeaus, and the phylogenetic position of the Lipotyphla. J. Mol. Evol. 41:952–957.[ISI][Medline]

    Lopez, J. V., S. Cevario, and S. J. O’Brien. 1996. Complete nucleotide sequence of the domestic cat (Felis catus) mitochondrial geneome and a transposed mtDNA tandem repeat (Numt) in the nuclear genome. Genomics 33:229–246.

    Maddison, W. P., and D. R. Maddison. 1993. MacClade 3.0. Sinauer, Sunderland, Mass.

    Martens, P. A., and D. A. Clayton. 1979. Mechanism of mitochondrial DNA replication of the light-stranded origin of replication. J. Mol. Evol. 135:327–351.

    Naylor, G. J. P., and W. M. Brown. 1997. Structural biology and phylogenetic estimation. Nature 388:527–528.

    ———. 1998. Amphioxus mitochondrial DNA, chordate phylogeny, and the limits of inference based on comparisons of sequences. Syst. Biol. 47:61–76.[ISI][Medline]

    Novacek, M. J. 1990. Morphology, paleontology, and the higher clades of mammals. Pp. 507–543 in H. H. Genoways, ed. Current mammalogy. Vol. 2. Plenum Press, New York.

    Penny, D., and M. Hasegawa. 1997. The platypus put in its place. Nature 387:549–550.

    Penny, D. L., M. Hasegawa, P. J. Waddell, and M. D. Hendy. 1999. Mammalian evolution: timing and implications from using the LogDeterminant transform for proteins of differing amino acid composition. Syst. Biol. 48:76–93.[ISI][Medline]

    Poe, S. 1998. Sensitivity of phylogeny estimation to taxon sampling. Syst. Biol. 47:18–31.[ISI][Medline]

    Posada, D., and K. A Crandall. 1998. MODELTEST: testing the model of DNA substitution. Bioinformatics 14:817–818.

    Prothero, D. R. 1994. The Eocene-Oligocene transition: paradise lost. Columbia University Press, New York.

    Radinsky, L. B. 1982. Evolution of skull shape in carnivores 3. The origin and early radiation of the modern carnivore families. Paleobiology 8:177–195.

    Saccone, C., A. Reyes, and G. Pesole. 1998. Complete mitochondrial DNA sequence of the fat dormouse, Glis glis: further evidence of rodent paraphyly. Mol. Biol. Evol. 15:499–505.[Abstract]

    Simon, C., F. Frati, A. Beckenbach, B. Crespi, H. Liu, and P. Flook. 1994. Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann. Entomol. Soc. Am. 87:651–701.[ISI]

    Simpson, G. G. 1945. The principles of classification and a classification of the mammals. Bull. Am. Mus. Nat. Hist. 85:1–350.[ISI]

    Sullivan, J., and D. L. Swofford. 1997. Are guinea pigs rodents? The importance of adequate models in molecular phylogenetics. J. Mamm. Evol. 4:77–86.

    Swofford, D. L. 1999. PAUP*. Phylogenetic analysis using parsimony (*and other methods). Version 4b2. Sinauer, Sunderland, Mass.

    Swofford, D. L., G. J. Olsen, P. J. Waddell, and D. M. Hillis. 1996. Phylogenetic inference. Pp. 407–514 in D. M. Hillis, C. Moritz, and B. K. Mable, ed. Molecular systematics. Sinauer, Sunderland, Mass.

    Waddell, P. J., Y. Cao, J. Hauf, and M. Hasegawa. 1999. Using novel phylogenetic methods to evaluate mammalian mtDNA, including amino acid-invariant sites-log det plus site stripping, to detect internal conflicts in the data, with special reference to the positions of hedgehog, armadillo, and elephant. Syst. Biol. 48:31–53.[ISI][Medline]

    Waddell, P. J., N. Okada, and M. Hasegawa. 1999. Towards resolving the interordinal relationships of placental mammals. Syst. Biol. 48:1–5.[ISI][Medline]

    Waits, L. 1996. A comprehensive molecular study of the evolution and genetic variation of bears. Ph.D. thesis, University of Utah, Salt Lake City.

    Xu, S., A. Janke, and U. Arnason. 1996. The complete mitochondrial DNA sequence of the greater Indian rhinoceros, Rhinoceros unicornis, and the phylogenetic relationship among Carnivora, Perrisodactyla and Artiodactyla (+Cetacea). Mol. Biol. Evol. 13:1167–1173.[Abstract]

    Xu, X., and U. Arnason. 1994. The complete mitochondrial DNA sequence of the horse, Equus cabellus: extensive heteroplasmy of the control region. Gene 148:357–362.

    ———. 1996. A complete sequence of the mitochondrial genome of the western lowland gorilla. Mol. Biol. Evol. 13:691–698.[Abstract]

    ———. 1997. The complete mitochondrial DNA sequence of the white rhinoceros, Ceratotherium simum, and comparison with the mtDNA sequence of the Indian rhinoceros, Rhinoceros unicornis. Mol. Phylogenet. Evol. 7:189–194.

    Yang, Z. 1996. Among site variation and its impact on phylogenetic analysis. Trends Ecol. Evol. 11(9):367–372.

    Zardoya, R., and A. Meyer. 1996. Phylogenetic performance of mitochondrial protein-coding genes in resolving relationships among vertebrates. Mol. Biol. Evol. 13:933–942.[Abstract/Free Full Text]

Accepted for publication October 12, 1999.