Bayesian Models of Episodic Evolution Support a Late Precambrian Explosive Diversification of the Metazoa

Stéphane Aris-Brosou1 and Ziheng Yang

Department of Biology, University College London, England

Correspondence: E-mail: z.yang{at}ucl.ac.uk.


    Abstract
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results and Discussion
 Acknowledgements
 Literature Cited
 
Multicellular animals, or Metazoa, appear in the fossil records between 575 and 509 million years ago (MYA). At odds with paleontological evidence, molecular estimates of basal metazoan divergences have been consistently older than 700 MYA. However, those date estimates were based on the molecular clock hypothesis, which is almost always violated. To relax this hypothesis, we have implemented a Bayesian approach to describe the change of evolutionary rate over time. Analysis of 22 genes from the nuclear and the mitochondrial genomes under the molecular clock assumption produced old date estimates, similar to those from previous studies. However, by allowing rates to vary in time and by taking small species-sampling fractions into account, we obtained much younger estimates, broadly consistent with the fossil records. In particular, the date of protostome–deuterostome divergence was on average 582 ± 112 MYA. These results were found to be robust to specification of the model of rate change. The clock assumption thus had a dramatic effect on date estimation. However, our results appeared sensitive to the prior model of cladogenesis, although the oldest estimates (791 ± 246 MYA) were obtained under a suboptimal model. Bayes posterior estimates of evolutionary rates indicated at least one major burst of molecular evolution at the end of the Precambrian when protostomes and deuterostomes diverged. We stress the importance of assumptions about rates on date estimation and suggest that the large discrepancies between the molecular and fossil dates of metazoan divergences might partly be due to biases in molecular date estimation.

Key Words: Divergence dates • Metazoa • molecular clock • Bayes


    Introduction
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results and Discussion
 Acknowledgements
 Literature Cited
 
The sudden appearance of numerous metazoan phyla in the early Cambrian (Knoll and Carroll 1999) suggests an explosive morphological diversification, but the timing of the divergence of the Metazoa is controversial (Valentine, Jablonski, and Erwin 1999). Paleontological evidence indicates this divergence occurred approximately 600 MYA (Valentine, Jablonski, and Erwin 1999; Conway Morris 2000), whereas molecular studies support much earlier dates (Wray, Levinton, and Shapiro 1996; Feng, Cho, and Doolittle 1997; Bromham et al. 1998; Wang, Kumar, and Hedges 1999). One hypothesis is that animals were small and unlikely to fossilize before the Cambrian explosion (Cooper and Fortey 1998), which is supported by the recent find of a crustacean in early Cambrian strata (511 MYA) (Siveter, Williams, and Waloszek 2001). However, in the absence of uncontroversial animal fossils older than the Ediacaran biota (Brasier 1998; Budd and Jensen 2000; Rasmussen et al. 2002), deep divergence dates may only be reasonably approached by indirect comparative studies.

Molecular date estimates vary widely among studies (Knoll and Carroll 1999), which indicates a possible problem when averaging results over genes evolving at different rates. However, all studies to date suggest that the basal divergence between protostomes and deuterostomes occurred more than 700 MYA (Bromham et al. 1998). This early-origin hypothesis is supported by the analysis of a large number of genes (Wray, Levinton, and Shapiro 1996; Feng, Cho, and Doolittle 1997; Gu 1998; Wang, Kumar, and Hedges 1999). However, these studies are all based on the molecular clock hypothesis (Zuckerkandl and Pauling 1965), with lineages violating the clock removed through sequential relative rate tests (Wu and Li 1985). The power of such tests was questioned (Bromham et al. 2000), and it has been suggested that violation of the clock could have drastic effects on date estimation (Ayala, Rzhetsky, and Ayala 1998; Bromham and Hendy 2000; Yoder and Yang 2000).

In this paper, we take a Bayes approach (Thorne, Kishino, and Painter 1998) to examine the effect of models of rate evolution on estimates of metazoan divergence dates. By analyzing 11 nuclear genes and 11 mitochondrial genes, we show that the Bayes analysis under the molecular clock gives date estimates comparable to those reported in previous molecular studies (Wray, Levinton, and Shapiro 1996; Feng, Cho, and Doolittle 1997; Gu 1998; Wang, Kumar, and Hedges 1999). However, when the clock assumption was relaxed, we obtained younger date estimates, largely consistent with the fossil records.


    Materials and Methods
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results and Discussion
 Acknowledgements
 Literature Cited
 
The Bayesian Framework
Divergence times (T) and evolutionary rates (R) are assigned prior distributions. Their joint posterior probability given data X is given by p(R,T | X) = p(X | R,T) p(R | T) p(T)/p(X) (Thorne, Kishino, and Painter 1998). The tree topology is assumed known and fixed (Nielsen 1995 [see Supplementary Material online]). The likelihood p(X | R,T) was computed assuming the HKY85 + {Gamma} nucleotide substitution model with eight rate classes (Hasegawa, Kishino, and Yano 1985; Yang 1994). The transition-transversion rate ratio {kappa} and the among-site rate variation parameter {alpha} were integrated out over uniform priors U(0, 1000) and U(0, 100), respectively. Base frequencies were estimated using their empirical frequencies in the data and fixed in the analysis.

Prior Models of Rate Change
Rate change is accommodated by a stochastic process in which the rate of a descendent branch follows a statistical distribution centered on the rate of the ancestral branch (Thorne, Kishino, and Painter 1998). The variance of the distribution is given by s2 = {sigma}2{Delta}t, where {Delta}t is the time duration of the branch. Thus rates evolve in an autocorrelated manner from ancestor to descendent, and branches far apart on the phylogeny are likely to have different rates. Parameter {sigma}2 determines how variable the rates are: a small {sigma}2 means a clock-like tree, whereas a large {sigma}2 means highly variable rates (fig. 1). Here, we explore two statistical distributions to model rate change over branches: the exponential distribution (EXP) and the Ornstein-Uhlenbeck process (OUP), which is a stationary Gaussian process (Aris-Brosou and Yang 2002). Because of the stationarity property, no systematic trend in the rate, either upward or downward, is assumed. EXP has no hyperparameters since the mean of the distribution is fixed at the rate of the ancestral branch. OUP has two hyperparameters: {sigma}2, which controls rate variation among branches, and ß, which is a friction term. Both {sigma}2 and ß were integrated out, assuming vague prior distributions: a gamma distribution with mean of 15 and variance of 25 for {sigma}2 and a log-normal distribution with mean log(0.5) and variance of 0.75 for ß. These values were chosen according to results from an Empirical Bayes study of the 18S rRNA gene (Aris-Brosou and Yang 2002), although optimal values may differ among genes.



View larger version (12K):
[in this window]
[in a new window]
 
FIG. 1. Effect of the variance parameter {sigma}2 in the model of rate change on the distribution of rates for branches. The rate ri of a branch is drawn from a distribution that is centered around the rate rA for the ancestral branch. When {sigma}2 is small, ri is close to rA, so that the rates of the different branches are very similar and the tree is clock-like (A and B). Conversely, a large {sigma}2 allows more variation in the branch-specific rates (C and D)

 
We used the likelihood ratio test (LRT) and the posterior Bayes factor (PBF) to test the molecular clock assumption. In the LRT, twice the log-likelihood difference between the clock and nonclock models is compared with a {chi}2 distribution with n 2 degrees of freedom for a tree of n species. PBF (Aitkin 1991) is the ratio of the likelihood L averaged over the posterior distribution under each model: PBF1,2 = E{theta}|X[L1] / E{theta}|X[L2]. PBF was also used to compare models of rate change. This measures the weight of evidence in the sample in favor of model 1 against model 2, with values greater than 20, 100, and 1,000, meaning strong, very strong, and overwhelming evidence in favor of model 1 (Aitkin 1991). The logarithm of PBF is used in the rest of the text. We note that PBF is highly controversial but used it because it is easy to calculate from the Markov chain Monte Carlo (MCMC) runs. For the data sets analyzed in this study, the LRT of the clock and PBF led to the same conclusions.

Prior Model of Speciation
The Bayes approach also requires modeling speciation events to specify the prior distribution of divergence times. This distribution was specified by a generalized birth-death process, allowing for species sampling (Yang and Rannala 1997). The birth-death process tends to generate trees with long internal branches. Taking incomplete species sampling into account generates more realistic trees (Yang and Rannala 1997). To avoid overparameterization, the time for the root is fixed at 1 (Rannala 2002; Aris-Brosou and Yang 2002), so that divergence times at other internal nodes are relative to that at the root. Hyperparameters of the model of speciation were integrated out: birth and death rates follow uniform prior distributions on (0, 15) and (0, 5), respectively. The sampling fraction {rho} follows a uniform prior distribution on (0, 0.001) unless otherwise stated. After running the MCMC, a calibration date was used to rescale the relative divergence times into absolute divergence times measured in MYA.

Approximation of the Posterior Distribution
For each gene, the marginal posterior distributions p(T | X) and p(R | X) were approximated using an MCMC algorithm based on the Metropolis-Hastings sampler. The algorithm was described in Aris-Brosou and Yang (2002) and is similar to that of Thorne, Kishino, and Painter (1998). Let {theta} be the parameters of the model (divergence times, branch-specific rates, all the hyperparameters of the prior distributions for rates and times, {kappa} and {alpha}). At each step of the Markov chain, a new state {theta}* is proposed, differing from the current state {theta} by the update of one single parameter. We used normal proposal densities centered on the current parameter values. The proposed state is then accepted with probability min{1, p({theta}* | X) / p({theta} | X)}. For the data sets analyzed in table 1, the 50,000 first steps of the chain were discarded (burn-in), and each chain was then sampled every 500 steps until 10,000 samples were collected. Convergence was checked by running four additional shorter chains. For the five runs, we analyzed time-series outputs for each parameter and checked consistency of the estimates across the different runs. Inferences were based on the median of the marginal distributions of parameters drawn from the longest run. The computation involved was intensive and took several months on a 30-processor Pentium III Beowulf cluster. The computer program implementing the models described is available at http://statgen.ncsu.edu/stephane/. A simulation study is conducted to examine the performance of the method (see Results and Discussion).


View this table:
[in this window]
[in a new window]
 
Table 1 Statistics of the Genes Analyzed and Divergence Dates of Two Major Clades (MYA).

 
DNA Sequences
We analyzed 11 nuclear genes and 11 mitochondrial genes, with an average number of 26 taxa and 1,388 nucleotides per gene (table 1). Nucleotide sequences were obtained from GenBank for 11 nuclear genes (18S rRNA, actin, {alpha}-tubulin, ß-tubulin, calreticulin, catalase, elongation factor 1 [EF-1], histone H1, heat shock protein 70 [Hsp70], protein kinase c [Pkc], troponin c) and 11 mitochondrial genes (cytochrome c oxidase [Cox] subunits I, II and III, cytochrome b [Cyt b], nicotinamide adenine dinucleotide dehydrogenase [ND] subunits 1 to 6 and 4L). For GenBank accession numbers of the sequences, see online Supplementary Material. These genes were chosen because (1) they are represented extensively across the Metazoa, (2) they cover both protostome-deuterostome and echinoderm-chordate splits, and (3) at least one fossil calibration point at an age at least 300 MYA is available. Alignments were performed with ClustalW version 1.8 (Thompson, Higgins, and Gibson 1994) and adjusted manually. For protein-coding genes, all three codon positions were included in the analyses, with among-site rate variation accounted for using a discrete gamma model with eight rate classes (Yang 1994).

Calibration Points
The phylogenetic tree for each gene is rooted using either a land plant (Arabidopsis), a fern (Polypodium for the 18S rRNA gene), or a fungus (Schizosaccharomyces for the troponin c gene). To reflect the most basal split (Parazoa-Eumetazoa), a diploblastic animal (Cnidaria) is included in the analysis whenever possible. To reduce errors associated with calibration points, only fossil-based dates were used, as in Bromham et al. (1998) (in): Collembola-Pterygota, 390 MYA; Aranaea-Scorpionida, 405 MYA; Coelacanth-Dipnoi/Tetrapoda, 418 MYA; Osteichthyes-Dipnoi/Tetrapoda, 428 MYA; Asteroidea-Echinoidea, 500 MYA; Agnata-Gnathostoma, 510 MYA; Arachnida-Merostomata, 520 MYA; Cephalochordata-Chordata, 530 MYA. Each gene has between one and eight calibration points (table 1). When several calibration points were used, the median of the estimated divergence times was used as the final estimate.

We focused on two key evolutionary transitions: the protostome-deuterostome (PD) divergence, which marks the appearance of "higher Metazoa" (Eumetazoa), and the echinoderm-chordate (EC) divergence, as it predates the origin of the vertebrates.


    Results and Discussion
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results and Discussion
 Acknowledgements
 Literature Cited
 
Computer Simulation to Examine the Performance of the Bayes Algorithm
We conducted computer simulations to examine the performance of the Bayes MCMC algorithm (fig. 2). One hundred replicate data sets were generated for three trees (fig. 2AC), each of which includes eight in-group taxa and one out-group taxon. Sequences, each of 1,000 nucleotides, were generated under the JC69 substitution model (Jukes and Cantor 1969) using the program Evolver from the PAML package (Yang 1997). Data sets were then analyzed under the same substitution model, assuming either a Bayesian molecular clock (with all branches having the same rate) or the exponential model of rate change. For each MCMC analysis, the first 20,000 samples of the chain were discarded as burn-in, after which 500 samples were collected for inference, sampling once every 1,000 iterations. Convergence was checked by time series plots for divergence times sampled from the posterior distribution. The results are presented in figure 2.



View larger version (19K):
[in this window]
[in a new window]
 
FIG. 2. Performance of the Bayesian molecular clock model and the exponential model of rate change in divergence date estimation. Data sets were simulated using topologies 1 to 3 with branch lengths indicated by the scale bars (A–C) and were analyzed under the Bayesian molecular clock (D–F) model and under the exponential model of rate change (G–I). The "true" divergence times are indicated by vertical gray broken lines. The clock holds in topologies 1 and 2, while in topology 3, a terminal branch shows a twofold rate acceleration. The density plots (D–I) represent the posterior distributions of the node times on a relative scale (0 = present; 1 = root time)

 
Topology 1 conforms to a perfect molecular clock, and both the Bayesian clock model and the EXP model of rate change performed well. Topology 2 also conforms to the clock, but the shape of the tree reflected in the relative divergence times is very different from that expected under the birth-death process model with species sampling. For this tree, the Bayes analysis assuming the clock performed well, but the exponential model of rate change tended to overestimate divergence times and underestimate substitution rates for terminal branches (fig. 2H). In topology 3, the molecular clock is violated, with a terminal branch having a rate twice as high as all other branches. When the clock was incorrectly assumed (fig. 2F), divergence dates t13 and t15 are seriously overestimated. Note that in the true tree topology, there are only two distinct divergence dates (i.e., t14 = t15, and the four most recent nodes have the same age). In contrast, the exponential model of rate change performed much better. The four recent nodes had dates all close to the truth, while t14 and t15 were similar as well although slighted biased upwards (fig. 2I). Note that posterior distributions under the EXP model are wider than under the clock model for all three trees, because the rate-change model involves more parameters and incorporates more variability in the data.

Divergence Dates Under the Bayesian Molecular Clock
The molecular clock hypothesis is tested using two approaches: the LRT and the posterior Bayes factor (PBF). For every gene studied, the clock hypothesis was strongly rejected by the LRT (P < 0.001 for all the genes) and the PBF (table 1). Date estimates of the PD divergence under the clock are above 700 MYA for most genes (table 1), with a median (1st to 3rd quartiles) of 1,090 MYA (1,411 to 841) (i.e., before the Vendian [fig. 3A]). These results are in agreement with previous molecular studies (Wray, Levinton, and Shapiro 1996; Feng, Cho, and Doolittle 1997; Gu 1998; Wang, Kumar, and Hedges 1999). Different genes produced substantially different estimates for the PD divergence, ranging from approximately 500 MYA for calreticulin to approximately 1,990 MYA for cox1 (fig. 3A).



View larger version (27K):
[in this window]
[in a new window]
 
FIG. 3. Posterior distributions of the divergence time between protostomes and deuterostomes. Eleven nuclear genes (solid lines) and 11 mitochondrial genes (broken lines) were analyzed under three models of rate change: (A) the Bayesian molecular clock, (B) the exponential distribution, (C) the Ornstein-Uhlenbeck process

 
Divergence Dates Under Models of Rate Change
The molecular clock assumption was relaxed by modeling rate change over time either by the exponential model (EXP) or by the Ornstein-Uhlenbeck process (OUP). Both processes are stationary, and neither assumes a trend in the evolutionary rate. In particular, unlike the model of Bromham and Hendy (2000), they do not posit fast early rates. Table 1 shows that EXP and OUP gave similar estimates, and the PD divergence is dated at 579 MYA (607 to 554) and 582 MYA (617 to 556), respectively (see also figure 3B and C). Those estimates are broadly consistent with paleontological data (Valentine, Jablonski, and Erwin 1999). The estimates for the PD divergence are close to the minimum estimate of 588 MYA obtained by Bromham and Hendy (2000), assuming elevated rates around the root. Note that the only difference between the clock and nonclock analyses in table 1 lies in the assumption about substitution rates while all other assumptions are the same. The effect of relaxing the molecular clock is thus dramatic.

Effect of the Prior Assumptions on Date Estimates
Before drawing a conclusion, we examine the effects of several factors. First, we assess the influence of the prior by running the MCMC with no data, that is, by fixing f(X|T,R) = 1 in the MCMC, so that the Markov chain converges to the prior distribution. For instance, with 36 species and the four calibration points used for the mitochondrial genes, prior divergence times under the clock and under models of rate change were very similar, with tclock(PD) = 526 (95% credible set (CS): 482 to 611) and tOUP(PD) = 526 (95% CS: 482 to 609). Although these dates are much younger than the dates estimated under the clock, they are close to the estimates under models of rate change. The results suggest that the prior for divergence times may have some influence on date estimates under models of rate change. This appears to be the same pattern as seen in the simulation study (fig. 2H).

Second, the exponential model used is crude and unrealistic, whereas the Ornstein-Uhlenbeck process is more complex. Table 1 shows that the latter model fitted the data better than the former for most genes (median log PBFOE = 16). Both models clearly outperformed the clock model (e.g., log PBFOC >> 10 [table 1]) and, more importantly, gave consistent estimates of divergence dates for the PD and the EC splits (table 1 and fig. 3B and C). Therefore, our date estimates for these two splits appear robust to the specification of the model of rate change and may not be improved by more complex models.

Third, the sampling fraction in the birth-death process with species sampling is known to affect the shape of the tree and thus may influence divergence date estimation. Our analysis (table 1 and fig. 3) assumes that the sampling fraction has a uniform prior distribution U(0, {rho}up) with the upper bound {rho}up = 0.001. Larger {rho}up values did not greatly affect estimates under the clock, but gave older estimates under models of rate change: if we assume that {rho}up = 0.5, the PD split was estimated at approximately 791 ± 246 MYA. However, such a large {rho}up was strongly rejected by the posterior Bayes factor (log PBF0.001/0.5 = 24.11). On the other hand, assuming a {rho}up smaller than that used in table 1 (0.0001) gave an estimate of approximately 561 ± 144 MYA for the PD split, which was favored by the Bayes factor (log PBF0.001/0.0001 = –6.27). We note that this latter estimate is very close to the one found for a sampling fraction an order of magnitude larger ({rho}up = D0.001). Thus the date estimates are stable over a reasonable range of {rho}up values.

Lastly, we did not account for uncertainty in the tree topology, although controversy exists concerning the metazoan phylogeny (Knoll and Carroll 1999). We note however that plausible topologies gave similar speciation date estimates in previous studies (Yoder and Yang 2000).

Episodic Molecular Evolution
It is of interest to examine whether the Cambrian explosion, as recorded by the fossils, has been preceded by a burst of molecular evolution, as suggested recently (Bromham and Hendy 2000). Figure 4 summarizes the estimates of relative rates against time from the exponential model of rate change. We define elevated relative rates as those greater than the 95th percentile of the distribution of relative rates over branches and over the sampled genes (rightmost panel of figure 4). High relative rates occur mainly between approximately 640 MYA (late Riphean) and approximately 420 MYA (Silurian). The average relative rate is almost twice as large during this period (1.37) than either before (0.71) or after (0.62) it. Elevated rates are mainly for branches before the PD, the EC, and the Agnatha-Gnathostoma (jawless and jawed vertebrates) divergences. The last two are contiguous and may belong to a single long period of elevated rates. Because of our limited sampling around the Parazoa-Eumetazoa split, it is difficult to detect any such burst at that time. It is remarkable that these bursts of evolution have concerned most of the 22 genes analyzed. We consider mutation rate variation to be an unlikely explanation for such evolutionary rate acceleration at the base of the metazoan, and a more sensible alternative is changed selective pressure. The large-scale rate acceleration might correspond to major duplication events, leading to the relaxation of selective constraints and higher evolutionary rates (Pollard and Holland 2000; Miyata and Suga 2001). Other lineages with high estimated rates (fig. 4) belong to the invertebrates, but here high rates are more gene-dependent. Subsequent "bursts" of evolution (<400 MYA) are smaller in magnitude and mainly concern, at least in our restricted species sampling, parasites. Our date estimates suggest that divergences have been explosive and happened in a relatively short period of time, as the PD and EC splits were found to be separated by an average of 63 MY. The environmental elements that could have triggered the Cambrian explosion remain unclear (Knoll and Carroll 1999), but our molecular date estimates are compatible with an origination of the Metazoa at about the Varanger ice age (~620 to 580 MYA), hereby renewing interest into possible refugia during a snowball Earth (Hyde et al. 2000).



View larger version (15K):
[in this window]
[in a new window]
 
FIG. 4. Relative rates of evolution plotted against the estimated divergence dates for the twenty-two genes studied. For each gene the relative rate is calculated by the estimated rate for the branch divided by the average rate for the gene. The x-axis indicates the age of the descendent node of the branch. The two dashed horizontal lines correspond to the median and the 95th percentile of the distribution, shown on the right, of relative rates over branches and genes. Four branches leading to the following divergences are indicated as: Parazoa–Eumetazoa split (x); PD split ({blacksquare}); EC split (+); split basal to vertebrates (•). A schematic geological scale is also given (topmost)

 
Consistently with other studies (Ayala, Rzhetsky, and Ayala 1998; Bromham and Hendy 2000; Yoder and Yang 2000), our results highlight the importance of assumptions about evolutionary rates in divergence date estimation. By relaxing the unrealistic molecular clock assumption, we obtained date estimates for Metazoan divergences much closer to fossil data than previous molecular estimates, most of which are based on the clock assumption. A very similar pattern has been reported recently by Adkins, Walton, and Honeycutt (2003), who estimated divergence dates among rodents using a large data set of more than 4,600 aligned nucleotide sequences. Those authors found that divergence dates among rodents estimated under a global molecular clock were compatible with previous molecule-based dates estimated under similar conditions, exceeding fossil-based dates. However, when a relaxed molecular clock was applied, estimated divergence dates were highly compatible with the fossil record.

We note that estimation of ancient divergence dates when the molecular clock is violated is problematic, and our date estimates should be confirmed by analyzing more genes. In particular, our approach of averaging over calibration points and over multiple genes is inferior to a simultaneous analysis of multiple genes under the constraints of multiple fossil calibrations (Thorne, Kishino, and Painter 1998; Thorne and Kishino 2002). We emphasize the sensitivity of date estimation to assumptions about rates, and suggest that development of powerful estimation methodologies and accumulation of more gene sequences will eventually resolve the issue of Metazoan divergences.


    Acknowledgements
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results and Discussion
 Acknowledgements
 Literature Cited
 
We thank Hirohisa Kishino, Jeffrey Thorne, Peter Holland, Lindell Bromham, and Andrew Rambaut for discussions, and we thank two referees for critical comments. This study was supported by a Biotechnological and Biological Sciences Research Council grant to Z.Y.


    Footnotes
 
1 Present address: Bioinformatics Research Center, North Carolina State University, Raleigh. Back

Manolo Gouy, Associate Editor Back


    Literature Cited
 TOP
 Abstract
 Introduction
 Materials and Methods
 Results and Discussion
 Acknowledgements
 Literature Cited
 

    Adkins, R. M., A. H. Walton, and R. L. Honeycutt. 2003. Higher-level systematics of rodents and divergence time estimates based on two congruent nuclear genes. Mol. Phylogenet. Evol. 26:409-420.[CrossRef][ISI][Medline]

    Aitkin, M. 1991. Posterior Bayes factors. J. R. Statist. Soc. B 53:111-143.[ISI]

    Aris-Brosou, S., and Z. Yang. 2002. The effects of models of rate evolution on estimation of divergence dates with a special reference to the metazoan 18S rRNA phylogeny. Syst. Biol. 51:703-714.[CrossRef][ISI][Medline]

    Ayala, F. J., A. Rzhetsky, and F. J. Ayala. 1998. Origin of the metazoan phyla: molecular clocks confirm paleontological estimates. Proc. Natl. Acad. Sci. USA 95:606-611.[Abstract/Free Full Text]

    Brasier, M. D. 1998. From deep time to late arrival. Nature 395:547-548.[CrossRef][ISI]

    Bromham, L. D., and M. D. Hendy. 2000. Can fast early rates reconcile molecular dates with the Cambrian explosion? Proc. R. Soc. Lond. B Biol. Sci. 267:1041-1047.[CrossRef][ISI][Medline]

    Bromham, L., D. Penny, A. Rambaut, and M. D. Hendy. 2000. The power of relative rates tests depends on the data. J. Mol. Evol. 50:296-301.[ISI][Medline]

    Bromham, L., A. Rambaut, R. Fortey, A. Cooper, and D. Penny. 1998. Testing the Cambrian explosion hypothesis by using a molecular dating technique. Proc. Natl. Acad. Sci. USA 95:12386-12389.[Abstract/Free Full Text]

    Budd, G. E., and S. Jensen. 2000. A critical reappraisal of the fossil record of the bilaterian phyla. Biol. Rev. Camb. Philos. Soc. 75:253-295.[CrossRef][ISI][Medline]

    Conway Morris, S. 2000. The Cambrian "explosion": slow-fuse or megatonnage? Proc. Natl. Acad. Sci. USA 97:4426-4429.[Abstract/Free Full Text]

    Cooper, A., and R. Fortey. 1998. Evolutionary explosions and the phylogenetic fuse. Trends Ecol. Evol. 13:151-156.[CrossRef][ISI]

    Feng, D. F., G. Cho, and R. F. Doolittle. 1997. Determining divergence times with a protein clock: update and reevaluation. Proc. Natl. Acad. Sci. USA 94:13028-13033.[Abstract/Free Full Text]

    Gu, X. 1998. Early metazoan divergence was about 830 million years ago. J. Mol. Evol. 47:369-371.[ISI][Medline]

    Hasegawa, M., H. Kishino, and T. Yano. 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22:160-174.[ISI][Medline]

    Hyde, W. T., T. J. Crowley, S. K. Baum, and W. R. Peltier. 2000. Neoproterozoic ‘snowball Earth’ simulations with a coupled climate/ice-sheet model. Nature 405:425-429.[CrossRef][ISI][Medline]

    Jukes, T. H., and C. R. Cantor. 1969. Evolution of protein molecules. Pp. 21–32 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York.

    Knoll, A. H., and S. B. Carroll. 1999. Early animal evolution: emerging views from comparative biology and geology. Science 284:2129-2137.[Abstract/Free Full Text]

    Miyata, T., and H. Suga. 2001. Divergence pattern of animal gene families and relationship with the Cambrian explosion. Bioessays 23:1018-1027.[CrossRef][ISI][Medline]

    Nielsen, C. 1995. Animal evolution: interrelationships of the living phyla. Oxford University Press, Oxford, U.K.

    Pollard, S. L., and P. W. Holland. 2000. Evidence for 14 homeobox gene clusters in human genome ancestry. Curr. Biol. 10:1059-1062.[CrossRef][ISI][Medline]

    Rannala, B. 2002. Identifiability of parameters in MCMC Bayesian inference of phylogeny. Syst. Biol. 51:754-760.[CrossRef][ISI][Medline]

    Rasmussen, B., S. Bengtson, I. R. Fletcher, and N. J. McNaughton. 2002. Discoidal impressions and trace-like fossils more than 1200 million years old. Science 296:1112-1115.[Abstract/Free Full Text]

    Siveter, D. J., M. Williams, and D. Waloszek. 2001. A phosphatocopid crustacean with appendages from the Lower Cambrian. Science 293:479-481.[Abstract/Free Full Text]

    Thompson, J. D., D. G. Higgins, and T. J. Gibson. 1994. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22:4673-4680.[Abstract]

    Thorne, J. L., and H. Kishino. 2002. Divergence time and evolutionary rate estimation with multilocus data. Syst. Biol. 51:689-702.[CrossRef][ISI][Medline]

    Thorne, J. L., H. Kishino, and I. S. Painter. 1998. Estimating the rate of evolution of the rate of molecular evolution. Mol. Biol. Evol. 15:1647-1657.[Abstract/Free Full Text]

    Valentine, J. W., D. Jablonski, and D. H. Erwin. 1999. Fossils, molecules and embryos: new perspectives on the Cambrian explosion. Development 126:851-859.[Abstract/Free Full Text]

    Wang, D. Y., S. Kumar, and S. B. Hedges. 1999. Divergence time estimates for the early history of animal phyla and the origin of plants, animals and fungi. Proc. R. Soc. Lond. B Biol. Sci. 266:163-171.[CrossRef][ISI][Medline]

    Wray, G. A., J. S. Levinton, and L. H. Shapiro. 1996. Molecular evidence for deep Precambrian divergences. Science 274:568-573.[Abstract/Free Full Text]

    Wu, C. I., and W. H. Li. 1985. Evidence for higher rates of nucleotide substitution in rodents than in man. Proc. Natl. Acad. Sci. USA 82:1741-1745.[Abstract]

    Yang, Z. 1994. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J. Mol. Evol. 39:306-314.[ISI][Medline]

    Yang, Z. 1997. PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci. 13:555-556.[Medline]

    Yang, Z., and B. Rannala. 1997. Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method. Mol. Biol. Evol. 14:717-724.[Abstract]

    Yoder, A. D., and Z. Yang. 2000. Estimation of primate speciation dates using local molecular clocks. Mol. Biol. Evol. 17:1081-1090.[Abstract/Free Full Text]

    Zuckerkandl, E., and L. Pauling. 1965. Evolutionary divergence and convergence in proteins. Pp. 97–166 in V. Bryson, and H. J. Vogel, eds. Evolving genes and proteins. Academic Press, New York.

Accepted for publication June 7, 2003.