Molecular Claims of Gondwanan Age for Australian Agamid Lizards are Untenable
http://www.100md.com
分子生物学进展 2004年第11期
* Centre for Evolutionary Biology and Biodiversity, School of Earth and Environmental Sciences, University of Adelaide, and Natural Sciences Division, The South Australian Museum, Adelaide, Australia
E-mail: Andrew.Hugall@adelaide.edu.au.
Abstract
A recent mtDNA study proposes a surprisingly deep (150 MYA) divergence between SE Asian and Australasian agamid lizards, consistent with ancient Gondwanan vicariance rather than dispersal across the Indonesian Archipelago. However, the analysis contains a fundamental error: use of rates of molecular evolution inferred from uncorrected sequence divergence to put a time frame on a tree with branch lengths greatly elongated by complex likelihood and rate-smoothing models. Furthermore, this date implies that basal splits within agamids occurred implausibly early, at least 300 MYA (100 Myr before the first fossil lizards and coincident with the earliest fossil reptiles). Analyses of the mtDNA data using more appropriate methods and new information from nuclear (c-mos) sequences suggest a much more recent divergence between SE Asian and Australian agamids (around 30 MYA). Using two fossil boundary dates, bootstrapping the c-mos data gives a 95% confidence interval for this divergence time that is sufficiently recent (14–41 MYA) to exclude an ancient Gondwanan vicariance and is more consistent with Miocene over-water dispersal. As with the mtDNA, the c-mos data implies implausibly old basal divergences among agamids if a Gondwanan age is assumed for the Australasian clade. The analyses also highlight how methods for creating ultrametric trees (especially nonparametric rate smoothing) can greatly modify branch lengths and, thus, always require internal calibrations. The errors associated with inferred dates in the previous study (inferred through parametric bootstrapping) were also unjustifiably low, as this method only considers stochasticity in the substitution model and ignores much larger sources of uncertainty, such as variation in character sampling, tree topology, and calibration accuracy.
Key Words: molecular clock ? rate-smoothing ? ultrametric trees ? Agamidae ? Wallace line ? Gondwanan biogeography
Introduction
Most of the extant Australasian biota is widely accepted to have originated from two ancestral sources. Australasia had land connections with other Gondwanan continents (Antarctica and South America) until the early Paleocene, 55 MYA at the latest (Audley-Charles 1987; Woodburne and Case 1996; Metcalfe 2001). There followed a period of isolation as the Australasian plate detached from Antarctica and drifted northwards; the isolation ended perhaps as early as 30 MYA, as the tectonic collision with Southeast Asia formed a filter-bridge allowing colonization of terrestrial organisms from the north (Hall 1998; but see Discussion). The hiatus defines a dichotomy between two distinct components of the Australasian biota: (1) ancient lineages that were distributed across Gondwana and became isolated in Australasia and (2) recent lineages that represent invasions from the north, when Australia collided with the Asian plate. Clearly, for any endemic Australasian radiation, these two routes predict very different distributions regarding nearest extralimital relatives (South America versus Asia) and very different divergence times from those relatives (55 MYA versus 30 MYA).
The 64 or more species of agamid lizards ("dragons") constitute a sizeable component (13.5%) of the Australian lizard fauna and traditionally have been interpreted to have arrived by dispersal from the north (e.g., Cogger and Heatwole 1981; Schwaner et al. 1985; Greer 1989). This interpretation, however, has been challenged by a recent analysis of mtDNA data that inferred deep (150 MYA) divergence dates between Australasian agamids and their sister group (the SE Asian mainland Physignathus cocincinus) that could only be consistent with an ancient vicariance during the early breakup of Gondwana (Schulte, Melville, and Larson 2003; see also Moody [1980] and Macey et al. [2000]).
Schulte, Melville, and Larson (2003) used an approximately 1,600-bp segment of mtDNA to construct the phylogeny of Australasian and a sparser sampling of other agamids. Branch lengths were reconstructed using a parameter-rich maximum-likelihood (ML) model and then, because of significant rate variation, further modified using nonparametric rate smoothing (NPRS) (Sanderson 1997) to obtain an ultrametric tree. These reconstructed branch lengths were translated into absolute time using a calibration rate of 0.65% change per Myr along each lineage, calculated using uncorrected differences in Laudakia agamids (Macey et al. 1998), and also consistent with other ectothermic vertebrates (Weisrock et al. 2001). The divergence between the Australasian radiation and P. cocincinus was thus dated at around 150 MYA. The problem here is that a tree with highly "stretched" branches constructed under one model has been assigned absolute dates using a calibration rate inferred from (and applicable to) a tree with much shorter branches reconstructed under a very different model.
Here, we reanalyze the mtDNA data of Schulte, Melville, and Larson (2003) using more appropriate methods and taxon sampling. The results all suggest divergence dates around four to five times shallower (30 MYA) than proposed, with associated errors that are much greater than proposed. This more recent divergence is corroborated by analysis of new nuclear gene (c-mos) sequences.
Mitochondrial Data
We repeated the analysis of Schulte, Melville, and Larson (2003), firstly using their data matrix (TreeBase matrix N1326), tree, and model (see figure 1) in PAUP* (Swofford 2002). The GTR-gamma-inv (GTRgi) model adopted stretches the observed branch length between the Australasian radiation and P. cocincinus from 19% (uncorrected) to around 42%. NPRS as implemented in TreeEdit (Rambaut and Charleston 2001) further stretches the branch length to 88% divergence (type 1 method, "weight rate difference across root") or 95% divergence (type 3 method, "weight rate differences at all nodes with mean"). That is, the original uncorrected difference has been doubled by the GTRgi model and doubled again by NPRS. These greatly extended branches (inflated divergences) are then translated into time using an uncorrected divergence rate of 0.65% per Myr per lineage to give 150 Myr (see figure 1).
FIG. 1.— Agamid mtDNA data set N1326 and trees with inferred branch lengths, all drawn to the same scale, for relative dating of the Australasian-Leiolepis and Australasian-Physignathus cocincinus splits (arrowed). (A) Jukes-Cantor; (B) Jukes-Cantor made ultrametric with MLK constraint; (C) GTRgi model with MLK constraint; (D) GTRgi made ultrametric with NPRS type 3 (as implemented in TreeEdit). Most taxon names have been omitted to simplify the figure. Ratio of depths of Leiolepis split versus depth of P. cocincinus split: MLK = 2.16, NPRS = 1.52 (thus, NPRS increases the relative depth of the latter split by a factor of 1.42 compared with MLK). Schulte, Melville, and Larson (2003).
However, using the same calibration rate on the rest of the tree dates the dichotomy between the clade containing Leiolepis and the clade containing the Australasian taxa at around 250 MYA (227 MYA using type 1 NPRS, 254 MYA using type 3 NPRS, and 310 MYA using MLK constraint instead of NPRS ). This is more than 100 Myr before the oldest unequivocal fossil agamid, much older than even the 190-Myr putative acrodont described by Evans, Prasad, and Manhas (2002), and very soon after the first unequivocal reptiles (around 288 MYA [Lee 1999]). It is now a common procedure to test the fit of the MLK clock constraint using a likelihood ratio test (e.g., Felsenstein 1981) and, if significantly worse, to make the tree ultrametric using NPRS. In the face of substantial accurately measurable rate heterogeneity, NPRS may be preferred to MLK for the ultrametric transformation (Sanderson 1997, 2002); however, attention needs to be paid to the effect it has on the shape and size of the tree. In the above analysis, NPRS greatly increases the overall length of the tree and also alters the ratio of the depths of the nodes in question, compared with the effect of the MLK transformation. Here, NPRS lowers the ratio of deeper (Leiolepis) to shallower (P. cocincinus) node distances, "fattening" the shape of the tree (see figure legends).
The results of Schulte, Melville, and Larson (2003) are erroneous because an ML+NPRS tree with greatly lengthened, reconstructed branch lengths has been related to absolute time using a divergence rate inferred in other studies using uncorrected sequence divergence. The obvious way (Rambaut and Bromham 1998; Benton and Ayala 2003) to avoid this error is to employ internal calibration(s), so that branch lengths (and inferred rates) subtending the calibration node(s) are modified in exactly the same fashion as branch lengths (and inferred rates) throughout the rest of the tree. Macey et al. (2000) analyzed agamid phylogeny using the same region as Schulte, Melville, and Larson (2003). The data set of Macey et al. (2000) is more appropriate to testing the questions posed in Schulte, Melville, and Larson (2003) than the data set in the latter study for two reasons: (1) taxon sampling is much more balanced, with many more non-Australasian agamids included, and (2) there are multiple geologically dated internal calibration points within Laudakia that are all highly consistent with each other and broadly consistent with other studies (Macey et al. 1998; Schulte, Melville, and Larson 2003). The divergence used to calibrate the present tree is the isolation of the Iranian plateau Laudakia (L. microlepis, L. erythrogastra, and L. caucasia) from their nearest northern relatives (identified as the L. himalayana clade by Macey et al. [2000]), which Macey et al. (1998) suggested was associated with a tectonic event of 10 MYA (the Pamir-TienShan uplift). Similarly, use of all the possible Laudakia calibration points simultaneously (e.g., using penalized-likelihood rate smoothing [Sanderson 2002]) would also give virtually identical results, given they are located on the same region of the tree and highly consistent with each other.
Using Macey et al.'s (2000) alignment (http://Systbiol url, matrix macey-99-59) and tree (their figure 8), with the GTRgi model, we employ the same procedure—unconstrained GTRgi tree made ultrametric by (1) MLK and (2) NPRS (fig. 2). Analyses with MrBayes, and with other models such as HKYg, yield very similar topologies and branch lengths. For instance, using the most complex available model (GTRgi) increases the P. cocincinus branch length over HKYg by a factor of about 1.08 and over JC by about 2.4 (much the same as with the mtDNA analysis above).
FIG. 2.— mtDNA data set (http://Systbiol url, matrix macey-99-59) and tree with inferred branch lengths, all drawn to the same scale, for relative dating of the Australasian-Leiolepis split, Australasian-Physignathus cocincinus split, and the 9-MYA calibration within Laudakia (arrowed). (A) GTRgi model unconstrained; –lnL 48004.05. (B) GTRgi made ultrametric with MLK constraint; –lnL 48211.92. (C) GTRgi made ultrametric with NPRS type 1 (as implemented in TreeEdit). Most taxon names have been omitted to simplify the figure. Ratio of depths of Leiolepis split versus depth of P. cocincinus split: MLK = 2.05, NPRS = 1.82 (thus, NPRS here increases the relative depth of the latter split by a factor of 1.13 compared with MLK). Macey et al. (2000).
The likelihood ratio test identifies the molecular clock model as a worse fit than the unconstrained model (P<<0.001), suggesting the need to account for apparent rate heterogeneity. If the internal calibration (Laudakia divergence discussed above) is fixed at 10 MYA, as proposed by Macey et al. (1998), then the Australasian-Physignathus cocincinus divergence is dated at 20 MYA by type 1 NPRS, 27 MYA by type 3 NPRS, and 28 MYA by MLK. The divergence between Leiolepis and the clade containing the Australasian taxa is dated at 36 MYA by type 1 NPRS, 52 MYA by type 3 NPRS, and 58 MYA by MLK. These dates are reasonably consistent with the fossil record, plate tectonics, and extant distributions (see Discussion). However, if we assume the Australasian-Physignathus cocincinus split is around 150 MYA, as suggested by Schulte, Melville, and Larson (2003), then the divergence between Leiolepis and the Australasian taxa becomes implausibly ancient (274 MYA using type 1 NPRS, 288 MYA using type 3 NPRS, and 307 MYA using MLK), approximately coinciding with the earliest fossil record for reptiles. Macey et al. (2000) acknowledged that uncorrected sequence divergence between Australasian agamids and P. cocincinus is only around twice that between clades of Laudakia separated by an assumed 10 Myr. They reconciled the low divergence between Australasian agamids and P. cocincinus with an early divergence date by suggesting that extensive saturation was causing divergences to level off after approximately 20 Myr, but this is inconsistent with deeper nodes in the data set exhibiting much greater uncorrected sequence divergence. Again, NPRS greatly increases branch lengths and lowers the ratio of deeper to shallower node distances.
Nuclear Data
Nuclear genes are often considered more useful than mitochondrial genes for dating older divergences because of generally slower substitution rates and, thus, a deeper time window before severe saturation occurs (e.g., Glazko and Nei 2003). To provide an independent test of the above issues, we use approximately 950 bp of the nuclear proto-oncogene c-mos for a broad cross section of acrodont (agamid plus chamaeleon) taxa, rooted with an iguanid and scincid (see figure 3 for GenBank accession numbers). Primers were developed and made available by Kathleen M. Saint: primer G303 (coding strand 5' to 3'); AATTATGCCATCMCCTMTTCC, primer G155 (non-coding 5'–3'); TGCTACTAWAGCYYTCCAGC. Concerns about base composition bias in c-mos apply to comparisons involving teiids versus other squamates and are not problematic within agamids (see Harris [2003]).
FIG. 3.— Nuclear c-mos trees with branch lengths, drawn to the same scale. (A) HKYg model unconstrained. –lnL 4515.09. Numbers beside clades denote bootstrap proportions based on 300 pseudoreplicates. GenBank accession numbers below taxon names. (B) HKYg-MLK model with topology and node depth bootstrap variances for selected clades (from 300 pseudoreplicates). –lnL 4530.78. The graph curves indicate the relative frequency of node depths in the bootstraps, for each of the five successive outgroups to the Australasian clade. Ratio of depths of Australasian-Leiolepis split versus depth of Australasian-P. cocincinus split: MLK = 3.16, type 1 NPRS = 2.59 (thus, NPRS increases the relative depth of the latter split by a factor of 1.22 compared with MLK).
Trees and branch lengths were constructed using maximum likelihood (ML) employing the model selected by hierarchical likelihood-ratio tests (as implemented in Modeltest [Posada and Crandall 1998]). The ML tree topology (fig. 3) is very robust. It remains constant across other ML models (HKY-inv to GTR-gamma) and equally weighted parsimony and if relevant mtDNA sequences (from Macey et al. [2000] and Schulte, Melville, and Larson [2003]) are added to the analysis. It is also consistent with the trees derived from a smaller (375-bp) section of c-mos (Harris 2003) and very similar to the mtDNA trees discussed above (Macey et al. 2000; Schulte, Melville, and Larson 2003) (fig. 2), differing only in a minor shift in the position of Leiolepis. Although beyond the scope of this study, we note that the data suggest that chameleons are sister to, rather than nested within, agamids (fide Macey et al. 2000).
There was little stretching of branches using the HKYg or more complex models over the observed differences, suggesting minimal problems with saturation. Limited taxon sampling may hide some saturation among the basal splits; however, the low amount of inferred multiple hits in complex models in the c-mos data suggests this may only amount to a small underestimate. This contrasts with the high inferred sequence divergence and great elongation of basal branches in the mtDNA analyses above, which indicate likely problems with saturation. The observed average pairwise difference across the Australasian-P. cocincinus is 3.5% (uncorrected p-distance), increasing to 4.2% with an HKYg-MLK model and to 6.6% with HKYg-NPRS (type 1) model; much less than for the mtDNA. The likelihood penalty of imposing a molecular clock assumption on the c-mos tree is slight (lnL = 15.7, chi-square test, P = 0.003 with df = 13) and recovers the same tree as the unconstrained search. Given the minor deviation from the molecular clock and low divergence levels, the MLK transformation could well be suitable and even preferable to NPRS (Sanderson 1997, 2002). For this reason, we focus on the MLK analysis, but NPRS dates are very similar (see below). Here, NPRS again increases branch lengths and lowers the ratio of deeper to shallower node distances. This pattern is also apparent in other studies (e.g., Barraclough and Vogler 2001; Martin et al. 2004) and so may be a general effect.
For four well-supported (>90% bootstrap) nodes, variance in branch lengths (and thus divergence dates) was assessed using the nonparametric bootstrap. This involves, for each bootstrap replicate, using the HKYg-MLK model to search for the optimal topology then calculating the depth of the most recent common ancestor of a given set of taxa (Baldwin and Sanderson 1998). Parametric bootstrapping of branch lengths (as employed by Schulte, Melville, and Larson [2003]) gives low variances because it measures only the inherent stochasticity in the chosen substitution model (Felsenstein 2004). In contrast, nonparametric bootstrapping of branch lengths captures errors caused by variation in character sampling and tree topology, which are additional and usually much larger sources of error. The curves showing bootstrap variance in branch lengths are shown in figure 3. Broadly, the 95% confidence interval usually includes 10% to 20% variation in divergence-time estimates, orders of magnitude larger than the standard errors indicated by Schulte, Melville, and Larson (2003) for the mitochondrial data. Even these variance estimates are overly conservative because they ignore such factors as errors in estimating divergence linearity ("stretching") in MLK or NPRS methods and calibration errors (which would uniformly stretch or compress the entire tree). The latter continues to be unjustifiably ignored in molecular clock studies (as emphasized by Graur and Martin [2004]; see also Hedges and Kumar [2004]). However, the magnitude of such fossil calibration errors is difficult to quantify, as it is itself an amalgam of several errors: (1) how closely the first known fossil of a group approaches the (necessarily earlier) true origin of the group, (2) uncertainties in the relative and absolute geochronological dating techniques used to assign an age to the stratum containing the fossil, and (3) uncertainties in the taxonomic assignment of the fossil (i.e., which node or clade the fossil is assumed to "date").
For the c-mos tree, if the Australasian-P. cocincinus split is assumed to be 150 MYA (fide Schulte, Melville, and Larson 2003), the split between Leiolepis and the clade containing the Australasian taxa gets pushed beyond 450 MYA, before the origin of reptiles and, in fact, of terrestrial vertebrates. In particular, the chameleon-agamid split is always at least three times the age of Australasian-P. cocincinus split. These findings are robust to lack of linearity and, hence, rule out any Gondwanan scenario.
There are no well-corroborated calibration points within the agamids for the c-mos tree, but there are two fossil dates to provide benchmarks for the most basal (acrodont-iguanid) split, with the caveat of assuming reasonable linearity. The earliest taxon that can be provisionally assigned to either side of this split is a likely stem acrodont from India, dated at around 190 MYA (Evans, Prasad, and Manhas 2002). However, this taxon is known only from jaw fragments, and the first unequivocal iguanians (priscagamids and Pristiguana) only occur much later (e.g., Evans 1993). The priscagamids have been rigorously shown to be stem acrodonts (Frost and Etheridge 1989) and first occur in Aptian-Albian deposits of Asia, around 110 MYA (Evans, Prasad, and Manhas 2002, and references therein). The phylogenetic relationships of other contemporaneous iguanians have not been properly analyzed, and they cannot yet be used for calibration. If the tentative, earlier date is used to calibrate the tree, the Australasian-P. cocincinus split is dated at 33 MYA (95% CI = 24–41 MYA), whereas the Australasian-Leiolepis split is dated at 106 MYA (95% CI = 90–126 MYA). If the later calibration date is used as an absolute minimum boundary, these divergences reduced to 19 MYA (95% CI = 14–24 MYA) and 62 MYA (95% CI = 52–73 MYA), respectively. For the c-mos tree with nonparametric rate smoothing (type 3), the two splits are dated at 50 MYA and 119 MYA (using the 190 MYA calibration) or at 29 MYA and 69 MYA (using the 110 MYA calibration). These dates are comparable with the revised mitochondrial DNA dates calculated above (MLK dates of 28 and 58 MYA, and NPRS type 3 dates of 27 and 52 MYA). Note that the Leiolepis split is approximately three times as old as the P. cocincinus split in the nuclear trees, compared with only twice as old in the mtDNA trees. The compression of the basal nodes in the mtDNA versus the c-mos trees suggests greater saturation effects in the mtDNA.
Discussion
Macey et al. (2000) used a mtDNA phylogeny of agamids to propose a Gondwanan age (>120 MYA) for the split between the SE Asian Physignathus cocincinus and the Australasian agamids, based on the monophyly of diverse Australasian radiation and citing suggestions that taxa could not have dispersed into Australasia before 10 MYA (Hall 1996; Metcalfe 1996). The incongruously close genetic distances between Physignathus cocincinus and the Australasian agamids was explained as an artifact of severe mtDNA saturation. Schulte, Melville, and Larson (2003), using a similar data set, estimated the age of this split at 150 MYA. These proposals have not yet been seriously questioned and have been cited in support of hypotheses from early fossils recently interpreted as agamids (Evans, Prasad, and Manhas 2002). However, these analyses are suspect because they use rates of molecular evolution inferred from uncorrected sequence divergence to put a time frame on a tree with branch lengths greatly elongated by complex likelihood and rate-smoothing models. Although we have not reanalyzed the varanid data of Schulte, Melville, and Larson (2003) in detail, we note their conclusion of deep (>100 MYA) divergences within the genus is affected by exactly the same methodological bias.
The more appropriate analyses of the agamid data here demonstrate that both the mitochondrial and nuclear sequence data imply a relatively recent divergence between Australasian agamids and the SE Asian P. cocincinus (between 14 and 41 MYA). These dates are sufficiently recent to exclude the Gondwanan vicariance hypothesis (Macey et al. 2000; Schulte, Melville, and Larson 2003), which suggested the split between Australasian agamids and P. cocincinus occurred before 120 MYA, when the last terrains split off from Gondwana and drifted north to Eurasia (Metcalfe 2001). Furthermore, for both the mtDNA and c-mos data sets, the proposed approximately 150-MYA split between P. cocincinus and Australasian agamids (Schulte, Melville, and Larson 2003; see also Macey et al. [2000]) implies implausibly ancient divergences in basal living agamids (300–400 MYA). These inferred divergence dates predate the fossil record for lizards and are as old as, or older than, the earliest fossil record for reptiles (Lee 1999; Graur and Martin 2004). Similary, their inferred dates (>110 MYA) for splits within Varanus imply more basal splits happened implausibly early (e.g., 300 MYA for the Lanthanotus-Varanus split, again coincident with the earliest fossil reptiles). Any proposed dates for particular nodes should be strongly suspect when they imply such unlikely dates for other splits within the tree.
Contrary to Schulte, Melville, and Larson (2003), the recent discovery of a 190-Myr-old putative acrodont (Evans, Prasad, and Manhas 2002) does not provide much support for their proposed dates. Using this very early fossil to calibrate the iguanid tree (fig. 3, c-mos) still only yields a divergence date for P. cocincinus and the Australasian agamids of 33 MYA, much more recent than the predicted 150 MYA. Finding a stem agamid that is 700 Myr old would approach the latter result; however, despite the problem with quantifying calibration uncertainty (see above), it seems safe to say such a scenario is unlikely. The "deep divergence" hypothesis also implies that every major agamid clade arose before the formation of Pangea in the Late Permian (250 MYA) and that multiple lineages from each clade should have been present on most continents. It thus has difficulty explaining the total absence of fossil and living agamids in both North and South America and the absence of agamid fossils in all continents before 210 MYA, even those regions with an excellent Permo-Triassic record (e.g., southern Africa and Russia). Finally, one might attempt to recast the Gondwanan vicariance model by assuming that the combined P. cocincinus-Australasian clade (rather than the Australasian clade alone) arose by vicariance, with P. cocincinus secondarily colonizing Eurasia. However, even this more inclusive clade appears too young (65 to 112 MYA using iguanid calibration bounds) to have formed through the early breakup of Gondwana, which requires it to be over 120 Myr old (see Schulte, Melville, and Larson [2003]).
Both the mtDNA and nuclear data indicate that the extant Australasian agamid radiation is between 14 and 41 Myr old (i.e., Oligo-Miocene). This age is consistent with the time frames for dispersal between Asia and Australasia suggested by plate tectonics (Cogger and Heatwole 1981). The gradual northward drift of the Australasian plate does not provide absolute constraints on the timing of dispersals (Hutchinson and Donnellan 1993; contra Hall [1996]). All that can be assumed is that during the Tertiary, dispersal between Asia and Australasia became progressively more likely with time as the plates approached one another. The latest studies suggest that short-range, island-hopping dispersal between Asia and Australasia could have occurred at least as long ago as 30 Myr (Metcalfe 2001), and longer-range dispersal was possible much earlier. The arboreal and/or amphibious habits of many acrodonts (including P. cocincinus and many basal Australasian agamids) would make them reasonable candidates for early dispersal through rafting. Schulte, Melville, and Larson (2003) suggested that agamids were not adept over-water dispersers; however, the three clades of iguanians (chamaeleons, agamids, and iguanids) all appear to contain forms that readily cross water barriers. Chameleons have been inferred to have crossed repeatedly between Madagascar and Africa (Raxworthy, Forstner, and Nussbaum 2002), and agamids and iguanids are found on numerous islands that have always been separated from the mainland by large water gaps (e.g., Moody 1980).
The dates estimated here are also compatible with the presence of agamids in Africa, Asia, and Australasia and their absence from the Americas. This is consistent with the relatively recent origin of agamids in Eurasia (where they potentially have their earliest fossil record at 190 MYA [Evans, Prasad, and Manhas 2002]), followed by dispersal into the Gondwanan landmasses of Africa, India, and Australasia when they drifted northwards. Since the breakup of Pangea, South America has never been connected to Eurasia, and North America connected only at high latitudes. The dates inferred here for basal divergences in agamids are also broadly consistent with the fossil record, placing the agamid radiation within the past 200 Myr, matching paleontological estimates of the timing of the lizard radiation (Estes 1983; Evans 1993).
Within agamids, the mtDNA data show large divergences, evidence of saturation effects, and considerable sensitivity to details of substitution model and rate-smoothing methods. In contrast, the large portion of c-mos appears less affected, performs well at retrieving basal nodes and giving stable estimates for deeper branch lengths in iguanians and should be a useful locus for other similar studies.
Acknowledgements
We thank Mark Hutchinson, Steve Donnellan, Remko Leijs, and two anonymous referees for discussion and comments; Kathleen M. Saint, Steve Donnellan and Mark Hutchinson for generously allowing access to their primer information and c-mos sequences; and the Australian Research Council for funding.
References
Audley-Charles, M. G. 1987. Dispersal of Gondwanaland: relevance to evolution of the angiosperms. Pp. 5–25 in T. C. Whitmore, ed. Biogeographical evolution of the Malay Archipelago. Oxford Monographs on Biogeography No. 4. Clarendon Press, Oxford.
Benton, M. J., and F. J. Ayala. 2003. Dating the tree of life. Science 300:1698–1700.
Baldwin, B. G., and M. J. Sanderson. 1998. Age and rate of diversification of the Hawaiian silversword alliance (Compositae). Proc. Nat. Acad. Sci. USA 95:9402–9406.
Barraclough, T. G., and A. P. Vogler. 2002. Recent diversification rates in North American tiger beetles estimated from a dated mtDNA phylogenetic tree. Mol. Biol. Evol. 19:1706–1716.
Cogger, H. G., and H. Heatwole. 1981. The Australian reptiles: origins, biogeography, distribution patterns and island evolution. Pp. 1333–1373 in A. Keast, ed. Ecological biogeography of Australia. W. Junk, The Hague.
Estes, R. 1983. The fossil record and early distribution of lizards. Pp. 365–398 in A. G. J. Rhodin and K. Miyata, eds. Advances in herpetology and evolutionary biology: essays in honour of Ernest E. Williams. Museum of Comparative Zoology, Cambridge, Mass.
Evans, S. E. 1993. Jurassic lizard assemblages. Revue de Paléobiologie (Suisse), Volume spéciale 7:55–65.
Evans, S. E., G. V. R. Prasad, and B. K. Manhas. 2002. Fossil lizards from the Jurassic Kota Formation of India. J. Vert. Paleont. 22:299–312.
Glazko, G. V., and M. Nei. 2003. Estimation of divergence time for major lineages of primate species. Mol. Biol. Evol. 20:424–434.
Graur, D., and W. Martin. 2004. Reading the entrails of chickens: molecular timescales of evolution and the illusion of precision. Trends Genet. 20:80–86.
Greer, A. E. 1989. The biology and evolution of Australian lizards. Surrey Beatty, Chipping Norton.
Hall, R. 1996. Reconstructing Cenozoic SE Asia. Pp. 153–184 in R. Hall and D. Blundell, eds. Tectonic evolution of Southeast Asia. Geological Society, London.
———. 1998. The plate tectonics of Cenozoic SE Asia and the distribution of land and sea. Pp. 99–131 in R. Hall and J. D. Holloway, eds. Biogeography and geological evolution of SE Asia. Backhuys, Leiden.
Harris, D. J. 2003 Codon bias variation in C-mos between squamate families might distort phylogenetic inferences. Mol. Phylogenet. Evol. 27:540–544.
Hedges, S. B., and S. Kumar. 2004. Precision of molecular time estimates. Trends Genet. 20:242–247.
Hutchinson, M. N., and S. C. Donnellan. 1993 Biogeography and phylogeny of the Squamata. Pp. 210–220 in C. J. Glasby, G. J. B. Ross, and P. L. Beesley, eds. Fauna of Australia, Vol 2A: Amphibia and Reptilia. AGPS Press, Canberra.
Felsenstein, J. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368–376.
———. 2004. Inferring phylogenies. Sinauer, Sunderland, Mass.
Lee, M.S.Y. 1999. Molecular clock calibrations and metazoan divergence dates. J. Mol. Evol. 49:385–391.
Macey, J. R., J. A. Schulte II, N. B. Ananjeva, A. Larson, N. Rastegar-Pouyani, M. Shammakov, and T. Papenfuss. 1998. Phylogenetic relationships among agamid lizards of the Laudakia caucasia species group; testing hypotheses of biogeographic fragmentation and an area cladogram for the Iranian plateau. Mol. Phylogenet. Evol. 10:118–131.
Macey, J. R., J. A. Schulte II, A. Larson, N. B. Ananjeva, Y. Wang, N. R. Pethiyagoda, N. Rastegar-Pouyani, and T. J. Papenfuss. 2000. Evaluating trans-Tethys migration: an example using Acrodont lizard phylogenetics. Syst. Biol. 49:233–256.
Martin, A. P., E. K. Costello, A. F. Meyer, D. R. Nemergut, and S. K. Schmidt. 2004. The rate and pattern of cladogenesis in microbes. Evolution 58:946–955.
Metcalfe, I. 1996. Precretaceous evolution of SE Asian terranes. Pp. 97–122 in R. Hall and D. Blundell, eds. Tectonic evolution of Southeast Asia. Geological Society, London.
———. 2001. Palaeozoic and Mesozoic tectonic evolution and biogeography of SE Asia-Australasia. Pp. 15–56 in I. Metcalfe, ed. Faunal and floral migrations and evolution in SE Asia-Australasia. Balkema Publishers, Exton, Pa.
Moody, S. 1980. Phylogenetic and historical biogeographical relationships of the genera in the Agamidae (Reptilia, Lacertilia). Doctoral dissertation, University of Michigan, Ann Arbor.
Posada, D., and K. Crandall, K. 1998. Modeltest: testing the model of DNA substitution. Bioinformatics 14:817–818.
Rambaut, A., and L. Bromham. 1998. Estimating divergence dates from molecular sequences. Mol. Biol. Evol. 15:442–448.
Rambaut, A., and M. Charleston. 2001. TreeEdit: phylogenetic tree editor. Version 1.0. Department of Zoology, University of Oxford. URL http://evolve.zoo.ox.ac.uk/software/TreeEdit/main.html
Raxworthy, C. J., M. R. J. Forstner, and R. Nussbaum. 2002. Chameleon radiation by oceanic dispersal. Nature 415:784–787.
Sanderson, M. J. 1997. A nonparametric approach to estimating divergence times in the absence of rate constancy. Mol. Biol. Evol. 14:1218–1231.
———. 2002. Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Mol. Biol. Evol. 19:101–109.
Schulte, J. A., J. Melville, and A. Larson. 2003. Molecular phylogenetic evidence for ancient divergence of lizard taxa on either side of Wallace's Line. Proc. R. Soc. Lond. B Biol. Sci. 270:597–603.
Schwaner, T. D., P. R. Baverstock, H. C. Dessauer, and G. A. Mengden. 1985 Immunological evidence for the phylogenetic relationships of Australian elapid snakes. Pp. 177–184 in G. C. Grigg, R. Shine, and H. Ehmann, eds. The biology of Australasian frogs and reptiles. Royal Zoological Society of NSW, Sydney.
Swofford, D. L. 2002. PAUP*: phylogenetic analysis using parsimony (*and other methods). Version 4. Sinauer Associates, Sunderland, Mass.
Weisrock, D. W., R. J. Macey, I. H. Ugurtas, A. Larson, and T. J. Papenfuss. 2001. Molecular phylogenetics and historical biogeography among Salamadrids of the ‘True’ Salamader clade: rapid branching of numerous highly divergent lineages in Mertensiella luschani associated with the rise of Anatolia. Mol. Phylogenet. Evol. 18:434–448.
Woodburne, M. O., and J. A. Case. 1996. Dispersal, vicariance, and the late Cretaceous to early Tertiary land mammal biogeography from South America to Australia. J. Mamm. Evol. 3:121–161.(Andrew F. Hugall* and Mic)
E-mail: Andrew.Hugall@adelaide.edu.au.
Abstract
A recent mtDNA study proposes a surprisingly deep (150 MYA) divergence between SE Asian and Australasian agamid lizards, consistent with ancient Gondwanan vicariance rather than dispersal across the Indonesian Archipelago. However, the analysis contains a fundamental error: use of rates of molecular evolution inferred from uncorrected sequence divergence to put a time frame on a tree with branch lengths greatly elongated by complex likelihood and rate-smoothing models. Furthermore, this date implies that basal splits within agamids occurred implausibly early, at least 300 MYA (100 Myr before the first fossil lizards and coincident with the earliest fossil reptiles). Analyses of the mtDNA data using more appropriate methods and new information from nuclear (c-mos) sequences suggest a much more recent divergence between SE Asian and Australian agamids (around 30 MYA). Using two fossil boundary dates, bootstrapping the c-mos data gives a 95% confidence interval for this divergence time that is sufficiently recent (14–41 MYA) to exclude an ancient Gondwanan vicariance and is more consistent with Miocene over-water dispersal. As with the mtDNA, the c-mos data implies implausibly old basal divergences among agamids if a Gondwanan age is assumed for the Australasian clade. The analyses also highlight how methods for creating ultrametric trees (especially nonparametric rate smoothing) can greatly modify branch lengths and, thus, always require internal calibrations. The errors associated with inferred dates in the previous study (inferred through parametric bootstrapping) were also unjustifiably low, as this method only considers stochasticity in the substitution model and ignores much larger sources of uncertainty, such as variation in character sampling, tree topology, and calibration accuracy.
Key Words: molecular clock ? rate-smoothing ? ultrametric trees ? Agamidae ? Wallace line ? Gondwanan biogeography
Introduction
Most of the extant Australasian biota is widely accepted to have originated from two ancestral sources. Australasia had land connections with other Gondwanan continents (Antarctica and South America) until the early Paleocene, 55 MYA at the latest (Audley-Charles 1987; Woodburne and Case 1996; Metcalfe 2001). There followed a period of isolation as the Australasian plate detached from Antarctica and drifted northwards; the isolation ended perhaps as early as 30 MYA, as the tectonic collision with Southeast Asia formed a filter-bridge allowing colonization of terrestrial organisms from the north (Hall 1998; but see Discussion). The hiatus defines a dichotomy between two distinct components of the Australasian biota: (1) ancient lineages that were distributed across Gondwana and became isolated in Australasia and (2) recent lineages that represent invasions from the north, when Australia collided with the Asian plate. Clearly, for any endemic Australasian radiation, these two routes predict very different distributions regarding nearest extralimital relatives (South America versus Asia) and very different divergence times from those relatives (55 MYA versus 30 MYA).
The 64 or more species of agamid lizards ("dragons") constitute a sizeable component (13.5%) of the Australian lizard fauna and traditionally have been interpreted to have arrived by dispersal from the north (e.g., Cogger and Heatwole 1981; Schwaner et al. 1985; Greer 1989). This interpretation, however, has been challenged by a recent analysis of mtDNA data that inferred deep (150 MYA) divergence dates between Australasian agamids and their sister group (the SE Asian mainland Physignathus cocincinus) that could only be consistent with an ancient vicariance during the early breakup of Gondwana (Schulte, Melville, and Larson 2003; see also Moody [1980] and Macey et al. [2000]).
Schulte, Melville, and Larson (2003) used an approximately 1,600-bp segment of mtDNA to construct the phylogeny of Australasian and a sparser sampling of other agamids. Branch lengths were reconstructed using a parameter-rich maximum-likelihood (ML) model and then, because of significant rate variation, further modified using nonparametric rate smoothing (NPRS) (Sanderson 1997) to obtain an ultrametric tree. These reconstructed branch lengths were translated into absolute time using a calibration rate of 0.65% change per Myr along each lineage, calculated using uncorrected differences in Laudakia agamids (Macey et al. 1998), and also consistent with other ectothermic vertebrates (Weisrock et al. 2001). The divergence between the Australasian radiation and P. cocincinus was thus dated at around 150 MYA. The problem here is that a tree with highly "stretched" branches constructed under one model has been assigned absolute dates using a calibration rate inferred from (and applicable to) a tree with much shorter branches reconstructed under a very different model.
Here, we reanalyze the mtDNA data of Schulte, Melville, and Larson (2003) using more appropriate methods and taxon sampling. The results all suggest divergence dates around four to five times shallower (30 MYA) than proposed, with associated errors that are much greater than proposed. This more recent divergence is corroborated by analysis of new nuclear gene (c-mos) sequences.
Mitochondrial Data
We repeated the analysis of Schulte, Melville, and Larson (2003), firstly using their data matrix (TreeBase matrix N1326), tree, and model (see figure 1) in PAUP* (Swofford 2002). The GTR-gamma-inv (GTRgi) model adopted stretches the observed branch length between the Australasian radiation and P. cocincinus from 19% (uncorrected) to around 42%. NPRS as implemented in TreeEdit (Rambaut and Charleston 2001) further stretches the branch length to 88% divergence (type 1 method, "weight rate difference across root") or 95% divergence (type 3 method, "weight rate differences at all nodes with mean"). That is, the original uncorrected difference has been doubled by the GTRgi model and doubled again by NPRS. These greatly extended branches (inflated divergences) are then translated into time using an uncorrected divergence rate of 0.65% per Myr per lineage to give 150 Myr (see figure 1).
FIG. 1.— Agamid mtDNA data set N1326 and trees with inferred branch lengths, all drawn to the same scale, for relative dating of the Australasian-Leiolepis and Australasian-Physignathus cocincinus splits (arrowed). (A) Jukes-Cantor; (B) Jukes-Cantor made ultrametric with MLK constraint; (C) GTRgi model with MLK constraint; (D) GTRgi made ultrametric with NPRS type 3 (as implemented in TreeEdit). Most taxon names have been omitted to simplify the figure. Ratio of depths of Leiolepis split versus depth of P. cocincinus split: MLK = 2.16, NPRS = 1.52 (thus, NPRS increases the relative depth of the latter split by a factor of 1.42 compared with MLK). Schulte, Melville, and Larson (2003).
However, using the same calibration rate on the rest of the tree dates the dichotomy between the clade containing Leiolepis and the clade containing the Australasian taxa at around 250 MYA (227 MYA using type 1 NPRS, 254 MYA using type 3 NPRS, and 310 MYA using MLK constraint instead of NPRS ). This is more than 100 Myr before the oldest unequivocal fossil agamid, much older than even the 190-Myr putative acrodont described by Evans, Prasad, and Manhas (2002), and very soon after the first unequivocal reptiles (around 288 MYA [Lee 1999]). It is now a common procedure to test the fit of the MLK clock constraint using a likelihood ratio test (e.g., Felsenstein 1981) and, if significantly worse, to make the tree ultrametric using NPRS. In the face of substantial accurately measurable rate heterogeneity, NPRS may be preferred to MLK for the ultrametric transformation (Sanderson 1997, 2002); however, attention needs to be paid to the effect it has on the shape and size of the tree. In the above analysis, NPRS greatly increases the overall length of the tree and also alters the ratio of the depths of the nodes in question, compared with the effect of the MLK transformation. Here, NPRS lowers the ratio of deeper (Leiolepis) to shallower (P. cocincinus) node distances, "fattening" the shape of the tree (see figure legends).
The results of Schulte, Melville, and Larson (2003) are erroneous because an ML+NPRS tree with greatly lengthened, reconstructed branch lengths has been related to absolute time using a divergence rate inferred in other studies using uncorrected sequence divergence. The obvious way (Rambaut and Bromham 1998; Benton and Ayala 2003) to avoid this error is to employ internal calibration(s), so that branch lengths (and inferred rates) subtending the calibration node(s) are modified in exactly the same fashion as branch lengths (and inferred rates) throughout the rest of the tree. Macey et al. (2000) analyzed agamid phylogeny using the same region as Schulte, Melville, and Larson (2003). The data set of Macey et al. (2000) is more appropriate to testing the questions posed in Schulte, Melville, and Larson (2003) than the data set in the latter study for two reasons: (1) taxon sampling is much more balanced, with many more non-Australasian agamids included, and (2) there are multiple geologically dated internal calibration points within Laudakia that are all highly consistent with each other and broadly consistent with other studies (Macey et al. 1998; Schulte, Melville, and Larson 2003). The divergence used to calibrate the present tree is the isolation of the Iranian plateau Laudakia (L. microlepis, L. erythrogastra, and L. caucasia) from their nearest northern relatives (identified as the L. himalayana clade by Macey et al. [2000]), which Macey et al. (1998) suggested was associated with a tectonic event of 10 MYA (the Pamir-TienShan uplift). Similarly, use of all the possible Laudakia calibration points simultaneously (e.g., using penalized-likelihood rate smoothing [Sanderson 2002]) would also give virtually identical results, given they are located on the same region of the tree and highly consistent with each other.
Using Macey et al.'s (2000) alignment (http://Systbiol url, matrix macey-99-59) and tree (their figure 8), with the GTRgi model, we employ the same procedure—unconstrained GTRgi tree made ultrametric by (1) MLK and (2) NPRS (fig. 2). Analyses with MrBayes, and with other models such as HKYg, yield very similar topologies and branch lengths. For instance, using the most complex available model (GTRgi) increases the P. cocincinus branch length over HKYg by a factor of about 1.08 and over JC by about 2.4 (much the same as with the mtDNA analysis above).
FIG. 2.— mtDNA data set (http://Systbiol url, matrix macey-99-59) and tree with inferred branch lengths, all drawn to the same scale, for relative dating of the Australasian-Leiolepis split, Australasian-Physignathus cocincinus split, and the 9-MYA calibration within Laudakia (arrowed). (A) GTRgi model unconstrained; –lnL 48004.05. (B) GTRgi made ultrametric with MLK constraint; –lnL 48211.92. (C) GTRgi made ultrametric with NPRS type 1 (as implemented in TreeEdit). Most taxon names have been omitted to simplify the figure. Ratio of depths of Leiolepis split versus depth of P. cocincinus split: MLK = 2.05, NPRS = 1.82 (thus, NPRS here increases the relative depth of the latter split by a factor of 1.13 compared with MLK). Macey et al. (2000).
The likelihood ratio test identifies the molecular clock model as a worse fit than the unconstrained model (P<<0.001), suggesting the need to account for apparent rate heterogeneity. If the internal calibration (Laudakia divergence discussed above) is fixed at 10 MYA, as proposed by Macey et al. (1998), then the Australasian-Physignathus cocincinus divergence is dated at 20 MYA by type 1 NPRS, 27 MYA by type 3 NPRS, and 28 MYA by MLK. The divergence between Leiolepis and the clade containing the Australasian taxa is dated at 36 MYA by type 1 NPRS, 52 MYA by type 3 NPRS, and 58 MYA by MLK. These dates are reasonably consistent with the fossil record, plate tectonics, and extant distributions (see Discussion). However, if we assume the Australasian-Physignathus cocincinus split is around 150 MYA, as suggested by Schulte, Melville, and Larson (2003), then the divergence between Leiolepis and the Australasian taxa becomes implausibly ancient (274 MYA using type 1 NPRS, 288 MYA using type 3 NPRS, and 307 MYA using MLK), approximately coinciding with the earliest fossil record for reptiles. Macey et al. (2000) acknowledged that uncorrected sequence divergence between Australasian agamids and P. cocincinus is only around twice that between clades of Laudakia separated by an assumed 10 Myr. They reconciled the low divergence between Australasian agamids and P. cocincinus with an early divergence date by suggesting that extensive saturation was causing divergences to level off after approximately 20 Myr, but this is inconsistent with deeper nodes in the data set exhibiting much greater uncorrected sequence divergence. Again, NPRS greatly increases branch lengths and lowers the ratio of deeper to shallower node distances.
Nuclear Data
Nuclear genes are often considered more useful than mitochondrial genes for dating older divergences because of generally slower substitution rates and, thus, a deeper time window before severe saturation occurs (e.g., Glazko and Nei 2003). To provide an independent test of the above issues, we use approximately 950 bp of the nuclear proto-oncogene c-mos for a broad cross section of acrodont (agamid plus chamaeleon) taxa, rooted with an iguanid and scincid (see figure 3 for GenBank accession numbers). Primers were developed and made available by Kathleen M. Saint: primer G303 (coding strand 5' to 3'); AATTATGCCATCMCCTMTTCC, primer G155 (non-coding 5'–3'); TGCTACTAWAGCYYTCCAGC. Concerns about base composition bias in c-mos apply to comparisons involving teiids versus other squamates and are not problematic within agamids (see Harris [2003]).
FIG. 3.— Nuclear c-mos trees with branch lengths, drawn to the same scale. (A) HKYg model unconstrained. –lnL 4515.09. Numbers beside clades denote bootstrap proportions based on 300 pseudoreplicates. GenBank accession numbers below taxon names. (B) HKYg-MLK model with topology and node depth bootstrap variances for selected clades (from 300 pseudoreplicates). –lnL 4530.78. The graph curves indicate the relative frequency of node depths in the bootstraps, for each of the five successive outgroups to the Australasian clade. Ratio of depths of Australasian-Leiolepis split versus depth of Australasian-P. cocincinus split: MLK = 3.16, type 1 NPRS = 2.59 (thus, NPRS increases the relative depth of the latter split by a factor of 1.22 compared with MLK).
Trees and branch lengths were constructed using maximum likelihood (ML) employing the model selected by hierarchical likelihood-ratio tests (as implemented in Modeltest [Posada and Crandall 1998]). The ML tree topology (fig. 3) is very robust. It remains constant across other ML models (HKY-inv to GTR-gamma) and equally weighted parsimony and if relevant mtDNA sequences (from Macey et al. [2000] and Schulte, Melville, and Larson [2003]) are added to the analysis. It is also consistent with the trees derived from a smaller (375-bp) section of c-mos (Harris 2003) and very similar to the mtDNA trees discussed above (Macey et al. 2000; Schulte, Melville, and Larson 2003) (fig. 2), differing only in a minor shift in the position of Leiolepis. Although beyond the scope of this study, we note that the data suggest that chameleons are sister to, rather than nested within, agamids (fide Macey et al. 2000).
There was little stretching of branches using the HKYg or more complex models over the observed differences, suggesting minimal problems with saturation. Limited taxon sampling may hide some saturation among the basal splits; however, the low amount of inferred multiple hits in complex models in the c-mos data suggests this may only amount to a small underestimate. This contrasts with the high inferred sequence divergence and great elongation of basal branches in the mtDNA analyses above, which indicate likely problems with saturation. The observed average pairwise difference across the Australasian-P. cocincinus is 3.5% (uncorrected p-distance), increasing to 4.2% with an HKYg-MLK model and to 6.6% with HKYg-NPRS (type 1) model; much less than for the mtDNA. The likelihood penalty of imposing a molecular clock assumption on the c-mos tree is slight (lnL = 15.7, chi-square test, P = 0.003 with df = 13) and recovers the same tree as the unconstrained search. Given the minor deviation from the molecular clock and low divergence levels, the MLK transformation could well be suitable and even preferable to NPRS (Sanderson 1997, 2002). For this reason, we focus on the MLK analysis, but NPRS dates are very similar (see below). Here, NPRS again increases branch lengths and lowers the ratio of deeper to shallower node distances. This pattern is also apparent in other studies (e.g., Barraclough and Vogler 2001; Martin et al. 2004) and so may be a general effect.
For four well-supported (>90% bootstrap) nodes, variance in branch lengths (and thus divergence dates) was assessed using the nonparametric bootstrap. This involves, for each bootstrap replicate, using the HKYg-MLK model to search for the optimal topology then calculating the depth of the most recent common ancestor of a given set of taxa (Baldwin and Sanderson 1998). Parametric bootstrapping of branch lengths (as employed by Schulte, Melville, and Larson [2003]) gives low variances because it measures only the inherent stochasticity in the chosen substitution model (Felsenstein 2004). In contrast, nonparametric bootstrapping of branch lengths captures errors caused by variation in character sampling and tree topology, which are additional and usually much larger sources of error. The curves showing bootstrap variance in branch lengths are shown in figure 3. Broadly, the 95% confidence interval usually includes 10% to 20% variation in divergence-time estimates, orders of magnitude larger than the standard errors indicated by Schulte, Melville, and Larson (2003) for the mitochondrial data. Even these variance estimates are overly conservative because they ignore such factors as errors in estimating divergence linearity ("stretching") in MLK or NPRS methods and calibration errors (which would uniformly stretch or compress the entire tree). The latter continues to be unjustifiably ignored in molecular clock studies (as emphasized by Graur and Martin [2004]; see also Hedges and Kumar [2004]). However, the magnitude of such fossil calibration errors is difficult to quantify, as it is itself an amalgam of several errors: (1) how closely the first known fossil of a group approaches the (necessarily earlier) true origin of the group, (2) uncertainties in the relative and absolute geochronological dating techniques used to assign an age to the stratum containing the fossil, and (3) uncertainties in the taxonomic assignment of the fossil (i.e., which node or clade the fossil is assumed to "date").
For the c-mos tree, if the Australasian-P. cocincinus split is assumed to be 150 MYA (fide Schulte, Melville, and Larson 2003), the split between Leiolepis and the clade containing the Australasian taxa gets pushed beyond 450 MYA, before the origin of reptiles and, in fact, of terrestrial vertebrates. In particular, the chameleon-agamid split is always at least three times the age of Australasian-P. cocincinus split. These findings are robust to lack of linearity and, hence, rule out any Gondwanan scenario.
There are no well-corroborated calibration points within the agamids for the c-mos tree, but there are two fossil dates to provide benchmarks for the most basal (acrodont-iguanid) split, with the caveat of assuming reasonable linearity. The earliest taxon that can be provisionally assigned to either side of this split is a likely stem acrodont from India, dated at around 190 MYA (Evans, Prasad, and Manhas 2002). However, this taxon is known only from jaw fragments, and the first unequivocal iguanians (priscagamids and Pristiguana) only occur much later (e.g., Evans 1993). The priscagamids have been rigorously shown to be stem acrodonts (Frost and Etheridge 1989) and first occur in Aptian-Albian deposits of Asia, around 110 MYA (Evans, Prasad, and Manhas 2002, and references therein). The phylogenetic relationships of other contemporaneous iguanians have not been properly analyzed, and they cannot yet be used for calibration. If the tentative, earlier date is used to calibrate the tree, the Australasian-P. cocincinus split is dated at 33 MYA (95% CI = 24–41 MYA), whereas the Australasian-Leiolepis split is dated at 106 MYA (95% CI = 90–126 MYA). If the later calibration date is used as an absolute minimum boundary, these divergences reduced to 19 MYA (95% CI = 14–24 MYA) and 62 MYA (95% CI = 52–73 MYA), respectively. For the c-mos tree with nonparametric rate smoothing (type 3), the two splits are dated at 50 MYA and 119 MYA (using the 190 MYA calibration) or at 29 MYA and 69 MYA (using the 110 MYA calibration). These dates are comparable with the revised mitochondrial DNA dates calculated above (MLK dates of 28 and 58 MYA, and NPRS type 3 dates of 27 and 52 MYA). Note that the Leiolepis split is approximately three times as old as the P. cocincinus split in the nuclear trees, compared with only twice as old in the mtDNA trees. The compression of the basal nodes in the mtDNA versus the c-mos trees suggests greater saturation effects in the mtDNA.
Discussion
Macey et al. (2000) used a mtDNA phylogeny of agamids to propose a Gondwanan age (>120 MYA) for the split between the SE Asian Physignathus cocincinus and the Australasian agamids, based on the monophyly of diverse Australasian radiation and citing suggestions that taxa could not have dispersed into Australasia before 10 MYA (Hall 1996; Metcalfe 1996). The incongruously close genetic distances between Physignathus cocincinus and the Australasian agamids was explained as an artifact of severe mtDNA saturation. Schulte, Melville, and Larson (2003), using a similar data set, estimated the age of this split at 150 MYA. These proposals have not yet been seriously questioned and have been cited in support of hypotheses from early fossils recently interpreted as agamids (Evans, Prasad, and Manhas 2002). However, these analyses are suspect because they use rates of molecular evolution inferred from uncorrected sequence divergence to put a time frame on a tree with branch lengths greatly elongated by complex likelihood and rate-smoothing models. Although we have not reanalyzed the varanid data of Schulte, Melville, and Larson (2003) in detail, we note their conclusion of deep (>100 MYA) divergences within the genus is affected by exactly the same methodological bias.
The more appropriate analyses of the agamid data here demonstrate that both the mitochondrial and nuclear sequence data imply a relatively recent divergence between Australasian agamids and the SE Asian P. cocincinus (between 14 and 41 MYA). These dates are sufficiently recent to exclude the Gondwanan vicariance hypothesis (Macey et al. 2000; Schulte, Melville, and Larson 2003), which suggested the split between Australasian agamids and P. cocincinus occurred before 120 MYA, when the last terrains split off from Gondwana and drifted north to Eurasia (Metcalfe 2001). Furthermore, for both the mtDNA and c-mos data sets, the proposed approximately 150-MYA split between P. cocincinus and Australasian agamids (Schulte, Melville, and Larson 2003; see also Macey et al. [2000]) implies implausibly ancient divergences in basal living agamids (300–400 MYA). These inferred divergence dates predate the fossil record for lizards and are as old as, or older than, the earliest fossil record for reptiles (Lee 1999; Graur and Martin 2004). Similary, their inferred dates (>110 MYA) for splits within Varanus imply more basal splits happened implausibly early (e.g., 300 MYA for the Lanthanotus-Varanus split, again coincident with the earliest fossil reptiles). Any proposed dates for particular nodes should be strongly suspect when they imply such unlikely dates for other splits within the tree.
Contrary to Schulte, Melville, and Larson (2003), the recent discovery of a 190-Myr-old putative acrodont (Evans, Prasad, and Manhas 2002) does not provide much support for their proposed dates. Using this very early fossil to calibrate the iguanid tree (fig. 3, c-mos) still only yields a divergence date for P. cocincinus and the Australasian agamids of 33 MYA, much more recent than the predicted 150 MYA. Finding a stem agamid that is 700 Myr old would approach the latter result; however, despite the problem with quantifying calibration uncertainty (see above), it seems safe to say such a scenario is unlikely. The "deep divergence" hypothesis also implies that every major agamid clade arose before the formation of Pangea in the Late Permian (250 MYA) and that multiple lineages from each clade should have been present on most continents. It thus has difficulty explaining the total absence of fossil and living agamids in both North and South America and the absence of agamid fossils in all continents before 210 MYA, even those regions with an excellent Permo-Triassic record (e.g., southern Africa and Russia). Finally, one might attempt to recast the Gondwanan vicariance model by assuming that the combined P. cocincinus-Australasian clade (rather than the Australasian clade alone) arose by vicariance, with P. cocincinus secondarily colonizing Eurasia. However, even this more inclusive clade appears too young (65 to 112 MYA using iguanid calibration bounds) to have formed through the early breakup of Gondwana, which requires it to be over 120 Myr old (see Schulte, Melville, and Larson [2003]).
Both the mtDNA and nuclear data indicate that the extant Australasian agamid radiation is between 14 and 41 Myr old (i.e., Oligo-Miocene). This age is consistent with the time frames for dispersal between Asia and Australasia suggested by plate tectonics (Cogger and Heatwole 1981). The gradual northward drift of the Australasian plate does not provide absolute constraints on the timing of dispersals (Hutchinson and Donnellan 1993; contra Hall [1996]). All that can be assumed is that during the Tertiary, dispersal between Asia and Australasia became progressively more likely with time as the plates approached one another. The latest studies suggest that short-range, island-hopping dispersal between Asia and Australasia could have occurred at least as long ago as 30 Myr (Metcalfe 2001), and longer-range dispersal was possible much earlier. The arboreal and/or amphibious habits of many acrodonts (including P. cocincinus and many basal Australasian agamids) would make them reasonable candidates for early dispersal through rafting. Schulte, Melville, and Larson (2003) suggested that agamids were not adept over-water dispersers; however, the three clades of iguanians (chamaeleons, agamids, and iguanids) all appear to contain forms that readily cross water barriers. Chameleons have been inferred to have crossed repeatedly between Madagascar and Africa (Raxworthy, Forstner, and Nussbaum 2002), and agamids and iguanids are found on numerous islands that have always been separated from the mainland by large water gaps (e.g., Moody 1980).
The dates estimated here are also compatible with the presence of agamids in Africa, Asia, and Australasia and their absence from the Americas. This is consistent with the relatively recent origin of agamids in Eurasia (where they potentially have their earliest fossil record at 190 MYA [Evans, Prasad, and Manhas 2002]), followed by dispersal into the Gondwanan landmasses of Africa, India, and Australasia when they drifted northwards. Since the breakup of Pangea, South America has never been connected to Eurasia, and North America connected only at high latitudes. The dates inferred here for basal divergences in agamids are also broadly consistent with the fossil record, placing the agamid radiation within the past 200 Myr, matching paleontological estimates of the timing of the lizard radiation (Estes 1983; Evans 1993).
Within agamids, the mtDNA data show large divergences, evidence of saturation effects, and considerable sensitivity to details of substitution model and rate-smoothing methods. In contrast, the large portion of c-mos appears less affected, performs well at retrieving basal nodes and giving stable estimates for deeper branch lengths in iguanians and should be a useful locus for other similar studies.
Acknowledgements
We thank Mark Hutchinson, Steve Donnellan, Remko Leijs, and two anonymous referees for discussion and comments; Kathleen M. Saint, Steve Donnellan and Mark Hutchinson for generously allowing access to their primer information and c-mos sequences; and the Australian Research Council for funding.
References
Audley-Charles, M. G. 1987. Dispersal of Gondwanaland: relevance to evolution of the angiosperms. Pp. 5–25 in T. C. Whitmore, ed. Biogeographical evolution of the Malay Archipelago. Oxford Monographs on Biogeography No. 4. Clarendon Press, Oxford.
Benton, M. J., and F. J. Ayala. 2003. Dating the tree of life. Science 300:1698–1700.
Baldwin, B. G., and M. J. Sanderson. 1998. Age and rate of diversification of the Hawaiian silversword alliance (Compositae). Proc. Nat. Acad. Sci. USA 95:9402–9406.
Barraclough, T. G., and A. P. Vogler. 2002. Recent diversification rates in North American tiger beetles estimated from a dated mtDNA phylogenetic tree. Mol. Biol. Evol. 19:1706–1716.
Cogger, H. G., and H. Heatwole. 1981. The Australian reptiles: origins, biogeography, distribution patterns and island evolution. Pp. 1333–1373 in A. Keast, ed. Ecological biogeography of Australia. W. Junk, The Hague.
Estes, R. 1983. The fossil record and early distribution of lizards. Pp. 365–398 in A. G. J. Rhodin and K. Miyata, eds. Advances in herpetology and evolutionary biology: essays in honour of Ernest E. Williams. Museum of Comparative Zoology, Cambridge, Mass.
Evans, S. E. 1993. Jurassic lizard assemblages. Revue de Paléobiologie (Suisse), Volume spéciale 7:55–65.
Evans, S. E., G. V. R. Prasad, and B. K. Manhas. 2002. Fossil lizards from the Jurassic Kota Formation of India. J. Vert. Paleont. 22:299–312.
Glazko, G. V., and M. Nei. 2003. Estimation of divergence time for major lineages of primate species. Mol. Biol. Evol. 20:424–434.
Graur, D., and W. Martin. 2004. Reading the entrails of chickens: molecular timescales of evolution and the illusion of precision. Trends Genet. 20:80–86.
Greer, A. E. 1989. The biology and evolution of Australian lizards. Surrey Beatty, Chipping Norton.
Hall, R. 1996. Reconstructing Cenozoic SE Asia. Pp. 153–184 in R. Hall and D. Blundell, eds. Tectonic evolution of Southeast Asia. Geological Society, London.
———. 1998. The plate tectonics of Cenozoic SE Asia and the distribution of land and sea. Pp. 99–131 in R. Hall and J. D. Holloway, eds. Biogeography and geological evolution of SE Asia. Backhuys, Leiden.
Harris, D. J. 2003 Codon bias variation in C-mos between squamate families might distort phylogenetic inferences. Mol. Phylogenet. Evol. 27:540–544.
Hedges, S. B., and S. Kumar. 2004. Precision of molecular time estimates. Trends Genet. 20:242–247.
Hutchinson, M. N., and S. C. Donnellan. 1993 Biogeography and phylogeny of the Squamata. Pp. 210–220 in C. J. Glasby, G. J. B. Ross, and P. L. Beesley, eds. Fauna of Australia, Vol 2A: Amphibia and Reptilia. AGPS Press, Canberra.
Felsenstein, J. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368–376.
———. 2004. Inferring phylogenies. Sinauer, Sunderland, Mass.
Lee, M.S.Y. 1999. Molecular clock calibrations and metazoan divergence dates. J. Mol. Evol. 49:385–391.
Macey, J. R., J. A. Schulte II, N. B. Ananjeva, A. Larson, N. Rastegar-Pouyani, M. Shammakov, and T. Papenfuss. 1998. Phylogenetic relationships among agamid lizards of the Laudakia caucasia species group; testing hypotheses of biogeographic fragmentation and an area cladogram for the Iranian plateau. Mol. Phylogenet. Evol. 10:118–131.
Macey, J. R., J. A. Schulte II, A. Larson, N. B. Ananjeva, Y. Wang, N. R. Pethiyagoda, N. Rastegar-Pouyani, and T. J. Papenfuss. 2000. Evaluating trans-Tethys migration: an example using Acrodont lizard phylogenetics. Syst. Biol. 49:233–256.
Martin, A. P., E. K. Costello, A. F. Meyer, D. R. Nemergut, and S. K. Schmidt. 2004. The rate and pattern of cladogenesis in microbes. Evolution 58:946–955.
Metcalfe, I. 1996. Precretaceous evolution of SE Asian terranes. Pp. 97–122 in R. Hall and D. Blundell, eds. Tectonic evolution of Southeast Asia. Geological Society, London.
———. 2001. Palaeozoic and Mesozoic tectonic evolution and biogeography of SE Asia-Australasia. Pp. 15–56 in I. Metcalfe, ed. Faunal and floral migrations and evolution in SE Asia-Australasia. Balkema Publishers, Exton, Pa.
Moody, S. 1980. Phylogenetic and historical biogeographical relationships of the genera in the Agamidae (Reptilia, Lacertilia). Doctoral dissertation, University of Michigan, Ann Arbor.
Posada, D., and K. Crandall, K. 1998. Modeltest: testing the model of DNA substitution. Bioinformatics 14:817–818.
Rambaut, A., and L. Bromham. 1998. Estimating divergence dates from molecular sequences. Mol. Biol. Evol. 15:442–448.
Rambaut, A., and M. Charleston. 2001. TreeEdit: phylogenetic tree editor. Version 1.0. Department of Zoology, University of Oxford. URL http://evolve.zoo.ox.ac.uk/software/TreeEdit/main.html
Raxworthy, C. J., M. R. J. Forstner, and R. Nussbaum. 2002. Chameleon radiation by oceanic dispersal. Nature 415:784–787.
Sanderson, M. J. 1997. A nonparametric approach to estimating divergence times in the absence of rate constancy. Mol. Biol. Evol. 14:1218–1231.
———. 2002. Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Mol. Biol. Evol. 19:101–109.
Schulte, J. A., J. Melville, and A. Larson. 2003. Molecular phylogenetic evidence for ancient divergence of lizard taxa on either side of Wallace's Line. Proc. R. Soc. Lond. B Biol. Sci. 270:597–603.
Schwaner, T. D., P. R. Baverstock, H. C. Dessauer, and G. A. Mengden. 1985 Immunological evidence for the phylogenetic relationships of Australian elapid snakes. Pp. 177–184 in G. C. Grigg, R. Shine, and H. Ehmann, eds. The biology of Australasian frogs and reptiles. Royal Zoological Society of NSW, Sydney.
Swofford, D. L. 2002. PAUP*: phylogenetic analysis using parsimony (*and other methods). Version 4. Sinauer Associates, Sunderland, Mass.
Weisrock, D. W., R. J. Macey, I. H. Ugurtas, A. Larson, and T. J. Papenfuss. 2001. Molecular phylogenetics and historical biogeography among Salamadrids of the ‘True’ Salamader clade: rapid branching of numerous highly divergent lineages in Mertensiella luschani associated with the rise of Anatolia. Mol. Phylogenet. Evol. 18:434–448.
Woodburne, M. O., and J. A. Case. 1996. Dispersal, vicariance, and the late Cretaceous to early Tertiary land mammal biogeography from South America to Australia. J. Mamm. Evol. 3:121–161.(Andrew F. Hugall* and Mic)