Estimating metazoan divergence times with a molecular clock

Materials and Methods

Cloning. Total RNA from 17 taxa was prepared from live animals by using a one-step TRIzol method (GIBCO/BRL) or RNAzol (Leedo Medical Laboratories, Houston). Taxa were purchased from Marine Biological Laboratory (Woods Hole, MA; Nucula proxima, Stylochus sp., Obelia sp., and Metridium senile), Gulf Specimen Aquarium and Marine Biological Supply (Panacea, FL; Encope michenlini, Eucidaris tribuloides, and Modiolus americanus), or Charles Hollahan (Santa Barbara, CA, Dendraster excentricus, Strongylocentrotus purpuratus, Mytilus edulis, and Mytilus califorianus). Saccoglossus kowalevskii clones, Monosiga brevicollis cDNA, Antedon mediterrania cDNA, Asterina miniata cDNA, and Priapulis caudatis animals were kind gifts of John Gerhart (Harvard University, Cambridge, MA), Nicole King and Sean Carroll (University of Wisconsin, Madison), Ina Arnone (Stazione Zoologica Anton Dohrn, Napoli, Italy), Veronica Hinman and Eric Davidson (California Institute of Technology, Pasadena, CA), and Graham Budd (University of Uppsala, Uppsala, Sweden), respectively. Ptychodera flava, Chaetopterus sp., Enallagma aspersum, Lestes congener, and Clypeatula cooperensis were already in the collections of K.J.P. and M.A.M. cDNA synthesis was performed with RETROscript (Ambion, Austin, TX) following the manufacturer's instructions by using 1–2 μg of total RNA as indicated above.

Partial fragments of seven nuclear-encoded genes were PCR amplified, cloned, and sequenced by using standard techniques: aldolase (200 aa), triosephosphate isomerase (217 aa), phosphofructokinase (175 aa), methionine adenosyltransferase (348 aa), elongation factor 1-α (418 aa), ATP synthase β chain (430 aa), and catalase (264 aa) (we were unable to amplify all seven from the choanoflagellate, M. califorianus, and S. kowalevskii). These genes were chosen because they had previously been shown to support the monophyly of Ecdysozoa or did not significantly support an alternative arrangement (9, 22) and/or had shown potential clock-like behavior (23). We stress that no molecule or region of a molecule was excluded from the analysis, and the successful amplification and cloning of only these seven (of 12 tested) proved tractable from this diversity of taxa using standard techniques. Gene-specific primers (sequences available on request) and 1 or 10 μl of cDNA plus the TaqPlus Precision PCR system (Stratagene) were mixed and used in touchdown style PCR. PCR fragments of the predicted sizes were excised, purified (Qiagen, Valencia, CA), ligated at 16°C overnight into the pGEM-T-Easy vector according to the manufacturer's instructions (Promega), and electroporated into DH10B cells. Clones containing the correct insert size were sequenced on an ABI373 model sequencer. Sequences were edited, translated, and aligned by using macvector, Ver. 7.0 (Genetics Computer Group).

Phylogenetic Analyses. Dipteran, vertebrate, and plant [Arabidopsis (mustard weed) and Oryza (rice)] sequences were searched by using blast, all significant hits were downloaded, and the inferred amino acid sequences of each gene were analyzed. The topology of each individual gene as deduced by neighbor-joining suggests that each is a case of “many-to-many orthologues” (21). The 50-gene data set was compiled from previous studies (9, 10); the sea urchin genes for this data set were acquired from the Sea Urchin Genome Project web site (http://sugp.caltech.edu). Distance methods used mega, Ver. 2.1, with pairwise deletion (24), and both the Poisson correction and γ distance models [the parameter scaling molecular substitution rate heterogeneity, α, ranged from ∞ (= Poisson distributed) to 0.28]; maximum likelihood (ML) used quartet puzzle, Ver. 5.0 (25) or paml (26). ML analyses used the Jones et al. (27) matrix of amino acid substitution, allowing the analysis to estimate the parameter for substitution rate heterogeneity; all amino acid substitution models gave effectively the same tree (analyses not shown). Relative rates tests used the output of quartet puzzle. Bootstrap values were derived by using 1,000 replications, and 1,000 puzzling steps were performed. Analyses of covariance were performed by using sas, Ver. 6.1 (SAS Institute, Cary, NC).

Divergence Estimates. Date estimates for uncalibrated nodes in phylogenies were derived by using r8s, Ver. 1.5 (M. J. Sanderson, http://ginger.ucdavis.edu/r8s). This software uses multiple calibration points to derive estimates of uncalibrated nodes by using various algorithms. All algorithms for estimating divergence times gave very similar results, and so we report only those derived by the Langley–Fitch likelihood method. Confidence intervals for divergence dates are based on the curvature of the likelihood surface (6) as implemented in r8s.

Results

Rate Heterogeneity Between Vertebrates and Insects. To first ask whether a significant molecular rate difference exists between vertebrates and dipterans, as suggested by comparative genomics (21), we assembled a data set consisting of the concatenation of 50 different nuclear-encoded protein sequences (7,613 aa) taken from the analyses of Wang et al. (9) and Nei et al. (10). The ML analysis of this data set is shown in Fig. 1. As expected, the correct topology is realized and strongly supported. We applied standard relative rates tests to examine whether rate differences exist between the vertebrate and dipteran lineages, using Arabidopsis as the outgroup. All pairwise tests for differences between the two lineages were significant (all P < 0.005), indicating strong rate heterogeneity between the lineages leading to vertebrates and dipterans. However, all pairwise relative rates tests comparing taxa within these lineages (e.g., comparing fish and mouse with Drosophila as the outgroup) were not significant (all P > 0.05), suggesting no rate heterogeneity within the two lineages.

If the bird/mammal divergence (300 Ma) is used to calibrate the tree, we find that Osteichthyes arose 445 Ma, an estimate congruent with the paleontological record (20) (Fig. 1 A Upper). In addition, we find that LCB arose 1,015 Ma, which is close to previous estimates (e.g., refs. 9 and 10). Nonetheless, we find that dipterans, which are near the apex of the insect tree (19, 28), arose 565 Ma, almost 200 million years (Myr) before the first appearance of insects in the fossil record (19). If instead the dipteran divergence is used to calibrate the tree, then the vertebrate divergences are far too shallow, with amniotes originating in the Early Cretaceous (125 Ma) and Osteichthyes originating during the Early Jurassic (185 Ma) (Fig. 1 A Lower).

Vertebrates Significantly Slowed Their Rate of Molecular Evolution. Because the presence of dipterans in the precambrian and the absence of bony fish in the Paleozoic are both untenable, these taxa must differ substantially in the rate of molecular evolution of the included sequences, as suggested by the relative rates tests. Although it is possible that insects increased their rate of molecular evolution with respect to vertebrates, possibly because of their faster generation time (20), it seems as likely that vertebrates slowed their rate of evolution with respect to dipterans.

To distinguish between these two alternatives, we analyzed seven of the 50 sequences discussed above from 14 invertebrate taxa from all three major bilaterian groups chosen specifically to maximize the number of calibration points across the Phanerozoic (Lowest Ordovician through Miocene) within a known phylogeny: Deuterostomia (five echinoderm calibrations), Ecdysozoa (two additional insect calibrations), and Spiralia (three bivalve calibrations); these calibration points are numbered 1–11 on Figs. 2 and 3 and are listed in Table 1. The seven different sequences were concatenated (2,052 aa) and analyzed with ML and minimum evolution; both analyses were accurate and, for most nodes, precise (Figs. 2 and 3). Furthermore, analyses of this data set gave qualitatively similar results when compared to the 50-gene data set (Figs. 1 and 2).

The regression analysis of calibration dates to distance derived from the ML analysis is also shown in Fig. 2. The dipteran divergence (node 10) is not an outlier from the regression line for the invertebrate calibration points (nodes 1–11); analysis of covariance among the lines generated from echinoderms (red), insects (blue), or spiralians (green) shows that neither the slope (F_3,5 = 1.35, P > 0.35) nor the elevation (F_2,7 = 3.10, P > 0.10) is significantly different among these taxa. However, analysis of covariance for the regressions of vertebrate vs. invertebrate calibration points shows that the vertebrate regression is displaced significantly above the invertebrate regression (F_1,8 = 35.81, P < 0.0001) (Fig. 2). This result demonstrates that rather than evolutionary rates increasing in insects (21), molecular evolutionary rates significantly decreased in vertebrates before the origin of crown-group Osteichthyes.

Molecular Clock Estimates of Metazoan Divergence Times. To estimate the origination date of Bilateria, as well as several other invertebrate divergences, we used the r8s software package to analyze the concatenated seven-gene data set for 23 in-group taxa (Fig. 3). Using the option in r8s that the calibration points are fixed, we estimate that LCB (node F, Fig. 3) evolved between 573 and 656 Ma, depending on the specified value of the rate heterogeneity parameter (Table 1). However, because the calibration points derived from the fossil record are paleontological minima (i.e., the first occurrence of a recognizable member of a total group), the estimates derived from these points must be minima as well. Estimating maxima for divergences is difficult (3). Nonetheless, if we specify in r8s that the calibration points are variable, then estimates scale linearly with the error; e.g., if a 10% error is associated with the calibration points, then all estimated dates are increased 7–12%. To ask whether the size of the data set changes the estimate for LCB, we added the sea urchin S. purpuratus to the 50-gene data set and reanalyzed the data without the vertebrates. Calibrating the resulting ML tree with the dipteran divergence gives an estimate of 541 Ma for the last common ancestor of S. purpuratus and Diptera, which is equivalent to LCB.

Finally, to address which γ parameter estimate (Table 1) might be more accurate, we had r8s estimate the origin of crown-group echinoderms (node 5) using both α = ∞ and α = 0.28. The predicted value based on the fossil record of echinoderms (Fig. 4) is between 485 (the minimum value based on the first occurrence of crinoids, Table 1) and 525 Ma (the maximum value that corresponds to the first occurrence of echinoderm skeletal material in the fossil record; ref. 29). Using α = ∞, r8s estimated the age of crown-group echinoderms at 508 ± 12 Ma, whereas with α = 0.28, r8s estimated the age at 527 ± 12 Ma. Thus, the simpler model of molecular evolution gives an estimate more consistent with the fossil record.

Discussion

Rate Heterogeneity Between Vertebrates and Invertebrates. Our data suggest that, inconsistent with most molecular clock estimates but consistent with paleontological predictions (e.g., refs. 13 and 14), bilaterians do not have a significant precambrian evolutionary history. The deep precambrian estimates for LCB derived from analyses that use nuclear protein-encoding genes calibrated to the vertebrate fossil record (e.g., refs. 9 and 10) are clearly artifacts associated with the significant rate reduction in the molecular evolution of the vertebrate genome.

Martin and Palumbi (30) argued that much of the rate heterogeneity that exists between species can be accounted for by differences in metabolic rates. Although no single factor can fully account for rate variation (30), we note that a difference in metabolic rates is unlikely to be the primary explanation in this case because the teleost fish Danio rerio is evolving at a similar rate to that of the two endothermic amniotes, and the molecular clock estimate of the origin of Osteichthyes is concordant with the vertebrate fossil record (Fig. 1 A). Given the metabolic rate difference between amniotes and fishes (31), we would expect a greater disparity between the clock estimate and the vertebrate fossil record if this were a significant driving factor. One possible explanation for the rate heterogeneity seen between vertebrates and invertebrates is the duplication of the vertebrate genome (12), which occurred sometime between the last common ancestor of cephalochordates and vertebrates and the last common ancestor of Osteichthyes (32). Despite the fact that gene duplication events are thought to increase rather than decrease the rate of molecular evolution (33), a genome duplication event would (at least initially) increase the number of interactors for each protein, potentially slowing the rate of molecular evolution across the entire genome (34).

Phylogenetic Considerations. Our protein tree (Fig. 3) gives the correct topology, where known (Table 1), and finds support for clades such as Ambulacraria (Echinodermata + Hemichordata), Spiralia, and Ecdysozoa, as well as Protostomia. This is now the third independent molecular data set supporting the monophyly of these clades, because all four are found with 18S rDNA analyses and Hox gene duplications, in addition to many newly elucidated characters (reviewed in ref. 35). We are not able to recover a monophyletic Deuterostomia when vertebrates are included into the analysis (Fig. 2), possibly because of the relatively few number of genes analyzed (36) and the pronounced rate heterogeneity detected with the vertebrate sequences (Fig. 2). Nonetheless, an analysis of the 50-gene data set, which includes both the sea urchin and the three vertebrates, results in a monophyletic but weakly supported Deuterostomia (not shown).

Tempo and Mode of Early Animal Evolution. Although the use of molecular clocks to infer divergence times is fraught with difficulties (3, 4), this analysis fulfills the suggested requirements of Shaul and Graur (37) for a molecular clock analysis: (i) the use of multiple primary calibration points; (ii) the accommodation of rate variation; and (iii) the calculation of confidence intervals associated with the estimates. Interestingly, both our analysis and the analyses of Aris-Brosou and Yang (11, 12) conclude that LCB evolved ≈570 Ma and split from cnidarians somewhere between 600 and 630 Ma (Fig. 4). Moreover, both analyses agree that the last common ancestor of protostomes evolved ≈550 Ma, and both analyses agree that phylum-level splits within Spiralia and Deuterostomia occurred 520–530 Ma. The congruence between these two clock studies and the fossil record is striking, although possibly not surprising given that both studies use multiple calibration points scattered across both phylogeny and time and account for rate heterogeneity either by removing vertebrates from the analysis, as was done here, or by using a Bayesian approach to account for rate change across lineages (11, 12). Because of this congruence, the Cambrian explosion must reflect, at least in part, the diversification of bilaterian phyla.

We extrapolate the origin of total-group Bilateria to be ≈615 Ma, almost 45 Myr before the appearance and rapid diversification of the crown group and ≈60 Myr before their first unequivocal appearance in the fossil record. This cryptic precambrian history suggests that these stem-group bilaterians must have been “micrometazoans,” because, although benthic (38, 39), they were seemingly incapable of leaving trace fossils (see also ref. 40). The bilaterian trace fossil record would commence only after the invention of pattern formation mechanisms that potentiated the evolution of larger body size in multiple animal clades near the end of Neoproterozoic (41, 42). Whether the origin of bilaterians or any other metazoan group can be triggered by environmental perturbations such as “snowball Earth” (43, 44) remains highly speculative at the moment, given the uncertainty about the exact number and ages of Neoproterozoic glaciations (45). Nonetheless, if the “Marinoan” glaciation interval is younger than 590 Ma (46), then an increase in bilaterian body size might have been facilitated by the greater productivity of the marine ecosystem after the glacial melt (47, 48). This increase in both body size and planktonic productivity would then allow for the evolution of broadcast spawning and external fertilization and, by ≈525 Ma, planktotrophic development. The absence of precambrian planktonic metazoans is consistent with our Early Cambrian estimate for the origin of the last common ancestor of living cnidarians, a population of animals whose life cycle lacked a medusa stage and thus was entirely benthic (49). The development of this new planktonic food web (47), coupled with the evolution of a dispersal stage and the reappearance of exposed continental shelf, may have provided the environmental stimuli necessary for the rapid evolution of disparate bilaterian body plans and ultimately the Cambrian explosion itself.