Fossil-informed biogeographic analysis suggests Eurasian regionalization in crown Squamata during the early Jurassic

Ranges & areas

Bird (2003) provided a high-resolution boundary dataset for 54 major, minor, and microplates. These are highly correlated (and causally linked) with terrestrial zoogeographic regions in many respects (Holt et al., 2013; Ficetola, Mazel & Thuiller, 2017), and by extension, previous global area classifications for groups such as amphibians (Bossuyt et al., 2006; Pyron, 2014a; Vilhena & Antonelli, 2015). First, we aggregated these 54 boundaries to the 8 major plates: Africa, Antarctica, Australia, Eurasia, India, North America, Pacific, and South America. We then separated three minor plates that intersect continental areas with endemic or transitional faunas: Arabia, Caribbean, and Sunda (Fig. 1; Appendix S1). We included these as distinct plates for area coding so that their faunas were not artifactually linked to the adjacent continental plates, which had formed distinct areas prior to the emergence of those three landmasses. Alternatively, Caribbean species occurring on recently emerged landmasses such as the Bahamas would be characterized as “North American,” which does not reflect reality in the Mesozoic.

These 11 plate categories describe most of the major qualitative biogeographic regionalizations (see Lomolino, Lomolino & Lomolino, 2010) and of squamate diversity (see Vitt & Caldwell, 2014), while also reflecting an objective paleogeographic reference. Antarctica contains no extant species and very few (primarily marine) fossils and was therefore omitted from further analyses (Legendre et al., 2020). Additionally, the Pacific plate currently contains substantial landmasses (e.g., Baja California, eastern Melanesia, and southern New Zealand) for which the biogeographic history of their squamate fauna is linked closely to the proximal zoogeographic region (e.g., the Nearctic, Oceania). Therefore, we assigned all Pacific species to their nearest major continental plate. We also opted to lump Madagascar in with Africa for this analysis. Usually, Madagascar is counted as its own plate, as it has a unique tectonic history (e.g., Bossuyt et al., 2006). However, the plate has remained relatively close to Africa since the middle Jurassic, so we argue that treating it as a separate unit is less redundant for our analysis given its exceptional level of squamate endemism, unlike groups such as frogs. This has the added benefit of decreasing the complexity of our models. This resulted in nine major global tectonic ecoregions encompassing nearly all extant, described squamate diversity (Fig. 1; Appendix S1).

To determine species occupancy in the nine regions, we first intersected the 10,064 polygon range shapefiles from the Roll et al. (2017) dataset with the plate boundaries for the 11 major plates, removed Antarctica and reassigned the Pacific species. We inspected these assignments for obvious errors, such as recent human-mediated transplant between continents and a few geometry problems. We then matched these species to the 9,755 taxa in the Tonini et al. (2016) phylogeny, for which 9,569 matched natively and 13 were matched to synonyms, for a total of 9,582 species coded using the Roll et al. (2017) maps. We downloaded range-map shapefiles for a further 64 seasnake taxa (Elapidae: Aipysurus, Emydocephalus, Ephalophis, Hydrelaps, Hydrophis, Kerilia, Kolpophis, Laticauda, Parahydrophis, and Thalassophis) from the 5 February 2020 update of the IUCN RedList ( http://www.iucn.org/) and classified them similarly. The 109 remaining species were verified manually using the Reptile Database (http://www.reptile-database.org/) or estimated from the literature. We thereby classified all 9,755 species in the Tonini et al. (2016) phylogeny to the nine plates (Appendix S1). Species were limited in their ranges to a maximum of four plates, as almost no species has a range exceeding four. The exceptions were a few species from the marine sea snake genus Hydrophis. We constrained those that exceed four plates to the ancestral range of Hydrophis (Australia, Eurasia, India, and Sunda) as estimated from the literature (Ukuwela et al., 2017). This was done to reduce the amount of both computational and analytical complexity (see Matzke, 2014).

No biographic coding scenario can account for all possible processes, and some recent Cenozoic patterns may be obscured by our scheme. Examples include the boundary between southeastern Eurasia and Sundaland or between Sundaland and Australasia (Oceania) and the corresponding faunal “lines” discussed by Simpson (1977). Similar criticisms may be leveled at the grouping of most of Central America with Jamaica, Puerto Rico, and the Lesser Antilles on the Caribbean plate, while grouping Cuba and Hispaniola on the North American plate (Fig. 1). The biogeographic history of Caribbean herpetofauna is complex, and this paleogeographical approach accounts for only a portion of recent historical processes (Rosen, 1975; Hedges, 1982; Crother & Guyer, 1996). Regardless, we suggest that this framework is a solid foundation for understanding early squamate biogeography. We anticipate future researchers will refine and revise these classifications using other quantitative methods such as network-based bioregionalization (Vilhena & Antonelli, 2015), incorporating spatial occurrence data and paleogeography to estimate finer-scale processes of recent dispersal and vicariance. Ultimately, we recorded 2,104 species in Africa, 344 in Arabia, 1,529 in Australia, 951 in the Caribbean, 1,202 in Eurasia, 702 in India, 1,173 in North America, 1,953 in South America, and 1,326 in Sundaland. This sums to 1,898 species in landmasses from Laurasian North America and 2,271 in Laurasian Eurasia, and 6,381 in Gondwanan continents and subcontinents.

Note that with nine areas and four allowed in a range, there are 256 possible states—too many to be visualized individually. Only 66 of these were occupied by squamates in our estimates. For visualization purposes, we primarily present our results summarized into four major synthetic, post-hoc ranges based on the ancestral estimates. These are Gondwana (Australia, Africa, Arabia, India, and South America), Laurasia (Eurasia, Caribbean, North America, and Sunda), Pangaea (any Gondwanan + Laurasian range), and Northern Pangaea (Laurasia + Africa). These areas are descriptive and not meant to be exclusive of each other or restricted to landmasses from which their names originated. Consequently, if a clade is reconstructed to have a “Laurasian” origin, that does not necessarily mean the clade originated during the Jurassic while Laurasia was still a united supercontinent; rather, it means that the lineage originated in an area arising from the paleocontinent (i.e., North America, Eurasia, Sunda, and the Caribbean), even if the lineage postdates fragmentation in the Cenozoic. If an ancestral range reconstructs areas belonging to both Gondwana and Laurasia, then it is considered to have a “Pangaean” distribution. The full 9-area results can be seen in Appendix S1. Therefore, while we collapse many Cenozoic biogeographic patterns into “Gondwana” and “Laurasia” for ease of illustration (Fig. 2), complex patterns among the 9 areas are present, particularly since the K–Pg boundary (Appendix S2).

Ancestral range estimation

We used the BioGeoBEARS package (Matzke, 2018) in R 3.6.0 to test several candidate models of biogeographic inference incorporating dispersal, extinction, cladogenesis, and founder-event speciation (Matzke, 2012). We tested the commonly used biogeographic models DEC (Ree & Smith, 2008), DIVA (Ronquist, 1997), and BAYAREA (Landis et al., 2013) with and without the addition of the free ‘j’ parameter, which allows for founder-event-speciation as a possible explanation for range expansion during cladogenesis (Matzke, 2014). Despite criticisms surrounding the ‘j’ parameter (Ree & Sanmartín, 2018), recent results demonstrate the validity of this approach (Matzke, 2021) when modeled and interpreted appropriately, as we are mindful to do here. We did not evaluate other possible parameters such as ‘x’ and ‘n’ (dispersal probability as a function of physical or ecological distance; see (Dam & Matzke, 2016)), or trait-based dispersal models (Klaus & Matzke, 2020), as these are less relevant to the questions here and difficult to optimize for a dataset of this size.

Ancestral range estimates typically produce partial support for multiple possible models. For example, an ancestral range of “Africa” may have 70% support, with 25% for “Africa + Arabia,” and 5% for “Arabia” alone for a given node. Consequently, “Africa” as an ancestral area occurs in 95% of the model space for that node. As we had nine geographic areas included in our analysis, we tested differences in total contribution to the ancestral range estimation from each area. A higher total probability means that the area is included across a greater proportion of estimated potential ancestral ranges. Using this method, we calculated the proportions for all nine areas for 11 major nodes (Squamata, Unidentata, Episquamata, Toxicofera, Dibamidae, Gekkota, Scincomorpha, Lacertoidea, Anguimorpha, Iguania, Serpentes; Fig. 3). From this, we can identify geographic areas of outsized importance in early squamate evolution.

Limitations on biogeographical analyses are not usually related to size of the tree used, but rather the number of areas in the model (Matzke, 2012; Landis et al., 2013). We used the 9,755-taxon tree from Tonini et al. (2016), pruned to the 5,415 species which had molecular data. This includes most described, extant genera, over 50% of all known, extant species, and covers the full range of biogeographic diversity in the group. While some of the calibration for the phylogeny may be updated in the future (see Simões & Pyron, 2021; Title et al., 2024), that concern is beyond the scope of our current study. We also considered the use of the backbone topology of 5,415 species with molecular data superior to the use of the fully sampled posterior distribution with imputed species. The latter would have required us to integrate ancestral area estimates across a sample of trees, a procedure of unclear statistical and biological validity (NJ Matzke, pers. comm.). We compared model fit using the likelihood-ratio test along with AICc to estimate relative likelihoods and select the best-fit model from the set of candidate models (Wagenmakers & Farrell, 2004). Future analyses may incorporate correspondence classes and fossil geographic occurrences to qualify area connectivity due to plate tectonics more finely in ancestral range estimation (e.g., Landis, 2017), particularly for trees including extinct or fossil lineages (e.g., Pyron, 2017; Simoes et al., 2018).

Fossil area constraints

In order to explore the intricate biogeographic history of Squamata, we must first interpret their early fossil record (see Evans, 2003) with respect to recent advances in phylogenetic understanding (Simões & Pyron, 2021). We compiled all verified records of Jurassic and early Cretaceous species from Pangaean landmasses, ranging from the Toarcian to Barremian (Table 1). Recent literature suggests that ancestral range estimation using extant-only datasets may fail to accurately reconstruct ancestral ranges for ancient nodes (Silvestro et al., 2016; Wisniewski, Lloyd & Slater, 2022). To address this issue, we ran models utilizing constraints from geographic occurrences derived from the Jurassic and early Cretaceous fossil record, referred to as “nodal constraints” hereafter. Nodal constraints match fossils of a specific clade to the ancestral lineage for that clade in the phylogeny. Once the fossil taxon has been matched to a node, the node is constrained so that its ancestral range must include the range of which the fossil is located. For instance, if fossil A, located in North America, is used to constrain node 3, then the ancestral range for node 3 must contain “North America” along with any other estimated areas.

This method hinges on the proper identification of fossils and assumes robust phylogenies, which are not always available. For instance, within Squamata, the lizard Jeddaherdan aleadonta was originally thought to be an iguanid from the middle Cretaceous of Africa. However, upon reexamination, morphological evidence suggested that Jeddaherdan was from the late Quaternary and was a member of the genus Uromastyx (Vullo et al., 2023). To help address this uncertainty, we chose fossil taxa based on three consistent criteria: (1) the taxon must be from the Jurassic Period or early Cretaceous, (2) the taxon must be close in temporal proximity to a branch in the squamate phylogeny occurring in the Jurassic or early Cretaceous, and (3) the taxon must be confidently assigned to the specific clade designated by the node based on an explicit phylogenetic or character-based hypothesis. This excludes earlier, stem-squamates (i.e., Megachirella; Simoes et al., 2018), as they can only be assigned to the base of Squamata, negating the purpose of an ancestral range estimation, or to earlier, non-squamate groups (which would be the rhychocephalian Sphenodon) which would be outside of the questions possible using molecular phylogenies of Squamata.

Four taxa met these criteria for use as nodal constraints. The first was the early stem snake Eophis underwoodi, discovered in the Kirtlington Cement Works Quarry, in Oxfordshire, United Kingdom from the Bathonian epoch (Middle Jurassic, 167.7–164.7 Ma) (Caldwell et al., 2015). While Eophis is considered stem-Serpentes (Caldwell et al., 2015), it is temporally closer to the base of Toxicofera. As all stem and crown Serpentes are toxicoferans, we therefore used this fossil occurrence to constrain the Toxicoferan node to include Eurasia. The next was the iguanid Introrsisaurus pollicidens, located in the Guimarota Mine and Alcobaca formation, Portugal from the Kimmeridgian epoch (Late Jurassic, 155.7–150.8 Ma; (Hoffstetter, 1967)). We used this to constrain the Iguania + Anguimorpha node to include Eurasia. The third was Dalinghosaurus longidigitus (Paleoanguimoprha), found in the Yixian Formation in the early Cretaceous of Liaoning, China (Barremian, 130–125.5 Ma; (Evans & Wang, 2005)). We used this to constrain Paleoanguimorpha branch to include Eurasia. The final fossil was Meyasaurus diazromerali, an early member of Laterata located in the Calizas de la Huérguina Formation in Cuenca, Spain (Early Cretaceous, 129.4–126.3 Ma; (Evans & Barbadillo, 1997)). We used this to constrain the Laterata node to the Eurasian plate. A complete list of the fossil history of Squamata in the Jurassic is included in Table 1.

Biogeographic history

Based on the fossil-constrained analysis described above, DEC+J received overwhelming support (AIC_wt ≈1; Table 2). We use this model as our best estimate of squamate biogeography. Ancestral-range estimates (see Appendix S1 for full results) show support for a localized Pangaean origin (Africa, Australia, Eurasia, and Sunda) of extant Lepidosauria (Squamata + Rhynchocephalia), which is congruent with Gondwanan Triassic presence of rhynchocephalians (e.g., Bonaparte & Sues, 2006) and Laurasian Triassic presence of squamates (Simoes et al., 2018). Subsequently, the node subtending extant squamates is estimated to have the same ancestral range of Pangaea (Africa, Australia, Eurasia, and Sunda) ∼180 Ma, some 68 million years after divergence from the rhynchocephalians, approximately coinciding with the beginning split between Gondwana and Laurasia in the Early Jurassic (Fig. 3).

The earliest diverging squamate lineage (Dibamia) exhibits a subsequent contraction from Pangaea into Laurasia (North America, Eurasia, and Sundaland) during the long ghost lineage between their ∼180 Ma divergence from the remainder of Squamata and their most recent common ancestor (MRCA) ∼73 Ma, possibly as the final portions of the Sundaland plate had collided with mainland Eurasia (Metcalfe, 1998). We estimate an ancestral range for Gekkota identical to that of Lepidosauria and Squamata (Africa, Australia, Eurasia, and Sunda). Most of the later-diverging lizard lineages in the Jurassic show an entirely Eurasian distribution, including Episquamata, Unidentata, Toxicofera, Laterata, Anguimorpha, and Iguania (see Burbrink et al., 2020 for clade definitions).

In the early Cretaceous, most stem lineages were regionalized to Laurasia, though we estimate a broad ancestral range for Scincomorpha, which shows a northern Pangaean distribution (Africa and Eurasia). Crown Serpentes also evolved in the Cretaceous, and the estimated ancestral range for extant snakes is also Pangaean (Africa, Eurasia, South America, and Sunda). Within Serpentes, we estimate an African origin for Scolecophidia and a South American origin for Alethinophidia, showing the first examples of Gondwanan regionalization in Squamata. In contrast with widespread Laurasian regionalization in the Jurassic and early Cretaceous, Gondwanan endemic lineages of extant taxa are established by the end Cretaceous in Gekkota, Iguania, Laterata, Scincomorpha, and Serpentes. These arose through several different mechanisms: via contraction from a broad Pangaean range in Scolecophidia, Alethinophidia, and Iguania; contraction from a Northern Pangaean range in Scincomorpha, some Gekkota, and colubroid snakes; or dispersal from Laurasia in some Gekkota and Laterata (Fig. 3).

In the Cenozoic, we estimate fewer transitions between ancestral areas designated by Laurasia, Gondwana, and northern Pangaea, and observe an increasing frequency of dispersal events within Old World and New World landmasses, with relatively infrequent instances of dispersal between the Old and New Worlds (Fig. 4B). Corresponding to its estimated status as a major locality of origination for Squamata (Table 3), Eurasia continued to act as a source of diversity for squamates throughout the Cenozoic, with numerous dispersals into Africa, India, and Sunda. Similarly, Sunda also acts as a major source area, with frequent dispersals into Eurasia, along with several events into Australia and India. Interestingly, both Eurasia and Sunda also show the highest rates of extirpation (i.e., lineage contraction from an area), with Eurasia having nine lineages contracting from the continental plate and Sunda having two. We also estimated one extinction each in Africa, North America, and India (Fig. 4B).

Ancestral range estimates for Mesozoic squamate lineages received insignificant contributions from the Arabian, Caribbean, Indian, and South American plates (Table 3). Some of these such as the Caribbean were not land positive in their current form in the Mesozoic, suggesting that these plates were not important areas of origination and diversification of squamate lineages. In contrast, it seems likely that South America and India had diverse endemic squamate faunas that are not represented among extant lineages but may be present undiscovered in the fossil record. For the root node of Squamata, the model shows significant contributions from the African, Eurasian, and Sunda plates as part of the ancestral range. These three continental plates are clearly important for the early evolution of Squamata across the Mesozoic (see Evans, 2003; Simões & Pyron, 2021) and should be an active area of research in the study of ancient squamates, particularly the phylogenetic placement of fossil taxa from those regions.

Biogeographic models

We estimate overwhelming support for a DEC+J model (Table 2), reflecting strong evidence for long-distance dispersal and colonization (particularly during the Cenozoic), a widely implicated mode of biogeographic range expansion for squamate lineages moving among former Pangaean landmasses currently or previously separated by great distances (Longrich et al., 2015). Recent authors (Ree & Sanmartín, 2018) raised potential issues with biogeographic inference using models (specifically DEC-type models) with founder-event speciation. However, their primary critique lay with cladogenesis parameters, which they claimed are problematic in being modeled probabilistically without respect to time in DEC-type models, an issue they further claim is compounded when additional weight is assigned these events via the ‘j’ parameter. We note that our best-fit model (DEC+J) does not seemingly over-estimate the jump parameter in respect to the dispersal, as the dispersal parameter (d = 0.0015) was estimated to be an order of magnitude higher than the parameter for the jump dispersal (j = 0.008; Table 2). Matzke (2021) showed that these concerns are likely unfounded for most empirical analyses.

As noted, we employed a geographic strategy of area coding based on present-day plate-tectonic boundaries, which are broadly congruent with modern zoogeographic regions (e.g., Holt et al., 2013). This approach has distinct advantages and limitations. Chief among the benefits is its discreteness, objectivity, and unambiguousness, and it can be replicated easily in other taxa for comparative analyses. This eliminates subjective uncertainty regarding where to draw lines between terrestrial ecoregions like the Nearctic and Tropical Middle America (e.g., Kreft & Jetz, 2010; Pyron, 2014a; Pyron, 2014b; Vilhena & Antonelli, 2015). Disadvantages include fixing boundaries between terrestrial ecoregions to their paleogeographic origins, which may not correspond with empirical evidence of the recent ecologically mediated biogeographic processes affecting related species in those regions (Mucina, 2018).

We attempt to avoid implying excessive certainty in our model-based reconstructions. While support for the most-likely ancestral range is high at nearly all internal nodes (Fig. S3; Table 2), we nevertheless concede that there are several factors that may confound such inferences. The first is the existence and phylogenetic placement of extralimital fossil taxa, species that are outside the modern range of the clade, which therefore does not reflect its historical distributions (Head, 2021). The second is topological uncertainty, such as the placement of Dibamia or relationships within Iguania (Burbrink et al., 2020). Third is variation in divergence-time estimates for nodes across the tree (see Jones et al., 2013; Simões & Pyron, 2021). While there are disparate estimations between phylogenomic (Irisarri et al., 2017) and total-evidence data (Pyron, 2017) especially for groups such as Gekkota, it is unlikely to affect the ancestral range estimation. Where differential divergence times can potentially impact results is in dispersal versus vicariance estimations, which is not something we attempted to capture in our study. Fourth are conceptual deficiencies in all existing biogeographic models, which include inadequate parameterization of extinction (see Matzke, 2014), and incomplete integration of correlated diversification and range evolution (see Goldberg, Lancaster & Ree, 2011; Goldberg & Igić, 2012; Caetano, O’Meara & Beaulieu, 2018). Recent results suggest that fully identifiable models for the latter may not be possible (see Louca & Pennell, 2020).

Regardless, we are confident that our results represent a substantial advance in presenting at least a partially informed view of the early biogeographic history of squamates, one that proposes several distinct hypotheses that can be tested further. Future studies could incorporate other parameters, such as trait-based dispersal (Klaus & Matzke, 2020), dispersal as a function of ecological or geographic distance between areas (Dam & Matzke, 2016). Body size (e.g., Feldman et al., 2016), limb reduction (Vidal et al., 2008), and the temperate/tropical transition (Pyron, 2014b) all potentially impact dispersal and are avenues worth pursuing.

Finally, a recent methodology aimed at incorporating both fossil geographic occurrences and complex geographic interaction across time has been used for ancestral range estimations (Landis, 2017), and could be implemented in future studies. To do this effectively, more Jurassic squamate fossils need to be included in morphological matrices for a better understanding of their phylogenetic relationships (see Conrad, 2008; Wiens et al., 2010; Gauthier et al., 2012; Reeder et al., 2015; Pyron, 2017; Simoes et al., 2018). Paramacellodid lizards in particular will likely have a strong bearing on these estimates, given their Jurassic and early Cretaceous distribution across Pangaea and uncertain placement within Squamata (see Bittencourt et al., 2020). There is a well-known bias in the fossil record towards species in Laurasia (Vilhena & Smith, 2013; Valenzuela-Toro & Pyenson, 2019; Croft & Lorente, 2021), which may have impacted our results here, as all squamates available were from Laurasia in the Jurassic and early Cretaceous. In particular, the results presented here suggesting a Eurasian regionalization for Squamata would be falsified by fossil occurrence of relevant lineages in former Gondwanan continents such as South America, Australia, India, or Antarctica. Prospecting and research on these locations could greatly benefit future studies.

Origin and diversification of extant Squamata

The phylogenetic relationships and present-day distributions of living squamates, along with well-constrained fossil taxa from the Jurassic and early Cretaceous, contain signal suggesting that the earliest squamate lineages (represented by the basal branches of the phylogeny of extant species) were localized along the plates that formed the coastline of the paleo-ocean, the Tethys Sea (Zhu, Zhao & Zhao, 2022; Fig. 3). Continuing into the Jurassic, a strongly supported pattern emerges of Eurasian origin for many major groups of squamates as Gondwana and Laurasia continued to break apart (Fig. 3). This geographical regionalization persisted at least until the mid to late Cretaceous, when Gondwanan radiations were established in Serpentes for Alethinophidia and Scolecophidia, the two major lineages of snakes (Fig. 3).

Later in the Cretaceous, Teiidae and later-diverging Gekkota also are reconstructed to have Gondwanan ancestral ranges (Fig. 3, for exact ranges see Appendix S2). These localizations formed via contraction from Pangaean or Northern Pangaean in Gekkota and Serpentes in contrast to Teiidae, which resulted from a jump dispersal from Laurasia to Gondwana (Fig. 3). This is coincident with and likely related to a substantial period of tectonic vicariance as the paleocontinents continued to diverge (Bird, 2003). It is worth noting that four of the nodes described above were constrained utilizing fossils and may seem circular. However, fossils are the only concrete evidence of the presence of a group in an area. Therefore, we argue that the estimated ancestral ranges are empirically valid and supported, since we are employing the increased certainty that comes along with fossil geographic occurrences as empirical data.

At more recent timescales, our results are congruent with most studies using molecular time trees to infer biogeographic histories to the extent that the area codings are comparable. For more detail, we refer to previous works on Dibamia (Townsend, Leavitt & Reeder, 2011), Gekkota (Agarwal et al., 2014; Carranza et al., 2000; Carranza & Arnold, 2006; Gamble et al., 2011; Gamble et al., 2008; Šmíd et al., 2013), Scincidae (Honda et al., 2000; Carranza & Arnold, 2003; Whiting et al., 2006; Brandley et al., 2011; Barley et al., 2015; Pereira & Schrago, 2017), Teiidae (Giugliano, Collevatti & Colli, 2007), Amphisbaenia (Longrich et al., 2015; Vidal et al., 2008), Lacertidae (Arnold, Arribas & Carranza, 2007; Tamar et al., 2016), Anguimorpha (Macey et al., 1999; Vidal et al., 2012), Iguania (Macey et al., 2000; Noonan & Sites Jr, 2010; Okajima & Kumazawa, 2010; Grismer et al., 2016), and in snakes, Typhlopoidea (Adalsteinsson et al., 2009; Vidal et al., 2010), Booidea (Graham Reynolds, Niemiller & Revell, 2014; Noonan & Chippindale, 2006a), and Colubroidea (Burbrink & Lawson, 2007; Alfaro et al., 2008; Wüster et al., 2008; Guo et al., 2012; Chen et al., 2013; Chen et al., 2017).

Subsequent studies may take a more focused perspective on individual lineages to enhance our understanding of the complex geographic history of Squamata. For example, colubroid snakes and Old World skinks have a vast, complex, and rather recent biogeographic history that deserves more attention than given here or in recent studies (see Cadle, 1985; Greer, 1970 for some early hypotheses and discussion). Similarly, little recent attention has been paid to teiid, lacertid, or anguimorph lizard biogeography at the continental or global scale since updated molecular phylogenies have become available (i.e., Tonini et al., 2016).

A final question of great paleoecological and biogeographic interest is the directionality or asymmetry of dispersal between the landmasses analyzed here. Recent studies have shown that a variety of taxa experienced higher rates of dispersal into North America from Asia across Beringia than the opposite (the Cenozoic Beringian Dispersal Hypothesis; Guo et al., 2012; Jiang et al., 2019). Similarly interesting results might be obtained by examining rates of interchange between North and South America (see Estes & Báez, 1985 for early speculation), between Africa and western Eurasia (e.g., Georgalis, Villa & Delfino, 2016; Tamar et al., 2016), between India and Eurasia (e.g., Agarwal et al., 2014; Datta-Roy et al., 2012; Grismer et al., 2016), and between Australasia and eastern Eurasia (see Oliver & Sanders, 2009) during the Cenozoic. A great deal remains to be learned about the biogeographic history of Squamata, and our results will provide a robust foundation for productive future investigations.

Preliminary assessment of our Cenozoic results (Fig. 4A) reveals relatively few lineages (only 463 out of 5,415 species, <10%) with a recent Pangaean (i.e., cosmopolitan) distribution. Regarding Cenozoic changes in distribution, there are very few range transitions from Laurasia to Gondwana (only 143 instances) and Gondwana to Laurasia (only 124 instances), suggesting a low rate of recent jump dispersal between the two former supercontinents. Accordingly, the network analysis (Fig. 4B) reveals a much clearer pattern of dispersal in the Cenozoic, which clearly shows a separation between clusters of dispersal within the New World (South America, North America, and the Caribbean) and Old World (Arabia, Australia, Africa, Eurasia, Sunda, and India), and very few events between the Old and New Worlds. This contrasts with the Gondwana-Laurasian division seen in our data during the Mesozoic (Fig. 3).

Additional fossil sampling, particularly in the Southern Hemisphere, will likely increase our knowledge of early Gondwanan squamate evolution. These data may broaden the ancestral range estimate to ‘Pangaea’ (see San Mauro et al., 2005 for a similar situation in Amphibia), completely alter it to ‘Laurasia’ or a component thereof (see (Zhang et al., 2013) for a similar result in plants), or a Gondwanan estimate with additional earlier dispersals into Laurasian landmasses (see (Gardner, Surya & Organ, 2019) for a similar example in early tetrapodomorphs). Recent Antarctic discoveries of an Eocene frog (Mörs, Reguero & Vasilyan, 2020) and late Cretaceous mosasaur material (Martin, 2006; Legendre et al., 2020) suggests the potential for a rich fossil history on this continent, which is almost entirely unknown in studies of most organisms (see Noonan & Chippindale, 2006a). Although Antarctica poses no significance for our study (as no modern squamates occupy the continent), other terrestrial tetrapod groups and marine squamates have been found in Antarctica (Goin & Goin, 1972; Rozadilla et al., 2016; Estrella et al., 2019), suggesting that with future exploration, potential terrestrial Antarctic squamates would expand our understanding of the biogeographic patterns on the continent.