MBE Advance Access originally published online on June 15, 2006
Molecular Biology and Evolution 2006 23(10):1832-1851; doi:10.1093/molbev/msl039
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Research Article |
Testing the Molecular Clock: Molecular and Paleontological Estimates of Divergence Times in the Echinoidea (Echinodermata)
,1
,1,1,1
* Department of Palaeontology, The Natural History Museum, London, United Kingdom; Bioinformatics Laboratory, The National University of Ireland, Maynooth, County Kildare, Ireland; and Department of Zoology, The Natural History Museum, Cromwell Road, London, United Kingdom
E-mail: a.smith{at}nhm.ac.uk.
 |
Abstract |
|---|
 TOP Abstract IntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
The phylogenetic relationships of 46 echinoids, with representatives from 13 of the 14 ordinal-level clades and about 70% of extant families commonly recognized, have been established from 3 genes (3,226 alignable bases) and 119 morphological characters. Morphological and molecular estimates are similar enough to be considered suboptimal estimates of one another, and the combined data provide a tree that, when calibrated against the fossil record, provides paleontological estimates of divergence times and completeness of their fossil record. The order of branching on the cladogram largely agrees with the stratigraphic order of first occurrences and implies that their fossil record is more than 85% complete at family level and at a resolution of 5-Myr time intervals.
Molecular estimates of divergence times derived from applying both molecular clock and relaxed molecular clock models are concordant with estimates based on the fossil record in up to 70% of cases, with most concordant results obtained using Sanderson's semiparametric penalized likelihood method and a logarithmic-penalty function. There are 3 regions of the tree where molecular and fossil estimates of divergence time consistently disagree. Comparison with results obtained when molecular divergence dates are estimated from the combined (morphology + gene) tree suggests that errors in phylogenetic reconstruction explain only one of these. In another region the error most likely lies with the paleontological estimates because taxa in this region are demonstrated to have a very poor fossil record. In the third case, morphological and paleontological evidence is much stronger, and the topology for this part of the molecular tree differs from that derived from the combined data. Here the cause of the mismatch is unclear but could be methodological, arising from marked inequality of molecular rates. Overall, the level of agreement reached between these different data and methodological approaches leads us to believe that careful application of likelihood and Bayesian methods to molecular data provides realistic divergence time estimates in the majority of cases (almost 80% in this specific example), thus providing a remarkably well-calibrated phylogeny of a character-rich clade of ubiquitous marine benthic invertebrates.
Key Words: molecular clock • rates of evolution • echinoids
|
Introduction |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
Congruence between independent sources of data is one of the most persuasive arguments when assessing the reliability of a phylogenetic hypothesis (De Queiroz et al. 1995
; Cunnigham 1997). Although some striking differences initially existed between morphological and molecular-based phylogenies, many of these are gradually being resolved as characters are added or reassessed, denser sampling is carried out, new fossils are found, and better methods of analysis adopted. For example, morphological and molecular data initially pointed to different sister-group relationships for the whales (Cetacea) within mammals. This conflict has now been resolved by the discovery of new fossils (Gingerich 2005), which demonstrated that previous morphological hypotheses were in error. Conversely, Xia et al. (2003) have shown that a bird–mammal linkage in 18S ribosomal RNA (rRNA) (which conflicted with morphological evidence) could be explained by misalignment of sequences, inappropriate generalization of base-frequency parameters over the whole tree, and poor sequence quality. The comparison of morphological and molecular data sets provides an important cross-check on the reliability of results.
In recent years, attention has shifted to another aspect where concordance between morphological and molecular estimates has been hard to achieve, namely, the dating of divergence times. Here workers seem to be more polarized in their outlook, either critical of the reliability and accuracy of molecular clock methods to date divergence times (e.g., Benton 1999; Rodriguez-Trelles et al. 2002; Benton and Ayala 2003) or dismissive of the quality of the fossil record (e.g., Easteal 1999). Until now, debate has logically focused on areas where the two approaches give the most divergent results, specifically the radiation of the metazoan phyla and on the Mesozoic origins of birds and of mammals (see Smith and Peterson 2002; Aris-Brosou and Yang 2003; Donoghue and Smith 2003; Peterson et al. 2004; Pisani et al. 2004; Blair and Hedges 2005; Ho et al. 2005; Peterson and Butterfield 2005; Welch et al. 2005). In all 3 cases, however, the fossil record is far from satisfactory. The fossil record of birds, for example, is miserable by comparison to almost all other fossil groups, with approximately two-thirds of fossil species from the Mesozoic still represented only by a single specimen (Fountaine et al. 2005), implying that the fossil bird record is still very poorly sampled. Similarly, the pre-Cambrian fossil record of bilaterian metazoans (Eumetazoans) is scant or nonexistent. In order to advance this debate, molecular clock methods need to be more widely tested in groups where the fossil record is considerably better and the phylogeny unequivocally established. Apart from the pioneering work of Pérez-Losada et al. (2004) on barnacles (Thoracica), no studies have looked at the reliability and consistency of molecular clock methods using marine invertebrate groups with a good fossil record.
Sea urchins (Echinoidea) are a diverse group of marine invertebrate deuterostomes (Schultz 2005; Smith 2006; Smith et al. 2004). Their multielement skeleton is complex, providing a large number of phylogenetically informative characters. This skeleton also preserves well, creating a rich fossil record that has been the focus for much paleontological research. Importantly, echinoid taxonomy, which was largely established in the great monographic works of Mortensen (1928, 1935, 1940, 1943a, 1943b, 1948a, 1948b, 1950, 1951), is based almost exclusively on skeletal characters and so is equally applicable to living and fossil taxa.
Today there are some 900 extant species distributed in about 50 families and 14 orders. Previous phylogenetic studies (e.g., Smith 1988; Smith et al. 1992; Littlewood and Smith 1995; Lessios et al. 1999, 2001; Jeffery et al. 2003; Lee 2003; Stockley et al. 2005) have generally found good levels of congruence between morphological and molecular estimates of relationship. With their good fossil record and well-established phylogenetic relationships, echinoids should provide a model system against which to examine the performance of molecular methods of dating.
Here we compare molecular and paleontological estimates of divergence times for Echinoidea as an empirical cross-check on the reliability of the methods and assumptions. We do this by 1) constructing phylogenies from both molecular data and traditional morphological data to arrive at the best-supported tree, 2) estimating the quality of the echinoid fossil record by calibrating this tree against the observed record of first occurrences, 3) estimating divergence times using a molecular clock method and a variety of relaxed molecular clock models applied to the molecular data, and 4) quantifying the match between the observed times of appearance of clades in the fossil record and the results derived from molecular data.
|
Materials and Methods |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
Taxa Included
To construct our phylogeny, we compiled morphological information and gene sequence data for 46 genera (listed in table 1) with representatives from 28 families and 13 of the 14 orders of living echinoid. Where sequence data for multiple species of the same genus existed, a strict consensus sequence was constructed. Because we restricted our analyses to only those regions of the genes that could be unambiguously aligned across all genera, sequences of congeneric species were effectively identical. However, the 2 species of Araeosoma were retained separately as these showed modest amounts of sequence divergence.
|
To compare divergence times based on the fossil record and molecular data, we selected 1 taxon from each family (asterisked in table 1). These were chosen, after examining our initial molecular phylogenetic analysis (fig. 2), to avoid taxa showing anomalously long or short branches, thereby creating a "partially linearized tree" (sensu Takezaki et al. 1995
).
|
Morphological data were rooted by outgroup comparison using the fossil Archaeocidaris, a late stem-group echinoid. Data on this taxon were taken from a remarkably well-preserved Carboniferous species (Lewis and Ensom 1982
). Molecular data were rooted using a combination of representatives from each of the other 4 classes of echinoderm as outgroup: the starfish Asterias, the ophiuroid Ophiocanops, the crinoid Antedon, and 2 holothurians, Cucumaria and Psychopetes. In regions where the ingroup and outgroup sequences were too divergent for meaningful alignment, outgroup sequences were scored as unknown to avoid spurious rooting.
Morphological Characters
Morphological characters were compiled from published analyses, notably from Smith (1988), Littlewood and Smith (1995), Jeffery et al. (2003), and Stockley et al. (2005). In total, 119 characters were scored, 29 of which are multistate (see Supplementary Material online). Two multistate characters (characters 6 and 8 in supplementary table 1, Supplementary Material online) reflect a clear and unambiguous ontogenetic sequence of character states and were, therefore, treated as ordered; the remainder was left unordered. The great majority of characters relate to skeletal features of the adult, which forms the basis for the classification of the group. Only 3 characters relate to the larval skeleton as larval characteristics have been shown to be more homoplastic than characters based on adult skeletal morphology (Smith et al. 1996; Smith and Littlewood 1997). Characters were obtained from direct observation or from modern descriptions in the literature and scored on the basis of the states shown in the species for which sequence data were available. The complete character listing and data matrix are provided as supplementary data (Supplementary Material online).
Molecular Characters
Three rRNA genes, 2 nuclear (18S small subunit [18S] and 28S large subunit [28S]) and 1 mitochondrial (16S large subunit [16S]), were sequenced. These were selected to encompass a range of different evolutionary rates and have been used successfully in a wide range of metazoan phylogeny studies aimed at resolving divergences over the past 250 Myr, including echinoderms (e.g., Littlewood and Smith 1995; Littlewood et al. 1997; Jeffery et al. 2003; Winchell et al. 2004; Stockley et al. 2005). Approximately 630 bp from the 3' end of the 16S gene, 1,250 bp from the 5' end of the 28S gene, and the entire 18S gene (ca. 1,800 bp) were sequenced.
For details of tissue selection and DNA extraction and amplification methods see Littlewood and Smith (1995) and Stockley et al. (2005). Sequence data were obtained for both forward and reverse reads. An initial multiple sequence alignment was made using MacClade (Maddison DR and Maddison WP 2001) followed by alignment of sequences by eye. Areas of high variability, for which no reliable alignment across the different orders could be made, were excluded from further data analysis. In total, out of 3,990 bp, 3,226 were alignable and 449 were phylogenetically informative. All sequences are lodged with GenBank/European Bioinformatics Institute under accession numbers listed in table 1. The aligned data matrix is available as supplementary data (Supplementary Material online).
Analytical Methods
Phylogenetic analysis of the morphological data was carried out with maximum parsimony as the optimality criterion, using the Macintosh version of PAUP* (4.0b10 [Altivec]) (Swofford 2002). Because of the large number of taxa included, we used a heuristic search method, with 1,000 random additional replicates and tree bisection reconnection branch swapping. Node support was tested by bootstrapping with 250 replicates and by clade decay analysis (Bremer 1994).
All 3 genes were combined for phylogenetic and molecular clock analyses. The program Modeltest (version 3.06) (Posada and Crandall 1998) was used to analyze each data set and produce an appropriate nucleotide substitution model. We used the general time reversible (GTR) + 
+ I model (rates set to gamma, with 6 rate categories). Bayesian inference analyses were conducted using a separate GTR + + I model for each data partition independently and also for the combined 3-gene analysis, thus allowing separate estimates for each model parameter per data set. Bayesian analyses were performed using MrBayes (Huelsenbeck and Ronquist 2003). The number of generations permitted was 5,000,000 with 4 chains, and the 50% majority rule consensus tree was constructed from the non–burn-in trees. Maximum likelihood (ML) analyses were implemented using subtree pruning regrafting (Hordijk and Gascuel 2005) under the best-fitting model (GTR + + I). Support for the nodes in the ML tree was estimated using the bootstrap (100 replicates).
Combined morphological and molecular data were analyzed using both parsimony and Bayesian methods. The likelihood model developed for morphological data by Lewis (2001) was used for the morphological partition, and separate, unlinked GTR + + I models were used for the molecular partitions.
Data Congruence Tests
The appropriateness of combining the morphological and molecular data sets was tested using the partition homogeneity test (Farris et al. 1994), as implemented in PAUP*. A Templeton (1983) test of data heterogeneity was performed on the trees that were obtained from morphological and molecular data, to determine whether they could be considered suboptimal estimates of the same underlying topology. We also used the recommended approach of Wiens (1998) of comparing support levels for nodes in the morphological and molecular analyses. We identified problematic areas as nodes where morphological and molecular data pointed to incongruent groupings with strong support, as indicated by bootstrap proportions of >70% or posterior probabilities of >95%.
Stratigraphic Congruence Tests
Benton's (1995, 2001) Relative Completeness Index (RCI) was used to measure the fit of stratigraphic data to cladogram topology. This measures the amount of missing range that must be added to make stratigraphic record fit the phylogeny and has the advantage over Huelsenbeck's (1994) Stratigraphic Consistency Index or the Spearman rank correlation proposed by Norrell and Novacek (1992) in that it takes account of the relative size of the mismatch over the entire tree (Hitchin and Benton 1997). We calculated RCI for the trees derived from the combined morphological and molecular data (fig. 3), using both parsimony and Bayesian analysis. We divided the fossil record into 5-Myr intervals and, for each sister-group pairing, used the earliest occurrence of either to establish the minimum time of origin of both sister groups. The geological ages of the earliest fossil representative of each clade at the family level or above included in this analysis are listed in table 2. Supporting evidence for the paleontological dating of nodes is provided in Appendix.
|
|
|
Where a family is the only representative of a larger clade, the oldest member of the more inclusive clade is given. For example, the Cidaridae is the only included family-level representative of the subclass Cidaroidea, sister group to all other living echinoids (fig. 1). Whereas the fossil record of Cidaridae is not much older than 150 Myr, stem-group members of the Cidaroidea extend back to 255 Myr. In such cases, we always use the earliest representative of a member of the most inclusive taxonomic group to date the divergence.
|
Molecular Estimates of Divergence Times
All molecular estimates of divergence times were calculated for the 28 family representatives selected (asterisked in table 1) plus outgroups, as explained in Materials and Methods. The likelihood ratio test (LRT) was used to test the null hypothesis that the data evolved under a molecular clock. This was done using Modeltest in "LRT calculator mode," after having estimated the likelihood for the molecular tree (fig. 2) under the GTR +
+ I model, both imposing a molecular clock and after having removed this assumption. As the molecular clock assumption could be rejected (P < 0.00001 with 29 degrees of freedom), divergence time estimates were calculated from the combined 16S, 18S, and 28S rRNA sequences using a variety of relaxed molecular clock approaches. We used the fully parametric Bayesian method of Thorne et al. (1998), the nonparametric rate smoothing (NPRS), and the semiparametric penalized likelihood (PL) methods of Sanderson (2002). However, we also used a dating method imposing a strict molecular clock, the Langley–Fitch method (LF; Langley and Fitch 1974), to evaluate its performance for comparison. The Bayesian method of Thorne et al. (1998) was implemented using the software Multidivtime (Thorne and Kishino 2002), whereas NPRS, PL, and LF analyses were performed using r8s version 1.70 (Sanderson 2004). The r8s version 1.70 allows divergence times under PL to be estimated using either a logarithmic- or an additive-penalty function, and divergence times for the internal nodes were obtained using both types of penalty function (see Sanderson 2004). Optimal smoothing parameters for the PL and log-PL analyses were obtained by cross-validation (Sanderson 2002, 2004).
Both r8s and Multidivtime need, as input information, a tree topology with associated branch lengths. For the r8s analyses, branch lengths for the ML molecular tree (fig. 2) were calculated via likelihood using PAUP 4b10, under the best-fitting substitution model (GTR + + I), which had been reestimated using Modeltest after removing all the taxa that were not selected for the molecular clock analyses. For comparison, in the Bayesian analyses branch lengths were estimated with the software Estbranches (which is part of the Multidivtime package) under the F84 + substitution model (as suggested in the Estbranches manual). Software limitations did not allow the implementation of more complex substitution model (i.e., GTR + + I), and the specific parameters for the F84 + model were estimated using phylogenetic analysis by maximum likelihood (PAML) (Yang 1997). Each Bayesian analysis was run 4 times, comparing the divergence times obtained in order to estimate whether convergence was reached. For each analysis, the Markov Chain Monte Carlo chain was run for 1,000,000 cycles, sampling every 100 cycles. The first 200,000 cycles were considered burn-in.
For all Bayesian divergence time estimates, 95% credibility intervals were calculated. For PL, NPRS, and the LF analyses, confidence intervals around the estimated divergence times were obtained by bootstrapping (Sanderson 2004). Thousand bootstrapped data sets were generated using Seqboot (Felsenstein 2004), and for each bootstrapped data set, the branch lengths of the input tree were reestimated using likelihood (see also above). The results of these analyses were sets of 1,000 bootstrap trees with associated branch lengths, with the trees in each set having the same topology and bootstrapped branch lengths. Divergence times were then obtained, and, taking each of these 1,000 bootstrap trees and each node in the input tree, the standard deviation (SD) of its estimated age was calculated. The 95% confidence intervals around each clade's estimated age were approximated as (X ± 2SD), where X is the estimated age of a given node. All the computational steps in the estimation of the bootstrap confidence intervals for the penalized likelihood analyses were automated using several PERL scripts written by D.P.
In the PL, NPRS, and LF analyses the ingroup node was fixed at 265 Myr, which is 10 Myr prior to the first appearance of the oldest recognizable member of the crown group in the fossil record. For the Bayesian analyses of divergence times, we set a date for this node of 265 Myr with a SD of 1, whereas the prior date of the root node was set to 480 Myr, which represents the earliest occurrence of members of 2 of our outgroups, the Asteroidea and Ophiuroidea (Dean 2005). Two sets of analyses were then run, the first without any internal constraints (with the exclusion of the ingroup node) and a second in which 4 local calibration points were enforced. These 4 points were set as minimal divergence estimates, and the basal dichotomy (ingroup node) was taken as a fixed point. The 4 internal calibrations were selected from across the tree topology to provide constraints on local rate variation. The following local calibration points were employed: the first occurrence of Pedinoida (fig. 5, node 4) at 210 Myr, the first occurrence of Scutellina (fig. 5, node 11) at 55 Myr, the first occurrence of paleopneustid spatangoids (fig. 5, node 15) at 95 Myr, and the first occurrence of Temnopleuridae (fig. 5, node 25) at 45 Myr (all dates from Smith 2006).
|
Measuring Congruence of Divergence Estimates
For comparing the different methods of estimating divergence times, we use the following measure. Each node where we have an independent molecular and paleontological estimate of divergence time was given a score between 0 and 2: 2 if the fossil date lies within 1 SD of the mean of the molecular estimate, 1 if it lies between 1 and 2 SD of the molecular estimate, and 0 when the paleontological estimate falls outside 2 SDs of the molecular estimate. Our overall measure of congruence is then simply the score summed over all nodes divided by the total number of nodes x 2; the higher the score (maximum 1 and minimum 0), the better the agreement achieved between molecular and paleontological estimates of divergence on our tree.
|
Results |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
Phylogenetic Relationships
Parsimony analysis of morphological data found 4,605 equally parsimonious solutions of tree length (TL), 249 steps with a consistency index (CI) of 0.60 and retention index 0.88. A semistrict consensus of these trees (fig. 1) shows good resolution in all parts of the tree except among members of the Camarodonta. Some nodes are well supported, especially within the irregular echinoids, but there are a number of weakly supported areas.
The ML analyses of molecular data identified a topology very similar, but not identical, to that supported by morphological data (fig. 2). Of the few differences, only 1 taxon, Fellaster, is placed with strong support at different positions in the 2 rival trees. Other conflicting placements are weakly supported in either one or both of the trees. The morphological and molecular trees passed Templeton's test. Although the trees based on molecular data could not have been produced by the morphological data, the consensus tree of equally parsimonious trees from the morphological analysis was not significantly different in terms of TL from the molecular ML tree at P > 0.05 (number of differences = 31, rank sums = 331.0, P = 0.07 for Templeton's test, and P = 0.15 for Winning sites test). Therefore, the 2 sources of data, morphological and molecular, can be considered to be suboptimal estimates of the same underlying topology.
Because different nodes are strongly supported in the 2 trees, a combined data set arguably provides the most appropriate way to combine the strengths of the 2 data sets. Parsimony analysis of the combined data found 6 trees whose strict consensus is shown (fig. 3, left-hand side). Bootstrap and Bremer support values were moderate to high for most branches. The Bayesian analysis of the combined data produces a tree that is very similar to the parsimony tree, except for its placement of the echinothurioids (Araeosoma) and the pedinoid Caenopedina (fig. 3, right-hand side). The topology in this part of the tree is closer to that derived from analysis of molecular data alone.
Stratigraphic Completeness Estimate
Both Bayesian and parsimony trees derived from the combined morphological and molecular data (fig. 3) were calibrated against the known fossil record to calculate the minimum amount of fossil record that must be missing. For the parsimony tree (table 2, fig. 4), the duration of all branch lengths (observed and inferred) implied by the combined morphological and molecular tree is 3,530 Myr, of which missing (ghost) lineages implied at family level constitute approximately 10% of the total duration (360 Myr inferred minimal time missing). Furthermore, almost half the perceived mismatch arises from relationships within the Camarodonta, where relationships are least well resolved. Almost identical results were obtained using the Bayesian combined-data tree (not shown).
Molecular Divergence Estimates
An ML analysis of the molecular data matrix for the 28 family representatives (plus outgroups) generated the tree shown in figure 5, and this was used to estimate molecular divergence times. When only the basal node of the ingroup was fixed, estimated divergence times varied considerably according to the specific method applied (table 3, fig. 6). The LF method failed to estimate the great majority (>80%) of paleontological nodes correctly, for the most part greatly underestimating divergence dates. The PL approach with additive-penalty function also performed poorly. However, the other 3 methods (NPRS, Bayesian, and log-PL) performed more or less equally well, having 17–18 of the 26 internal nodes congruent with paleontological estimates (i.e., the paleontological estimate lies within 2 SDs of the molecular estimate). Overall, the PL approach using logarithmic-penalty function gave the closest match to paleontological estimates.
|
|
|
|
When the additional 4 local calibration points were enforced, all molecular estimation methods performed reasonably well (table 4), having 17–18 of the 26 internal nodes congruent with paleontological estimates (figs. 5 and 7). However, in these conditions PL with additive-penalty function and LF greatly improved their rates of success, whereas the performance of NPRS was distinctly poorer compared with when only the ingroup node was constrained. Focusing on the number of nodes where the molecular estimate encompasses the paleontological estimate within 1 SD of the mean, the log-PL method outperforms other methods with NPRS performing the poorest.
Surprisingly, under these conditions LF performed better than all the relaxed methods we applied except for log-PL. This implies that relaxed methods should not be considered a priori better than methods for dating divergence times imposing a molecular clock. Methods imposing a single, global, rate of evolution could return biased results, but the same is true for relaxed models as correctly noted, for example, by Welch et al. (2005) and Ho et al. (2005).
It is also interesting to note that log-PL and the Bayesian method of Thorne et al. (1998) were the most consistent, correctly estimating a high proportion of divergence times, no matter whether internal constraints different from the ingroup node were enforced. However, log-PL always performed better than the method of Thorne et al. (1998). It is, however, difficult to evaluate what causes this difference. The method of Thorne et al. (1998), as it is currently implemented, can only use simple models (e.g., F84 + ), but it build these models into the likelihood calculations, whereas Sanderson's (2004) software (r8s) relies on trees with branch lengths that must be previously estimated using other software (e.g., PAUP) and does not build the models into the likelihood calculations. However, in this way, it allows using more complex models and hence more accurate branch length estimations. If the use on simple and perhaps suboptimal model will be discovered to be the cause of the difference in performance between the 2 methods, then we should expect that the Bayesian method of Thorne et al. (1998) would improve as more models are integrated into it.
In any case, when the molecular estimates and paleontological estimates for each node are regressed, there is a strong and highly significant correlation with r2 = 0.91 (for log-PL estimates). For the 8 cases where there is significant mismatch, paleontological data mostly underestimate the date of divergence.
Irrespective of which specific method was applied, mismatch between molecular and morphological estimates of divergence time is confined to the same 3 areas of the cladogram, 2 nodes near the base of the tree (the divergence of Euechinoidea; fig. 5, nodes 1 and 2), 4 nodes within the Clypeasteroida (fig. 5, nodes 8 and 11–13), and 2 within the Camarodonta (fig. 5, nodes 23 and 24).
When the tree derived from the combined-data analysis (fig. 3) was used as the base for calculating divergences instead of the molecular tree, the overall congruence between molecular and paleontological estimates decreases (table 5), and the camarodont and clypeasteroid nodes that were incongruent in the previous analyses are also incongruent in this analysis. However, the basal nodes now show a better fit, with molecular and paleontological estimates for nodes 1 and 2 now congruent and in close agreement.
|
|
Discussion |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
There has been much recent discussion about the validity of the different methods used to estimate divergence times from molecular data but few empirical tests in groups with a good fossil record. The study of Pérez-Losada et al. (2004)
tested a variety of different methods on the phylogeny of barnacles and their relatives but examined only a small number of nodes. In this study, we have adopted a similar approach but applied across a much larger number of paleontologically well-dated nodes. The fossil record and molecular-based divergence dates provide effectively independent estimates for the 26 internal nodes on our family-level tree. Furthermore, we specifically used the topology derived from only molecular data (rather than the combined morphological and molecular data tree) in order to set a more stringent test of the accuracy of the molecular approach.
It has been suggested that approaches that can accommodate rate variation in their inference procedures generate significantly more realistic results (e.g., Kishino et al. 2001; Sanderson 2002; Thorne and Kishino 2002; Pérez-Losada et al. 2004). However, whether these methods truly return better results seem strongly dependent on the model used to estimate rate variation, on the calibration points used (Welch et al. 2005), and, in the case of the Bayesian methods, on the prior on-divergence times chosen. For example, Aris-Brosou and Yang (2003) obtained high concordance with the fossil record from their molecular clock analysis; however, the major determinant of their divergence times was not their data, it was their prior (Ho et al. 2005; Welch et al. 2005). This also seems to be the case in our study, which found that although the LF method performed poorly when only the root is fixed, the inclusion of internal calibration points surprisingly increased its performance, whereas this same operation negatively affected NPRS. Furthermore, it is interesting to note that with our data, when LF behaved poorly, it underestimated most of the divergence time, contradicting the claim of Benton and Ayala (2003) that molecular clock methods by default overestimate divergence times. The correspondence we obtain between paleontological and molecular estimates of divergence using both the fully parametric Bayesian method of Thorne et al. (1998) and the semiparametric PL methods of Sanderson (2002) is reassuringly strong. In approximately 70% of cases, morphological and molecular estimates are congruent, that is, the paleontological estimate falls within the 95% confidence error bars of the molecular estimate. The use of internal fossil calibration points made little difference for log-PL and the Bayesian method of Thorne et al. (1998), and this consistency of performance is encouraging, suggesting that these methods should probably be favored over others. Overall, for our data, log-PL was best at recovering known divergence times, although employing better substitution models to estimate branch lengths in the Bayesian approach of Thorne et al. (1998) might be expected to improve the success rate of this method.
The few nodes that consistently show a mismatch between paleontological and molecular divergence estimates are restricted to 3 parts of the tree (fig. 5). These affect basal nodes 1 and 2 and nodes within the Echinoida clade (nodes 23 and 24) and clypeasteroid clade (nodes 8 and 11–13).
There are 3 reasons why disagreement between paleontological and molecular estimates of divergence might arise: inaccuracy of the phylogenetic reconstruction being used, incompleteness of the fossil record, or methodological problems in the way molecular estimates are derived. In order to test the first possibility, we calculated molecular divergence dates using the tree constructed from the combined morphological and molecular data (fig. 3), which we take as our best-supported estimate of phylogenetic relationships. This tree differs most obviously from the molecular tree (fig. 2) in the relationships of the basal nodes and the monophyly of clypeasteroids. Using the combined-data tree reduces the correspondence between molecular and paleontological estimates of divergence in some parts of the tree but improves the fit of basal nodes 1 and 2 considerably. Consequently, the poor match between paleontological and molecular estimates of divergence times seen in the 2 basal nodes in the molecular analyses could be the result of there being a suboptimal arrangement of basal taxa in the molecular tree. Significantly, however, mismatch between paleontological and molecular divergence estimates in the 2 other regions of the tree remains when using the best overall supported tree. In these cases, therefore, it is unlikely that the mismatch is simply a problem of inaccurate phylogenetic reconstruction.
To check whether inadequacies in the fossil record could explain the mismatch, we calibrated the combined morphology and gene tree against the observed fossil record (fig. 4). Although this implies that the fossil record of echinoids is relatively well sampled overall, almost half the missing record identified is concentrated in the Camarodonta, suggesting that this group is seriously underrepresented in the fossil record. The reason for this is obvious. It is notoriously difficult to assign fossil echinoids to specific families within the Camarodonta because Mortensen (1943a, 1943b) established his taxonomy on the basis of structures (details of pedicellarial anatomy) that are rarely preserved in fossils. Consequently, it is likely that the fossil record in this part of the tree is better than that suggested by our data and that it is simply a lack of diagnostic morphological characters that creates the apparent mismatch. The fossil record may well be there, it is just that we cannot currently place many fossil taxa with any degree of confidence. On this part of the tree, molecular estimates of divergence are always deeper than paleontological estimates (figs. 5–7
). Consequently, the observed mismatch of molecular and paleontological estimates of divergence at nodes 23 and 24 is probably due to inadequacies of the fossil record, rather than due to errors in the molecular-based estimate.
For the remaining 4 nodes (nodes 8 and 11–13), the mismatch between paleontological and molecular estimates of divergence times is less easily explained. All these nodes refer to divergences within clypeasteroids, a morphologically complex and intensively studied clade (Durham 1955; Kier 1982; Mooi 1987, 1990). Furthermore, they have a highly distinctive synapomorphy (large numbers of pores perforating all ambulacral plates) that makes them immediately recognizable in the fossil record, even from small fragments. The group primarily inhabits shallow water, living in environments that are well represented in the fossil record. Consequently, the chance that crown-group clypeasteroids have a deep, hidden fossil record stretching back into the Early Cretaceous seems improbable from what we know of the echinoid fossil record.
Clypeasteroid divergence times estimated from molecular data may be afflicted by two potentially significant problems. First, the topology of this part of the molecular tree differs from that supported by the combined analysis of morphological and molecular data that we take as our best estimate of the phylogeny (cf. figs. 3 and 5). Specifically, the 2 cassiduloid families group within the clypeasteroids in the molecular tree (fig. 5) rather than as their sister group (fig. 3). However, even when the combined-data tree is used as the model to estimate divergence times (log-PL analysis, table 5), the paleontological and molecular-based estimates of divergence time remain strongly incongruent. Errors in the topology of the tree cannot, therefore, explain why molecular and paleontological estimates of divergence times are so incongruent in this part of the tree.
A second potential problem arises from the very uneven rate of molecular evolution shown by the cassiduloids and clypeasteroids on the combined-data tree. Whereas cassiduloids show 0.0051–0.0057 substitutions per site, all 5 clypeasteroids have rates about 4 times faster (0.018–0.026 substitutions per site). This may be a large amount of rate heterogeneity for molecular methods to accommodate. We predict that denser taxon sampling together with two or more local tie points (minima) might help generate a closer match between paleontology and molecular data in this part of the tree.
In conclusion, because morphological and molecular data both point to closely similar phylogenetic relationships among echinoid clades (with the exception of the clypeasteroid Fellaster) and we estimate that their fossil record is relatively complete, we have the first real opportunity to compare the accuracy of molecular and paleontological methods of estimating divergence times empirically. In approximately 70% of nodes tested, paleontological and molecular methods give congruent estimates of divergence dates using methods that allow for rate variation over the tree. Although this is pleasantly reassuring, it is not as good as we were hoping for. For those nodes where paleontological and molecular estimates of divergence time are incongruent, paleontological data mostly underestimate divergence times. In some cases, the problem clearly lies with the paleontological data because of the poor preservation potential of key diagnostic characteristics by which paleontologists recognize those clades. A second problem may be that the molecular tree being used to establish divergence times is suboptimal in its arrangement of certain branches. When these problems are discounted, only a small proportion of nodes (ca. 15%) show a significant mismatch between molecular and paleontological estimates where the error may be the fault of our molecular techniques. Although care is still needed in selecting calibration points when using molecular data to estimate divergence times, our study demonstrates that, so long as a realistic model of rate variation is applied (see Welch et al. 2005), modern parametric and semiparametric approaches that assume rate heterogeneity can and do generate realistic divergence time estimates in the great majority of cases.
|
Supplementary Material |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
Figures 6 and 7, in color, and the full description of morphological characters and character states together with the data matrix and the aligned sequences used in our analyses given as supplementary tables 1, 2, and 3 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).
|
Appendix |
|---|
|
TOPAbstractIntroductionMaterials and MethodsResultsDiscussionSupplementary MaterialAppendixAcknowledgementsReferences |
|---|
Supporting evidence for the paleontological dating of nodes in this study is given in this section.
Paleontological dates for nodes 1–26 in figures 5–7 are fixed as follows.
- Root: Cidaroidea–Euechinoidea Divergence. These 2 clades have fundamentally different lantern and perignathic girdle structures, both derived compared with the arrangement seen in Paleozoic late stem-group archaeocidarids. Kier (1984a) recognized that this divergence probably occurred in the Permian, and Smith and Hollingworth (1990) showed that the well-dated Late Permian (Kazanian) "Miocidaris" keyserlingi possessed cidaroid synapomorphies in both the structure of its lantern and perignathic girdle. As the oldest demonstrable member of the Cidaroidea, M. keyserlingi sets the minimum time of crown-group divergence at 255 Myr.
- 1: Echinothurioida–Acroechinoidea Divergence. Echinothurioida are the only extant clade of echinoids whose tests have remained imbricate, and, consequently, they have the poorest fossil record of any of the groups considered here (see Smith and Wright 1990
). This divergence is, therefore, dated by the first occurrence of a member of the Acroechinoidea, which has, as its primary synapomorphies, a fully tessellate test, the lack of multiple peristomial ambulacral plates, and the presence of true ambulacral plate compounding. "Hemipedina" hudsoni (Kier 1977), though relatively poorly known, has a solid, nonimbricate test plating and compound ambulacral plating characteristic of an Acroechinoidea. It comes from the Norian of Oman and dates the origin of Acroechinoidea at around 220 Myr.
- 2: Aulodont–[Echinacea + Irregularia] Divergence. The major synapomorphy distinguishing the Echinacea and Irregularia from aulodonts, is their possession of keeled, as opposed to grooved, teeth. Isolated teeth are not uncommon from the Middle Triassic (Carnian) St Cassian beds, but all are grooved (Vadet 1999
). The lantern of the Late Triassic (Rhaetian) Diademopsis serialis is known in detail (Smith 1981) and is aulodont in structure. The earliest known keeled teeth come from the Lower Jurassic (Pliensbachian; Smith 1981; Markel 1978). However, Kier (1977, p. 33) argued that the Late Triassic (Rhaetian) Pseudodiadema silbinense (which is more correctly placed in the genus Stereopyga) represented the oldest echinacean on the basis of general test characteristics. These occurrences suggest that this divergence had definitely occurred by 195 Myr and probably by 210 Myr.
- 3: Diadematidae–[Aspidodiadematidae + Pedinidae] Divergence. The pairing of Aspidodiadematidae and Pedinidae as sister groups to the exclusion of Diadematidae in the molecular analysis is unexpected. Traditionally, the Diadematidae and Aspidodiadematidae have been grouped together in the order Diadematoida and distinguished from members of the order Pedinoida by their crenulated tubercles (e.g., Mortensen 1940
; Smith 1981). No morphological synapomorphy exists uniting the molecular clade Aspidodiadematidae + Pedinidae, and so the date of this divergence is set by the oldest member showing synapomorphies of any 1 of the 3 constituent families. The Upper Triassic (Rhaetian) Diademopsis serialis is generally placed in the Pedinidae (Smith 1981, 2006) as is the Lower Jurassic (Pliensbachian) Hemipedina, which displays the derived dicyclic pedinid apical disc. Synapomorphies of the Aspidodiadematidae include their distinctive apical disc plating, which is only very loosely connected to the corona, and their distinctive sphaeridial pit arrangement. Both occur in the Lower Jurassic (Pliensbachian) Gymnotiara varusense (Smith 2006). Divergence must have occurred before 195 Myr and probably by 210 Myr, the date used here.
- 4: Pedinoida–Aspidodiadematidae Divergence. See Discussion for node 4. The oldest member of the Aspidodiadematidae is G. varusense, from the Pliensbachian at 195 Myr (Smith 2006
). The oldest member of the Pedinoida is Diademopsis serialis from the Rhaetian, at 210 Myr.
- 5: Echinacea–Irregularia Divergence. Both these groups have keeled teeth but differ in their apical disc structure: Echinacea having a normal endocyclic apical disc with 5 gonopores in contrast to the more derived apical disc arrangement of primitive irregular echinoids in which the periproct is displaced to the posterior behind the posterior genital plates and genital plate 5 has lost its gonopore. The oldest putative echinacean is the Rhaetian Stereopyga silbinense (redescribed by Kier 1977
). However, there is nothing to definitely place this as a derived Echinacea rather than a stem-group Echinacea + Irregularia. Cladistic analysis consistently identifies the Sinemurian Jesionekechinus hawkinsi as the earliest member of the Irregularia (Smith and Anzalone 2000; Barras forthcoming), establishing divergence at 200 Myr.
- 6: Echinoneoida–Microstomata Divergence. Among irregular echinoids, Microstomata is more derived than its sister group Echinoneoida in having evolved specialized ambulacral zones in the form of aboral petals and adoral phyllodes (Barras forthcoming). The Lower Jurassic (Toarcian) Galeropygus sublaevis has distinct phyllodes and a subanal groove and establishes the date of this split as no younger than 185 Myr.
- 7: Atelostomata–Neognathostomata Divergence. Both atelostomates and neognathostomates have long, well-documented stem groups. The synapomorphies that distinguish basal atelostomates are mostly associated with apical disc arrangement, atelostomates having a "stretched" apical disc with ocular plates 2 and 4 intercalated between the anterior and posterior pair of genital plates, as opposed to the more compact ethmophract disc structure of early neognathostomates (Smith 1981
; Barras forthcoming). The earliest echinoid to show the derived atelostomate disc plating is Aulacopygus caudatus from the Late Bajocian (Smith 2006), setting this divergence as no later than 175 Myr.
- 8: Fellaster–[Cassiduloida + Scutellina] Divergence. As discussed in the text, the position taken by Fellaster in the molecular tree is in strong contradiction to that in traditional taxonomies based on morphology (Mortensen 1948
b; Smith 1981; Mooi 1987, 1990). Traditional taxonomies place Fellaster as a member of the Clypeasteroida along with its sister group the Scutellina, whereas the molecular topology implies that Clypeasteroida are diphyletic, with the 2 cassiduloids as closer to members of the Scutellina than Fellaster. A further complication is that the Cassiduloida, as currently recognized, is a paraphyletic group (Smith 2001). Cladistic analysis of morphological data suggests that the cassiduloid families Cassidulidae and Echinolampadidae form a clade, whereas a third cassiduloid family, Apatopygidae (not included in the molecular analysis), represents the closest living sister group to Fellaster and the Scutellina (Smith 2001). This latter grouping is based on a derived pattern of ambulacral plating seen in apatopygids and undoubted stem-group Clypeasteroida and implies that the Clypeasteroid–[Cassidulidae + Echinolampadidae] split can be traced back in the Late Cretaceous (Nucleopygus angustatus being the earliest recognizable representative of the Apatopygidae + Clypeasteroida total group). Given the molecular topology, divergence at this node is thus placed at 105 Myr.
- 9: Cassiduloida–Scutellina Divergence. As for node 8, above.
- 10: Cassidulidae–Echinolampadidae Divergence. The phylogenetic relationships of cassiduloids remain problematic because of the small number of convincing synapomorphies available from morphology (Smith 2001
). Convincing cassidulids (Rhyncholampas) and echinolampadids (Vologesia) are, however, present and clearly differentiated by the Maastrichtian (70 Myr; Smith and Jeffery 2000) with echinolampadids extending back to the Santonian (Hungaresia ovum; Smith 2006). This places their divergence at no younger than 85 Myr.
- 11: Scutelliformes–Laganiformes Divergence. The most paleontologically useful synapomorphies characterizing these 2 clades are 1) the single terminal interambulacral plate unique to laganines and 2) the complex arrangement of buttressing internal meshwork in the Scutellina (e.g., Mooi 1987
). Both Sismondia and Echinocyamus have interambulacra that terminate in a single plate, and both extend back to the Eocene. However, the very oldest members of each genus (the Lower Eocene Echinocyamus gurnahensis and Sismondia logotheti) are not sufficiently known to confirm that their interambulacra show the laganiform arrangement (e.g., Roman and Struogo 1994). The oldest scutelliform, Eoscutum doncieuxi, is also of Lower Eocene and has the apomorphic dense internal buttressing of a scutelliform (Roman 1990; Smith 2006). Divergence is thus taken as having occurred around 50–55 Myr.
- 12: Echinocyamidae–Laganidae Divergence. Echinocyamidae are small laganiforms with interambulacral zones ending adapically in a single plate and with internal buttressing composed of radial buttresses only. Unlike other laganiform families, there is no certainty that this is a monophyletic clade, and it could include pedomorphic forms of other laganiforms. Small laganids might be extremely difficult to distinguish from taxa currently placed in the Echinocyamidae. Forms attributed to Echinocyamus extend back into the Eocene (Kier 1968
), with E. gurnahensis being the oldest (Roman and Struogo 1994). Laganidae are larger forms with a well-developed internal skeleton of both concentric partitions and adradial bars. The highly stellate nature of their basicoronal plates, with their radially directed points, is also a synapomorphy. One stem-group branch of the Laganidae shows an additional synapomorphy of having pseudocompound plating in their petals and is found as far back as the Middle Eocene (Kier 1968, 1980). The 2 groups were thus clearly separated by the Middle Eocene with a divergence probably in the Early Eocene, at 50–55 Myr.
- 13: Mellitidae–Astroclypidae Divergence. Mellitidae and Astriclypeidae are sister taxa, united by having a well-developed microcanal system internally, spines and tubercles on the oral surface clearly differentiated into food-gathering and locomotory areas, and the periproct opening in the first postbasicoronal interambulacral plate. Mellitidae is the more derived taxa in possessing an anal lunule (Smith 2006
). Although both clades possess ambulacral notches or lunules, these are constructed - 1: Echinothurioida–Acroechinoidea Divergence. Echinothurioida are the only extant clade of echinoids whose tests have remained imbricate, and, consequently, they have the poorest fossil record of any of the groups considered here (see Smith and Wright 1990