• E-mail
• Print
• Comment
• Font Size
• Digg
• del.icio.us
• Discuss article

Molecular Phylogenetics of Squamata: The Position of Snakes, Amphisbaenians, and Dibamids, and the Root of the Squamate Tree

Posted on: Thursday, 2 December 2004, 03:00 CST

Abstract.-

Squamate reptiles (snakes, lizards, and amphisbaenians) serve as model systems for evolutionary studies of a variety of morphological and behavioral traits, and phylogeny is crucial to many generalizations derived from such studies. Specifically, the traditional dichotomy between Iguania (anoles, iguanas, chameleons, etc.) and Scleroglossa (skinks, geckos, snakes, etc.) has been correlated with major evolutionary shifts within Squamata. We present a molecular phylogenetic study of 69 squamate species using approximately 4600 (2876 parsimony-informative) base pairs (bp) of DNA sequence data from the nuclear genes RAG-1 (~2750 bp) and c-mos (~360 bp) and the mitochondrial ND2 region (~1500 bp), sampling all major clades and most major subclades. Under our hypothesis, species previously placed in Iguania, Anguimorpha, and almost all recognized squamate families form strongly supported monophyletic groups. However, species previously placed in Scleroglossa, Varanoidea, and several other higher taxa do not form monophyletic groups. Iguania, the traditional sister group of Scleroglossa, is actually highly nested within Scleroglossa. This unconventional rooting does not seem to be due to long-branch attraction, base composition biases among taxa, or convergence caused by similar selective forces acting on nonsister taxa. Studies of functional tongue morphology and feeding mode have contrasted the similar states found in Sphenodon (the nearest outgroup to squamates) and Iguania with those of Scleroglossa, but our findings sviggest that similar states in Sphenodon and Iguania result from homoplasy. Snakes, amphisbaenians, and dibamid lizards, limbless forms whose phylogenetic positions historically have been impossible to place with confidence, are not grouped together and appear to have evolved this condition independently. Amphisbaenians are the sister group of lacertids, and dibamid lizards diverged early in squamate evolutionary history. Snakes are grouped with iguanians, lacertiforms, and anguimorphs, but are not nested within anguimorphs. [Amphisbaenia; Dibamidae; DNA; Iguania; lizards; long-branch attraction; mitochondrial; nuclear; phylogeny; Scleroglossa; Serpentes, Squamata.]

Evolutionary biologists often seek generalities about evolutionary processes from detailed studies of particular model systems, and squamate reptiles have provided a large number of such systems (e.g., Huey et al., 1983; Vitt and Pianka, 1994). An accurate squamate phylogeny is crucial to studies of morphological, behavioral, and life-history variation because phylogeny is a key part of comparative methodology (Miles and Dunham, 1993). For example, herpetological studies of foraging mode and prey chemical discrimination (Cooper, 1995; Perry, 1999), demographic tactics (Clobert et al, 1998), and home-range variation (Perry and Garland, 2002) have all explicitly incorporated phylogeny into their testing framework to insure appropriate, independent comparisons. The evolution of squamate tongue morphology and chemoreception abilities is cited as a prime example of the importance of history in present- day distribution patterns and ecology (Schwenk and Wagner, 2001; Vitt et al., 2003). A subtle shift in prey-prehension technique is thought to have allowed major changes in tongue morphology and chemosensory abilities (Cooper, 1995; Schwenk, 1993; Schwenk and Wagner, 2001). This shift coincided with the Early Jurassic split between the two major squamate clades, Iguania and Scleroglossa, and allowed the scleroglossans to exploit a variety of new habitats and foraging modes unavailable to iguanians, such that scleroglossans now predominate on a global scale, over 200 million years (my) later (Schwenk and Wagner, 2001; Vitt ct al., 2003). All of these inferences are heavily dependent on a correct rooting of the squamate tree, which itself is dependent on comparison of character states in outgroup taxa.

The nearest extant outgroup to squamates consists only of the superficially lizard-like tuataras from New Zealand, the only remaining members of a once much more widespread and diverse group (Evans et al., 2001; Reynoso, 2000). Squamata itself is a diverse assemblage including all reptiles commonly called lizards plus three limbless groups: snakes, amphisbaenians, and dibamid lizards. Table 1 gives a summary of current squamate classification and includes the higher taxon names used in this paper. Recent morphological studies (Estes et al, 1988; Lee, 1998; Lee and Caldwell, 2000; Reynoso, 1998) agree on some groupings of families into higher taxa (e.g., Anguimorpha, Iguania, Scleroglossa), and on the basal dichotomy between Iguania and Scleroglossa, but the phylogenetic relationships among many higher squamate taxa remain uncertain. Historically, the limbless clades have been particularly difficult to place based on morphology because the limbless condition eliminates many characters, although utilizing fossil taxa can sometimes help (e.g., Lee and Caldwell, 1998; Zaher and Rieppel, 1999). Furthermore, because limblessness is often associated with a fossorial lifestyle, cranial morphology in these animals is also often radically altered from that of nonburrowing squamates (Lee, 1998).

Snakes are by far the most ecologically diverse and familiar of these limbless groups, with over 2900 species occupying a variety of terrestrial, arboreal, fossorial, and aquatic habitats on all major land masses except Antarctica (Pough et al., 2004). Amphisbaenians are a much more homogeneous, completely fossorial group with major radiations in South America and northern Africa (extending into southern Europe), and two small clades confined to North America. One of these, the Bipedidae, is unique among amphisbaenians in its retention of forelimbs. Dibamid lizards are a small, poorly known, completely fossorial group with a curiously disjunct distribution (Southeast Asia/Suncla Shelf and northeastern Mexico), which suggests that this group was previously more widespread. Although snakes (Rieppel, 1983; Rieppel, 1985), amphisbaenians (Gans, 1978), and dibamids (Greer, 1985) each exhibit characters that might place them phylogenetically outside all other squamates, most authors now agree that all three are probably nested within lizards. Estes et al. (1988) designated these three groups as "Scleroglossa incertae sedis," and many morphological studies (e.g., Greer, 1985; Hallermann, 1998; Lee, 1998; Lee and Caldwell, 2000; Reynoso, 1998; Rieppel, 1984; Wu et al., 1996) have specifically addressed placement of these taxa. A common feature of most recent morphological studies (e.g., Lee, 1998; Lee and Caldwell, 2000; Reynoso, 1998; Wu et al., 1996) is the grouping of at least two of these limbless groups in a clade.

TABLE 1. Phylogenetic taxonomy of sequamntes based on morphology.a

Pew molecular studies have addressed higher-level relationships within Squama ta on a broad scale. Most studies concerned with suprafamilial relationships have had limited outgroup taxon sampling (e.g., Ast, 2001; Donnellan et al., 1999; Macey et al., 1997b, 1999, 2000; Odierna et al., 2002; Saint et al., 1998; Whiting et al., 2003) and were thus not designed to produce a comprehensive higher- level squamate phylogeny, although several specific points have been clarified through these molecular studies. Donellan et al. (1999) support Kluge's (1987) conclusion that the Australian Pygopodidae (another limbless group) is closely related to the Australian diplodactyline geckos. Ast (2001) reports monophyly of Varanoidea (Varanidae + Helodcrmatidac), but only relative to anguids and anniellids, because Xenosaurus and SMmscmrHS were not included as outgroups. Whiting et al. (2003) find strong support for the New World Xantusiidac as the sister taxon of the African Cordylidae.

Harris (2003) and Harris et al. (2001; 1999) use sequence data from the nuclear proto-oncogene c-mos to investigate higher squamate relationships (ultimately 162 sequences representing all major squamate clades except Dibamidae were analyzed). Many of the higher- level relationships recovered in these studies conflict with those of morphological studies. However, the gene fragment used for these studies is only approximately 360 bp, and most basal relationships are weakly supported. Vidai and Hedges (2004) use this same c-mos fragment along with approximately 500 bp of the protein-coding nuclear gene RAG-1 to investigate relationships among major snake taxa as well as the position of snakes within Squamata. This study, which includes representatives of all recognized squamate families, finds strong bootstrap support for snakes in a phylogenetic position outside of Anguimorpha, contradicting several morphological studies (e.g., Lee, 1998; Lee and CaIdwell, 2000; McDowcll and Bogert, 1954). However, most other basal squamate relationships are not well supported.

Here we present a phylogenetic study of Squamata using three independent molecular data sets. Along with mitochondria! data from the ND2 region, we have utilized the same c-mos fragment as the previous studies, and we have collected data from almost t\he entire length of RAC-1.

MATERIALS AND METHODS

Taxon Sampling

Rhynchocephalia is traditionally considered the closest outgroup to Squamata, although some molecular evidence (Hedges and Poling, 1999) suggests the arrangement (Rhynchocephalia, (Testudines, Archosauria)). We have therefore included representatives from all three of these taxa as outgroups. Within the ingroup, all recognized major squamate clades (i.e., lizard families, amphisbaenians, dibamids, and snakes) are represented, as well as many major subclades. In diverse families (or other equivalent taxa), we attempted to sample species from both sides of the most basal divergence, as inferred from morphological studies and/or previous molecular work, as follows: Agamidae (Frost and Etheridge, 1989; Macey et al., 2000), Chamaeleonidae (Klaver and Bohme, 1986; Townsend and Larson, 2002, unpublished data), Gekkonidae (Donnellan et al., 1999; Kluge, 1987), Amphisbaenia (Kearney, 2003), Scincidae (Greer, 1970; Whiting et al., 2003), Cordylidae (Lang, 1991; Odierna et al., 2002), Serpentes (Heise et al, 1995; Rieppel, 1988), and Anguidae (Gauthier, 1982; Macey et al, 1999). For families with uncertain intrafamilial relationships, e.g., Iguanidae (Frost and Etheridge, 1989), we tried to sample all major subclades to assure that the deepest divergence was spanned. A total of 69 ingroup species were sampled for RAG-1 and the mitochondrial fragment, and 44 ingroup species (including all major clades) were sampled for c- mos. Almost all RAG-1 and mitochondrial sequences were collected from the same individuals; however, many of the c-mos sequences are from previously published studies, and exact species matches were not always possible. For combined analyses in which species representing particular higher taxa were not the same across all data partitions, we have labeled the resulting trees with the name of the higher taxon. For example, draconine agamids are represented by Calotes calotes and C. versicolor in the RAG-1 and c-mos data sets, respectively, and their concatenated sequences are labeled "Draconinae" in the combined analyses. see Appendix for museum and GenBank accession numbers for all specimens.

Laboratory Protocols

Genomic DNA was extracted from muscle, liver, or skin tissue (stored either frozen or in 70% to 95% ethanol) using DNEasy Tissue Extraction Kits (Qiagen, Inc.) and stored in AE buffer. Mitochondrial polymerase chain reaction (PCR) products were amplified from genomic DNA using an initial denaturation at 95C for 2 min, then a denaturation at 95C for 35 s, annealing at 5OC for 35 s, and extension at 7OC for 150 s with 4 s added to each successive extension cycle for 30 cycles. Nuclear genomic DNA was originally amplified using the touchdown protocol of Groth and Barrowclough (1999). All PCR products were purified on 1.3% low-melt agarose gels and reamplified under the same conditions used to amplify the mitochondrial genomic DNA, except that the annealing temperature was reduced to 45C. Promega Taq polymerase (Promega, Inc., Madison, Wisconsin) was used for all amplifications. Reamplified products were purified on 0.8% high-melt agarose gels, and template extracted using Viogene Gel Extraction Kits (Viogene, Inc., Taipei, Taiwan) and sequenced in both directions using ABI Prism Dye Terminator Cycle Sequencing Ready Reaction Kits with AmpliTaq DNA Polymerase (Perkin Elmer, Norwalk, Connecticut) following the manufacturer's instructions. Sequencing products were analyzed with ABI373 or 377 Automated Sequencers (Applied Biosystems, Foster City, California) or an MJ Research BaseStation (MJ Research, San Francisco, California). Multiple overlapping PCR fragments were sequenced for the mitochondrial and RAG-1 partitions, which helped guard against PCR contamination problems. When certain taxa fell in unexpected places in some trees (e.g., Heloderma), multiple individuals or closely related taxa were at least partially sequenced when possible to confirm our findings.

Primers used in this study can be found on the Systematic Biology website (http://systematicbiology.org/). Mitochondrial sequences cover an approximately 1500 bp region corresponding to positions 4419 to 5933 on the human mitochondrial genome (Anderson et al., 1981), including portions of AfDJ (subunit 1 of NADH dehydrogenase) and COl (subunit I of cytochrome c oxidase), the 8 tRNA genes for glutamine, isoleucine, methionine, tryptophan, alanine, asparagine, cysteine, and tyrosine, the entire ND2 (subunit 2 of NADH dehydrogenase) gene, and the stem-and-loop structure representing the origin for light-strand replication (O]J. RAG-1 sequences cover an approximately 2800-bp coding region corresponding to positions 84 to 3126 on the published chicken RAG-1 gene (Carlson et al., 1991). C-mos sequences cover an approximately 374-bp coding region corresponding to positions 513 to 888 on the human c-mos gene (Watson et al., 1982). Forty-four mitochondrial and 35 c-mos sequences were obtained from GenBank (see Appendix), and 28 mitochondrial and 12 c-mos sequences were newly generated for this study, as were 71 of the 73 RAG-1 sequences. All mitochondrial GenBank sequences were previously generated in our laboratory, and were checked against the corresponding sequences in their original alignments to verify their identity. When possible, c-mos sequences from taxa closely related to those used in our study were also downloaded from GenBank for comparison to help detect mislabeled sequences (no mistakes were found).

Alignments and Phylogcnetic Analyses

Sequences were edited and assembled using SeqMan II (DNASTAR, Inc., Madison, Wisconsin). Alignments of protein-coding regions were performed on amino acid translations using Clustal X (Thompson et al., 1997) at a variety of gap-opening and gap-extension penalties. For pairwise augments, gap-opening penalties were set to 10, 20, and 35 with respective gap-extension penalties of 0.1, 0.45, and 0.75. Corresponding multiple-alignment penalties were 10, 15, and 20 (gap- opening) and 0.1, 0.2, and 0.3 (gap-extension). Regions for which alignments differed between the three suites of settings were excluded from all analyses. Genes coding for tRNAs were aligned manually following the structural models of Kumazawa and Nishida (1993). Length-variable loops that could not be confidently aligned were excluded from all analyses. All gaps were treated as missing data. There are several opinions on combining data from different sources for phylogenetic analysis. One is that all available data should be included in any analysis (Kluge, 1989), and a second is that partitions should always be analyzed separately, with congruence between partitions inferred as strong support for particular relationships (Miyamoto and Fitch, 1995). Finally, a third alternative is to test for congruence between data partitions, then combine them if the test is passed (e.g., Bull et al., 1993; Farris et al., 1994; Huelsenbeck and Bull, 1996; Rodrigo et al., 1993). However, it is difficult to take a strictly formulaic approach to this problem. For example, as discussed by Wiens (1998), significant global tests of data incongruence do not indicate whether the incongruence is spread throughout the data or is concentrated in specific parts, and therefore not performing combined analyses based solely on the results of these tests seems inappropriate. Furthermore, recently developed methods of mixed- model analysis (e.g., Nylander et al., 2004) may make this issue less relevant (barring horixontal transfer and other similar evolutionary events). Our strategy was to perform both separate and combined analyses, and, in the case of conflicting topologies, to examine characteristics (rates of evolution, relative branch lengths, etc.) of each data partition to find potential explanations for the conflict (see Results for discussion of specific cases).

To explore the possibility of heterogeneous selective pressures on the protein-coding nuclear genes (sec Results), we used DnaSP (Rozas and Rozas, 1999) to calculate ratios of synonymous substitutions per synonymous site (K^sub s^) to nonsynonymous substitutions per nonsynonymous site (K^sub a^) for all possible pairwise taxon comparisons. Average K^sub a^/K^sub s^ ratios were then calculated within each major clade as well as among clades. Average K^sub a^/K^sub s^ ratios significantly greater than one (as determined by t-tests) indicate directional selection in at least some of the species/clades compared, whereas ratios significantly less than one indicate stabilizing or purifying selection (see Messier and Stewart, 1997 for a more detailed discussion).

The model of evolution and all maximum-likelihood (ML) parameters were estimated for each data set individually using hierarchical likelihood-ratio tests as implumunLcd in Modeltest (fosoda and Crandall, 1998). Maximum-likelihood analyses were conducted using the heuristic search option of PAUP* (Swofford, 1998) and a neighbor- joining tree as a starting tree for TBR branch swapping. Computational limitations precluded the use of nonparametric bootstrapping under the likelihood criterion.

Baycsian analyses were performed using MrBayes 3.0b3 (Ronquist and Huelsenbeck, 2003) under the same model used for the corresponding likelihood analyses. One major concern with combining separate data sets is that evolutionary models may differ substantially between them (see Huelsenbeck et al., 1996). Version 3 of MrBayes allows parameter values to be estimated separately (under potentially different evolutionary models) for different data partitions, and this has the potential to alleviate this problem (Nylander et al., 2004). However, as more parameters are estimated, the potential for loss of statistical power increases. The magnitude of this problem is not fully understood, but it seems likely that excessive partitioning of the data could create its own probl\ems. We have therefore used multiple partitioning schemes in our combined analyses and compared their effects on topology and branch support (see Results). For all Dayesian analyses, four incrementally heated Markov chains were started from random trees and run for 4,000,000 generations each. The effect of heating the chains is to increase the probability of acceptance of new parameter-value propositions; this flattens the landscape somewhat, allowing the chains to cross valleys and to explore trccspacc more effectively. Chains were sampled every 1000 generations to ensure that the samples were independent. Through inspection of the likelihood scores and model parameters in the output file, we determined the number of generations required for stabilization and discarded all trees obtained prior to stabilisation as burnin. Two independent analyses were conducted for each data set, and their resulting topologies, posterior clade probabilities, and log-likelihood (InL) values at stationarity were compared to prevent drawing the posterior distributions from local optima. Trees from the posterior distribution were imported into PAUP* (Swofford, 1998) and, after excluding the burn-in, a majority-rule consensus tree was constructed showing relative occurrences (i.e., the posterior probabilities) of all nodes in the tree.

Maximum-parsimony (MP) analyses were performed using PAUP* (Swofford, 1998) under the heuristic search option with 100 random- addition replicates. Nonparametric bootstrap resampling was applied to assess heuristic support for individual nodes (Felsenstein, 1985b) using 1000 bootstrap pseudoreplicates with 25 random additions per pseudoreplicate. Branch-support (decay) indices (Bremer, 1994) were calculated as heuristic support measures for all resolved internal branches of the tree using the "Decay Index PAUP File" feature of MacClade (Maddison and Maddison, 2000). DeBry (2001 ) showed that the variance among significant decay indices within a single tree can be large, and that rigorous interpretation of decay values must take branch lengths into account. However, in the absence of explicit statistical testing of each node, we feel that in many cases the decay index is still useful as a rough guide to relative levels of support (especially once bootstrap values reach their maximum of 100). As an indicator of relative homoplasy among data sets, retention indices (Farris, 1989) were also calculated.

Testing Alternative Topologies

Statistical support for individual nodes was assessed using two separate nonparametric tests, the parsimony-based Wilcoxon signcd- ranks (Templeton) test (Felsenstein, 1985a; Templeton, 1983) and the likelihood-based SH test (Shimodaira and Hasegawa, 1999). Felsenstein (1985a) showed that one-tailed probabilities for the Templeton test are close to the exact probabilities and that use of two-tailed probabilities is always conservative. Consistent with these findings, the two-tailed version of this test is generally conservative relative to alternative tests that ask whether character data statistically discriminate alternative phylogenetic topologies (e.g., Lee, 2000; Townsend and Larson, 2002). Bonferroni corrections for multiple tests (Rice, 1989) are very conservative, and were not applied to this already conservative test. Goldman et al. (2000) commented that the Templeton test is appropriate only when all trees being tested are specified a priori, because by using the best (MP) tree derived from the data at hand, the test is potentially biased to be less conservative. The magnitude of this potential bias is unknown.

Goldman et al. (2000) suggested using instead the SH test (Shimodaira and Hasegawa, 1999), which uses a resampling method to overcome this potential bias and also makes corrections for multiple comparisons. Theoretically, this test requires that all possible topologies be compared simultaneously, an obvious impossibility with data sets of more than a few taxa. Buckley et al. (2001) suggested restricting the set of possible topologies to only those reasonably likely to be the true topology, but even this is impractical with most data sets. In testing particular nodes, we conducted both of these tests as follows. First, constraint trees containing only a single resolved node were constructed using MacClade (Maddison and Maddison, 2000). Next, for the Templeton test, the shortest trees not containing this node were found using PAUP* (Swofford, 1998), and these trees were then compared to the shortest unconstrained tree using the "Tree Scores" option of PAUP* (Swofford, 1998). An analogous procedure (using likelihoods instead of tree lengths) was followed for the SH test. This use of the SH test reduces to an appropriately "centered" KH test (Goldman et al., 2000). Because neither test was performed under technically perfect conditions, borderline-significant results should be interpreted with caution.

Our RAG-1 tree roots at one of two relatively long branches (MP roots it at Dibamus, ML at Gekkonidae; see Figure 1 and TreeBASE website). To test the hypothesis that long-branch attraction (LBA; Felsenstein, 1978) might cause an aberrant rooting, we followed Wiens and Hollingsworth's (2000) implementation of the parametric bootstrapping method of Huelsenbeck (1997). We first rerooted our RAG-1 ML topology at Iguania (the morphological root; see Results) and reoptimized all branch lengths and other model parameters on this tree. Next, we used Seq-Gen (Rambaut and Crassly, 1997) to simulate 100 data sets on this topology, with sequence length equal to that of the RAG-1 data set. MP and ML analyses were performed on each of these data sets, and a tally was kept for each optimality criterion of the number of correct and incorrect roofings (relative to the simulated topology), as well as the number of times the tree rooted specifically at either geckos or Dibamus. As a general rule with this type of test, if parsimony tends to root the tree incorrectly but likelihood does not, this suggests that LBA is a potential problem. In this specific case, if either analysis tends to recover a tree rooted at geckos and/or Dibamus, the original rooting from the real RAG-1 data would be highly suspect. Because unconstrained ML analyses were not computationally feasible, intrafamilial relationships were constrained to match those from the original RAG-1 ML analysis, but interfamilial relationships were free to vary.

Results from the mitochondrial analysis suggest that LBA might occur between two specific clades at the end of long internal branches. We therefore performed a similar study with the mitochondrial data and topology, this time separating the two suspicious long branches in the model tree used for the simulations. In this case, if parsimony tends incorrectly to join the long branches while likelihood does not, this result suggests that LBA is a problem.

The most consistent well-supported difference between our topology and the topologies found in all recent morphological studies concerns the placement of the squamate root. If the morphological rooting is affected by misleading convergence between two or more taxa, the characters supporting this rooting might be found to be concentrated in one anatomical area. This would not necessarily be true for all scenarios but if, for example, there were convergence in feeding morphology between sphenodontids and iguanians (and iguanians did not actually branch off early in squamate history), we might expect that skull and jaw characters might be overrepresented in the list of characters supporting a rooting at Iguania. Likewise, if convergence is a problem with the molecular rooting, we might expect its supporting characters to be concentrated in one particular genie functional domain. To identify the source of conflict on this point, we performed separate tests using Lee and Caldwell's (2000) morphological data (including fossils) and our RAG-1 data. First, using Lee and Caldwell's (2000) taxa, we constructed a tree congruent with our RAG-1 ML topology. We then made a second topology by rerooting this tree at Iguania (in accordance with morphological hypotheses). After excluding 27 characters identified by Lee (1998) as potentially correlated to a fossorial existence (this exclusion does not reduce support for the morphological rooting), we mapped Lee and Caldwell's (2000) morphological characters onto each of these trees and identified those that required more steps on the tree with the molecular rooting. Using the anatomical divisions given in Lee and Caldwell (2000), we then performed chi-square tests to determine if one or more anatomical regions were overrepresented in the list of characters opposing the molecular rooting (see Harshman et al., 2003, for a similar use of this test). Performing an exactly analogous test on the molecular data is difficult, because some functional domains of the gene are known to be more variable than others (Willett et al, 1997), and the variable regions might be expected to contain proportionately more informative characters than the conserved regions. As a proxy, we divided the RAG-2 gene into two regions, the more highly variable 5' one-third of the gene that codes for protein-binding sites, and the more highly conserved 3' two-thirds of the gene responsible for target-site recognition and DNA binding (reviewed in Willett et al., 1997). These regions were analyzed separately to check for congruence among their respective topologies.

RESULTS

Phylogenetic Results from Nuclear Genes

All complete, aligned data files (with excluded positions marked as such), along with trees from all individual data sets, are available on the TreeBASE website (http://www.treebase.org/treebase/ ). The RAG-1 MP, ML, and Bayesian topologies are all very similar, and all nodes receiving high heuristic support from the parsimony analysis (bootstrap >90%) also have posterior probabil\ities >95% in the Bayesian analysis (Fig. 1). The c-mos ML and Bayesian topologies are very similar, and both analyses recover all moderately to highly supported nodes (bootstrap >80%) from the c-mos parsimony analysis (Fig. 2). Although the c-mos data set contains fewer species, all major clades from the RAG-1 analysis are still represented. The topology of the c-mos MP strict consensus tree is largely compatible with the RAG1 topology, although many deeper relationships are not resolved (Fig. 2).

The model parameters estimated for the RAG-1 and c-mos genes are similar (although Modeltest chose the simpler HKY model for c-mos, probably due to the relatively short length of this data set), and the two genes have very similar levels of divergence among squama tes (Table 2, Figs. 1 and 2). Sequence data from these two genes were combined and analyzed as a single nuclear data set (Fig. 3). For clarity, results of these combined analyses are detailed here, with references to individual analyses (Figs. 1 and 2; see TreeBASE website for individual MP and ML/Bayesian topologies) as necessary. Parsimony, Bayesian, and likelihood topologies from the combined RAG- 1 and c-mos data are largely congruent with each other as well as with corresponding trees from each data set analyzed singly. RAG-1/ c-mos parsimony and Bayesian support values are at least as high as those from the RAG-1 data alone, and often substantially higher (Figs. I and 3). Combined Bayesian analyses were run both with parameter values estimated separately for the RAG-1 (GTR+l+G) and c- mos (HKY+I+G) partitions (Fig. 3) and also as a single data partition under a GTR+I+G model (results not shown). Topologies were identical between these analyses, and posterior probabilities were very similar as well.

TABLE 2. Properties of character variation for all protein coding genes (analyzed by first, second, and third codon positions), plus tRNA genes of the mitochondrial genome.

In all analyses, when more than one subclade is represented, monophyly of almost all recognized families is recovered with strong support (the one exception is a paraphyletic Agamidae found in ML and Bayesian analyses) (Fig. 3). Note that because Estes et al. (1988) defined their taxa so that taxon names would remain stable, it is technically impossible that Scincidae, for example, could be nonmonophyletic (i.e., only its taxon composition can change). For brevity, however, we will make reference to monophyly, paraphyly, etc. of these taxa throughout this paper, with the understanding that we are actually referring to the groups placed in these taxa by Estes et al. (1988) at the time of their definition. Traditional suprafamilial groups recovered with strong support in the RAG-1/c- mos analyses include Acrodonta, Iguania, Anguimorpha (also characterized by a one-codon insertion at positions 128 to 130 in the aligned RAG-1 data set), and Teioidea. Interestingly, several nontraditional relationships are also recovered with strong support.

Tests of Phylogenetic Rooting

Most striking is the absence of a monophyletic ScIeroglossa as the sister taxon of Iguania (Fig. 3). Instead, the deepest divergence is between Dibamus and all other squamates, and Iguania occupies a highly nested position in the tree. Both a paraphyletic Scleroglossa and highly nested Iguania are contradicted by substantial morphological evidence (e.g., Lee and Caldwell, 2000). Furthermore, this specific conflict between molecules and morphology can be resolved by simply rerooting our nuclear topology at Iguania. We therefore considered four scenarios that might have led to an incorrect rooting caused by misleading sequence convergence between nonsister taxa in the nuclear analyses: heterogeneous base composition among taxa/ heterogeneous selection pressures, phylogcnetic randomization of outgroup sequences with respect to ingroup sequences, and long-branch attraction.

FIGURE 1. RAG-J data, ML phylogram (CTR-U+G model; -InL = 50519.59; A = 0.3007, C = 0.2242, G = 0.2254, T = 0.2497; AC = 1.3332, AG = 4.7011, AT = 0.9186, CG = 0.8644, CT = 5.7274, GT = 1.0; 1 = 0.3352; G = 1.7108). Asterisks indicate branches that receive a posterior probability of 95% or greater in the Bayesian analysis. MP bootstrap proportions >70% (above branches) and decay indices (below branches) are provided for all nodes congruent between analyses based on the two optimality criteria. Numbers to the right denote major taxonomic units as follows: 1. Chamaeleonidae; 2. Agamidae; 3. Iguanidae; 4. Anguidae; 5. Helodermatidae; 6. Xenosauridae; 7. Varanidae; 8. Shinisauridae; 9. Serpentes; 10. Lacertidae; 11. Amphisbaenia; 12. Teiidae; 13. Gymnophthalmidae; 14. Scincidae; 15. Xantusiiclae; 16. Cordylidae; 17. Dibamidae; 18. Gekkonidae. Outgroup branches with hatch marks have been shortened.

FIGURE 2. C-mos data, ML phylogram (HKY+I+C model, -InL = 5431.52; A = 0.2712, C = 0.2527, C = 0.2231, T = 0.2530; Ti/Tv = 2.4357; I = 0.3084; G - 3.5785). Asterisks indicate branches that receive a posterior probability of 95% or greater in the Baycsian analysis. MP bootstrap proportions >70% (above branches) and decay indices (below branches) are provided for all nodes congruent between analyses based on the two optimality criteria. Annotations and numbering as in Figure 1.

FIGURE 3. Combined RAG-1 and c-mos data, ML phylogram (GTR+I+G model; -InL = 44562.45; A = 0.3022, C = 0.2212, G = 0.2243, T = 0.2523; AC = 1.4034, AG = 4.9417, AT = 0.9147, CG = 0.9256, CT = 5.9064, GT = 1.0; I = 0.3543; G = 2.0844). Asterisks indicate branches that receive a posterior probability of 95% or greater in the Bayesian analysis. MP bootstrap proportions >70% (above branches) and decay indices (below branches) are provided for all nodes congruent between analyses based on the two optimality criteria. A = Amphibolurinae; S = Scincinae. Other annotations as in Figure 1.

Many authors have documented the potentially misleading effects of heterogeneous base composition among taxa in phylogenetic studies (e.g., Lockhart et al., 1994; Steel et al., 1993; Tarrio et al., 2000; Tarrio et al., 2001 ). Chi-square tests for homogeneity of base frequencies on the RAG-1 data show that only third positions are significantly heterogeneous (P < 0.001). We therefore used Lockhart et al/s (1994) LogDet transformation (with an invariant sites parameter) to correct for base frequency heterogeneity (see Tarrio et al., 2001 for comments on this test) in a minimum evolution (ME) analysis of the full R/4G-7 data set. This analysis still finds Cekkonidae as the sister taxon of all other squamates (a Scincoidea-DibamHS clade is the next to diverge), and Iguania is still highly nested. Furthermore, this same basic topology is also found from MP and ML analyses using only first and second codon positions.

Harris (2003) found substantial heterogeneity in c-mos third codon positions among squamates, specifically high GC content among some teiids (the most basal ingroup taxon in his analysis) and his outgroup taxa. In our study, base frequencies at third positions of c-mos codons are not significantly heterogeneous (P = 0.1664 for all taxa, P = 1.000 for ingroup only). However, third-position CC content of our outgroup taxa (average 63.2%) is markedly higher than that of the ingroup taxa (average 41.5%), and Tarrio et al. (2000) suggested that this situation could cause incorrect rooting. However, ME analysis of LogDet-corrected c-mos data places Gekkonidae near the base of the tree and finds a highly nested Iguania, as do MP and ML analyses using only first and second codon positions. All of these topologies are similar to those derived from MP and ML analyses of the full c-mos data set (Fig. 2).

Heterogeneous selection pressure affecting the genes used for phylogenetic inference is another potentially confounding factor. If two or more nonsister lineages undergo a period of similar selection that is divergent with respect to the remaining lineages, then parallel or convergent amino acid replacements (in the case of a proteincoding gene) in these nonsister lineages could be problematic for phylogenetic analyses. However, K^sub a^/K^sub b^ values are remarkably uniform within and among clades for both the RAG-1 and c- mos genes. The test clades used for each gene include lguanidae, Acrodonta, Gekkonidae, Serpentes, Anguimorpha, Lacertiformes (including Amphisbaenia), and Scincoidca (including Xantusiidae). For both RAG-1 and c-wos, K^sub a^/K^sub s^ values for both within- clade comparisons (RAG-1 average 0.13 [0.12-0.15]; c-mos average 0.30 [0.23-0.54]) and between-clade comparisons (.RAG-I average 0.11 [0.10-0.12]; c-mos average 0.24 [0.16-0.33J) are significantly different from 1. For both data sets, all except within-clade snake comparisons are highly significant even after Bonferroni correction for multiple tests (Rice, 1989). These results indicate that both genes are under strong stabilizing selection, both within and among clades. K^sub a^/K^sub s^ values significantly greater than 1 between particular test clades or groups of test clades would indicate a shift in selective regimes, even if all within-clade values remained small (Messier and Stewart, 1997), and this would be a cause for concern. However, because all values were small, no evidence exists for heterogeneous selection that might mislead results of the nuclear analyses by causing convergence in protein structure between nonsister taxa.

Graham et al. (2002) showed that if the nearest outgroups are extremely divergent such that phylogenetic information is essentially randomized with respect to the ingroup taxa, trees will tend to root incorrectly on long terminal ingroup branches (this is an extreme example of LBA). However, even with uninformative characters excluded, average uncorrected Sphenodon-ingroup distance is only 0.358, which is well below the 0.75 expected from random DNA data. Average ML distances betwee\n ingroup taxa and turtles and crocodylids, respectively, are generally less than 5% higher than Sphenodon-ingroup distances (birds are somewhat more divergent). Analyses repeated with all possible outgroup combinations always give the same ingroup topology, with support values very similar to those obtained in the original analysis.

Finally, both Dibamus and Gekkonidae sit at the ends of relatively long branches in the RAG-1 analysis (Fig. 1). Therefore, even though the outgroup sequences are definitely not random with respect to ingroup sequences, LBA needs to be considered as a potential cause of our unconventional rooting. Although support is high for the RAG-1 branching order, LBA is related to the phenomenon of statistical inconsistency, and higher support values for an incorrect topology are predicted from theory as the amount of data increases (Felsenstein, 1978). One simple test for LBA is to remove the suspicious long branches and to see if the remaining topology is stable. We therefore sequentially removed (without replacement) the sister taxa of the remaining squamates from the RAG-2 data set (i.e., first Gekkonidae, then D/iamus, etc.; see Fig. 1) and then analyzed each of these modified data sets with both parsimony and likelihood. In each of these analyses, Iguania remained the most nested group. We also tried excluding Acrodonta from all these analyses, thus making the branch to lguanidae even longer, but the topology remained congruent with those from the original analyses.

Simulation studies likewise suggest that LBA is not responsible for our rooting. Analysis of simulated data sets modeled on the RAG- ? ML tree rerooted at Iguania never produced a rooting at Gekkonidae and /or Dibamus (Fig. 4A); this was true even when the branches to Iguania and Dibamus were artificially shortened and lengthened, respectively, to accentuate any tendency toward rooting at Dibamus (Fig. 4B). In all of these simulations, parsimony failed to recover the modeled topology more often than likelihood; however, this is because parsimony often incorrectly rooted the tree at Serpentes, another clade subtended by a relatively long branch (Fig. 1).

Analyses performed to localize the signal for the molecular and morphological rootings in the RAG-1 and morphological (Lee and Caldwell, 2000) data, respectively, found that support for the alternative roofings is not concentrated in one particular subset of either data set. Chi-square tests show that the morphological characters favoring the morphological rooting are randomly distributed among anatomical regions (Table 3), and separate analyses performed on sequence from the two broad functional domains of the RAG-1 gene likewise both recover the molecular rooting found with the full RAG-7 data.

FIGURE 4. Results from two sets of 100 parametric simulations designed to detect long-branch attraction that, if present, might have caused an incorrect rooting at geckos or Dibamus in the RAG-1 analyses. T+G = Teiidae/Gymnophthalmidae; L+A = Lacertidae/ Amphisbaenia; X+C = Xantusiidae/Cordylidae. (A) In the model tree for the first set, branch lengths are simply reoptimized on the RAC- I ML topology rerouted at Iguania. (B) The second model tree is identical, except that the branch leading to Iguania is multiplied by 0.33, and the branch to Oibamus is multiplied by 1.5 (see text). For clarity, outgroups are not shown. For both MP and ML analyses, the table shows the percentage of times the tree correctly rooted at Iguania as well as the number of incorrect roofings, both at the potential problem taxa and at various other points in the tree. In both sets of simulations, incorrect footings at Serpentes in the MP analyses accounted for most of the difference in accuracy compared to the ML analyses. As a control, a third set of 100 simulated data sets was constructed using as a model the exact RAG-1 ML topology, which is rooted at geckos (Fig. 1). For these data sets, MP and ML recover the modeled topology 57 and 87 times, respectively. Thirty- four of the 43 incorrect topologies recovered under MP root the tree at Oibamus.

Phylogenetic Positions of Limbless Taxa

Snakes, amphisbaeniarts, and Dibamus each are placed in separate parts of the tree, and alternative hypotheses placing any two of these as sister taxa are statistically rejected (Table 4). Analyses based on the two different optimality criteria disagree on the exact placement of snakes, although it is clear that they are nested well within squamates (Fig. 3). Parsimony places snakes as the sister taxon of Lacertiformes (including amphisbaenians) with weak support (bootstrap of 52). However, likelihood recovers a clade containing snakes, anguimorphs, and iguanians, and this arrangement is strongly supported by Bayesian results (Fig. 3). A sister-taxon relationship between snakes and Varanidae is statistically rejected (Table 4).

Inclusion of amphisbaenians within the traditional Lacertiformes (Lacertidae + Teioidea) is statistically supported (Table 4), and heuristic support is strong for a sister-taxon relationship between lacertids and amphisbaenians (Fig. 3). This latter relationship may be supported by a structural character as well. Gallotia (a lacertid) has a seven-codon deletion at c-mos positions 220 to 240, and all sampled amphisbaenians share an overlapping eight-cod on deletion at positions 217 to 240, suggesting that the original deletion was simply extended by one codon in amphisbaenians. Harris et al. (1999) reported a seven-codon deletion in this general region for two gekkonines. Although the alignment in this area is not completely unambiguous, alignments made with Clustal X (Thompson et al., 1997) at a variety of gap penalties ( see Materials and Methods) suggest that the lacertid and gekkonine deletions do not involve the same codon positions. Furthermore, forcing the gekkonine and lacertid deletions to coincide requires two separate, smaller amphisbaenian deletions instead of the one found at all Clustal gap- penalty settings used.

TABLE 3. Chi-square tests of the random distribution among character partitions of morphological characters(a) favoring the morphological rooting (at lguania) of the molecular topology(b) over the molecular rooting (at Dibninus) of the molecular topology. SS = 7.54; P = 0.479.

TABLE 4. Results of Wilcoxon signcd-ranks (Tcmpleton) and SH tests of topology.

Relationships within Amphisbaenia are strongly supported. The amphisbaenian family Rhincuridac is not represented in the combined RAG-1 and c-mos data set because of problems amplifying the c-mos fragment. However, a Templeton test performed on the RAG-1 data alone provides strong support (P < 0.0001) for monophyly of the other three amphisbaenian families to the exclusion of RAincura (Townscnd, 2002). Furthermore, in the combined analysis, Trogonophidae and Amphisbaenidae form a well supported clade exclusive of Bipedidae (Fig. 3).

In MP, ML, and Bayesian analyses, Dibamus is the sister taxon of a clade containing all other squamates, and geckos are the second group to diverge from the ancestral squamate lineage. Both parsimony and Bayesian measures strongly support grouping the remaining squamates to the exclusion of geckos and Diiamus (Fig. 3).

Other Well-Supported Clades from the Nuclear Analyses

Strong support is found for Xantusiidae as the sister taxon of Cordylidae, and for the placement of HeIodermatidae within a Xenosaurus-Anguidae clade to the exclusion of Varanidae, which is often considered the sister group of helodermatids (Fig. 3). Shinisaurus and Varanidae form a clade, and a sister-taxon relationship between xenosaurus and Shinisaurus (the traditional Xenosauridac) is statistically rejected (Table 4).

Within Scincidae, phylogenetic positions of the two limbless subfamilies are well supported. Acontinae is the sister group to a clade containing all other skinks, and Feylininae is closely related to African scincines (actually nested within this group; see Fig. 1). Monophyly of African and North American scincines is not supported (Fig. 3).

Relationships within Gekkonidae are well supported. Pygopodinae is the sister taxon of Diplodactylinae (Fig. 3), and this relationship is further supported by a shared one-codon deletion in the RAG-1 data set at positions 125 to 127. Tcrafoscmcus and Sphaerodactylinae form the sister group of Gekkoninae with high bootstrap and Bayesian support (Fig. 3). MP recovers Eublepharinae as the sister taxon of (TenzfoscmCMS 4- Sphaerodactylinae + Cekkoninae), but the bootstrap value is <70%. However, ML and Bayesian analyses recover this same relationship, and Bayesian support is high (Fig. 3). Furthermore, independent support for this arrangement comes from a shared fourcodon deletion at positions 95 to 106 in the KAG-I data set. No other sampled gekkonids have any deleted bases in this region, and the surrounding amino acid sequence is conserved across geckos, making alignment unambiguous.

Mitochondrial-DNA Analyses

The mtDNA MP strict consensus of four trees is unresolved at many deeper nodes, but resolved portions of the tree are largely compatible with the mtDNA ML topology (Fig. 5). The Bayesian consensus topology is similar to the ML topology, and both analyses recover all moderately to highly supported nodes (bootstrap >80%) from the parsimony analysis, with one exception within geckos (see below).

Chamaeleonidae, Agamidae, Acrodonta, Iguanidae, Anguimorpha, Serpentes, Scincidae, Amphisbaenia, and Teioidea all receive moderate to high parsimony bootstrap support, and all of these clades except Agamidae receive high Bayesian support (Fig. 5). Not all relationships within Amphisbaenia could be evaluated because mtDNA sequence could not be obtained from Rhineura. However, relationships among the remaining amphisbaenian families are strongly supported, and mirror exactly the results of the nuclear analysis (Fig. 3). Furthermore, the- lacertid-amph\isbaenian clade identified in RAG-1 and c-mos analyses is once again recovered with moderate parsimony and high Bayesian support (Fig. 5).

Within Gekkonidae, the mtDNA ML and Bayesian analyses find the same topology as all nuclear analyses, except that Eublepharinae (represented by Eublepharus turkmenicus) is the sister taxon of all other gekkonids. The mtDNA MP tree differs in finding moderate bootstrap support (81) for a sister-taxon relationship between Gekkoninae and Sphaerodactylinae, although this result is not supported by statistical tests (see Table 4). Interestingly, in the ML analysis the branches leading to gekkonines and sphaerodactylines are each roughly twice as long as the branch leading to Teratoscincus (Fig. 5). This result suggests that LBA may account for the mitochondrial parsimony gekkonine/sphaerodactyline clade, which is at odds with all other analyses in this study. Mitochondrial data agree with the nuclear data on the nesting of feylinine skinks within African scincines, and a sister-taxon relationship between xantusiids and cordylids.

In disagreement with the nuclear analyses, moderate (parsimony) to strong (Bayesian) support is found for a snake-acrodont clade (Fig. 5), although this result is not supported by the nonparametric statistical tests (Table 4). The branches subtending each of snakes and acrodonts are much longer than most other branches of similar depth in the tree (Fig. 5). Interestingly, both of these taxa have gene rearrangements associated with unusual positions of mitochondrial replication origins (Kumazawa and Nishida, 1995; Macey et al., 1997a), a situation that might be related to an increased rate of molecular evolution. The long snake and acrodont branches, combined with strong support for a monophyletic Iguania from the nuclear data (Fig. 3, Table 4), suggests that LBA might be causing this very unorthodox arrangement.

Results from parametric bootstrapping simulations support the LBA hypothesis (Fig. 6). When snakes and anguimorphs are constrained to be sister taxa (a more traditional scenario) in 100 simulated data sets, equalweights parsimony correctly recovers this clade in only 25% of replicates, whereas in 62% of replicates parsimony incorrectly recovers a snake-acrodont clade, as in the original analysis of the real mitochondrial data. In contrast, ML recovers the correct snake-anguimorph clade 78% of the time, and incorrectly recovers a snake-acrodont clade only 14% of the time.

In a second, more extreme deviation from the original mitochondrial-based topology, 100 data sets are simulated in which acrodonts and iguanids form a monophyletic Iguania as the sister taxon of a clade containing snakes and anguimorphs, a topology compatible with all well-supported nodes from the nuclear analysis (Fig. 3). Results from this analysis further support a role for LBA in the mitochondrial results. Equal-weights parsimony recovers the correct topology only 12% of the time, and 46% of the analyses incorrectly place snakes as the sister taxon of acrodonts. Meanwhile, likelihood recovers the correct topology 56% of the time, and incorrectly recovers a snake-acrodont clade in only 5% of the simulation replicates.

Average ML-corrected distances between ingroup taxa are nearly five times as large for the mitochondrial data as they are for the RAG-1 data (Table 2). Likewise, although the mitochondrial data set has only slightly more than half the number of characters of the RAG- 1 data set (and only one less species) (Table 2), the mitochondrial MP tree is approximately 45% longer than the corresponding RAG-1 tree (17237 and 9900 steps, respectively). Thus, the mitochondrial data are more likely to show saturation at more basal nodes, perhaps explaining the lack of resolution and poor support at deeper levels (Fig. 5). For this reason, nodes at the deepest levels of the squamate tree are probably best assessed from the nuclear data alone. However, the high levels of homoplasy in the mitochondrial data should not affect all levels of the tree equally; indeed, at more shallow levels, there is strong support for many clades. Although we know that the nuclear and mitochondrial data sets conflict in some areas (e.g., the placement of snakes), there is no reason to assume that these partitions are wholly incongruent, and we therefore combine all data sets for a final analysis. Bayesian runs were performed using one, two (nuclear and mitochondrial), and three (RAG-I, c-mos, and mitochondrial; Fig. 7) data partitions. Topologies were identical and support values were very similar for all partitioning schemes.

Monophyly of Agamidae has moderate (parsimony) to strong (Bayesian) support (Fig. 7), but is not supported by nonparametric tests (Table 4). Support is statistically significant (Table 4) for lacertid-amphisbaenian, trogonophid-amphisbaenid, and xantusiid- cordylid (SH test only) clades. As in the mitochondrial analysis, MP places snakes as the sister taxon of acrodonts (Fig. 7A). Interestingly, the ML analysis, which should be more resistant to LBA, instead places snakes as the sister group of (Anguimorpha + Iguania) (Fig. 7B). Monophyly of Iguania (Acrodonta + lguanidae) is supported by a posterior probability >95% in the Bayesian analysis, which agrees with well supported MP and Bayesian results from the nuclear analyses (Fig. 3).

FlGURE 5. Mitochondrial data, ML phylogram (GTR+I+G model; -InL = 63989.06; A = 0.4150, C = 0.3394, G = 0.0593, T = 0.1863; AC = 0.3588, AG = 2.5905, AT = 0.4988, CG = 0.3158, CT = 2.4827, GT = 1.0; I = 0.0951, G = 0.5652). Asterisks indicate branches that receive a posterior probability of 95% or greater in the Bayesian analysis. MP bootstrap proportions >70% (above branches) and decay indices (below branches) are provided for all nodes congruent between analyses based on the two optimality criteria. To facilitate discussion of potential LBA (see text), long internal branches leading to snakes and acrodonts are in bold, and pertinent higher taxa are labeled. Numbering as in Figure I.

FIGURE 6. Results from two sets of 100 parametric simulations designed to detect potential long-branch attraction in the mitochondrial analysis. For all simulations, the long-branch taxa (snakes and acrodonts) are separated in the model trees by making snakes the sister taxon of anguimorphs. (A) In the model tree for the first set (top), Iguania is not monophyletic. (B) In the model tree for the second set (bottom), Iguania is monophyletic (see text). For clarity, only the relevant portions of the model trees are shown; the remainder of each tree is identical to the original mitochondrial ML topology of Figure 5. The tables show the number of times incorrect and correct topologies are recovered under MP and ML, respectively. As a control, 100 simulated data sets were constructed using as a model the exact ML topology of Figure 5, in which snakes and acrodonts are sister taxa. For these control data sets, MP and ML each recover the model topology in 98% of the analyses.

Figure 8 summarizes molecular support for higherlevel phylogenetic relationships within Squamata. In this figure, we have not constructed a consensus of all trees recovered from the three data sets used in this study. Rather, we present all nodes receiving both strong (>95%) bootstrap and Bayesian support from either the combined nuclear or the mitochondrial data set, and which are not contradicted with similar levels of support in the other data set. All nodes not meeting the above criteria are collapsed, thus giving a conservative estimate of well-supported squamate relationships.

DISCUSSION

Paraphyly of Scleroglossa and Evolution of Squamate Feeding

The most important discrepancy between our results and the morphological hypotheses is our strong statistical rejection of the hypothesis that taxa traditionally included in Scleroglossa form a monophyletic group. This grouping (though not the taxon name) dates to Camp's (1923) study, and it is supported by numerous osteological and soft-tissue characters (Estes et al., 1988; Schwenk, 1988 and subsequent authors), as well as behavioral characters related to prey prehension (Schwenk and Throckmorton, 1989). The basal branches of the ingroup are all relatively small for all genes. However, although we do not specifically test the cost of moving the root across each of these branches individually, our results clearly reject its placement a t Iguania, as congruence with the morphology would require.

FIGURE 7. Results of combined R/lG-2/omos/mtDNA analyses. To highlight its change in position under the two different oplimality criteria, the snake clade is in bold. Other relevant higher taxa are labeled. A = Amphibolurinae; S = Scincinae. (A) Strict consensus of four most parsimonious trees (L = 22234, RI = 0.415). MP bootstrap proportions >70% are shown above branches and decay indices are in bold below branches. (B) ML phylogram (GTR+I+G model; -InL = 93061.53458; A = 0.3456, C = 0.2690, G = 0.1506, T = 0.2348; AC = 1.7122, AG = 3.6083, AT = 0.9770, CG = 0.6230, CT = 5.0995, GT = 1.0; I = 0.2362, G = 0.7308). Asterisks indicate branches that receive a posterior probability of 95% or greater in the Bayesian analysis.

Our results suggest reinterpretation of studies that have used comparative methodology to contrast ScIeroglossa and Iguania. For example, Schwenk (1993) found a fundamental difference in tongue morphology and prey-prehension technique between iguanian (lingual prehension) and scleroglossan (jaw prehension) lizards. Schwenk (1986) reported that the tongue of Sphenodon (a lingual feeder) shares many features with iguanid lizards, including muscle-fiber architecture and hyobranchial-foretongue coupling. Based on these similarities, along with independent evidence for a basal dichotomy within squamates between Iguania and ScIeroglossa (Estes et al., 1988), Schwenk (1986) concluded that Sphenodon and iguanians e\xhibit the ancestral squamate (and lepidosaurian) condition. This inference of the ancestral condition is problematic, however, because jaw prehension is widespread in the closest outgroups to lepidosaurs (birds, turtles, and crocodilians). Schwenk (1989) cites examples of lingual prehension in some of these groups as evidence that it is the ancestral state; however, given the difficulty in comparing the highly modified feeding apparatus between these distantly related groups (Schwenk, 1988), this conclusion may be unwarranted.

Under Schwenk's (1986) scenario, the common ancestor to Scleroglossa evolved a fundamentally different feeding system and associated tongue morphology. Nonherbivorous iguanians are generally considered ambush predators with little ability to detect chemical cues from prey items (e.g., Cooper, 1995), whereas scleroglossans are often actively foraging lizards that tongue-flick to collect chemical cues from prey items (although several exceptions exist; see Perry, 1999). Release of the tongue from its prey-prehension duties is thought to have allowed this new role to evolve, whereas the functional constraints imposed by lingual prey-prehension presumably have prevented most iguanians and Sphenodon from developing olfactory capabilities to the same extent (Schwenk, 1993).

FIGURE 8. Summary of higher-level squamate phylogenetic relationships well supported by molecular data. Branches with any type of bar are supported by MP bootstraps and Bayesian posterior probabilities ≥95% in the combined nuclear analysis. Solid bars denote branches also supported by bootstraps and posterior probabilities ≥95% in the mitochondrial analysis. Hatched bars denote branches supported by posterior probabilities (but not bootstraps) ≥95% in the mitochondrial analysis. Open bars denote branches not congruent with any mitochondrial topology but which are also not strongly contradicted (by the above support criteria) in the mitochondrial analyses. Note that Rltineurn (Amphisbaenia) was not included in the analyses from which this figure was derived.

However, even if lingual prehension is assumed to be the ancestral lepidosaurian condition, it is possible that the similar feeding behavior and tongue morphology of Sphenodon and iguanians represent homoplasy rather than homology. As mentioned by Schwenk (1986), several authors (e.g., Cans, 1983; see also Wu, 1994) have noted that Sphenodon is not a basal but rather a highly nested taxon within Rhynchocephalia, a once widespread group that included a diversity of body plans and lifestyles, including long-legged terrestrial forms, long-bodied obligately aquatic forms, and specialized herbivores (Evans et al., 2001; Reynoso, 2000). Although Sphenodon is almost certainly the closest living relative to squamates, considering its character states ancestral for Squama ta, especially when the characters involve largely soft-tissue anatomy and behavior, is problematic.

FIGURE 9. Evolution of feeding mode in lepidosaurs. White branches indicate lingual prehension, black branches jaw prehension, and the gray branch is equivocal. (A) Under a traditional monophyletic Scleroglossa, produced here by simply rerouting our nuclear topology, lingualfeeding arose once in a common ancestor of Sphenodon and squamates, and was lost in an ancestor to Scleroglossa. Note that we have adopted Schwenk's (1986) assumption that lingual prehension is ancestral for lepidosaurs (but see text). (B) Under our nuclear topology, feeding mode is more labile, with lingual feeding arising at least twice, once either in the lineage leading to Spheiiodim or in a common ancestor of Splienodan and squamates (allowing uncertainty in the outgroup designations), and once in an ancestor to Iguania. Dotted lines indicate ambiguity in the position of the first acquisition of lingual feeding. Data from Schwenk (2000).

Schwenk and Wagner (2001) used suites of characters associated with both lingual-prehension and jawprehension modes of feeding to illustrate their cvolutionarily stable configuration (ESC) concept, arguing that the phylogcnetic stability of lingual feeding across a variety of habitats and lifestyles is evidence of a complex, integrated system. Internal selection for maintenance of the entire functional system results in only rare transitions from one system to another. In the example discussed here, the strong interdependence among components of the lingual-prehension feeding mode is thought to have led to its persistence in virtually all iguanians, regardless of habitat, diet, or other ecological variation. Only when the components of this system were somehow decoupled in the common ancestor to scleroglossans could jaw prehension and its associated olfactory and behavioral traits evolve.

Under our phylogenetic hypothesis, iguanians and Sphenodon (or some possibly distant ancestor to Sphen-odon) are inferred to have acquired lingual prey-prehension techniques independently (Pig. 9). Beca-use food prehension techniques, tongue musculature, and chemosensory ability are unknown for rhyncho-cephalians other than Sphenodon, this scenario is only slightly less parsimonious than the traditional view. Al-though similarities in muscle fiber and connective-tissue architecture between Sphenodon and iguanians may be explained most parsimoniously by symplesiomorphy (Schwenk, 1986), if lepidosaurian feeding systems truly are highly integrated and constrained, tongue morphol-ogy could evolve to be markedly similar in unrelated groups adopting the same feeding mode. Indeed, the ap- parent lability of feeding mode and tongue morphology is exemplified in the recovery (with each of the three gene regions) of an iguanian- anguimorph-snake clade, a group that represents the extremes of these traits within squamates.

Phylogenetics of Limblessness in Squamates

Our nuclear data stati

Source: Systematic Biology

More News in this Category