The Jurassic rise of squamates as supported by lepidosaur disparity and evolutionary rates

Arnau Bolet; Thomas L Stubbs; Jorge A Herrera-Flores; Michael J Benton

doi:10.7554/eLife.66511

. 2022 May 3;11:e66511. doi: 10.7554/eLife.66511

The Jurassic rise of squamates as supported by lepidosaur disparity and evolutionary rates

Arnau Bolet ^1,^2,^✉, Thomas L Stubbs ², Jorge A Herrera-Flores ², Michael J Benton ²

Editors: Min Zhu³, George H Perry⁴

PMCID: PMC9064307 PMID: 35502582

Abstract

The squamates (lizards, snakes, and relatives) today comprise more than 10,000 species, and yet their sister group, the Rhynchocephalia, is represented by a single species today, the tuatara. The explosion in squamate diversity has been tracked back to the Cretaceous Terrestrial Revolution, 100 million years ago (Ma), the time when flowering plants began their takeover of terrestrial ecosystems, associated with diversification of coevolving insects and insect-eating predators such as lizards, birds, and mammals. Squamates arose much earlier, but their long pre-Cretaceous history of some 150 million years (Myr) is documented by sparse fossils. Here, we provide evidence for an initial radiation of squamate morphology in the Middle and Late Jurassic (174–145 Ma), and show that they established their key ecological roles much earlier than had been assumed, and they have not changed them much since.

Research organism: Other

Introduction

Lepidosaurs, currently mainly represented by squamates, are one of the most species-rich tetrapod clades (Uetz et al., 2019), only rivaled by birds in terms of diversity. Evidence points to an explosion in squamate biodiversity 100 million years ago (Ma), in the Cretaceous (Lloyd et al., 2008; Longrich et al., 2012; Mongiardino Koch and Gauthier, 2018), corresponding to the time when flowering plants were diversifying and restructuring terrestrial ecosystems. However, the origin of squamates at least 250 Ma (Jones et al., 2013; Simões et al., 2018) poses two challenges: was that 150 million years (Myr) of poorly documented history real or evidence of a poor fossil record; and when did squamates acquire their current range of ecological adaptations? These uncertainties contrast somewhat with more confident knowledge of the early radiations of birds (Lee et al., 2014; Wang and Lloyd, 2016; Benson et al., 2014) and mammals (Close et al., 2015).

The squamate fossil record is patchy, especially through the Triassic to Early Cretaceous interval (252–100 Ma) when fossils are sparse and incomplete (Cleary et al., 2018). The earliest unambiguously identified squamate fossils date from the Middle and Late Jurassic (174–145 Ma), and among them are forms that can be assigned to major modern clades of squamates, including both lizards and snakes (Evans, 1998; Evans, 2003; Caldwell et al., 2015), but many are isolated jaws and skull bones of difficult identification. Further, their rarity suggests diversity was not high.

The mid-Cretaceous shows an increase in abundance and diversity of squamates (Gauthier et al., 2012), linked to the Cretaceous Terrestrial Revolution (KTR), which triggered an outburst of terrestrial life, including major new clades, such as angiosperms, as well as ferns, hornworts, liverworts among plants, currently highly diverse insect groups, including cockroaches, termites, many groups of beetles, bugs, the wasp, ant, and bee lineage, and the butterfly and moth lineage. These rich new supplies of plants and insects provided food for expanding clades of insect-eaters, including spiders, birds, mammals, and lizards, and even perhaps some dinosaur groups (Lloyd et al., 2008; Doyle, 2008; Benton, 2010; Meredith et al., 2011; Cardinal and Danforth, 2013). It could be plausible to identify this as the time when squamate ecological adaptation expanded, but is this truly the case?

Here, we explore morphospace distribution, disparity (morphological richness), and evolutionary rates of lepidosaurs to understand these important stages in the first three-quarters of squamate history. Species richness is hard to document with confidence in the face of such a patchy fossil record (Cleary et al., 2018), although when combined with phylogenomic data, the relative timing of origins of major modern clades can be identified (Jones et al., 2013; Simões et al., 2018). Sparse fossil data can be used, on the other hand, to document disparity, even though extreme forms may be absent, and the morphospace may be poorly filled. Phylogenetic comparative methods are used here to explore whether the available fossil record before the Late Cretaceous suggests that the importance of the Jurassic in the rise of squamates has been underestimated, hinting at a cryptic diversity hidden behind an impoverished fossil record.

Results

Our dated phylogeny of lepidosaurs (Figure 1, Figure 1—figure supplement 1), based on a morphological tree constrained by phylogenomic evidence (Figure 1—figure supplements 2–4, see Figure 1—figure supplements 5–7 for results of an unconstrained analysis), shows that the clade originated around the Permian-Triassic boundary, and that by the Mid-Triassic it was represented by different extinct groups. The Rhynchocephalia diversified in the Mesozoic, but reverted to a single species subsequently, the living tuatara (Sphenodon punctatus) from New Zealand. The Squamata, on the other hand, show a step of diversification through the Jurassic, as the main modern clades emerged, and then a further diversification in the mid-Cretaceous, perhaps linked to the KTR. The ranges of ages for nodes differ depending on the method used (Figure 1—figure supplement 1), with the Hedman method (Figure 1) yielding the oldest divergence dates for nodes, the MBL (Minimum Branch Length) method yielding the youngest ones, and the equal method being intermediate. According to these results, it was regarded as necessary to consider all three methods of dating when performing the evolutionary rate analyses in order to discard the possibility that results were biased by the selected dating method.

The dated phylogenetic tree only describes the outline of the origins of squamate biodiversity, but does not map species numbers or, importantly, the range of morphology, and presumably the range of adaptation, reflecting ecological impact, of the group. Using a large morphological dataset (Conrad, 2018), covering 201 species of living and fossil lepidosaurs scored for 836 skeletal morphological characters, we analyzed disparity for lepidosaurs through time and tracked changes to morphospace occupation and major expansions. Stacked temporal morphospaces (Figure 2) show that rhynchocephalians and squamates occupy mutually exclusive morphospaces. Stepping up through time, from bottom to top of the stack, shows how lepidosaur morphospace expanded, not gradually, but marked by a single major step. At first, the total morphospace is small, formed by stem lepidosaurs and rhynchocephalians in the Triassic, and exclusively by rhynchocephalians in the Early Jurassic. Note that the lack of Early Jurassic stem lepidosaurs is artificial because they are found again (and possibly for the last time) in the Middle Jurassic as represented by Marmoretta. Then, with the addition of squamates in the Middle to Late Jurassic bin, morphospace expands to five or six times the area – the limits are established by rhynchocephalians, generalized lizards, anguimorph lizards, and snakes, each occupying a separate area of morphospace. We coin the term Jurassic Morphospace Expansion (JME) for the event related to this sudden increase in morphospace, which is interpreted as evidence of the initial radiation of the total group Squamata. Note that the coincident loss of stem lepidosaurs does not result in a modification of the morphospace hull because of the central position of these taxa. This morphospace configuration remains remarkably stable from the Late Jurassic through to the present, with only subtle increases in morphospace occupation and in the density of points inside the envelope, notably in the mid Cretaceous coinciding with the KTR and the consequent recorded increase in diversity.

Figure 2. — Blue circles correspond to squamates, with the blue scincid silhouette indicating the position of generalized lizards, the yellow varanid indicating the position of anguimorphs, and the violet snake the position of snakes (and other limbless squamates). Green triangles correspond to rhynchocephalians (green *Sphenodon* silhouette). Red crosses correspond to stem lepidosaurs (red kuehneosaur silhouette). For additional plots of morphospace occupation through time, see Figure 2—figure supplements 1–6. JME: Jurassic Morphospace Expansion; KTR: Cretaceous Terrestrial Revolution.

Figure 2—figure supplement 1. — Blue circles correspond to squamates, with the blue scincid silhouette indicating the position of generalized lizards, the yellow varanid indicating the position of anguimorphs, and the violet snake the position of snakes (and other limbless squamates). Green triangles correspond to rhynchocephalians (green *Sphenodon* silhouette). Red crosses correspond to stem lepidosaurs (red kuehneosaur silhouette). For additional plots of morphospace occupation through time, see Figure 2—figure supplements 1–6. JME: Jurassic Morphospace Expansion; KTR: Cretaceous Terrestrial Revolution.

At this point, we should comment on the form-function relationship. It is well understood that form (skeletal morphology) does not always equate to function (Wainwright et al., 2005), with many functions sometimes performed by organisms of apparently similar morphology, or many different morphologies capable of performing a single function. However, here we calibrate the morphospace by mapping living taxa of known function and ecology onto the fossil time slices. This means we can mark (Figures 2 and 3) the rhynchocephalian (cluster 1, in green) pole as dietary generalists with robust jaws and tongue prey prehension like the modern Sphenodon; the generalized lizard (cluster 2, in blue) pole as diverse insect-eaters, like modern skinks; the anguimorph (cluster 3, in yellow) pole as active foragers with tendency to carnivory; and the snake (cluster 4, in violet) pole as limbless predators that feed mainly on other vertebrates and, to a lesser degree, invertebrates. Note that these clusters are loosely based on the ones recovered using the R ‘pamk’ function of the fpc package (see Figure 3—figure supplement 6, Source code 3). Admittedly, this is a simplification because, just as an example, rhynchocephalians contain dietary specialists (durophagous, piscivorous, etc.) and examples of forms adapted to swimming (like pleurosaurs), and the same occurs with specific clades of squamates (e.g., herbivore iguanians). Although a more precise ecomorphological classification would potentially provide more information on the distribution of ecologies through time, it has proven impractical for the current dataset, also in agreement with previous results (e.g., Simões et al., 2020). However, our classification serves the purpose of depicting that the extremes of morphospace had been achieved by the Middle to Late Jurassic (Figures 2 and 3D). Regarding squamates, it is worth noting that the first members of the group sampled in our analyses are Late Jurassic in age, but Middle Jurassic forms are known, just happen to be not included in the current morphological dataset because they are too incomplete.

The illustrated morphospaces (Figures 2, 3A and B, Figure 2—figure supplement 2, Figure 3—figure supplement 1) represent the first two major axes of variation, and there could be additional morphospace expansions along the other main axes: this is not the case (Figure 2—figure supplements 3 and 5; Figure 3—figure supplements 2 and 3). Further, the story does not change when the post-Cretaceous time bin is divided into Paleogene and Neogene time slices (Figure 2—figure supplements 2; 4 and 6). In a plot of total disparity (i.e., the sum of variances [SoVs] across all morphospace axes) of lepidosaurs through geological time (Figure 4A), the two peaks of elevated disparity (Middle–Late Jurassic and mid–Late Cretaceous) are clear. These summary data also confirm the much higher total disparity of squamates than rhynchocephalians (Figure 4—figure supplement 1) and among the former, a higher disparity of snakes and anguimorphs (and less clearly mosasaurs) than the remaining main groups of squamates among which dibamids, lacertids, and the extinct group of ardeosaurs (sensu lato) are the ones with the lowest disparity (Figure 4—figure supplement 2).

Figure 4—figure supplement 1. — (A) Temporal disparity patterns (bootstrapped and rarefied within bin sum of variances for all axes). For additional plots of disparity patterns, see Figure 4—figure supplements 1–4. (B) Evolutionary rates through time in epoch-scale bins. Black solid line corresponds to results for all taxa, blue solid line for lizards, and green solid line for rhynchocephalians plus stem lepidosaurs (all according to the constrained phylogeny). Dashed gray line corresponds to results for all taxa and unconstrained phylogeny. The curves represent averages from 25 iterations of each analysis using randomly selected trees dated with the Hedman method. For additional plots of evolutionary rates, see Figure 4—figure supplements 5–22. JME: Jurassic Morphospace Expansion; KTR: Cretaceous Terrestrial Revolution; ETr: Early Triassic; MTr: Middle Triassic; LTr: Late Triassic; EJ: Early Jurassic; MJ: Middle Jurassic; LJ: Late Jurassic; EK: Early Cretaceous; LK: Late Cretaceous; Pg: Paleogene; Ng-Q: Neogene-Quaternary.

The ‘two-peak’ pattern is also identified through a study of evolutionary rates among lepidosaurs (Figure 4B). Maximum-likelihood analyses of rates of morphological character evolution clearly show that the highest rates occur in the Late Jurassic (coinciding with the peak in disparity and roughly coinciding with the observed expansion of morphospace), and there is a lower peak in the mid-Cretaceous (see Figure 4—figure supplements 5–7). The high peak is for squamates (see Figure 4—figure supplements 8–10 for squamates), not rhynchocephalians (see Figure 4—figure supplements 11–13), and in our opinion all three observations (morphospace expansion, increase of disparity, and fast evolutionary rates) are linked and support the existence of the JME event. Rhynchocephalians showed their highest rates of evolution in the Triassic, and those rates declined substantially through time. On the other hand, squamates show a further step-up in rates during the Neogene, the past 23 Myr. The peak in evolutionary rates towards the present is particularly acute when using the MBL method, in what we regard as an artifact due to the condensing of divergence dates towards the present. These evolutionary rates are not the result of choosing a specific phylogenetic context – the black line represents summed rates for all lepidosaurs using the phylogenomic constrained trees, and the dashed gray line the rates from unconstrained morphological trees, with iguanians and fossorial species in traditional positions. This alternative version yields a similar general result, where the Late Jurassic peak is again clearly recovered (see Figure 4—figure supplements 14–22). In both cases, most parsimonious trees (MPTs) were randomly selected and dated five times according to each of the three dating methods, and evolutionary rates were calculated for each of the resulting dated trees. The curve shown represents average evolutionary rates among those calculated using the randomly selected MPT’s and the equal method (see Figure 4—figure supplements 5–22 for results according to all methods). The Late Jurassic peak in evolutionary rates is thus also robust to changes in the particular (randomly selected) point inside the stratigraphical range of a given fossil, and to different methods of dating the trees, including ‘equal’ and ‘Hedman’ methodologies, but not the ‘MBL’ method (but see the discussion for a possible explanation).

Discussion

In exploring the nature of the ‘early burst’ in squamate disparity, we wanted to understand how the different clades occupied morphospace. Such an early establishment of squamate morphospace has never been documented in the few studies dealing with early radiations of squamates (e.g., Simões et al., 2020). Our finding that morphospace dimensions were established as early as the Middle or Late Jurassic is a counterintuitive result because it substantially predates the apparent increase in species richness that is usually tied to the KTR. Also, the high morphological rates recovered by a similar approach in Simões et al., 2020 applied to a different dataset of lepidosaurs plus a wide array of other diapsids (Simões et al., 2018) are in their case not correlated to a high disparity. In the latter study, morphospace was not plotted through time, so our time series (Figure 2) or chronophylomorphospace (Figure 3D) have no counterpart in their results.

It is also interesting to observe so little addition to squamate morphospace after its establishment; the great expansion in species numbers up to 10,000 today has happened partly by minor expansions of the total morphospace envelope, but mainly by packing ever more species inside the existing morphospace area. Although we would expect that a more fine-tuned grouping of ecomorphotypes would reveal minor changes in the specific portions of occupied morphospace, it is unlikely that such improvement in resolution would change the recovered outer limits of morphospace.

The total morphospace (Figures 2 and 3A) confirms the central location of stem lepidosaurs, and that rhynchocephalians and squamates explored distinct morphospace throughout their evolution (see their spread along PCO2). Martínez et al., 2021 reported that rhynchocephalians seem to present a lower proportion of morphological traits not shared with other diapsids than squamates do, which they interpreted as supporting the traditional view of rhynchocephalians as retaining a more plesiomorphic morphotype than squamates. Although our dataset is not comparable in the sense that it does not contain such a wide array of non-lepidosaurian diapsids, our morphospace, showing stem lepidosaurs in an intermediate position between squamates and rhynchocephalians, does not seem to support a plesiomorphic morphology for the latter. Our results are in line with current reinterpretations of many supposedly plesiomorphic traits of rhynchocephalians as actually derived in the context of Lepidosauromorpha. Regarding the distribution of squamates groups in particular, in our 2D phylomorphospace (Figure 3B) the squamate clades that lie furthest from the centroid are anguimorphs and marine mosasaurs (forming cluster 3) as well as amphisbaenians (worm lizards), dibamids (blind skinks), and snakes (forming cluster 4, associated with the limbless morphotype). A more detailed distribution of clades through the morphospace is represented in Figure 3—figure supplement 4.

Simões et al., 2020 reported a comparable morphospace for lepidosaurs, where rhynchocephalians and squamates do not overlap in morphospace, and squamates are mainly divided into limbed (lizards) and limbless (snakes and amphisbaenians) taxa. However, their large sample of non-lepidosaur diapsids forced a relatively low sample of lepidosaurs (in comparison to our study) that resulted in a loss of resolution for lepidosaur and squamate morphospace. It is worth noting that our morphospace is also similar to the one recovered by Watanabe et al., 2019 as based on skull geometric morphometrics, in identifying a cluster of snakes (and other limbless taxa like amphisbaenians and dibamids), anguimorphs, and mosasaurs at another pole, and a poor differentiation between other lizard groups (what we regard as ‘generalized lizards’).

The 3D phylomorphospace (Figure 3C, Supplementary files 4 and 5) shows how additional variation along PCO3 reveals the separation of several lizard morphotypes. This morphospace is easier to interpret in the 3D interactive plot (Supplementary file 5), which more clearly separates amphisbaenians and dibamids from snakes, the latter occupying an intermediate position between the two former clades and anguimorphs along PCO3. The inclusion of PCO3 slightly improves separation among ‘generalized lizards,’ although some clades persist as mixed groups in morphospace (e.g., scincids and gekkotans). However, PCO3 separates limbed (e.g., Eumeces and Acontias) from limbless (Typhlacontias and Feylinia) scincids. There is also a superposition between iguanians and the portion of anguimorphs that is closer to the middle part PCO3, which seem to represent less specialized anguimorphs (e.g., Elgaria) than those at the edge (varanid-like forms). In our opinion, all available studies (including ours) fail in fine-tuning ecomorphological groups of squamates, hinting at a problem that seems to be shared to both types of source data (discrete morphological characters and morphometric data). Ours is, however, the first study to track lepidosauromorph morphospace changes through time.

The chronophylomorphospace (Figure 3D) highlights the Mid–Late Jurassic expansion of morphospace, linked with the first radiation of the novel ‘snake-like’ morphology (see also Appendix 3—figure 1), but also that of clearly predatory forms as represented by Dorsetisaurus. In our interpretation, the first event (JME) is tracking the initial radiation of total group Squamata, when the clade radiated into its main components, as revealed by the branching timing of the dated trees, the primary expansions seen in morphospace plots, the Late Jurassic peak in disparity, and the Late Jurassic peak in evolutionary rates. This is particularly true for the trees dated using the Hedman and equal methods, but in the case of the MBL method the signal is overprinted by the high concentration of short branches close to the present time, which result in artificially higher rates for that period, and relatively lower Mesozoic peaks. Besides the MJE, a second event, which fits well with the timing of the KTR, would be coincident with the radiation of the constituent crown groups of Squamata as revealed by a limited expansion and infilling of morphospace plus the record of coincident peaks of disparity and evolutionary rates. Whether this second event represents an actual event of diversification among squamates or is the result of a greatly improved fossil record remains unclear.

Our study offers a new perspective on the early evolution of the major clade Squamata and the other groups of lepidosauromorphs. It benefits from current phylogenomic evidence on phylogeny, as well as fossil data on the timings of events and the expansion of skeletal morphologies and disparity. Our results show that, although the first assemblages of lizards (and possibly snakes) in the Middle–Late Jurassic are not particularly diverse or abundant, the basic structure of the present morphospace distribution had already been achieved (Figures 2 and 3D, Appendix 3—figure 1). This is independent of the interpretation of the affinities of a given taxon because the points in morphospace do not change with changes in topology (only the branches uniting them in phylomorphospaces). Finding support for this early burst of disparity and associated rapid evolutionary rates was rather unexpected, especially so long before the KTR – a reported key driver of squamate evolution (Gauthier et al., 2012) – and before a good fossil record is documented.

Further, we confirm that these distributions in morphospace, marking broad ecological and functional groupings, were remarkably stable for the subsequent 150 Myr, through to the present day (Figure 2, Figure 2—figure supplements 1–6). In other words, the range of adaptations in the current huge diversity of squamate species tracks back very deep in Earth history, some 60 Myr before the KTR. The only observable changes from then on correspond to a slight expansion of the edges of the occupied morphospace, and a notable increase in the density of points filling this morphospace. We acknowledge, however, that the recovered structure represents a simplification that only corresponds to groups according to general bauplans (e.g., limbed vs. limbless morphotypes) and, to a lesser degree, adaptations achieved by specific clades like, for example, anguimorphs. It is thus not possible to track finer ecomorphologies like, for example, adaptations of snakes to different environments (marine, fossorial, or ground-dwelling), which would likely add some variability in the form of shifts in the occupied morphospace through time. Although this is possibly related, in part, to the use of a phylogenetic morphological matrix that was constructed to capture the deep phylogenetic relationships of the constituent groups inside Squamata, it is worth noting that our morphospace is not too different from one recovered from a geometric morphometrics approach (Watanabe et al., 2019), suggesting that this poor resolution is not entirely explained by this procedural choice.

As a final note on squamate morphospace distribution and evolutionary rates, the results presented here also differ from recently published studies dealing with dentition shape, jaw size disparity (as informed by geometric morphometrics Herrera-Flores et al., 2021a), and body size as a continuous character (Herrera-Flores et al., 2021b). In the latter, even though the divergence times for most clades were in line with our results (they applied the Hedman method to their dataset), their results on evolutionary rates greatly differ from the ones presented here. Our study here yields consistently higher rates for squamates than rhynchocephalians, when the opposite trend was recovered in reference (Herrera-Flores et al., 2021b). Moreover, our results show a trend of decreasing evolutionary rates for rhynchocephalians through time, whereas the opposite was recovered by Herrera-Flores et al., 2021b through the Mesozoic, with a marked increase across the Jurassic–Cretaceous boundary. We regard these striking differences as related to the radically different sources of information used in both studies. Although differences in body size can be used to track shifts in evolution and, accordingly, to hint at macroevolutionary patterns, they do not need to be necessarily related to the same processes explored here. Both mentioned studies (Herrera-Flores et al., 2021b; Herrera-Flores et al., 2021a) dealt with particular aspects of lepidosaur evolution (body size and dentition/diet), whereas the results presented here are derived from an approach that considered many different ecomorphological aspects (as many as can be reflected in a morphological matrix that includes osteological characters for the entire skeleton, as well as soft-tissue characters). In addition, the focus on Mesozoic taxa in the aforementioned studies make the results obtained here, which include a good representation of extant taxa, difficult to compare.

It is important to consider whether results could represent bias in the fossil record. It is well understood that the Mesozoic fossil record of squamates is patchy, including some very poorly sampled time intervals (Cleary et al., 2018; Evans, 2003). As a counter to this concern, we note that the occupied total squamate morphospace in the Middle–Late Jurassic is just slightly smaller than that for the Late Cretaceous or the Paleogene to extant time bin (Figure 2), even though the two latter have yielded much higher sample sizes of specimens (Appendix 1—figure 5) that are also anatomically more complete (Appendix 1—figures 1–4). In particular, note that the morphospaces through geological time, from the Middle Jurassic onwards, are not much smaller than the morphospaces occupied by the represented sample among 10,000-strong extant squamates. Therefore, we have either identified more or less the correct extent of morphospace for the Middle to Late Jurassic, despite the poor fossil record at that time, or that with much richer finds from that time interval occupied morphospace was even larger than we identify here. This would then enhance our interpretation of an early burst in squamate morphology and function. We posit that our conclusions regarding morphospace expansion are not affected by sampling; we predict that new fossil finds in the future will mostly fit inside the demarcated area of occupied morphospace.

Our interpretation here of the long-term stable morphospace occupation by squamates is compelling because the apparent increase in species richness through time, even if partially influenced by bias in the fossil record (Cleary et al., 2018), is not linked to a great increase in occupied morphospace. The observed expansion in occupied morphospace and rapid evolutionary rates coincides not only with the first presence of squamates in the fossil record, but also with the time when many crown groups are first recorded (e.g., scincoids, anguimorphs, and likely gekkotans and snakes), in the Middle–Late Jurassic (Evans, 2003; Caldwell et al., 2015; Estes, 1983; Evans, 2008; Figure 1). The Jurassic expansion of squamates is further supported by (1) the fact that all three main squamate morphological groups in morphospace (clusters 2–4) are already present in the Middle–Late Jurassic bin (Figure 2); (2) bootstrapped and rarefied measures of disparity through time (Figure 4A) present a peak roughly corresponding to the Late Jurassic for all lepidosaurs and for squamates alone; and (3) the evolutionary rates calculations also show peak in the Late Jurassic (Figure 4B).

This explosive adaptive radiation of squamates in the Middle–Late Jurassic situates the dates of origin of major clades in line with current phylogenomic analyses (Mulcahy et al., 2012; Zheng and Wiens, 2016; Pyron, 2017; Burbrink et al., 2020), mostly into the Jurassic, except for some groups in the Hedman trees, the divergence ages of which are even older, placed in the Triassic. It is worth noting that other key tetrapod groups also radiated ecologically in the Jurassic, namely, illustrated by the rapid diversification of paravians, the clade including birds and related small, feathered theropods with elongate wing-like arms (Lee et al., 2014; Puttick et al., 2014; Brusatte et al., 2015) and the expansion of early mammalian clades (Close et al., 2015). This predates the second ecological expansion of these three major clades, accounting for more than 95% of the modern biodiversity of tetrapods, which happened in the mid-Cretaceous in association with the KTR (Lloyd et al., 2008; Doyle, 2008; Benton, 2010; Meredith et al., 2011; Cardinal and Danforth, 2013), when diversity, and probably also abundance, exploded in line with the new food resources on land. Squamates remained at low diversity through the Triassic (where unambiguous fossils are yet to be recovered, but supposed to be present), Jurassic, and Early Cretaceous (Jones et al., 2013; Caldwell et al., 2015; Evans, 2008), and species richness seems to have risen massively during the KTR some 100 Ma, but the morphological expansion had already happened some 60 Myr earlier, in the Middle to Late Jurassic. This seems to fit an already identified pattern where the main diapsid groups present a long chronological lag between the initial phenotypic radiation of the group and its subsequent taxonomic diversification (Simões et al., 2020; Close et al., 2019) or, alternatively, it is related to a failure of an impoverished fossil record to reveal the true diversity achieved in pre-Late Cretaceous times.

What was happening in the Middle Jurassic that could have triggered squamate morphological diversity? (1) The supercontinent Pangaea began to split into precursors of the modern continents; (2) temperatures rose sharply for a short time; (3) gymnosperm plants diversified; and (4) various insect groups (e.g., mayflies, crickets, cockroaches, bugs, cicadas) diversified. All these factors may have had a role in driving some aspects of the early burst of squamate disparity, and they all require further investigation. This early radiation of squamates had been previously inferred (Evans, 1998; Evans, 2003) on the basis of a crude interpretation of the fossil record and tree topologies. However, it is the first time this issue has been approached with quantitative methods involving such an array of diverse points of view (phylogeny, dating, fossil record, morphospace, disparity, and evolutionary rates). Moreover, most of the Jurassic forms are very difficult to classify, and many of them have been reinterpreted since this was proposed (e.g., Marmoretta, see Griffiths et al., 2021). Our methods, however, do not necessarily rely on the achieved identification of each fossil because they feed on the morphological information stored in the character matrix, and not the specific topology derived from its analysis. Thus, the method used accounts for possible shifts in the phylogenetic position recovered for each form.

In their later evolution in the Mesozoic, all living clades of squamates diversified rapidly through the KTR. In addition, new and short-lived squamate groups arose in the Late Cretaceous, such as the terrestrial borioteiioids and the marine mosasaurs and relatives, but they disappeared, together with non-avian dinosaurs and other groups of diapsids at the end of the Cretaceous. The non-survival of such groups emphasizes the importance of the origin of the key modern clades in the Middle Jurassic and the establishment of their key ecomorphological adaptations – these then proved robust to various crises, including the end-Cretaceous mass extinction. Our integrative study here, incorporating current phylogenomic analyses of relationships of squamate clades with current fossil data, and novel computational methods in disparity and evolutionary rates, provides a synthetic narrative of the origin of one-third of modern tetrapod biodiversity, the Squamata.

Although morphological matrices might not be ideal for macroevolutionary inferences because they were specifically built for inferring phylogenetic relationships, they are handy in that they represent readily available sources of information, and they allow the mixture of taxa for which ecology is known (extant taxa) and taxa for which it can only be inferred (fossils). Moreover, results actually show that some ecomorphological signal is present in such datasets. Although the recovery of a limbless cluster of taxa might seem trivial, in fact it shows that the ecomorphological signal is overprinting the phylogenetic signal in that case because otherwise snakes and amphisbaenians would cluster with their respective closer clades (anguimorphs and iguanians for the former, lacertids for the latter). There are many other examples of this convergence in morphospace that can be interpreted as related to ecological niche convergence, for example, xantusiids and gekkotans, two groups that are not closely related phylogenetically but greatly overlap in our morphospace. The tight clustering of multiple groups close to the centroid does not help in interpretation, but the overall morphospace distribution shares many similarities with the niche plots reported by reference (Pianka et al., 2017) according to extant taxa scored for five niche dimensions. Their Figure 4 perfectly shows that a mixture of niche conservatism (phylogenetically close taxa tend to occupy similar niches) and niche convergence (distantly related species with similar ecomorphology tend to cluster together) occurs. Although represented taxa are not completely comparable (Pianka et al.’s dataset includes only lizards, lacking snakes and amphisbaenians among squamates, and also rhynchocephalians), and even the lizards sampled are different at the genus or species levels, similarities between our morphospace plot and their niche plot for extant groups include (1) anguimorphans are in both cases the most differentiated group, far from a much more populated cluster of taxa around the centroid that includes most of the rest of lizards; (2) this centered cluster includes scincoids, lacertoids, gekkotans, and some iguanians in our plot, whereas in Pianka’s plot gekkotans and most iguanians overlap outside this cluster, in the opposite direction of anguimorphans along PC1; and (3) small teiids overlap lacertids in both cases, but large teiids (Dracaena and Tupinambis in our case, Tupinambis in the Pianka et al.’s plot) are closer to anguimorphs.

Macroevolutionary studies can be strongly influenced by an array of potential biases that sometimes compromise results to variable degrees. Several potential issues have been identified through the design, development, and review of this study, ranging from sampling to methodological and interpretative factors. Moreover, methods are quickly evolving and can be quickly displaced by more refined approaches or criticized in their use or misuse. We have made an effort to consider as many variables as possible by assessing multiple potential resolutions for the phylogenies (constrained vs. unconstrained), specific changes in topology (by randomly selecting multiple MPTs), dating (by using three different methods, and randomly dating each fossil tip multiple times, accounting for geological range uncertainty), and by using multiple metrics and time bins when necessary. Other factors, like the possibility that our results are biased by the nature of the fossil record and how it conditions effective sampling across different time bins, are difficult to circumvent. We think, however, that if the poor fossil record is affecting results, it is most probably undermining the effect of the JME because (1) we have not been able to include any of the known Middle Jurassic squamate fossils and (2) the Late Jurassic contains a low number of samples compared to the Early and Late Cretaceous.

Among other studies that have emphasized the importance of the KTR in squamate evolution is Lafuma et al., 2021. In a study of origins and losses of tooth complexity across the clade, they found that tooth complexity first increased in the Late Jurassic, although it is regarded as marginal until the KTR. This increase in tooth complexity is apparent when its distribution through the Jurassic and Cretaceous is analyzed, but the possibility that the change in the quality of the fossil record might be enhancing the much more complete sample occurring in the Late Cretaceous is not discussed. Further, lumping the diverse morphologies of unicuspid teeth into a single category is potentially problematic if carnivores and insectivores are to be considered as distinct styles of predators. Another interesting but ignored result of that study is that, besides presenting a Cretaceous turnover (speciation/extinction) peak coinciding with the KTR, there is a previous peak mostly coinciding with the Jurassic–Cretaceous boundary. In any case, it seems that there is a shared pattern to our results, where innovations are initially explored in the Jurassic and then fully exploited in the Cretaceous, coinciding with the KTR.

To the uncertainty generated by the incomplete fossils that can only be scored for a minor portion of the morphological characters, we face the added problem of unscorable characters. In a morphological matrix of characters for lepidosaurs, this is not a minor issue because there is a long list of characters that cannot be scored for multiple groups, for example, characters related to limbs in limbless taxa, or characters related to structures only found in the snake skull. Even though the use of inapplicable characters has been discouraged (Gerber and Ruta, 2019), we think that simply removing them from the analyses is not the best solution, just as it would not be for a phylogenetic analysis.

Finally, results that directly depend on the estimation of time-calibrated branch lengths, such as the calculation of evolutionary rates, should be treated as preliminary because they must be validated under the use of more robust methods of time calibration, such as those that incorporate molecular data alongside the fossil record, as well as employ more realistic models of diversification such as the fossilized birth–death prior. In this sense, an ongoing study (work in progress) aims to analyze the present matrix and other datasets by using Bayesian tip dating under relaxed morphological clocks as described in Zhang and Wang, 2019. This allows us to calculate phylogeny and estimate divergence times and evolutionary rates while accounting for their uncertainties, and allow the use of both morphological and combined (morphological plus molecular data) matrices. This additional study should help clarify if the signal recovered in this work is reliable, or on the contrary it is biased by the chosen methodology. Meanwhile, the results presented here question the alternative view that regards the great diversification of squamates as occurring in the mid-Cretaceous, coinciding with the KTR. The first half of the Mesozoic has a great potential for unveiling the key milestones in the evolutionary history of lepidosauromorphs in general, but also of squamates in particular. Current reanalyses of classic material and the description of new specimens and taxa are already displacing the focus from the Late Cretaceous to the first half of the Mesozoic and are expected to provide insights on the issue presented here.

Materials and methods

Taxa and character data

The data source for all morphological character and taxon data analyses is the morphological data matrix of Conrad, 2018, reduced in our study to 201 species of living and fossil lepidosaurs scored for 836 skeletal morphological characters. We used this data matrix because it is by far the most extensive in terms of taxa and characters.

Phylogeny and timescaling

Phylogenetic analyses were performed in TNT 1.5 (Goloboff and Catalano, 2016). The settings for the unconstrained analysis are the same as in the original publication (Conrad, 2018) (ratchet and drift options activated, except that we set analyses to 100 replicates instead of 200). An alternative version of the phylogeny was obtained after constraining the general relationships recovered in molecular studies for those groups that present discrepancies in their position in morphological analyses, among others the sister group relationship of Iguania to the rest of crown squamates and the grouping of limb-reduced and limbless forms in the called ‘fossorial’ group (including dibamids, snakes, amphisbaenians, and limbless skinks), which is the result of convergences and clearly do not form a monophyletic group. For this, we randomly chose one of the MPTs recovered in the first analysis and forced the topology of phylogenomic studies for extant clades by defining the monophyly of the main extant groups according to Pyron et al., 2013. We set up fossils, which account for more than half of the taxa comprising the matrix, as floaters, so they could freely move around the tree. In both cases (constrained and unconstrained analysis), the resulting MPTs were exported to PAUP (Howard et al., 2002), where consensus trees were calculated. We produced the time trees for illustration (Figure 1, Figure 1—figure supplement 1) and rates calculations using fossil data to date origins of clades and time calibrated the trees in Paleotree v. 3.3.0 (Bapst, 2012) and using the Hedman method (see below).

Morphological disparity

All disparity and macroevolutionary analyses were performed in R (R Development Core Team, 2013). For disparity analyses, the pipeline started with the calculation of a pairwise morphological distance from the original character data using the package Claddis and maximum observable rescaled distances (MORD; Lloyd, 2016). The pairwise distances data was then subject to principal coordinates analysis (PCO) to identify the major axes of morphological variation. The resulting ordination matrix was used to plot morphospace based on PCOs 1–3. This morphospace was combined with a single topology (dated using the same method) to illustrate phylomorphospace and a chronophylomorphospace. We also plotted morphospace occupation in temporal bins. Finally, we used both pre-ordination (weighted mean pairwise distance, WMPD) and post-ordination (SoV, calculated in DispRity, Guillerme and Poisot, 2018) metrics to calculate global disparity, disparity in specific groups, as well as disparity through time. We also calculated completeness and sampling across the different time bins for comparisons with disparity results. We used various packages in R for plotting, namely, Plotly (Sievert, 2018), ggplot (Wickham, 2016), Geomorph (Adams et al., 2019), Claddis (Lloyd, 2016), and Phytools (Revell, 2012).

Morphological evolutionary rates

Rates of morphological evolution were analyzed using maximum-likelihood methods applied to the discrete skeletal character dataset and a range of phylogenetic trees. We used the DiscreteCharacterRate function from the R package Claddis and ran calculations for five of the unconstrained MPTs and five of the constrained MPTs, separately. We used a modified version of the code from Moon and Stubbs, 2020. The methodology first seeks to identify rate heterogeneity across the whole tree and then highlights branches or temporal bins with significant rate deviations (notably fast or slow) using likelihood ratio tests (Lloyd, 2016). To ensure rate results are consistent, the different topologies were dated multiple times (in our case, five dating replicates for each of the five randomly selected trees, for both unconstrained and constrained trees). We also repeated this for three dating methodologies, using the ‘equal’ method (Brusatte et al., 2008), ‘minimum branch length’ approach (Laurin, 2004), using the R functions from Lloyd et al., 2016 and a whole-tree extension of the Bayesian Hedman algorithm (Hedman, 2010). The Hedman node-dating approach uses Bayesian statistics, incorporating probability distribution constraints based on successive outgroup taxa ages (Hedman, 2010). We calculated per-bin evolutionary rates in two sets of time bins, one corresponding to geological stages and one corresponding to equal 10 Myr bins. To illustrate the rates results, we use ‘spaghetti plots’ showing individual lines for each combination of tree and dating (25 individual lines), as well as an average line, and also highlighting iterations and bins with significantly fast and slow evolutionary rates (Figure 4—figure supplements 5–22). In the main Figure 4B, we present summaries of these analyses.

Plots of Cramér coefficients

We used Cramér coefficients (Appendix 1—figure 7, Appendix 2—figure 2) to show correspondence between characters and PCO axes (Kotrc and Knoll, 2015; Nordén et al., 2018). See Appendix 4 for more details.

R scripts are available as Source code 1 (for morphospace and disparity), Source code 2 (evolutionary rates), Source code 3 (morphospace clusters), and Source code 4 for alternative analysis without integument and Cramér values, and Source code 5 for plotting clade-colored consensus trees.

Acknowledgements

We thank Graeme Lloyd, Benjamin Moon, Armin Elsler, and Guillermo Navalón for help with the scripts and discussion. Funding for AB comes from a Newton International Fellowship (NF170464) and a Juan de la Cierva Incorporación fellowship (IJC2018-037685-I) of the Spanish Government, and by the CERCA Programme of the Generalitat de Catalunya. This is part of the ERC Innovation Advanced Grant to MJB. (ERC 788203). We thank Min Zhu, George Perry, and two anonymous reviewers for comments on an early version of this manuscript.

Appendix 1

Appendix 1—figure 7. — Characters do not appear in the original order because they have been grouped by anatomical regions. Inside each region, characters on the left are significant (according to their p-values, in green when significant, in red when not significant).

Appendix 2

Appendix 3

Appendix 4

Extended methods

The morphological matrix (Conrad, 2018) used for phylogenetic analyses included in this work feeds upon previously published matrices for osteological characters (Conrad, 2008; Conrad et al., 2012; Smith, 2009; Gauthier et al., 2012), soft tissue (Schwenk, 1988), and salivary compounds (Fry et al., 2006). A few preliminary runs of disparity and evolutionary rates in alternative matrices (e.g. Simões et al., 2018) revealed that the results (not shown) were compromised by the poor sampling of squamates at key time bins (e.g., Late Jurassic), and we abandoned their use. Moreover, during the process of review of our article, Simões et al., 2020 published their own study on evolutionary rates and disparity in reptiles based on their own matrix (2018), and an additional work (Martínez et al., 2021) deals with the description a new stem lepidosaur and its position in morphospace. Results of those studies are compared to our own in the results section below. Regarding the final version of the matrix (Conrad, 2018) used here, it should be noted that running the original file as it was published does not yield the results reported in the corresponding paper. One problem relates to the comment in the text (p. 618, line 6) regarding deactivating characters 236, 242, and 364. This should not be done because these characters were evidently removed from the final matrix, and this is the reason why it has only 836 characters instead of the 839 stated in the text. We also decided to delete a few taxa for various reasons. Some terminal taxa correspond to additional specimens of a same taxon (e.g., Slavoia, Globaura, or Eoxanta). Others, like the fragmentary Tikiguania (Datta and Ray, 2006), are no longer relevant since it has been reinterpreted as a non-Triassic (probably Quaternary) agamid lizard (Hutchinson et al., 2012), and in any case it is not particularly complete or informative. Regarding Ardeosaurus, the matrix contains more forms than those reported in the original study (Conrad, 2018), so we deleted the taxon named ArdeosauruscfBRE. Finally, we deleted the taxon labeled ChineseScincoid because it is not clear which specimen is intended. A second problem concerns the list of ordered characters. Apparently, this list remained unchanged in the published matrix after the deletion of characters 236, 242, and 364. Accordingly, ordering of characters from 236 onwards needs to be changed, and all these modifications have been applied to our published version of the matrix, so this file (Nexus file in Source code 1) is ready to be used. Finally, we suggest that some aspects of the matrix may benefit from a detailed revision (e.g., the scorings of some characters, character ordering, addition of taxa belonging to poorly sampled taxonomical groups) if phylogenetic results are the aim of future works.

Several different time bin schemes are used throughout the analyses. Two of them use stratigraphic bins corresponding to the Triassic, Early Jurassic, Middle–Late Jurassic, Early Cretaceous, Late Cretaceous, and Paleogene-present (time scheme 1) or to the Triassic, Jurassic, Early Cretaceous, Late Cretaceous, Paleogene–Pliocene and Pleistocene-present (time scheme 2). Other time schemes correspond to bins with a specified constant length (e.g., 10 Myr or 24 Myr).

Analyses of disparity and morphospace occupation (see Source code 1 for the code used and source files) are based on the same matrix as the phylogeny. However, we removed Bharatagama rebbanensis and Adriosaurus microbrachis from those analyses because their inclusion resulted in incalculable distances. We used the MORD metric because of reported problems (Conrad, 2008) with the Geometric Euclidean Distance (GED) in matrices with a high percentage of missing information. We used the uncorrected option for the cmdscale function in Claddis, but we also checked that results were comparable using the corrected version (results not shown). We report the stacked temporal plots for PCO1–PCO6 and both stratigraphical schemes but note that the rest of the axes can be easily plotted upon minimal editing of the code. A single randomly selected topology (among the MPTs resulting from the constrained phylogenetic analysis) dated using all three methods (‘Hedman,’ ‘equal,’ and ‘MBL’) was selected for plotting phylomorphospaces using the function plotGMPhyloMorphoSpace of the package Geomorph (Adams et al., 2019). We used colors for different groups in separate plots: (1) major lepidosaur groups (stem lepidosaurs, rhynchocephalians, and squamates); (2) more restricted groups (e.g., iguanians, anguimorphs, etc.); and (3) toxicoferan vs. non-toxicoferan squamates. In order to be able to plot the chronophylomorphospace, we re-calculated principal coordinates using the function MorphMatrix2PCoA in Claddis applying the following options: distance.method = "MORD", transform.proportional.distances = "arcsine_sqrt", correction = "none", estimate.allchars = FALSE, and estimate.tips = FALSE. Then we edited the Claddis package function ChronoPhylomMorphospace so that plot parameters can be easily changed. The chrono.subsets function of DispRity (Guillerme and Poisot, 2018) was used to split morphospace along two different time bin schemes. The first scheme contains equally long (24 Myr) time bins, whereas the second one corresponds to the stratigraphic time scheme 1 described above. The function boot.matrix was used to bootstrap the data 100 times, and rarefy to n = 6. The SoV was calculated for the two different time bin schemes with the function DispRity. For both time schemes, we plotted the disparity of all taxa and of squamates alone. WMPD was calculated in Claddis (Lloyd, 2016) and bootstrapped 100 times. Analyses were performed for all taxa and for squamates alone. The clusters used in Figure 3D are based on the result of applying the pamk R function of the package fpc, as plotted in Figure 3—figure supplement 6 (data for the first two axes used).

In order to assess changes in sampling across different time bins, we plotted both the sample number (number of taxa in each time bin) and completeness (percentage of characters recorded for each taxon, grouped by time bin). This was applied to the two stratigraphic time schemes (time scheme 1 and 2). Note that for those time bins with one or no taxa sampled completeness and disparity measures are not calculated.

Evolutionary rate analyses were run for three subsets of the morphological matrix: (1) all taxa, (2) all squamates, and (3) all rhynchocephalians and stem lepidosaurs (see Source code 2 for the code used and source files). We ran analyses twice, one for five randomly selected MPTs of the constrained analysis, and the other for five randomly selected MPTs with the unconstrained analysis as the source. We dated each of these five MPTs five times with three different methods: ‘equal’ method (Brusatte et al., 2008), ‘minimum branch length’ approach (Laurin, 2004) using the R functions from Lloyd, 2016, and a whole-tree extension of the Bayesian Hedman algorithm (Hedman, 2010). The Hedman node-dating approach uses Bayesian statistics, incorporating probability distribution constraints based on successive outgroup taxa ages (Hedman, 2010). Successive outgroup taxa that are both more ‘basal’ and predate the first lepidosaurs were required to date the nodes close to, and including, the root. Occurrence dates of the following outgroup taxa were used: Weigeltisaurus jaekeli, Eunotosaurus africanus, Lanthanolania ivakhnenkoi, Orovenator mayorum, Petrolacosaurus kansensis, Anthracodromeus longipes, Hylonomus lyelli, Palaeomolgophis scoticus, Casineria kiddi, Ossirarus kierani, Tulerpeton curtum, Ymeria denticulata, and Ichthyostega stensioi (ages are provided in the supplementary code). For illustrating the results of the per-bin evolutionary rates (according to the two time sets, one corresponding to geological stages and one corresponding to equal 10 Myr bins), we used ‘spaghetti plots’ showing individual lines for each combination of tree and dating (25 individual lines), as well as an average line, and also highlighting iterations and bins with significantly fast (red triangle) and slow (blue rhomb symbol) evolutionary rates.

R scripts are available as Source code 1 (for morphospace and disparity), Source code 2 (evolutionary rates), Source code 3 (morphospace clusters), Source code 4 for alternative analysis without integument and Cramér values and Source code 5 for plotting clade-colored consensus trees.

Appendix 5

Extended phylogenetic analyses results

Results of phylogenetic analyses are explained here in detail for two main reasons: (1) the constrained analysis is completely new, as such constraints have not been applied to analyses of this matrix before; (2) a detailed explanation of the results of the unconstrained analysis was not provided by Conrad, 2018, who did not figure complete consensus trees because the subject of interest for the article was the position and interrelationships of ardeosaurs and related forms; (3) for both analyses the composition of our sampled taxa is different from the original because we removed a few taxa from the matrix; and (4) we changed the ordering of a few characters (see comment above).

Constrained analysis

Note that the phylogenetic relationships between extant clades were forced according to phylogenomic results, prior to this analysis. In contrast, the interrelationships between taxa forming these clades and, more importantly, the position of all fossil taxa, which were treated as floaters and could thus move to any possible position in the tree (Figure 1—figure supplements 1–4), are results of the present analysis. The strict consensus tree (Figure 1—figure supplement 4) has poor resolution because of the poorly defined position of a few rogue taxa, so we selected the majority rule tree (Figure 1—figure supplement 2) to discuss the results of the phylogeny. Note, however, that five randomly selected MPTs were used for evolutionary rate analyses, being one of them randomly selected for constructing phylomorphospace and chronophylomorphospace.

Outgroup and non-lepidosaur lepidosauromorphs

We followed Conrad, 2018 in considering Pamelina polonica as the outgroup. The list of stem lepidosaurs is completed with Icarosaurus siefkeri and Kuehneosaurus latus (sister taxa, considered kuehneosaurs in most works) and Marmoretta oxoniensis, which is sister to crown lepidosaurs.

Rhynchocephalians

Rhynchocephalians form a well-supported monophyletic group containing all forms typically regarded as such but, unexpectedly, with Scandensia ciervensis as their sister taxon. The latter had been previously regarded as a stem squamate (e.g., Evans and Barbadillo, 1998), as a stem ‘scleroglossan’ (e.g., Conrad, 2008; Bolet and Evans, 2011) or in a more derived position (Bolet and Evans, 2011), but never as closely related to rhynchocephalians. Note, however, that we have noticed that Conrad, 2018 scored the dentition of Scandensia as triangular (character 212, state 1), when it would be best coded as straight and pointed (state 0). This, together with the persistent notochord (character 230, state 0), a character state that is most probably convergent with rhynchocephalians, and a few other characters, is possibly behind this unexpected (and presumably unreliable) result. Accordingly, we regard the previously reported squamate affinities of Scandensia as unquestioned. The few MPTs that do not recover Scandensia as sister to rhynchocephalians place it in a more typical position among stem squamates. Our results also differ from specific studies on rhynchocephalians (e.g., Herrera-Flores et al., 2017) in that Gephyrosaurus, Planocephalosaurus, and Diphydontosaurus usually (most MPTs) form a monophyletic group, instead of forming a paraphyletic assemblage on the stem of ‘derived rhynchocephalians.’ Among the latter, two main monophyletic groups are recovered: the first contains Kallimodon (a sapheosaur), Priosphenodon (an eilenodontine), Pamizinsaurus, Palaeopleurosaurus, Pleurosaurus (pleurosaurids), and Bharatagama; the second one contains Sphenodon and Cynosphenodon (sphenodontids) and Ankylosphenodon, Polysphenodon, and a paraphyletic Clevosaurus (clevosaurs).

Stem squamates

Bavarisaurus macrodactylus, Huehuecuetzpalli mixtecus (forming a monophyletic group), and Hoyalacerta sanzi are recovered as stem squamates in a high percentage (97%) of the MPTs. B. macrodactylus had previously been recovered as a stem scincogekkonomorph in the original analysis of this matrix (Conrad, 2018), but all three taxa have been previously suggested to be stem squamates (Reynoso, 1998; Evans and Barbadillo, 1999; Evans et al., 2006).

Dibamidae

We forced the position of dibamids as the sister to the rest of crown squamates, but it is worth mentioning that none of the fossils included in the analysis is recovered as related to dibamids in the majority rule consensus tree, rendering the branch leading to them as one of the longest in the tree (together with that of the sole extant rhynchocephalian Sphenodon). A minor portion of the MPTs, however, recovered the Cretaceous taxa Sineoamphisbaena and Polrussia on the stem of Dibamidae.

Gekkota

Norellius nyctisaurops and Gobekko cretacicus are recovered as stem gekkotans, but Hoburogekko suchanovi forms part of the crown in our analysis because the extant Pygopus lepidopus and eublepharids are recovered as consecutive sister taxa of the remaining crown gekkotans (including Hoburogekko).

Scincoidea

Paramacellodus oweni and the paramacellodid from Utah (but not other paramacellodids in the classic view, see below) are recovered as stem scincoids. Tepexisaurus tepexii is recovered as a stem xantusiid, coinciding with other studies (Gauthier et al., 2012). This genus was originally considered a stem scincoid (Reynoso and Callison, 2000), but note that in the corresponding phylogenetic analysis xantusiids were recovered as lacertoids rather than scincoids. Our results are coincident with most previous studies (e.g., Gauthier et al., 2012) in that Palaeoxantusia is a crown xantusiid. Pseudosaurillus becklesi is nested within globaurids, a Campanian radiation on the stem of Scincidae plus Cordyliformes and containing Globaura, Bainguis, Eoxanta, Parmeosaurus, Hymenosaurus, Myrmecodaptria, and Slavoia in our results. Palaeolacerta is recovered (in just 68% of the trees) as sister to cordyliformes, but its position is rather unstable as it is recovered in several distant positions.

Lacertoidea

None of the fossil taxa included in our analysis is recovered as a stem lacertoid. Dracaenosaurus is confirmed as a crown lacertid (e.g., Čerňanský et al., 2017). Besides the position of Liushusaurus as a stem teioid (an unusual position for this taxon), the fossil record of this group in our analysis is limited to fossils of species reaching the present (e.g., Tupinambis teguixin).

Serpentes

An unexpected result of the present analysis is the position of Jucaraseps, Sineoamphisbaena, and Cryptolacerta in the stem of Serpentes in most trees. None of these has ever been claimed to be related to snakes. The phylogenetic relationships of Jucaraseps are not clear, but it has been regarded as a ‘scincogekkonomorphan’ (Bolet and Evans, 2012). The position of Sineoamphisbaena has been problematic as it was described as related to amphisbaenians and later reinterpreted as related to macrocephalosaurs (e.g., Kearney, 2003). Cryptolacerta presents a similar case: it was described as a stem amphisbaenian (Müller et al., 2011), but it has been regarded as more closely related to lacertids (e.g., Longrich et al., 2015). Portugalophis lignites and Parviraptor estesi are recovered as Jurassic crown snakes (Caldwell et al., 2015) in all trees. Dinilysia patagonica and Najash rionegrina, on one side, and Haasiophis terrasanctus and Pachyrhachis problematicus, on the other, are two separate lineages in our results.

Anguimorpha

Different versions of Conrad’s matrix (Conrad, 2008; Conrad, 2018; Conrad et al., 2011) are consistent in recovering as anguimorphs a few taxa that had been rarely referred to this group of lizards. In our majority rule tree, these include the clade formed by Meyasaurus diazromerali, Eolacerta robusta, Ornatocephalus metzleri, Yabeinosaurus tenuis, and Becklesius hoffstetteri. The latter has been usually regarded as a paramacellodid (e.g., Estes, 1983). Regarding the crown, Xenosauridae includes the extant Xenosaurus grandis as well as the fossil Exostinus serratus and Restes rugosus. The next clade to diverge is Shinisauridae, including the extant Shinisaurus crocodilurus and the fossil taxa Bahndwivici ammoskius, Merkurosaurus ornatus, and, less expectedly, Dalinghosaurus longidigitus and Parasaniwa wyomingensis. Anguidae are represented by the extant Elgaria coerulea and the fossil glyptosaur Peltosaurus granulosus. The rest of the anguimorphs form a monophyletic group consisting of Dorsetisaurus purbeckensis, Eosaniwa koehni, and a clade formed by mosasaurs and more advanced ‘varanoids,’ replicating the typical morphological structure of the tree. This is because although a sister group relationship between iguanians and anguimorphs was forced, the molecular topology of the constituent clades of anguimorphs was not. Mosasaurs include Eonatator sternbergi and Aigialosaurus with Paravaranus angustifrons as sister taxon of both, on one side, and Adriosaurus plus Proplatynotia longirostrata, on the other. ‘Varanoids’ consist of Parviderma inexacta, helodermatids (the extant Heloderma plus Saniwides mongoliensis, Palaeosaniwa canadensis, Estesia mongoliensis, Gobiderma pulchrum, Paraderma bogerti, and Primaderma nessovi), and varanids that, in our analysis, include the extant Lanthanotus borneensis, Varanus varius, and Psammosaurus griseus, plus the fossil Telmasaurus grangeri and a paraphyletic Saniwa.

Iguania

One of the most interesting results of the constrained analysis is the position of both ardeosaurs (sensu lato) and borioteiioids, on the stem of Iguania. This is most likely the result of forcing iguanians into a more crownward position (sister to anguimorphs), which somehow seems to drag both groups that were stem scleroglossans in the original analysis of Conrad, 2018 and our unconstrained analysis here with them. If this is correct, then ardeosaurs would be filling a gap between the earliest known anguimorphs and the earliest known iguanians. Ardeosaurs in our analysis consist of Ardeosaurus brevipes, Schoenesmahl dyspepsia, and Chometokadmon fitzingeri. A monophyletic Eichstaettisaurus is the next to diverge and is thus very closely related to the group above, but rendering Ardeosauridae paraphyletic if Eichstaettisaurus is included. Note that Conrad, 2018 considered Ardeosauridae, Eichstaettisauridae, and Bavarisauridae (as well as borioteiioids) as separate monophyletic clades on the stem of ‘Scleroglossa.

Borioteiioids form in our results a monophyletic group consisting of Polyglyphanodon sternbergi, Erdenetesaurus robinsonae, Tianyusaurus zhengi, Darchansaurus estesi, Gilmoreteius ferrugenous, Gilmoreteius chulsanensis, Adamisaurus magnidentatus, Cherminsaurus kozlowskii, Gobinatus arenosus, and Tchingisaurus multivagus. Among crown iguanians, the two main groups (the paraphyletic pleurodont iguanians and monophyletic Acrodonta) have radiations in the Campanian. On the one hand, Gobiguania includes the following pleurodont iguanians: Zapsosaurus sceliphros, Anchaurosaurus gilmorei, Ctenomastax parva, Temujinia ellisoni, and Saichangurvel davidsoni; on the other hand, the stem of Acrodonta is formed by priscagamids (Priscagama gobiensis, Phrynosomimus asper, Mimeosaurus crassus, Gladidenagama semiplena, and Flaviagama dzerzhinskii) and Arretosaurus ornatus.

Unconstrained analysis

Trees resulting from the unconstrained phylogenetic analysis (Figure 1—figure supplements 5–7) roughly match the results of the original analysis of this matrix (Conrad, 2018), mainly if we take into account that we removed a few taxa, and that a detailed comparison is not possible because this author did not publish a complete tree, but a strict consensus tree with most major clades collapsed instead. K. latus and I. siefkeri (kuehneosaurs) are sister taxa. However, in contrast to the results of the constrained analysis described above, M. oxoniensis is sister to rhynchocephalians instead of being sister to crown lepidosaurs. The phylogenetic relationships among rhynchocephalians are less well resolved and, besides this, the group formed by S. punctatus and Cynosphenodon huizachalensis is nested within non-clevosaur advanced rhynchocephalians, instead of grouping with clevosaurs.

Stem Squamata

Hoyalacerta sanzi and Huehuecuetzpalli mixtecus, but not Bavarisaurus macrodactylus (see constrained analysis above), are recovered on the stem of Squamata.

Iguania

The unconstrained version of the analysis recovers the typical position in morphological analyses of iguanians as sister to the rest of the squamates (e.g., Gauthier et al., 2012; contra Simões et al., 2018). Neither ardeosaurs nor borioteiioids are recovered on the stem of iguanians. Gobiguanians, with the same composition as in the constrained analysis, are here recovered on the stem of Iguania, instead of nested within crown iguanians. Crown iguanians are recovered again as containing a paraphyletic group of pleurodont iguanians and a monophyletic Acrodonta, the latter with Arretosaurus in its stem, and priscagamids as sister clade.

Stem ‘Scleroglossa’

The stem of ‘Scleroglossa’ is formed by borioteiioids and ardeosaurs (sensu lato), coinciding with the original results of Conrad, 2018. The latter are paraphyletic, formed by a first monophyletic group formed by S. dyspepsia, A. brevipes, and the unnamed genus corresponding to PMUR58; a second monophyletic group formed by B. macrodactylus, S. ciervensis, Eichstaettisaurus gouldi; and, finally, Eichstaettisaurus schroederi.

Gekkotans

Eoxanta lacertifrons and N. nyctisaurops are recovered as stem gekkotans. In contrast to the constrained analysis, Gobekko joins Hoburogekko as a Cretaceous crown gekkotan.

‘Autarchoglossa’

The group containing non-gekkotan ‘scleroglossans’ differs slightly from previous results in that teiioids are sister to a group formed by scincoids and anguimorphs. Scincoids include, besides the typical members of the group (scincids, cordyliforms, and xantusiids), globaurids (on the stem of scincidae + cordyliforms), dibamids, amphisbaenians, and snakes. The three latter form a convergent group inside Scincidae, which is why Conrad, 2008 coined the term Scincophidia.

Note: The unexpected and rather unlikely position of Jucaraseps, Sineoamphisbaena, and Cryptolacerta on the stem of Serpentes (Figure 1—figure supplement 2) in the constrained phylogenetic analysis seems to be incongruent with morphospace (see Figure 3—figure supplement 1) because these taxa cluster with lizards instead of snakes. This seems to hint at some problem with these taxa in particular that for some reason (possibly convergence in some cranial characters related to fossoriality or semi-fossoriality) are attracted to snakes in the phylogeny, but are plotted closer to lizards in morphospace.

Disparity and morphospace occupation

Plots of morphospaces with all taxa labels for PCO1–PCO6 are reported here (Figure 3—figure supplements 1–3), but plots for additional axes can be easily obtained by modifying the code to get the desired PCOs. Figure 3—figure supplement 4 shows hulls by taxonomic groups. A simplification of these groups into rhynchocephalians, generalized lizards, anguimorphs, and snakes and other limbless squamates in Figures 2 and 3b is based on the results of applying the pamk function (see resulting plot in Figure 3—figure supplement 6). Figure 3B and Figure 3—figure supplement 4 show that iguanians, and especially chameleons, are the squamate group that is closest to rhynchocephalians. Stem lepidosaurs are situated between rhynchocephalians and squamates, but closer to the latter. Besides iguanians, the cluster containing generalized lizards is formed by scincoids (scincids, cordyliforms, and xantusiids), gekkotans, lacertids, teiioids, and four fossil clades: ardeosaurs (and the possibly related eichstaettisaurs), borioteiioids, globaurids, and paramacellodids. A third cluster contains anguimorphs (including mosasaurs, which are contained in this clade in our phylogenetic results), and a fourth cluster consists of all clades of limbless taxa (snakes, amphisbaenians, and dibamids). Of these, snakes are situated furthest from the centroid.

We also provide the 2D phylomorphospace plot (Figure 3—figure supplement 5) that corresponds to the morphospace colored by squamates, rhynchocephalians, and lepidosaurs (Figure 3A), and a 3D phylomorphospace plot (Appendix 3—figure 2, Supplementary file 6) where toxicoferan vs. non-toxicoferan squamates are represented. This plot shows that overlap between both groups is present, but rather limited. In order to complement the temporal morphospace stack of Figure 2, we provide stacks for additional PCOs and time scheme 2 (Figure 2—figure supplements 4–6). The sudden increase of morphospace in the Middle–Late Jurassic and stability through time until the present day is confirmed in all these additional plots.

Regarding measures of disparity, the SoV (Figure 4A) presents two highs in the Late Jurassic and mid Cretaceous, and marked drops in the earliest Cretaceous and around the K-Pg boundary, a point from which disparity remains stable and intermediate between the low disparity of the Triassic and the high disparity of the Late Jurassic and mid Cretaceous peaks. The WMPD plotted by taxonomic group shows that squamates have a much higher disparity than either rhynchocephalians or stem lepidosaurs (Figure 4—figure supplement 1). For the WMPD by less inclusive taxonomical groups (Figure 4—figure supplement 2), if we arbitrarily set a high disparity for values above 0.4, low disparity for values below 0.35, and intermediate disparity for values between these values, we see that scincids, cordyliforms, snakes, mosasaurs, and anguimorphs show high disparity; stem lepidosaurs, xantusiids, lacertids, dibamids, borioteiioids, and ardeosaurs (including the possibly related eichstaettisaurs) show low disparity; and for rhynchocephalians, gekkotans, amphisbaenians, and iguanians disparity is intermediate. This same measure plotted through time (Figure 4—figure supplements 3 and 4) fails to recover the Middle–Late Jurassic peak on disparity, although disparity for this time bin is higher than those of the Triassic and Early Jurassic, similar to that of the Early Cretaceous, and only slightly lower than that of the Late Cretaceous. Note, however, that this measure has been bootstrapped but not rarefied, and is thus more sensitive to uneven sampling.

Evolutionary rates

Results for the stratigraphical time scheme (top plot in each figure of Figure 4—figure supplements 5–22, with bins corresponding to the Early, Middle, and Late Triassic; Early, Middle, and Late Jurassic; Early and Late Cretaceous; and Paleogene and Neogene) recover similar results as the 10 Myr long time bins (bottom plot of each figure), although, as expected, some resolution is lost. Results for the constrained phylogeny including all taxa and the ‘Hedman’ method (Figure 4—figure supplement 5) recover the highest peaks of evolutionary rates by the Late Jurassic and the Neogene-present time bins, with a moderate peak corresponding to the Late Cretaceous. The same analysis using the equal method (Figure 4—figure supplement 6) recovers the same two highest peaks in the Late Jurassic and Neogene-present, and two additional peaks, in the Late Triassic and Late Cretaceous, are only slightly lower. Using the ‘MBL’ method (Figure 4—figure supplement 7) instead produces an extremely high peak in the Neogene-present, possibly as the result of the placement of many nodes towards modern times according to the specific procedure of the method for establishing the age of a node. The other three peaks are present, but only in the Late Triassic and Late Cretaceous of the 10 Myr time bin scheme some trees show significantly high evolutionary rates (the rest being nonsignificant). Results for squamates alone using the ‘Hedman’ method (Figure 4—figure supplement 8) are very similar to those of all taxa, again with the highest peak in the Late Jurassic, a slightly lower peak in the Neogene-present, and an even lower (although still significant) peak in the Late Cretaceous. Roughly the same results are recovered when using the ‘equal’ method (Figure 4—figure supplement 9), but not when using the ‘MBL’ method (Figure 4—figure supplement 10), which again boosts the Neogene-present evolutionary rates to such a high level that the Late Jurassic and Late Cretaceous peaks are no longer significant. Rhynchocephalians (plus stem lepidosaurs) present high evolutionary rates mostly in the Late Triassic in all three dating methods (Figure 4—figure supplements 11–13), and in all cases there is a rather constant trend of decreasing evolutionary rates towards the present.

Results for all analyses of the unconstrained topology (Figure 4—figure supplements 14–22) are very similar to those for the constrained topology. The presence and prevalence (except for the MBL method) of the Late Jurassic peak in all iterations of the complete dataset and of squamates alone are confirmed. Slight differences include a less clear peak in the Late Cretaceous either because a lower number of trees present significantly high rates or because the peak moves to the Early Cretaceous (e.g., Figure 4—figure supplement 15, top and bottom, respectively). Regarding rhynchocephalians (plus stem lepidosaurs), results of the unconstrained analyses (Figure 4—figure supplements 19–22) are almost identical to those of the analyses of the constrained phylogeny, although this was expected because most changes in topology affecting divergence times (and thus branch lengths and related evolutionary rates) are concentrated in squamates.

Sampling, selection of axes, and Cramér values

Because the number of taxa and their relative completeness (number of scored characters against the total number of characters) may have an influence on disparity and evolutionary rates results, we plotted these statistics for the entire dataset in time schemes 1 and 2 (Appendix 1—figures 1–5). Time scheme 1 shows that the Middle–Late Jurassic has one of the lowest taxon counts (only slightly surpassing that of the Triassic and Early Jurassic) and the lowest completeness of all time bins (Appendix 1—figures 1, 3 and 5A). The Late Cretaceous is the time bin with greatest completeness if the one including the extant taxa (the Paleogene-present) time bin is excluded and presents also a high taxon count. For time scheme 2 (Appendix 1—figures 2, 4 and 5B), the Jurassic is retained as one of the two time bins with the lowest completeness. The Paleogene–Pliocene bin (used only in this second time scheme) records a decrease of completeness after the Late Cretaceous bin. In order to discard an important influence of characters with a high degree of missing scores (either because they are rarely preserved in fossils or because they represent unscorable characters in a good portion of the taxa sampled in the matrix), we performed an additional analysis where only those characters scored for more than 40% of taxa were included. This procedure removed an important number of characters, but resulted in a matrix containing less than 25% of total missing data, a threshold considered safe in terms of completeness. Interestingly, resulting plots of morphospace (not shown) are completely comparable to those reported for the complete matrix, suggesting that characters showing a high degree of incompleteness are barely contributing to the distribution of taxa in morphospace of the first axes.

It is worth noting that although results for all PCOs are not figured, they can easily be obtained applying small edits to the code. In any case, we show in Appendix 1—figure 6 that the variance explained by additional axes is very low in relation to the first ones.

Finally, we report Cramér coefficients in order to show the correlation between characters and the first two PCOs. Because the chosen methodology (calculation of PCoA) transforms the distribution of character states into a distance matrix, it is not possible to calculate PC loadings as is usually done with continuous data and PCA. Cramér coefficients have been used as an alternative to PC loadings (e.g., Kotrc and Knoll, 2015; Nordén et al., 2018), although most papers using PCoA skip this step entirely (e.g., Simões et al., 2020; Martínez et al., 2021). We generated a first plot (Appendix 1—figure 7) where we show all characters first grouped by anatomical regions (e.g., preorbital region, postorbital region, palate, braincase, etc.) and then, inside every region, characters that are significant are grouped on the left (green columns of the second and fourth plots, sometimes barely visible), whereas nonsignificant characters are situated to the right (red columns of the same plots). This is a good way to see all data, clearly showing that characters with a high Cramér value are concentrated in the integument/myology anatomical region. However, because most of these characters cannot be scored for fossils, we decided to recalculate the distance matrix, and then plot the new morphospace without including all those characters that were likely missing in all fossils (e.g., those related to skin or myology), and leaving only those that had some potential to be preserved in fossils (e.g., those related to osteoderms), besides those of osteology. The resulting morphospace (Appendix 2—figure 1) is very similar to that of the full dataset. This seems to be indicating that despite presenting the highest Cramér coefficients (meaning that they correlate well with the distribution in the corresponding axis), they are not essential to recover the morphospace we are discussing. This probably has to do with the high number of sampled characters, which dilute the weight of any given character. In order to focus discussion to a smaller number of characters, we plotted Cramér values for this second (reduced) dataset, which does not include these soft-tissue characters, only for those characters that presented a Cramér value >0.4 (Appendix 2—figure 2). Note, however, that characters with lower Cramér values are still contributing to the distribution, mainly for those taxa that cannot be scored for characters with higher Cramér values. Characters with a Cramér value >0.4 are predominantly concentrated in the vertebral column anatomical region (nine characters for both PCO1 and PCO2) corresponding to characters related to vertebral morphology. The rest of the regions are much less widely represented in PCO1, ranging from four characters related to the braincase to one related to the dentition or postorbital region, or even without representation like the hindlimbs. For PCO2, characters related to the vertebral column are equaled by characters related to the antorbital region (nine in both cases), being the rest of regions represented by a much smaller number of characters, from 0 (e.g., palate, forelimbs, dentition, or hindlimbs), 1 (postorbital, lower jaw), to 2 (e.g. integument – osteoderms or braincase). Note that, despite the fact that all these characters are contributing to the recovered distribution of taxa in morphospace, only those that present a p-value<0.05 can be regarded as significant (marked with an asterisk in the plot). These characters, which are more reliable for interpretation, include in the case of PCO1 the following list of characters: 260 (coracoid anterior emargination), 304 (dorsal compound osteoderms), 306 (ventral compound osteoderms), 420 (facial notch in the crista prootica), 422 (ventrolateral margin of the paroccipital process of the otoccipital), 736 (Hyoid cornu), and 791 (compound supraorbital scale osteoderms in the orbit). In the case of PCO2, the list is formed by the following characters: 49 (jugal, postorbital branch), 54 (quadratojugal), 204 (prearticular, crest with imbedded angular process), 237 (presacral vertebrae, length of transverse processes), 291 (egg teeth), 304 (presence/absence of dorsal compound osteoderms), 541 (frontal(s), subolfactory process fusion, and obliteration of midline suture), 672 (opisthotic/otoccipital, contribution to the posterior auditory foramen), and 800 (dorsal vertebra, midline inter-zygosphenoidal spur). Characters like, for instance, 49 or 54 are clearly contributing to the separation of rhynchocephalians and squamates along PCO2. Interpreting which characters are contributing to the separation of groups along PCO1 appears much more complicated and is not attempted further here. The conclusion would be that, although the characters related to the integument–myology are the ones with highest Cramér coefficients, they are not essential to produce a stable morphospace. Among the rest, it seems that the vertebral column first, and then regions like the braincase, dentition, and postorbital region would contain character states distributions that fit best with PCO1, whereas PCO2 would be more tightly related to the vertebral column again, and the antorbital region. However, as said above, the number of characters is so high that it is unlikely that a given character or set of characters is responsible to a great degree of the recovered morphospace distribution. Instead, the morphospace represents a culmination on multiple changing character scores, some with complex distributions. This is partially expected given that the dataset is composed of characters from across the whole skeleton, all of which change across the major morphotypes.

Appendix 6

Extended discussion

Phylogeny

Phylogenetic results for all main groups are reported above, but there are a few specific results that are worth discussing here. One is that forcing the molecular constraints results in a crownward movement of borioteiioids and ardeosaurs, which become stem iguanians. This is important in filling a gap between the earliest anguimorphs and the earliest iguanians. This placement also provides added evidence of the presence in the Jurassic of Toxicofera, a clade that is usually considered as highly derived, and suggesting that an important part of the evolutionary history of squamates occurred before the end of that period. The three main toxicoferan clades would be represented by Jurassic forms, namely, P. estesi and P. lignites as snakes (according to Caldwell et al., 2015, and results herein), Dorsetisaurus as an anguimorph, and ardeosaurs (sensu lato) as iguanians. Difficulties in the placement of borioteiioids, as reflected in their unstable position among different phylogenies in previous works, are potentially related to a possible convergent nature of similarities between teioids and toxicoferans. Note, however, that the exact phylogenetic position of problematic taxa (e.g., the identification of Parviraptor and Portugalophis as snakes) is irrelevant to morphospace and disparity discussions because topology is not considered in the construction of the morphospace, just used to illustrate phylomorphospaces.

Morphospace

Results of our disparity, morphospace occupation, and evolutionary rate analyses are discussed here in the context of all available evidence, including patterns of diversification as informed from the fossil record and current phylogenies. Phylogenetic results are not discussed in greater detail because they were only meant to provide a phylogenetic framework for evolutionary rate analyses (and a specific topology to be used in plots of phylomorphospace). Morphospace distribution of points and associated measures of disparity do not rely on phylogeny and, as expressed in the Materials and methods section, we have incorporated phylogenetic uncertainty into our analyses by conducting separate analyses for constrained and unconstrained phylogenies, and by using multiple randomly selected most parsimonious trees. Although the constrained phylogeny has been chosen to illustrate the results of our analyses, it is worth noting that our results are robust to changes in topology. A conflicting point regarding the evolutionary history of squamates is that the acceptance of the molecular topology requires that the numerous morphological similarities between iguanians and Sphenodon (e.g., Estes et al., 1988) are the product of reversals and convergences (Losos et al., 2012) and should expectedly change the timing of the events in the evolutionary history of the group. This similarity between iguanians and rhynchocephalians is evident in the 2D and 3D plots of morphospace (e.g., Figure 3B, Figure 3—figure supplement 4, Supplementary file 2), where they are plotted close to each other in morphospace. The 3D plot (Figure 3C) shows iguanians (especially chameleons) among all squamates as the closest group to rhynchocephalians. Note that striking similarities between rhynchocephalians and iguanians in soft tissue and osteology, such as an apparently conserved morphology of the tongue, the vomeronasal organ, and closely placed cranial bones (Conrad, 2008; Conrad et al., 2011), must be convergences if molecular studies are correct in placing iguanians with anguimorphs (Mongiardino Koch and Gauthier, 2018). Of note, removing rhynchocephalians and coloring toxicoferan vs. non-toxicoferan squamates (Appendix 3—figure 2, Supplementary file 6) results in a good separation of both groups, except in the contact region between them, where some overlap occurs. This suggests that some morphological support for Toxicofera is present in the dataset, even if unconstrained phylogenetic analyses fail to recover monophyly for the group.

The JME event

The results obtained here, implying the existence of a previously unidentified event triggering an early increase in disparity of squamates linked to high evolutionary rates by the Late Jurassic at the latest, have profound implications for interpreting evolutionary dynamics in lepidosaurs. Moreover, it not only changes the focus for understanding the main radiation of lizards (and lepidosaurs) from the Cretaceous to the Jurassic, but also adds to current discussions on the importance of Jurassic events in shaping Mesozoic tetrapod assemblages as a whole. Understanding the processes behind extremely successful clades is key not only to acquiring a more complete picture of past and present biodiversity, but also to help in the prediction of future trends for vulnerable portions of the tree of life. The study of squamates (lizards, amphisbaenians, and snakes), as one of the largest tetrapod clades that dominate modern landscapes yet one of the least understood, is not trivial in this regard. This limited knowledge on the early evolutionary history of the group is, in part, linked to a fossil record that is poor for a great part of its early history (Evans, 2003), as well as the fact that other clades of vertebrates have received greater attention for various reasons. This admittedly poor and uneven fossil record limits our understanding of the timing, mode, and reasons behind the diversification of the clade. As an example, a recent study (Cleary et al., 2018) on the diversity of lepidosaurs, based on generic occurrences along the Mesozoic and Paleogene fossil record, showed an apparent low diversity for the group for the greatest part of the Mesozoic, and then a sudden peak in the Late Cretaceous. However, the same study highlights that the available data for the greatest part of the Mesozoic is too poor to provide confident conclusions. This is in line with Evans, 2008, who stated that mid-Cretaceous squamates showed increased diversity, although they claimed that it was not possible to determine if this was a real Cretaceous trend or the result of a more complete record in the Cretaceous. At the same time, this impoverished early fossil record could be partly behind the delay in the first observed record of many clades and, according to this, in the calculated divergence ages and tied evolutionary rates of many key lineages.

From the Early Jurassic onwards, mesic Laurasian deposits are characterized by the presence of fish, lissamphibians, crocodiles, turtles, and choristoderes, in what Evans, 2008 describes as representing lowland, freshwater lagoonal, or wetland deposits. According to these authors, no Triassic/Early Jurassic deposit has yielded an equivalent assemblage, which instead usually contain rhynchocephalians. The fact that rhynchocephalians and squamates are recorded together in some Jurassic and Cretaceous deposits but always in unequal proportions (Evans, 1995) led some authors (e.g., Evans, 2008) to the conclusion that these differences in composition may be related to differences in ecology, rather than being purely taphonomic. We argue that the previously detected (but possibly artificial) increase in diversity related to the improved fossil record of the Late Cretaceous has overemphasized the importance of the events occurring in this period regarding the radiation of squamates. Although the KTR might be a real event that triggered diversification among different clades, it is not the first time that the fit of a clade radiation in this Late Cretaceous event is questioned. Just as an example, Lloyd et al., 2008 suggested that a sampling bias might be behind the apparent increase in diversification of dinosaurs through the Late Cretaceous, in a case similar to the one presented herein.

Another example showing that oversampling in a particular time bin might result in an exaggerated signal is the high rates recovered for the last 10 Myr time bin, an issue that sees its most exaggerated version in the evolutionary rates calculated using the MBL method. The latter method tends to place divergence times towards the present, mainly when the number of extant taxa is high, as a result of an effect equivalent to the pull of the recent described for diversity analyses (Raup, 1979; Sepkoski, 2016). This concentration of nodes in this time bin results in short branches that are necessarily correlated to generally higher evolutionary rates. This is because, if we assume a constant rate of changes, such a high concentration of short branches is expected to yield higher evolutionary rates. The opposite might be true for extremely long branches, which could be behind the low evolutionary rates recovered for extant taxa lacking closely related forms (e.g., the taxon Sphenodon in Herrera-Flores et al., 2017). The fact that the Late Jurassic peak in evolutionary rates appears in all the versions of the analysis performed, and in most of them being the highest peak, is interpreted here as proof of a true signal for high evolutionary rates because this is not a particularly highly sampled period (in terms of number of taxa, or of their completeness, see Appendix 1—figures 1, 3 and 5A).

Other considerations regarding our results concern groups that have a record that is much younger than their expected origin, or that have fragmentary and/or early records not included in the data matrix used (which is biased towards more complete specimens informative for phylogenetic analyses). For example, some Jurassic lizards have been interpreted as iguanians (e.g., Bharatagama, but see results in the original analysis of the matrix and herein), but the first clear record of the clade might be as young as the Late Cretaceous. If anguimorphs were present in the Late Jurassic (Dorsetisaurus), then iguanians, as their sister taxon, should be equally old. As expressed above, the recovery of ardeosaurs and borioteiioids as stem iguanians is interesting in this regard, as, if the position of the former is confirmed, it would fill a great portion of the gap between the first recorded anguimorphs and the first recorded iguanians. Other groups, like dibamids, have a much younger fossil record (if any) than would be expected from their sister group position to the remaining squamates. Other examples exist, but the point here is that a future discovery of stem members of clades with a poor fossil record is expected to result in even higher evolutionary rates (and possibly increased disparity) especially in early time bins. Current patterns of diversification and the apparently poorly sampled Middle and Late Jurassic, point to these intervals as the ones with a greatest potential for providing new evidence of an early radiation of squamates as new forms are discovered.

Another point that needs to be considered is that our analysis is biased towards more complete specimens. This is the result of the initial selection of taxa (Conrad, 2018), based on those scored for a maximum amount of codifiable characters in order to be useful in phylogenetic analysis. According to this selection procedure, the first sampled member of a group is not always the earliest taxon (or specimen) of that group to appear in the fossil record highlighting, again, that the morphospace occupation (and disparity) in Jurassic bins is probably underestimated. The possible influence of the ‘Lagerstätten effect’ in the current results is not investigated here, but deserves attention and will be assessed in future studies.

A second time bin scheme (scheme 2) was used to explore the possible effect of the K-Pg boundary. The results show a very slightly decreased morphospace occupation for the Paleogene–Pliocene time bin. However, because sampling and completeness for this time bin is smaller than for the previous time bin (Appendix 1—figures 2, 4 and 5B), we interpret this decrease as possibly an artifact of an impoverished fossil record and/or sample selection. This is related, in our particular case, to the limited number of well-preserved fossils of this age (Messel and the Green River formation being among the few exceptions), rather than a comprehensively poor fossil record. Even if we accepted the results as related to the true signal, occupied morphospace after the K-Pg boundary is not severely reduced as would have been expected if the K-Pg extinction had strongly affected squamates, as previously suggested (e.g., Longrich et al., 2012), but exploring the effect of this event would require shorter bins in order to track the post K-Pg event recovery, and much more complete Paleogene sampling. Even in terms of diversity alone, statements that lizards underwent a mass extinction at the K-Pg boundary are hard to sustain, considering that the number of clades that completely disappeared during this extinction is extremely low. Moreover, the number of extant clades that were already present in the Late Cretaceous (and thus, survived the extinction) is notable. Regarding the measures of disparity, the SoV and WMPD through time both show a decrease in disparity for the Paleogene–Pliocene time bin, although again we argue that this lower relative disparity might be related to a low sample size for this time bin, which hampers the recording of representatives of forms that would be necessary to get the true extent of morphospace and associated disparity values.

Evolutionary rates seem to provide stronger evidence for a decrease in evolutionary rates through the entire Paleogene. There is little diversification through this time period, which results in a higher proportion of long branches that, at the same time, can result in lower evolutionary rates. It is not clear, thus, if the Paleogene represented a period of poor diversification and low evolutionary rates or if, on the contrary, we are looking at an additional case of decrease in this measure related to poor sampling. Even if the Paleogene was interpreted as a time of slightly reduced morphospace and lower evolutionary rates, all seem to recover in the last time bin. This shows long-term stability for the morphospace occupied by lepidosaurs in general (and that of squamates in particular). Several specific features of the distribution of the different points in post K-Pg morphospace occupation in our results explain the apparent stability recorded: (1) the presence of the extant Sphenodon is important in maintaining the vertex corresponding to rhynchocephalians. It is not difficult to imagine an alternative scenario where this taxon (and thus all rhynchocephalians) had gone extinct, causing the removal of an important pole of lepidosaur morphospace. It has to be noted, however, that the large empty space between rhynchocephalians and squamates is gradually filled by iguanians (specially by chameleons) from the Cretaceous onwards, and thus the lack of rhynchocephalians would not result in such a large loss of occupied morphospace as might have been the case in earlier periods. (2) The clades that appear in post K-Pg assemblages do not expand the occupied morphospace. Instead, they fill empty spaces inside the envelope of already occupied morphospace or they directly overlap with clades already present.

Extended conclusions

According to our results, the age for crown Squamata would be Middle Jurassic (‘MBL’ method), or even earlier, in the Late Triassic (‘equal’ method), or the Middle Triassic (‘Hedman’ method). The age provided by other studies (e.g., Jones et al., 2013) is around the Triassic–Jurassic boundary, thus suggesting an intermediate age. In summary, the results provided by our disparity and evolutionary rate analyses strongly suggest that squamates underwent a great adaptive radiation by the Late Jurassic at the latest. Disparity analyses show a great change in occupied lepidosaur morphospace by the Middle–Late Jurassic, linked to the appearance of squamates in the fossil record. The Late Jurassic peak in lepidosaur evolutionary rates is present in all analyses, showing that these results are robust to changes in composition of the dataset and topology of the phylogeny. Importantly, recognized uncertainties regarding the phylogeny of squamates were considered by randomly sampling among different MPTs (uncertainties at a small scale, concerning a small number of labile taxa), but also confronting the pure morphological topology against the topology from constraining extant taxa to the molecular position. The Late Jurassic peak is recovered in all cases, suggesting that these results are robust to these changes in topology, and thus more reliable than if only supported by a few of the trees. Both the temporal range of each fossil and alternative methods in dating nodes have been considered. Because disparity results do not depend on a phylogeny (except when plotting a phylomorphospace), our results reveal a strong signal for this Middle–Late Jurassic squamate adaptive radiation, independent of which topology or dating for the trees is preferred. There has been a recent shift towards more integrative methods such as Bayesian clock estimates of divergence times and morphological evolutionary rates that are not used in this work because it was devised before these were widely used. It is worth noting that (1) our divergence times cover the ranges of divergence times recovered by other works that have applied Bayesian methods (e.g., Simões et al., 2020) and (2) our own preliminary results after applying the same Bayesian methods to the present matrix (work in progress) are largely congruent with the conclusions presented here.

It has been said (Gauthier et al., 2012) that rhynchocephalians dominated the fossil record of Lepidosauria during the Triassic and Jurassic, ‘with the squamate branch becoming abundant, in a classic pattern evolutionary relay, only much later during the Cretaceous.’ However, these authors refer to abundance, presumably in the fossil record, which in the case of squamates coincides with an increase of diversity at the genus level (see how Appendix 1—figure 5 reflects this increase in sampled taxa in the morphological matrix used for this study). This Cretaceous diversification, coinciding with the KTR, can be interpreted as a secondary radiation of the crown into modern clades, and is thus not strictly related to the primary radiation we describe in the Middle–Late Jurassic. Another point that favors our interpretations is that dated phylogenies, and no matter what method is used to date them, suggest that many of the main squamate clades had already diverged by the Late Jurassic (Figure 1—figure supplement 1). Because we know that the fossil record is fragmentary and uneven, it is even possible that the addition of new taxa/fossils could move this radiation backwards, but it is unlikely that an improved fossil record would result in a movement of this radiation towards the Cretaceous. The discoveries of the last decades from the Jurassic of China (e.g., Sullivan et al., 2014) regarding dinosaurs (especially feathered avian dinosaurs), but also mammals and lizards show that our knowledge of the timing of events can radically change with the reporting of new deposits containing key fossils, and have a tendency to move backwards in time. The Early Cretaceous of China has shown that squamates from this age already present extreme adaptations (an example would be Xianglong, the earliest known gliding squamate described by Li et al., 2007), supporting the view that important events in the radiation of squamates occurred before the Late Cretaceous. Accordingly, we suggest that more attention should be paid to events occurring in the Jurassic that could have triggered the diversification and disparity expansion of squamates.

The Jurassic sees a first major diversification in terms of taxa, but also of morphotypes. As explained above, all four major morphotypes of lepidosaurs (rhynchocephalians, generalized lizards, anguimorphs, and the limbless forms as represented by snakes) were already present by the Middle Jurassic. Finally, it seems that the KTR might have had a moderate influence on the evolutionary history of squamates, likely overinterpreted by an enhanced fossil record from the middle–Late Cretaceous. Most of the extant clades of squamates were already present in the Mesozoic (many of them as early as the Middle–Late Jurassic), and differences between Mesozoic and post-Mesozoic assemblages are the result of the mixture of the extinction of a low number of clades at the K-Pg boundary, and the rise of a few new clades (e.g., lacertids, amphisbaenians) along the Paleogene. Considering that the groups of lepidosaurs containing the taxa with a largest body size (mosasaurs, borioteiioids, and most rhynchocephalians) seem to have been most affected by the K-Pg extinction, small size might be one of the reasons for a greater survival of squamates, although providing statistical evidence for this is beyond the scope of our work.

Funding Statement

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Contributor Information

Arnau Bolet, Email: arnau.bolet@icp.cat.

Min Zhu, Institute of Vertebrate Paleontology and Paleoanthropology, Chinese Academy of Sciences, China.

George H Perry, Pennsylvania State University, United States.

Funding Information

This paper was supported by the following grants:

Royal Society NF170464 to Arnau Bolet.
Ministerio de Ciencia, Innovación y Universidades IJC2018-037685-I to Arnau Bolet.
European Commission ERC 788203 to Michael Benton.
Natural Environment Research Council NE/I027630/1 to Michael J Benton.

Additional information

Competing interests

No competing interests declared.

Author contributions

Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Visualization, Writing - original draft, Writing – review and editing.

Methodology, Validation, Writing – review and editing.

Conceptualization, Validation, Writing – review and editing.

Additional files

Supplementary file 1. Interactive plot of 3D morphospace (PCO1, PCO2, and PCO3), colors according to high taxonomical groups.

elife-66511-supp1.zip^{(1.4MB, zip)}

Supplementary file 2. Interactive plot of 3D morphospace (PCO1, PCO2, and PCO3), colors according to low taxonomical groups.

elife-66511-supp2.zip^{(1.4MB, zip)}

Supplementary file 3. Interactive plot of 3D morphospace (PCO1, PCO2, and PCO3), colors according to time bins.

elife-66511-supp3.zip^{(1.4MB, zip)}

Supplementary file 4. Interactive plot of 3D phylomorphospace (PCO1, PCO2, and PCO3), colors according to high taxonomical groups.

elife-66511-supp4.zip^{(160.3KB, zip)}

Supplementary file 5. Interactive plot of 3D phylomorphospace (PCO1, PCO2, and PCO3), colors according to low taxonomical groups.

elife-66511-supp5.zip^{(160.5KB, zip)}

Supplementary file 6. Interactive plot of squamate 3D morphospace (PCO1, PCO2, and PCO3).

Black spheres correspond to toxicoferans, white spheres correspond to non-toxicoferans.

elife-66511-supp6.zip^{(157.3KB, zip)}

Transparent reporting form

elife-66511-transrepform1.pdf^{(322.9KB, pdf)}

Source code 1. R script and source files for morphospace and disparity analyses.

elife-66511-code1.zip^{(362.5KB, zip)}

Source code 2. R script and source files for evolutionary rate analyses.

elife-66511-code2.zip^{(80.6KB, zip)}

Source code 3. R script and source files for plotting morphospace clusters.

elife-66511-code3.zip^{(136.6KB, zip)}

Source code 4. R script and source files for plotting Crámer values.

elife-66511-code4.zip^{(5.7MB, zip)}

Source code 5. R script and source files for plotting consensus trees.

elife-66511-code5.zip^{(13.5KB, zip)}

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files. More specifically, all necessary files (Nexus files and source tables) are included in the Source code 1–4.

References

Adams D, Collyer M, Kaliontzopoulou A. geomorph: an r package for the collection and analysis of geometric morphometric shape data. Methods in Ecology and Evolution. 2019;1:12035. doi: 10.1111/2041-210X.12035. [DOI] [Google Scholar]
Bapst DW. Paleotree: an R package for paleontological and phylogenetic analyses of evolution. Methods in Ecology and Evolution. 2012;3:803–807. doi: 10.1111/j.2041-210X.2012.00223.x. [DOI] [Google Scholar]
Benson RBJ, Campione NE, Carrano MT, Mannion PD, Sullivan C, Upchurch P, Evans DC. Rates of dinosaur body mass evolution indicate 170 million years of sustained ecological innovation on the avian stem lineage. PLOS Biology. 2014;12:e1001853. doi: 10.1371/journal.pbio.1001853. [DOI] [PMC free article] [PubMed] [Google Scholar]
Benton MJ. The origins of modern biodiversity on land. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences. 2010;365:3667–3679. doi: 10.1098/rstb.2010.0269. [DOI] [PMC free article] [PubMed] [Google Scholar]
Bolet A, Evans SE. New material of the enigmatic Scandensia, an Early Cretaceous lizard from the Iberian Peninsula. Spec. Pap. Palaeontol. 2011;86:99–108. doi: 10.1111/J.1475-4983.2011.01074.x. [DOI] [Google Scholar]
Bolet A, Evans SE. A tiny lizard (Lepidosauria, Squamata) from the Lower Cretaceous of Spain. Palaeontology. 2012;55:491–500. doi: 10.1111/j.1475-4983.2012.01145.x. [DOI] [Google Scholar]
Brusatte SL, Benton MJ, Ruta M, Lloyd GT. Superiority, competition, and opportunism in the evolutionary radiation of dinosaurs. Science (New York, N.Y.) 2008;321:1485–1488. doi: 10.1126/science.1161833. [DOI] [PubMed] [Google Scholar]
Brusatte SL, O’Connor JK, Jarvis ED. The Origin and Diversification of Birds. Current Biology. 2015;25:R888–R898. doi: 10.1016/j.cub.2015.08.003. [DOI] [PubMed] [Google Scholar]
Burbrink FT, Grazziotin FG, Pyron RA, Cundall D, Donnellan S, Irish F, Keogh JS, Kraus F, Murphy RW, Noonan B, Raxworthy CJ, Ruane S, Lemmon AR, Lemmon EM, Zaher H. Interrogating Genomic-Scale Data for Squamata (Lizards, Snakes, and Amphisbaenians) Shows no Support for Key Traditional Morphological Relationships. Systematic Biology. 2020;69:502–520. doi: 10.1093/sysbio/syz062. [DOI] [PubMed] [Google Scholar]
Caldwell MW, Nydam RL, Palci A, Apesteguía S. The oldest known snakes from the Middle Jurassic-Lower Cretaceous provide insights on snake evolution. Nature Communications. 2015;6:5996. doi: 10.1038/ncomms6996. [DOI] [PubMed] [Google Scholar]
Cardinal S, Danforth BN. Bees diversified in the age of eudicots. Proceedings. Biological Sciences. 2013;280:20122686. doi: 10.1098/rspb.2012.2686. [DOI] [PMC free article] [PubMed] [Google Scholar]
Čerňanský A, Bolet A, Müller J, Rage J-C, Augé M, Herrel A. A new exceptionally preserved specimen of Dracaenosaurus (Squamata, Lacertidae) from the Oligocene of France as revealed by micro-computed tomography. Journal of Vertebrate Paleontology. 2017;37:e1384738. doi: 10.1080/02724634.2017.1384738. [DOI] [Google Scholar]
Cleary TJ, Benson RBJ, Evans SE, Barrett PM. Lepidosaurian diversity in the Mesozoic-Palaeogene: the potential roles of sampling biases and environmental drivers. Royal Society Open Science. 2018;5:171830. doi: 10.1098/rsos.171830. [DOI] [PMC free article] [PubMed] [Google Scholar]
Close RA, Friedman M, Lloyd GT, Benson RBJ. Evidence for a Mid-Jurassic Adaptive Radiation in Mammals. Current Biology. 2015;25:2137–2142. doi: 10.1016/j.cub.2015.06.047. [DOI] [PubMed] [Google Scholar]
Close RA, Benson RBJ, Alroy J, Behrensmeyer AK, Benito J, Carrano MT, Cleary TJ, Dunne EM, Mannion PD, Uhen MD, Butler RJ. Diversity dynamics of Phanerozoic terrestrial tetrapods at the local-community scale. Nature Ecology & Evolution. 2019;3:590–597. doi: 10.1038/s41559-019-0811-8. [DOI] [PubMed] [Google Scholar]
Conrad JL. Phylogeny And Systematics Of Squamata (Reptilia) Based On Morphology. Bulletin of the American Museum of Natural History. 2008;310:1–182. doi: 10.1206/310.1. [DOI] [Google Scholar]
Conrad JL, Ast JC, Montanari S, Norell MA. A combined evidence phylogenetic analysis of Anguimorpha (Reptilia: Squamata. Cladistics : The International Journal of the Willi Hennig Society. 2011;27:230–277. doi: 10.1111/j.1096-0031.2010.00330.x. [DOI] [PubMed] [Google Scholar]
Conrad JL, Balcarcel AM, Mehling CM. Earliest example of a giant monitor lizard (Varanus, Varanidae, Squamata. PLOS ONE. 2012;7:e41767. doi: 10.1371/journal.pone.0041767. [DOI] [PMC free article] [PubMed] [Google Scholar]
Conrad JL. A new lizard (Squamata) was the last meal of Compsognathus (Theropoda: Dinosauria) and is A holotype in A holotype. Zoological Journal of the Linnean Society. 2018;183:584–634. doi: 10.1093/zoolinnean/zlx055. [DOI] [Google Scholar]
Datta PM, Ray S. Earliest lizard from the late triassic (Carnian) of India. Journal of Vertebrate Paleontology. 2006;26:795–800. doi: 10.1671/0272-4634. [DOI] [Google Scholar]
Doyle JA. Integrating Molecular Phylogenetic and Paleobotanical Evidence on Origin of the Flower. International Journal of Plant Sciences. 2008;169:816–843. doi: 10.1086/589887. [DOI] [Google Scholar]
Estes R. Sauria Terrestria, Amphisbaenia (Handbuch Der Palaoherpetologie 10A. Stuttgart: Gustav Fischer Verlag; 1983. [Google Scholar]
Estes R, Queiroz KD, Gauthier J. In: Phylogenetic Relationships of the Lizard Families. Estes R, Pregill G, editors. Stanford: Stanford University Press; 1988. Phylogenetic relationships within Squamata; pp. 119–281. [Google Scholar]
Evans SE. Lizards: evolution, early radiation and biogeography. Proceedings of the 6th symposium on Mesozoic terrestrial ecosystems and Biota. Short pap.1995. [Google Scholar]
Evans SE. Crown Group Lizards (Reptilia, Squamata) from the Middle Jurassic of the British Isles. Palaeontographica Abteilung A. 1998;250:123–154. doi: 10.1127/pala/250/1998/123. [DOI] [Google Scholar]
Evans SE, Barbadillo LJ. An unusual lizard (Reptilia: Squamata) from the Early Cretaceous of Las Hoyas, Spain. Zoological Journal of the Linnean Society. 1998;124:235–265. doi: 10.1111/j.1096-3642.1998.tb00576.x. [DOI] [Google Scholar]
Evans SE, Barbadillo LJ. A short-limbed lizard from the Lower Cretaceous of Spain. Special Papers in Palaeontology. 1999;60:73–85. doi: 10.1111/j.1475-4983.2012.01145. [DOI] [Google Scholar]
Evans SE. At the feet of the dinosaurs: the early history and radiation of lizards. Biological Reviews of the Cambridge Philosophical Society. 2003;78:513–551. doi: 10.1017/s1464793103006134. [DOI] [PubMed] [Google Scholar]
Evans SE, Raia P, Barbera C. The Lower Cretaceous lizard genus Chometokadmon from Italy. Cretaceous Research. 2006;27:673–683. doi: 10.1016/j.cretres.2006.03.004. [DOI] [Google Scholar]
Evans SE. In: The Skull of Lepidosauria. Gans C, Gaunt AS, editors. Society for the Study of Amphibians and Reptiles; 2008. Biology of the Reptilia; pp. 1–347. [Google Scholar]
Fry BG, Vidal N, Norman JA, Vonk FJ, Scheib H, Ramjan SFR, Kuruppu S, Fung K, Hedges SB, Richardson MK, Hodgson WC, Ignjatovic V, Summerhayes R, Kochva E. Early evolution of the venom system in lizards and snakes. Nature. 2006;439:584–588. doi: 10.1038/nature04328. [DOI] [PubMed] [Google Scholar]
Gauthier JA, Kearney M, Maisano JA, Rieppel O, Behlke ADB. Assembling the Squamate Tree of Life: Perspectives from the Phenotype and the Fossil Record. Bulletin of the Peabody Museum of Natural History. 2012;53:3–308. doi: 10.3374/014.053.0101. [DOI] [Google Scholar]
Gerber S, Ruta M. Use and misuse of discrete character data for morphospace and disparity analyses. Palaeontology. 2019;62:305–319. doi: 10.1111/pala.12407. [DOI] [Google Scholar]
Goloboff PA, Catalano SA. TNT version 1.5, including a full implementation of phylogenetic morphometrics. Cladistics : The International Journal of the Willi Hennig Society. 2016;32:221–238. doi: 10.1111/cla.12160. [DOI] [PubMed] [Google Scholar]
Griffiths EF, Ford DP, Benson RBJ, Evans SE, Ruta M. New information on the Jurassic lepidosauromorph Marmoretta oxoniensis. Papers in Palaeontology. 2021;7:2255–2278. doi: 10.1002/spp2.1400. [DOI] [Google Scholar]
Guillerme T, Poisot T. dispRity: A modular R package for measuring disparity. Methods in Ecology and Evolution. 2018;9:1755–1763. doi: 10.1111/2041-210X.13022. [DOI] [Google Scholar]
Hedman MM. Constraints on clade ages from fossil outgroups. Paleobiology. 2010;36:16–31. doi: 10.1666/0094-8373-36.1.16. [DOI] [Google Scholar]
Herrera-Flores JA, Stubbs TL, Benton MJ, Ruta M. Macroevolutionary patterns in Rhynchocephalia: is the tuatara ( Sphenodon punctatus ) a living fossil? Palaeontology. 2017;60:319–328. doi: 10.1111/pala.12284. [DOI] [Google Scholar]
Herrera-Flores JA, Elsler A, Stubbs TL, Benton MJ, Benson R. Slow and fast evolutionary rates in the history of lepidosaurs. Palaeontology. 2021a;65:12579. doi: 10.1111/pala.12579. [DOI] [Google Scholar]
Herrera-Flores JA, Stubbs TL, Benton MJ. Ecomorphological diversification of squamates in the Cretaceous. Royal Society Open Science. 2021b;8:201961. doi: 10.1098/rsos.201961. [DOI] [PMC free article] [PubMed] [Google Scholar]
Howard TJ, Swofford JB, Wagner DL, Sherman S, Lehman GA. Quality of life after bilateral thoracoscopic splanchnicectomy: long-term evaluation in patients with chronic pancreatitis. Journal of Gastrointestinal Surgery. 2002;6:845–852. doi: 10.1016/s1091-255x(02)00123-3. [DOI] [PubMed] [Google Scholar]
Hutchinson MN, Skinner A, Lee MSY. Tikiguania and the antiquity of squamate reptiles (lizards and snakes. Biology Letters. 2012;8:665–669. doi: 10.1098/rsbl.2011.1216. [DOI] [PMC free article] [PubMed] [Google Scholar]
Jones MEH, Anderson CL, Hipsley CA, Müller J, Evans SE, Schoch RR. Integration of molecules and new fossils supports a Triassic origin for Lepidosauria (lizards, snakes, and tuatara. BMC Evolutionary Biology. 2013;13:208. doi: 10.1186/1471-2148-13-208. [DOI] [PMC free article] [PubMed] [Google Scholar]
Kearney M. The phylogenetic position of Sineoamphisbaena hexatabularis reexamined. Journal of Vertebrate Paleontology. 2003;23:394–403. doi: 10.1671/0272-4634(2003)023[0394:TPPOSH]2.0.CO;2. [DOI] [Google Scholar]
Kotrc B, Knoll AH. A morphospace of planktonic marine diatoms. I. Two views of disparity through time. Paleobiology. 2015;41:45–67. doi: 10.1017/pab.2014.4. [DOI] [Google Scholar]
Lafuma F, Corfe IJ, Clavel J, Di-Poï N. Multiple evolutionary origins and losses of tooth complexity in squamates. Nature Communications. 2021;12:6001. doi: 10.1038/s41467-021-26285-w. [DOI] [PMC free article] [PubMed] [Google Scholar]
Laurin M. The evolution of body size, Cope’s rule and the origin of amniotes. Systematic Biology. 2004;53:594–622. doi: 10.1080/10635150490445706. [DOI] [PubMed] [Google Scholar]
Lee MSY, Cau A, Naish D, Dyke GJ. Dinosaur evolution. Sustained miniaturization and anatomical innovation in the dinosaurian ancestors of birds. Science (New York, N.Y.) 2014;345:562–566. doi: 10.1126/science.1252243. [DOI] [PubMed] [Google Scholar]
Li PP, Gao KQ, Hou LH, Xu X. A gliding lizard from the Early Cretaceous of China. PNAS. 2007;104:5507–5509. doi: 10.1073/pnas.0609552104. [DOI] [PMC free article] [PubMed] [Google Scholar]
Lloyd GT, Davis KE, Pisani D, Tarver JE, Ruta M, Sakamoto M, Hone DWE, Jennings R, Benton MJ. Dinosaurs and the Cretaceous Terrestrial Revolution. Proceedings. Biological Sciences. 2008;275:2483–2490. doi: 10.1098/rspb.2008.0715. [DOI] [PMC free article] [PubMed] [Google Scholar]
Lloyd GT. Estimating morphological diversity and tempo with discrete character-taxon matrices: implementation, challenges, progress, and future directions. Biological Journal of the Linnean Society. 2016;118:131–151. doi: 10.1111/bij.12746. [DOI] [Google Scholar]
Lloyd GT, Bapst DW, Friedman M, Davis KE. Probabilistic divergence time estimation without branch lengths: dating the origins of dinosaurs, avian flight and crown birds. Biology Letters. 2016;12:20160609. doi: 10.1098/rsbl.2016.0609. [DOI] [PMC free article] [PubMed] [Google Scholar]
Longrich NR, Bhullar BAS, Gauthier JA. Mass extinction of lizards and snakes at the Cretaceous-Paleogene boundary. PNAS. 2012;109:21396–21401. doi: 10.1073/pnas.1211526110. [DOI] [PMC free article] [PubMed] [Google Scholar]
Longrich NR, Vinther J, Pyron RA, Pisani D, Gauthier JA. Biogeography of worm lizards (Amphisbaenia) driven by end-Cretaceous mass extinction. Proceedings of the Royal Society B. 2015;282:20143034. doi: 10.1098/rspb.2014.3034. [DOI] [PMC free article] [PubMed] [Google Scholar]
Losos JB, Hillis DM, Greene HW. Evolution. Who speaks with a forked tongue? Science (New York, N.Y.) 2012;338:1428–1429. doi: 10.1126/science.1232455. [DOI] [PubMed] [Google Scholar]
Martínez RN, Simões TR, Sobral G, Apesteguía S. A Triassic stem lepidosaur illuminates the origin of lizard-like reptiles. Nature. 2021;597:235–238. doi: 10.1038/s41586-021-03834-3. [DOI] [PubMed] [Google Scholar]
Meredith RW, Janečka JE, Gatesy J, Ryder OA, Fisher CA, Teeling EC, Goodbla A, Eizirik E, Simão TLL, Stadler T, Rabosky DL, Honeycutt RL, Flynn JJ, Ingram CM, Steiner C, Williams TL, Robinson TJ, Burk-Herrick A, Westerman M, Ayoub NA, Springer MS, Murphy WJ. Impacts of the Cretaceous Terrestrial Revolution and KPg extinction on mammal diversification. Science (New York, N.Y.) 2011;334:521–524. doi: 10.1126/science.1211028. [DOI] [PubMed] [Google Scholar]
Mongiardino Koch N, Gauthier JA. Noise and biases in genomic data may underlie radically different hypotheses for the position of Iguania within Squamata. PLOS ONE. 2018;13:e0202729. doi: 10.1371/journal.pone.0202729. [DOI] [PMC free article] [PubMed] [Google Scholar]
Moon BC, Stubbs TL. Early high rates and disparity in the evolution of ichthyosaurs. Communications Biology. 2020;3:68. doi: 10.1038/s42003-020-0779-6. [DOI] [PMC free article] [PubMed] [Google Scholar]
Mulcahy DG, Noonan BP, Moss T, Townsend TM, Reeder TW, Sites JW, Wiens JJ. Estimating divergence dates and evaluating dating methods using phylogenomic and mitochondrial data in squamate reptiles. Molecular Phylogenetics and Evolution. 2012;65:974–991. doi: 10.1016/j.ympev.2012.08.018. [DOI] [PubMed] [Google Scholar]
Müller J, Hipsley CA, Head JJ, Kardjilov N, Hilger A, Wuttke M, Reisz RR. Eocene lizard from Germany reveals amphisbaenian origins. Nature. 2011;473:364–367. doi: 10.1038/nature09919. [DOI] [PubMed] [Google Scholar]
Nordén KK, Stubbs TL, Prieto-Márquez A, Benton MJ. Multifaceted disparity approach reveals dinosaur herbivory flourished before the end-Cretaceous mass extinction. Paleobiology. 2018;44:620–637. doi: 10.1017/pab.2018.26. [DOI] [Google Scholar]
Pianka ER, Vitt LJ, Pelegrin N, Fitzgerald DB, Winemiller KO. Toward a Periodic Table of Niches, or Exploring the Lizard Niche Hypervolume. The American Naturalist. 2017;190:601–616. doi: 10.1086/693781. [DOI] [PubMed] [Google Scholar]
Puttick MN, Thomas GH, Benton MJ. High rates of evolution preceded the origin of birds. Evolution; International Journal of Organic Evolution. 2014;68:1497–1510. doi: 10.1111/evo.12363. [DOI] [PMC free article] [PubMed] [Google Scholar]
Pyron RA, Burbrink FT, Wiens JJ. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evolutionary Biology. 2013;13:1–56. doi: 10.1186/1471-2148-13-93. [DOI] [PMC free article] [PubMed] [Google Scholar]
Pyron RA. Novel Approaches for Phylogenetic Inference from Morphological Data and Total-Evidence Dating in Squamate Reptiles (Lizards, Snakes, and Amphisbaenians. Systematic Biology. 2017;66:38–56. doi: 10.1093/sysbio/syw068. [DOI] [PubMed] [Google Scholar]
R Development Core Team . Vienna, Austria: R Foundation for Statistical Computing; 2013. http://www.R-project.org/ [Google Scholar]
Raup DM. Biases in the fossil record of species and genera. Bulletin of Carnegie Museum of Natural History. 1979;13:85–91. [Google Scholar]
Revell LJ. phytools: an R package for phylogenetic comparative biology (and other things. Methods in Ecology and Evolution. 2012;3:217–223. doi: 10.1111/j.2041-210X.2011.00169.x. [DOI] [Google Scholar]
Reynoso VH. Huehuecuetzpalli mixtecus gen. et sp. nov: a basal squamate (Reptilia) from the Early Cretaceous of Tepexi de Rodríguez, Central México. Philosophical Transactions of the Royal Society of London. Series B. 1998;353:477–500. doi: 10.1098/rstb.1998.0224. [DOI] [Google Scholar]
Reynoso VH, Callison G. A new scincomorph lizard from the Early Cretaceous of Puebla, Mexico. Zoological Journal of the Linnean Society. 2000;130:183–212. doi: 10.1111/j.1096-3642.2000.tb01629.x. [DOI] [Google Scholar]
Schwenk K. In: Phylogenetic Relationships of the Lizard Families. Estes R, Pregill G, editors. Stanford: Stanford University Press; 1988. Comparative morphology of the lepidosaur tongue and its relevance to squamate phylogeny; pp. 569–597. [Google Scholar]
Sepkoski JJ. A factor analytic description of the Phanerozoic marine fossil record. Paleobiology. 2016;7:36–53. doi: 10.1017/S0094837300003778. [DOI] [Google Scholar]
Sievert C. ggplot2; 2018. https://plotly-r.com [Google Scholar]
Simões TR, Caldwell MW, Tałanda M, Bernardi M, Palci A, Vernygora O, Bernardini F, Mancini L, Nydam RL. The origin of squamates revealed by a Middle Triassic lizard from the Italian Alps. Nature. 2018;557:706–709. doi: 10.1038/s41586-018-0093-3. [DOI] [PubMed] [Google Scholar]
Simões TR, Vernygora O, Caldwell MW, Pierce SE. Megaevolutionary dynamics and the timing of evolutionary innovation in reptiles. Nature Communications. 2020;11:3322. doi: 10.1038/s41467-020-17190-9. [DOI] [PMC free article] [PubMed] [Google Scholar]
Smith KT. Eocene Lizards of the Clade Geiseltaliellus from Messel and Geiseltal, Germany, and the Early Radiation of Iguanidae (Reptilia: Squamata. Bulletin of the Peabody Museum of Natural History. 2009;50:219–306. doi: 10.3374/014.050.0201. [DOI] [Google Scholar]
Sullivan C, Wang Y, Hone DWE, Wang Y, Xu X, Zhang F. The vertebrates of the Jurassic Daohugou Biota of northeastern China. Journal of Vertebrate Paleontology. 2014;34:243–280. doi: 10.1080/02724634.2013.787316. [DOI] [Google Scholar]
Uetz P, Freed P, Hošek J. The Reptile Database. 2019. [February 19, 2020]. http://www.reptile-database.org
Wainwright PC, Alfaro ME, Bolnick DI, Hulsey CD. Many-to-One Mapping of Form to Function: A General Principle in Organismal Design? Integrative and Comparative Biology. 2005;45:256–262. doi: 10.1093/icb/45.2.256. [DOI] [PubMed] [Google Scholar]
Wang M, Lloyd GT. Rates of morphological evolution are heterogeneous in Early Cretaceous birds. Proceedings of the Royal Society B. 2016;283:20160214. doi: 10.1098/rspb.2016.0214. [DOI] [PMC free article] [PubMed] [Google Scholar]
Watanabe A, Fabre AC, Felice RN, Maisano JA, Müller J, Herrel A, Goswami A. Ecomorphological diversification in squamates from conserved pattern of cranial integration. PNAS. 2019;116:14688–14697. doi: 10.1073/pnas.1820967116. [DOI] [PMC free article] [PubMed] [Google Scholar]
Wickham H. Ggplot2: Elegant Graphics for Data Analysis. Cham: Springer-Verlag; 2016. [DOI] [Google Scholar]
Zhang C, Wang M. Bayesian tip dating reveals heterogeneous morphological clocks in Mesozoic birds. Royal Society Open Science. 2019;6:182062. doi: 10.1098/rsos.182062. [DOI] [PMC free article] [PubMed] [Google Scholar]
Zheng Y, Wiens JJ. Combining phylogenomic and supermatrix approaches, and a time-calibrated phylogeny for squamate reptiles (lizards and snakes) based on 52 genes and 4162 species. Molecular Phylogenetics and Evolution. 2016;94:537–547. doi: 10.1016/j.ympev.2015.10.009. [DOI] [PubMed] [Google Scholar]

eLife. doi: 10.7554/eLife.66511.sa0

Editor's evaluation

Min Zhu ¹

This article presents an evaluation of the macroevolutionary history of squamates (lizards, snakes, and relatives) and is relevant to evolutionary biologists and paleontologists interested in this group. The ‘early burst’ of disparity in squamates demonstrates that squamates established their morphospace range much earlier than had been assumed, and the long-term stable morphospace occupation ever since.

eLife. doi: 10.7554/eLife.66511.sa1

Decision letter

Editor: Min Zhu¹

Reviewed by: Min Zhu²

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Decision letter after peer review:

Thank you for submitting your article "The Jurassic rise of squamates as revealed by lepidosaur disparity and evolutionary rates" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Min Zhu as Reviewing Editor and Reviewer #1, and the evaluation has been overseen by George Perry as the Senior Editor.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The present work should be of interest to evolutionary biologists concerned with macroevolution as well as taxonomic specialists. The main conclusions are supported by the data with caution, and they complement and refine some other recent results. The authors go to considerable lengths to rule out alternative explanations. Overall, while the conclusion of this manuscript is novel, many concerns addressed by the reviewers (below) should be considered for revisions and improvement. The manuscript could be published without the need for new analyses, although we would like to see more, especially how rates are affected by enforcing divergence times more in line with molecular trees. If new analyses are not performed, then the manuscript should be toned down to account for the many sources of uncertainty that permeate the present analyses, and acknowledge the ways in which the results depend on methodological decisions that could be improved upon.

(1) The discussion of ecology (lines 95 and following) is oversimplified (as the authors recognise), perhaps inevitably so.

The authors recognise that the "ecology" of the clusters discovered by the authors are difficult to generalise (what is the ecology of limbless forms that unites scolecophidians and boas?). Given that, we would consider striking "ecology" from the title and focus just on disparity.

(2) While we believe that the authors have done sufficient due diligence to rule out potential confounding factors, it would still be worthwhile to be more explicit about one limitation, namely the issue of missing data and inapplicable characters. For instance, this paper could be brought up:

Gerber, S. (2019). Use and misuse of discrete character data for morphospace and disparity analyses. Palaeontology, 62(2), 305-319.

We are not convinced that a broad-scale expunging of "inapplicable" characters (say) is the best way to solve the issue (one thinks of limb character in limbless forms), but a more explicit acknowledgment of this limitation would make the paper more compelling.

(3) One weakness is the issue of completeness, which has been raised prominently in recent years. The authors approach this by way of examining Cramér coefficients and throwing out characters that are generally missing for fossils (like those of integument). While this is not unreasonable, and results in disparity plots that are broadly comparable to analyses in which these characters are present, there remain other characters (especially those related to the limbs), this problem was less thoroughly addressed than others.

(4) One of the striking aspects of the results is that the interpreted expansion of morphospace in both the Jurassic and Cretaceous occurs before the interpreted spike in evolutionary rates (not exactly "roughly coincident" as stated). This seemingly counterintuitive result is not addressed.

(5) As someone who has performed principal coordinate analyses on lepidosaur morphological datasets, one reviewer can confirm that the first few dimensions invariably separate limb-bearing from limbless forms, as well as rhynchocephalians from squamates (in the case of this analysis, anguimorphs seem to also contribute to this axis). This is exactly what is found in Figure 1b, but also previously in Watanabe et al., (2019) and Simoes et al., (2020). As such, separating these four groups (rhynchocephalians, snakes, anguimorphs and "generalized lizards", as defined in lines 95-104 by the authors) explains a large proportion of total variance, which loads prominently on the first axes but is likely to contribute to some extent to others. Given that the matrix being used was "designed to identify the relationships among the major squamate clades" (Conrad 2018), it is no surprise that the same pattern emerges in this study. The "Jurassic morphospace expansion" (JMA) proposed here is nothing more than a period of time when a morphospace occupied by only one of these groups (rhynchocepalinas) gives way to one occupied by all four. This result is not novel (i.e., does not depend on any discovery made in this study), it is just presented here in a novel and quantitative way. In fact it is already stated in some of the first lines of the introduction that these morpho-groups are known to already be present in the Middle and Late Jurassic (lines 47-49), and that a similar pattern had already been mentioned in the literature (line 238). Furthermore, whether this does represent a true event of expansion or an artifact of the fossil record is also unclear given the poor fossil record of the clade during this period, as noted by the authors. Furthermore, it ultimately depends on many methodological decisions not entirely validated (such as temporal binning, choice of distance metric, etc), as well as character and taxon sampling decisions that were not optimized for the purpose of macroevolutionary inference, and might therefore supported biased patterns, not just during the Jurassic but also later on.

(6) The discovery of this peak in disparity if further linked to high rates of morphological evolution. Note however that disparity and rates are not necessarily linked, and high rates are not necessary to support a JME. However, this result is entirely dependent on the topology and divergence times supported by this study. The methods of inference and calibration employed here are much more simplistic than those of previous studies, and allow for a much poorer integration of uncertainty into subsequent analyses. Results differ crucially from those of all previous studies, including a much faster early diversification of squamates, and much older ages for some crown clades such as snakes, anguimorphs and iguanians. This leads to a concentration of branches with high levels of morphological change in the mid to late Jurassic, ultimately responsible for the peak in rates found. The methods employed and the way in which they are presented once again do not rule out the possibility that this is an artifact of the way the study was set up, and a major reason why no such increase in rates was discovered by previous studies.

(7) There is a lack of consistency throughout the manuscript. For example, the JME label applied to each of the four panels of Figure 1 seems to represent four different periods of time. The topology in Figure 1A was calibrated using a different method that the rates shown in Figure 1D, making it difficult to draw conclusions. While a staggering 51 supplementary figures are presented, many of these are unedited and thus highly cryptic, others such as the many ordination plots are unlabeled and thus do not contribute much.

Some aspects of this manuscript are very interesting, such as the position of extinct lineages under topologies constrained to show different relationships among major extant clades. However, its main claim is somewhere between trivial (a change from 1 to 4 morpho-groups somewhere in the second half of the Jurassic) and potentially unsubstantiated (a coincident peak in rates).

The current narrative, emphasizing a completely new and unexpected view on squamate diversification might be true, but it rests on many assumptions. A more realistic interpretation should focus on the existence of four morpho-groups within squamates (which could even be tested with this data, see below), and their timing of origin. The manuscript starts by posing the question "why are some groups highly species-rich and others are not?" (line 16). This question is never answered, and in fact the manuscript does not seek to do so or use data that could. But that kind of exaggeration permeates the rest of the manuscript, and the reviewer thinks it will be benefit from a more realistic approach to the data at hand.

(8) One reviewer thinks the paper should validate some of their methodological choices better, or even better improve upon them. For example, two disparity metrics are employed (one of which does not even show a JME), but no reason is given why these were chosen. Results from other post-ordination metrics should be presented as well, along with some validation of why certain metrics where favored. This is implemented by Guillerme et al., 2020 Eco Evo and employs some of the same programs already used. Furthermore, a large part of the findings relies on time-calibrated analyses that are poorly explained. For example, although 3 methods of post hoc scaling of topologies are explored, only 1 topology is presented as a supplementary figure. Details on the methods themselves are missing, as it is unclear what time information went into the calibration and how it was gathered and treated. Ultimately however, branch lengths scaling seems to be resulting in unlikely ages for many clades. The reviwer believes results need to be repeated using a topology calibrated using molecular data, or at the very least including secondary node constraints based on other published molecular studies.

One aspect that could also help improve this paper is an exploration of morpho-functional groups within morphospace. Currently, traditional views are simply repeated (see above regarding the four groups of lepidosaurs) and forced onto the results rather than extracted from them. What the data at hand could help establish (possibly along other datasets that are said to have been dropped due to poor sampling of Jurassic fossils) is how many morphological clusters of taxa can be found (using k-means or expectation maximization algorithms for example), and when in time do those originate based on the topologies at hand (not just the one used here but those previously published as well).

(9) The earliest fossil records of the crown squamates are the Middle Jurassic. All these records (for some constituent major clades) in the Middle Jurassic should be clearly stated and/or cited in the text. As one reviewer thinks, this is why the authors separated Early Jurassic from Middle-Late Jurassic in their analyses. The criteria for this separation should be discussed. Otherwise, the readers would doubt why you did not separate Early-Middle Jurassic from Late Jurassic. Figure 1a, which shows the earliest record of the crown squamates later than the Middle Jurassic, is misleading in this regard.

(10) The manuscript should be edited carefully following the author guidelines of eLife. Some of the discussions in the supplementary information text can be moved to the main text to gain a better readership. The Discussion and especially Extended Discussion would benefit from some further subdivision. Some of the supplementary figures (Extended Data Figures 1-10), which are not well demonstrated, can be deleted or replaced with a batch file. Overall, the supplementary figures should be re-organized to have a uniform format. The authors should put more effort into the supplementary figures, many of which are unintelligible.

Other claims that need to be revisited:

(1) Line 132 – "These evolutionary rates are not dependent on phylogeny": They are very dependent on phylogeny. This analysis shows they are not dependent on topology (molecular vs. morphological) but both of those are likely to share very similar branch lengths for all branches in common, making the result still strongly dependent on phylogeny.

(2) Lines 172-173 – "The second event, which fits well with the timing of the KTR, would be coincident with the radiation of the constituent crown groups of Squamata": Figure S1, which is the only time-calibrated topology shown in its entirety, shows crown snakes, anguimorphs, iguanians, lacertoids and skinks originating in the Jurassic, and radiating well before the KTR.

3) Lines 189-191 – "The only observable changes from then on correspond to a slight expansion of the edges of the occupied morphospace, and a notable increase in the density of points filling this morphospace": This is entirely contingent on analyses using a matrix that was tuned to capture the deep phylogeny of the clade, and is something that should be acknowledged.

4) Line 271 – "basal position of Iguania": This and other instances of the use of the word basal should be corrected, living groups cannot be basal to other living groups.

(5) The reviewer has not been able to find which dimensions were considered to estimate the disparity curve shown in Figure 1C. Please add this information to the caption.

6) The reviewer would consider striking some comments from the body of the paper that seem more like personal annotations ("We expected that….").

7) L.24: "the Middle and Late Jurassic (174-145 Ma)" is clearer than "the second half of the Jurassic (170-145 Ma)" to the readers, and also consistent with that in Line 48.

(8) L.82: Figure 1b-d as a separate figure?

9) L.90: For "pan-Squamata", do you mean the total group Squamata?

10) L.108, 171, 212: "mid-late Jurassic" Ambiguous here、Middle and Late Jurassic?

(11) L.145：'stem lepidosaurs' for 'ancestral lepidosaurs'

(12) L. 151: Extended Data Figure 11 as a main-text figure?

(13) L.172, 175, 'late', upper case.

(14) Extended Data Figures 23, 24, 43-47, 'early', 'middle', 'late', upper cases.

Reviewer #1:

This study reveals the explosive adaptive radiation of squamates in the Middle and Late Jurassic, situating the dates of origins of their constituent major clades. The "early burst" of disparity in squamates demonstrates that squamates established their morphospace range much earlier than had been assumed, and the long-term stable morphospace occupation ever since. Overall, the results are innovative and significant to evolutionary biology of lepidosaurs.

Reviewer #2:

The evolution of disparity, as distinct from species diversity, is a central question in macroevolution. The authors apply recently developed methods keyed to widely available, increasingly large-scale data sets (phylogenetic character matrices) to examine when major changes in body form occurred in the evolution of lizards and snakes, one of the most species-rich living vertebrate groups. They conclude that major changes in body form -- and an expansion of the parameter space occupied by lizards and snakes -- occurred tens of millions of years earlier than anticipated by most previous studies.

Strengths:

One of the major strengths of the disparity method (as is well recognised) is the fact that it is based on widely available data sets (phylogenetic character matrices), rather than the bespoke morphometric data sets traditionally employed in the study of disparity. This makes it much more widely applicable.

One of the major strengths of the paper is the lengths to which the authors go to rule out alternative explanations. For instance, the analyses of evolutionary rates were conducted under different assumptions of tree topology and divergence times to rule out a dependency on that factor (a potentially serious problem in squamates).

Weaknesses:

The discussion of ecology (lines 95 and following) is oversimplified (as the authors recognise), perhaps inevitably so.

One weakness is the issue of completeness, which has been raised prominently in recent years. The authors approach this by way of examining Cramér coefficients and throwing out characters that are generally missing for fossils (like those of integument). While this is not unreasonable, and results in disparity plots that are broadly comparable to analyses in which these characters are present, there remain other characters (especially those related to the limbs), this problem was less thoroughly addressed than others.

One of the striking aspects of the results is that the interpreted expansion of morphospace in both the Jurassic and Cretaceous occurs before the interpreted spike in evolutionary rates (not exactly "roughly coincident" as stated). This seemingly counterintuitive result is not addressed.Reviewer #3:

Bolet et al., explore the evolutionary history of lepidosaur morphology. Using a previously published dataset sampling a wide range of living and extinct lineages, they analyze both morphospace occupation and rates of evolution through time. By combining a series of phylogenetic and macroevolutionary methods, they reveal an event of morphospace expansion (coupled with increased evolutionary rates) during the second half of the Jurassic. This contrasts with previous studies, which had recovered a predominant signal of diversification and ecological innovation during the mid-Cretaceous, coinciding with the radiation of other major animal and plant lineages and a large-scale restructuring of terrestrial ecosystems. While their conclusion is novel, the manuscript is not based on any new morphological or paleontological data, and employs relatively simple methodologies that can potentially determine the outcome. My main concerns are the following:

(1) As someone who has performed principal coordinate analyses on lepidosaur morphological datasets, I can confirm that the first few dimensions invariably separate limb-bearing from limbless forms, as well as rhynchocephalians from squamates (in the case of this analysis, anguimorphs seem to also contribute to this axis). This is exactly what is found in Figure 1b, but also previously in Watanabe et al., (2019) and Simoes et al., (2020). As such, separating these four groups (rhynchocephalians, snakes, anguimorphs and "generalized lizards", as defined in lines 95-104 by the authors) explains a large proportion of total variance, which loads prominently on the first axes but is likely to contribute to some extent to others. Given that the matrix being used was "designed to identify the relationships among the major squamate clades" (Conrad 2018), it is no surprise that the same pattern emerges in this study. The "Jurassic morphospace expansion" (JMA) proposed here is nothing more than a period of time when a morphospace occupied by only one of these groups (rhynchocepalinas) gives way to one occupied by all four. This result is not novel (i.e., does not depend on any discovery made in this study), it is just presented here in a novel and quantitative way. In fact it is already stated in some of the first lines of the introduction that these morpho-groups are known to already be present in the Middle and Late Jurassic (lines 47-49), and that a similar pattern had already been mentioned in the literature (line 238). Furthermore, whether this does represent a true event of expansion or an artifact of the fossil record is also unclear given the poor fossil record of the clade during this period, as noted by the authors. Furthermore, it ultimately depends on many methodological decisions not entirely validated (such as temporal binning, choice of distance metric, etc), as well as character and taxon sampling decisions that were not optimized for the purpose of macroevolutionary inference, and might therefore supported biased patterns, not just during the Jurassic but also later on.

(2) The discovery of this peak in disparity if further linked to high rates of morphological evolution. Note however that disparity and rates are not necessarily linked, and high rates are not necessary to support a JME. However, this result is entirely dependent on the topology and divergence times supported by this study. The methods of inference and calibration employed here are much more simplistic than those of previous studies, and allow for a much poorer integration of uncertainty into subsequent analyses. Results differ crucially from those of all previous studies, including a much faster early diversification of squamates, and much older ages for some crown clades such as snakes, anguimorphs and iguanians. This leads to a concentration of branches with high levels of morphological change in the mid to late Jurassic, ultimately responsible for the peak in rates found. The methods employed and the way in which they are presented once again do not rule out the possibility that this is an artifact of the way the study was set up, and a major reason why no such increase in rates was discovered by previous studies.

(3) There is a lack of consistency throughout the manuscript. For example, the JME label applied to each of the four panels of Figure 1 seems to represent four different periods of time. The topology in Figure 1A was calibrated using a different method that the rates shown in Figure 1D, making it difficult to draw conclusions. While a staggering 51 supplementary figures are presented, many of these are unedited and thus highly cryptic, others such as the many ordination plots are unlabeled and thus do not contribute much.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "The Jurassic rise of squamates as supported by lepidosaur disparity and evolutionary rates" for further consideration by eLife. Your revised article has been evaluated by George Perry (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

See comments from Reviewer #3.

Reviewer #1:

The manuscript is well revised to address the comments from the reviewers.

Reviewer #3:

I congratulate the authors on the effort put towards revising their manuscript, and after carefully reading their response and new version of the manuscript I have just a handful of minor suggestions that I would like to see implemented, or more strongly stated, in the manuscript before the final version is published.

(1) Editing some of the supplementary files shown would not take that much time at all, and it would make them much more readable than they are now. Reducing the size of the label text in the trees now shown as supplementary figures to figure 1, and color coding them by clade, does not represent too much work. Same holds true for color coding by clade the dots in the many supplementary PCO plots.

(2) The authors state in their response that "Methods in some other recent studies include extensive genomic data as well, but such data is only available for modern taxa and can only change the overall tree topology, but not the dating of deep divergences, as here." Genomic data can most definitely change a lot more than just topology, and in fact do modify the dating of deep divergences. The methods implemented currently for time-scaling make a number of simplified decisions regarding the calculation of branch lengths that come from the fact that they do not incorporate any quantifiable (morphological or molecular) sense of the degree of divergence occurred between successive nodes. Methods directly employing genomic data have access to this amount of divergence, and can therefore not only incorporate it into the analysis but also eschew some of the simplistic assumptions of methods such as equal or mbl (which remains true regardless of the latter being very common in the literature). It is wrong of the authors to assume that time calibration relies exclusively on the use of the right fossils (although this is certainly crucial), and there are good reasons why approaches incorporating genomic data to divergence time estimation are unambiguously better than those that don't, even when molecular data is available only for extant terminals. I strongly believe that this caveat needs to be made explicit in the main text of the manuscript, and that it is not a problem that disappears just because several different methods explored agree with each other.

3) In a similar vein, I would suggest the authors to explicitly state that some results that directly depend on the estimation of time-calibrated branch lengths (such as the calculation of evolutionary rates) necessitate a confirmation with more robust methods of time calibration, such as those that incorporate molecular data alongside the fossil record, as well as employ more realistic models of diversification such as the fossilized birth-death prior.

eLife. 2022 May 3;11:e66511. doi: 10.7554/eLife.66511.sa2

Author response

Essential revisions:

The present work should be of interest to evolutionary biologists concerned with macroevolution as well as taxonomic specialists. The main conclusions are supported by the data with caution, and they complement and refine some other recent results. The authors go to considerable lengths to rule out alternative explanations. Overall, while the conclusion of this manuscript is novel, many concerns addressed by the reviewers (below) should be considered for revisions and improvement. The manuscript could be published without the need for new analyses, although we would like to see more, especially how rates are affected by enforcing divergence times more in line with molecular trees.

We agree with the reviewers in that the divergence ages provided by some methods (mainly the Hedman method) are older than those reported in other studies. However, the use of the Hedman method upon a different dataset (e.g. Herrera-Flores et al., 2021a) yielded similar ages. In that case, the Hedman method was the only method considered, when we have instead performed the analyses according to two additional dating methods that yield more recent divergence ages estimates. According to this, it cannot be considered that results are tied to the old divergence times resulting from applying the Hedman method, because they are recovered when using the equal method too. Our results using the MBL method are admittedly younger, but we think that the MBL suffers from some sort of pull of the recent effect. We also performed a preliminary Bayesian analysis to see if divergence times were in line with those from the parsimony methods just noted. Surprisingly, the ages recovered using MrBayes (preliminary results provided for evaluation, but not implemented in the manuscript) are more in line with those of the Hedman method (see the attached MCCT tree) than with the equal or MBL methods, which both provide younger ages. We are planning to repeat analyses forcing specific time divergences according to molecular studies, but this approach will require a careful assessment of the ages of fossils in our dataset that could represent older records than the ages of given clades in molecular analyses.

If new analyses are not performed, then the manuscript should be toned down to account for the many sources of uncertainty that permeate the present analyses, and acknowledge the ways in which the results depend on methodological decisions that could be improved upon.

(1) The discussion of ecology (lines 95 and following) is oversimplified (as the authors recognise), perhaps inevitably so.

The authors recognise that the "ecology" of the clusters discovered by the authors are difficult to generalise (what is the ecology of limbless forms that unites scolecophidians and boas?). Given that, we would consider striking "ecology" from the title and focus just on disparity.

This failure to recover fine-tuned ecomorphological clusters is a recurrent problem in studies of lepidosaurs, and seems to be independent of the type of dataset (geometric morphometrics in Watanabe et al., 2018 or Herrera-Flores et al., 2021a vs discrete morphological characters like in Simões et al., 2020, Martínez et al., 2021 or the present study) or differences in sampling (a broader sampling of non-lepidosaur diapsids of e.g. Simões et al., 2020 vs our broader sampling of lepidosaurs, specially Mesozoic squamates). In our opinion, the lack of resolution among ecological clusters is related to the fact that the latter are superimposed on general bauplans that sometimes have a greater weight in the general distribution in morphospace. Just as an example, snakes and other limbless taxa are plotted together, but resolution is not good enough to separate fossorial vs marine vs surface dwelling taxa inside this main cluster. The truth is that we are looking for broad-scale diversification events – e.g., we are more interested in the timing of the event that led from a limbed lizard-like squamate to the earliest limbless morphotype or from a lizard-like headed form to a snake-like headed form, than to subsequent changes once a bauplan is achieved that might have occurred multiple times, like adaptations to fossoriality or to a marine environment inside a particular group. A more detailed study of particular groups might be interesting in many ways, but it is unlikely to provide information regarding the initial diversification of squamates, which is the focus of our paper.

Although we agree with the reviewers that the term ecology would be setting perspectives not reached in the text, the term “ecology” is not used in the paper title or any subtitles.

(2) While we believe that the authors have done sufficient due diligence to rule out potential confounding factors, it would still be worthwhile to be more explicit about one limitation, namely the issue of missing data and inapplicable characters. For instance, this paper could be brought up:

Gerber, S. (2019). Use and misuse of discrete character data for morphospace and disparity analyses. Palaeontology, 62(2), 305-319.

We are not convinced that a broad-scale expunging of "inapplicable" characters (say) is the best way to solve the issue (one thinks of limb character in limbless forms), but a more explicit acknowledgment of this limitation would make the paper more compelling.

We agree on the importance of the limitation of missing and inapplicable data, which on the other hand is a recurrent problem identified in squamate phylogenetic data matrices (also for conducting phylogenetic analyses in particular). According to Gerber (2019), the inclusion of inapplicable characters is not advisable, but we concur with the reviewers that simply removing the inapplicable characters does not seem to be the solution; in fact, an alternative approach circumventing this problem still needs to be devised. To test the impact of missing data we performed an additional analysis where we removed all characters that presented a percentage of missing data above 40% (note that because missing and unscorable entries were both assigned a question mark in the original matrix, they cannot be treated differently). In any case, it seems that other works (e.g. Martínez et al., 2021) ended up considering both missing and inapplicable characters as ? in the morphospace analyses, so a separation of ? and – states does not seem to be necessary at this stage. The procedure described above reduced the matrix from its original 836 characters to just 184, but contained a total amount of missing data just below the 25% of missing data considered safe for interpretation. The morphospace resulting is different regarding the position of specific taxa, but still recovers the same main morphospace groups in the same positions as in the full analysis. This confirms that our main results are reliable despite the huge amount of missing data. We suspect (but have not investigated in detail yet) that those characters with a huge amount of missing data were barely contributing to the first axes, explaining the weak effect of their removal in morphospace for these main axes. We have added a comment on the limitation as discussed by Gerber (2019) where needed but we have not included this additional analysis in the results (we are open to doing so, however, if the reviewers or editors find it necessary). Note that one of the performed analyses in the original manuscript already focussed only on osteological characters, which had a similar effect of boosting completeness (mainly among fossils). Another potential problem, the express removal of autapomorphies and/or the supposed lack of collection of autapomorphies seems not to apply to the current matrix because multiple autapomorphies were identified in Conrad (2018), either in the form of autapomorphic binary characters or one or more states in multistate characters, and they were not deactivated. Regarding the true incompleteness of fossils, this cannot be avoided without deleting incomplete fossils. We are reluctant to do so because we think they still are the only source of true (not inferred) morphological information in specific points of the fossil record.

(3) One weakness is the issue of completeness, which has been raised prominently in recent years. The authors approach this by way of examining Cramér coefficients and throwing out characters that are generally missing for fossils (like those of integument). While this is not unreasonable, and results in disparity plots that are broadly comparable to analyses in which these characters are present, there remain other characters (especially those related to the limbs), this problem was less thoroughly addressed than others.

Again, we cannot disagree with the reviewer regarding the importance of completeness and unscorable characters (which effectively act as missing data). However, true incompleteness is inherent in the use of fossils, and removing those characters that concentrate the greatest amount of missing data barely affects results. The percentage of missing data across the entire matrix could be reduced by removing fossils, but even when they are incomplete, fossils provide information regarding which sets of character states were present at specific points in time. This is precisely why they were included originally, and removing them only seems advisable when they are revealed as unstable taxa that are recovered in different points of the tree (although this situation was considered by using multiple randomly selected most parsimonious trees). In the absence of fossils, or of most fossils, our only information on trait distribution in the past would be based on ancestral states estimates, and actual fossil data has been shown repeatedly to be preferable to ancestral data inferred only from modern taxa.

For each variable we want to consider, we are forced to repeat every other analysis performed in order to confirm that any difference in the results is related to this variable in particular. Considering that the time elapsed for some of the performed analyses can reach a few weeks (because the morphological matrix is very large), there is a limit to the number of options one can reasonably check. We consider the issue of incompleteness (via partially preserved taxa in the case of fossils, or unscorable characters in the case of e.g. snakes) a valid concern and we have added mentions to the problem in the text, but we feel this cannot be easily addressed at this stage. As mentioned above we performed, however, a complementary analysis where we deleted all characters that showed more than 40% of missing entries. This results in a strongly reduced matrix that, when analysed, yields very similar results of morphospace distribution. This is interesting in two ways: (1) it seems that those characters with great amounts of missing data contribute little to distributions in morpohospace along the main axes and (2) characters related to the postcranium concentrate the vast majority of characters removed because they are unknown in many fossils, and they are unscorable in limbless taxa. This result is also interesting in showing that the separation of limbless taxa (e.g. snakes and amphisbaenians) is strongly supported by skull characters too, probably due to strong skull adaptations of these taxa to fossoriality and specialised feeding, with independence of the lack of limbs.

4) One of the striking aspects of the results is that the interpreted expansion of morphospace in both the Jurassic and Cretaceous occurs before the interpreted spike in evolutionary rates (not exactly "roughly coincident" as stated). This seemingly counterintuitive result is not addressed.

Actually, the new figure 4 shows that the Jurassic peak in morphospace and evolutionary rates are coincident (note that the disparity peak is placed in the midpoint of a 24-Myr long bin). We admit that the Cretaceous peak in disparity occurs earlier than the Cretaceous peak in evolutionary rates, but the manuscript focussed on discussing the importance of the Jurassic peaks. It is worth noting that the Jurassic peak in evolutionary rates in the late Jurassic is indeed as coincident as it can be according to the disparity metric binning, which shows its peak in the time bin spanning from 168 to 144 Ma. This bin spans a very minor portion of the latest Middle Jurassic (the latest ~ 5 Myr), the entire Late Jurassic (~18 Myr), and the earliest Early Cretaceous (less than 2 Myr). The middle point of this span, which is the one used to graphically represent the value obtained, is 156 Ma, thus well inside the Late Jurassic. Morphospace occupation represents Middle and late Jurassic taxa together, which results in a representation that suggests that morphospace occupation occurred earlier than the disparity or evolutionary rates peak. Note, however, that in the case of squamates, all sampled taxa are Late Jurassic in age because Middle Jurassic taxa were not sampled due to incompleteness of the fossils. One might expect that the inclusion of Middle Jurassic taxa would move all peaks towards the Middle Jurassic but, given the limitations of the fossil record and forced incongruences in time binning, we regard the results as congruent enough to infer that they all hint to the same event.

We acknowledge that the situation is not ideal, but all these issues are common to all such deep-time studies and easy solutions are not available to remove such uncertainties. Because the requirements for each metric and plot are different, we decided to maximise resolution rather than consistency across bins (at the cost of lower resolution). Moreover, and as expressed by the reviewers in other parts of the review, there is no need for evolutionary rates and disparity to be tightly related, although if they are (and we consider that this is a real possibility), they are more likely to be related to a specific event.

5) As someone who has performed principal coordinate analyses on lepidosaur morphological datasets, one reviewer can confirm that the first few dimensions invariably separate limb-bearing from limbless forms, as well as rhynchocephalians from squamates (in the case of this analysis, anguimorphs seem to also contribute to this axis). This is exactly what is found in Figure 1b, but also previously in Watanabe et al., (2019) and Simoes et al., (2020).

Yes, and this was already acknowledged in the text. Also note that we have shown that even if the first dimensions separate limb-bearing from limbless forms, this is not necessarily related to the presence vs absence of limbs, but possibly to strong modifications of the skull, explaining why Watanabe et al., (2019) were able to get this separation with a dataset comprising only skulls (without postcranial information).

As such, separating these four groups (rhynchocephalians, snakes, anguimorphs and "generalized lizards", as defined in lines 95-104 by the authors) explains a large proportion of total variance, which loads prominently on the first axes but is likely to contribute to some extent to others. Given that the matrix being used was "designed to identify the relationships among the major squamate clades" (Conrad 2018), it is no surprise that the same pattern emerges in this study. The "Jurassic morphospace expansion" (JMA) proposed here is nothing more than a period of time when a morphospace occupied by only one of these groups (rhynchocepalinas) gives way to one occupied by all four.

This result is not novel (i.e., does not depend on any discovery made in this study), it is just presented here in a novel and quantitative way.

Precisely, this is acknowledged when Evans (2003) is mentioned. However, going from the crude interpretation of the fossil record and phylogeny to quantitative analyses and finding support for this overlooked idea (most works tend to talk about the Cretaceous, not the Jurassic, as the key period for squamate diversification) is important.

In fact it is already stated in some of the first lines of the introduction that these morpho-groups are known to already be present in the Middle and Late Jurassic (lines 47-49), and that a similar pattern had already been mentioned in the literature (line 238). Furthermore, whether this does represent a true event of expansion or an artifact of the fossil record is also unclear given the poor fossil record of the clade during this period, as noted by the authors. Furthermore, it ultimately depends on many methodological decisions not entirely validated (such as temporal binning, choice of distance metric, etc), as well as character and taxon sampling decisions that were not optimized for the purpose of macroevolutionary inference, and might therefore supported biased patterns, not just during the Jurassic but also later on.

We understand and, to a certain degree share, the concerns of the reviewers regarding the potential impact on methodological decisions. We have tried to take all these variables into account by repeating analyses according to different topologies, temporal binnings, distance metrics, dating methods, etc. We ultimately decided to apply the temporal binnings that optimised resolution for each analysis, even if in the process we sacrificed some consistency when plotting the results. For disparity, we chose two widely used metrics: the weighted mean pairwise distance (WMPD), which is a pre-ordination matrix; and the sum of variances (SoV) which is a post-ordination matrix, in order to assess the effect of choosing one or the other type. Regarding character and taxon sampling, most morphological data matrices used in these kinds of studies were initially built for conducting phylogenetic analyses. The same observation would apply to other macroevolutionary studies, including the mentioned work of Simões et al., (2020) but also Martínez et al., (2021), and many others concerning other groups. Regarding taxon sampling, we also took this into account when we chose the Conrad (2018) matrix, which is the most complete in terms of both extant taxa and fossils. We refer the reviewers to Smith et al., (2021), who state that using morphological matrices intended for inferring phylogeny as the source data for disparity analyses is acceptable as long as the tree is balanced in terms of sampling. We think that the matrix used here presents a very good balance between fossils and extant taxa, as well as across different clades and through time. We made some initial runs of the analyses with other smaller matrices (Gauthier et al., 2012; Simões et al., 2018), but we decided not to use them because the poor sample of fossils outside the late Cretaceous apparently biased results. We would expect that increasing the number of Mesozoic taxa by choosing Conrad’s (2018) matrix would contribute to a greater resolution and a more reliable result. The same applies to the number of characters. It would be hard to justify the use of any other matrix containing a lower number of taxa and characters, unless it was implied that a given matrix had issues regarding the quality of the chosen characters.

That being said, and acknowledging that there are not many alternatives to the choices made, we recognise that some patterns might be influenced by the unevenness of the quality of the fossil record. This is already stated in the text, and is an unavoidable effect of the nature of the fossil record. However, using the alternative (ignoring the fossil record and inferring morphologies from extant taxa) introduces, in our opinion, an even greater bias, among other things because control on the exact time when a given morphological character state appears is lost.

Although all the issues reported by the reviewers are potential sources of bias, and this is acknowledged throughout the paper, only those choices consistently leading to a concentration of time divergences in the mid-late Jurassic are potentially biasing patterns towards the result we present in our manuscript. We have shown that applying Bayesian inferences (new results provided in this review, but not incorporated to the manuscript) does not yield strikingly different results in dating the trees (i.e. Bayesian results are inside the range of dates recovered by the other three methods used). A higher amount of divergence times in the Jurassic (compared to other studies) can be also explained by the proportionally larger amount of Jurassic sampled taxa (in contrast to other matrices). Other studies have used a very limited number of Jurassic and Early Cretaceous taxa, and this invariably leads to younger divergence times for key clades– no matter which method is used for dating, removing key taxa, specifically among the oldest ones, is expected to move divergence times towards the present. There is another important observation. All these datings feed upon the current knowledge of the fossil record. We have already expressed our opinion that our sample of Mesozoic squamates is better than that of other studies, but another concern seems to be the quality of the fossil record itself and associated completeness. Note, however that (1) this is unavoidable because the fossil record is invariably incomplete and, more importantly (2) improvement of knowledge of the fossil record is likely to support our results because any lineage newly found in the Jurassic has the potential to move divergence times towards older ages. The opposite would necessarily imply massive reinterpretations of the affinities of many pre-Late Cretaceous forms. Even in the latter case, the morphospace occupation would not change because it is not tied to a given phylogenetic interpretation.

Regarding the mention of the fact that not all methodological decisions have been entirely validated, we commit ourselves to fully explore alternatives to such methodological choices in an ongoing study that will complement the present manuscript (but we provide preliminary results to the reviewers in order to justify our claims regarding the fact that many of the proposedly conflictive factors have a weak influence upon results).

(6) The discovery of this peak in disparity if further linked to high rates of morphological evolution. Note however that disparity and rates are not necessarily linked, and high rates are not necessary to support a JME. However, this result is entirely dependent on the topology and divergence times supported by this study.

Although evolutionary rates as calculated depend on the topology and dating of the trees, it is important to note that differences in topology were indeed considered in two complementary ways. The main results refer to the constrained analysis, where a backbone constraint was used to force the general interrelationships among extant groups, but a second set of results feeds on a second set of trees resulting from an unconstrained analysis, so both schemes have been considered. Moreover, for each of the two schemes, several MPT’s were randomly selected, and the curve of evolutionary rates considers when evolutionary rates are consistently recovered among different topologies and dating of the trees. Regarding the latter, it has been explained above that the three methods give different divergence times for key clades, but the high evolutionary rates in the late Jurassic remain a constant result, except for the mbl method. In the latter, the concentration of nodes close to the present day results in increased evolutionary rates for these associated short branches. It is worth noting, however, that the Jurassic and Cretaceous peaks are still recovered, they just become less significant because of the massive boost reached by branches near the present. Also, we can now state that the use of alternative methods for dating the trees (Bayesian analyses) are consistent with our results in terms of dating (they are intermediate between those of the Hedman and Equal methods), and in terms of which branches have higher rates of evolution. This is in part related to the fact that the branch leading to snakes, which is invariably situated in the Jurassic, has the highest evolutionary rates of the entire tree by far (together with other branches inside this same clade). This same result regarding high evolutionary rates in the branch leading to snakes was identified in Simões et al., (2020) although in their case Jurassic taxa were not sampled, what presumably could have moved branches (and associated high evolutionary rates) towards earlier times.

The methods of inference and calibration employed here are much more simplistic than those of previous studies, and allow for a much poorer integration of uncertainty into subsequent analyses. Results differ crucially from those of all previous studies, including a much faster early diversification of squamates, and much older ages for some crown clades such as snakes, anguimorphs and iguanians.

Older ages than usually found for crown clades are only recovered in the figured tree (Figure 1), dated with the Hedman method. Note, however, that evolutionary rates were calculated for all three methods of dating, which in the other cases (e.g. equal) do not yield significantly older ages for crown clades, and still recover high evolutionary rates for the late Jurassic (see explanation above). Moreover, as explained above, we have now used Bayesian methods for dating the trees, demonstrating (see attached files) that the Bayesian trees are inside the range of ages recovered under other methods. Methods in some other recent studies include extensive genomic data as well, but such data is only available for modern taxa and can only change the overall tree topology, but not the dating of deep divergences, as here.

This leads to a concentration of branches with high levels of morphological change in the mid to late Jurassic, ultimately responsible for the peak in rates found. The methods employed and the way in which they are presented once again do not rule out the possibility that this is an artifact of the way the study was set up, and a major reason why no such increase in rates was discovered by previous studies.

This is assuming that all methods used recover older divergence ages than other studies, which is not the case. As explained above, the ages we get for the equal method are younger than those we got using Bayesian methods upon the current dataset, and they also recover the Jurassic peak. In any case, and even if we were getting slightly older ages than other methods, it is important to note that because they are usually fragmentary, most of the earliest fossils of clades are not included in the matrix, and thus divergence times recovered in all these studies are probably slightly younger than they would be if we included all available fossils. We are considering including these fossils as calibration points in the Bayesian analysis dating, what would allow to use their temporal information without decreasing the resolution of resulting trees.

In any case, we agree with the reviewers in that the recovered high evolutionary rates in the late Jurassic are not necessary to support the observed Jurassic Morphospace Expansion. These are linked in the paper because appear to be roughly contemporary (more on that later). Again, evidence of the fact that we were aware of the importance of properly dating the trees is reflected in our use of three methods (minimum branch length, equal, and Hedman) for this purpose, and comparing their results regarding evolutionary rates. Note that the divergence times vary from one method to the other, but the late Jurassic peak is consistently recovered, except for the mbl method, that we consider is stretching nodes towards the present because of the high proportion of extant taxa. We recognise that since we designed the study alternative methods for dating the trees have appeared and been used. We think that it is unlikely that using any of these new methods would dramatically change the results because the range of recovered datings is already wide and, in fact, we preliminarily demonstrate that using Bayesian inference makes little difference. In our opinion, failure of previous studies in recovering the Late Jurassic peak in evolutionary rates is more related to their poor sampling of Jurassic taxa than to their use of a more sophisticated dating method (see below, however, regarding our decision to provide a statement that our results will be further validated by an ongoing analysis using Bayesian methods).

(7) There is a lack of consistency throughout the manuscript. For example, the JME label applied to each of the four panels of Figure 1 seems to represent four different periods of time. The topology in Figure 1A was calibrated using a different method that the rates shown in Figure 1D, making it difficult to draw conclusions.

We acknowledge a lack of consistency, with two different sources: one is related to the conscious decision to favour resolution across the different measures and plots above consistency. This results in minor disarrangements in the exact time when an observation is represented. Because of this lack of consistency and the inherent error margin resulting of the use of relatively wide bins, we consider that a single event might be the cause of all these patterns. Figure 1 (former figure 1A) no longer shows the JME, but we have represented evolutionary rates and disparity (Figure 4A and B) in a way that shows that both results show a peak that is coincident in time. Regarding the methods used in figures 1A and 1D (now Figure 1 and Figure 3D), there was an error in the label, as they were both calculated under the Hedman method. Note also that results for all methods are reported in the supplementary data.

While a staggering 51 supplementary figures are presented, many of these are unedited and thus highly cryptic, others such as the many ordination plots are unlabeled and thus do not contribute much.

We think that the fact that the figure legends of the supplementary figures are placed in a separate list might make interpretation more difficult, so we are open to moving each legend below the corresponding figure, and also try to add some information that make them more readable. If the reviewers are asking us to label individual points of plots, this is very difficult for 200 points in each plot, but it was done in supplementary figure 8, which can be used as a guide to know the position in other plots (in the case of PCO1-PCO2, in figures 11,12, 14-19). Finally, it must be said that one can easily get any label plotted by slightly modifying the R code and running it again.

Some aspects of this manuscript are very interesting, such as the position of extinct lineages under topologies constrained to show different relationships among major extant clades. However, its main claim is somewhere between trivial (a change from 1 to 4 morpho-groups somewhere in the second half of the Jurassic) and potentially unsubstantiated (a coincident peak in rates).

Although we acknowledge that an early (Jurassic) diversification of squamates had been already pointed out (e.g. Evans, 2003, already mentioned in our text), this claim had never been tested using numerical methods implying disparity of evolutionary rates. Besides this, we challenge the widespread idea that not much occurred in the evolution of squamates until the apparent radiation of crown groups in the mid-Cretaceous. Moreover, an important observation is that any Jurassic fossil newly added to the dataset (because it was formerly not sampled or because it is newly discovered) is likely to increase support of the proposed diversification event by (1) moving divergence dates to older times; (2) proving a given clade was present in the Jurassic and (3) likely increasing fast evolutionary rates for that period. The opposite (this is, demonstrating that the Jurassic event is an artefact of the methodology used) would probably require that many of the Jurassic fossils were reinterpreted as rather basal forms that were not related to the diversification of the crown, and/or using a dating method that gives very young divergence times for crown clades even when Jurassic taxa are included, but divergence times cannot be younger than fossils in the clades concerned! We aim to provide additional evidence for our claims in the complementary manuscript we are preparing. The fact that other recent studies using similar methodologies have failed to recover some of the patterns presented here or have ignored them, shows that it is indeed interesting to present and discuss them. In contrast to Martínez et al., (2021), who recovered stem lepidosaurs overlapping in morphospace with rhynchocephalians for PC1-PC2, in our results the overlap is between stem-lepidosaurs and squamates for both PC1-PC2 and PC3-PC4. This seems to contradict the assumption by Martínez et al., that rhynchocephalians would present a conservative morphology. This supposed plesiomorphic morphology of rhynchocephalians has been previously claimed, but current studies suggest that at least some of the characters that were supposed to support it are controversial or their polarity was simply reversed (e.g. the presence of acrodont teeth or of a complete lower temporal bar).

The current narrative, emphasizing a completely new and unexpected view on squamate diversification might be true, but it rests on many assumptions. A more realistic interpretation should focus on the existence of four morpho-groups within squamates (which could even be tested with this data, see below), and their timing of origin. The manuscript starts by posing the question "why are some groups highly species-rich and others are not?" (line 16). This question is never answered, and in fact the manuscript does not seek to do so or use data that could. But that kind of exaggeration permeates the rest of the manuscript, and the reviewer thinks it will be benefit from a more realistic approach to the data at hand.

We agree that the manuscript does not seek to answer that particular question, which has been removed, but it investigates the timing in which squamates became morphologically diverse, and if this process was fast (fast evolutionary rates) or not. These are important questions that need to be addressed in order to attempt questions with a wider scope.

8) One reviewer thinks the paper should validate some of their methodological choices better, or even better improve upon them. For example, two disparity metrics are employed (one of which does not even show a JME), but no reason is given why these were chosen. Results from other post-ordination metrics should be presented as well, along with some validation of why certain metrics where favored. This is implemented by Guillerme et al., 2020 Eco Evo and employs some of the same programs already used.

We refer the reviewers to Gerber (2019):

“For the construction of disparity curves, I do not recommend carrying out analyses from the morphospace ordination (‘post-ordination disparity’ in Lloyd 2016) but instead obtain disparity measures directly from the distance matrix to minimize the impact of missing data”.

Regarding the two metrics employed, we chose two of the most widely used metrics, one of them is a pre-ordination metric (WMPD) and the other one (SoV) is a post-ordination metric. Gerber (2019) recommended the calculation of a pre-ordination metric to minimize the impact of missing data, so we have added a disclaimer regarding the possibility that the SoV results might be affected by missing data. We are reluctant to calculate additional postordination metrics following Gerber’s (2019) advice.

Furthermore, a large part of the findings relies on time-calibrated analyses that are poorly explained. For example, although 3 methods of post hoc scaling of topologies are explored, only 1 topology is presented as a supplementary figure. Details on the methods themselves are missing, as it is unclear what time information went into the calibration and how it was gathered and treated. Ultimately however, branch lengths scaling seems to be resulting in unlikely ages for many clades. The reviwer believes results need to be repeated using a topology calibrated using molecular data, or at the very least including secondary node constraints based on other published molecular studies.

Regarding the dating of the trees, we selected one of the dated trees (that resulting from the Hedman method) to illustrate the full results. Note, however, that we performed separate runs of evolutionary rate analyses including the trees resulting from the two other methods, and for each method, several randomly selected trees are used and dated multiple times. According to this methodology, the dating choice is actually considered when interpreting the results, and uncertainties related to specific topologies too. We have added now a figure of a randomly selected dated tree where the range of ages for each node according to the three methods is shown.

The reviewers state that there are no details for the methods used. The mbl and the equal method have been widely used (as implemented in the Paleotree package), and we have added a sentence in the text referring the readers to the package documentation. Regarding the Hedman method, we have added a sentence explaining details on the methodology. The time information was gathered from the Paleobiology database as FAD (first appearance datum) and LAD (last appearance datum) for each sampled taxon. This is the standard way to introduce time information in the methods used for time-calibration (mbl, equal and Hedman).

As expressed elsewhere in our answers to the reviewers, we are preparing a manuscript for submission as a Research Advance to eLife that will deal with alternative methodologies (e.g. Zhang and Wang, 2019) for dating the trees that preliminarily support our results. It is possible to use the same matrix (Conrad 2018) and run a Bayesian analysis that reports a single dated (majority rule tree) and the corresponding evolutionary rates. Again, we will calculate results for both a constrained and unconstrained phylogenetic analysis. Moreover, we will use the version of the method that deals with combined morphological and molecular datasets in order to test if the matrix of Gauthier et al., (2012) recovers some signal for the proposed event in the middle-late Jurassic. The latter approach will be interesting in that the possibility that the Jurassic peak in evolutionary rates is also detectable using smaller matrices will be tested. In the meantime, we have included a explicit statement explaining that the relevant conclusions require additional supporting data that will be published in an ongoing study.

One aspect that could also help improve this paper is an exploration of morpho-functional groups within morphospace. Currently, traditional views are simply repeated (see above regarding the four groups of lepidosaurs) and forced onto the results rather than extracted from them.

Although a more elaborate approach was considered in early stages of this study, it has some limitations. The main problem is that half the sampled taxa are fossils, and as such, their ecomorphology can only be inferred. Thus, one of the possibilities would be to assign an ecomorphological group to each of the extant taxa and use the plotted position for fossils to infer their ecomorphology. An alternative approach, to assign inferred ecomorphology to fossils according to literature, renders the argument circular. The strict approach (assigning ecomorphology for those taxa for which it is confidently known, that is, extant taxa) strongly reduces the number of truly sampled taxa (to around half of the original) and relies on the assumption that the dataset is capable of separating discrete ecomorphological groups and then correctly assign fossils to them. We have doubts regarding the viability of the methodology because, at least on some occasions, the phylogenetic and ecomorphological signals are too strongly mixed to provide a clear output at the species level. As an example, the phylogenetic signal is strong enough to hold marine rhynchocephalians and marine squamates (mosasaurs) separate. The great overlap of “generalised” lizards is another example of the limitations of the proposed exploration. Please note, however, that this same limitation has been previously encountered and is likely to derive of the way characters are distributed in lepidosaur morphological matrices. We suspect that the widely claimed high degree of homoplasy among squamates that is apparently behind the lack of correspondence between molecular and morphological phylogenetic analysis might be behind this poor resolution too. Apparently, morpho-functional signal is strong enough to force artificial groups like the limbless taxa in morphological analyses, but the phylogenetic signal is strong enough to recover many of the true clades (according to molecular data) too. Although the exploration of morphospace proposed by the reviewers is definitely interesting, we think it cannot be easily achieved, and we knowingly simplified discussion so the considered morphospace groups could be tracked through time.

What the data at hand could help establish (possibly along other datasets that are said to have been dropped due to poor sampling of Jurassic fossils) is how many morphological clusters of taxa can be found (using k-means or expectation maximization algorithms for example), and when in time do those originate based on the topologies at hand (not just the one used here but those previously published as well).

We thank the reviewers for the suggestion of the use of k-means and expectation maximization algorithms, which is very interesting, and it will be explored in our ongoing follow-up study. As an alternative approach, we have used the R function pamk of the fpc package in order to get statistically supported clusters. We have added a plot to the supplementary figures showing that when the function is set to group points into 2-4 clusters, the groups we get are almost identical to the ones we used in former figure 2B (now Figure 3b). Allowing for a greater number of groups (8 or 16) gives a much larger number of clusters (not shown) that, on the other hand, gradually increase in similarity to the groups observed in the extended data Figure 11 (this is, phylogenetic groups rather than ecomorphological groups).

(9) The earliest fossil records of the crown squamates are the Middle Jurassic. All these records (for some constituent major clades) in the Middle Jurassic should be clearly stated and/or cited in the text.

A mention to the fact that the Middle and Late Jurassic contains the earliest unambiguous squamates, and also members of the crown is stated in lines 53-56. No more detail is provided because these forms are in many cases not included in the dataset, mostly because they are too incompletely preserved and tend to decrease the resolution of the phylogenetic analyses. Middle Jurassic forms are not sampled, but if anything their inclusion would favour our interpretations.

As one reviewer thinks, this is why the authors separated Early Jurassic from Middle-Late Jurassic in their analyses. The criteria for this separation should be discussed. Otherwise, the readers would doubt why you did not separate Early-Middle Jurassic from Late Jurassic. Figure 1a, which shows the earliest record of the crown squamates later than the Middle Jurassic, is misleading in this regard.

Regarding the fossil record of squamates in Figure 1A, it is important to clarify that the represented record is not the global record of the clades, but the recorded range according to the distribution of sampled taxa in the matrix used. This is the reason why squamates are not recorded until the late Jurassic. This was (and is) explicitly stated in the caption of Figure 1a (now Figure 1).

As we acknowledged above, Middle Jurassic squamates do exist, but they are in almost every case too fragmentary to be included (and they were not in the original data matrix). In this sense, it made more sense to join two time bins that actually contained squamates (even if those from the Middle Jurassic are not sampled), than merging the Middle Jurassic with the Early Jurassic. The chosen binning also resulted in a more balanced distribution of time. Considering that our peak in evolutionary rates and disparity occur in the Late Jurassic, the time binning (early+middle and separate late Jurassic) proposed by the reviewers would fit best our interpretations. However, the inclusion of Middle Jurassic fossils would potentially move the event to the Middle Jurassic.

(10) The manuscript should be edited carefully following the author guidelines of eLife. Some of the discussions in the supplementary information text can be moved to the main text to gain a better readership. The Discussion and especially Extended Discussion would benefit from some further subdivision.

We have added some more weight to the main discussion, and have placed some subdivisions, mainly in the extended discussion.

Some of the supplementary figures (Extended Data Figures 1-10), which are not well demonstrated, can be deleted or replaced with a batch file.

We have replaced Extended Figure 1 with a figure of a randomly selected MPT dated with the Hedman method and showing the ranges of the ages according to all three methods. We are not sure what the reviewers mean with the expression “which are not well demonstrated”, so we retain all supplementary figures until the reviewers clarify which can safely be deleted.

Overall, the supplementary figures should be re-organized to have a uniform format. The authors should put more effort into the supplementary figures, many of which are unintelligible.

We have improved the supplementary figure captions where possible to introduce information that helps in interpreting them. Considering that these are just supplementary figures, and that are the result of the use of multiple packages in R, providing a uniform format is beyond something that can be easily done. However, we hope that the fact that each figure will appear together with its own caption in the published version will help in its interpretation.

Other claims that need to be revisited:

(1) Line 132 – "These evolutionary rates are not dependent on phylogeny": They are very dependent on phylogeny. This analysis shows they are not dependent on topology (molecular vs. morphological) but both of those are likely to share very similar branch lengths for all branches in common, making the result still strongly dependent on phylogeny.

We revise this. We probably should have used the term topology here (we have changed the text now). What we were trying to express was that discrepancies in topology (both at the level of large groups like we see in the constrained vs unconstrained results, and at the species level representing minor changes in topology) were considered by including multiple topologies in the analyses.

(2) Lines 172-173 – "The second event, which fits well with the timing of the KTR, would be coincident with the radiation of the constituent crown groups of Squamata": Figure S1, which is the only time-calibrated topology shown in its entirety, shows crown snakes, anguimorphs, iguanians, lacertoids and skinks originating in the Jurassic, and radiating well before the KTR.

Revised. We are trying to convey a distinction here between the radiation of the crown group Squamata into their corresponding main groups (occurring in the middle-late Jurassic, when the first representatives are found) and the observed main radiation of each of these groups, which seems to be delayed until the KTR. Regarding the latter, we have argued that this is likely an artefact of the fossil record, that suddenly improves in the middle-late Cretaceous.

(3) Lines 189-191 – "The only observable changes from then on correspond to a slight expansion of the edges of the occupied morphospace, and a notable increase in the density of points filling this morphospace": This is entirely contingent on analyses using a matrix that was tuned to capture the deep phylogeny of the clade, and is something that should be acknowledged.

Revised. The problems of using matrices originally built for phylogenetic analysis have been discussed above. However, we want to note that the same main morphotypes could have appeared much later (e.g. in the Late Cretaceous), and they would not be supporting the JME, but the KTR.

(4) Line 271 – "basal position of Iguania": This and other instances of the use of the word basal should be corrected, living groups cannot be basal to other living groups.

We have replaced this expression by “sister-group relationship to the rest of squamates” in this case, and the appropriate expression in the rest of cases (in the supplementary text).

(5) The reviewer has not been able to find which dimensions were considered to estimate the disparity curve shown in Figure 1C. Please add this information to the caption.

Disparity in former Figure 1C (now Figure 4a) corresponds to the Sum of Variances and thus was calculated from the ordinated data, considering all axes data. We have added a statement in the text.

(6) The reviewer would consider striking some comments from the body of the paper that seem more like personal annotations ("We expected that….").

We have modified these comments following the reviewer advice.

(7) L.24: "the Middle and Late Jurassic (174-145 Ma)" is clearer than "the second half of the Jurassic (170-145 Ma)" to the readers, and also consistent with that in Line 48.

We use “the Middle and Late Jurassic” consistently now.

(8) L.82: Figure 1b-d as a separate figure?

Following the reviewer’s advice, and in order to improve the readability of the figures, we provide now Figures 1b and 1c,d as separate figures (new figures 2 and 4, respectively), and former Figure 1a has its own figure (new figure 1) now.

(9) L.90: For "pan-Squamata", do you mean the total group Squamata?

Yes, we have replaced the term (although pan-Squamata has been formally erected, see Gauthier et al., 2020, in Phylonyms, a companion to the PhyloCode)

(10) L.108, 171, 212: "mid-late Jurassic" Ambiguous here、Middle and Late Jurassic?

We use Middle to Late Jurassic now.

(11) L.145：'stem lepidosaurs' for 'ancestral lepidosaurs'

We have replaced ancestral by stem.

(12) L. 151: Extended Data Figure 11 as a main-text figure?

This could be easily done, and we are open to it if the reviewers and editors consider it necessary. However, we wanted to point out that this plot is redundant with former Figure 2B (now Figure 3B). The morphospace in both plots is the same, the differences are that in the main figure (former Figure 2B) the phylogenetic branches are shown and only the main clusters are labelled. In Extended data figure 11 we provided hulls, but the points are the same as in former Figure 2B (new Figure 3B)

(13) L.172, 175, 'late', upper case.

Done.

(14) Extended Data Figures 23, 24, 43-47, 'early', 'middle', 'late', upper cases.

Done.

Reviewer #2:

Weaknesses:

The discussion of ecology (lines 95 and following) is oversimplified (as the authors recognise), perhaps inevitably so.

We certainly recognise that a more detailed discussion on the ecology of groups appears as one of the expected outputs. We face, however, two main limitations here. One is that, even when you look at total morphospace, it is complicated to interpret ecology beyond the main four recognised groups. This is a recurrent problem in ecomorphological studies of squamates, and seems to be independent of the type of dataset (geometric morphometrics vs discrete morphological characters) or differences in sampling. In our opinion, the lack of resolution among ecological clusters is related to the fact that the latter are superimposed on general bauplans that sometimes have a greater weight in the general distribution. We plotted the ecology of taxa in the morphological plot (not shown) in order to see if we were able to make additional interpretations. The main problem is that there are local adaptations to different ecological roles inside of each clade, resulting in an overlap of ecologies that hinder interpretation. Note, however, that we do not regard this as a limitation of the methodology itself, but rather a problem derived from the complexity of the diversification patterns among squamates. This could be regarded as a result in itself and discussed further, but we decided to focus on surveying the timing of broad-scale diversification events, the main of which would have occurred during the Jurassic.

One weakness is the issue of completeness, which has been raised prominently in recent years. The authors approach this by way of examining Cramér coefficients and throwing out characters that are generally missing for fossils (like those of integument). While this is not unreasonable, and results in disparity plots that are broadly comparable to analyses in which these characters are present, there remain other characters (especially those related to the limbs), this problem was less thoroughly addressed than others.

We acknowledge that inapplicable characters are a source of uncertainty in the data, but there is not much that can be done to fix this issue. According to Gerber (2019), the inclusion of inapplicable characters is not advisable, but we concur with the reviewers in that simply removing the inapplicable characters does not seem to be the solution. Moreover, some of the performed analyses were designed to minimise the influence of incompleteness (e.g. the calculation of a pre-ordination metric like the WMPD). We have added explanations regarding the limitations of the methodology used (or of the matrix itself) where corresponding in the text, and we will explore this and other limitations in a subsequent study. We also run an additional analysis where completeness was raised by deleting characters containing more than 40% of missing entries. This resulted in a much smaller dataset, but results of morphospace are completely comparable to those of the full dataset.

One of the striking aspects of the results is that the interpreted expansion of morphospace in both the Jurassic and Cretaceous occurs before the interpreted spike in evolutionary rates (not exactly "roughly coincident" as stated). This seemingly counterintuitive result is not addressed.

We have provided additional explanations for this apparent contradiction in the text and in the detailed responses to reviewers. In the case of comparing the measure of disparity (SoV) and evolutionary rates, the peaks are actually coincident, as it is shown in the new Figure 4 (note that the layout of Figure 4A has been modified for consistency with the time binning of Figure 4B).

Reviewer #3:

Bolet et al., explore the evolutionary history of lepidosaur morphology. Using a previously published dataset sampling a wide range of living and extinct lineages, they analyze both morphospace occupation and rates of evolution through time. By combining a series of phylogenetic and macroevolutionary methods, they reveal an event of morphospace expansion (coupled with increased evolutionary rates) during the second half of the Jurassic. This contrasts with previous studies, which had recovered a predominant signal of diversification and ecological innovation during the mid-Cretaceous, coinciding with the radiation of other major animal and plant lineages and a large-scale restructuring of terrestrial ecosystems. While their conclusion is novel, the manuscript is not based on any new morphological or paleontological data, and employs relatively simple methodologies that can potentially determine the outcome. My main concerns are the following:

(1) As someone who has performed principal coordinate analyses on lepidosaur morphological datasets, I can confirm that the first few dimensions invariably separate limb-bearing from limbless forms, as well as rhynchocephalians from squamates (in the case of this analysis, anguimorphs seem to also contribute to this axis). This is exactly what is found in Figure 1b, but also previously in Watanabe et al., (2019) and Simoes et al., (2020). As such, separating these four groups (rhynchocephalians, snakes, anguimorphs and "generalized lizards", as defined in lines 95-104 by the authors) explains a large proportion of total variance, which loads prominently on the first axes but is likely to contribute to some extent to others. Given that the matrix being used was "designed to identify the relationships among the major squamate clades" (Conrad 2018), it is no surprise that the same pattern emerges in this study. The "Jurassic morphospace expansion" (JMA) proposed here is nothing more than a period of time when a morphospace occupied by only one of these groups (rhynchocepalinas) gives way to one occupied by all four. This result is not novel (i.e., does not depend on any discovery made in this study), it is just presented here in a novel and quantitative way. In fact it is already stated in some of the first lines of the introduction that these morpho-groups are known to already be present in the Middle and Late Jurassic (lines 47-49), and that a similar pattern had already been mentioned in the literature (line 238).

We agree that the results regarding this point are not completely novel (and this is made clear in the text), but our study confirms previous assumptions based on the crude interpretation of the fossil record and, more importantly, it contradicts recent studies using a similar approach that recovered a much later diversification of squamates. It thus contributes to an unsettled debate.

Furthermore, whether this does represent a true event of expansion or an artifact of the fossil record is also unclear given the poor fossil record of the clade during this period, as noted by the authors.

We agree with the reviewer that the fossil record is patchy and, as such, it is just providing partial information. To what degree this is biasing results it is unknown, but we want to note that our choice of morphological matrix and methods reflect our concerns regarding the fossil record issue by (1) including the largest number of fossils available (including many of the late Jurassic forms lacking in other data matrices) and (2) by limiting our conclusions to the information extracted from the data. Moreover, the poor fossil record in the Jurassic and the lack of middle Jurassic sampled squamates (due to incompleteness) are, if anything, weakening the signal for an early burst of diversification of squamates. Adding middle Jurassic taxa would move divergence times towards the past, and thus reinforce our claim of an early event in squamate diversification. It is hard to see, however, how an improved fossil record (or sampling) could move divergence times and associated evolutionary rates towards the mid Cretaceous. Moreover, if the poor fossil record is affecting results, the use of a lower amount of fossils in other studies should be even more problematic. To sum up (1) if the poor fossil record is a problem, it is even more prominent for other studies including a low number of fossils; (2) it is more likely to keep finding older and older representatives of the main groups of squamates than the opposite situation, this is, that we reinterpreted the oldest current examples of such groups in what would be the only way to move the origin of clades towards the present. In our opinion, the poor Jurassic fossil record is more likely to be underrepresenting the importance of Jurassic taxa in the establishment of lizard morphospace limits. Conversely, the enhanced fossil record of the Middle-Late Cretaceous could be artificially inflating measures of disparity, evolutionary rates or diversity for that period. According to this, we think it is worth considering the possibility that the Jurassic signal we get in our manuscript represents a true event.

Furthermore, it ultimately depends on many methodological decisions not entirely validated (such as temporal binning, choice of distance metric, etc), as well as character and taxon sampling decisions that were not optimized for the purpose of macroevolutionary inference, and might therefore supported biased patterns, not just during the Jurassic but also later on.

We agree that any methodological decision can have an impact on the results, and we have tried to provide many ways of auditing the effect of one choice as possible by trying multiple options. We have controlled the effect of topology by including multiple topologies of two different schemes: one is the result of a pure morphological phylogenetic analysis; and one is the result of a constrained analysis according to molecular studies; we have provided two measures of disparity; we have calculated measures of disparity and evolutionary rates according to different time binnings; we calculated evolutionary rates for all lepidosaurs, and for squamates and rhynchocephalians (plus stem-lepidosaurs) separately. The time spent in some of the analyses is measured in weeks, and it is thus prohibitive to keep adding variables. We aim at testing the robustness of our results for many of the points raised by the reviewers in a subsequent paper using complementary methodologies, including the use of tip dating to infer both divergence times and evolutionary rates while accounting for their uncertainties in a coherent Bayesian statistical framework.

(2) The discovery of this peak in disparity if further linked to high rates of morphological evolution. Note however that disparity and rates are not necessarily linked, and high rates are not necessary to support a JME. However, this result is entirely dependent on the topology and divergence times supported by this study. The methods of inference and calibration employed here are much more simplistic than those of previous studies, and allow for a much poorer integration of uncertainty into subsequent analyses. Results differ crucially from those of all previous studies, including a much faster early diversification of squamates, and much older ages for some crown clades such as snakes, anguimorphs and iguanians. This leads to a concentration of branches with high levels of morphological change in the mid to late Jurassic, ultimately responsible for the peak in rates found. The methods employed and the way in which they are presented once again do not rule out the possibility that this is an artifact of the way the study was set up, and a major reason why no such increase in rates was discovered by previous studies.

We are working on applying a more integrated analysis on the same matrix in order to test the robustness of results. We want to note, regarding the fact that other studies have not recovered fast rates in the Jurassic, that it is unlikely that you can get a concentration of branches with high levels in that period if most of the fossils included are late Cretaceous in age (with a minimum sampling of Early Cretaceous and Jurassic forms). No matter which approach for dating the trees you are using, results will be highly affected by a poor sampling in a key period. A Late Jurassic peak in evolutionary rates is recovered in all the three dating methods used, and only in the third (mbl) the peak appears as low in comparison to the extremely high evolutionary rates recovered in branches concentrated close to the present(we have provided details on why this method can be producing biased results in our responses to the reviewers).

(3) There is a lack of consistency throughout the manuscript. For example, the JME label applied to each of the four panels of Figure 1 seems to represent four different periods of time. The topology in Figure 1A was calibrated using a different method that the rates shown in Figure 1D, making it difficult to draw conclusions. While a staggering 51 supplementary figures are presented, many of these are unedited and thus highly cryptic, others such as the many ordination plots are unlabeled and thus do not contribute much.

Regarding the different method used in figures 1A and 1D, in fact there was an error in the caption of former figure 1D. The method used for that figure was, indeed, the Hedman method and, as such, it is congruent with that of figure 1A. Lack of consistency between Figure 1C and 1D was the result of differences in the plot layout, that have been corrected now (new Figures 4A and 4B).

Some aspects of this manuscript are very interesting, such as the position of extinct lineages under topologies constrained to show different relationships among major extant clades. However, its main claim is somewhere between trivial (a change from 1 to 4 morpho-groups somewhere in the second half of the Jurassic) and potentially unsubstantiated (a coincident peak in rates).

Although it might be true that the identification of the JME might not be surprising to some (based on a crude interpretation of the fossil record), validating this view is not a trivial point. Moreover, many of the claims of numerical approaches recently applied in other works have pointed to the KTR as the main event in the diversification of squamates. Challenging this view is important.

Regarding the initial statement that reads: [the manuscript] “employs relatively simple methodologies that can potentially determine the outcome”, before starting the process of reviewing our manuscript, we followed the reviewers advice and explored the use of other approaches like the Bayesian tip dating under relaxed morphological clocks in order to calculate divergence times and evolutionary rates while accounting for their uncertainties (as implemented by e.g. Zhang and Wang, 2019). One of the main issues we encountered is that, if we report evolutionary rates for the majority rule tree, a great polytomy at the base of Squamata precludes any interpretation of the results. This is critical because this portion of the tree is where we expect to get the high evolutionary rates, at least according to our results from other methods. There is no simple solution for this, other than changing the sampling (by adding or removing taxa) in the hope that we get a more resolved topology. This problem of poor resolution of the tree was circumvented in our methodology by using a subsample of randomly selected most parsimonious trees as the source for evolutionary rates instead of a consensus tree. The reported figures of evolutionary rates of our manuscript account for this variation by showing when high or low evolutionary rates are consistently recovered across all topologies. We are aiming at investigating how the points raised by the reviewers regarding methodology affect results in a forthcoming paper. We are in position, however, to say that this alternative methodology provides divergence times that are within those retrieved by the methodologies used in our manuscript (more specifically, between those of the Hedman and the equal methods). We have shared some of these provisional results with the reviewers.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

See comments from Reviewer #3.

Reviewer #3:

I congratulate the authors on the effort put towards revising their manuscript, and after carefully reading their response and new version of the manuscript I have just a handful of minor suggestions that I would like to see implemented, or more strongly stated, in the manuscript before the final version is published.

(1) Editing some of the supplementary files shown would not take that much time at all, and it would make them much more readable than they are now. Reducing the size of the label text in the trees now shown as supplementary figures to figure 1, and color coding them by clade, does not represent too much work. Same holds true for color coding by clade the dots in the many supplementary PCO plots.

(2) The authors state in their response that "Methods in some other recent studies include extensive genomic data as well, but such data is only available for modern taxa and can only change the overall tree topology, but not the dating of deep divergences, as here." Genomic data can most definitely change a lot more than just topology, and in fact do modify the dating of deep divergences. The methods implemented currently for time-scaling make a number of simplified decisions regarding the calculation of branch lengths that come from the fact that they do not incorporate any quantifiable (morphological or molecular) sense of the degree of divergence occurred between successive nodes. Methods directly employing genomic data have access to this amount of divergence, and can therefore not only incorporate it into the analysis but also eschew some of the simplistic assumptions of methods such as equal or mbl (which remains true regardless of the latter being very common in the literature). It is wrong of the authors to assume that time calibration relies exclusively on the use of the right fossils (although this is certainly crucial), and there are good reasons why approaches incorporating genomic data to divergence time estimation are unambiguously better than those that don't, even when molecular data is available only for extant terminals. I strongly believe that this caveat needs to be made explicit in the main text of the manuscript, and that it is not a problem that disappears just because several different methods explored agree with each other.

(3) In a similar vein, I would suggest the authors to explicitly state that some results that directly depend on the estimation of time-calibrated branch lengths (such as the calculation of evolutionary rates) necessitate a confirmation with more robust methods of time calibration, such as those that incorporate molecular data alongside the fossil record, as well as employ more realistic models of diversification such as the fossilized birth-death prior.

We thank again the expert reviewer who made additional comments on our reviewed version of the manuscript. We found them highly valuable, and we have followed them by making the following modifications on the manuscript:

1. We have modified the morphospace plots where all tips are labeled and points are now colored according to clades. We have also replaced the complete taxon names by abbreviations, to avoid (or at least minimise) label overlapping.

2. We have modified the consensus trees, now showing colors according to clade

3. We have replaced the corresponding code for morphospace by an updated version including the code for the new morphospace plots.

4. We have added a new source file for the code for plotting the new (coloured) consensus trees.

5. Following the reviewer’s advice, we have added an explicit statement highlighting that some results that directly depend on the estimation of time-calibrated branch lengths (such as the calculation of evolutionary rates) need a confirmation with more robust methods of time calibration, such as those that incorporate molecular data alongside the fossil record, as well as employ more realistic models of diversification such as the fossilized birth-death prior. This way we address the comments 2 and 3 of the reviewer.

6. We have updated the manuscript text by removing the previous track changes, and now contains only the few modifications made according to the changes described above (mainly related to the captions corresponding to modified figures)

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supplementary file 1. Interactive plot of 3D morphospace (PCO1, PCO2, and PCO3), colors according to high taxonomical groups.

elife-66511-supp1.zip^{(1.4MB, zip)}

Supplementary file 2. Interactive plot of 3D morphospace (PCO1, PCO2, and PCO3), colors according to low taxonomical groups.

elife-66511-supp2.zip^{(1.4MB, zip)}

Supplementary file 3. Interactive plot of 3D morphospace (PCO1, PCO2, and PCO3), colors according to time bins.

elife-66511-supp3.zip^{(1.4MB, zip)}

Supplementary file 4. Interactive plot of 3D phylomorphospace (PCO1, PCO2, and PCO3), colors according to high taxonomical groups.

elife-66511-supp4.zip^{(160.3KB, zip)}

Supplementary file 5. Interactive plot of 3D phylomorphospace (PCO1, PCO2, and PCO3), colors according to low taxonomical groups.

elife-66511-supp5.zip^{(160.5KB, zip)}

Supplementary file 6. Interactive plot of squamate 3D morphospace (PCO1, PCO2, and PCO3).

Black spheres correspond to toxicoferans, white spheres correspond to non-toxicoferans.

elife-66511-supp6.zip^{(157.3KB, zip)}

Transparent reporting form

elife-66511-transrepform1.pdf^{(322.9KB, pdf)}

Source code 1. R script and source files for morphospace and disparity analyses.

elife-66511-code1.zip^{(362.5KB, zip)}

Source code 2. R script and source files for evolutionary rate analyses.

elife-66511-code2.zip^{(80.6KB, zip)}

Source code 3. R script and source files for plotting morphospace clusters.

elife-66511-code3.zip^{(136.6KB, zip)}

Source code 4. R script and source files for plotting Crámer values.

elife-66511-code4.zip^{(5.7MB, zip)}

Source code 5. R script and source files for plotting consensus trees.

elife-66511-code5.zip^{(13.5KB, zip)}

Data Availability Statement

[bib1] Adams D, Collyer M, Kaliontzopoulou A. geomorph: an r package for the collection and analysis of geometric morphometric shape data. Methods in Ecology and Evolution. 2019;1:12035. doi: 10.1111/2041-210X.12035. [DOI] [Google Scholar]

[bib2] Bapst DW. Paleotree: an R package for paleontological and phylogenetic analyses of evolution. Methods in Ecology and Evolution. 2012;3:803–807. doi: 10.1111/j.2041-210X.2012.00223.x. [DOI] [Google Scholar]

[bib3] Benson RBJ, Campione NE, Carrano MT, Mannion PD, Sullivan C, Upchurch P, Evans DC. Rates of dinosaur body mass evolution indicate 170 million years of sustained ecological innovation on the avian stem lineage. PLOS Biology. 2014;12:e1001853. doi: 10.1371/journal.pbio.1001853. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib4] Benton MJ. The origins of modern biodiversity on land. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences. 2010;365:3667–3679. doi: 10.1098/rstb.2010.0269. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib5] Bolet A, Evans SE. New material of the enigmatic Scandensia, an Early Cretaceous lizard from the Iberian Peninsula. Spec. Pap. Palaeontol. 2011;86:99–108. doi: 10.1111/J.1475-4983.2011.01074.x. [DOI] [Google Scholar]

[bib6] Bolet A, Evans SE. A tiny lizard (Lepidosauria, Squamata) from the Lower Cretaceous of Spain. Palaeontology. 2012;55:491–500. doi: 10.1111/j.1475-4983.2012.01145.x. [DOI] [Google Scholar]

[bib7] Brusatte SL, Benton MJ, Ruta M, Lloyd GT. Superiority, competition, and opportunism in the evolutionary radiation of dinosaurs. Science (New York, N.Y.) 2008;321:1485–1488. doi: 10.1126/science.1161833. [DOI] [PubMed] [Google Scholar]

[bib8] Brusatte SL, O’Connor JK, Jarvis ED. The Origin and Diversification of Birds. Current Biology. 2015;25:R888–R898. doi: 10.1016/j.cub.2015.08.003. [DOI] [PubMed] [Google Scholar]

[bib9] Burbrink FT, Grazziotin FG, Pyron RA, Cundall D, Donnellan S, Irish F, Keogh JS, Kraus F, Murphy RW, Noonan B, Raxworthy CJ, Ruane S, Lemmon AR, Lemmon EM, Zaher H. Interrogating Genomic-Scale Data for Squamata (Lizards, Snakes, and Amphisbaenians) Shows no Support for Key Traditional Morphological Relationships. Systematic Biology. 2020;69:502–520. doi: 10.1093/sysbio/syz062. [DOI] [PubMed] [Google Scholar]

[bib10] Caldwell MW, Nydam RL, Palci A, Apesteguía S. The oldest known snakes from the Middle Jurassic-Lower Cretaceous provide insights on snake evolution. Nature Communications. 2015;6:5996. doi: 10.1038/ncomms6996. [DOI] [PubMed] [Google Scholar]

[bib11] Cardinal S, Danforth BN. Bees diversified in the age of eudicots. Proceedings. Biological Sciences. 2013;280:20122686. doi: 10.1098/rspb.2012.2686. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib12] Čerňanský A, Bolet A, Müller J, Rage J-C, Augé M, Herrel A. A new exceptionally preserved specimen of Dracaenosaurus (Squamata, Lacertidae) from the Oligocene of France as revealed by micro-computed tomography. Journal of Vertebrate Paleontology. 2017;37:e1384738. doi: 10.1080/02724634.2017.1384738. [DOI] [Google Scholar]

[bib13] Cleary TJ, Benson RBJ, Evans SE, Barrett PM. Lepidosaurian diversity in the Mesozoic-Palaeogene: the potential roles of sampling biases and environmental drivers. Royal Society Open Science. 2018;5:171830. doi: 10.1098/rsos.171830. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib14] Close RA, Friedman M, Lloyd GT, Benson RBJ. Evidence for a Mid-Jurassic Adaptive Radiation in Mammals. Current Biology. 2015;25:2137–2142. doi: 10.1016/j.cub.2015.06.047. [DOI] [PubMed] [Google Scholar]

[bib15] Close RA, Benson RBJ, Alroy J, Behrensmeyer AK, Benito J, Carrano MT, Cleary TJ, Dunne EM, Mannion PD, Uhen MD, Butler RJ. Diversity dynamics of Phanerozoic terrestrial tetrapods at the local-community scale. Nature Ecology & Evolution. 2019;3:590–597. doi: 10.1038/s41559-019-0811-8. [DOI] [PubMed] [Google Scholar]

[bib16] Conrad JL. Phylogeny And Systematics Of Squamata (Reptilia) Based On Morphology. Bulletin of the American Museum of Natural History. 2008;310:1–182. doi: 10.1206/310.1. [DOI] [Google Scholar]

[bib17] Conrad JL, Ast JC, Montanari S, Norell MA. A combined evidence phylogenetic analysis of Anguimorpha (Reptilia: Squamata. Cladistics : The International Journal of the Willi Hennig Society. 2011;27:230–277. doi: 10.1111/j.1096-0031.2010.00330.x. [DOI] [PubMed] [Google Scholar]

[bib18] Conrad JL, Balcarcel AM, Mehling CM. Earliest example of a giant monitor lizard (Varanus, Varanidae, Squamata. PLOS ONE. 2012;7:e41767. doi: 10.1371/journal.pone.0041767. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib19] Conrad JL. A new lizard (Squamata) was the last meal of Compsognathus (Theropoda: Dinosauria) and is A holotype in A holotype. Zoological Journal of the Linnean Society. 2018;183:584–634. doi: 10.1093/zoolinnean/zlx055. [DOI] [Google Scholar]

[bib20] Datta PM, Ray S. Earliest lizard from the late triassic (Carnian) of India. Journal of Vertebrate Paleontology. 2006;26:795–800. doi: 10.1671/0272-4634. [DOI] [Google Scholar]

[bib21] Doyle JA. Integrating Molecular Phylogenetic and Paleobotanical Evidence on Origin of the Flower. International Journal of Plant Sciences. 2008;169:816–843. doi: 10.1086/589887. [DOI] [Google Scholar]

[bib22] Estes R. Sauria Terrestria, Amphisbaenia (Handbuch Der Palaoherpetologie 10A. Stuttgart: Gustav Fischer Verlag; 1983. [Google Scholar]

[bib23] Estes R, Queiroz KD, Gauthier J. In: Phylogenetic Relationships of the Lizard Families. Estes R, Pregill G, editors. Stanford: Stanford University Press; 1988. Phylogenetic relationships within Squamata; pp. 119–281. [Google Scholar]

[bib24] Evans SE. Lizards: evolution, early radiation and biogeography. Proceedings of the 6th symposium on Mesozoic terrestrial ecosystems and Biota. Short pap.1995. [Google Scholar]

[bib25] Evans SE. Crown Group Lizards (Reptilia, Squamata) from the Middle Jurassic of the British Isles. Palaeontographica Abteilung A. 1998;250:123–154. doi: 10.1127/pala/250/1998/123. [DOI] [Google Scholar]

[bib26] Evans SE, Barbadillo LJ. An unusual lizard (Reptilia: Squamata) from the Early Cretaceous of Las Hoyas, Spain. Zoological Journal of the Linnean Society. 1998;124:235–265. doi: 10.1111/j.1096-3642.1998.tb00576.x. [DOI] [Google Scholar]

[bib27] Evans SE, Barbadillo LJ. A short-limbed lizard from the Lower Cretaceous of Spain. Special Papers in Palaeontology. 1999;60:73–85. doi: 10.1111/j.1475-4983.2012.01145. [DOI] [Google Scholar]

[bib28] Evans SE. At the feet of the dinosaurs: the early history and radiation of lizards. Biological Reviews of the Cambridge Philosophical Society. 2003;78:513–551. doi: 10.1017/s1464793103006134. [DOI] [PubMed] [Google Scholar]

[bib29] Evans SE, Raia P, Barbera C. The Lower Cretaceous lizard genus Chometokadmon from Italy. Cretaceous Research. 2006;27:673–683. doi: 10.1016/j.cretres.2006.03.004. [DOI] [Google Scholar]

[bib30] Evans SE. In: The Skull of Lepidosauria. Gans C, Gaunt AS, editors. Society for the Study of Amphibians and Reptiles; 2008. Biology of the Reptilia; pp. 1–347. [Google Scholar]

[bib31] Fry BG, Vidal N, Norman JA, Vonk FJ, Scheib H, Ramjan SFR, Kuruppu S, Fung K, Hedges SB, Richardson MK, Hodgson WC, Ignjatovic V, Summerhayes R, Kochva E. Early evolution of the venom system in lizards and snakes. Nature. 2006;439:584–588. doi: 10.1038/nature04328. [DOI] [PubMed] [Google Scholar]

[bib32] Gauthier JA, Kearney M, Maisano JA, Rieppel O, Behlke ADB. Assembling the Squamate Tree of Life: Perspectives from the Phenotype and the Fossil Record. Bulletin of the Peabody Museum of Natural History. 2012;53:3–308. doi: 10.3374/014.053.0101. [DOI] [Google Scholar]

[bib33] Gerber S, Ruta M. Use and misuse of discrete character data for morphospace and disparity analyses. Palaeontology. 2019;62:305–319. doi: 10.1111/pala.12407. [DOI] [Google Scholar]

[bib34] Goloboff PA, Catalano SA. TNT version 1.5, including a full implementation of phylogenetic morphometrics. Cladistics : The International Journal of the Willi Hennig Society. 2016;32:221–238. doi: 10.1111/cla.12160. [DOI] [PubMed] [Google Scholar]

[bib35] Griffiths EF, Ford DP, Benson RBJ, Evans SE, Ruta M. New information on the Jurassic lepidosauromorph Marmoretta oxoniensis. Papers in Palaeontology. 2021;7:2255–2278. doi: 10.1002/spp2.1400. [DOI] [Google Scholar]

[bib36] Guillerme T, Poisot T. dispRity: A modular R package for measuring disparity. Methods in Ecology and Evolution. 2018;9:1755–1763. doi: 10.1111/2041-210X.13022. [DOI] [Google Scholar]

[bib37] Hedman MM. Constraints on clade ages from fossil outgroups. Paleobiology. 2010;36:16–31. doi: 10.1666/0094-8373-36.1.16. [DOI] [Google Scholar]

[bib38] Herrera-Flores JA, Stubbs TL, Benton MJ, Ruta M. Macroevolutionary patterns in Rhynchocephalia: is the tuatara ( Sphenodon punctatus ) a living fossil? Palaeontology. 2017;60:319–328. doi: 10.1111/pala.12284. [DOI] [Google Scholar]

[bib39] Herrera-Flores JA, Elsler A, Stubbs TL, Benton MJ, Benson R. Slow and fast evolutionary rates in the history of lepidosaurs. Palaeontology. 2021a;65:12579. doi: 10.1111/pala.12579. [DOI] [Google Scholar]

[bib40] Herrera-Flores JA, Stubbs TL, Benton MJ. Ecomorphological diversification of squamates in the Cretaceous. Royal Society Open Science. 2021b;8:201961. doi: 10.1098/rsos.201961. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib41] Howard TJ, Swofford JB, Wagner DL, Sherman S, Lehman GA. Quality of life after bilateral thoracoscopic splanchnicectomy: long-term evaluation in patients with chronic pancreatitis. Journal of Gastrointestinal Surgery. 2002;6:845–852. doi: 10.1016/s1091-255x(02)00123-3. [DOI] [PubMed] [Google Scholar]

[bib42] Hutchinson MN, Skinner A, Lee MSY. Tikiguania and the antiquity of squamate reptiles (lizards and snakes. Biology Letters. 2012;8:665–669. doi: 10.1098/rsbl.2011.1216. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib43] Jones MEH, Anderson CL, Hipsley CA, Müller J, Evans SE, Schoch RR. Integration of molecules and new fossils supports a Triassic origin for Lepidosauria (lizards, snakes, and tuatara. BMC Evolutionary Biology. 2013;13:208. doi: 10.1186/1471-2148-13-208. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib44] Kearney M. The phylogenetic position of Sineoamphisbaena hexatabularis reexamined. Journal of Vertebrate Paleontology. 2003;23:394–403. doi: 10.1671/0272-4634(2003)023[0394:TPPOSH]2.0.CO;2. [DOI] [Google Scholar]

[bib45] Kotrc B, Knoll AH. A morphospace of planktonic marine diatoms. I. Two views of disparity through time. Paleobiology. 2015;41:45–67. doi: 10.1017/pab.2014.4. [DOI] [Google Scholar]

[bib46] Lafuma F, Corfe IJ, Clavel J, Di-Poï N. Multiple evolutionary origins and losses of tooth complexity in squamates. Nature Communications. 2021;12:6001. doi: 10.1038/s41467-021-26285-w. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib47] Laurin M. The evolution of body size, Cope’s rule and the origin of amniotes. Systematic Biology. 2004;53:594–622. doi: 10.1080/10635150490445706. [DOI] [PubMed] [Google Scholar]

[bib48] Lee MSY, Cau A, Naish D, Dyke GJ. Dinosaur evolution. Sustained miniaturization and anatomical innovation in the dinosaurian ancestors of birds. Science (New York, N.Y.) 2014;345:562–566. doi: 10.1126/science.1252243. [DOI] [PubMed] [Google Scholar]

[bib49] Li PP, Gao KQ, Hou LH, Xu X. A gliding lizard from the Early Cretaceous of China. PNAS. 2007;104:5507–5509. doi: 10.1073/pnas.0609552104. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib50] Lloyd GT, Davis KE, Pisani D, Tarver JE, Ruta M, Sakamoto M, Hone DWE, Jennings R, Benton MJ. Dinosaurs and the Cretaceous Terrestrial Revolution. Proceedings. Biological Sciences. 2008;275:2483–2490. doi: 10.1098/rspb.2008.0715. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib51] Lloyd GT. Estimating morphological diversity and tempo with discrete character-taxon matrices: implementation, challenges, progress, and future directions. Biological Journal of the Linnean Society. 2016;118:131–151. doi: 10.1111/bij.12746. [DOI] [Google Scholar]

[bib52] Lloyd GT, Bapst DW, Friedman M, Davis KE. Probabilistic divergence time estimation without branch lengths: dating the origins of dinosaurs, avian flight and crown birds. Biology Letters. 2016;12:20160609. doi: 10.1098/rsbl.2016.0609. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib53] Longrich NR, Bhullar BAS, Gauthier JA. Mass extinction of lizards and snakes at the Cretaceous-Paleogene boundary. PNAS. 2012;109:21396–21401. doi: 10.1073/pnas.1211526110. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib54] Longrich NR, Vinther J, Pyron RA, Pisani D, Gauthier JA. Biogeography of worm lizards (Amphisbaenia) driven by end-Cretaceous mass extinction. Proceedings of the Royal Society B. 2015;282:20143034. doi: 10.1098/rspb.2014.3034. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib55] Losos JB, Hillis DM, Greene HW. Evolution. Who speaks with a forked tongue? Science (New York, N.Y.) 2012;338:1428–1429. doi: 10.1126/science.1232455. [DOI] [PubMed] [Google Scholar]

[bib56] Martínez RN, Simões TR, Sobral G, Apesteguía S. A Triassic stem lepidosaur illuminates the origin of lizard-like reptiles. Nature. 2021;597:235–238. doi: 10.1038/s41586-021-03834-3. [DOI] [PubMed] [Google Scholar]

[bib57] Meredith RW, Janečka JE, Gatesy J, Ryder OA, Fisher CA, Teeling EC, Goodbla A, Eizirik E, Simão TLL, Stadler T, Rabosky DL, Honeycutt RL, Flynn JJ, Ingram CM, Steiner C, Williams TL, Robinson TJ, Burk-Herrick A, Westerman M, Ayoub NA, Springer MS, Murphy WJ. Impacts of the Cretaceous Terrestrial Revolution and KPg extinction on mammal diversification. Science (New York, N.Y.) 2011;334:521–524. doi: 10.1126/science.1211028. [DOI] [PubMed] [Google Scholar]

[bib58] Mongiardino Koch N, Gauthier JA. Noise and biases in genomic data may underlie radically different hypotheses for the position of Iguania within Squamata. PLOS ONE. 2018;13:e0202729. doi: 10.1371/journal.pone.0202729. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib59] Moon BC, Stubbs TL. Early high rates and disparity in the evolution of ichthyosaurs. Communications Biology. 2020;3:68. doi: 10.1038/s42003-020-0779-6. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib60] Mulcahy DG, Noonan BP, Moss T, Townsend TM, Reeder TW, Sites JW, Wiens JJ. Estimating divergence dates and evaluating dating methods using phylogenomic and mitochondrial data in squamate reptiles. Molecular Phylogenetics and Evolution. 2012;65:974–991. doi: 10.1016/j.ympev.2012.08.018. [DOI] [PubMed] [Google Scholar]

[bib61] Müller J, Hipsley CA, Head JJ, Kardjilov N, Hilger A, Wuttke M, Reisz RR. Eocene lizard from Germany reveals amphisbaenian origins. Nature. 2011;473:364–367. doi: 10.1038/nature09919. [DOI] [PubMed] [Google Scholar]

[bib62] Nordén KK, Stubbs TL, Prieto-Márquez A, Benton MJ. Multifaceted disparity approach reveals dinosaur herbivory flourished before the end-Cretaceous mass extinction. Paleobiology. 2018;44:620–637. doi: 10.1017/pab.2018.26. [DOI] [Google Scholar]

[bib63] Pianka ER, Vitt LJ, Pelegrin N, Fitzgerald DB, Winemiller KO. Toward a Periodic Table of Niches, or Exploring the Lizard Niche Hypervolume. The American Naturalist. 2017;190:601–616. doi: 10.1086/693781. [DOI] [PubMed] [Google Scholar]

[bib64] Puttick MN, Thomas GH, Benton MJ. High rates of evolution preceded the origin of birds. Evolution; International Journal of Organic Evolution. 2014;68:1497–1510. doi: 10.1111/evo.12363. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib65] Pyron RA, Burbrink FT, Wiens JJ. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evolutionary Biology. 2013;13:1–56. doi: 10.1186/1471-2148-13-93. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib66] Pyron RA. Novel Approaches for Phylogenetic Inference from Morphological Data and Total-Evidence Dating in Squamate Reptiles (Lizards, Snakes, and Amphisbaenians. Systematic Biology. 2017;66:38–56. doi: 10.1093/sysbio/syw068. [DOI] [PubMed] [Google Scholar]

[bib67] R Development Core Team . Vienna, Austria: R Foundation for Statistical Computing; 2013. http://www.R-project.org/ [Google Scholar]

[bib68] Raup DM. Biases in the fossil record of species and genera. Bulletin of Carnegie Museum of Natural History. 1979;13:85–91. [Google Scholar]

[bib69] Revell LJ. phytools: an R package for phylogenetic comparative biology (and other things. Methods in Ecology and Evolution. 2012;3:217–223. doi: 10.1111/j.2041-210X.2011.00169.x. [DOI] [Google Scholar]

[bib70] Reynoso VH. Huehuecuetzpalli mixtecus gen. et sp. nov: a basal squamate (Reptilia) from the Early Cretaceous of Tepexi de Rodríguez, Central México. Philosophical Transactions of the Royal Society of London. Series B. 1998;353:477–500. doi: 10.1098/rstb.1998.0224. [DOI] [Google Scholar]

[bib71] Reynoso VH, Callison G. A new scincomorph lizard from the Early Cretaceous of Puebla, Mexico. Zoological Journal of the Linnean Society. 2000;130:183–212. doi: 10.1111/j.1096-3642.2000.tb01629.x. [DOI] [Google Scholar]

[bib72] Schwenk K. In: Phylogenetic Relationships of the Lizard Families. Estes R, Pregill G, editors. Stanford: Stanford University Press; 1988. Comparative morphology of the lepidosaur tongue and its relevance to squamate phylogeny; pp. 569–597. [Google Scholar]

[bib73] Sepkoski JJ. A factor analytic description of the Phanerozoic marine fossil record. Paleobiology. 2016;7:36–53. doi: 10.1017/S0094837300003778. [DOI] [Google Scholar]

[bib74] Sievert C. ggplot2; 2018. https://plotly-r.com [Google Scholar]

[bib75] Simões TR, Caldwell MW, Tałanda M, Bernardi M, Palci A, Vernygora O, Bernardini F, Mancini L, Nydam RL. The origin of squamates revealed by a Middle Triassic lizard from the Italian Alps. Nature. 2018;557:706–709. doi: 10.1038/s41586-018-0093-3. [DOI] [PubMed] [Google Scholar]

[bib76] Simões TR, Vernygora O, Caldwell MW, Pierce SE. Megaevolutionary dynamics and the timing of evolutionary innovation in reptiles. Nature Communications. 2020;11:3322. doi: 10.1038/s41467-020-17190-9. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib77] Smith KT. Eocene Lizards of the Clade Geiseltaliellus from Messel and Geiseltal, Germany, and the Early Radiation of Iguanidae (Reptilia: Squamata. Bulletin of the Peabody Museum of Natural History. 2009;50:219–306. doi: 10.3374/014.050.0201. [DOI] [Google Scholar]

[bib78] Sullivan C, Wang Y, Hone DWE, Wang Y, Xu X, Zhang F. The vertebrates of the Jurassic Daohugou Biota of northeastern China. Journal of Vertebrate Paleontology. 2014;34:243–280. doi: 10.1080/02724634.2013.787316. [DOI] [Google Scholar]

[bib79] Uetz P, Freed P, Hošek J. The Reptile Database. 2019. [February 19, 2020]. http://www.reptile-database.org

[bib80] Wainwright PC, Alfaro ME, Bolnick DI, Hulsey CD. Many-to-One Mapping of Form to Function: A General Principle in Organismal Design? Integrative and Comparative Biology. 2005;45:256–262. doi: 10.1093/icb/45.2.256. [DOI] [PubMed] [Google Scholar]

[bib81] Wang M, Lloyd GT. Rates of morphological evolution are heterogeneous in Early Cretaceous birds. Proceedings of the Royal Society B. 2016;283:20160214. doi: 10.1098/rspb.2016.0214. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib82] Watanabe A, Fabre AC, Felice RN, Maisano JA, Müller J, Herrel A, Goswami A. Ecomorphological diversification in squamates from conserved pattern of cranial integration. PNAS. 2019;116:14688–14697. doi: 10.1073/pnas.1820967116. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib83] Wickham H. Ggplot2: Elegant Graphics for Data Analysis. Cham: Springer-Verlag; 2016. [DOI] [Google Scholar]

[bib84] Zhang C, Wang M. Bayesian tip dating reveals heterogeneous morphological clocks in Mesozoic birds. Royal Society Open Science. 2019;6:182062. doi: 10.1098/rsos.182062. [DOI] [PMC free article] [PubMed] [Google Scholar]

[bib85] Zheng Y, Wiens JJ. Combining phylogenomic and supermatrix approaches, and a time-calibrated phylogeny for squamate reptiles (lizards and snakes) based on 52 genes and 4162 species. Molecular Phylogenetics and Evolution. 2016;94:537–547. doi: 10.1016/j.ympev.2015.10.009. [DOI] [PubMed] [Google Scholar]

PERMALINK

The Jurassic rise of squamates as supported by lepidosaur disparity and evolutionary rates

Arnau Bolet

Thomas L Stubbs

Jorge A Herrera-Flores

Michael J Benton

Roles

Abstract

Introduction

Results

Figure 1. Lepidosaur— phylogeny, morphospace, disparity, and evolutionary rates.

Figure 1—figure supplement 2. Majority rule consensus tree for the constrained analysis of lepidosaurs.

Figure 1—figure supplement 3. Adams consensus tree for the constrained analysis of lepidosaurs.

Figure 1—figure supplement 4. Strict consensus tree for the constrained analysis of lepidosaurs.

Figure 1—figure supplement 5. Majority rule consensus tree for the unconstrained analysis of lepidosaurs.

Figure 1—figure supplement 6. Adams consensus tree for the unconstrained analysis of lepidosaurs.

Figure 1—figure supplement 7. Strict consensus tree for the unconstrained analysis of lepidosaurs.

Figure 2. Morphospace occupation through time.

Figure 2—figure supplement 1. Morphospace of lepidosaurs through time for PCO1 and PCO2.

Figure 2—figure supplement 2. Morphospace of lepidosaurs through time for PCO1 and PCO2.

Figure 2—figure supplement 3. Morphospace of lepidosaurs through time for PCO3 and PCO4.

Figure 2—figure supplement 4. Morphospace of lepidosaurs through time for PCO3 and PCO4.

Figure 2—figure supplement 5. Morphospace of lepidosaurs through time for PCO5 and PCO6.

Figure 2—figure supplement 6. Morphospace of lepidosaurs through time for PCO5 and PCO6.

Figure 3. Lepidosaur morphospace.

Figure 3—figure supplement 1. Lepidosaur morphospace for PCO1 and PCO2 showing abbreviated taxa labels and colors according to clades (see color legend in Figure 1—figure supplement 2).

Figure 3—figure supplement 2. Lepidosaur morphospace for PCO3 and PCO4 showing abbreviated taxa labels and colors according to clades (see color legend in Figure 1—figure supplement 2 and abbreviations in Figure 3—figure supplement 1).

Figure 3—figure supplement 3. Lepidosaur morphospace for PCO5 and PCO6 showing abbreviated taxa labels and colors according to clades (see color legend in Figure 1—figure supplement 2 and abbreviations in Figure 3—figure supplement 1).

Figure 3—figure supplement 4. Lepidosaur morphospace for PCO1 and PCO2.

Figure 3—figure supplement 5. Lepidosaur phylomorphospace in 2D for PCO1 and PCO2.

Figure 3—figure supplement 6. Plot of morphospace and clusters after applying the pamk function of the fpc package.

Figure 4. Disparity and evolutionary rates through time.

Figure 4—figure supplement 1. Disparity (weighted mean pairwise distance) by high taxonomical group.

Figure 4—figure supplement 2. Disparity (weighted mean pairwise distance) by low taxonomical group.

Figure 4—figure supplement 3. Disparity (weighted mean pairwise distance) through time for lepidosaurs and time scheme 1.

Figure 4—figure supplement 4. Disparity (weighted mean pairwise distance) through time for lepidosaurs and time scheme 2.

Figure 4—figure supplement 5. Evolutionary rates for lepidosaurs according to the constrained phylogeny as represented by five most parsimonious trees (MPT) dated five times using the Hedman method.

Figure 4—figure supplement 6. Evolutionary rates for lepidosaurs according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the equal method.

Figure 4—figure supplement 7. Evolutionary rates for lepidosaurs according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the MBL method.

Figure 4—figure supplement 8. Evolutionary rates for squamates according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the Hedman method.

Figure 4—figure supplement 9. Evolutionary rates for squamates according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the equal method.

Figure 4—figure supplement 10. Evolutionary rates for squamates according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the MBL method.

Figure 4—figure supplement 11. Evolutionary rates for stem lepidosaurs plus rhynchocephalians according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the Hedman method.

Figure 4—figure supplement 12. Evolutionary rates for stem lepidosaurs plus rhynchocephalians according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the equal method.

Figure 4—figure supplement 13. Evolutionary rates for stem lepidosaurs plus rhynchocephalians according to the constrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the MBL method.

Figure 4—figure supplement 14. Evolutionary rates for lepidosaurs according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the Hedman method.

Figure 4—figure supplement 15. Evolutionary rates for lepidosaurs according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the equal method.

Figure 4—figure supplement 16. Evolutionary rates for lepidosaurs according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the MBL method.

Figure 4—figure supplement 17. Evolutionary rates for squamates according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the Hedman method.

Figure 4—figure supplement 18. Evolutionary rates for squamates according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the equal method.

Figure 4—figure supplement 19. Evolutionary rates for squamates according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the MBL method.

Figure 4—figure supplement 20. Evolutionary rates for stem lepidosaurs and rhynchocephalians according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the Hedman method.

Figure 4—figure supplement 21. Evolutionary rates for stem lepidosaurs and rhynchocephalians according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the equal method.

Figure 4—figure supplement 22. Evolutionary rates for stem lepidosaurs and rhynchocephalians according to the unconstrained phylogeny as represented by five most parsimonious trees (MPTs) dated five times using the MBL method.

Discussion

Materials and methods

Taxa and character data

Phylogeny and timescaling

Morphological disparity

Morphological evolutionary rates

Plots of Cramér coefficients

Acknowledgements

Appendix 1

Appendix 1—figure 1. Completeness percentage by time bin for the complete dataset (all lepidosaurs and stem), for time scheme 1.

Appendix 1—figure 2. Completeness percentage by time bin for the complete dataset (all lepidosaurs and stem) for time scheme 2.

Appendix 1—figure 3. Completeness percentage by time bin for squamates and time scheme 1.

Appendix 1—figure 4. Completeness percentage by time bin for squamates and time scheme 2.

Appendix 1—figure 5. Plot of the number of taxa sampled for each time bin in (A) time scheme 1 and (B) time scheme 2.

Appendix 1—figure 6. Percentage of variance explained by each axis.

Appendix 1—figure 7. Plot of Cramér coefficients for PCO1 and PCO2 for the full dataset.

Appendix 2

Appendix 2—figure 1. Morphospace (PCO1 and PCO2) when myological/integument (except for osteoderm characters) are deactivated.

Appendix 2—figure 2. Plot of Cramér coefficients for PCO1 and PCO2 for the dataset without myological/integument (except for osteoderm characters), limited to characters with a coefficient <0.4.

Appendix 3

Appendix 3—figure 1. Chronophylomorphospace for squamates.

Appendix 3—figure 2. Phylomorphospace in 3D for PCO1, PCO2, and PCO3 for squamates only.

Appendix 4

Extended methods

Appendix 5

Extended phylogenetic analyses results