Squamates (lizards, snakes, and their kin such as amphisbaenians, or “worm lizards”) represent the world's most diverse clade of terrestrial vertebrates with ∼11,000 described extant species, representing key components in many of the world's most diverse ecosystems. With an evolutionary history dating back at least to the Middle Triassic at 242 Ma, the squamate Tree of Life also features numerous diverse but extinct branches, with hundreds of fossil species found all over the world. Despite their biological relevance both today and in the geological past, there remains a centuries-old controversy on how the major lineages of squamates are related to each other, with a direct impact on studies in ecology, evolution, paleontology, toxinology, and other fields. Here, we provide a historical overview of this long research tradition, from 19th century naturalists to 21st century phylogenomics, with special emphasis on several recent advances over the last two decades. These insights have had a dramatic effect on our understanding of the squamate Tree of Life and clarify several possible future research agendas. We provide an integrative perspective derived from genomics, morphology, and the fossil record and propose several points of synthesis in our current knowledge of broadscale squamate evolution and systematics. Key topics of interest include dating the origin and early evolution of lizards, the phylogenetic origin of snakes, the evolution of venom, recent agreements between morphological and molecular squamate evolutionary trees, genomic patterns of evolution, and the integration of morphological and molecular data sets. We conclude by providing perspectives on possible advancements in the field, directing researchers to promising future lines of investigation that are necessary to further expand our synthetic knowledge of squamate evolution.
Squamates (lizards, snakes and amphisbaenians) are the most diverse known clade of extant terrestrial vertebrates with ca. 11,000 species, nearly twice as that of mammals (Uetz and Hošek, 2021). With an evolutionary history spanning more than 242 million years (Simões et al., 2018), squamates have inhabited nearly all terrestrial and aquatic environments from the distant past to today, including all of the world's oceans during the Cretaceous, and from Antarctica to the Arctic Circle during warmer periods (Pianka and Vitt, 2003; Vitt and Caldwell, 2013). An accurate understanding of the squamate Tree of Life is thus a crucial basis for reconstructing the evolutionary history of a considerable component of the world's biodiversity. A robust phylogenetic tree is therefore the necessary starting point towards a broader understanding of the evolvability of functional traits (e.g., Wiens et al., 2006), distribution of extinction risk (e.g., Tonini et al., 2016), origin of life history strategies (e.g., Pyron and Burbrink, 2014), and myriad other macroevolutionary dynamics (e.g., Title and Rabosky, 2017).
Although we are now in the fourth decade of quantitative phylogenetic analyses of squamates (beginning with Estes et al., 1988; see below) and their various subgroups, numerous questions remain unanswered. Morphological data sets differ strongly in the placement of key fossil groups such as Mosasauria and Paramacellodidae (see Pyron, 2017) and relationships among the major extant families. Molecular data sets differ in the placement of key groups such as Dibamidae; the relationships of Iguania, Anguimorpha, and Serpentes; and relationships within Iguania (see Burbrink et al., 2020). Perhaps most famously, many morphological and molecular data sets appear fundamentally irreconcilable with respect to the placement of taxa such as Iguania (Losos et al., 2012). Although recent genomic and morphological data sets resolve some of these issues (Simões et al., 2018; Singhal et al., 2021), numerous branches are still unresolved. These lingering questions about the squamate Tree of Life have affected our ability to answer fundamental questions about squamate systematics and evolutionary biology.
First, what is the best phylogenetic hypothesis for Squamata, including both extinct and extant lineages? Are the various traditional morphological hypotheses driven by widespread homoplasy (e.g., Losos et al., 2012)? Are evolutionary processes such as paralogy, incomplete lineage sorting, and horizontal gene transfer or hybridization driving conflicting topologies in genome-scale data (e.g., Siu-Ting et al., 2019)? Could other sources of morphological (e.g., soft tissue; Siegel et al., 2011) or molecular (e.g., synteny mapping; Sevillya and Snir, 2019) data offer alternative perspectives on species relationships? Could both morphological and molecular estimates be confounded by analytical problems such as gene tree error (Roch and Warnow, 2015) or character construction and coding (Simões et al., 2017b)? Can refinement of characters (Simões et al., 2018) and integration of ontogeny (Wiens et al., 2005) help to alleviate the apparent molecular-morphological discordance?
Second, once a higher level resolution of the squamate Tree of Life is achieved, what can we learn about the evolutionary history of snakes, lizards, and worm-lizards? What are the earliest evolving lineages of squamates and their closest reptile relatives? When did crown squamates originate? Do snakes have a marine, terrestrial, or fossorial ancestry (Vidal and Hedges, 2004; Lee, 2005a)? Did venom originate only once or multiple times across lizards and snakes (e.g., Fry et al., 2006)? The questions that such an estimate can help answer are nearly endless.
Here, we review the recent history and future directions of research into the squamate Tree of Life. We detail many of the most contentious relationships, the data and methods used to arrive at various results, and the conclusions reached by recent authors regarding phylogenetic topologies and evolutionary histories in the group. We then offer a synthesis and perspective on where knowledge currently stands for these contentious branches. Finally, we detail potentially fruitful avenues for future investigations. We believe that promise exists in both data and methods to drastically increase our knowledge of lizard and snake relationships and that these will provide grounds for an increasingly broadscale synthetic view of squamate evolution.
A BRIEF OVERVIEW OF SQUAMATE DIVERSITY
Extant Lizards and Snakes
The phylogenetic diversity and taxonomic history of extant squamates (as well as extinct lineages; see below) is staggeringly complex, with a bewildering array of clade names, definitions, and transfers of lineages between taxa as phylogenetic understanding increased (see Conrad, 2008; Vitt and Caldwell, 2013). Recounting these names and their evolving composition over time is too substantial an undertaking to present here. Instead, we will briefly outline the major clade names and their composition from the recent study of Burbrink et al. (2020), who analyzed higher level relationships within Squamata with genome-scale data. Those authors provide a baseline understanding of phylogeny and taxonomy for a series of higher level clades for which consensus has emerged since the formalization of the “molecular” hypothesis of squamate relationships (see Vidal and Hedges, 2005, 2009; Pyron et al. 2013). Many of these taxa (e.g., Serpentes, the snakes) have been recognized in some form or another since initial historical classifications were first erected. We also comment on the diversity and composition of each group (see Vitt and Caldwell, 2013; Uetz and Hošek, 2021).
Although species, genera, and even entire families and subfamilies have been moved into or out of these taxa based on morphological and molecular analyses over the years, we present them here as defined by most recent molecular phylogenies (e.g., Vidal and Hedges, 2005; Pyron et al., 2013; Burbrink et al., 2020). Serpentes is perhaps the most well characterized, as few extant snakes have ever been misclassified as lizards or vice versa. This group was historically classified (see below) into the blindsnakes and threadsnakes (Scolecophidia); boas, pythons, and allies like pipe-snakes, sunbeam snakes, and shieldtails (Henophidia); and advanced snakes (Caenophidia), comprising many of the most common species worldwide and all medically significant venomous species. Scolecophidia and Henophidia are most likely not monophyletic, although Caenophidia and Caenophidia+“Henophidia” (Alethinophidia) are strongly supported in most recent phylogenies (see Burbrink et al., 2020; Singhal et al., 2021). Snakes represent ∼3,900 species (Uetz and Hošek, 2021) occurring on every continent except Antarctica and in nearly every terrestrial and aquatic environment, spanning a massive range of morphological and ecological diversity.
Iguania is a diverse clade of ∼2,000 species of primarily diurnal, terrestrial and arboreal lizards including charismatic and well-known taxa such as Anolis, chameleons and agamas, iguanas, and other diverse assemblages. Whereas most iguanian families (Pleurodonta) inhabit Neotropical regions in the Americas, acrodontans (agamids and chameleons) inhabit Old World continents. The exception to that are some pleurodont families with species in Madagascar (Opluridae) and the Galapagos and Fiji (Iguanidae), which give puzzling indicators of either Gondwanan vicariant origin or long-distance oceanic dispersal (Noonan and Chippindale, 2006; Noonan and Sites, 2010). Surprisingly, the most abundant fossil record of extinct pleurodontan iguanians occurs in East Asia (e.g., Gao and Norell, 2000) whereas the oldest known acrodontans came from what is today South America and Africa (Simões et al., 2015b; Apesteguía et al., 2016), further confounding our understanding of their biogeographic history.
Several apparently “primitive” or plesiomorphic morphological traits, such as a fleshy, lobular tongue used for lingual prey prehension, led historical authors and early morphological systematists to place Iguania as the earliest diverging lineage of squamate (see below; Losos et al., 2012). Relationships among the many pleurodontan iguanian families have proven intractable to estimate even from genome-scale data sets, apparently because of a rapid evolutionary radiation with short intervening branch lengths (Townsend et al., 2011; Singhal et al., 2021).
Another group of widely recognized lizards has often been referred to as Anguimorpha (see below), but was referred to as Anguiformes (Conrad, 2006) by Burbrink et al. (2020). This includes the famous monitors (Varanidae) such as the Komodo dragon; the venomous Gila monster and beaded lizards (Helodermatidae); the Mexican knob-scaled lizards (Xenosauridae); the Chinese crocodile lizard (Shinisauridae); the American legless lizards (Anniellidae); the galliwasps (Diploglossidae); and the alligator lizards, glass lizards, and slow worms (Anguidae). Despite a relatively small number of species (∼250) compared with other diverse squamate clades, anguiforms nonetheless contain a number of exceptionally distinct species that are ingrained in the public consciousness. The group is primarily tropical in distribution, with only a few species in temperate regions, but occurring on every continent except Antarctica. Some of the oldest known records of the group come from the Late Cretaceous of North America and the Early Cretaceous of the UK—represented by Dorsetisaurus (Hoffstetter, 1967; Prothero and Estes, 1980). Anguids in particular become more abundant during the Eocene in North America, including the now extinct and heavily armored Glyptosaurinae (Sullivan, 1979; Estes, 1983).
Scincoidea or Scinciformata (see Vidal and Hedges 2005, 2009; Pyron et al. 2013) represents a distinct clade of terrestrial, arboreal, or fossorial lizards, typically presenting hard shiny dorsal scales or forms of plate armor, in the families Scincidae (skinks), Gerrhosauridae (plated lizards), Cordylidae (girdled lizards), and Xantusiidae (night lizards). Gerrhosaurids, cordylids, and xantusiids contain only a few dozen species each and are restricted to southern Africa (Cordylidae), Africa and Madagascar (Gerrhosauridae), or the desert southwest of North America, Cuba, and tropical Central America (Xantusiidae). In contrast, skinks are hyperdiverse (∼1,700 species) and distributed worldwide in nearly every region and community containing squamates, with exceptional ecomorphological diversity and habitat variability. Fossils of scincoids are relatively rare in the fossil record, but early scincoids (or close relatives) may comprise one of the first squamate clades to show up in the fossil record (see below—Paramacellodidae).
Lacertoids or Laterata represents a diverse group (∼1,000 species) of primarily terrestrial, diurnal, highly active, and conspicuous lizards, along with an enigmatic clade of burrowing, limb-reduced species. Lacertidae (wall or true lizards) are a common feature of squamate communities across Europe, Africa, the Middle East, and Asia. Gymnophthalmidae (spectacled lizards) are ecologically widespread in Central and South America, including savannas, deserts, and montane and lowland rainforests, and are typically smaller, more nocturnal, and more cryptic or fossorial than the other terrestrial lacertoids. Teiids are a widespread component of typically arid- or semiarid-adapted squamate assemblages across North and South America, but including some rainforest and semiaquatic species, ranging from small (∼10 cm total length) to large (.1 m) body sizes.
Nested within Lacertoidea as the likely sister lineage to Lacertidae (see below), are the enigmatic Amphisbaenia, the “worm lizards.” These are limb-reduced or limbless species, usually with stout bodies and heads adapted for a fossorial burrowing lifestyle, with a highly fragmented and disjunct distribution across Florida, Central America, the Caribbean, northern South America, sub-Saharan Africa, the Iberian Peninsula and adjacent northern Africa, Madagascar, and the Arabian Peninsula (Vitt and Caldwell, 2013). Their phylogenetic placement has been contentious across various morphological and molecular analyses through time (see below). The early fossil record of amphisbaenians includes extremely well preserved Eocene and Oligocene rhineurids within the genus Spathorhynchus (Berman, 1973, 1977; Müller et al., 2016). However, the fossil record of amphisbaenians before the Eocene includes several purported intermediate forms that have subsequently been reassigned to other clades by morphological phylogenetic analyses, such as Sineoamphisbaena (Wu et al., 1993; Kearney, 2003; Gauthier et al., 2012) and Slavoia (Conrad, 2008; Tałanda, 2016; Simões et al., 2018), both from the Late Cretaceous of Mongolia.
Gekkota represents a diverse lineage of highly distinct and well-known species, consisting of more than 2,000 species. Gekkotans are classified into several extant families (Diplodactylidae, Carphodactylidae, Eublepharidae, Sphaerodactylidae, Phyllodactylidae, and Gekkonidae), typically distinguished from all other squamates by their reduced irregular granular scales and expanded toe pads that allow many species to scale smooth vertical surfaces and even hang upside down by Van Der Waals forces generated by subdigital lamellae (Autumn et al., 2002). Females typically only lay 1–2 eggs per clutch, in contrast to the wide range of fecundity observed in other squamate lineages. Gecko species are found worldwide, with their highest diversities in tropical and warm desert regions. Nested within geckos are the flap-footed lizards (Pygopodidae), a small clade (∼35 species) of limb-reduced or limbless forms of highly specialized geckos found in Australia and New Guinea.
In both morphological and molecular phylogenies, Gekkota is typically estimated as an early-diverging lineage, subject only to the placement of Iguania and Dibamia. The gekkotan fossil record is relatively scarce, with the most informative fossils coming from the Late Cretaceous of Mongolia (Daza et al., 2014), but it includes the oldest adhesive toepads preserved in amber, indicating this remarkable adaptation had already evolved at least by 100 Ma (Arnold and Poinar, 2008). Stem gekkotans also display adaptations for scansoriality, despite not yet having adhesive toepads on the basis of their hand morphology (Daza et al., 2014; Simões et al., 2017a).
Finally, Dibamia (Dibamidae; blind skinks) is a small group (∼25 species) of highly derived limb-reduced or limbless fossorial species found primarily in mesic temperate or tropical habitats from Southeast Asia to western New Guinea (Dibamus), with a single species (Anelytropsis) in the deserts of northeastern Mexico. They are primarily insectivorous, with reduced eyes, and typically lay a single egg with a hard calcified shell, unlike the leathery eggs of most other squamates. Their reduced body form, with convergence on a burrowing body shape and associated adaptations in skull morphology, makes their phylogenetic affiliations uncertain by morphological data (see below). Most molecular analyses place them as an early-branching lineage of Squamata, either the earliest diverging lineage, the sister lineage of Gekkota (with Gekkota+Dibamia forming the earliest diverging squamate lineage), or as the sister lineage to all other squamates excluding Gekkota. As with Iguania, rapid evolutionary radiations with short internodes and long descendant branches make decisive resolution of these early squamate nodes difficult, even with genome-scale data (Singhal et al., 2021). The complete absence of a fossil record of dibamids makes the matter of understanding the origin of dibamids even more complex, because it prevents access to the early morphology of the group and calibration points for its origin.
Key Fossil Taxa
Fossils are known throughout the entire evolutionary history of Squamata, from early stem taxa in the Triassic (Figs. 1A–C) to Pleistocene ancestors of extant species. (For broad previous reviews, see Estes, 1983; Evans, 2003; Caldwell, 2019.) Many extinct species are easily assignable to crown lineages such as Serpentes (snakes) (Figs. 1D, E) or Gekkota (geckos) and provide clear insight into the evolutionary history of those groups (Daza et al., 2014; Caldwell, 2019). However, several diverse enigmatic groups of extinct species form entirely distinct evolutionary branches within crown squamates, but of uncertain phylogenetic affinity to extant lineages. Among the several fossil groups of squamates, we here provide a stronger focus on three clades that have had a large effect on reconstructing the squamate Tree of Life, either because of their age or taxonomic richness: Paramacellodidae, Borioteiioidea (or Polyglyphanodontidae), and Mosasauria.
Paramacellodids (Figs. 1F–I) are perhaps the least understood groups of fossil squamates. They primarily appear to be terrestrial lizards but have an extremely incomplete fossil record, mostly composed of fragmentary and disarticulated fossil remains that complicates any further inference of ecomorphology. However, this poor fossil record is a consequence of the age of the group: They were mostly abundant between the Middle Jurassic and the end of the Early Cretaceous (Bittencourt et al., 2020), an age from which the terrestrial vertebrate fossil record is scarce globally (Close et al., 2019) and furthered biased against small-bodied vertebrates (Evans, 2003). On the other hand, the age of paramacellodids makes them an incredibly important group, because they represent the oldest recognizable clade in the evolutionary history of squamates to diversify and occupy multiple continents (Bittencourt et al., 2020). Their fossils have been found in the Middle Jurassic of Europe and possibly Morocco; the Late Jurassic of North America, Europe, Africa, Central and East Asia (China); the Early Cretaceous of Europe, Africa, North America, South America, Central Asia (Siberia), and East Asia (Japan and Mongolia); and the Late Cretaceous of Europe (Waldman and Evans, 1994; Nydam and Cifelli, 2002; Folie and Codrea, 2005; Haddoumi et al., 2016; Bittencourt et al., 2020).
Owing to the fragmentary nature of most fossils, the primary diagnostic feature uniting paramacellodids lies in their dentition: teeth labiolingually expanded at their bases, along with lingually concave tooth apices (Figs. 1G–I). These characteristics for dental morphologies are relatively rare within squamates and only found in conjunction in paramacellodids. Another highly relevant phylogenetic character from this lineage is the frequent presence of body osteoderms (Fig. 1F). This osteological feature has been one of the most important characters for placing paramacellodids as closely allied to scincoids in most phylogenies (see below). It is likely that the discovery of more complete specimens of paramacellodids, or the redescription of poorly known but fairly complete taxa (e.g., Sharovisaurus; Hecht and Hecht, 1984), will be key to illuminate the placement of this group in the squamate Tree of Life.
From the Late Cretaceous, another taxonomically and geographically diverse group of squamates has drawn a lot of attention for reconstructing the squamate Tree of Life. This group has been named and rediagnosed several times in recent decades, including Polyglyphanodontinae (Estes, 1983), Borioteiioidea (Nydam et al., 2007), Polyglyphanodontidae (Conrad, 2008), and Polyglyphanodontia (Gauthier et al., 2012) (Figs. 1J–L), as more species have been described and anatomical knowledge has increased. Historically, borioteiioids have typically been linked to teiioids, initially as a subfamily of Teiidae (Estes, 1983) on the basis of gross morphological similarity. Later estimates from quantitative phylogenetic analyses place the group as the sister lineage to Teiidae (Conrad, 2008) or Teiioidea (Nydam et al., 2007; Simões et al., 2018) on the basis of several coded similarities between their dentition and that of teiioids, such as the heavy deposits of cementum at tooth bases, the presence of labiolingually expanded teeth, and the semicircular shape of the resorption pits (Nydam et al., 2007; but see below for a distinct hypothesis proposed by Gauthier et al., 2012).
Borioteiioids were geographically widespread throughout North America and East Asia during the Late Cretaceous and may be an ancestral lineage involved in the dispersal of teiioids to South America (Nydam et al., 2007). Their known fossils include individuals with a wide variety of dental morphologies, including species with a molarlike multicuspid dentition resembling that of mammals, as found in Peneteius (Nydam et al., 2000). Another example of extreme adaptation in the group is found in Polyglyphanodon sternbergi, by far the most completely known borioteiioid from North America with nearly perfectly articulated individuals (Figs. 1J, K), including an ontogenetic series and possible evidence of sexual dimorphism (Simões et al., 2016). Polyglyphanodon had an extreme level of labiolingual expansion of its dentition, with tiny denticles on the occlusal surface that were apparently highly adapted for cropping vegetation (Nydam and Cifelli, 2005). Polyglyphanodon and the Chinese borioteiioid Tianyusaurus further possessed an extremely unusual feature among squamates—the presence of a complete lower temporal bar forming a double temporal fenestration of the skull (Figs. 1l, 2, char. 3). This complete lower temporal bar is a condition reminiscent of early diapsid reptiles, such as archosaurs, that was reacquired at least twice independently within borioteiioids, but never in any known other group of living or extinct squamate (Simões et al., 2016). Owing to the completeness of several specimens of borioteiioids from North America and Asia, particularly Polyglyphanodon and Gilmoreteius, borioteiioids represent one of the very few clades of extinct squamates for which their morphology is extremely well characterized and generally well understood.
The third and final group of interest for the discussions below is Mosasauria (Figs. 1M, N), a diverse group (∼80 valid species) of mostly marine lizards with an extraordinary range in size, from ∼5 cm in length in the dolichosaurid Aphanizocnemus libanensis (Dal Sasso and Pinna, 1997) up to an estimated 17 m in Mosasaurus hoffmannii (Lingham-Soliar, 1995)—a size range spanning at least two orders of magnitude. Mosasaurians include early evolving lineages with only a few described species, such as Dolichosauridae (e.g., Dolichosaurus and Pontosaurus) and Aigialosauridae (e.g., Aigialosaurus) (DeBraga and Carroll, 1993; Campbell Mekarski et al., 2019). These taxa are generally considered to be the forerunners of mosasaurs—the giant aquatic lizards that occupied all of the world's oceans by the end of the Cretaceous (Polcyn et al., 2014; Jiménez-Huidobro et al., 2017). Mosasaurs (Mosasauridae) have traditionally been subdivided into two main groups, the Mosasaurinae and Russellosaurina (Bell, 1997). The first phylogenetic analysis of the group using Bayesian inference suggested the latter might not be monophyletic (Simões et al., 2017d), but it was later sustained by relaxed clocks (Madzia and Cau, 2020). Interestingly, several recent phylogenetic analyses of mosasaur relationships suggest an independent acquisition of aquatic adaptations in the group, with at least two independent land-to-water evolutionary transitions (Caldwell and Palci, 2007; Simões et al., 2017d).
Mosasaurians have played a central role in controversies over morphological analysis of squamate relationships, specifically owing to their proposed sister group relationship to snakes, or as an early branching lineage of squamates (see below). Of special interest to broader squamate relationships is a better understanding of the morphology and phylogeny of early mosasaurians that, historically, have been far less studied than mosasaurids. Some key taxa in that regard are the dolichosaurids Pontosaurus and Adriosaurus (Lee and Caldwell, 2000; Pierce and Caldwell, 2004; Caldwell, 2006), which along with the aigialosaurid Aigialosaurus (DeBraga and Carroll, 1993; Caldwell and Dutchak, 2006) represent the best preserved early mosasaurians. Deeper insights into those specimens from high-resolution computed tomography scanning, hopefully supplemented by discoveries of new well-preserved specimens, will play a key role in resolving the placement of mosasaurs.
The Precladistic Era
Attempts to understand the squamate Tree of Life have a long and complex history, including simple pre-evolutionary classifications in the 18th century from gross anatomy, complex morphologically based taxonomies informed by evolutionary theory in the 19th century, proto-cladistic graphical representations of phylogeny in the early 20th century, and statistical inferences in recent decades from morphological and molecular data sets. In the pre-evolutionary era, we can mark the start of squamate systematics (as with most animal groups) by the work of Linnaeus (1758), with Laurenti (1768), Daudin (1802), Oppel (1811), Cuvier (1817), Duméril and Bibron (1839), Fitzinger (1843), and Gray (1845) expanding our knowledge of species richness and erecting increasingly sophisticated classifications on the basis primarily of morphological features, both internal and external. In the evolutionary era, Cope (1864), Boulenger (1885), Nopcsa (1923), and many others began to expand and refine these classifications in an evolutionary and quasi-phylogenetic approach.
These early works are too numerous and expansive in scope, and their varying taxonomic concepts and nomenclatural changes too complex to summarize easily here. Yet, from them, we derive many of the clade names and natural groupings that are still recognized today. The extremely detailed anatomical investigations of these early authors quickly led to the recognition of robustly supported “clades” defined by morphological “synapomorphies,” even if these concepts were not yet developed or fixed in the minds of investigators at the time. Thus, many genera, subfamilies, and families were erected on the basis of clear anatomical affinities in the 18th and 19th centuries, and many of these hypotheses have stood the test of time as both morphological and molecular systematics have increased in precision and accuracy. In contrast, processes such as ecomorphological convergence, particularly for limb-reduced and burrowing body forms, made the early allocation of many such taxa difficult and uncertain. Similarly, hyperdiverse groups without clearly hierarchical diagnostic morphological characters (e.g., colubroid snakes, geckos, and skinks) proved challenging to classify in a coherent way below the “family” or “superfamily” level. Thus, compelling resolution of many branches has only been possible in recent decades with molecular data, yielding many phylogenetic novelties (see Pyron et al., 2013). Similarly, clade names such as “Sauria” for lizards persist even in the present day, despite the fact that Sauria has long been known to be paraphyletic with respect to Serpentes and Amphisbaenia; both snakes and worm lizards are simply derived lizards.
Two major syntheses in the 20th century (Camp, 1923; Underwood, 1967) represent the first systematic attempts to unify available anatomical data and classify lizards and snakes in a comparative framework with explicit notions of common ancestry and ancestor-descendant relationships among groups. Whereas the latter postdates the English edition of Hennig (1966) by one year, both are precladistic in method and represent the culmination of premodern thought in squamate systematics on the basis of morphological examination. Taxa were grouped and arranged by an implicitly subjective understanding of shared characters, but generally without formal recognition of this principle or quantitative analysis thereof. Camp (1923) is by far the most influential of all of the precladistic phylogenetic studies, depicting several aspects of squamate relationships that would dominate morphology-based research in squamate evolution for the next century. For instance, he proposed an early dichotomy between a group formed by iguanians and geckos (Ascalabota) and all other squamates, followed by a dichotomy between Ascalabota and autarchoglossans (composed of Scincomorpha and Anguimorpha)—many names and clade compositions that would be later formalized by the first statistical phylogenies that used morphological data (see more below). Later authors used such essays and the morphological characters illustrated in them as the basis for cladistic analyses of the relationships within groups such as “primitive” snakes (Rieppel, 1979) and geckos (Kluge, 1987).
A revolutionary classification of snakes introduced by Underwood (1967) was based on extensive anatomical documentation of a nearly exhaustive accounting of both soft tissue and osteological characters. In doing so, he provided the first systematic overview of snake diversity since the 19th century works of Cope and Boulenger. He documented characteristics of systematic relevance in extraordinary detail and breadth, providing diagnoses or character-based descriptions for most taxa down to the subfamily level, and in many cases to the tribe or genus. Furthermore, he explicitly recognized the link between such taxonomic arrangements, the characters on which they are based, and the evolutionary history and series of evolutionary transitions they implied. The enormity of this work and its towering importance for nearly all subsequent studies of snake classification are difficult to overstate, but we will here only note that it formalized the tripartite infraordinal arrangement of snakes into the “lower” snakes, the blindsnakes and allies (Scolecophidia) and the boas, pythons, and allies (Henophidia), and the “higher” or advanced snakes (Caenophidia), such as colubrids, dipsadids, elapids, and viperids.
Despite the extensive amount of research in herpetological systematics during the middle of the 20th century and the rapid expansion of “cladistic” systematics after the English publication of Hennig (1966), few studies examined any squamate groups by explicit quantitative analyses of morphological character matrices in the decades after Camp (1923). As noted above, only a few researchers investigated small groups such as “primitive” snakes (Rieppel, 1979) and geckos (Kluge, 1987). In contrast, no workers attempted an explicit analysis of Squamata as a whole for 65 years, until the first comprehensive phylogeny of the group was estimated by Estes et al. (1988), marking the modern period of squamate systematics a mere 33 years ago.
Early Morphological Phylogenies
The first attempt at a comprehensive phylogenetic analysis of broadscale squamate relationships was conducted by Estes et al. (1988). The analysis included 148 morphological characters and 19 extant squamate families as tips (= operational taxonomic units) analyzed by maximum parsimony—an approach followed by most subsequent morphological analyses. Squamates were diagnosed by several characters, including fused premaxillae, fused parietals, opisthotics fused to exoccipitals, loss of palatine teeth, procoelous vertebrae, loss of posterior trunk (thoracolumbar) intercentra, and presence of an anterior coracoid emargination. (See Fig. 2 for a brief sample of morphological characters that have been historically used to diagnose squamates and several squamate clades.) Within squamates, the results of Estes et al. (1988) recovered a broad-level structure very similar to the highly influential precladistic classification of Camp (1923) (Fig. 3A). Specifically, the earliest branching event in the squamate tree would be that between Iguania (diagnosed in part by fleshy, lobular tongues; frontals strongly constricted between orbits; postfrontals reduced or lost; and fingerlike angular process on retroarticular process) and Scleroglossa (i.e., all other squamates, diagnosed in part by scaly, forked tongues; prey prehension exclusively mediated by jaws rather than tongue; postfrontals forked medially; enlarged cephalic scales; loss of the dorsal process of squamosal; alar process of prootic elongated; and clavicles strongly angulated, curving anteriorly; e.g., Fig. 2, chars. 1, 4, 5, and 11).
Within scleroglossans, Estes et al. (1988) estimated a major dichotomy between gekkotans (diagnosed in part by the loss of the lacrimals and postorbitals; reduction of the jugals and postfrontals; and persistent notochordal canal in adults, with amphicoelic vertebrae in some families) and autarchoglossans (diagnosed by only three synapomorphies: loss of jugal-squamosal contact on supratemporal arch; presence of dermal bone rugosities; and presence of an m. rectus abdominis lateralis muscle). Autarchoglossans were subdivided into two main lineages: Scincomorpha and Anguimorpha. Scincomorphs included xantusiids, lacertids, teiioids (teiids and gymnophthalmids), cordyloids (cordylids and gerrhosaurids), and scincids, whereas anguimorphs comprised anguids, varanids, helodermatids, and xenosaurids (Fig. 3A).
In the study by Estes et al. (1988), all sampled elongate and limb-reduced or legless taxa except pygopodids—snakes, amphisbaenians, and dibamids—were “clumped” together within anguimorphs. However, this was considered a problematic result because the group was primarily supported by characters associated with limb reduction, which were thought to be convergent. Therefore, the final and preferred phylogenetic tree reported by the authors did not indicate the placement of those three limb-reduced clades. Pygopodids, on the other hand, were recovered as the sister clade to gekkotans instead of clustering with other limb-reduced clades, most likely driven by soft tissue characters and gecko-specific cranial characters that were able to overcome the signal from limbless characters.
A series of subsequent modifications and additions to this initial morphological data set followed in the years afterwards, including the first studies to add fossil taxa to the original all-extant data set (e.g., Wu et al., 1996; Evans and Barbadillo, 1998; Evans and Chure, 1998; Caldwell, 1999). All of those studies except Caldwell (1999) found a very similar basic configuration related to the works of Camp (1923) and Estes et al. (1988). Importantly, the latter studies provided the first phylogeny-based systematic classification of important extinct lineages, such as the placement of Borioteiioidea (= Polyglyphanodontidae) as a sister clade to either Teiioidea or Lacertoidea (Teiioidea+Lacertidae), and paramacellodids as either a sister taxon to or a member of Scincoidea (Scincidae+Cordyloidea) (Evans and Barbadillo, 1998; Evans and Chure, 1998). However, the relationships of other fossil species varied considerably across some of those analyses, such as the placement of the oldest articulated fossil lizards known at the time from the Late Jurassic of Germany, including Eichstaettisaurus, Bavarisaurus, and Ardeosaurus (Evans and Barbadillo, 1998; Evans and Chure, 1998).
The largest reassessment of the original matrix of Estes et al. (1988) was provided by Caldwell (1999), who revised the characters from the original study on the basis of their interdependence and also removed soft tissue characters (which could not be scored for fossils), resulting in the exclusion or merger of nearly one-third of the original characters. New fossil taxa were also included, such as the oldest known fossil snake at the time (Dinilysia patagonica) and Mosasauria (dolichosaurids, aigialosaurids, and mosasaurids). The results still estimated iguanians as the earliest diverging squamate crown clade along with a monophyletic Anguimorpha and Scincomorpha. However, gekkotans were found to be the sister lineage of scincomorphs instead of autarchoglossans. Importantly, mosasaurs were found to be the sister lineage of snakes, with Mosasauria+Serpentes representing an early-branching lineage in squamate evolution.
After Estes et al. (1988), the first new morphological data set aimed at a comprehensive squamate phylogeny was provided by Lee (1998). After a considerable reevaluation of previous data sets and inclusion of a large number of new observations, this matrix totaled 230 characters scored for 22 taxa. This data set and its subsequent modifications and expansions (Lee and Caldwell, 2000; Lee, 2005a,b) provided (at that time) the most substantial shift from the classical framework previously established by Camp (1923) and Estes et al. (1988) (Fig. 3B).
This phylogeny was the first to estimate scincomorphs to be paraphyletic: Scincids and cordyloids were more closely related to anguimorphs (clade Diploglossa) instead of forming a clade with lacertoids (lacertids+teiioids). This novel clade was diagnosed in part by the retroarticular process twisted distally, broad distally and with a smooth dorsal surface; osteoderms present both dorsally and ventrally on the body and over the skull table; and caudal transverse processes converging distally. Furthermore, Lee (1998) recovered Gekkota and Xantusiidae as sister lineages and within the same clade as Dibamidae and Amphisbaenia—clade Nyctisaura, a clade name resurrected from Cope (1900)—diagnosed in part by the absence of a postorbital, an occipital condyle weakly or strongly bipartite; medial closure of the Meckelian canal by the dentary; and absence of the angular as a discrete element. Finally, this study also found Mosasauria to be the sister taxon to early-diverging snakes within the clade Pythonomorpha (diagnosed by several synapomorphies, including the quadrate being suspended mostly by the supratemporal; strongly developed parietal downgrowths; stapedial footplate surrounded by the prootic and opisthotic [but see for the homology of this character Rieppel and Kearney, 2002]; supraoccipital sutured to parietal; mobile contact between angular and splenial [contributing to an intramandibular joint]; recumbent replacement teeth; reduction of median premaxillary teeth; and presence of zygosphenes and zygantra on the vertebral neural arches; e.g., Fig. 2, char. 9).
The resurrection of Pythonomorpha, also a name first proposed by Cope (1869), for grouping mosasaurs and snakes was also recovered by other authors around the same time (Scanlon, 1996; Caldwell, 1999), including slightly earlier versions of the Lee (1998) data set (Lee, 1997a,b). This is perhaps the most controversial result among morphological phylogenies because it led to a series of subsequent debates on the marine, terrestrial, or fossorial origin of snakes and their placement within the squamate Tree of Life (e.g., Rieppel and Zaher, 2000; Tchernov et al., 2000; Conrad, 2008; Gauthier et al., 2012; see below under “The Origin of Snakes”). Those controversial (although not necessarily unexpected) results were subsequently maintained by analyzing soft tissue data alone (Lee, 2000), and some recent morphological studies (see below). It is perhaps surprising that the placement of snakes and mosasaurs as sister taxa within squamates (both groups never having been tested together in a phylogenetic analysis before the late 1990s in the aforementioned studies) sparked more controversy than the proposed paraphyly of scincomorphs (a clade that had been established for nearly 75 years) and the subsequent placement of cordylids and scincids with anguimorphs, in direct conflict with the traditional classifications of Camp (1923) and Estes et al. (1988). This situation underscores the effect that forces such as an organism's charisma or even its penetration into popular culture and religion can have on the broader effect of scientific discovery and research agendas, regardless of how much they actually contradict previously established hypotheses.
Large-scale Morphological Data Sets
Two decades after the first quantitative phylogenetic analysis of squamate relationships (Estes et al., 1988), new morphological data sets encompassing an increasingly larger number of characters and taxa began to emerge. These data reflected a growing body of evidence during the late 1990s and early 2000s regarding the advantages of increased taxon and character sampling in phylogenetic inference (Graybeal, 1998; Wiens, 2004), specifically regarding the benefits of including fossil taxa (Wiens, 2003, 2006). These matrices also represented further attempts to “test” the morphological signal in squamates (i.e., the strength of the morphological data for driving broad-level squamate phylogenetic relationships), given the contrasting results obtained by the first emerging molecular phylogenies, which contested all previous morphology-based hypotheses of squamate evolution (Townsend et al., 2004; Vidal and Hedges, 2005; see below).
The first large-scale morphological data set for squamates was that of Conrad (2008), including 223 taxa and 363 primarily osteological characters. This data set and its subsequent expansions and modifications, including combined evidence analyses (e.g., Wiens et al., 2010; Simões et al., 2015a; Conrad, 2017; Pyron, 2017) still represent the most extensive taxonomic sampling with morphological data available today. The overall structure of the tree in all studies analyzing the Conrad (2008) matrix alone is similar to that of Estes et al. (1988), including the early divergence of iguanians and Scleroglossa, with a monophyletic Autarchoglossa, Scincomorpha, and Anguimorpha (Fig. 3C). However, the much broader taxonomic sampling of extant and fossil taxa had the dual effects of breaking long branches and introducing more terminals with observed character state reversals. Thus, considerably fewer characters were found to diagnose each major clade unambiguously. For instance, total group Scleroglossa (or Scincogekkonomorpha under Conrad's terminology) was diagnosed by only five synapomorphies (compared with 27 by Estes et al., 1988). Similar to previous morphology-based phylogenies focused on iguanian relationships (Frost and Etheridge, 1989), Conrad (2008) obtained very low resolution and stability among the different clades of iguanians, a limiting aspect also observed later in Gauthier et al. (2012). Snakes, amphisbaenians, dibamids, limb-reduced skinks (Feylinidae and Acontidae) were clustered in a limb-reduced clade: Scincophidia (except pygopodids). Within snakes, scolecophidians were the earliest evolving clade of snakes, with subsequent divergence of extinct snake taxa (e.g., Dinilysia, Pachyrhachis, and Wonambi), and finally caenophidians+henophidians. The overall structure of Conrad's (2008) parsimony tree was later sustained under both maximum parsimony and Bayesian inference (Wiens et al., 2010).
The much greater sampling of fossil taxa by Conrad (2008) relative to all other studies also provided the first assessment of some fossil families in a comprehensive estimate. For instance, necrosaurids were recovered as early-diverging varanoids, whereas several fossil species from the Late Cretaceous of Mongolia (e.g., Gobiderma, Estesia, and Paraderma) were recovered as closely related to the Gila monster (helodermatids) in a new clade named Monstersauria. Importantly, Mosasauria was recovered as deeply nested within varanoids and only distantly related to snakes, in sharp contrast to most of the work done during the previous decade that established a sister group relationship between snakes and mosasaurs (Caldwell and Lee, 1997; Lee, 1998; Caldwell, 1999; Lee and Caldwell, 2000; Lee, 2005b; but see Rieppel and Zaher, 2000).
The next step towards the expansion of character and taxon sampling for squamate phylogeny was provided by Gauthier et al. (2012), with 192 taxa and 610 morphological characters, including highly detailed character descriptions and illustrations. Similar to Conrad (2008), despite a considerable increase in the number of taxa and characters, the overall tree structure of extant squamate families under maximum parsimony reinforced the framework first introduced by Camp (1923) and Estes et al. (1988) (but different from Lee  and its subsequent expansions) and contested results obtained from molecular data sets (Fig. 3D). However, Bayesian inference revealed a significant lack of resolution by finding a large polytomy among the major iguanian and scleroglossan clades. As with nearly all previous morphology-based squamate phylogenies, most limb-reduced taxa (except pygopodids) formed a clade that was inclusive of snakes, amphisbaenians, and dibamids, plus Anniella (a limb-reduced anguid) and the fossil taxon Sineoamphisbaena. This assemblage was called a “fossorial group,” although the vast majority of snake species are not fossorial (see Wallach et al., 2014).
Gauthier et al. (2012) included a much larger sampling of snake taxa than any previous studies and found contrasting results concerning early snake phylogeny. Whereas Cretaceous fossil snakes such as Dinilysia and Najash were estimated as the earliest evolving snakes under maximum parsimony, scolecophidians were found as the earliest evolving snakes under Bayesian inference. Also in contrast to previous studies, Gauthier et al. (2012) found new and controversial placements for some fossil lizard lineages, most notably borioteiioids and mosasaurians. The two groups were estimated as non-scleroglossan squamates, suggesting that they are the earliest diverging squamate lineages after the initial divergence of Iguania, and thus successive sister lineages to Scleroglossa. The subsequent expansion of their data set by Reeder et al. (2015) with the addition of 81 soft tissue characters (totaling 691 morphological characters), yielded similar results when analyzed alone under all tested optimality criteria (maximum parsimony, maximum likelihood, and Bayesian inference).
After nearly 100 years of attempts to establish relationships among the major squamate lineages, dozens of different analyses and data sets have yielded a few consistent but occasionally opposing paradigms with the use of morphological data. First, most studies agree that many squamate groups long recognized by prephylogenetic workers (going back to Linnaeus or earlier) are monophyletic, including snakes, geckos, skinks, and several fossil lineages, including Mosasauridae, Borioteiioidea, Paramacellodidae, Madtsoiidae, and Glyptosaurinae. A primary challenge, however, has been identifying how those different groups are related to each other. Another common point of agreement is in the placement of Iguania as the earliest diverging group of crown squamates. There is also widespread agreement on a close relationship between anguids, helodermatids, xenosaurids, and varanoids, comprising Anguimorpha. Most studies also concur on a close relationship between paramacellodids and scincoids (either within scincoids or the sister clade), as well as the placement of Mosasauria within anguimorphans, regardless of whether snakes were found as the sister group to Mosasauria or elsewhere.
The most considerable sources of disagreement between morphological hypotheses stem from two major groups of competing data sets and analyses: one with the basic tree structure first introduced by Camp (1923) and maintained by Estes et al. (1988), Conrad (2008), and Gauthier et al. (2012) and the one introduced by Lee (1998) and further promoted by Lee and Caldwell (2000), among others. Whereas the former is in agreement on the monophyly of Scincomorpha and the placement of snakes in the “fossorial” clade, including dibamids and amphisbaenians, the latter suggests the placement of scincoids in a clade with Anguimorpha and the placement of snakes as the sister lineage to Mosasauria. Another and less-recognized source of discordance between morphological data sets is the conflict regarding the placement of geckos. Whereas most morphological phylogenies recovered gekkotans as the sister clade to autarchoglossans, diverging early in squamate evolutionary history subsequent to iguanians (e.g., Estes et al., 1988; Conrad, 2008; Gauthier et al., 2012), other morphological data sets recovered gekkotans in a clade with Mosasauria and snakes (Caldwell, 1999; Simões et al., 2017b), close to borioteiioids (Lee, 2005a), or possibly forming a clade with borioteiioids or amphisbaenians, dibamids, and xantusiids, or a combination of these groups (Lee, 2000, 2005a,b).
Similarly, the controversy over the placement of snakes as either the sister lineage of Mosasauria or in a clade with other limb-reduced taxa sparked a wave of debate over the origin of snakes, including marine (Lee, 1997a, 1998, 2005a,b), terrestrial (Vidal and Hedges, 2004), or fossorial (e.g., Rieppel and Zaher, 2000; Gauthier et al., 2012; Hsiang et al., 2015) hypotheses for the ecological driver of early snake evolution. The “snake origins” debate further boosted another discussion on the reliability of morphological characters in the face of widespread apparent homoplasy, particularly those associated with limb reduction in squamates (Lee, 1998, 2005a,b). Characters highly correlated with limb reduction, miniaturization of the skull, and fossorial adaptations seem invariably to “attract” limb-reduced taxa to cluster in morphological data sets, unless each limb-reduced taxon is included and analyzed separately from other limb-reduced clades or character weighting is applied (Lee, 1998, 2005b). Yet, few if any researchers actually “believe” that these taxa form a natural group, particularly when it includes species such as legless skinks, which are almost indisputably highly derived and deeply nested in Scincidae, for instance (see Wiens and Lambert  for more discussion of this point).
Subsequent studies eventually called into question the entire conceptual and methodological approach taken during the construction of nearly all previous squamate morphological data sets. Upon a reassessment of more than 1,000 morphological characters used by earlier authors, Simões et al. (2017b) found that ∼35%–45% of all characters used in previous analyses violated basic principles of morphological character construction, including not accounting for logical and biological dependencies among characters (Wilkinson, 1995; Strong and Lipscomb, 1999), not keeping multiple exclusivity among states of a single character (Sereno, 2007; Brazeau, 2011), violating principles of character coding that avoid biases introduced by inapplicable or missing data (Strong and Lipscomb, 1999; Brazeau, 2011), and not treating continuous variables as such by introducing arbitrarily delimited character discretization (Rae, 1998; Goloboff et al., 2006). Upon reevaluating the largest available squamate morphological data sets at the time (Conrad, 2008; Gauthier et al., 2012) and recoding or removing problematic characters, Simões et al. (2017b) found results in substantial contrast to the original studies. A reanalysis of the matrix of Conrad (2008) resulted in a tree with strong similarities to that of Lee (1998) and subsequent expansions thereof, by estimating scincids and cordyloids as more closely related to anguimorphs (Diploglossa). In a similar reanalysis of Gauthier et al. (2012), Mosasauria moved from an early-diverging position in the squamate tree (an outgroup to Scleroglossa) to form a clade along with the limb-reduced squamates, including snakes.
Following stricter criteria for morphological character construction and an expanded taxonomic sampling of rhynchocephalians, early lepidosaurs, and other diapsid lineages, Simões et al. (2018) built a new morphological data matrix inclusive of the major stem and crown extinct and extant lineages of squamates. The result was the first ever agreement of the morphological signal with the molecular signal (see below) for several key nodes in the squamate Tree of Life (i.e., gekkotans at the base of the tree and iguanians within a monophyletic Toxicofera; Fig. 4). Numerous fossil lineages were included, with borioteiioids recovered in their classical placement as closely related to teiids (specifically, the sister group to Teiioidea) and paramacellodids as the sister group to skinks. Mosasaurians were found to be the sister taxon to snakes within Toxicofera, and several previously overlooked fossil lepidosaurs were found to represent some of the oldest known squamate fossils, expanding the fossil record of the group to the Middle Triassic. This result offers great promise for unifying the phylogenetic signal found in different ontological partitions of squamate biology (see Losos et al., 2012), perhaps the final remaining frontier for understanding the squamate Tree of Life.
History of the Molecular Hypothesis
Comprehensive molecular analyses of squamate phylogeny were relatively late in arrival, only beginning in earnest several years into the 21st century. In the 1990s, however, a few studies began to suggest that molecular data may exhibit substantially different phylogenetic signals than traditional morphological hypotheses (see above). In a study of mitochondrial transfer RNAs (tRNAs), Kumazawa and Nishida (1995) sampled a gecko, a skink, and an iguanian, and found that either the skink or the gecko represented the earliest divergence in squamates in various analyses, but never the iguanian. In a data set sampling the mitochondrial gene ND4 and associated tRNAs for several snakes, an anguimorph, a gecko, two iguanians, and a skink, Forstner et al. (1995) found that the skink was the sister lineage to the remaining taxa. A reanalysis of this data set by Macey and Verma (1997: 277) confirmed support for this placement, noting that “no topology is consistent with the squamate phylogenetic hypotheses derived from morphological data.” Finally, Gorr et al. (1998: 481) noted that hemoglobin sequences did not place Serpentes closer to Anguimorpha than Iguania as hypothesized by Estes et al., (1988) and others, stating “neither single chain . . . nor any of the a–β tandem alignment trees . . . are consistent with this hypothesis.”
Starting with Saint et al. (1998), a series of studies (e.g., Harris et al., 1999; Harris, 2003) began to assemble increasingly well-sampled data sets comprising fragments of the single-copy, intron-free nuclear protein–coding gene CMOS for a broad array of squamate species. In essentially all analyses, Iguania, Anguimorpha, and Serpentes formed a clade, and Teiioidea was the earliest diverging squamate lineage. As Harris (2003: 540) noted, “such a relationship is completely anomalous to estimations of relationships based on morphology or the fossil record.” Vidal and Hedges (2004) subsequently combined parts of these data sets with sequence fragments from a second protein-coding nuclear gene (RAG1), sampling nearly all squamate families recognized at the time. Node support varied, but their analyses estimated a topology with Dibamidae as the earliest diverging squamate lineage, followed by geckos; skinks and relatives; teiids, lacertids, and amphisbaenians; and a clade comprising Anguimorpha, Iguania, and Serpentes (Fig. 5). Despite the extraordinary incongruence with previous morphological hypotheses, Vidal and Hedges (2004) were primarily concerned with snake origins and thus did not offer much additional commentary on the implications of their results for understanding higher level squamate relationships.
The landmark study of Townsend et al. (2004) was the first multilocus data set (CMOS+RAG1) to address explicitly higher level squamate relationships with relatively comprehensive taxon sampling. They corroborated Vidal and Hedges (2004) in cementing the “molecular” hypothesis, in which dibamids or geckos are the earliest diverging lineages, followed by skinks and relatives; teiids, lacertids, and amphisbaenians; and, again, a clade comprising Anguimorpha, Iguania, and Serpentes (Fig. 5). They noted “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” (Townsend et al., 2004: 735) and suggested that the morphological topology was likely due to homoplasy. These results were quickly corroborated by Vidal and Hedges (2005), sampling nine nuclear protein–coding genes (albeit for fewer taxa), who introduced a tree-based nomenclature for these groups. In particular, Scinciformata comprises Scincidae, Cordylidae, and Xantusiidae, whereas Laterata comprises teioids, lacertids, and amphisbaenians. Finally, the name Toxicofera was applied to the clade of Anguimorpha, Iguania, and Serpentes. Like Townsend et al., (2004) Vidal and Hedges (2005) also suggested that phenotypic homoplasy was the most likely explanation for the unorthodox results compared to the morphological hypothesis.
A flood of subsequent research yielded highly similar results across a variety of data sets, including full mitochondrial genomes (Kumazawa, 2007), dozens of protein-coding nuclear genes (Wiens et al., 2012), thousands of species sampled for a smaller set of mitochondrial and nuclear genes (Pyron et al., 2013), thousands of ultra-conserved element (UCE) loci (Streicher and Wiens, 2017), thousands of orthologous loci generated from RNA transcriptomes (Irisarri et al., 2017), hundreds of anchored hybrid-enrichment (AHE) loci (Burbrink et al., 2020), and a combination of UCEs, AHE, and nuclear protein–coding loci (Singhal et al., 2021). This drastic overhaul of traditional, primarily morphology-based classification mirrors similar situations that occurred in other groups, such as lissamphibians (see Frost et al., 2006). Thus, the “molecular” hypothesis of squamate relationships is now broadly accepted by most systematists. Questions remain, however, regarding the incongruence with existing morphological data sets and the many poorly supported branches observed across the myriad molecular studies referenced above.
Molecular Data Sets and Uncertainty
As noted, consistency across molecular data sets is remarkable in terms of support for the “molecular” hypothesis of a nested Iguania (Fig. 5). The signal for the molecular hypothesis is not marker specific but is found across both the nuclear and mitochondrial genomes in both coding and noncoding regions. However, support is not universal across these loci. For example, seemingly artifactual results (paraphyletic Iguania and early-branching Serpentes) from some whole-mitochondria alignments are apparently the result of long-branch attraction and an ancient episode of convergent adaptive evolution (Castoe et al., 2009). Similarly, the BovB LINE family of transposable elements has experienced extraordinary levels of ectoparasite-mediated horizontal transfer between mammals and squamates (see Pasquesi et al., 2018, and references therein), and its phylogeny is essentially uninterpretable for understanding squamate relationships.
Almost as remarkable as the consistency across studies for the molecular hypothesis is their consistency in reflecting the same regions of uncertainty in the tree, with a set of poorly supported branches being common to most analyses (see Burbrink et al., 2020; Singhal et al., 2021, and references therein). Along the backbone of the tree, two major nodes merit discussion. The first is the most recent common ancestor (MRCA) of Squamata, where the earliest branch is either Gekkota, Dibamia, or both, forming a clade that is the sister lineage of remaining squamates. The second is the resolution of Toxicofera: whether Serpentes is the sister lineage to Iguania, Anguimorpha, or both.
Within higher level clades, four major groups stand out. The first are pleurodont iguanians, where relationships among the numerous families are famously intractable (see Townsend et al., 2011). The second is monophyly of Amphisbaenia: whether Lacertidae is more closely related than Rhineuridae to the remainder of amphisbaenians (Wiens et al., 2010; Streicher and Wiens, 2017). The third is monophyly of Scolecophidia (blindsnakes), where Anomalepididae is sometimes the sister lineage of all other snakes exclusive of Typhlopidae+Leptotyphlopidae, or of all other snakes (e.g., Burbrink et al., 2020). Fourth, within Caenophidia, relationships among the numerous families of booid and colubroid snakes often receive strong support for vastly different topologies (see Pyron et al., 2014; Reynolds et al., 2014; Streicher et al., 2016, and references therein) mirroring the situation in pleurodont iguanians.
Most of these instances (except Amphisbaenia) appear to stem from the well-known problem of rapid evolutionary radiations and long descendant branches. These scenarios can be difficult if not impossible to resolve because of a lack of informative sites accumulating during brief periods of divergence, as well as frequent homoplasy arising along deep stem lineages (Rokas and Carroll, 2006). Phylogenomic data sets alone may not be enough to overcome these various sources of nonphylogenetic signal and other artifacts, which may stem from nonorthology, model error for gene trees or species trees, taxon sampling, and so on (Philippe et al., 2011). In contrast, non-monophyly of Amphisbaenia is apparently due to long-branch attraction and is seemingly “corrected” by increased taxon sampling in molecular data sets (see below and Pyron et al., 2013).
Advancements and Controversies
Recent attempts to understand the distribution and variation of phylogenetic signal in phylogenomic data sets for squamates have corroborated these hypothesized explanations for topological discordance and low support but have unfortunately offered limited resolution, although with some promising potential avenues for future investigation. McMahan et al. (2015) suggested that since the branch leading to the sole outgroup Sphenodon was so long, the signal for the root of Squamata was most likely obliterated. Thus, the rooting of the ingroup may essentially be random. They then claimed to show that rerooting the molecular tree with Iguania is supported by more molecular synapomorphies than when the tree is rooted with Teiioidea or Dibamidae.
Harrington et al. (2016) showed that this procedure erroneously conflated branch support with branch length. The number of molecular synapomorphies counted along a branch is an estimate of the length of that branch in substitutions, not a measure of character support. Thus, McMahan et al. (2015) were in effect arguing that the longest ingroup branch of any clade should be the root. Furthermore, their counts were obscured by collapsing ingroup clades to a single basal dichotomy and by considering only two possible roots. When the full tree is analyzed under their procedure, both Teiioidea and Serpentes represent “longer” root branches than Iguania. Additionally, no known gene tree places Iguania at the root, and none of the dozens of published molecular phylogenies have obtained this result. If the Sphenodon rooting was “random,” one would presumably expect occasional molecular support for a basal Iguania or a non-Gekkota/Dibamia root in at least some trees.
Another attempt to explain the differences between (most of) the morphological and molecular phylogenetic data sets (Koch and Gauthier, 2018) raised similar objections, claiming to show that the nested placement of Iguania was due to systematic bias and noise in molecular data—in essence, long-branch attraction of Iguania and Serpentes. They purported to demonstrate this by removing “fast-evolving” sites, which collapsed support for Toxicofera and eventually all backbone nodes as more characters were removed. Subsequently, the authors showed that individual genes did not decisively support Toxicofera with significance tests of topology. They suggested that because 85% of 46 genes from Reeder et al. (2015) included more than one topological alternative in the confidence set, Toxicofera was not supported decisively by the molecular data. Finally, they offered an analysis of phylogenetic informativeness (PI) of the 46 genes from Reeder et al. (2015). These showed younger peaks in PI than the ages of the root branches of interest, concluding that these data sets thus had inadequate power to resolve early divergences.
Owing to the recency of this criticism, to our knowledge it has so far not been critically evaluated in the literature. However, some immediate limitations are apparent in their conclusions, including:
1) The authors did not provide any direct evidence that “fast-evolving” sites were artifactually biased by any specific mechanism, beyond showing that snakes and iguanians have higher rates of molecular evolution and a slight Adenine-Thymine (AT) bias. Thus, they have shown only that removing the most potentially informative sites reduces phylogenetic resolution, which is a truism that does not bear on the placement of Iguania. Importantly, it has been recognized that the branch leading to snakes also exhibits exceptionally high rates of morphological evolution (Watanabe et al., 2019; Simões et al., 2020b) which are in fact proportionately faster than molecular rates of evolution (Simões et al., 2020b).
2) The lack of decisiveness for individual genes is a theoretical expectation for short internodes because of incomplete lineage sorting that is typically resolved by the use of species tree methods (Edwards, 2009) and thus does not reflect on the inclusive phylogenetic signal of the constituent genes. In contrast, well-sampled data sets of stochastically discordant gene trees are expected to converge on the correct species tree with high confidence, even when a low proportion of individual genes matches the true topology (Edwards et al., 2007; Edwards, 2009).
3) Although the root nodes of Squamata do not exactly match the peak of PI in terms of timing, they are still close to it, with drastically higher PI than later in the tree. Under the authors' logic, it should thus be difficult to estimate recent divergences in Squamata with these data; yet, as noted earlier, morphological and molecular data are highly congruent for many younger nodes and estimated relationships. It is also unclear, under the authors' formulation, how it would ever be possible to estimate more ancient divergences with such data (e.g., Hugall et al., 2007), in that they come before the peak of PI.
4) Although the authors briefly mention the 4,178 UCE data set of Streicher and Wiens (2017), they only tested the signal of the 46 loci of Reeder et al. (2015). Thus, they did not address how the stochastic factors that they claim drive the putatively erroneous estimate of Toxicofera could be so consistently and widely distributed across both the nuclear and mitochondrial genome and thousands of coding and noncoding loci, without even occasional support for a basal Iguania being found among thousands of independent gene trees.
More recently, two studies took a systematic approach to examining locus-level variation in the phylogenetic signal in squamate genome-scale data sets, one using 394 AHE loci (Burbrink et al., 2020) and a combined probe set (Singhal et al., 2021) consisting of 38 nuclear protein–coding loci from Wiens et al. (2012), 5,052 UCEs from Streicher and Wiens (2017), and 372 AHE loci from Singhal et al. (2017). Both were consistent in recovering a reassuring set of explanatory variables for discordance among loci. Generally speaking, longer loci with more informative sites are more congruent in supporting the molecular hypothesis and Toxicofera. Many loci that do not support these nodes instead exhibit incomplete lineage sorting, which is resolved by the use of coalescent species tree methods and subsequently favors Toxicofera. Finally, some discordance is due to short internodes in the species tree, for which few informative sites have accumulated in any locus. Additionally, some gene tree error is driven by noise arising from a nonphylogenetic signal such as loci with few informative sites, model violations such as rate heterogeneity or saturation, and alignment issues such as low sequencing or assembly quality or incorrect ortholog identification (see Phillippe et al. 2011). However, these low-quality loci do not support any alternative topology strongly, and their removal only increases support for the molecular hypothesis and Toxicofera, while no locus supports Scleroglossa or a basal Iguania.
Concomitantly, both Burbrink et al. (2020) and Singhal et al. (2021) reported continuing uncertainty regarding numerous nodes. Chief among these is the placement of Dibamia, which represents a very short internode subtending a very long branch; conditions that make resolution extremely difficult (Philippe et al., 2011) and may continue to be intractable regardless of the number of loci sampled until other types of data are examined (see below). Similarly, the resolution of Toxicofera remains indecisive, although most recent analyses are congruent in estimating an Iguania-Anguimorpha clade. The same is true of scolecophidian monophyly, where many recent data sets are congruent in recovering Anomalepididae as the sister lineage to all other snakes exclusive of Typhlopoidea. Like Dibamia, both Toxicofera and Scolecophidia represent short internodes with long descendent branches that may require alternative data sources to resolve conclusively.
Another contentious point is monophyly of Amphisbaenia. Burbrink et al. (2020) did not sample Rhineuridae (likely the earliest branch in the group), but both Streicher and Wiens (2017) and Singhal et al. (2021) included it and estimated weak support for amphisbaenian paraphyly. This is an apparent case of long-branch attraction that has been putatively solved through increased taxon sampling in other molecular data sets (see below). Taxon sampling may be a factor in many of the examples discussed here, because most genome-scale data sets thus far have only sampled dozens or a few hundred representatives from the ∼11,000 species total diversity of Squamata.
Finally, the problematic radiations of booid snakes and pleurodont iguanians were shown to be in the anomaly zone (Singhal et al., 2021), where short internodes predominantly produce discordant gene trees and cause nearly insurmountable biases for estimating the species tree (Degnan and Rosenberg, 2006). Although sampling a greater number of informative loci should eventually allow resolution in these scenarios (Linkem et al., 2016), that level has clearly not been reached, even with thousands of markers in squamates. As few as three short internodes in rapid succession can produce as many as 45 anomalous gene trees (Rosenberg and Tao, 2008), and Burbrink et al. (2020) and Singhal et al. (2021) find at least four and nine short internodes in Pleurodonta, respectively. Thus, whole-genome data or other character types (e.g., synteny mapping; see below) may be needed to provide any conclusive resolution for these branches. Neither Burbrink et al. (2020) nor Singhal et al. (2021) sampled enough colubroid or elapoid snake families to make robust inferences regarding those historically difficult groups (see Pyron et al., 2014), but similar processes seem likely to be at play, with numerous short branches during rapid Cenozoic radiations in advanced snakes (Pyron et al., 2011; Zaher et al., 2019).
COMBINED EVIDENCE OF SQUAMATE RELATIONSHIPS
Early Searches for Congruence
The phylogeny of Estes et al. (1988) cemented the “morphological” hypothesis in the minds of most herpetological systematists as the best quantitative estimate of squamate relationships. Concomitantly, as molecular data increasingly gained ascendancy in the field of systematics during the late 1990s and early 2000s, such data sets provided spectacular new resolutions of phylogeny and taxonomy for numerous squamate groups. Studies were published on dozens of families and genera, with examples such as Gymnophthalmidae (Pellegrino et al., 2001), rat snakes in the genus Elaphe sensu lato (Utiger et al., 2002), Scincidae (Whiting et al., 2003), Iguania (Sites et al., 1996; Frost et al., 2001), and dwarf boas (Wilcox et al., 2002), representing merely a random sample of papers from this period with .100 citations on Google Scholar as of February 2021. Anticipation in the field for a widely sampled phylogeny of Squamata thus began to increase rapidly, holding the potential promise of merging morphological and molecular data matrices to provide a new, unified view of lepidosaur evolution.
At the same time, molecular trees such as those of Saint et al. (1998) began to accumulate in the literature, and rumors of preliminary results from several labs that would result in publications such as Townsend et al. (2004) and Vidal and Hedges (2004) began to circulate among workers in the field. Thus, many herpetologists gradually became cognizant that the “molecular” tree might end up being significantly different from the “morphological” tree. With this dawning awareness in the background, the United States National Science Foundation (NSF) awarded approximately US$2.4 million in 2004 to an international team of researchers under the “Assembling the Tree of Life” program to investigate squamate origins with molecular and morphological data and estimate a total evidence phylogeny. Titled “The Deep Scaly Project: Resolving Squamate Phylogeny Using Genomic and Morphological Approaches,” this project aimed at a grand synthesis. As the official NSF abstract stated (DEB-0334966):
Many critical questions in squamate evolution remain unresolved, such as identification of the most primitive lineage of squamates, the origin of snakes, and the relationships of venomous snakes to other snake lineages. An international team of eight investigators from diverse institutions (Brigham Young University, Field Museum of Natural History, San Diego State University, State University of New York-Stony Brook, University of Adelaide, University of Texas-Austin, and Yale University) will collaborate to resolve squamate relationships. Anatomical data from living and fossil forms will be combined with DNA sequences from 50 genes for 142 representative squamate species. Anatomical data will be obtained using traditional methods and new high resolution X-ray scanning techniques. DNA data will be generated by incorporating new tools and databases from recent vertebrate genome projects. Computer modeling will be used to determine how data from molecular and fossil studies can best be combined to reconstruct evolutionary trees.
As the numerous publications referenced above appeared between 2004 and 2010, it became increasingly apparent that an illuminating evolutionary consensus of the various data sets would not be forthcoming. After 8 years of attempting to find a congruent explanation linking the molecular and morphological data sets, Daza (2014: 341) notes that: “. . . the two teams diverged and published independently, and each group developed a different phylogeny based on their own independent analyses of a data set that was composed of the largest data partitions for squamate reptiles ever produced (Gauthier et al., 2012; Wiens et al., 2012).”
In many ways, this fundamental divergence is still among the primary questions for researchers in squamate systematics. As Losos et al. (2012) asked: Who speaks with a forked tongue? Which data set is “lying” about the true order of branching events? They note (Losos et al., 2012: 1429): “When two phylogenies are fundamentally discordant, at least one data set must be misleading. . . . We are left with a conundrum. The molecular data imply an astonishing pattern of morphological homoplasy and suggest very limited knowledge of the functional link between structures and lifestyle; if convergence is so pervasive, what faith can we have in the placement of fossil taxa for which no molecular data are available? Conversely, morphology implies a pattern of molecular evolution that has yet to be explained.”
It is from this basis, a fundamentally puzzling and apparently intractable clash of phylogenetic signals from what should be epistemologically equivalent data sets, that all subsequent attempts at integration and consilience have proceeded. These include early attempts published before the main Deep Scaly data sets (Lee, 2005a, 2009), continuations of the Deep Scaly project (e.g., Reeder et al. 2015), comparisons of multiple existing morphological data sets including Conrad (2008) and Gauthier et al. (2012) by Pyron (2017), and de novo construction of morphological matrices in an attempt to understand and eliminate bias in character ontology and coding (Simões et al., 2017b,, 2018). With this recent historical background in mind, we now examine some of these results in more detail.
Results, Successes, and Limitations
To our knowledge, the first study to combine morphological and molecular data as independent partitions to assess broadscale relationships in squamates was that of Lee (2005a). In that study, the molecular data set from Vidal and Hedges (2005) was combined with an updated version of a previous morphological data set (Lee and Caldwell, 2000). Given the relatively small size of the molecular data set (especially in terms of informative sites), the combined evidence tree analyzed under maximum parsimony greatly resembled the framework obtained by analyzing the morphological data set of Lee (1998) and its subsequent expansions alone. Specifically, the earliest split in the squamate tree represented by that between iguanians and scleroglossans, a monophyletic Nyctisaura (Xantusiidae+Gekkota+Dibamidae+Amphisbaenia), and Diploglossa (Scincoidea+Anguimorpha), in which early mosasaurians represent the sister taxa to Serpentes. Some important novelties however included the placement of borioteiioids and lacertoids as the earliest diverging branches within Scleroglossa, before the dichotomy between Nyctisaura and Diploglossa.
Despite the initial dominance of the morphological signal over the molecular signal, this scenario soon reversed with the incorporation of larger molecular data sets into combined evidence studies of squamate phylogeny. After combining nine nuclear and one mitochondrial loci (totaling 8,086 bp) from both Townsend et al. (2004) and Vidal and Hedges (2005), Lee (2009) found a result that was largely congruent with previous molecular phylogenies of squamates. Ever-growing molecular data sets increasingly tended to override most signal available from morphological partitions. The molecular component of combined evidence data sets increased to 22 nuclear loci (totaling 15,794 bp) analyzed with the 363 morphological characters from Conrad (2008) by Wiens et al., (2010), and then to 46 loci (45 nuclear and 1 mitochondrial, totaling 35,673 bp) with 691 morphological characters (Reeder et al., 2015). These studies consistently found the estimated trees to reflect the overall structure of the molecular signal, albeit with widely variable placements of the fossil species and extinct lineages.
Similar effects were observed in subsequent studies implementing smaller molecular data sets, including Pyron (2017) with six nuclear loci (totaling 8,710 bp) and a morphological partition of either 363 or 610 characters from Conrad (2008) and Gauthier et al. (2012), and Simões et al. (2018) with 16 loci (13 nuclear and 3 mitochondrial, totaling 11,532 bp) and 347 newly assembled morphological characters (Fig. 6). However, we note that the number of changes from the morphological to the combined evidence tree in Simões et al. (2018) is somewhat reduced because the morphological tree from the latter already recovers relationships more similar to the molecular tree (Fig. 4). Furthermore, Wiens et al. (2010) and Pyron (2017) noted cases in which the addition of morphological data altered the molecular hypotheses, as well, indicating that even relatively small morphological partitions can still exert a decisive influence on larger molecular matrices.
Generally, the overwhelming signal from molecular data in nearly all combined evidence analyses has an indirect effect even on the placement of extinct lineages. The placement of fossil taxa can only be directly influenced by morphological data. However, molecular data influence the relationships between extant lineages for which morphological data are also available. Thus, molecular data can thereby indirectly influence the placement of fossil taxa by altering the likelihood of various arrangements of branches relative to the morphological matrix. The reverse is also true in theory, although such effects are typically observed to a far lesser extent in most empirical analyses (see Wiens et al.. 2010; Pyron, 2017).
For instance, Huehuecuetzpalli mixtecus (Early Cretaceous, Mexico) is commonly estimated as one of the earliest diverging branches of the squamate tree, as an outgroup to iguanians in morphological phylogenies (Conrad, 2008; Gauthier et al., 2012; but see Simões et al., 2018). This results from Huehuecuetzpalli sharing some morphological characters with iguanians, such as a dorsal process of the squamosal and a reduced postfrontal located anteriorly to the postorbital (Fig. 2, chars. 1, 4). However, with the incorporation of molecular data, the combined evidence trees place iguanians within Toxicofera, thus “dragging” Huehuecuetzpalli with them inside Toxicofera and as the sister taxon to iguanians owing to their shared morphological attributes (Wiens et al., 2010; Pyron, 2017; but see Reeder et al., 2015). In a similar manner, the morphological signal linking Mosasauria to snakes and anguimorphs drives mosasaurs to be nested within Toxicofera either as the sister taxon to snakes (Lee, 2009; Reeder et al., 2015; Pyron, 2017; Simões et al., 2018) or anguimorphs (Wiens et al., 2010; Simões et al., 2020b). Finally, extinct subclades of crown families (e.g., glyptosaurines, a subclade of anguids) are found to be nested within anguids (inside Anguimorpha) and, therefore, inside Toxicofera (Wiens et al., 2010; Reeder et al., 2015; Pyron, 2017). We predict that such extinct groups within major crown clades (e.g., glyptosaurines) are highly unlikely to have their phylogenetic placement determined independently by the morphological signal in the presence of sufficiently large (e.g., ca. .8–10 kb) molecular sampling of extant species—except in the case of a rampant homoplasy-driven morphological signal, such as the clustering of limb-reduced taxa.
The major exception to the observed pattern described above is in the case of fossil taxa that do not have strong support from the morphological data to be the sister taxon to (or a subclade of) an extant clade. For instance, the controversial fossil lizard Sineoamphisbaena hexatabularis (Late Cretaceous, Mongolia) has been previously proposed to be the sister lineage of amphisbaenians (see above; Wu et al., 1996; Gauthier et al., 2012), but other studies have placed it as a sister lineage to or a member of Borioteiioidea (Kearney, 2003; Lee, 2005b; Conrad, 2008) or estimated borioteiioids and Sineoamphisbaena as the successive sister lineages to amphisbaenians (Lee, 2009). The lack of consensus on the phylogenetic affinities of Sineoamphisbaena is further evidenced by the fact that all proposed systematic placements for this taxon are usually met with extremely low topological support (e.g., Bremer index of 1 in maximum parsimony analyses). As such, the inclusion of molecular data and the nearly inevitable rearrangement of the relationships among extant families will likely affect and shift those controversial fossil taxa to potentially novel phylogenetic positions compared with the morphological tree (e.g., Lee, 2009; Wiens et al., 2010; Reeder et al., 2015; Pyron, 2017). However, as argued above, such rearrangement is most likely a result of their poor phylogenetic signal rather than the actual capabilities of the morphological data to drive their phylogenetic placement independently of the molecular data. Such paradoxes can only be solved by improving character construction and increasing character sampling in the situations where this is possible, given fragmentary remains.
Other fossil species with similarly non-consensual and poorly supported placement among different analyses include: Ardeosaurus (Late Jurassic, Germany), the Early Cretaceous Calanguban, Tijubina, and Olindalacerta (Brazil); Tepexisaurus (Mexico); Scandensia, Pedrerasaurus, and Jucaraseps (Spain); Dalinghosaurus (China); and Slavoia, Globaura, and Eoxanta (Late Cretaceous, Mongolia) (Conrad, 2008; Bolet and Evans, 2012; Gauthier et al., 2012; Simões et al., 2015a, 2017a). These species are also the ones that tend to reduce topological resolution and support (and most likely, tree accuracy) in combined evidence analyses and total evidence dating. We argue that algorithms for rogue taxon identification should be used to avoid wild-cards (see Pyron 2017), as evidence suggests they have a considerable negative effect for estimating accurate tree topologies under a variety of scenarios, both by increasing the overall amount of missing data and by destroying consensual support for optimal topologies (Luo et al., 2020; Vernygora et al., 2020).
SYNTHESIS IN SQUAMATE EVOLUTION
The Origin of Squamates
Understanding the roots of squamate diversification, including the stem divergence from its sister clade (rhynchocephalians) and the placement of lepidosaurians in the greater context of diapsid reptile evolution, has remained elusive for nearly a century (Broom, 1925; Carroll, 1975, 1977). By far, the greatest limitations have been imposed by a patchy fossil record of squamates from the Triassic and Jurassic, mostly comprising fragmentary jaw elements during the Jurassic (e.g., Waldman and Evans, 1994) and a complete absence of Triassic squamates until very recently. Another major factor was the lack of understanding of broadscale diapsid relationships: Which taxa should act as the outgroup to lepidosaurs, and which clades make up Lepidosauromorpha? Finally, the overwhelming disagreement between most morphological and molecular hypotheses concerning the phylogeny of crown squamates (see above) has hindered a consensual approximation of the relationships between extant lineages, and thus our understanding of recent squamate evolution.
The recent discovery of the oldest known squamate from the Middle Triassic of Italy, followed by the first major agreement between morphological and molecular datasets concerning important parts of the squamate tree, provides the first modern piece of evidence towards elucidating squamate origins (Simões et al., 2018). By including a large number of other diapsid lineages, Simões et al. (2018) determined that many diapsid reptile clades previously proposed as early-diverging lepidosauromorphs, and thus supposedly the closest relatives to lepidosaurs (kuehneosaurids, younginiforms, sauropterygians, among others), actually fall in other parts of the diapsid Tree of Life. This latter conclusion was more thoroughly corroborated by the inclusion of all major clades of diapsids in a single data set by Simões et al. (2018), although other studies with more limited taxon sampling already hinted towards this conclusion over the past 20 years. For instance, kuehneosaurids, gliding reptiles from the Triassic that were once thought to represent the earliest lizards (Robinson, 1962; Estes, 1983), and later to be non-squamate lepidosauromorphs (Benton, 1985; Gauthier et al., 1988a; Motani et al., 1998), have been recovered as non-lepidosauromorphs in the previous large-scale phylogenies inclusive of most diapsid lineages (Müller, 2004; Hill, 2005; but see Chen et al., 2014). The same was detected for younginiforms, terrestrial and semiaquatic Paleozoic reptiles (Müller, 2004; Hill, 2005), and more recently, for sauropterygians, semiaquatic and aquatic reptiles (Hill, 2005; Chen et al., 2014).
The growing current consensus is that the early part of the lepidosauromorph Tree of Life consisted of only a few stem taxa before the divergence between rhynchocephalians and squamates (e.g., Palaeagama, Sophineta, and Vellbergia, spanning the latest Permian [∼252 Ma] to the early Middle Triassic [∼230 Ma]) (Carroll, 1975, 1977; Evans and Borsuk-Białynicka, 2009; Sobral et al., 2020). The identification of Megachirella (Figs. 1A–C) as the oldest stem squamate has now also been further corroborated by expansions and reanalyses of the largest available diapsid-squamate data set (Simões et al., 2018; Garberoglio et al., 2019a; Bittencourt et al., 2020; Simões et al., 2020b; Sobral et al., 2020) by a new data set focused on sphenodontians (Simões et al., 2020a) and by an independently built data set (Scheyer et al., 2020). Given that the oldest known fossil squamate (Megachirella) and rhynchocephalian (cf. Diphydontosaurus) come from the Middle Triassic at ∼242 Ma and ∼230 Ma, respectively, we are now able to narrow down the divergence between lizards and tuataras to the vicinity of the Permian–Triassic extinction. Indeed, this has been recovered by independent data sets implementing phylogenomic (Irisarri et al., 2017; Burbrink et al., 2020) and total evidence dating methods (Simões et al., 2018, 2020). Yet, it remains to be established whether the split between lizards and tuataras happened before or after the Permian–Triassic mass extinction. The latter will depend on even more precise estimates of divergence times and additional fossil material with the results being of great consequence for understanding the early drivers of squamate evolution.
The radiation of crown squamates also finds consensual support on several grounds regarding the placement of gekkotans as the earliest evolving crown clade of squamates, along with dibamids (but see more below in Perspectives). Most notably, the early divergence of gekkotans in the squamate tree has been recovered by essentially all molecular data sets to date and finds agreement in the latest morphological hypotheses (Simões et al., 2018; Simões et al., 2020b). Importantly, this hypothesis also finds support in the fossil record, where the oldest known gekkotans come from the Middle Jurassic (Caldwell et al., 2015), and articulated stem gekkotans are known since the Late Jurassic (Simões et al., 2017a). In contrast, the oldest putative iguanian fossils are represented by fragmentary elements from the Early Cretaceous and more complete specimens only in the Late Cretaceous (Gao and Nessov, 1998; see also review by Simões et al., 2017c). The fossil record may also support an earlier divergence of scincoids (the next major branch in the squamate tree following the consensual trees), because paramacellodids (Figs. 1F–I) also have a fossil record since the Middle Jurassic (Waldman and Evans, 1994) and achieved a near cosmopolitan distribution by the Early Cretaceous (Bittencourt et al., 2020).
Furthermore, it is important to recognize that the placement of gekkotans as part of the earliest radiation of crown squamates has substantial support from a morphological perspective. Gekkotans possess a variable number of features that are common to early evolving diapsid reptiles from the Triassic, including early lepidosaurs and most fossil sphenodontians, but which are absent in all other squamates. Those include the presence of amphicoelic vertebrae and the persistence of a notochordal canal in adults (both in gekkonines and diplodactylines, and also in stem squamates and Eichstaettisaurus), a perforated stapes (in eublepharines and in dibamids), the metotic fenestra undivided externally, although already divided medially (gekkotans and xantusiids), and the presence of paired premaxillae (also in stem squamates and Eichstaettisaurus) (Simões et al., 2018) (e.g., Fig. 2, chars. 7, 8). Those anatomical attributes have historically been suggested as the result of evolutionary reversals within gekkotans (Estes et al., 1988). However, as molecular hypotheses had already suggested (see above and Fig. 5), with further support by recent morphological phylogenies (Simões et al., 2018; Simões et al., 2020b), those attributes actually represent plesiomorphic traits in squamates. On the other hand, purported characters shared by iguanians and sphenodontians are, in fact, the result of convergent evolution, which is further clarified by morphological aspects of stem squamates, especially Megachirella and Huehuecuetzpalli, that are shared with squamate outgroups and at least some gekkotans but are absent in all other squamates, such as the presence of amphicoelic vertebrae and paired premaxillae. Therefore, the inclusion of closely related outgroup taxa to the early divergence of squamates in morphological phylogenies, such as early-evolving sphenodontians and early lepidosaurs, is fundamental to retrieve a more appropriate estimate of character evolution than the historically more common and simplistic approach of including only the extant Sphenodon punctatus, which is separated from squamates by a branch of ∼250 million years.
Phylogenetic Dating of Squamata
Attempts to date the origin of major squamate lineages accompanied early estimates of the “molecular” hypothesis of squamate phylogeny (e.g., Vidal and Hedges, 2005). Those authors estimated a particularly old age (240 Ma; 221–251 Ma range) for the MRCA of crown Squamata in the Middle Triassic. Subsequent authors began to refine these estimates (Mulcahy et al., 2012; Jones et al., 2013) by integrating more complex models for molecular clock estimation and the incorporation of fossil data. These studies, and numerous subsequent estimates, have been remarkably consistent in estimating the MRCA of crown Squamata in the early Jurassic: ∼180–190 Ma. The Timescale of Life database (TimeTree, 2020) records at least 22 studies estimating the age of this node, ranging from 173 Ma (Gamble et al., 2015) to the 240 Ma of Vidal and Hedges (2005), with a median of 190 Ma. However, these estimates are not independent of each other in terms of molecular or fossil data, relying as they do on similar loci and interpretations of existing fossil for use as node age constraints. Thus, these estimates are subject to the myriad concerns and considerations that affect secondary node age calibrations (Parham et al., 2012).
These issues can be partially alleviated by total evidence dating that directly employs fossils as terminal taxa, incorporating morphological as well as molecular substitutions with direct observations of speciation, extinction, and divergence times (Pyron, 2011; Ronquist et al., 2012). However, such methods still suffer from incomplete modeling of relevant processes and adequate associated priors (e.g., Arcila et al., 2015) and are thus still in their infancy. Only two sets of studies have attempted total evidence dating of Squamata to date, those of Pyron (2017) and Simões et al. (2018) [and subsequent expansions of the latter (Garberoglio et al., 2019a; Simões et al., 2020b)], incorporating both molecular and morphological matrices (Table 1).
Divergence times (and 95% highest posterior density when available) estimated from recent phylogenomic and total evidence dating relaxed clock Bayesian inference analyses inclusive of the major lineages of Squamata. Results reported below were obtained from the authors' preferred tree for methodological reasons (see original work for further details). Not applicable (N/A) indicates no distinction between the crown clade and a total clade (absence of stem fossils in the inferred tree topologies).
Using the “standard” Markov K (Mk) models available for morphological evolution under the fossilized birth–death process in MrBayes (see Ronquist et al., 2012), Pyron (2017) recovered an estimate of 190 Ma for the MRCA of crown Squamata. This remarkable congruence with the DNA-only studies mentioned above is notable because he included only a single calibration on the root of Lepidosauria (238–250 Ma) from Jones et al. (2013). All other temporal calibration came from the observed ages of the fossil tips, with no internal node calibrations, and yet still recovered the same estimate of ∼190 Ma as previous studies. Incorporating a more complex model for morphological evolution (MkA), he estimated a slightly younger date of 187 Ma, which is consistent with previous DNA-only studies showing that underparameterized models or undersampled data sets can inflate divergence time estimates (Schenk and Hufford, 2010; Mulcahy et al., 2012).
Simões et al. (2018) employed essentially the same analytical strategy, but with a data set containing significantly more fossil taxa, particularly stem lineages. They tested calibrating the tree based on both tip-dated fossils only, as well as internal node plus tip age calibrations. They recovered a slightly older, Late Triassic age for crown Squamata at ∼210 Ma, and ∼257 Ma for the MRCA of all squamates (total group Squamata). Garberoglio et al. (2019a) essentially repeated the analyses of Simões et al. (2018) with the addition of more fossil and extant snakes, yielding nearly identical results. However, in the latest expansion of this data set, Simões et al. (2020b) employed approaches to avoid overestimating divergence times owing to deep root attraction in total evidence dating (Ronquist et al., 2016) with the use of two distinct software packages. The best performing model from Simões et al. (2020b) estimates times for the age of the MRCA of crown squamates congruent with the preferred analysis of Pyron (2017) implementing the Mka model at ∼186.5 Ma. Additionally, the age for the MRCA of all squamates was brought down to ∼250 Ma (immediately after the Permian–Triassic mass extinction) and the timing for the split between squamates and rhynchocephalians at ∼259 Ma.
Thus, the prospects for total evidence dating to converge widely accepted estimates for squamate origins seem promising, especially when appropriate models are used to handle both morphological and molecular data. Chief among these are large-scale morphological matrices sampling more fossil taxa with well-constructed characters that are comparable in extant species (e.g., Simões et al., 2017b). Another crucial avenue of investigation will be more accurately parameterized models for morphological evolution (e.g., Harrison and Larsson, 2015; Bapst et al., 2016; Matzke and Wright, 2016; Wright et al., 2016; Pyron, 2017; Simões et al., 2020a). Together, these total evidence dating approaches hold great promise for simultaneously resolving the topology of early squamate groups, as well as estimating their timescale.
One of the most biologically intriguing developments in our recent understanding of squamate evolution has been the functional assemblage of salivary proteins (i.e., venoms) in Toxicofera. Human knowledge and fear of snake venoms is prehistoric and may have contributed to the evolution of primate vision systems (Isbell, 2006). Knowledge of venom in helodermatid lizards is much more recent (Woodson, 1947). Anecdotally, belief that various other lizard species are venomous, including anguids, skinks, geckos, chameleons, and pleurodontans, is widespread in folklore around the world. Squamate venoms generally consist of ordinary and abundant somatic proteins that undergo subsequent biochemical modifications that increase their toxicity (Hargreaves et al., 2014a). This process likely results either from genomic duplication and physiological recruitment from body tissues into modified salivary glands in the jaws or previous occurrence of these proteins in the salivary glands and subsequent restriction and localization into these glands therein. But how many times has this occurred?
Recruitment of toxins was previously thought to represent a more recent evolutionary event restricted to the classically “venomous” snakes in the families Elapidae and Viperidae, with additional convergent origins in a few medically significant colubrids. However, a series of biochemical analyses concurrent with the development of molecular hypotheses for squamate origins found widespread presence of “venom” proteins across snakes (Fry et al., 2003). Subsequently, and contemporaneous with the expanded results of Vidal and Hedges (2005), Fry et al. (2006) announced the discovery of apparently functional venom genes, venom glands, and venoms across Anguimorpha, Iguania, and Serpentes, comprising Toxicofera (gk. “those who bear toxins”), or the “venom clade.”
This discovery implies a single early origin of “venom” in the MRCA of Toxicofera ∼170 Ma, with varying subsequent selection for toxicity and delivery apparatus. Although most of the .4,600 species in this clade are not medically significant to humans, bioactive salivary proteins may thus play a much larger role in prey capture than previously thought. For instance, the notably noxious and medically significant bite of the Komodo dragon (Varanus komodoensis) was long thought to be related to toxic oral bacteria obtained from eating carrion. It is now better understood to result (at least in part) from toxicoferan venoms (Fry et al., 2009), representing a second major lineage of medically significant lizards (Koludarov et al., 2017); even the acrodont iguanian Pogona, a common species in the pet trade, has a venom-secreting glandular apparatus in its jaws (Fry et al., 2006). In contrast, most other lizard species appear to have lost these apparatus and toxicity entirely.
However, this interpretation of an early single origin of “venom” with a unique purpose for prey capture and a wide variety of subsequent evolutionary outcomes in terms of toxicity and complexity of apparatus is not without controversy. Broader sampling of species and tissue types reveals widespread presence of many “venom”type genes in nonfunctional roles across nontoxic structures (Hargreaves et al., 2014b). Even within snakes, many proposed toxins are present across a variety of tissues, including salivary glands, and may thus have been restricted to those glands during independent origins of venom (Hargreaves et al., 2014a). There is also significant disagreement about the structure, function, and terminology of venoms and their apparatus within and among species (see Weinstein et al., 2012 and subsequent replies).
Thus, it may be more accurate to think of the necessary genomic and physiological substrate for venoms (a particular suite of physiologically active proteins) to have evolved early in the history of Toxicofera. In contrast, the several lineages with generally agreed “venoms” and associated apparatus may represent multiple independent parallel origins from this shared substrate. Resolving this question is complex, requiring multiple lines of evidence, including genomics, biochemistry, physiology, histology, functional anatomy, phylogeny, and paleontological evidence. Well-defined ontologies regarding exactly what we consider “venom” and “venomous” are also a necessary complement to these data. Regardless, venom systems in several groups of toxicoferan squamates are among the most consequential and charismatic evolutionary innovations in the Tree of Life. A necessary component to this understanding, however, is an accurate resolution of the Squamate Tree of Life itself, with which to estimate the evolutionary history of these processes.
The Origin of Snakes
The placement of snakes within the squamate Tree of Life and the origins of the extreme morphological adaptations that characterize the snake body plan have long been a hot topic of debate among evolutionary biologists. Morphological data have long been considered limited in their capacity to understand the phylogenetic origins of snakes, owing to the rampant occurrence of apparent homoplasies uniting limb-reduced and skull-miniaturized taxa, discussed in detail above. Although molecular data have provided considerable support for the placement of snakes with iguanians and anguimorphans within Toxicofera, they have not been able to indicate their exact systematic placement within the latter.
Despite limitations of their own, many studies combining morphological and molecular data with extensive sampling of extant and fossil taxa have placed snakes within Toxicofera, and as the sister clade to Mosasauria, forming the clade Pythonomorpha (Reeder et al., 2015; Pyron, 2017; Simões et al., 2018; Garberoglio et al., 2019a). Accordingly, Pythonomorpha is the sister clade to Iguania+Anguimorpha (Reeder et al., 2015; Simões et al., 2018; Garberoglio et al., 2019a). Given the ambiguous molecular support for the placement of snakes within toxicoferans, morphological data thus play a considerable role in inducing the latter configuration. As such, the choice of the morphological matrix will affect the placement of snakes within toxicoferans, and variations on the latter pattern may occur (e.g., Pyron, 2017). Therefore, the continual assessment and reevaluation of morphological characters and fossil specimens remains critical.
In the last decade, important advances have been made for understanding the homology of components of the snake phenotype, aspects that have historically had an important effect in the placement of snakes on the basis of morphological data. Those aspects include: homology of the “postorbital element,” now recognized as the jugal bone instead of a postorbital or postfrontal (Palci and Caldwell, 2013; Garberoglio et al., 2019a); homology of the cervical and caudal peduncles, now recognized as homologous to the hypapophyses and haemapophyses of other squamates; and the absence of a crista circumfenestralis enclosing the stapes in the earliest snake fossils (Palci and Caldwell, 2013; Garberoglio et al., 2019a,b). New and more complete fossils of Cretaceous snakes have also continually expanded our understanding of the acquisition of various features characteristic of the body plan of crown snakes (Garberoglio et al., 2019a). Those findings have provided increasing morphological support for the (Pythonomorpha [Iguania+Anguimorpha]) configuration of Toxicofera.
However, the sister-group relationship between snakes and mosasaurs estimated by many recent analyses creates an additional problem in squamate evolution: a ghost lineage of nearly 40 million years exists between the snake–mosasaur split at ∼171 Ma (Simões et al., 2018; Garberoglio et al., 2019a) and the earliest mosasaurians in the fossil record ∼132.9–125 Ma (Campbell Mekarski et al., 2019). This gap suggests that either further rearrangements of this topological configuration will emerge on the basis of new fossil and morphological data, or that a major gap in the early fossil record of mosasaurians remains to be filled. Simões et al. (2020b) recently applied diversification priors that penalize excessively long ghost lineages (which negatively affect divergence time estimates), recovering mosasaurians as the sister lineage to anguimorphs instead of snakes, despite a strong morphological signal linking mosasaurians and snakes. This investigation indicated that available models are capable of efficiently penalizing such long ghost lineages, despite a strong signal in the data (see also Pyron, 2017). However, because the earliest members of the mosasaurian lineage were, at some point, most likely inhabiting a terrestrial environment, it is possible that such long ghost lineage could be the result of the poor fossil record of terrestrial squamates during most of the Mesozoic until the late Early Cretaceous (Cleary et al., 2018; Close et al., 2019). Further sampling of Late Jurassic and Early Cretaceous squamates is essential to elucidate this matter.
Despite an apparent trend towards the consensual placement of snakes in the Squamate Tree of Life, those new results provide little insight regarding the environmental drivers of the early evolution of the snake body plan. The long and heated debate between a fossorial versus aquatic hypotheses for the origin of snakes (Lee, 1998; Rieppel and Zaher, 2000; Lee, 2005a,b; Gauthier et al., 2012; Hsiang et al., 2015) inherently depends not only on understanding the sister group relationships between early snakes and other squamates, but also a good understanding of the environments where the earliest snakes first evolved and their community structure. This knowledge is currently lacking, in that the oldest known snake fragments come from a variety of environments such as coal swamps, mixed coastal lake and pond systems, and epicontinental islands (Caldwell et al., 2015), and little is known of their ecological context and co-occurrence with other ecologically relevant species. If consensual answers are ever to emerge to settle this debate, they will depend on an increased knowledge of the fossil record and paleoenvironments of mid-Jurassic snakes and other squamates, sources of data outside the strict confines of molecular or morphological systematics.
The phylogenetic origin of snakes has also stirred broad avenues of research towards understanding the origins of the extreme adaptations characterizing the snake body plan. Although limb reduction is widespread among squamates (Wiens et al., 2006), most groups retain at least partial development of their limb skeletons, with the major exception being snakes. Although some relatively early diverging snakes retain remnants of their pectoral girdle (e.g., “henophidians”, such as boas and pythons), the vast majority of snakes have lost their limbs entirely.
Recently, genomic studies have revealed fundamental clues about the genetic mechanisms responsible for the loss of limbs in snakes. The loss of limbs in squamates leading to the origin of snakes has been demonstrated to stem from mutations on an enhancer—the zone of polarizing activity regulatory sequence (ZRS)—of the sonic hedgehog gene (Shh). It was already known that this was a conserved regulatory region (Sagai et al., 2004) and that mutations on this enhancer can result into limb malformation in vertebrates (Sagai et al., 2005; Hill and Lettice, 2013), and so it would be expected that snakes would have a high number of mutations on the ZRS. Surprisingly, however, this enhancer has been observed to be extremely conserved across extant early-evolving snakes retaining remnants of the hindlimbs (e.g., boas and pythons) (Kvon et al., 2016). Only later-evolving snakes without limb remnants demonstrate a much higher degree of substitutions on this enhancer, suggesting a progressive loss of function of the ZRS in snake evolution (Kvon et al., 2016). This hypothesis of a more gradual evolutionary loss of limbs in snakes has been later supported by new fossil evidence and phylogenetic reconstructions suggesting that early snakes retained relatively well-developed hindlimbs for at least the first 80 million years since their origin (Garberoglio et al., 2019a). Combined, these results suggests that early snake evolution is not characterized by drastic changes on their limbs, but rather, on the skull (see also Caldwell, 2019). This idea finds further support from recently detected high rates of morphological evolution of the snake skull among the earliest branches of the snake evolutionary tree (Watanabe et al., 2019) and early snake fossils already possessing mandibular features extremely very similar to modern day taxa (Caldwell et al., 2015).
Phylogenetic Signal and Fossil Data
Despite the long-term debate between the role of morphology and molecules in systematics (e.g., Scotland et al., 2003; Wiens, 2004), this issue was mostly concerned with the quality of such data sets per se for systematic purposes (regardless of analytical procedure) or focused only on more traditional analytical approaches such as maximum parsimony. However, the increasing availability of phylogenomic data sets and the emergence of total evidence dating with relaxed clock Bayesian inference under rapidly improving analytical procedures, makes it important to assess the effect of and the best approach towards combining morphological (from fossils and extant taxa) and molecular data to reconstruct the squamate Tree of Life.
As illustrated in the sections above, evidence is ample for the limited role of morphological data and fossil taxa in driving estimated topological relationships in squamates when combined with moderate to large-sized molecular data—for either non-clock or relaxed clock total evidence dating. This situation may only be significantly altered by character weighting to balance out the overwhelmingly large number of molecular characters (Giribet, 2010), which has thus far not been attempted in large squamate data sets. However, even in the absence of character weighting, we argue that, despite such limitations, morphological data still play an important role for the placement of both extant and extinct taxa in the instances in which the molecular signal is weak or ambiguous (see Wiens et al., 2010; Pyron, 2017). A clear example comes from the ambiguous molecular signal described above concerning the earliest dichotomy in the squamate tree (dibamids and gekkotans) and over the relationship of the three major toxicoferan clades. In such situations, the morphological signal drives the placement of gekkotans as the earliest evolving crown squamate lineage relative to dibamids on the basis of the plesiomorphic features of some gekkotans in the context of early diapsid and early lepidosaur evolution, and which are not present (or inapplicable) in dibamids (see above; Simões et al., 2018).
This signal is so strong that even other morphological data sets without a broad sampling of early diapsids and lepidosaurs (including only a handful of sphenodontian taxa for outgroup comparison) also recover gekkotans as diverging earlier relative to dibamids when combined with molecular data, using both the Conrad (2008) and Gauthier et al. (2012) data sets tested by Pyron (2017). The morphological signal may also play an important role in the placement of snakes within Toxicofera given the ambiguous molecular signal, although this has yet to be formally tested. Furthermore, morphological data play an important role even in the presence of extremely large molecular data sets by estimating the relationships within entirely extinct lineages, such as mosasaurians, borioteiioids, and priscagamids, among many others, which will also affect their sister group relationships to extant clades (see Wiens et al., 2010; Pyron, 2017). Finally, the inclusion of both morphological and molecular data has a fundamental role in total evidence dating, by avoiding ad hoc hypotheses for the placement of fossil taxa and further avoiding a priori constraining of topological relationships (Ronquist et al., 2012; Zhang et al., 2016).
Another important point of consideration is that a great portion of the Tree of Life, especially in early-diverging branches, is represented only by fossils, with some accounts estimating that up to 99.9% of all species that ever lived are now extinct (Novacek and Wheeler, 1992). This observation implies an inherent property of phylogenetic trees of clades with extant representatives: The broader the taxonomic sampling, or the deeper in time we need to go to find the last common ancestor of all sampled species, the greater the likelihood of encountering entirely extinct lineages. As noted above, the phylogenetic relationships of these extinct lineages can thus only be informed by morphological (usually osteological) characters and data sets.
A fundamental corollary of this property is the prediction that the broader or deeper in time the taxonomic sampling (e.g., higher level phylogenies) extends, the greater the chances of creating artifactual sampling gaps if fossils are not included (e.g., compare Figs. 5 and 6). We call this property the “deep-time sampling bias.” Thus, many, if not most, “long branches” in extant-only trees (e.g., the branch leading to Sphenodon) are not actually long; they are simply undersampled, creating increasingly larger taxon sampling gaps irrespective of the sampling effort of extant taxa. Concomitantly, “breaking” these long branches by sampling additional fossil taxa is well known to exert a strong and presumably positive effect on phylogenetic inference and dating (Gauthier et al., 1988b; Wiens, 2004; Guillerme and Cooper, 2016; Luo et al., 2020).
In the case of squamates, the lineages representing most extant genera and families coalesce at some point in the last 75–100 million years (Pyron, 2017; Simões et al., 2018; Burbrink et al., 2020; Simões et al., 2020b), but the divergence between the major squamate higher clades is much more ancient, concentrated between 75 and 175 Ma (Table 1; Pyron, 2017; Simões et al., 2018; Burbrink et al., 2020; Simões et al., 2020b), sometimes marked by some long stem branches. As expected, the 100 million years between the origin of the higher squamate clades and the divergence of most extant families represent a timeframe densely interpermeated by extinct families (e.g., Mosasauria, Borioteiioidea, Paramacellodidae; Figs. 1J–N) and also many stem members of extant lineages (e.g., Eichstaettisaurus and Gobekko for gekkotans; Najash and Dinilysia for snakes; Figs. 1D, E). Fossil squamates ranging in age from ∼75 Ma to the time of the oldest known fossil squamate (242 Ma) are thus of particular relevance to sample several extinct lineages and break ancient long branches and avoiding deep fossil sampling biases.
Therefore, genomic data are most likely to provide highly accurate and robust phylogenetic hypothesis relationships among younger or more recent lineages given the recent radiation of the sampled species (e.g., family or subfamily level). However, higher level phylogenies have a much greater need for sampling of extinct lineages for a balanced taxonomic representativeness, thus necessitating a morphological partition for the inclusion of fossil taxa. Recent studies demonstrate the latter to be fundamental for improving accuracy for divergence time estimation (Guillerme and Cooper, 2016; Luo et al., 2020). Although such an effect was already known in some sense (e.g., Graybeal, 1998; Wiens, 2003), the particular example of the squamate Tree of Life throws it into sharp relief for broadscale phylogenetic estimation of temporally diverse, ancient groups.
Genomics is still in relative infancy for evolutionary and phylogenetic studies, and whole-genome data present a wealth of data for resolving relationships among lineages and clarifying the evolutionary origin of traits that have not yet been fully deployed for these purposes, particularly in squamates. For instance, analysis of gene order conservation (synteny mapping) represents a promising source of phylogenetic signal—essentially measuring the “morphology” of the genome. Analyses thereof have resolved contentious branches in prokaryotes (Shifman et al., 2014) and mammals (Luo et al., 2012). Other sophisticated methods of analyzing whole-genome data include topology-by-location scans, wherein genealogies are estimated for loci across the whole genome, showing complete patterns of gene tree variation (Li et al., 2019).
One can thereby highlight the empirical distribution of discordant gene trees across areas of the genome that are duplicated, introgressed, or vary in effective population size, giving more decisive estimates of the overall species tree (e.g., Ravinet et al., 2018). Such approaches might be particularly applicable to longer internal or terminal branches, such as the placement of Dibamia or amphisbaenian monophyly, by revealing which topologies are linked most closely to highly conserved regions of the genome (see Singhal et al., 2021).
Other sources of phylogenetic signal abundant in whole-genome data are mobile genetic elements such as SINEs and LINEs, which have long been noted for their potential applicability to difficult phylogenetic problems (Hillis, 1999), but rarely deployed at scale. Genomic scans can produce a wealth of such data as bycatch from other types of loci, offering a corroborative counterpoint to gene tree heterogeneity because of short internal branches and incomplete lineage sorting (Cloutier et al., 2019). Such approaches may be confounded by the complex diversity and history of such elements in squamate genomes (Pasquesi et al., 2018). Nevertheless, element-based approaches might be helpful in discriminating among topologies in the anomaly zone of rapid radiations, such as Pleurodonta and booid, colubroid, and elapoid snakes.
Genomic data are, of course, also indispensable in understanding the origin, evolution, and development of complex traits, in conjunction with detailed morphological and physiological observations. Myriad examples are present in the recent literature, such as the origins of visual and olfactory adaptations for prey capture (Perry et al., 2018), venom evolution (Schield et al., 2019), metabolic adaptations (Lind et al., 2019), sex chromosomes (Gamble et al. 2017), and viviparity (Gao et al. 2019). Genomic assemblies of squamate species are becoming increasingly common for snakes (Kerkkamp et al., 2016) and annotated drafts available for some lizards species within iguanians (Alföldi et al., 2011), anguimorphans (Song et al., 2015), and geckos (Xiong et al., 2016). As of 28 May 2020, more than 241 genome-scale data sets are present on GenBank, providing a wealth of data for both phylogenetic inference and comparative analysis.
Understanding the overall relationships among the numerous lineages of living and extinct squamates has been one of the major goals in herpetology and evolutionary biology for over a century. This long-lasting problem has generated fierce debates across disciplines towards understanding some of the most fundamental questions regarding the evolution of lizards, snakes, amphisbaenians, and tuataras. However, the last two decades have been especially promising toward clarifying some of the most disputed aspects of squamate evolution.
Squamate phylogenomics has revealed that nearly all genomic regions result in very similar species trees, with only a few points of standing disagreement. Importantly, agreement between morphological and molecular hypotheses is starting to emerge. Those hypotheses, along with combining morphological and molecular data in total evidence dating, have enabled the development of important points of increasingly consensual agreement and synthesis in squamate evolution. We propose that major points of synthesis in squamate evolution are represented by a much improved understanding of the early fossil record and origin of squamates from other reptile lineages—a remarkable convergence of estimates of diversification times for the major lineages of squamates from total evidence dating, monophyly of Toxicofera, and the placement of snakes as the sister lineage to iguanians+anguimorphans, irrespective of the placement of mosasaurians. Other emerging patterns, but which need further corroboration, include monophyly of amphisbaenians and paraphyly of scolecophidians.
The major standing problems in squamate phylogenetics that remain to be solved by larger scale phylogenomics and morphological data sets (most crucially, including more fossil lineages) are the internal relationships of pleurodont iguanians and colubroid snakes, the placement of dibamids, and the placement of many fossil families (e.g., mosasaurians and paramacellodids). The relationships of iguanians, colubroids, and dibamids will be resolved, one hopes, by genome-scale data sets, as well as possibly by the discovery and inclusion of new fossil taxa in morphological matrices. Estimating fossil relationships and more accurately dating and reconstructing deep branches in the squamate tree depend on expanding our currently scarce knowledge of the squamate fossil record during the first half of squamate evolution: from the Middle Triassic to the Early Cretaceous. Perhaps more than any other factor, this is the most essential source of data to understand not only early squamate taxonomy, but also to improve our ability to break long stem branches and reconstruct early patterns of squamate evolution and the time of origin of its major clades.
As evidence accumulates regarding the benefits of combining morphological and molecular data for inferring accurate time-trees, the future of estimating higher level evolutionary trees seems to be inherently dependent on accurate estimates of both morphological and molecular evolution. Therefore, efforts should be put into improving models and protocols to assemble and analyze such data sets, both independently as well as combined in total evidence dating. Finding agreement between both data types is an important starting point. It may be impossible to obtain a full agreement between morphological and molecular trees, because they represent fundamentally distinct dimensions of organismal evolution and are differentially affected by selective pressures and by fixed and plastic variation in development and ontogeny. However, at least partial levels of agreement are achievable and necessary to avoid conflicting phylogenetic signals that negatively affect timetree inference. The development of more complex models of morphological character evolution is paramount to a more biologically realistic modeling of the phenotype, comparable to the array of molecular substitution models currently available.
AMNH, American Museum of Natural History, New York, New York, USA; CM, Carnegie Museum of Natural History, Pittsburgh, Pennsylvania, USA; FHSM, Fort Hays Sternberg Museum, Hays, Kansas, USA; FMNH, Field Museum of Natural History, Chicago, Illinois, USA; MCZ, Museum of Comparative Zoology, Harvard University, Cambridge, Massachusetts, USA; MPCA, Museo Carlos Ameghino, Cipolletti, Río Negro Province, Argentina; NHMUK, Natural History Museum, London, United Kingdom; NMNH, National Museum of Natural History, Washington, D.C., USA; PZO, Museo Archeologico dell'Alto Adige, Bolzano (Bozen), Italy; TMP, Royal Tyrrell Museum, Drumheller, Alberta, Canada; UAMZ, University of Alberta Museum of Zoology, Edmonton, Alberta, Canada.
We thank several researchers with whom we have had thoughtful and intense debates regarding the squamate Tree of Life, despite not always holding the same views, especially M. W. Caldwell, R. Nydam, M. S. Y. Lee, J. Gauthier, J. C. Rage, J. J. Wiens, D. L. Rabosky, P. O. Title, S. Singhal, F. T. Burbrink, H. Zaher, and J. Conrad (R.I.P.). This research was funded in part by U.S. NSF grant DEB-1441719 to R.A.P. and an Alexander Agassiz Postdoctoral Fellowship to T.R.S. provided by the Museum of Comparative Zoology, Harvard University.