Molecular Biology and Evolution 19:2005-2021 (2002)
© 2002 Society for Molecular Biology and Evolution
Fossil Calibration of Molecular Clocks and the Divergence Times of Geminate Species Pairs Separated by the Isthmus of Panama
Naos Marine Laboratory, Smithsonian Tropical Research Institute, Republic of Panama and Molecular Systematics Laboratory, Natural History Museum of Los Angeles County
| Abstract |
|---|
|
|
|---|
Calibration of nucleotide sequence divergence rates provides an important method by which to test many hypotheses of evolution. In the absence of an adequate fossil record, geological events, rather than the first appearances of sister taxa in the geological record, are often used to calibrate molecular clocks. The formation of the Isthmus of Panama, which isolated the tropical western Atlantic and eastern Pacific oceans, is one such event that is frequently used to infer rates of nucleotide sequence divergence. Isthmian calibrations assume that morphologically similar "geminate" species living now on either side of the isthmus were isolated geographically by the latest stages of seaway closure 3.13.5 MYA. Here, I have applied calibration dates from the fossil record to cytochrome c oxidase-1 (CO1) and nuclear histone-3 (H3) divergences among six pairs of geminates in the Arcidae to test this hypothesis. Analysis of CO1 first and third positions yield geminate divergences that predate final seaway closure, and on the basis of CO1 first positions, times for all six geminates are significantly greater than 3.5 Myr. H3 sequences produce much more recent geminate divergences, some that are younger than 3.1 Myr. But H3-derived estimates for all arcid geminates are not significantly different from both 0 and 15 Myr. According to CO1, one of the two most divergent pairs, Arca mutabilis and A. imbricata, split more than 30 MYA. This date is compatible with the fossil record, which indicates that these species were morphologically distinct at least 1621 MYA. Across all CO1 nucleotide sites, divergence rates for arcids are slower than the rates reported for other taxa on the basis of isthmian calibrations, with the exception of rates determined from the least divergent species pair in larger surveys of multiple transisthmian pairs. Rate differences between arcids and some taxa may be real, but these data suggest that divergence rates can be greatly overestimated when dates corresponding to final closure of the Central American Seaway are used to calibrate the molecular clocks of marine organisms.
| Introduction |
|---|
|
|
|---|
"Geminate species" (Jordan 1908
Among the handful of taxa for which transisthmian molecular sequences have been collected for multiple species pairs, simultaneous isolation of all geminates seems an unlikely hypothesis (Bermingham and Lessios 1993
; Knowlton et al. 1993
; Bermingham, McCafferty, and Martin 1997
; Knowlton and Weigt 1998
; Lessios 1998
; Lessios, Kessing, and Pearse 2001
; Marko and Moran 2002
). Although rate variation can explain some differences in divergence among geminate pairs (Lessios 1979
; Bermingham and Lessios 1993
; Bermingham, McCafferty, and Martin 1997
), most patterns of varying molecular sequence divergence across the isthmus are explained by the hypothesis of staggered evolutionary divergence of geminates over time (Knowlton and Weigt 1998
; Lessios, Kessing, and Pearse 2001
). The fossil record, which does not provide evidence of when gene flow stopped between EP and Caribbean populations, indicates that extinctions associated with the uplift of the isthmus also accrued in a gradual rather than an abrupt fashion, starting more than 10 MYA in the Middle Miocene (Jackson et al. 1993
; Coates and Obando 1996
; L. S. Collins 1996
; Jackson, Jung, and Fortunato 1996
; Beu 2001
; Roopnarine 2001
; Vermeij 2001
). Similarly, bursts of origination of taxa (e.g., L. S. Collins 1996
), as early as the Late Miocene (5.210.4 Myr), suggest that early stages of seaway constriction may have presented a barrier to some taxa. Nevertheless, the assumption that single geminate pairs were isolated at the time of complete cessation of water exchange between EP and WA is still used routinely for clock calibrations (Baldwin et al. 1998
; Chenoweth et al. 1998
; Metz, Gomez-Gutierrez, and Vacquier 1998
; Schubart, Diesel, and Hedges 1998
; Donaldson and Wilson 1999
; Hellberg and Vacquier 1999
; Lessios et al. 1999
; Tringali et al. 1999
; Van Oppen, Willis, and Miller 1999
; Hellberg, Moy, and Vacquier 2000
; Quinteiro, Vidal, and Rey-Mendez 2000
; Hickerson and Ross 2001
; Stillman and Reeb 2001
).
On the other hand, some molecular data can be used to address the validity of assuming that gene flow was maintained until the time of final seaway closure. Because of an uninformative fossil record for most of the marine taxa investigated to date, studies of molecular evolution involving multiple geminates must assume that the least divergent species pair was isolated at the time of final seaway closure, providing the best estimate for rates of molecular evolution (Bermingham, McCafferty, and Martin 1997
; Knowlton and Weigt 1998
; Lessios 1998
; Lessios, Kessing, and Pearse 2001
). Ideally, calibration of clocks with organismal divergences in the fossil record, that are completely independent of isthmus formation, is necessary to determine when geminate pairs were isolated. Only one study to date has used calibration points from the fossil record to estimate the divergence of sister species distributed on both sides of the isthmus. This lone nonisthmian calibration indicates preisthmian divergence times of 5.38.5 Myr for three tropical American gastropod geminate pairs (data in Collins 1989
; for analyses see Cunningham and Collins 1994
; Collins et al. 1996
; Knowlton and Weigt 1998
).
Determining whether nonisthmian calibrations generally result in slower inferred rates of nucleotide substitution and correspondingly older divergence times for transisthmian pairs is fundamentally important for understanding the rates of molecular evolution and the history of biotic separation between the tropical WA and the EP oceans. Here, I analyzed mitochondrial cytochrome c oxidase-1 (CO1) and nuclear histone-3 (H3) sequences from six pairs of transisthmian arcid species, as part of a molecular phylogenetic analysis of the bivalve mollusc family Arcidae (arcids or ark shells). Strong morphological similarities have led to wide acceptance that these six arcid species pairs merit geminate status (Olsson 1961
, 574 p.; Abbott 1974
, 662 p.; Keen 1974
; Vermeij 1978
, 416 p.). Recent morphological analyses indicate that multivariate statistical analyses are required to distinguish living transisthmian arcids (Marko and Jackson 2001). Discriminant analyses of specimens of Arca from the fossil record also demonstrate that the shell shapes characterizing the living transisthmian species pair, A. mutabilis and A. imbricata, are both present in 16- to 21-Myr-old fossil deposits from the Early Miocene (Marko and Jackson 2001), suggesting that some arcid geminates were morphologically distinct long before final seaway closure. Using phylogenetic trees based on analyses of partial CO1 and H3 sequences, I have employed three fossil-based divergence times (i.e., nonisthmian events) to calibrate the arcid molecular clock and infer the absolute divergence dates for the six geminate pairs, including A. mutabilis and A. imbricata. Because most studies involving geminates assume genetic isolation at final seaway closure, the present data set has implications for the assumptions made and conclusions drawn from molecular clock studies that use geological events rather than the first appearances of taxa in the fossil record as calibration points.
| Methods |
|---|
|
|
|---|
Taxon Sampling
Between September 1997 and June 1999, specimens were collected at various localities in tropical America and elsewhere (table 1 ). An IndoWest Pacific representative of Glycimerus was used as an out-group in all analyses. Other molecular analyses suggest a close relationship between the Glycimeridae and the Arcidae (Distel 2000
|
DNA Extraction, Amplification, and Sequencing
Genomic DNA was extracted from each specimen by proteinase K digestion of tissue in 2x CTAB (2% hexadeclytrimethylamonium-bromide, 1.4 M NaCl, 20 mM EDTA, 100 mM Tris, 0.2% ß-mercaptoethanol, pH 8.0) for 3 h to overnight followed by two chloroformisoamyl alcohol (24:1) extractions, isopropanol precipitation, and two washes with 70% ethanol.
To amplify and sequence a 603-bp portion of CO1 and a 324-bp portion of the H3 gene, I used the primers of Folmer et al. (1994)
and W. F. Ponder and D. R. Lindberg (unpublished data, but see GenBank accession number AF033701 for primer sequences), respectively. For each gene, 12 µl of each genomic DNA extraction was used as a template in a 25-µl polymerase chain reaction (PCR) that consisted of 10 mM Tris, 50 mM KCl, 20 mM MgSO4, 0.6 mM dNTPs, 0.2 µl Taq polymerase, and 0.5 mM of each primer. For CO1, each reaction was subjected to 40 cycles of 94°C for 30 s, 4045°C for 30 s, and 6572°C for 1.52 min. Most CO1 amplifications required thermal cycler ramp speeds of 1°C/s, both down to and up from the annealing step. H3 sequences were amplified with 40 cycles of 94°C for 30 s, 5055°C for 30 s, and 72°C for 30 s to 1 min.
PCR amplifications that produced a single band of the expected size when visualized on an agarose gel were excised and separated from the agarose by a spin column (Qiagen). Some CO1 products were faintly visible on agarose gels, and were excised and reamplified. Products were then sequenced directly in both directions with fluorescently labeled dye-terminators (Applied Biosystems, Inc.). In total, 49 CO1 sequences were collected from 27 nominal species. A total of 23 H3 sequences were obtained from 23 of the nominal taxa for which CO1 was sequenced.
Phylogenetic Analyses
Because of an absence of insertions and deletions, both CO1 and H3 sequences were easily aligned. Phylogenetic trees were constructed with PAUP* 4.0 (Swofford 2001
) using maximum parsimony (MP) and maximum likelihood (ML). For MP searches, all nucleotides were weighted equally, searches were heuristic, and included 20 random additions of taxa, with tree-bisectionreconnection branch swapping. The reliability of trees obtained with MP was assessed with 500 nonparametric bootstrap replicates (Felsenstein 1985
).
Because ML branch length estimates and hypothesis tests of the molecular clock can vary as a result of substitution model assumptions (Yang 1999
; Buckley, Simon, and Chambers 2001
), before any ML searches, Modeltest (ver. 3.06) (Posada and Crandall 1998
, 2001
) was used to conduct hierarchical likelihood ratio tests (HLRTs) (Huelsenbeck and Crandall 1997
) as a method to choose the most appropriate substitution model for each of the complete CO1 sequences, complete H3 sequences, and combined data sets. Modeltest uses a neighbor-joining tree topology derived from Jukes-Cantor (Jukes and Cantor 1969
) genetic distances to infer all substitution model parameters before calculation of log-likelihood ratio tests. HLRTs resulted in the selection of the General Time Reversible model (Rodriguez et al. 1990
) with invariant sites and gamma-distributed rate heterogeneity (GTR + I + G) for each of the CO1 and combined data sets. A model in which transversions and transitions occur at different rates (but all transversions occur at an equal rate) with rate heterogeneity among sites (TrN + G) was selected for the H3 data set.
Statistical Tests of Molecular Clocks
I used likelihood ratio tests (LRTs) (Felsenstein 1981
) to test the null hypothesis of a Poisson-distributed clock. Because a clock was easily rejected for complete CO1 sequences (see Results), I subdivided the CO1 data set according to codon position and used HLRTs to choose substitution models for each subset of the data. The following models were chosen: (1) GTR with equal nucleotide frequencies and rate heterogeneity among sites (SYM + G) for first positions, (2) Felsenstein-81 model (Felsenstein 1981
) with rate heterogeneity among sites (F81 + G) for second positions, (3) a transition model (Posada and Crandall 1998
) with rate heterogeneity among sites (TIM + G) for third positions, and (4) Kimura's 3-parameter model (Kimura 1981
) with unequal nucleotide frequencies and rate heterogeneity among sites (K81uf + G) for first and second positions combined.
Calibration Points
Given the importance of fossil calibration points in analyses involving molecular clocks (Lee 1999
; Yoder and Yang 2000
), I chose three nodes that received relatively strong bootstrap support in analyses of the combined CO1-H3 data set. I also considered only calibration points in which the identities of fossils had been critically evaluated by several other authors.
The first calibration point (C1) was based on the split between the subfamilies Noetiinae and Striarcinae. The oldest representatives of the Noetiinae in the fossil record are from the Early Cretaceous (Upper Aptian) of Lebanon, aged 83118 Myr (Whitfield 1891
; Vokes 1946
; Cox et al. 1969
). The oldest species assigned to Striarcinae are known from the Late Cretaceous (Santonian), dated at 8386 Myr (Richards 1958
). Therefore, 83118 Myr was used as the range of dates for the C1 calibration node.
The second calibration point (C2) marks the divergence between the Anadarinae and the clade comprising the closely related barbatid subgenera Fugleria and Cucullaearca. The earliest member of the Anadarinae is from the Late Cretaceous (Early Campanian) in strata aged 7983 Myr (Conrad 1869
, 1870
; Sheldon 1917
; Reinhart 1935
; Richards 1958
). The subgenus Cucullaearca has a much older geologic history than does Fugleria (Reinhart 1937
), first known 8389 MYA in the Late Cretaceous (Early Senonian) of North Carolina (Stephenson 1923
; Cox et al. 1969
). Therefore, a range of 7989 Myr dates was considered for the C2 calibration point.
The third calibration point (C3) is the most recent and corresponds to the split between Anadara sensu stricto (a subgeneric group represented in the molecular phylogeny by A. tuberculosa and A. similis) and the subgenus Grandiarca (Anadara grandis). Both Anadara s.s. and Grandiarca first appeared and diversified in the Early Miocene 1623 MYA (Nelson 1870; Spieker 1922
; Reinhart 1943
; Noda 1966
, 1986
, 1991
; Eames 1967
; Cox 1969
). The anadarid subgenus Grandiarca is combined frequently with Larkinia (not represented in the phylogenies) by some authors (Cox et al. 1969
; but see Abbott 1974
; Keen 1974
), but because Larkinia also first appears in the Early Miocene, this taxonomic uncertainty does not affect the C3 calibration date.
Calculation of Geminate Divergence Times
For each CO1 and H3 data set, I calculated the divergence times of geminate pairs in two ways. First, divergence rates for each calibration point were estimated separately using the average sequence divergence between taxa in calibration clades. For each calibration point, I calculated sequence divergence rates using each of the upper and lower limits from the fossil record. Second, I used a simple weighted linear regression of time on sequence divergence, with the intercept forced through the origin to determine the relationship between time since divergence and sequence divergence (Hillis, Mable, and Moritz 1996
). Regressions models were then used to calculate the predicted divergence times of geminates. To determine the statistical significance of inferred divergences, I calculated the 95% confidence interval for each divergence time (i.e., prediction intervals, see Hillis, Mable, and Moritz 1996
). In the regression approach, I used the midpoint between the upper and lower bounds from the fossil record as the values of time for each calibration point. Sequence divergences were calculated patristically, and geminate divergences were corrected for intraspecific variation by calculating the net nucleotide divergence (Nei and Li 1979; Edwards and Beerli 2000
).
| Results |
|---|
|
|
|---|
Phylogenetic Analyses: CO1
A total of 603 nucleotide positions of CO1 were aligned for the 49 CO1 sequences. Across all sites, 422 positions were variable and 355 were parsimony informative. For complete CO1 sequences, both ML and MP recover nearly identical tree topologies (figs. 1 and 2 , respectively) with some exceptions. First, Acar is monophyletic under ML (fig. 1 ), whereas in the MP reconstruction this taxon does not form a clade because of the placement of Barbatia divaricata (fig. 2 ). Second, the position of B. cancellaria and the relationship between Acar and Arca also differ between the MP and ML trees. ML and MP analyses of the first and second CO1 positions also produced trees (not shown) highly similar to those in figures 1 and 2 .
|
|
Monophyly of CO1 sequences from each geminate species and each geminate species pair is strongly supported by MP bootstrap percentages (fig. 2 ). Although the genus Barbatia (comprising the subgenera Acar, Fugleria, Cucullaearca, and Barbatia s.s.) is polyphyletic, most other generic and subgeneric taxa erected on morphological grounds are monophyletic with respect to CO1 (i.e., Anadara, Arca, Arcopsis, Fugleria, and Cucullaearca, figs. 1 and 2 ). Nodes representing all three calibration points also were present in trees from ML and MP analyses of CO1 sequences (figs. 1 and 2 ). Although the C3 node received strong support (100% of MP bootstrap replicates, fig. 2 ) and the sister-clades descended from all three calibration points received good support (71%100%, fig. 2 ), the C1 and C2 nodes received only poor to moderate MP bootstrap support: no greater than 60% of bootstrap replicates included these two nodes (fig. 2 ). If the CO1 third positions are excluded, bootstrap percentages for the C1 and C2 nodes increase marginally to 64% and 62%, respectively (fig. 2 ).
For CO1 third positions, MP produces a tree with substantial phylogenetic structure (fig. 3A ), much of which is shared with the tree produced in analyses of complete sequences (59% of nodes of the tree in fig. 1 ). Use of ML in the third position, however, gains five additional nodes shared by the trees in figures 1 and 2 , including all three calibration nodes (fig. 3B ).
|
Phylogenetic Analyses: H3
For the nuclear H3 gene, 23 sequences of 361 nucleotides were aligned. A total of 100 sites were variable, 84 of which were parsimony informative. An MP tree based on H3 sequences (fig. 4A ) was similar to the CO1 trees. All included geminates formed sister taxa and were each well supported by bootstrap percentages (greater than 84%), except for the most divergent species pair, B. gradata and B. domingensis, which received poor support (less than 50%). All three calibration points also were present in the MP tree based on H3 sequences. The H3 tree differed from the CO1 topology in that Arca, Noetia, Arcopsis, and all but one member of the subgenus Acar did not form a monophyletic clade. Arca also was not monophyletic in the MP tree based on H3 sequences (fig. 4A ).
|
An ML tree based on H3 sequences produced a tree similar to the topology of the MP tree, except that the subgenus Cucullaearca was placed within Anadara (fig. 4B ); thus, the C2 calibration point was not present in this tree. In the ML tree, the B. gradata and B. domingensis geminate pair also did not form a sister pair (fig. 4B ).
Phylogenetic Analyses: Total Evidence
Phylogenetic analysis of the combined H3 and CO1 data sets produces a topology (fig. 5
) very similar to the CO1 and H3 trees built with MP. Both ML and MP searches of the combined data set yield identical tree topologies. All included geminates formed sister taxa and were well supported by MP bootstrap percentages, with the exception of the A. nux and A. chemnitzii pair. Total evidence MP analyses also yield more robust bootstrap support for the calibration nodes than do either the CO1 or the H3 MP trees. With all the nucleotides weighted equally, bootstrap percentages are 73, 72, and 100 for the C1, C2, and C3 nodes, respectively (fig. 5
). If CO1 third positions are excluded, bootstrap percentages increase to 76 and 88 for the C1 and C2 nodes, respectively (fig. 5
).
|
Tests of the Molecular Clock
Likelihood ratio tests comparing trees in which the clock was both relaxed and enforced are presented in table 2 . For each of the complete CO1 sequences, second positions, and first and second positions combined, a clock is easily rejected. But for each of the first and third positions, a clock could not be rejected. A clock also could not be rejected for complete H3 sequences.
|
Saturation
To investigate CO1 sequences for evidence of substitution saturation, I first plotted the number of nucleotide differences (uncorrected) for each codon position versus the corrected sequence divergences across all sites for ingroup taxa (fig. 6 ). Both first and second positions start to show the effects of multiple substitutions at sequence divergences >0.75, very close to the divergences observed at the C1 and C2 calibration nodes (fig. 6A and B ). Third positions accumulate differences much more rapidly, suggesting that saturation starts occur at divergences <0.25 across all sites, less than the average number of differences between sequences at the C3 node; at the C1 and C2 nodes, the third positions appear nearly completely saturated. Effects of saturation are less pronounced in the nuclear H3 sequences, even at third positions (fig. 7B ). Most H3 differences are at third positions, with only a fraction of differences at first positions (fig. 7 ). No H3 sequences exhibited differences at second positions.
|
|
I also compared uncorrected pairwise H3 differences with each of the CO1 differences at first positions and the corrected CO1 third position divergences. CO1 first position differences show little evidence of saturation until H3 differences approach
25 (fig. 8A
), a value slightly larger than the number of H3 differences at the C2 calibration node. Corrected CO1 third position divergences begin to plateau at the same point with respect to H3 differences (fig. 8B
). Therefore, the C1 calibration appears substantially biased because of the effects of multiple substitutions, with respect to both uncorrected first position differences and corrected third position divergences. Plots in which H3 divergences and CO1 first position differences are corrected for multiple substitutions are nearly identical to those in figure 8
.
|
Divergence Rates and Divergence Times of Geminates
Divergence rates are presented in table 3 . Across all CO1 sites the C2 calibration yielded the largest divergence rates, with complete sequences diverging at 1.21%/Myr. The C1 calibration yielded the smallest rates for all CO1 data sets, as was expected because of evidence of saturation. Excluding the C1 calibration, third positions diverge at a rate of 4.03%6.84%, whereas first positions evolve much more slowly, at rates of 0.15%0.46%/Myr. Not surprisingly, the H3 sequences evolve much more slowly than CO1.
|
Because CO1 sequence divergences at the C1 node are clearly affected by saturation at both first and third positions, CO1-derived divergence times based on only the C2 and C3 calibrations are listed in table 4 . Across table 4 the youngest inferred divergence time for any geminate pair is 9.94 Myr for Arcopsis solida and A. adamsi. Although the confidence interval is large, the divergence time for the Arcopsis pair is significantly greater than 3.5 Myr, the oldest date from the geologic record associated with final seaway closure. Divergence times based on CO1 first positions also indicate that all the other five geminate pairs split significantly earlier than 3.5 Myr. The divergence date for B. reeveana and B. candida is significantly greater than 10 Myr, and divergence times for the three other geminates are significantly greater than 20 Myr.
|
CO1 third positions yield similar results with respect to the point estimates for geminate divergence times, but the 95% confidence intervals indicate that third positions estimates are less precise (table 4 ). The reduced statistical power of the CO1 third position clock is best exemplified by the confidence limits associated with the split of the least divergent pair (B. reeveana and B. candida), bound by -1.1 and 24.1 Myr (table 4 ). Thus, using CO1 third positions, the hypothesis that geminates diverged at the time of final seaway closure is not rejected for the three relatively less divergent geminates (Arcopsis solida and A. adamsi, B. reeveana and B. candida, and B. illota and B. tenera). Isolation at the time of final seaway closure is, however, rejected easily for the three most divergent geminate pairs (Anadara nux and A. chemnitzii, Arca mutabilis and A. imbricata, and B. gradata and B. domingensis), all of which have inferred divergence dates significantly greater than 10 Myr.
Analysis of H3 sequences produces much younger divergence times for arcid geminate pairs (table 5 ). According to the H3 data set, three of the four geminates diverged less than 3 MYA. But, like CO1 third positions, predicted values for divergence times based on H3 sequences have very broad confidence limits: the 95% confidence intervals for the predicted divergence times for all geminates include both 0 and 15 Myr. Even the most divergent species pair (B. gradata and B. domingensis) has confidence limits for its divergence that range from -12.9 to 47.7 Myr for H3-derived estimates.
|
| Discussion |
|---|
|
|
|---|
Although the molecular clock potentially provides a temporal perspective on patterns of evolutionary change, calibration is often problematic for taxa that lack a fossil record. When appropriate fossil divergences are not available, biogeographic or geologic events are often used instead. Final closure of the Central American Seaway is the most widely used calibration point for marine taxa, in part because of its relative recency but also because it is one of the most intensively studied geologic events of the Cenozoic. Rather than assuming that final seaway closure 3.13.5 MYA was the cause of disruption of gene flow between geminate species, I have used calibration dates independent of isthmus formation to infer the divergence times for six geminate species pairs in the bivalve family Arcidae. Application of this approach to CO1 first-codon positions results in geminate divergence times that significantly exceed 3.5 Myr, allowing rejection of the hypothesis that all six arcid geminates diverged at the time of final seaway closure.
Calculation of the 95% confidence limits for the predicted divergence times, however, reveals that the error associated with these estimates is formidable. Even for CO1 first positions, which exhibit relatively greater precision than do either CO1 third positions or H3 sequences, the confidence intervals for most geminate divergences span a period similar to the entire time during which the isthmus was formed. The greater precision of CO1 first positions compared with either CO1 third positions or complete H3 sequences (and the much younger point estimates of divergence based on H3 sequences) may reflect the variation in rates of substitution among these different types of sites. CO1 first positions diverge more slowly than CO1 third positions but not nearly as slowly as nuclear H3 sites (table 3
). Over the temporal scale of interest, CO1 first positions may accumulate enough substitutions so that the number of observed differences is not only eclipsed by stochastic variation but also accrues quite a few changes so that saturation does not introduce significant error in the estimation of sequence divergences (fig. 6A
). Given that CO1 third positions diverge at rates greater than 4%/Myr (table 3
), estimation of sequence divergences among calibration taxa likely entails much greater error because of the effects of multiple substitutions at single-nucleotide sites (fig. 6C
) and, thus, much greater proportional deviations about the regression line used to estimate the divergence times of geminates. Greater variances associated with the relatively small numbers of substitutions observed at loci like H3 may explain why the confidence intervals for H3 estimates are so large (Hillis, Mable, and Moritz 1996
). Substitutions at the nuclear H3 gene are rare and may occur too infrequently to resolve divergence times less than 30 Myr. For example, the taxa descended from the C3 calibration (which split approximately 20 MYA) have accumulated only two differences, both of which occurred in the Anadara grandis lineage; no changes have occurred in the A. tuberculosa and A. similis clade over approximately 20 Myr. Therefore, on the basis of the C3 calibration node, approximately one nucleotide difference is expected every 10 Myr, too few to expect temporal resolution over the last 20 Myr.
Regardless of the causes of the variance in estimates among the data sets considered here, the relatively older divergence dates derived from CO1 sequences are most compatible with the direct evidence available from the fossil record. Discriminant analyses of recent and fossil specimens of Arca mutabilis and A. imbricata indicate that these taxa were morphologically distinct during the early Miocene, 1621 MYA (Marko and Jackson 2001). The general agreement between the fossil record and the CO1 molecular clock estimates for the divergence of the Arca geminate pair indicates that this split likely represents a speciation event that was not a consequence of isthmus formation. CO1 divergences also indicate that other geminate divergence may be equally ancient, likely unrelated to formation of the isthmus. For example, according to CO1 first positions, both Barbatia gradata and B. domingensis, and Anadara nux and A. chemnitzii likely split before the end of the Oligocene (23 Myr), more than 10 Myr before the earliest hypothesized disruptions of deep water exchange between the EP and the Caribbean.
But in the two least divergent geminate pairs (A. solida and A. adamsi, and B. illota and B. tenera), which exhibit divergence times derived from CO1 first positions as recent as 9.94 Myr, disruption of gene flow between the Caribbean and the EP could have been caused by a combination of climate change and oceanographic effects related to the rising isthmus. During the latest Middle Miocene (11.812.9 MYA), ocean depths over the isthmus were reduced significantly by uplift of the Panama sill (Duque-Caro 1990b
). This period is immediately followed by drops in sea level, the formation of distinctive benthic EP invertebrate assemblages not present in the WA, and the first dispersal of terrestrial mammals (8.09.3 Myr) between North and South America (Marshall et al. 1979
; Marshall 1985
; Webb 1985
; Duque-Caro 1990b
). Equatorial shifts in the latitudinal distributions of benthic foraminiferans in the northern EP at this same time suggest intensification of the southward flow of the California Current down the coast of Central America, which, combined with partial emergence of the isthmus, may have produced a transient circulation barrier between the EP and the Caribbean from 11.8 to 7.0 MYA (Duque-Caro 1990b
). Nearshore oceanographic conditions over the rising isthmus also may have limited the dispersal of organisms with planktonic larvae. Among snapping shrimp in the genus Alpheus, geminate pairs that exhibit the smallest transisthmian sequence divergences are those that live in mangroves, nearshore environments that were likely the very last habitats separated by the isthmus (Knowlton and Weight 1998
). Because neither Arcopsis nor Barbatia (Cucullaearca) geminates are restricted to these environments, these taxa could represent species pairs isolated by the relatively early oceanographic changes described above rather than by the final termination of shallow water flow 3.1 MYA.
Several potential confounding factors might also affect interpretations of the evolutionary significance of sequence divergences between transisthmian taxa. First, inadvertent amplification and sequencing of pseudogenes (mitochondrial genes transferred to the nuclear genome), which appear to be common in Alpheus (Williams and Knowlton 2001
), may confound interpretations of varying divergence times among arcid geminate species pairs. No stop codons or frameshift mutations were found in any sequences in this study, and no sequences appeared to be mixtures of more than one PCR product. Nevertheless, very recently transferred pseudogenes may exhibit none of these characteristics. Although it is impossible to rule out this hypothesis, the inclusion of pseudogene sequences among the six geminate pairs likely would not fundamentally affect the primary result of this study because pseudogenes are expected to evolve more slowly than the functional mtDNA copies (Arctander 1995
; Sunnucks and Hales 1996
; Lopez et al. 1997
; Williams and Knowlton 2001
) and would therefore result in apparently younger geminate divergences.
Second, incomplete taxon sampling may explain why some transisthmian taxa exhibit such larger divergences than do others. For example, the subgenus Acar, which includes one of the three most divergent geminate pairs (B. gradata and B. domingensis), is also represented in tropical America by two relatively rare (although morphologically distinct, Bartsch 1931
) EP species, Barbatia bailyi and B. rostae, not presently included in the phylogeny. Similarly, the anadarid subgenus Cunearca, which contains the A. nux and A. chemnitzii pair, also is represented in tropical America by several species not yet sampled (Olsson 1961
; Keen 1974
). For other subgeneric arcid groups, however, incomplete taxon sampling likely is not a potential explanation for deep divergences of hypothesized geminate pairs: with the exception of one species (Arca zebra), no other extant species of Fugleria, Cucullaearca, Arca, and Arcopsis is known from either the EP or the WA. Sampling of species from other geographic regions, such as the Central Pacific and IndoWest Pacific (e.g., Arca ventricosa, B. divaricata) and the tropical eastern Atlantic (B. plicata), is yet to reveal that species from other regions might be related more closely to one member of a hypothesized geminate pair. But because cryptic species are nearly ubiquitous among marine taxa (Knowlton 1993
), additional geographic sampling within nominal transisthmian taxa could reveal additional species. Morphometric analyses of shell shape among living populations (particularly for the subgenus Acar, see Marko and Jackson 2001) indicate that interpopulation differences within some arcid species rival the morphometric differences observed across the isthmus, further suggesting the potential of cryptic species within either the EP or the Caribbean.
Lastly, extinction is an important hypothesis that has yet to be investigated in detail with respect to the evolutionary significance of geminate species, and it can also explain why some hypothesized geminates have relatively large sequence divergences. For example, if one member of a geminate pair isolated at the time of final seaway closure was lost to extinction after formation of the isthmus, the next most closely related species found on the same side of the isthmus would appear to be the corresponding partner in a geminate pair (Schneider 1995
; Lessios 1998
). Given that at least 70% of some molluscan subgeneric groups were lost to extinction in the Caribbean after seaway closure (Stanley 1986
; Jackson et al. 1993
), transisthmian taxa that survived in both the EP and the Caribbean may be rare.
Among the arks, most of the hypothesized geminates are thin-shelled taxa found nestling in gaps and crevices on reefs and rocky shores (Arca, Barbatia, and Arcopsis), which suggests that taxa that produce thick shells and burrow in soft bottoms (Anadara and Noetia) are less likely to have survived on both sides of the isthmus after seaway closure. But given the relatively ancient divergence times suggested by the CO1 data sets, extinction of reef and rocky shore species may have been as great as for burrowing taxa. Some evidence from the fossil record suggests that differential extinction of geminates was important after isthmus formation. For the lone geminate pair for which fossil data have been analyzed, a small number of specimens corresponding to the shell shapes of both A. mutabilis and A. imbricata are known from some Pleistocene deposits in the Caribbean, indicating that both species may have survived for some time in the WA after final seaway closure (Marko and Jackson 2001). Because oceanographic conditions in the Caribbean were more similar to those in the EP before seaway closure (Coates and Obando 1996
; L. S. Collins 1996
), it is also possible that the surviving Caribbean members of geminate pairs could be restricted geographically to small areas of the WA that today are more "Pacific-like" with respect to temperature, seasonality, and biological productivity (Vermeij 1978
; Petuch 1981
; Vermeij and Petuch 1986).
Ancient divergence times for arcid bivalve species pairs currently separated by the isthmus suggest that some clocks calibrated with final seaway closure may overestimate the rates of molecular evolution and may therefore underestimate the divergences of other taxa (for discussions, Gorman, Kim, and Rubinoff 1976
; Lessios 1979
, 1998
; Vawter, Rosenblatt, and Gorman 1980
; Hart, Byrne, and Smith 1997
; Knowlton and Weigt 1998
). For example, across all sites the arks exhibit CO1 divergence rates of 0.7%1.2%/Myr (table 3
). In contrast, sea star geminates in the genus Oreaster, which differ by 17.2%, yield an isthmian rate of 5.0% (Hart, Byrne, and Smith 1997
). Relatively large transisthmian divergences are also reported for Penaeus (13%) and Eucidaris (8.5%11.1%), yielding rates of 4.3%/Myr and 2.8%3.7%/Myr, respectively (Baldwin et al. 1998
; Lessios et al. 1999
). Among taxa with a fossil record, an isthmian calibration for the gastropod genus Tegula yields a rate of 2.4% (Hellberg and Vacquier 1999
), which is twice as fast as the maximum arcid rate (1.2%). Differences in rates between arks and other taxa may be real, but consideration of the available fossil and phylogenetic information suggest that CO1 substitution rates probably are overestimated at least for Tegula. Although Norissia, the closest out-group for Tegula (Hickman and McLean 1990
), is only sporadically known from the Pleistocene (Hickman and McLean 1990
) and Pliocene (unpublished data), its plesiomorphic morphological features indicate a basal position in the Tegulinae (Hickman and McLean 1990
; Hellberg 1998
) and, therefore, an origin equal to or earlier than that of Tegula, whose unique synapomorphies first appear 15 MYA in the fossil record (Hickman and McLean 1990
; Hellberg 1998
). But the range of CO1 sequence divergences between Tegula and Norissia spans 11.0%21.9% (Hellberg 1998
), yielding divergence times of 4.69.1 MYA under the 2.4% isthmian rate. Therefore, either Tegula is not monophyletic with respect to Norissia or the isthmian calibration overestimates CO1 divergence rates by at least 40%. Using the assumption that Tegula and Norissia diverged 15 MYA, the divergence rates (0.7%1.5%) are more similar to those for arcids, yielding divergence times of 5.110.1 Myr for the Tegula geminate pair. This divergence is qualitatively similar to those obtained for other rocky-shore gastropods geminates (5.38.5 Myr), in which calibrations were from the fossil record rather than from the isthmus (T. M. Collins 1989
, 1996
; Cunningham and Collins 1994
).
The divergence rates inferred for arcid bivalve and Tegula sequences calibrated with the fossil record are most similar to the CO1 rates based on isthmian calibrations in which surveys of several transisthmian pairs were first completed and the least divergent transisthmian pair was used to infer the rate (assuming that these pairs split
3 MYA). For example, on the basis of the smallest sequence divergences observed among 15 pairs of snapping shrimp in the genus Alpheus, Knowlton and Weigt (1998)
infer a rate of 1.4%/ Myr for CO1. Similarly, the least divergent echinoids among eight transisthmian pairs (Meoma, 4.55%, Zigler et al., unpublished data, cited in Lessios, Kessing, and Pearse 2001
; Diadema, 4.57%, Lessios, Kessing, and Pearse 2001
) yield a rate of 1.5%/Myr. A rate of 1.2% across all sites also was inferred after comparison of 19 transisthmian fish lineages (Bermingham, McCafferty, and Martin 1997
). Directly testing the assumption that the least divergent geminate pairs provide the best estimates for these rates is difficult because these other taxa lack the perspective of the fossil record. For arks, application of this assumption to CO1 first positions results in divergence times for the calibration nodes that cannot clearly be accommodated by the fossil record. For example, if the least divergent species pair, A. solida and A. adamsi, is assumed to have a divergence time of 3.5 Myr, the age of the C2 node must be less than 28 Myr (on the basis of CO1 first positions), which is more than 50 Myr before the first appearance of the C2 calibration taxa in the fossil record.
In summary, molecular sequence data from arcid bivalves indicate that despite strong morphological similarities, isolation of geminate species did not necessarily stem from the latest stages of closure of the Central American Seaway. As previous studies have shown, deep differentiation and speciation of marine taxa often does not require impenetrable physical barriers to dispersal, and many marine species are of relatively recent origin (Burton and Lee 1994
; Palumbi 1994
; Duffy 1996
; Miya and Nishida 1997
; Palumbi et al. 1997
; Burton 1998
; Marko 1998
; Barber et al. 2000
). Therefore, the use of calibration dates associated with the final geological steps associated with vicariance events may be uncertain for some marine organisms, leading to inflated rates of molecular evolution and other organismal traits with which molecular sequence divergences are frequently compared. Although earlier divergences of geminates are known as potential sources of error in the calibration of geminate molecular clocks (Knowlton and Weigt 1998
), molecular and paleontological evidence indicates that for at least one pair (Arca mutabilis and A. imbricata), speciation may have predated isthmus formation by more than 10 Myr. Calculation of confidence intervals indicates that the expected variance in accumulation of substitutions also presents a considerable source of error in these estimates, likely requiring additional estimates of divergence from independent genes evolving at appropriate rates (Hillis, Mable, and Moritz 1996
) as well as comparison with the fossil record.
| Supplementary Material |
|---|
|
|
|---|
The CO1 sequences reported in this paper are deposited in GenBank under accession numbers AF25347694, AF34564147, and AF41682040. Histone-3 sequences are under accession numbers AF41684164.
| Acknowledgements |
|---|
|
|
|---|
This work was supported by postdoctoral fellowships from the Smithsonian Institute, the Friday Harbor Laboratories (University of Washington), and the W. M. Keck Foundation (grants awarded to the Molecular Systematics Laboratory, Natural History Museum of Los Angeles County). N. Knowlton and S. Williams collected specimens from Brazil and the Cape Verde Islands, and Tahitian samples were collected by C. Thacker and D. Geiger. I thank N. Knowlton, J. Jackson, S. Williams, W. Toller, H. Lessios, C. Cunningham, C. Thacker, D. Geiger, D. Kizirian, A. Moran, two anonymous reviewers, and the associate editor for discussion and comments on the manuscript.
| Footnotes |
|---|
Stephen Palumbi, Reviewing Editor
Abbreviations: CO1, cytochrome c oxidase 1; EP, eastern Pacific; HLRT, hierarchical likelihood ratio test; H3, histone-3; LRT, likelihood ratio test; ML, maximum likelihood; MP, maximum parsimony; NJ, Neighbor-joining; WA, western Atlantic. ![]()
Keywords: bivalves
calibration
cytochrome oxidase-1
histone-3
fossil record
molecular clock
molluscs
speciation dates
substitution models ![]()
Address for correspondence and reprints: Peter B. Marko, Department of Marine Sciences, 12-7 Venable Hall, CB# 3300, The University of North Carolina at Chapel Hill, Chapel Hill, North Carolina 27599-3300. pmarko{at}unc.edu ![]()
| References |
|---|
|
|
|---|
Abbott R. T., 1974 American seashells. 2nd edition Van Nostrand Reinhold Co., New York
Arctander P., 1995 Comparison of a mitochondrial gene and a corresponding nuclear pseudogene Proc. R. Soc. Lond. Ser. B Biol. Sci 262:13-19[Medline]
Baldwin J. D., A. L. Bass, B. W. Bowen, W. H. Clark, 1998 Molecular phylogeny and biogeography of the marine shrimp Penaeus. Mol. Phylogenet. Evol 10:399-407[ISI][Medline]
Barber P. H., S. R. Palumbi, M. V. Erdmann, M. K. Moosa, 2000 Biogeography: a marine Wallace's line? Nature 406:692-693[Medline]
Bartsch P., 1931 The West American mollusks of the genus Acar U.S. Nat. Mus. Proc 80:1-4.
Bermingham E., H. A. Lessios, 1993 Rate variation of proteins and mitochondrial DNA evolution as revealed by sea urchins separated by the Isthmus of Panama Proc. Natl. Acad. Sci. USA 90:2734-2738
Bermingham E., S. S. McCafferty, A. P. Martin, 1997 Fish biogeography and molecular clocks: perspectives from the Panamanian Isthmus Pp. 113128 in T. D. Kocher and C. A. Stepien, eds. Molecular systematics of fishes. Academic Press, New York
Beu A. G., 2001 Gradual Miocene to Pleistocene uplift of the Central American Isthmus: evidence from tropical American tonnoidean gastropods J. Paleontol 75:706-720
Buckley T. R., C. Simon, G. K. Chambers, 2001 Exploring among-site variation models in a maximum likelihood framework using empirical data: effects of model assumptions on estimates of topology, branch lengths, and bootstrap support Syst. Biol 50:67-86[ISI][Medline]
Burton R. S., 1998 Intraspecific phylogeography across the Point Conception biogeographic boundary Evolution 52:734-745[ISI]
Burton R. S., B.-N. Lee, 1994 Nuclear and mitochondrial gene geneaologies amd allozyme polymorphism across a major phylogeographic break in the copepod Tigriopus californicus Proc. Natl. Acad. Sci. USA 91:5197-5201
Chenoweth S. F., J. M. Hughes, C. P. Keenan, S. Lavery, 1998 When oceans meet: a teleost shows secondary intergradation at an Indian-Pacific interface Proc. R. Soc. Lond. Ser. B Biol. Sci 265:415-420
Coates A. G., J. A. Obando, 1996 The geologic evolution of the Central American isthmus Pp. 2156 in J. B. C. Jackson, A. F. Budd, and A. G. Coates, eds. Evolution and environment in tropical America. University of Chicago Press, Chicago
Collins L. S., 1996 Environmental changes in Caribbean shallow waters relative to the closing tropical American seaway Pp. 130167 in J. B. C. Jackson, A. F. Budd, and A. G. Coates, eds. Evolution and environment in tropical America. University of Chicago Press, Chicago
Collins T. M., 1989 Rates of mitochondrial DNA evolution in transisthmian geminate species Doctoral dissertation, Yale University, New Haven, Conn
Collins T. M., 1996 Molecular comparisons of transisthmian species pairs: rates and patterns of evolution Pp. 303334 in J. B. C. Jackson, A. F. Budd, and A. G. Coates, eds. Evolution and environment in tropical America. University of Chicago Press, Chicago







