Successive climate crises in the deep past drove the early evolution and radiation of reptiles.

Climate change-induced mass extinctions provide unique opportunities to explore the impacts of global environmental disturbances on organismal evolution. However, their influence on terrestrial ecosystems remains poorly understood. Here, we provide a new time tree for the early evolution of reptiles and their closest relatives to reconstruct how the Permian-Triassic climatic crises shaped their long-term evolutionary trajectory. By combining rates of phenotypic evolution, mode of selection, body size, and global temperature data, we reveal an intimate association between reptile evolutionary dynamics and climate change in the deep past. We show that the origin and phenotypic radiation of reptiles was not solely driven by ecological opportunity following the end-Permian extinction as previously thought but also the result of multiple adaptive responses to climatic shifts spanning 57 million years.

INTRODUCTION

A major challenge in macroevolution and macroecology is understanding the magnitude and duration of organismal response to rapid environmental changes. On ecological time scales, it is well known that organismal diversity (), population density (), morphology (, ), and life history traits—such as clutch size, age of sexual maturity, and longevity ()—can all be affected by fast shifts in the environment, such as those caused by anthropogenic climate change. However, inferring long-term organismal response to environmental perturbations on geological (or evolutionary) time scales is considerably more challenging. Fortunately, the fossil record can provide direct evidence of past environmental crises and how organisms were affected and recovered from such events. The current rate of global warming and species extinction only has parallels deep into the geological record (, ), thus making our understanding of global climatic shifts in the geological past especially informative for inferring its long-term consequences for organismal evolution.

One of the best measures of long-term biological response to a changing environment is the reconstruction of organismal evolutionary rate trajectories, either genetic or phenotypic. These trajectories can be used to understand the adaptive responses of organisms to a changing climate (–), as evolutionary theory predicts that phenotypic evolutionary rates are tightly connected to biotic (e.g., competition for resources) and abiotic (e.g., temperature) environmental variables at both local and global scales through deep geological time (, ). However, such an approach has rarely been applied to empirical study systems because of the limited availability of data or quantitative analytical tools to explore such questions at broad macroevolutionary scales. Hence, reconstructing the rate of phenotypic evolution across major environmental shifts in deep time is a challenging empirical task, although it provides unique insights into Earth-life system interactions over vast tracts of time.

Two of the biggest climatic crises in Earth’s history occurred at the end of the Permian. Those events were triggered by unprecedented rates of global warming culminating in two major mass extinctions—the end-Guadalupian [EGE; at ~261 million years (Ma ago)] and the Permian-Triassic mass extinctions (PTE; at 251.9 Ma ago) (, ). The PTE represents the largest mass extinction in the history of complex life, with an estimated loss of ~86% of species worldwide (). In contrast, the EGE has only been recognized as a major mass extinction in recent decades, but evidence suggests that it was similar in magnitude to other major mass extinctions in Earth’s history, including the Cretaceous-Paleogene mass extinction (). Both the EGE and PTE events reset entire ecosystems worldwide (e.g., extinction of several groups of early synapsids on land and major reef-building organisms in the oceans, such as brachiopods) and opened up ecological opportunities for surviving animal groups to undergo adaptive radiations during the subsequent geological period—the Triassic (, , , ).

Reptiles represent one such group that survived the Permian-Triassic environmental crises and adaptively radiated during the Triassic, occupying habitats and ecological roles previously dominated by other organisms, such as early synapsids—the evolutionary forerunners of mammals (, , , ). These factors make reptiles a rare study system to understand the impact and successful recovery of tetrapod faunas in the aftermath of cataclysmic environmental change. However, several aspects of how reptiles responded to these climatic crises remain unanswered (, ), owing to logistical challenges in collecting vast amounts of phenotypic data for fossil vertebrates and in reconstructing organismal evolutionary trajectories in the deep past. Furthermore, such questions are made even more complex by ongoing conflicts regarding fundamental aspects of the reptile tree of life. For instance, there is a strong disagreement between morphological and molecular phylogenies on the relationships of important lineages (e.g., turtles), which remain hotly debated (–). These conflicts hamper the reconstruction of a precise timeline of reptile evolution, and yet, they remain critical for accurately calibrating major events in their evolutionary history using either morphological or molecular data.

Here, we address these challenges by constructing a new morphological dataset inclusive of all major lineages of early reptiles, early synapsids, and their closest relatives (early amniotes). This new dataset represents the largest taxon sample ever assembled to assess this portion of the tree of life, being 40% larger than the previously largest dataset () and the first with the vast majority of all sampled species data based on primary data collection from direct specimen observation and micro–computed tomography scans. We explore this dataset with emerging techniques from Bayesian evolutionary analyses that combine species relationships, diversification dynamics, and evolutionary rates within an integrated and multivariate statistical framework (, –). We further combined our results with data from body sizes and global temperature data for over 140 Ma, from the Late Carboniferous to the Late Triassic, to provide a broad overview of the major biological responses toward climatic shifts at this critical interval in Earth’s history. This approach enabled us to shed new light on current conflicts in early amniote and reptile relationships and allowed us to establish a precise timeline for the reptile tree of life and the timing of the reptile-mammal split and to answer the following questions: (i) How do rates of phenotypic evolution in reptiles respond to major changes in temperature patterns during past climate change events? (ii) How can the detected long-term patterns inform on fundamental ecological and evolutionary processes, such as selective pressures and organismal adaptability, operating before, during, and after global climatic crises?

Our results provide the first instance of agreement between morphology- and molecular-based hypotheses of early amniote and reptile evolution and a substantially more precise estimate for the time of origin of their major evolutionary branches compared to previous studies. We shed light on the primary causes for the historical challenge of placing turtles into the reptile tree of life using either morphological or molecular datasets. Furthermore, we find a heretofore undetected and intimate relationship between long-term global mean temperature trends and rates of phenotypic (i.e., morphological) evolution in reptiles between the Permian and Triassic. Last, we discuss how reptiles responded adaptively to the end-Permian extinction events, which were marked by some of the most drastic environmental crises in Earth’s history, finding that ecological opportunity was not the sole driver for the adaptive radiation of reptiles following the PTE.

RESULTS AND DISCUSSION

A revised early amniote and reptile time tree of life

Phylogenetic analyses using distinct optimality criteria and models (see Materials and Methods) indicate a strong support for crown amniotes—a monophyletic Synapsida and Sauropsida (Reptilia or reptiles here)—but with the exclusion of captorhinids, protorothyridids, and araeoscelidians from reptiles. The latter groups are here inferred as stem amniotes, forming a clade that is sister to crown Amniota (Fig. 1 and figs. S1 to S10). This rather unorthodox relationship goes against most recent hypotheses of early amniote and reptile evolution (, –), but it is supported by numerous character-state transitions (Supplementary Text) and also recently inferred by another study that used a distinct dataset ().

Fig. 1.

Revised evolutionary tree and divergence times for the major groups of early amniotes, early synapsids, and reptiles.

Maximum compatible tree from relaxed morphological clock Bayesian inference analysis under the SFBD tree model with a single clock partitioning and best-fitting models (see Materials and Methods and table S11). Node values represent median ages; node violin plots (cyan) represent the distribution of the 95% highest posterior density intervals. Evolutionary branches leading to taxa originating after the Jurassic are omitted for simplicity (for full tree, see figs. S5 and S6).

The removal of araeoscelidians—long considered to be the hallmark of early diapsid evolution—from reptiles to the amniote stem has major impacts on reptile systematics and the recognition of Diapsida as a natural group (i.e., clade). Diapsids have historically been circumscribed as the group including araeoscelidians, younginiforms, and all their descendants to the exclusion of captorhinids and parareptiles [see Supplementary Materials and (, )]. Therefore, our results [and ()], where crown reptiles (Sauria) and younginiforms share a more recent common ancestor with “parareptiles” and synapsids than araeoscelidians, necessarily imply on the polyphyly of Diapsida. As more inclusive reptile groups, such as Sauropsida and Eureptilia, have been previously defined using Diapsida as a reference point—e.g., (, )—this necessitates a substantial revision of early reptile taxonomy, which will be addressed in greater detail elsewhere.

Furthermore, our results support the traditional hypothesis of varanopids as a subgroup of synapsids, not as early sauropsids as recently suggested (). We also find strong evidence for the paraphyly of groups traditionally classified as parareptiles as early-deriving sauropsids. Instead of a single group of early reptiles, parareptiles comprise at least four distinct radiations in the reptilian stem: Acleistorhinidae, Millerettidae, Procolophonia (Procolophonoidea + Pareiasauria), and Mesosauridae (Fig. 1). Of notable interest is the phylogenetic placement of acleistorhinids—a poorly known group of early Permian reptiles from North and South America—as the earliest reptilian clade to evolve. This result stems from numerous morphological features that acleistorhinids share with early synapsids (e.g., caseasaurs), including a large overlap between the lacrimal and the nasal process of the maxilla (the latter lying lateral to the lacrimal), a narial fossa on the premaxillary process of the maxilla, the presence of anterior caniniform teeth, and teeth on both the anterior coronoid and posterior coronoids (see Supplementary Text).

Last, our expansive taxonomic sampling of both early- and later-deriving reptiles provided the opportunity for a new look into one of the oldest problems in reptile evolution—the phylogenetic origin of turtles. Accurately inferring the placement of turtles among other groups of reptiles suffers from several hurdles, affecting both morphological and molecular data. These include limited taxonomic sampling in most previous morphological datasets [focused either on early reptiles (, –) or later-deriving reptiles (, , ), but not both simultaneously—see () for an exception], a large number of entirely extinct early reptile lineages for which no genomic data can ever be made available—nearly 80% of all groups sampled here (data file S1), and the extremely poor fossil record of stem turtles, with only a handful of known taxa from the Triassic (and none from the Permian if Eunotosaurus is excluded) (, ). The result has been a long history of conflict among existing morphological hypotheses and between those and molecular hypotheses. Previous morphology-based studies have inferred turtles to be a subgroup of parareptiles (closely related to the extinct Procolophonia) (, , ), lepidosauromorphs (the group inclusive of lizards, snakes, and tuataras) (, , ), or early-deriving diapsids (, ). On the other hand, molecular studies, from smaller sequence data to phylogenomics, strongly support turtles as the sister group to archosaurs among extant taxa (, , –), thus making turtles more closely related to crocodiles and birds than to any other reptile. However, even among molecular studies, several genes have been recognized as strongly supporting alternative placements of turtles, especially within lepidosauromorphs, indicating important areas of conflict within the molecular signal (, ).

Contrary to all previous morphological phylogenetic studies, we find turtles to be more closely related to archosaurs than other living reptiles—forming the clade Archelosauria (Fig. 1) ()—in agreement with the overwhelming majority of molecular studies over the past two decades (, , –). This agreement of morphological data with the dominant molecular hypothesis of turtle origin is considered here to be a result of our much broader taxonomic sampling, including all reptile groups to which turtles have been previously proposed to be more closely related, and new data provided by key fossil species (). The latter combined with emerging tip-dated Bayesian inference models for morphological data (see Materials and Methods) also provided considerable precision (narrow 95% highest posterior densities) for divergence time estimates for most nodes (Fig. 1 and fig. S8), even nodes closer to the root that are notoriously harder to estimate ().

It is important to note, however, that the node supporting Archelosauria has a low posterior probability, suggesting a pervasively strong conflict within the morphological signal. Only two unambiguous synapomorphies (estimated with maximum parsimony ancestral state reconstruction) are found to be supporting turtles with archosaurs: the presence of a sagittal crest on the supraoccipital and the absence of a humeral entepicondylar foramen—although the latter is conspicuous to archosauromorphs, has a low level of homoplasy (retention index = 0.85), and is occasionally used to identify isolated archosauromorph remains [e.g., ()]. A similar issue has been detected for molecular datasets and has been suggested to result from short time intervals separating the divergence of turtles from other major groups of reptiles, especially archosaurs (). This would result in less time available for enough changes to accumulate in some genes, with trees derived from these genes therefore containing less information to infer evolutionary relationships (), besides further complications arising from gene paralogy (). The former issue also affects the morphological data, as short time tree internodes are expected to reduce the likelihood of morphological changes to accumulate and to be detected on preserved fossils. Our divergence time estimates indicate an exceptionally short time span separating the base of Archelosauria and the next diverging node (only 1.9 Ma between median ages). This is one of the shortest time intervals between any internal nodes separating major clades in our results (Fig. 1 and fig. S8). We also find exceptionally short internal branches separating the divergence of the major groups of archosauromorphs (e.g., only 2.8 Ma between Archosauromorpha and the archosauriform + protorosaur node) (Fig. 1). Together, our substantially more precise timeline of reptile evolution provides the first morphological support for Archelosauria and empirical evidence that supports previous predictions of a short time interval between the origin of turtles and archosauromorphs, with the latter as the primary cause for internal conflicts in both molecular and morphological data to infer the phylogenetic origin of turtles.

Evolutionary homogeneity rather than mosaicism is predominant among early amniotes

Here, we used both single and partitioned relaxed clock Bayesian evolutionary analyses for examining the rates of morphological evolution in early amniote and reptile evolution. This fast-evolving technique has provided unique insights into deep time macroevolutionary dynamics across various portions of the tree of life in recent years—e.g., (, –). We find that early amniotes, including early synapsids and most early reptile lineages, evolved under relatively low rates of morphological change (≤1.0) (Fig. 2A, figs. S6 and S9 to S11, and table S1). The major exception is Archelosauria, with its major clades (e.g., archosauriforms, protorosaurs, marine reptiles, and turtles) undergoing accelerated evolutionary rates of up to ~2.5 times higher than clock base levels. This pattern is the same for rates measured from either cranial or postcranial subdivisions of their phenotype (Figs. 2A and 3, figs. S9 and S10, and tables S1 and S2).

Fig. 2.

Relative rates of evolution across subdivisions of the early amniote, early synapsid, and reptile phenotype.

(A) Mirrored trees depicting median relative rates of morphological evolution under the SFBD tree model using partitioned relaxed morphological clocks (see also figs. S9 and S10). (B) Kernel density plots depicting rate distributions for cranial and postcranial partitions for each major synapsid and reptile clade. Names in bold* indicate clades with significantly different rates between morphological partitions (i.e., mosaic evolution). Relative rates increase substantially in Archelosauria relative to all other groups of reptiles and early synapsids, but rate heterogeneity occurs in both slow- and fast-evolving evolutionary lineages.

Fig. 3.

Relative evolutionary rates by morphological region and across time for reptile groups under strong selective regimes between the middle Permian and Late Triassic.

(A) Relative cranial evolutionary rates. (B) Relative postcranial evolutionary rates. Horizontal black dashed line indicates the base of the clock rate. Protorosauria* includes a mixture of terrestrial (ground dwelling and arboreal) and semiaquatic taxa. Vertical red dashed lines indicate the time of major mass extinctions caused by global warming: the EGE and the PTE. Each data point represents a sampled species in the morphological dataset (depicted in Fig. 1 and figs. S1 to S10).

When comparing the rates between cranial and postcranial partitions within each major clade, most values are not significantly different (homogeneous) between morphological partitions (Fig. 2 and table S3). Significantly different (heterogenic) rates between morphological partitions are detected for only 7 of the 22 (34%) major groups assessed (Fig. 2). This result is unexpected, considering the quite different functional and adaptive demands acting upon cranial and postcranial morphological characters—e.g., sensory and trophic adaptations in the skull versus locomotory adaptations in the postcranium.

Such distinct functional demands among taxonomic and morphologically disparate animal groups should result in a stronger contrast between cranial and postcranial rates. Accordingly, among the seven clades evolving at heterogenic rates, three were evolving semiaquatic to fully aquatic morphotypes (younginiforms, thalattosaurs, and sauropterygians), and one includes a wide mixture of aquatic-, arboreal-, and ground-dwelling forms (protorosaurs), resulting in significantly faster rates of morphological change in their postcranium relative to the skull (Figs. 2B and 3). However, the overall pattern detected is that such cases of mosaic evolution—characterized by phenotypic traits evolving at substantially different rates—do not characterize the bulk of early amniote phenotypic evolution. This pattern suggests that strong selective pressures were acting with similar magnitude upon both cranial and postcranial characters, regardless of functionality, during this phase of early amniote evolution (see more below).

Rates of morphological evolution during the Permian-Triassic environmental crises

Comparing long-term geological trends in global mean temperature () to rates of morphological change in reptiles, from the Early Carboniferous to the Late Triassic, at step intervals of 1 Ma (Fig. 4, A to C), reveals a number of major evolutionary trends. During the Carboniferous, there was a steady overall increase in global mean temperature from 25°C until it peaked to around 30°C during the Late Carboniferous (ca. 310 Ma ago) and then slightly decreased back to 25°C toward the end of the Carboniferous (Fig. 4B). However, these long-term changes were not matched by major rates of temperature change (Fig. 4C), which stayed below 3°C of change at every 1 Ma. During this time, stem amniotes and early synapsids were undergoing a steady decrease in their evolutionary rates, while the earliest reptiles held stable at just below base rates (Fig. 4A and figs. S12 and S13).

Fig. 4.

Rates of morphological and global temperature change during early reptile evolution.

(A) Relative rates of morphological evolution (for whole body) across time for the major groups of reptiles. See fig. S6 for rates in a time tree and figs. S12 and S13 for results inclusive of early synapsids. SAC, Sakmarian-Artinskian Crisis. (B) Mean sea surface temperatures (SSTs) at every 1 Ma (black dots) and their respective 95% confidence interval (CI; blue vertical lines) [from data file S1 and ()]. Loess best-fit line across all mean SST values (black solid line) and its 95% CI (gray area) indicate geological long-term temperature trends. (C) Rates of temperature change (the total step change in global mean SST per million year) (from data file S1). Horizontal black dashed lines indicate critical temperature threshold of 3° and 5.2°C. Vertical red dashed lines indicate the time of major mass extinctions caused by global warming: the EGE and the PTE. The vertical orange dashed line marks a mass extinction caused by global cooling: the SAC.

The situation changes markedly during the Permian, which is governed by substantial rates of temperature change, with several points peaking above 3°C/Ma and two above the critical 5.2°C/Ma rate level—a recently recognized critical threshold common to all periods of climate change directly related to major mass extinctions (Fig. 4C and table S4) (). These climate change events were characterized by both large-scale positive and negative temperature fluctuations—exceptional periods of global warming and global cooling, respectively—with several data points and their respective 95% confidence interval (CI) either above or below the CI of the regression line (Fig. 4B). Such temperature alternation generated a low net temperature change during most of the Permian (illustrated by the mostly horizontal regression line), which might create a false impression of temperature stability. The big climate changes of the Permian start with the Sakmarian-Artinskian Crisis (SAC)—caused by a major glaciation event (Fig. 3C) ()—and mark the onset of a long-term pattern of phenotypic radiation in reptiles (Fig. 4A). During the latter part of the Permian, we observe the origin of several new groups of reptiles [Fig. 1 and (, ))] and increasingly higher evolutionary rates (Fig. 4A and figs. S6, S9, and S10), which continue to increase toward the end of the Permian concomitant with substantial increases in the global mean temperature during the EGE and PTE crises.

Subsequently, during at least the first 10 Ma of the Triassic, global mean temperatures were substantially higher compared to the Permian and the Carboniferous (Fig. 4B) and are further characterized by high rates of temperature change with links to global warming—one during the earliest Triassic and three in the Middle Triassic (Fig. 4C). Great fluctuations in CO2 levels and isotopic carbon over the course of 4 Ma during the Early Triassic indicate a long-sustained greenhouse crisis (, ). It is during this period that we see the highest rates of morphological change in early reptile evolution, characterized by the diversification of the archosauriform and sauropterygian body plans and, at a smaller rate scale, ichthyosauromorphs, protorosaurs, and rhynchosaurs (Figs. 2 and 3A). Despite some significant differences in rates of evolution between cranial and postcranial morphological subdivisions among those lineages, both rates are always above base rate levels (≥1; Fig. 3 and figs. S14 and S15). A major exception to this pattern is that depicted by lepidosauromorphs, which show steadily decreasing rates of morphological change from the latest Permian until the end of the Triassic (Figs. 3 and 4A), which is discussed further below.

It is only during the transition to the Late Triassic that temperature levels substantially decrease, although they are still higher compared to the Permian and Carboniferous (Fig. 4B). Furthermore, rates of temperature change decrease during the Late Triassic—and are inferior to those in the Early and Middle Triassic—until the late Norian (215 Ma ago). Similarly, during the Late Triassic, there is a steady decline of evolutionary rates for all reptiles, including archelosaurs, closer to the basal rate level (~1.0; Fig. 4A and figs. S12 to S15).

Stabilizing selection among early reptiles broken by archelosaurs

It has long been argued that stasis is the predominant phenotypic macroevolutionary mode in the fossil record, being interspersed with geologically short time intervals of fast phenotypic change (–)—a pattern that has formed the basis of the theory of punctuated equilibria (). Such long-term patterns of phenotypic stasis are suggested to be the product of stabilizing selection (, ). One factor contributing to stasis is a species’ ability to migrate to new habitats as an effort to secure optimal environmental and climatic conditions for their survival. Therefore, periods of time where lineages break away from stabilizing selection and stasis and undergo positive (directional) selection to new adaptive zones are expected to occur more rarely, such as during geologically fast changes in biotic (e.g., introduction of predators) or abiotic (e.g., temperature) conditions ().

On the basis of this theoretical foundation, rates of morphological evolution can provide key insights into the mode and strength of natural selection operating on distinct phenotypic regions (, –) and are the only source of information available to understand such processes in mostly or entirely extinct evolutionary lineages (, , ). Using our recent approach to infer selection mode from Bayesian clock rates (, ), we detect, according to theoretical expectation, widespread stabilizing selection in most groups of early amniotes, early synapsids, and reptiles (Fig. 5). This pattern is only broken by fast-evolving archelosaurs, indicating strong positive selection toward new morphotypes in both their cranial and postcranial partitions (Fig. 5).

Fig. 5.

Inferred mode of selection across subdivisions of the early amniote, early synapsid, and reptile phenotype.

Branches in shades of blue indicate rates significantly lower than the background rates, indicative of stabilizing selection (or stasis). Branches in shades of red indicate rates significantly higher than the background rates, indicative of positive (or directional) selection. Branches in gray indicate evolutionary rates not significantly different from the background rates (indicating neutral evolution)—following interpretation in (). Stasis is the most frequent pattern, broken almost exclusively by archelosaurs toward positive selection.

The time frame of such changes coincides with the marked environmental changes happening at the end of the Permian, responsible for the EGE and PTE, and in the Early to Middle Triassic (Fig. 4). Therefore, our findings suggest that a shift from evolutionary stasis to a strongly directional evolutionary regime by archelosaurs at the end of the Permian is associated with an adaptive response to those fast climatic shifts. Combined with ecological opportunity arising from the demise of several groups of early synapsids after the EGE and PTE (, , , ), climate change–driven adaptive evolution resulted in the rapid diversification of the vast diversity of reptile morphotypes that came to characterize worldwide ecosystems later on during the Triassic ().

Two clear exceptions to this pattern are the rates of evolution in lepidosauromorphs and millerettids, with both groups always displaying lower rates relative to the base rate (≤1) throughout the late Permian to Middle Triassic (Figs. 3 and 4A and figs. S12 and S13). This strong signal of phenotypic stasis and stabilizing selection in the face of climatic instability contrasts with the pattern observed for most other reptile groups living at the same time interval, suggesting a radically different adaptive strategy.

Body size and its impacts on the adaptive response to climate change

A likely explanation for the detected difference in evolutionary rates and mode of selection among early reptiles (i.e., archelosaurs versus lepidosauromorphs) relates to differences in body size. The phenomenon, known as the Lilliput effect (), has been suggested as the observed pattern of marked reduction in body sizes across the PTE in both terrestrial and marine faunas (, –), including lepidosauromorph-bearing sites (), although rarely quantified for terrestrial ecosystems. Small body sizes among ectothermic animals are expected to be favored by natural selection in the aftermath of environmental changes characterized by global warming, as smaller area-volume ratios make them better capable of heat exchange with the surrounding environment, especially at lower latitudes where average temperatures are higher compared to higher latitudes (, , ). Therefore, it is expected that smaller terrestrial animals tend to cope better with temperatures changes that can be lethal to larger bodied forms, thus being selected to maintain their preadapted morphotype—i.e., undergo stabilizing selection.

To test this hypothesis, we examined the body size variation across reptile clades and various paleolatitudes between the late Permian to Late Triassic (data file S1). Our results indicate a clear pattern of much greater disparity and mean body sizes in Archelosauria (especially as represented by archosauriforms, rhynchosaurs, and allokotosaurs) relative to lepidosauromorphs and millerettids (Fig. 6A, figs. S16 to S19, and tables S5 to S10). The pattern is sustained regardless of the body size measure estimate used (figs. S16 to S19). The much reduced body size of early lepidosauromorphs and millerettids compared to archelosaurs, including their presence in fossil assemblages marked by the Lilliput effect (), thus suggests that the maintenance of small body size was a key element to their survival across these environmental changes dominated by exceptionally high temperatures. The detected predominant mode of stabilizing selection in early lepidosauromorphs further indicates that several other elements of their skeletal morphology, not just overall body size, were also being adaptively constrained during this time frame.

Fig. 6.

Body size comparisons among terrestrial reptile clades and paleolatitudes.

(A) Boxplots of body size for terrestrial clades inferred as undergoing strong selective pressures (directional or stabilizing) for both phenotypic partitions (Fig. 5) between the late Permian and Late Triassic. Y axis depicts ln-transformed mandibular length as used for statistical analyses (see figs. S16 to S19 for additional plots also using femoral length, absolute values, and sampling up to the Middle Triassic). Dots, species data points; diamonds, medians; horizontal lines, means. (B) Kernel density plots depicting the distribution of body sizes for all terrestrial reptile clades by geographical region between the late Permian and Late Triassic (see also figs. S20 and S21 for additional plots also using femoral length, ln-transformed data, and sampling up to the Middle Triassic). The fields are transparent, and the darker region is the overlap between the two geozones. (C) Geographical distribution of body sizes among all terrestrial reptile clades (i.e., according to their paleolatitudes) between the late Permian and Late Triassic (see also fig. S22 for additional plots also using femoral length, ln-transformed data, and sampling up to the Middle Triassic). Dots, species data points. Whereas small-sized reptiles occur in both temperate and tropical zones, there is a maximum limit for reptile body sizes in tropical regions, creating a tropical size gap (B and C). Body size and latitude data are available in data file S1, summary statistics are available in tables S5 to S8, and pairwise Mann-Whitney tests for significant differences between major clades and paleolatitudes are reported in tables S9 and S10.

In addition, we detect that small-bodied groups were more common in tropical paleolatitudes, whereas temperate regions were able to sustain a much greater diversity of body sizes (Fig. 6B, figs. S20 and S21, and tables S7 and S8). Although relatively small body sizes predominate in both tropical and temperate regions, there seems to be a maximum limit for terrestrial amniote body sizes in the tropics between the late Permian and Late Triassic—creating a “tropical size gap” (Fig. 6C and fig. S22). Together, this provides strong quantitative evidence for the predominance of substantially smaller body sizes in clades undergoing phenotypic stasis compared to those undergoing fast rate changes and directional selection. It also indicates that large-bodied terrestrial taxa were exclusive to temperate regions, irrespective of clade membership, and therefore large body sizes required strong adaptive solutions to periods of fast global warming.

With relatively large body sizes since the Permian (data file S1), archelosaur survival and subsequent adaptive radiation in the Triassic must have required adaptive strategies to cope with the extreme temperatures during the Permian-Triassic transition. One such strategy was an overall reduction in body size among some lineages—e.g., protorosaurs, which have the smallest average body size among early archosauromorphs (Fig. 6A). A second strategy was to migrate to more temperature-stable aquatic environments (enabled by the specific heat capacity of water being orders of magnitude higher relative to land), such as represented by the evolution of marine reptiles from terrestrial ancestors during the latest Permian and Early Triassic—i.e., sauropterygians, thalattosaurs, and ichthyosauromorphs. This transition to the aquatic environment required several adaptive changes in their skeletal morphology, which is coupled with some of the fastest evolutionary rates within archelosaurs (Fig. 3 and figs. S14 and S15). Many of those lineages displayed large body sizes marked by elongated body plans (), some already established by the Early Triassic (), further suggesting that the marine environment also enabled the release of body size constraints among reptile groups (Fig. 6, B and C).

Last, a third adaptive response proposed here is the migration of large-bodied terrestrial archelosaurs (i.e., terrestrial archosauromorphs) to higher latitudes, where they could better cope with higher temperatures. Extremely hot temperatures coupled with major rates of temperature change would have made conditions in the tropics unbearable to reptiles within the range of body sizes displayed by archosauriforms, rhynchosaurs, and allokotosaurs (clade Crocopoda) between ~260 and 230 Ma ago (Figs. 2 and 3). This is exemplified by Early and Middle Triassic terrestrial vertebrate assemblages in temperate paleolatitudes—e.g., southern South America and South Africa, between 45°S and 60°S (Fig. 6C and data file S1)—in which large-bodied crocopod archosauromorphs (among other large-bodied vertebrates, such as therapsids) were abundant, whereas diminutive-sized vertebrates (including small reptiles) predominated in similarly aged tropical latitudes (Fig. 6C and data file S1)—e.g., Vellberg, in Germany () and Czatkowice, in Poland (). Migration to new habitats and new environmental conditions by geographically scattered reptile populations (survivors of the late Permian extinctions) would be expected to promote peripatric speciation concomitant with fast rates of phenotypic change (, ). Increased clade-wide lineage diversification () and phenotypic rates (Fig. 3) (, ) characterize the evolution of various groups of crocopod archosauromorphs during the Early and Middle Triassic. However, the scarcity of data for late Permian crocopods currently limits our ability to test this biogeographic hypothesis at the present time.

Long-term organismal response to climate change

Ecological opportunity provided by the demise of dominant terrestrial faunas (e.g., early synapsids) during the PTE at 251.9 Ma ago has been previously recognized as the primary driver of reptile adaptive radiation in the Triassic (, , , ). However, our analyses indicate that a rapidly changing climate during the late Permian to Middle Triassic also greatly affected reptile evolutionary dynamics. In particular, high rates of temperature change coupled with higher global mean temperatures during the early stages of the Triassic triggered a long-term adaptive response by reptiles that was sustained up to at least 15 Ma after the PTE. These included accelerated rates of morphological evolution among large-bodied archosauromorph reptiles, invasion of the marine realm by ichthyosauromorphs and sauropterygians, as well as maintenance of a small-bodied morphotype in lepidosauromorphs.

Specifically among terrestrial archosauromorphs, their pattern of body size and paleolatitudinal distribution indicates the preferential survival of large-bodied taxa in temperate regions (possibly enhanced by migratory events from lower latitudes) and their adaptive radiation in higher latitudes. However, the observed trends between global temperature fluctuations and evolutionary rates do not appear to be triggered by the PTE; it can be traced all the way back to the SAC event in the early Permian at 294 Ma ago (Fig. 4) or 42 Ma before the end-Permian extinction. This coincides, at least in part, with the time of origin of the major reptile groups (e.g., turtles, marine reptiles, archosauriforms, and lepidosauromorphs) in the middle to late Permian (Fig. 1). Together, these results suggest a previously unrecognized association between long-term patterns of global mean temperature and its downstream impact on global climate and evolutionary response by terrestrial faunas, including the origin of new clades, construction of new body plans, and geographical/latitudinal species distributions.

Furthermore, our findings shed light on the macroevolutionary processes responsible for the origin and diversification of major clades. The emergence of the uniquely diverse reptile fauna that is characteristic of the Triassic, from marine reptiles to chameleon-like drepanosaurs, true lizards, and the forerunners of dinosaurs, is not just the product of ecological opportunity or evolutionary serendipity after the end-Permian extinction events as previously thought (, , , ). Rather, their unique morphotypes also represent the result of a sustained long-term adaptive response to an exceptionally long period of climatic shifts that lasted for ca. 57 Ma, one of the longest in Earth’s history.

MATERIALS AND METHODS

Morphological phylogenetic dataset

A new morphological dataset was constructed to analyze the relationships among early amniotes, early synapsids, and sauropsids. We sampled between 3 and 10 species for each major evolutionary lineage, focusing on the oldest and most completely known species, totaling 125 species. All species were scored on the basis of personal observation of the most completely known specimens for each taxon, except for two taxa. To our knowledge, this is the first large-scale study of early amniote and reptile relationships, with sampling to the species level based on nearly 100% first-hand anatomical observations, which, in our experience, improves the amount and precision of data scoring for morphological phylogenetics.

Morphological characters were scored on the basis of the first-hand observations of organismal morphology, combined with a wide variety of literature sources including previous phylogenetic studies on early synapsids and early reptiles. All characters were constructed following guidelines given in a recent summary of best practices of character construction (). We also included 43 new characters based on the morphological variation that we observed among the studied specimens. The final dataset includes 348 morphological characters that were evenly sampled across all studied groups. We included 24 direct autapomorphic characters—sensu (), characters in which one state is observed in only one taxon, in the case of binary characters—and 17 characters expected to behave as indirect (homoplastic) autapomorphies (independently acquired by distantly related terminal taxa). Both types of autapomorphic characters, besides character state reversals in terminal taxa, are important to contribute to branch lengths at the tips, making branch lengths and evolutionary rates at both deep and shallow internal branches comparable to each other (). Accordingly, in our final results, we observed a considerable number of changes on terminal branches comparable to the number of changes observed among internal branches (data file S2).

Body size data

Most groups of early amniotes and reptiles represented here are entirely extinct and do not have well-preserved fossils that permit direct measures of body size. Therefore, we used two separate linear measurements to approximate their body size: one based on cranial data (mandibular length) and one postcranial measure (femoral length). As observed in our results (e.g., figs. S16 to S19), trends in body size difference among clades remain the same regardless of body measure, and we depict mandibular length herein (Fig. 6A), as it has the least amount of missing data.

We obtained digital measurements using specimen photographs with embedded scale bars for size calibration. Most of those measurements were conducted using high-resolution photographs taken by one of us (T.R.S.), whereas measurements for other specimens were obtained from their published images—all measurements and literature sources are provided in data file S1. Measurements were restricted to clades entirely (or mostly) terrestrial, thus excluding all clades composed by taxa living primarily in aquatic habitats for two primary reasons: (i) Tetrapods secondarily adapted to aquatic habitats are generally more anteroposteriorly elongated and have reduced limb proportions, making estimates of body size based either on mandibular length or femoral length not directly comparable to terrestrial clades (); (ii) clades primarily inhabiting aquatic habitats live in a more temperature-stable environment compared to terrestrial clades; and (iii) physical constraints linked to the viscosity of the medium, in this case water, and the reduced effects of gravity on phenotype due to buoyancy. In combination, these changes enabled the release of size constraints among aquatically adapted clades, resulting into a great disparity of body sizes in nearly all marine reptiles (, ).

The total sample size includes 90 specimens from 62 species, representing nine major terrestrial early reptile clades between the middle Permian and Late Triassic. Some of these species were not included in the phylogeny, as they were not personally examined by us and the size data were obtained from the literature alone. The final body size sample used for the analyses here include only species in clades undergoing strong selection (indicated in data file S1), totaling 51 species (data file S5). This final number represents the maximum number of species for which we could assign either the mandibular or femoral data based on the available literature and personal data collection, representing the largest available dataset for reptile body sizes between the late Permian and Late Triassic that we are aware of. However, we expect to continue to expand this dataset in the future for a more extensive sampling of some clades in the Late Triassic.

For species with two or more specimens with body size measures available, we calculated the species average to be used for subsequent analyses (BodySize_by Species; data file S1). For the final sample size used for each analysis, see tables S9 and S10.

Paleotemperature and paleolatitudinal data

We used a recently published large temperature dataset of sea surface temperatures (SSTs) based on oxygen isotopes for the last 450 Ma (). We chose this dataset because of its vast data sampling across different latitudes—and so mean temperatures reflect global mean estimates across tropical and temperate zones (as also sampled for reptile taxa)—and because it provides a uniform and comparable temperature measure across broad geological time scales. To our knowledge, there are no equivalent datasets measuring the atmospheric temperature between the Permian and Triassic, but these are expected to correlate closely with SST.

Paleolatitudinal data were obtained from (), which provides an up-to-date estimate of paleolatitudes based on the different paleogeographical projections and how these shifted across geological time. When estimates were not available for the location of some fossil species, these were complemented with the paleolatitudinal data available in the Paleobiology Database (paleobiodb.org) for these same species.

Morphological character partition

We conducted analyses in which all phylogenetic morphological characters are treated as a single partition and a set of analyses in which all characters were subdivided into morphological partitions. Partitions were initially detected following an automatic partitioning approach using partitioning around medoids (PAM), with the optimal partition number was set using the silhouette index (Si), followed by a graphic clustering approach [t-distributed stochastic neighbor embedding (t-SNEs)] to verify results by PAM + Si—as described in () and recently implemented in the R package EvoPhylo v.0.2 (). PAM + Si indicate the utilization of two partitions as the best, but the distribution of characters among those two partitions did not have a strong fit with the partitions detected by t-SNEs. Therefore, we attempted a second partitioning scheme, in which partitions were identified on the basis of anatomical regions, but limited to two partitions only following the results by PAM + Si (i.e., cranial and postcranial partitions). The latter had a much better agreement with the t-SNE–identified partitions (fig. S23), and thus, this was our partitioning scheme of choice for partitioned clock analyses.

For the latter, we assigned a separate morphological clock to each partition, thus allowing us to obtain separate estimates of character evolution for each partition (e.g., cranial versus postcranial partitions). This strategy is especially relevant to address how different body regions had distinct adaptive responses to a changing climate.

Character and clock models

We assessed the fit of different character and clock model combinations to our data using the stepping-stone sampling strategy to assess the marginal model likelihoods () that are used to calculate Bayes factors (BFs)—35 steps between 100 and 150 million generations, depending on the specific test. For the model of character evolution, we tested for the best-fitting probability distribution (gamma or lognormal) to model among-character rate variation (ACRV) in our dataset. For clock models, we tested between two classes of relaxed clock models: the uncorrelated independent gamma rate, and the autocorrelated continuous autocorrelated (TK02)—clock models, available in MrBayes (). Last, we tested the best-fitting distribution for the variance on the clock rate prior: An exponential distribution with mean values = 1.0 and 10 (default). Using the significance thresholds in (), we found a mild preference for a gamma ACRV model (BF = 4.54), a very strong support (all BFs >80) for the autocorrelated clock model (TK02), and a strong preference for the prior on clock variance following an exponential distribution with mean = 1.0 (BF = 42.36) (table S11). Those values were used for all clock analyses.

The starting value for the prior on the clock rate was given an informative prior as per the previous nonclock analysis—the median value for tree height in substitutions from posterior trees divided by the age of the tree based on the median of the distribution for the root prior (11.63/322.5 = 0.036). The mean of the lognormal distribution was given the value based on the nonclock tree estimate in a natural log scale: ln(0.036) = −3.3242. Last, we chose a broad SD around the mean (σ = 1.0).

Age calibrations

The minimum age of the root was given the maximum age for the oldest known fossil amniotes, which come from the Moscovian, middle Pennsylvanian (~315 Ma ago) (). The maximum age for crown amniotes (under a node-based definition: the split between synapsids and sauropsids) has been suggested to be at ca. 332 Ma ago on the basis of the age of the oldest known putative stem amniote fossils, Westlothiana (, ). More recent divergence times based on genome-wide datasets are somewhat younger than the previous, estimating a mean age for the crown amniote node of 316 Ma ago (95% CI = 299 to 330.15 Ma ago) () and 319.4 Ma ago (95% CI = 317.2 to 322.6 Ma ago) (). Therefore, to provide a conservative maximum possible age for the root, we used the inferred previous estimates for the soft maximum age for the crown amniote node (330 Ma ago). We used an offset exponential distribution for the prior on root age, which gives a higher sampling probability for values closer to the minimum age and a relatively low (but nonzero) probability of age values higher than 330 Ma ago (95% quantile = 360 Ma ago).

For initial analyses, all fossil calibrations (apart from the root node) were based on tip dates only, which account for the uncertainty in the placement of fossil taxa and avoids the issue of assigning maximum age constraints that are necessary for node-based age calibrations (). In addition, we avoided recently detected strong biases that can be introduced by point age calibrations on the age of the fossils () by using a uniform prior distribution on the age range of the stratigraphic occurrence of the fossils (data file S1). We implemented the fossilized birth-death (FBD) process for the tree model with a flat prior (uniform distribution: samples ∈ U[0,1]) on the parameters of relative extinction (turnover) and probability of sampling fossils, whereas the net diversification rate was sampled from an exponential distribution with mean = 10 (default in MrBayes).

Tree models

We tested the impact on our results from various tree modeling conditions by implementing both the constant rates FBD tree model and the skyline variant of the FBD model (SFBD) (). The SFBD tree model enables a more biologically realistic parameterization of the diversification process, as it enables shifts across net diversification, relative extinction, and relative fossilization parameters across the tree (, ). This is especially useful in the case of higher-level phylogenies with taxa sampled across the span of several million years, in which diversification dynamics are expected to be affected by global events such as mass extinctions and major shift in global temperature, among others. Incorporating the SFBD model in trees rich in fossil data across the temporal scales thus makes tree models more realistic, with inferred improvements on divergence time estimates for empirical datasets ().

Diagnostic parameters

Convergence of independent runs was assessed using the average SD of split frequencies (~0.01), potential scale reduction factors (≈1 for all parameters), and effective sample size for each parameter greater than 200 and was analyzed using Tracer v. 1.7.1 (). Furthermore, the evolutionary rate parameter estimates were contrasted against the branch duration to check for time scaling biases, which are expected to affect deterministic evolutionary rate estimates (, ). As expected from our probabilistic-based approach (), we found no direct correlation between branch duration and rate value (figs. S24 to S26). Analyses were conducted using the CIPRES Science Gateway v.3.3 ().

Strength of selection

Detecting the strength (and mode) of natural selection using evolutionary rates from phenotypic data has been conducted for decades—e.g., (, , ). When rates are significantly accelerated, it provides support for positive phenotypic selection in analogy with the ratio of non-synonymous to synonymous substitutions (dN/dS) in molecular evolution, whereas strongly decelerating rates represent an instance of stabilizing selection, stasis, or constraint (, ). Here, we used our adaptation of this principle to clock-based evolutionary rates () recently implemented in the R package EvoPhylo v.0.2 (), which we attribute the strength of selection upon distinct morphological partitions for each evolutionary lineage based on the number of SDs that rates for each partition and lineage (Δv) deviate from the mean background rate (μΔb). When the mean rate of evolution on a branch is greater than the background rate plus 1 SD (Δv > μΔb + σ), it is interpreted as support for positive selection. When the mean rate of evolution on a branch is less than the main background rate minus 1 SD (Δv > μΔb − σ), it is interpreted as support for stabilizing selection or stasis. Deviations higher than 3 SDs (>3σ) are considered to indicate strong support of adaptive evolution (either for positive or stabilizing selection). For further details, see Simões et al. ().

Macroevolutionary parameters

Macroevolutionary parameter inferences from each node of the posterior trees were extracted from output log files and imported into R to estimate summary statistics and CIs. Parameters depicted in the summary maximum compatible tree (MCT) were imported from tree files into R, where parameter statistics were completed using custom scripts and functions from the R package EvoPhylo v.0.2 ()—see data file S5. In all cases, we discarded the parameter values extracted from nodes of post–Late Triassic age (used here only for evolutionary tree inference phase—see above).

Statistical analysis

We performed pairwise statistical tests to assess whether the distribution of evolutionary rates across cranial and postcranial partitions and body sizes between time bins and latitudinal bins were significantly different for each major clade assessed here. The results for most clades failed tests of normality of data distribution (Shapiro-Wilk normality test and visual assessment of data distribution; Bartlett test of homogeneity of variances to assess homoscedasticity in the data). Therefore, we conducted nonparametric pairwise Wilcoxon rank sum (Mann-Whitney) for all analyses (reported in tables S9 and S10). All statistical tests used a significance level (α) = 0.05.