Non-technical Summary
Fossils give us a rare look into how life on Earth has changed over hundreds of millions of years, showing when species appeared, thrived, and disappeared. Using recently developed mathematical modeling, we studied an updated dataset of three groups of early relatives of mammals: Ophiacodontidae, Edaphosauridae, and Sphenacodontidae. Our results confirm that these animals diversified during most of the late Carboniferous, but waned thereafter. However, our new results differ in several major points compared with our previous study (published in 2024). Notably, the slowdown in their diversification seems to have started earlier than we previously suggested, at the end of the Carboniferous, rather than early in the Permian, and it appears to have coincided with a moderate (53% survival probability), previously unreported extinction event. Moreover, this group seems to have experienced a much more severe crisis in the mid-Kungurian, which left very few surviving lineages. This pattern may be explained by a climate instability that started about 299 million years ago and lasted throughout the early Permian. This was likely caused by intense volcanism that produced large lava and volcanic ash deposits. The results seem to have been a stagnating biodiversity between two brief extinction events. This is reminiscent of the “press–pulse model,” in which a prolonged buildup of stress (the press) is followed by a sudden, intense collapse (the pulse), but it differs by starting with a crisis. Thus, these early mammal relatives may have experienced a more complex series of events, consisting of two brief crises separated by a long period of stress.
Introduction
The Permian (298.9–251.9 Ma) witnessed a spectacular diversification of amniotes and of other terrestrial life forms (Laurin Reference Laurin2010). The synapsids (mammal total group) dominated the terrestrial faunas until the end of the Permian, when the largest mass extinction event (i.e., the extinction of many lineages over a short period) took place (e.g., Irmis and Whiteside Reference Irmis and Whiteside2012; Viglietti et al. Reference Viglietti, Benson, Smith, Botha, Kammerer, Skosan, Butler, Crean, Eloff and Kaal2021). After that crisis, reptiles would dominate the biosphere for the entire Mesozoic (Carroll Reference Carroll1988). However, important taxonomic turnover occurred in synapsids well before the end of the Permian. This may have result in part from an important climate change that occurred around the Carboniferous/Permian boundary. It consisted of global warming and aridification, as the Late Paleozoic Ice Age (LPIA) was waning (Montañez and Poulsen Reference Montañez and Poulsen2013; Richey et al. Reference Richey, Montañez, Goddéris, Looy, Griffis and DiMichele2020). In addition, up to four mass extinction events have been recognized in the Permian alone (Lucas Reference Lucas2017; Didier and Laurin Reference Didier and Laurin2021). In chronological order, the first one is in the Sakmarian or the Artinskian (Benton Reference Benton1985, Reference Benton1989; Brocklehurst et al. Reference Brocklehurst, Kammerer and Fröbisch2013), between 293.5 and 283.3 Ma. The second, near the end of the Kungurian (Lucas Reference Lucas2017; Brocklehurst Reference Brocklehurst2018), is often called “Olson’s extinction” (Sahney and Benton Reference Sahney and Benton2008; Brocklehurst Reference Brocklehurst2018, Reference Brocklehurst2020), and occurred around 274.4 Ma. The third, near the end of the Capitanian (Day et al. Reference Day, Ramezani, Bowring, Sadler, Erwin, Abdala and Rubidge2015; Arefifard and Payne Reference Arefifard and Payne2020; Rampino and Shen Reference Rampino and Shen2021; Fielding et al. Reference Fielding, Bryan, Crowley, Frank, Hren, Mays, McLoughlin, Shen, Wagner and Winguth2023), is sometimes called the “dinocephalian extinction event” (Lucas Reference Lucas2009) and dates from 259.5 Ma. The famous end-Permian event (Ward et al. Reference Ward, Botha, Buick, De Kock, Erwin, Garrisson, Kirschvink and Smith2005; Smith et al. Reference Smith, Rubidge, Van der Walt and Chinsamy-Turan2012; Viglietti et al. Reference Viglietti, Benson, Smith, Botha, Kammerer, Skosan, Butler, Crean, Eloff and Kaal2021), which was probably the most severe, occurred around 252 Ma. Recovery from that event may have extended into the Middle Triassic (Irmis and Whiteside Reference Irmis and Whiteside2012), possibly because other crises took place in the Early Triassic (Hochuli et al. Reference Hochuli, Sanson-Barrera, Schneebeli-Hermann and Bucher2016).
Of these four Permian mass extinction events, the one that resulted in the most spectacular faunal changes among synapsids is arguably the Kungurian event. It witnessed the demise of most lineages that originated during the first (Carboniferous and Cisuralian) evolutionary radiation of synapsids. These were ectothermic taxa (Faure-Brac and Cubo Reference Faure-Brac and Cubo2020; de Buffrénil et al. Reference de Buffrénil, de Ricqlès, Zylberberg, Padian, Laurin and Quilhac2021) with a sprawling posture (Gônet et al. Reference Gônet, Bardin, Girondot, Hutchinson and Laurin2023) that bore little resemblance to mammals. These include five main clades, which we list in their order of inferred appearance, based on the established phylogeny (Reisz Reference Reisz and Wellnhofer1986; Huttenlocker et al. Reference Huttenlocker, Singh, Henrici and Sumida2021). First, Caseasauria produced some fairly large herbivorous taxa. Second, varanopids included agile, typically midsized carnivores. Third, ophiacodontids were large predators and were long thought to have been amphibious. Fourth, the sail-backed edaphosaurids included the first known herbivorous amniotes. Fifth, the sphenacodontids, most of whom also had a sail, were the top predators on land (Romer and Price Reference Romer and Price1940; Reisz Reference Reisz and Wellnhofer1986). These synapsid taxa did not inhabit exactly the same habitats. Ophiacodontids, edaphosaurids, and to a lesser extent the sphenacodontids (collectively called “OES grade” hereafter) are believed to have inhabited lowlands (Romer and Price Reference Romer and Price1940; Olson Reference Olson1975). This is shown in reconstructions of these animals in their presumed habitats (Fig. 1), although this does not necessarily imply that they spent much time swimming (Felice and Angielczyk Reference Felice, Angielczyk, Kammerer, Angielczyk and Fröbisch2014; Laurin and de Buffrénil Reference Laurin and de Buffrénil2016). Caseasauria and Varanopidae inhabited more upland environments (Brink et al. Reference Brink, MacDougall and Reisz2019), a little farther from large bodies of standing water.
Two typical lowland limbed vertebrates from the late Pennsylvanian and Cisuralian: Ophiacodon mirus holding a prey (an embolomere, an aquatic to amphibious stegocephalian) in its mouth. Reconstruction by Dmitry Bogdanov (a cardiologist and paleo-artist living in Chelyabinsk, Russia, who kindly gave us permission to reproduce it) sent to us by Amin Khaleghparast (a biologist from Tehran, Iran). Lower-resolution versions of this image are available in the Wikimedia commons at: https://commons.wikimedia.org/wiki/File:Ophiacodon_mirus.jpg.
Figure 1. Long description
At the center foreground, Ophiacodon mirus is depicted in lateral view with robust limbs splayed to the sides, a broad head, and a long tail curving rightward. Its jaws are open, gripping an elongate, greenish embolomere whose body arches downward and left, with visible fin-like limbs and a laterally compressed tail. The Ophiacodon’s skin is textured with subtle striping and shading, and its claws are spread on the substrate. The background consists of horizontally layered, dark and olive logs or sediment, suggesting a swampy lowland environment. No text labels are present. The predation interaction is the focal point, with the prey’s head and upper body inside the predator’s mouth.
A limited diversity of caseasaurs and varanopids extended well into the Guadalupian, after the Cisuralian. In contrast, ophiacodontids and edaphosaurids are unknown in the fossil record after the Kungurian, and the last sphenacodontid occurs barely later, in the early Roadian (Laurin and Hook Reference Laurin and Hook2022). Thus, by the Guadalupian, the synapsid faunas were dominated by the endothermic therapsids (Faure-Brac and Cubo Reference Faure-Brac and Cubo2020). These taxa had a more erect posture than earlier synapsids (Blob Reference Blob2001), and are now represented by mammals (Sidor Reference Sidor2001; Angielczyk and Kammerer Reference Angielczyk, Kammerer, Zachos and Asher2018). This brief taxonomic overview of early synapsids follows the established view that varanopids are synapsids. This long-established idea (Romer and Price Reference Romer and Price1940; Didier and Laurin Reference Didier and Laurin2024) has been challenged by a few (Laurin and Piñeiro Reference Laurin and Piñeiro2018; MacDougall et al. Reference MacDougall, Modesto, Brocklehurst, Verrière, Reisz and Fröbisch2018; Ford and Benson Reference Ford and Benson2020), but not all (Benoit et al. Reference Benoit, Ford, Miyamae and Ruf2021; Simões et al. Reference Simões, Kammerer, Caldwell and Pierce2022; Jenkins et al. Reference Jenkins, Benson, Ford, Browning, Fernandez, Dollman, Gomes, Griffiths, Choiniere and Peecook2024) recent studies.
Important gaps remain in our knowledge of how the ectothermic synapsids became extinct and how they were replaced by therapsids. Most notably, the fossil record of continental vertebrates is scant in the late Kungurian and early Roadian (Laurin and Hook Reference Laurin and Hook2022). That record is also geographically disjunct, being best represented by the U.S. Southwest (especially Texas, New Mexico, and Oklahoma) in the Kungurian (Olson Reference Olson1962, Reference Olson1967), and in the Russian Platform in the Roadian (Golubev Reference Golubev2015). A controversy surrounding the age of the youngest North American Permian lithological units that yielded continental vertebrates has existed since Olson’s (Reference Olson1962, Reference Olson1967) pioneering works, but it seems to have been settled recently (Lozovsky Reference Lozovsky2005; Brocklehurst Reference Brocklehurst2020; Laurin and Hook Reference Laurin and Hook2022). Recent works (e.g., Forte et al. Reference Forte, Kustatscher, van Konijnenburg-van Cittert, Looy and Kerp2017; Spina et al. Reference Spina, Marchetti, Diez, Capezzuoli, Saleh, Barreiro, Kustatscher and Ronchi2026) showed that most of the available paleobotanic evidence that had been used by Laurin and Hook (Reference Laurin and Hook2022) is irrelevant to this issue, because it is compatible with a Kungurian or Roadian age (among others) of the San Angelo and Chickasha formations. However paleobotanical data formed only a minor part of the evidence provided by Laurin and Hook (Reference Laurin and Hook2022). The early Roadian age advocated by that study was based primarily on vertebrate biostratigraphy and lithostratigraphy, in addition to marine fossils (fusulinids and ammonoids).
We undertook a series of studies to reassess the pattern of biodiversity evolution of early synapsids to better understand the demise of ectothermic synapsids (Didier and Laurin Reference Didier and Laurin2021, Reference Didier and Laurin2024). Our last study on the evolution of biodiversity of early synapsids found no evidence of mass extinction events in the OES grade, but we did find a strong drop in cladogenesis rate in the early Cisuralian. This suggested that the OES grade simply declined gradually and did not support a Kungurian mass extinction event (Didier and Laurin Reference Didier and Laurin2024). Our study on Varanopidae suggested that this clade became extinct in the Capitanian event, which is unsurprising, because the last varanopids occur only a little before the suspected extinction horizon (Laurin and Didier Reference Laurin and Didier2025). It also suggested a previously unrecognized extinction event near the Carboniferous/Permian boundary. This is surprising, because according to Stanley and Yang (Reference Stanley and Yang1994), that temporal interval featured few extinctions.
Here, we assess the dynamics of decline and extinction in three clades of eupelycosaurs: Ophiacodontidae, Edaphosauridae, and Sphenacodontidae. We also compare our results with our recent findings on varanopids (Laurin and Didier Reference Laurin and Didier2025). The fact that varanopids appear to have inhabited more upland environments (Brink et al. Reference Brink, MacDougall and Reisz2019) raises the possibility that differences in preferred habitat between these taxa might explain differences in their diversification dynamics.
Moreover, we compare our results with those that we previously obtained on the OES grade, using an earlier implementation of the fossilized birth–death model (FBD) (Didier and Laurin Reference Didier and Laurin2024) in which the limits between time slices had to be placed a priori rather than estimated. The current version of our software to assess the skyline FBD (introduced in Laurin and Didier Reference Laurin and Didier2025) estimates the age of the boundaries between time slices, in addition to estimating the intensity of the crises. This difference in FBD implementation (see “Methods” for other, more technical innovations) yields new results about the evolution of biodiversity in the OES grade. Notably, our new analyses contradict our previous conclusion (Didier and Laurin Reference Didier and Laurin2024) that the OES grade was not affected by mass extinction events. We thus confirm the existence of an extinction event around the Carboniferous/Permian boundary that we first detected among varanopids (Laurin and Didier Reference Laurin and Didier2025) and provide strong evidence that Olson’s extinction (Lucas Reference Lucas2017; Brocklehurst Reference Brocklehurst2018) had a strong impact on the OES grade in the Kungurian.
Compared with previous studies by others, our work introduces several types of improvements. First, the stratigraphic resolution of the database used is finer than stage level (Laurin and Didier Reference Laurin and Didier2025) while avoiding time binning, a strategy that can create artifacts, as explicitly acknowledged by Brocklehurst (Reference Brocklehurst2018), who used alternative strategies to minimize this problem. Second, we use phylogeny-informed analyses rather than taxon counts (e.g., Irmis and Whiteside Reference Irmis and Whiteside2012), because the latter is known to be suboptimal (Bertrand et al. Reference Bertrand, Pleijel and Rouse2006). Third, our analyses are performed at the lineage (species) level, rather than on families (e.g., Benton Reference Benton1989) or genera (e.g., Day et al. Reference Day, Ramezani, Bowring, Sadler, Erwin, Abdala and Rubidge2015). Finally, our FBD-based method was especially designed to model biodiversity evolution using the fossil record in a phylogenetic context, whereas other implementations of the FBD (e.g., Heath et al. Reference Heath, Huelsenbeck and Stadler2014) and other analytical methods do not fully use the stratigraphic and phylogenetic information that fossils can provide.
Methods
We updated the absolute age of the 175 fossil occurrences of the 50 terminal taxa considered by Didier and Laurin (Reference Didier and Laurin2024) to reflect the geologic timescale (GTS from here on) published in September 2023 (Cohen et al. Reference Cohen, Harper, Gibbard and Car2023), because Didier and Laurin (Reference Didier and Laurin2024) had used an earlier timescale published by Gradstein and Ogg (Reference Gradstein, Ogg, Gradstein, Ogg, Schmitz and Ogg2020). In most cases, the differences between these two GTSs are small for the relevant chronostratigraphic units, namely the interval between late Pennsylvanian and early Guadalupian. The sole exception is the Kungurian/Roadian boundary, which was placed at 274.4 Ma in the 2020 GTS, but was moved to 273.01 Ma in the 2023 GTS. These records were extracted both from the Paleobiology Database (Peters and McClennen Reference Peters and McClennen2016) and from the primary literature (see Didier and Laurin [Reference Didier and Laurin2024] for the references).
We used the same tree population sample of 100 equiparsimonious trees as Didier and Laurin (Reference Didier and Laurin2024). Analyses were carried out by using our current implementation of the skyline FBD model, a model for which other implementations are available (Gavryushkina et al. Reference Gavryushkina, Welch, Stadler and Drummond2014; Zhang et al. Reference Zhang, Stadler, Klopfstein, Heath and Ronquist2016; do Rosario Petrucci et al. Reference do Rosario Petrucci, May and Heath2025). Under this model, lineages can undergo three types of events: speciation (cladogenesis), extinction, or fossilization with respective rates λ, μ, and ψ that are constant over a given number of time slices. Our current implementation, first used in our study of varanopid diversification (Laurin and Didier Reference Laurin and Didier2025), allows us to estimate both boundary times between the slices and the rates under user-defined constraints. For instance, one may constrain the fossilization rate to remain the same over all the slices, as we do here. We test models with up to three time slices and up to two mass extinction events. Didier and Laurin (Reference Didier and Laurin2024) and older studies had tested only simpler models with no more than two time slices and a single extinction event, if any. Data required are the tree topology (without the divergence times) of the studied taxa, which may be extinct or not, and the fossil ages, generally provided as stratigraphic intervals. Our software handles missing data by dealing with the situation where the fate of certain lineages is unknown (or not considered for some reason) from a time specific to each of these lineages.
Some technical features of our software changed from the version used in Laurin and Didier (Reference Laurin and Didier2025). In particular, the likelihood optimization now proceeds in two stages (rather than one), with a global search followed by a local refinement. This two-step approach is essential, because the maximum likelihood is not computed solely over the model parameters, but also over the fossil ages, which are themselves constrained by their stratigraphic intervals. We use the NLopt implementation of the evolutionary algorithm described in da Silva Santos et al. (Reference da Silva Santos, Goncalves and Hernandez-Figueroa2010) for the global stage, followed by the Subplex algorithm for the local optimization (Rowan Reference Rowan1990).
Results
Our new analyses support the existence of two mass extinction events, one of moderate intensity at the Carboniferous/Permian boundary and a more severe one in the mid-Kungurian. In this respect, the present results differ sharply from the model favored in our previous study on the OES grade, which featured no mass extinction event (Didier and Laurin Reference Didier and Laurin2024). In addition, our results also show stronger support for the best model, which receives an Akaike information criterion (AIC) weight of 0.95 (Table 1). In our previous studies that used the skyline FBD, the best model received AIC weights of 0.65 in Didier and Laurin (Reference Didier and Laurin2024) for the OES grade and 0.64 (reduced dataset) to 0.74 (full dataset) for varanopids (Laurin and Didier Reference Laurin and Didier2025). Thus, the most-supported model here (M3EE; Fig. 2), which includes three time slices and two partial extinction events, is our central scenario. The other models need not be considered, given their low support (AIC weight ≤ 0.05). That well-supported model (M3EE; Table 1) infers boundaries between time slices at 298.9 Ma, precisely at the Carboniferous/Permian boundary (placed at 298.9 Ma in the 2023 GTS used here), and at 278.9 Ma, in the mid-Kungurian (Figs. 2, 3).
Fit of the various fossilized birth–death (FBD) models showing that the data strongly support three time slices, a mild extinction event at the end of the Gzhelian and a severe mass extinction event in the mid-Kungurian. The number in the model name indicates the number of time slices modeled; “N” indicates that no mass extinction is modeled between two time slices; “E” indicates that a mass extinction is modeled between two time slices; models with three time slices include two letters that indicate if an extinction is modeled or not between the first and second, and between the second and third time slices (respectively). The parameters for each time slice are listed in the following order: cladogenesis (“speciation”), extinction, fossilization (here assumed to be constant), and (when a mass extinction is modeled at the end of the time slice) survival rate. AIC, Akaike information criterion.
Table 1. Long description
The table consists of seven columns labeled M1, M2N, M2E, M3NN, M3NE, M3EN, and M3EE, representing fossilized birth–death models with varying time slices and extinction events. The first row lists log likelihood values from minus 692.09 for M1 to minus 647.5 for M3EE. The second row shows the number of parameters, ranging from 3 for M1 to 11 for M3EE. The third row presents A I C values, decreasing from 1738.18 for M1 to 1665.01 for M3EE. The fourth row displays A I C weights, with M3EE having the highest at 0.9495, indicating strongest model support. Estimated parameters for the first time slice are listed for each model, including cladogenesis, extinction, fossilization, and survival rate when applicable. For example, M3EE has values 0.564, 0.278, 0.295, and 0.527. Estimated time boundaries for the first and second slices are provided, such as 299.02 for M2N and 279.00 for M3NN. Parameters for the second and third time slices are detailed, with values like 0.060, 0.060, 0.295, 0.113 for M3EE in the second slice and 0.000, 0.082, 0.295 for M3EE in the third slice. The table demonstrates that models with three time slices and extinction events, especially M3EE, are best supported by the data.
Models tested and parameter estimates. Following the “M” for “model,” the number indicates how many time slices are considered. One or two letters follow for models that include two or three time slices, respectively. The letters “E” and “N” indicate the presence or not (respectively) of a mass extinction event at the time boundary between two slices. For models with two slices, there is a single such boundary. For models with three time slices, the first letter specifies whether a mass extinction is allowed (E) or not (N) between the two oldest slices, and the second letter specifies whether a mass extinction event is allowed (E) or not (N) between the two youngest slices. Mass extinction events are represented with thick vertical lines. The best-supported model, with an Akaike information criterion (AIC) weight of 0.95, is model M3EE (Table 1).
Figure 2. Long description
From top to bottom, seven horizontal panels labeled M1, M2N, M2E, M3NN, M3NE, M3EN, and M3EE display model structures. Each panel contains one to three rectangular boxes, each with parameter values: lambda, mu, psi. For models with multiple time slices, boxes are connected by arrows; thick vertical lines indicate mass extinction events. For example, M2E and M3EE show thick red lines at time boundaries, denoting mass extinctions. Numeric values above arrows indicate the timing of transitions. The bottom section is a stratigraphic chart spanning 320 to 270 million years ago, divided into Carboniferous (Pennsylvanian) and Permian (Cisuralian, Guadalupian, Roadian) periods, with further subdivisions (e.g., Bashkirian, Asselian, Kungurian). The best-supported model, M3EE, includes three time slices, two mass extinction events, and parameter shifts at each boundary.
Timetree with modeled mass extinction events. This figure was generated using our computer program (https://github.com/gilles-didier/DateFBD). The probability densities of divergence times are shown in violet on one of the 100 most parsimonious trees used in our analyses; extinction time probability density for each lineage is shown in green. The fossil ages that maximize maximum likelihood under model M3EE are indicated by narrow orange vertical bars. Divergence and extinction time densities are computed from the single-slice model, without considering any mass extinction event (in which case the extinction distributions would be mixtures of a Dirac distribution at the time of the event and of another distribution similar to those plotted here—the first component of the mixture accounting for cases where the lineage is actually killed by the event, and the second for cases where it is not). The two vertical bars crossing the tree represent the mass extinction events supported by the most-supported model. Line transparency for extinction events is proportional to survival rate; hence, the most severe extinction corresponds to the darkest, most obvious line. The outlines of several taxa were downloaded from PhyloPic (https://www.phylopic.org/). These include, from top to bottom, Ophiacodon uniformis by Dmitry Bogdanov; Edaphosaurus pogonias by Matt Celeskey (Attribution 3.0 Unported); Haptodus garnettensis by Dmitry Bogdanov, vectorized by T. Michael Keesey; Tetraceratops insignis by Nobu Tamura, modified by T. Michael Keesey; Secodontosaurus obtusidens (CC0 1.0 Universal Public Domain Dedication); Sphenacodon ferox by Dmitry Bogdanov (Attribution-ShareAlike 3.0 Unported); and Dimetrodon by Alessio Ciaffi (Attribution 4.0 International). The Ophiacodon, Haptodus, Tetraceratops, and Sphenacodon illustrations were published under the Attribution-ShareAlike 3.0 Unported license.
Figure 3. Long description
Anchored at the left, the horizontal axis marks time from 320 to 260 million years ago, subdivided into Carboniferous and Permian periods with labeled stages. The main panel is a rightward-branching phylogenetic tree, with each tip labeled by taxon names in yellow-green text. Blue-violet curves along branches represent divergence time probability densities, while green curves indicate extinction time probability densities. Narrow orange vertical bars mark fossil age estimates for each lineage. Brown horizontal bars show fossil stratigraphic intervals. Two vertical red lines cross the tree, labeled as mass extinction events, with the darker line indicating a more severe event. Line transparency for extinction events varies by survival rate. Taxon silhouettes are arrayed vertically on the right, corresponding to Ophiacodon, Edaphosaurus, Haptodus, Tetraceratops, Secodontosaurus, Sphenacodon, and Dimetrodon, matching the tree's major clades: Ophiacodontidae, Edaphosauridae, and Sphenacodontidae. The legend at lower left explains the color coding for divergence time, extinction time, fossil intervals, and fossil age estimates. The geological timescale at the bottom is color-coded by period and stage, aligned with the time axis above.
The Carboniferous/Permian extinction event recovered by the best-supported model (M3EE) is mild, with a survival probability of 0.527 (Table 1). Its presence is strongly supported, to the extent that the best-supported model without that event (M3NE) has a low AIC weight of 0.050. That moderately severe crisis also marks a distinct shift in diversification dynamics in the OES grade. Before that crisis, in the first modeled time slice, the difference (r) between cladogenesis (0.564) and extinction rates (0.278), yields a positive and fairly strong (0.286) net diversification rate. By contrast, after that crisis, in the second modeled time slice that represents most of the Cisuralian, the net diversification is nil (r = 0.000). This results from a drastic drop in the rate of cladogenesis to only 0.006, which is not fully compensated by a strong, but smaller drop in extinction rate to 0.006.
The mid-Kungurian extinction event that affected the OES grade in model M3EE appears to be much more severe (Figs. 2, 3), with a survival rate of only 0.113. It is also better supported than the Carboniferous/Permian extinction event, given that the best-supported model without the Kungurian event (M3EN) barely reaches an AIC weight of 0.0001. After that crisis, in the third time slice that represents the late Kungurian and at least the early Roadian, the background rate of biodiversity loss seems to have accelerated. This results from a sharp drop in the rate of cladogenesis (0.000) and a slight increase in the extinction rate, to 0.082 (Table 1), which yields a negative net diversification rate (r) of −0.082, as opposed to 0.000 in the previous time slice. However, given that very few lineages in the OES grade are represented in that third time slice, this estimate is poorly constrained. These taxa include minimally Dimetrodon angelensis, from the San Angelo Formation (Olson Reference Olson1962), which probably dates from the early Roadian but could conceivably extend into the latest Kungurian (Laurin and Hook Reference Laurin and Hook2022). However, several other taxa from the OES grade may also extend into this time slice (Fig. 3). These additional taxa potentially include lineages whose latest known fossil occurrences fall entirely or partly (given the uncertainty over their exact age) in that time slice: six species of Dimetrodon (D. grandis, D. giganhomogenes, D. loomisi, D. limbatus, D. macrospondylus, and D. dollovianus), two species of Edaphosaurus (E. cruciger and E. pogonias), and three species of Ophiacodon (O. retroversus, O. uniformis, and O. major). Other OES grade taxa that may extend into this third time slice are those whose latest occurrence predates that time a bit, but whose extinction time probability densities (Fig. 3) do not allow us to rule out the possibility that they extended into that time slice. Among these we find the following taxa: Varanosaurus acutirostris, Secodontosaurus obtusidens, Dimetrodon natalis, and Tetraceratops insignis. The latest records of these taxa probably date from between 283 and 280 Ma, but they may well have become extinct after 278.9 Ma.
Discussion
Comparisons with Previous FBD-based Studies of Early Synapsid Diversification
Several differences with our previous study on biodiversity evolution in the OES grade can be explained by the fact that we (Didier and Laurin Reference Didier and Laurin2024) previously only compared support for models with one or two time slices and at most a single extinction event between the two time slices. Furthermore, the position of boundaries between time slices had to be determined a priori, and these boundaries were placed between successive geological stages. Thus, in Didier and Laurin (Reference Didier and Laurin2024), the best-supported model (with an AIC weight of 0.65) had a limit between time slices at the Asselian/Sakmarian boundary, one geological stage later than found here. However, the second best-supported model in that study, with an AIC weight of 0.22 (about a third of the weight of the best-supported model) placed the limit between the two time slices at the Carboniferous/Permian boundary (298.9 Ma), as in our present estimate. This congruence is better than expected, given the age uncertainty that is associated with all fossil ages (the average range of possible ages of the 175 fossil occurrences is about 3.6 Myr). Similarly, the age and the exact stratigraphic position of the Carboniferous/Permian boundary in continental sediments has long been debated and remains uncertain, including in the basins from which the bulk of our OES data originates (Utting and Lucas Reference Utting and Lucas2010).
The boundary between the second and third time slices, in the mid-Kungurian (278.9 Ma), is marked by a fairly severe mass extinction event (survival probability of 0.113) in the OES grade (Figs. 2, 3). That extinction event could not be detected for methodological reasons in our previous study on varanopids (Laurin and Didier Reference Laurin and Didier2025) or on the OES grade (Didier and Laurin Reference Didier and Laurin2024). However, the existence of an extinction event around that time was first inferred long ago (Olson Reference Olson, Silver and Schultz1982). It was subsequently nicknamed “Olson’s extinction” (Sahney and Benton Reference Sahney and Benton2008), and its existence is supported by a few more recent studies (Ruta and Benton Reference Ruta and Benton2008; Brocklehurst et al. Reference Brocklehurst, Kammerer and Fröbisch2013; Brocklehurst Reference Brocklehurst2018, Reference Brocklehurst2020). Our analytical method, which simultaneously maximizes likelihood of the model and of the age of each fossil within the specified age interval, as was already done in Didier and Laurin (Reference Didier and Laurin2024), can push the ages of several fossils toward one of the specified age bounds. This happened here for several taxa (listed in the “Results”) for which the latest occurrence could date from between 279 and 277 Ma. These occurrences are here estimated to have occurred between 278.86 and 279 Ma, an interval that ends at the inferred mid-Kungurian crisis (Fig. 3). Note that this result has some independent stratigraphic support, because as described in “The Pulse,” the late Kungurian in both Texas and Oklahoma includes a large stratigraphic gap (of unknown exact duration) in the fossil record (Laurin and Hook Reference Laurin and Hook2022). That crisis may have been especially severe among synapsids, because it was partly responsible for the fairly quick, partial replacement of the presumably ectothermic Permo-Carboniferous synapsid clades by the apparently endothermic therapsids (de Ricqlès Reference de Ricqlès1972; Faure-Brac and Cubo Reference Faure-Brac and Cubo2020; de Buffrénil et al. Reference de Buffrénil, de Ricqlès, Zylberberg, Padian, Laurin and Quilhac2021).
The diversification dynamics inferred here for the OES grade shows a dramatic drop in the rate of cladogenesis and extinction between the first two time slices, as in Didier and Laurin (Reference Didier and Laurin2024). However, we now infer two mass extinction events (see “Results”) and a nil net diversification rate in the second time slice, rather than the slow decline inferred in our previous work. The discrepancy between these studies is due partly to the estimated boundaries (rather than fixed a priori) between time slices, which could explain the detection of a mid-Kungurian crisis. Didier and Laurin (Reference Didier and Laurin2024) could not test that hypothesis, because extinction events were constrained to occur at stage boundaries. Another important methodological difference is the additional assumption made here that the fossilization rate is constant over the entire diversification period, as in Laurin and Didier (Reference Laurin and Didier2025), whereas it was allowed to vary between time slices in Didier and Laurin (Reference Didier and Laurin2024). This assumption is motivated by the difficulties encountered when maximizing the likelihood while allowing both the limit ages and the fossilization rates to vary. These technical difficulties generate instability in our estimates of limit ages, at least with datasets of moderate size such as ours (50 terminal taxa and 175 fossil occurrences). Constraining the fossilization rate to remain constant resolves this issue. Despite these numerical differences, both sets of estimates result in significantly positive diversification in the first time slice. Likewise, in the second time slice, both sets of estimates indicate a significant slowdown of diversification, although both studies differ in detail. Our new estimates of stagnating biodiversity over time (Fig. 2), with an r of about 0.000, contrasts with the biodiversity decline (r of about −0.078) inferred by Didier and Laurin (Reference Didier and Laurin2024), who had estimated cladogenesis and extinction rates of 0.021 and 0.099, respectively. This difference in rate of biodiversity loss between both studies is readily explained by the fact that our new estimate incorporates a severe mass extinction in the Kungurian. In contrast, the impossibility of modeling such an event in our previous study resulted in a steeper rate of biodiversity loss as an alternative way to explain why the OES grade vanished from the fossil record around the Kungurian/Roadian boundary.
Comparisons with our recent results on varanopids are much easier, because Laurin and Didier (Reference Laurin and Didier2025) used our current implementation of the FBD that estimated the boundary between time slices, and we tested for the presence of up to two mass extinction events. Only two main caveats hamper these comparisons. The first is that because of the smaller size of the analyzed varanopid datasets (16 or 21 taxa with 31 or 36 fossil occurrences, as opposed to 50 taxa and 175 fossil occurrences here), no more than two time slices were modeled. A third time slice was allowed in some models, but its parameters were set arbitrarily, because there were no varanopid fossils in that time slice. The second caveat is that the fossil record of varanopids extends up to the Capitanian, whereas Ophiacodontidae and Edaphosauridae disappear from that record in the Kungurian, and Sphenacodontidae probably extends a bit longer into the early Roadian (Laurin and Hook Reference Laurin and Hook2022). Thus, there might be genuine differences in diversification dynamics between Varanopidae and the OES grade in the Pennsylvanian and Cisuralian, and comparisons are simply impossible in the Guadalupian, given that the OES grade left a negligible fossil record in that epoch.
With these caveats in mind, we note that for Varanopidae, Laurin and Didier (Reference Laurin and Didier2025) estimated the boundary between the first two time slices close to the Carboniferous/Permian boundary, either synchronous with our current estimate (298.9 Ma) or about 2 Myr later (296.94 Ma). This was based on the two best-supported models in each varanopid datasets (full and reduced), which resulted in four estimated models. This confirms a significant change in diversification of all these eupelycosaurs around the Carboniferous/Permian boundary, if varanopids are indeed eupelycosaurs (see “Introduction”).
Likewise, both studies find support for a moderate extinction event around the Carboniferous/Permian boundary, which as far as we know, was not previously reported (at least for vertebrates) in the literature. Our inferred survival rate of 0.527 for the OES grade translates into an extinction rate of 0.473, which is less severe than in the “big five” (or six) mass extinction events. The intensity of the latter is difficult to compare with our results, because the reported extinction rates are typically expressed for taxa of higher taxonomic levels (e.g., Labandeira Reference Labandeira2005), but McGhee et al. (Reference McGhee, Clapham, Sheehan, Bottjer and Droser2013: table 3) reported that the six most severe extinction events had 39–83% extinction rates for genera. This would translate into specific extinction rates greater than 60%, according to the rarefaction curve proposed by Stanley (Reference Stanley2016: fig. 5). In any case, this OES extinction rate is apparently about as high as in Varanopidae, for which Laurin and Didier (Reference Laurin and Didier2025) estimated a survival rate of 0.747 on the exhaustive dataset (21 taxa), and only 0.358 for the reduced dataset (16 taxa). However, given the small sample size for Varanopidae, which is constrained by the small size of that clade, these numbers must be interpreted with caution. Our new analyses significantly strengthen the case for a moderate Carboniferous/Permian extinction event, given that the dataset of our OES grade dataset is about three times larger than the varanopid dataset.
Both sets of taxa seem to have stagnated or declined slowly in the Cisuralian (and, in the case of Varanopidae, until their final demise in the Capitanian). These rates of nil or negative diversification are very similar in both groups: r = 0.000 in the OES grade and r = −0.013 in Varanopidae. According to the models obtained here and by Laurin and Didier (Reference Laurin and Didier2025), the fact that varanopids may not have been affected by the mid-Kungurian crisis could explain their much greater longevity.
The mid-Kungurian crisis apparently affected the OES grade much more than Varanopidae (Figs. 2, 3). The OES grade appears to have experienced a fairly severe crisis in the mid-Kungurian (a survival rate of only 0.113 translates into an extinction rate of 0.887) and to have become extinct shortly thereafter. Indeed, the last known lineage of that grade is Dimetrodon angelensis from the latest Kungurian to (more likely) early Roadian San Angelo Formation (Olson Reference Olson1962; Laurin and Hook Reference Laurin and Hook2022). Of course, Therapsida had probably started their evolutionary radiation by then, and they are rooted in the OES grade in all recent phylogenies (e.g., Sidor Reference Sidor2001; Brocklehurst et al. Reference Brocklehurst, Ford and Benson2022; Matamales-Andreu et al. Reference Matamales-Andreu, Kammerer, Angielczyk, Simões, Mujal, Galobart and Fortuny2024), but their physiology may have differed rather sharply from that of the OES grade (de Ricqlès Reference de Ricqlès1972; Faure-Brac and Cubo Reference Faure-Brac and Cubo2020; de Buffrénil et al. Reference de Buffrénil, de Ricqlès, Zylberberg, Padian, Laurin and Quilhac2021). On the contrary, we found no evidence of an extinction event in Varanopidae in the Kungurian, although the low sample size of fossil occurrences in that clade and the methodological caveats noted earlier make that conclusion tentative (Laurin and Didier Reference Laurin and Didier2025). The absence of a Kungurian extinction event among varanopids, if genuine, might be linked to their “somewhat more upland” environment, which they shared with Caseasauria, which also survived that event.
Note that the meaning of “upland” in vertebrate paleontology (at least for Permo-Carboniferous stegocephalians) has never been rigorously defined, as far as we know. For instance, the Bromacker locality of the Tambach Formation (Sakmarian or Artinskian) experienced periodic flooding, so streams and at least ephemeral bodies of standing water were present, but no large, perennial lake or river appears to have been present. Yet Bromacker is reputedly the best example of an “upland” locality in the Cisuralian (Eberth et al. Reference Eberth, Berman, Sumida and Hopf2000). Unsurprisingly, this concept has recently been criticized by paleobotanists (Bashforth et al. Reference Bashforth, DiMichele, Eble, Falcon-Lang, Looy and Lucas2021). In this context, “upland” mostly means “farther from large, perennial water bodies” and does not necessarily imply high altitude, mountains, or even hills, even though the concept encompasses such habitats, at least potentially (Olson Reference Olson1975). Nevertheless, this difference in habitat may explain the different fates of the OES grade compared with varanopids, especially if these differences in habitat resulted in different diets.
Comparisons with Previous Studies on Diversification in the Late Paleozoic
The change in diversification dynamics detected by the three recent FBD-based diversification studies of early synapsids (Didier and Laurin Reference Didier and Laurin2024; Laurin and Didier Reference Laurin and Didier2025; and the present one) around the Carboniferous/Permian boundary had not been discussed in the literature previously, but at least one recent study of marine metazoans found a similar pattern. Namely, Fan et al. (Reference Fan, Shen, Erwin, Sadler, MacLeod, Cheng, Hou, Yang, Wang and Wang2020) described “a previously unrecognized great biodiversification event” that they called “the Carboniferous/Permian Biodiversification Event.” As in our study, their species-level biodiversity curve (their fig. 1A) shows a diversification restricted to the Carboniferous (with a slight dip near the mid-Carboniferous that may be a modest extinction event) and a sudden drop in diversity in the earliest Permian. As in our new results, the Permian data of Fan et al. (Reference Fan, Shen, Erwin, Sadler, MacLeod, Cheng, Hou, Yang, Wang and Wang2020) suggest no long-term biodiversity trend in most of the Cisuralian.
This similarity is surprising given the fact that the data used by Fan et al. (Reference Fan, Shen, Erwin, Sadler, MacLeod, Cheng, Hou, Yang, Wang and Wang2020) come from marine organisms (without vertebrate data; the diversification events described in the “Results” are especially visible among brachiopods and foraminifera) and that they were largely derived from Chinese sections. By contrast, our dataset is restricted to a fairly small number of continental amniotes, chiefly from what is now western North America and Europe, which were then parts of equatorial Pangea, rather far from the blocks that now form China (Kent and Muttoni Reference Kent and Muttoni2020). Thus, despite strong taxonomic, ecological, and geographic differences between these datasets, both studies (Fan et al. Reference Fan, Shen, Erwin, Sadler, MacLeod, Cheng, Hou, Yang, Wang and Wang2020 and the present one) point to a drop in diversity around the Carboniferous/Permian boundary. This point clearly deserves additional scrutiny.
Among the land plant communities, the Carboniferous/Permian boundary coincides with an increase in the relative abundance of xeromorphic, drought-tolerant floras alternating with hygromorphic, wetland taxa. This marks the appearance of so-called mixed floras (e.g., Falcon-Lang et al. Reference Falcon-Lang, Lucas, Kerp, Krainer, Montañez, Vachard, Chaney, Elrick, Contreras and Kurzawe2015; DiMichele et al. Reference DiMichele2014, Reference DiMichele, Hotton, Looy and Hook2019; Marchetti et al. Reference Marchetti, Forte, Kustatscher, DiMichele, Lucas, Roghi, Juncal, Hartkopf-Fröder, Krainer and Morelli2022). That floral change may well have affected insects by potentially driving faunal turnover. This is indicated by the high values of insect origination and extinction rates estimated during the Upper Pennsylvanian, followed by an early Permian minimum and a subsequent renewed increase (Labandeira Reference Labandeira2005; Nicholson et al. Reference Nicholson, Mayhew and Ross2015). Given the dominance of phytophagy, fungivory, and detritivory in Paleozoic insect communities (Rainford and Mayhew Reference Rainford and Mayhew2015), shifts in floral composition from hygromorphic to xeromorphic taxa may have acted as primary driver of the observed insect dynamics, as observed during other large-scale biotic transitions (Jouault et al. Reference Jouault, Nel, Perrichot, Legendre and Condamine2022). Indeed, the detritivore-rich insect communities of the Middle Pennsylvanian were complemented in the Upper Pennsylvanian by the emergence of new herbivorous taxa. This is shown by the first occurrence of fossilized galls on petioles and tunnel-and-gallery networks on tree ferns (Labandeira and Phillips Reference Labandeira and Phillips1996, Reference Labandeira and Phillips2002), which indicates the establishment of more complex trophic interactions. Did these floral changes negatively affect the OES grade and varanopids?
A Press–Pulse Extinction Event or a Pulse–Press–Pulse Event?
The observed biodiversity pattern is reminiscent of (but not identical to) the general pattern of the press–pulse model proposed by Arens and West (Reference Arens and West2008). The press–pulse model predicts a prolonged decline (the press) followed by a brief crisis (the pulse). The pattern that we infer for the OES grade differs slightly from that, because instead of a press, we see a stagnating biodiversity lasting about 20 Myr (from about 299 to 279 Ma). While a 20 Myr stagnation in biodiversity might not seem like a strong signal, this contrasts rather sharply with the rapid diversification that took place in the Pennsylvanian, and this suggests that the environment (biotic or abiotic) was less favorable for the OES grade and varanopids. Furthermore, Didier and Laurin (Reference Didier and Laurin2024) had inferred negative diversification, which matches the expected press. This is followed by a fairly severe extinction event (a survival probability of 11.3% translates into an extinction rate of 88.7%), which certainly qualifies as a pulse. However, we also observe a first pulse (with a survival probability of 52.7%), weaker than the second one, just before the press (Figs. 2, 3). Because of this, the pattern that we observe could better be called a “pulse–press–pulse.”
In the original proposal of this model, Arens and West (Reference Arens and West2008) envisioned the press (a slow, long-lasting biodiversity decline) to generally be caused by a phase of intense volcanism and the pulse (a brief dramatic biodiversity loss) to be caused by the impact of a large asteroid. However, Arens and West (Reference Arens and West2008) indicated that other events might cause similar patterns; for instance, climate and sea-level changes could cause the press, whereas marine anoxic incursions could cause the pulse (in the seas).
Origin of the First Pulse
A plausible cause of this mild, Carboniferous/Permian extinction event might be a warming event as recorded in δ18Oapt in conodonts from offshore sediments from South China (Chen et al. Reference Chen, Joachimski, Wang, Shen, Qi and Qie2016). The 1‰ decrease in the oxygen isotope record corresponds to a maximum temperature rise of 4°C, although the increase may have been more moderate if accompanied by ice sheet melting (Chen et al. Reference Chen, Joachimski, Wang, Shen, Qi and Qie2016). This is supported by the global sea-level curve during the Gzhelian (Haq and Schutter Reference Haq and Schutter2008), which shows a rise consistent with milder warming. This warming is also supported by vegetational changes (toward more arid-adapted land plants) that took place at the time (DiMichele et al. Reference DiMichele, Montañez, Poulsen and Tabor2009). This warming could be due to the first phase of the Tarim and Skagerrak-centered large igneous provinces (LIPs), which have been dated at ~300 Ma (Xu et al. Reference Xu, Wei, Luo, Liu and Cao2014) and 297.5 ± 3.8 Ma (Torsvik et al. Reference Torsvik, Smethurst, Burke and Steinberger2008), respectively.
The Press
In the Cisuralian, the press could plausibly be attributed to a trend toward a warmer and more arid (at least seasonally) climate. This trend could be linked to the end of the LPIA with the waning of the glaciation covering the Southern Hemisphere and an increase in atmospheric pCO2 levels from the early Cisuralian onward (e.g., Montañez and Poulsen Reference Montañez and Poulsen2013; Richey et al. Reference Richey, Montañez, Goddéris, Looy, Griffis and DiMichele2020). The transition was not continuous but rather stepwise, and plausibly diachronous from west to east (e.g., Cleal and Thomas Reference Cleal and Thomas2005; DiMichele et al. Reference DiMichele, Montañez, Poulsen and Tabor2009, 2014; Richey et al. Reference Richey, Montañez, Goddéris, Looy, Griffis and DiMichele2020). It probably lasted a long time; its earliest manifestation may be visible in the Kasimovian, a stage for which we have scant data on amniote diversity, probably because amniote diversification had begun only shortly before (Fig. 3). Sahney et al. (Reference Sahney, Benton and Falcon-Lang2010) argued that a rainforest collapse occurred at that time, resulting in dietary regime shifts among stegocephalians, and that this accelerated diversification, especially among amniotes. The initially high diversification of the OES grade that we detect here (as in Didier and Laurin Reference Didier and Laurin2024) is consistent with this interpretation. This climatic trend appears to have affected areas inhabited by the OES grade, because a flora adapted to a seasonally dry climate was reported in the Asselian of what is now New Mexico (Falcon-Lang et al. Reference Falcon-Lang, Lucas, Kerp, Krainer, Montañez, Vachard, Chaney, Elrick, Contreras and Kurzawe2015). As the climate continued becoming warmer and drier (well into the mid- to late Cisuralian), diversification of eupelycosaurs slowed dramatically. This apparently occurred rather suddenly around the end of the Carboniferous.
The coeval rise in pCO2 is possibly linked to volcanic events that took place mostly from the late Asselian to early Kungurian (Shellnutt et al. Reference Shellnutt, Bhat, Brookfield and Jahn2011; Xu et al. Reference Xu, Wei, Luo, Liu and Cao2014; Wei et al. Reference Wei, Yang, Cheng, Domeier, Wu, Kravchinsky, Zhou, Jiang, Wu and Huo2020). That activity was especially intense from the Sakmarian–Artinskian onward, with volcanic episodes such as the Tarim II, Tarim III, Panjal Traps, Qiangtang basalt, and Silicic Choiyoi igneous province (Ernst and Buchan Reference Ernst and Buchan2001; Richey et al. Reference Richey, Montañez, Goddéris, Looy, Griffis and DiMichele2020). Evidence from both western (e.g., Montañez et al. Reference Montañez, Tabor, Niemeier, Dimichele, Frank, Fielding, Isbell, Birgenheier and Rygel2007; Poulsen et al. Reference Poulsen, Pollard, Montañez and Rowley2007; Montañez and Poulsen Reference Montañez and Poulsen2013; Richey et al. Reference Richey, Montañez, Goddéris, Looy, Griffis and DiMichele2020) and central tropical Pangea low-latitude successions (e.g., Schneider et al. Reference Schneider, Körner, Roscher and Kroner2006, Reference Schneider, Lucas, Scholze, Voigt, Marchetti, Klein, Opluštil, Werneburg, Golubev and Barrick2020; Marchetti et al. Reference Marchetti, Forte, Kustatscher, DiMichele, Lucas, Roghi, Juncal, Hartkopf-Fröder, Krainer and Morelli2022) supports an intensification of the aridification trend and the following biotic turnover during the Artinskian warming episode (Forte et al. Reference Forte, Vallé and Kustatscher2023).
This aridification process is documented by pedotypes in the Appalachian Basin, as surveyed by Hembree (Reference Hembree2022) in the upper Pennsylvanian to lower Permian Conemaugh; Monongahela, and Dunkard groups in eastern North America; and ichnological and paleobotanical data from the Southern Alps, in northern Italy (Marchetti et al. Reference Marchetti, Forte, Kustatscher, DiMichele, Lucas, Roghi, Juncal, Hartkopf-Fröder, Krainer and Morelli2022). This includes plant macro- and microfossil assemblages supporting dry, seasonal, and semiarid conditions and a decline in diversity and abundance among anamniotic stegocephalians from the Artinskian–Kungurian onward (e.g., Forte et al. Reference Forte, Kustatscher, van Konijnenburg-van Cittert and Kerp2018a,Reference Forte, Kustatscher, Roghi and Pretob; Marchetti et al. Reference Marchetti, Forte, Kustatscher, DiMichele, Lucas, Roghi, Juncal, Hartkopf-Fröder, Krainer and Morelli2022; Vallé et al. Reference Vallé, Nowak, Kustatscher, Erkens, Roghi, Morelli, Krainer, Preto and Hartkopf-Fröder2023). This long-term aridification is further supported by a gradual temperature increase of a maximum of 2°C (Chen et al. Reference Chen, Joachimski, Shen, Lambert, Lai, Wang, Chen and Yuan2013) and a global decrease in river input into the oceans as recorded by strontium (87Sr/86Sr) (Korte et al. Reference Korte, Jasper, Kozur and Veizer2006; Wang et al. Reference Wang, Katchinoff, Garbelli, Immenhauser, Zheng, Zhang, Yuan, Shi, Wang and Planavsky2021) and lithium isotope (δ7Li) ratios (Wang et al. Reference Wang, Zhang, Liu, von Strandmann, Xiong, Zheng, Garbelli, Zhang, Yuan and Shen2024). A recently reported xerophytic Lopingian paleoflora from North America is also consistent with this inferred climatic trend (Hotton et al. Reference Hotton, Bercovici, Looy, Duijnstee, Chaney, Eble, Montañez, Bourquin, Cecil and Nelson2025).
The fact that we did not detect an extinction event during the Artinskian warming episode may simply be a constraint linked to the fact that we tested models with up to three time slices and two extinction events only. This constraint is imposed by the amount of data that we have (on 50 taxa and 175 fossil occurrences) and, to a lesser extent, by their associated uncertainty (phylogenetic and geochronological). However, given that our model estimated the position of the diversification dynamics shifts, our results suggest that these changes were more substantial near the Carboniferous/Permian boundary and in the Kungurian than in the Artinskian. In future studies, when we have more data (i.e., on additional amniote clades), we will be able to test more complex models, at least if these clades seem to follow similar diversification patterns.
The climatic trend marked by increasing seasonality and aridity could also reflect a northward drift of the main continental fossiliferous basins of what had been equatorial Pangea in the Cisuralian into an arid zone located at slightly higher (tropical to subtropical) paleolatitudes of both hemispheres in the Guadalupian (Kent and Muttoni Reference Kent and Muttoni2020). In that case, the inferred climatic change might be a local, rather than global, climatic trend. If so, as we recently suggested for varanopids (Laurin and Didier Reference Laurin and Didier2025), it is conceivable that representatives of the OES grade persisted in other regions that have been less intensively prospected for fossils. In fact, a possible example of such records has been documented in the Lopingian or Early Triassic Buena Vista Formation from Uruguay, which has yielded vertebrae that Piñeiro et al. (Reference Piñeiro, Verde, Ubilla and Ferigolo2003) argued belonged to a varanopid and, possibly, a sphenacodontid. At that time, the Buena Vista was located at a high, southern paleolatitude, outside the inferred arid latitudinal zones reconstructed by Kent and Muttoni (Reference Kent and Muttoni2020), so the climate should in any case have been humid enough to sustain early synapsid populations. Basal, ectothermic synapsids were not previously known in what is now South America, but such records were recently found (Angielczyk et al. Reference Angielczyk, Fröbisch, Kammerer, Smith, Marsicano, Pardo, Richter and Cisneros2026). This adds weight to the taxonomic interpretations of Piñeiro et al. (Reference Piñeiro, Verde, Ubilla and Ferigolo2003). However, the absence of the OES grade in the Guadalupian strata of two other high-latitude areas that have been intensively studied, namely the Russian Platform and the Karoo Basin, raises doubts about these putative late records of the OES grade. Caseids (Brocklehurst and Fröbisch Reference Brocklehurst and Fröbisch2017; Romano et al. Reference Romano, Brocklehurst and Fröbisch2018) and varanopids (Reisz and Berman Reference Reisz and Berman2001; Anderson and Reisz Reference Anderson and Reisz2004) are known from the middle Permian of Russia, and the Karoo Basin has yielded varanopids (Modesto et al. Reference Modesto, Smith, Campione and Reisz2011). All this suggests that the OES grade may well have been extinct by the late Wordian, or that the climate of these higher-latitude areas was unsuitable for the OES grade.
The Pulse
The cause of the pulse of extinctions in the mid-Kungurian (around 278 Ma) is probably a fast, global warming event that occurred at the end of the P2 glaciation, one of the glacial episodes of the LPIA. This mid-Kungurian interval is characterized by a global perturbation of the carbon cycle, evidenced by a prominent negative excursion (~3‰) in the carbon isotope record, hereafter referred to as the Kungurian carbon isotope excursion (KCIE; Liu et al. Reference Liu, Du, Jarochowska, Yan, Munnecke and Lu2018; Huang et al. Reference Huang, Sun, Chen, Ogg, Hou, Li, Li and Yong2025). This isotopic anomaly is observed in both inorganic and organic carbon reservoirs and across marine and terrestrial environments. Its widespread occurrences, documented in South and North China (Liu et al. Reference Liu, Jarochowska, Du, Vachard and Munnecke2017; Lu et al. Reference Lu, Zhou, Yang, Zhang, Shao and Hilton2021; Fang et al. Reference Fang, Zhang, Bai, Tang, Chao and Wei2023), Venezuela and North America (Laya et al. Reference Laya, Tucker, Gröcke and Perez-Huerta2013), South Africa (Van de Wetering et al. Reference Van de Wetering, Esterle, Golding, Rodrigues and Götz2019), and eastern Australia (Birgenheier et al. Reference Birgenheier, Frank, Fielding and Rygel2010), supports its global extent.
The KCIE coincides with peaks in Hg/total organic carbon ratios, which suggests a link to enhanced volcanic activity (Lu et al. Reference Lu, Zhou, Yang, Zhang, Shao and Hilton2021). Notably, LIP activity in Tarim III and Qiangtang (Tibet) (Zhang and Zhang Reference Zhang and Zhang2017), alongside contemporaneous arc volcanism in the Okhotsk-Taigonos and Choiyoi regions (Biakov Reference Biakov2012; Sato et al. Reference Sato, Llambías, Basei and Castro2015), has been dated to the time interval when our analyses place a mass extinction event (Fig. 4). However, carbon mass balance models indicate that volcanism alone cannot fully account for the magnitude of the KCIE (Liu et al. Reference Liu, Du, Jarochowska, Yan, Munnecke and Lu2018). One plausible supplementary source is the release of isotopically light carbon from destabilized methane hydrates during the melting of the Gondwanan ice sheet (Liu et al. Reference Liu, Du, Jarochowska, Yan, Munnecke and Lu2018; Huang et al. Reference Huang, Sun, Chen, Ogg, Hou, Li, Li and Yong2025).
One of the most recent Dimetrodon species (D. grandis, from the Kungurian Clear Fork Group; Reisz Reference Reisz and Wellnhofer1986) in front of a brightly colored sky, such as would have been produced by the intense volcanism that took place at the time. Even though much of the volcanism documented by the large igneous provinces (LIPs) was effusive, some explosive arc volcanism occurred in the Okhotsk-Taigonos and Choiyoi regions (Biakov Reference Biakov2012; Sato et al. Reference Sato, Llambías, Basei and Castro2015), and caldera eruptions occurred in the Southern Alps, in what is now northern Italy (Morelli et al. Reference Morelli, Bargossi, Mair, Marocchi and Moretti2007, Reference Morelli, Marocchi, Moretti, Bargossi, Gasparotto, De Waele, Klötzli and Mair2012; Marocchi et al. Reference Marocchi, Morelli, Mair, Klötzli and Bargossi2008; Willcock et al. Reference Willcock, Cas, Giordano and Morelli2013, Reference Willcock, Bargossi, Weinberg, Gasparotto, Cas, Giordano and Marocchi2015). Reconstruction by Dmitry Bogdanov, sky and digital manipulation by Roman Yevseyev (a photo editor and paleo-artist from Moscow, Russia), sent to us by Amin Khaleghparast, who suggested the colorful sky; all these individuals gave us permission to publish it.
Conodont δ¹⁸O records show the onset of global warming in the early to mid-Kungurian (Chen et al. Reference Chen, Joachimski, Shen, Lambert, Lai, Wang, Chen and Yuan2013). This warming, coupled with glacial melting, transgression, and the oxidation of released methane, likely contributed to widespread marine dysoxia and anoxia (Huang et al. Reference Huang, Sun, Chen, Ogg, Hou, Li, Li and Yong2025). Additionally, lithium isotope data—a proxy for silicate weathering intensity—reveal exceptionally high values during the Kungurian compared with the Asselian–Sakmarian and Capitanian–Changhsingian stages (Kalderon-Asael et al. Reference Kalderon-Asael, Katchinoff, Planavsky, Hood, Dellinger, Bellefroid, Jones, Hofmann, Ossa and Macdonald2021; Wang et al. Reference Wang, Zhang, Liu, von Strandmann, Xiong, Zheng, Garbelli, Zhang, Yuan and Shen2024). The transition from a glacial to a more temperate climate likely intensified runoff and, together with elevated atmospheric CO₂, enhanced incongruent weathering that explains these high δ7Li values.
In summary, the mid-Kungurian experienced a short-term but significant destabilization of the global carbon cycle, driven by volcanic emissions and possibly amplified by methane hydrate release. This provoked a rapid increase in pCO2 and temperature which enhanced the hydrological cycle and increased river runoff. While this event has been documented mainly in the marine record until now, our study reveals a notable terrestrial response as well (Figs. 2–4). This short-lived carbon cycle disruption corresponds to the start of a longer-term climatic shift away from glacial conditions, with subsequent glacial intervals (P3: late Kungurian to latest Roadian; P4: late Wordian to late Capitanian) reflecting more localized glaciation and reduced ice volume (Chen et al. Reference Chen, Joachimski, Shen, Lambert, Lai, Wang, Chen and Yuan2013).
However, the apparent Kungurian mass extinction event might not be directly linked to the abovementioned climatic disturbances. Instead, it may be local and may result from a northward drift of the areas that yielded fossils of the OES grade (see “The Press”). Indeed, the fact that the bulk of the data on the OES grade comes from southwestern United States raises the possibility that the apparent extinction event is a taphonomic artifact. The upper part of the Clear Fork Formation (considered a group in the older literature) and the basal part of the San Angelo Formation, both in Texas, are rich in evaporites and were not conducive to fossilization (Presley and McGillis Reference Presley and McGillis1982); they have yielded no vertebrate fossils. A coeval gap occurs in Oklahoma, as Olson (Reference Olson1967) noticed long ago. Altogether, these barren strata in the late Kungurian and possibly earliest Roadian in Texas amount to a thickness of about 130 m (Laurin and Hook Reference Laurin and Hook2022). Olson (Reference Olson1967: 91) reported that the coeval fossil gap in Oklahoma is represented by “several hundred to nearly 1,000 feet of section between the Hennessey beds and the fossiliferous Chickasha beds.” More fieldwork in late Kungurian continental basins would be useful to untangle possible taphonomic effects from a genuine mass extinction event at that time. Alternatively, similar analyses of biodiversity evolution on other, geographically more widespread clades of amniotes and other continental taxa in the late Kungurian could shed additional light on this matter.
At least two other studies on biodiversity evolution of early amniotes may be relevant to begin to address this issue. The first study assessed rates of morphological evolution in captorhinids (Brocklehurst Reference Brocklehurst2017). While not focused on biological crises, that study showed a clear, spectacular drop in morphological disparity of Captorhinidae in the late Kungurian to early Roadian. The second, more relevant study used tip dating to try to discriminate between two hypotheses (Olson’s gap or Olson’s extinction) to explain the stratigraphic distribution of Permo-Carboniferous amniotes, using three datasets (one on amniotes, one concentrating on captorhinids, and a last one on caseids). That study concluded that the analyses “overwhelmingly” supported the hypothesis of an extinction (Brocklehurst Reference Brocklehurst2017). Some of the data in these two studies could be affected by the same taphonomic factors as our data. However, both captorhinids and caseids extend beyond the early Roadian, and their Guadalupian distribution is documented in other areas, such as Russia (Maddin et al. Reference Maddin, Sidor and Reisz2008; Brocklehurst and Fröbisch Reference Brocklehurst and Fröbisch2017), China (Reisz et al. Reference Reisz, Liu, Li and Müller2011), Niger (de Ricqlès and Taquet Reference de Ricqlès and Taquet1982), and Morocco (Jalil and Dutuit Reference Jalil and Dutuit1996). We also note that Olson’s gap and Olson’s extinction, which are often presented as mutually exclusive explanations of the observed stratigraphic range of amniotes in the Kungurian and Roadian, could coexist, because a genuine extinction event could occur just before a gap in the fossil record. Such a coincidence would challenge paleontologists, but geochronological (especially biostratigraphic) studies appear to rule out the existence of Olson’s gap (Reisz and Laurin Reference Reisz and Laurin2001, Reference Reisz and Laurin2002; Lozovsky Reference Lozovsky2005; Benton Reference Benton2012; Brocklehurst et al. Reference Brocklehurst, Day, Rubidge and Fröbisch2017; Brocklehurst Reference Brocklehurst2020; Laurin and Hook Reference Laurin and Hook2022). Thus, while these studies do not fully address the taphonomic problem that may affect our OES dataset, they do support our tentative conclusion that an extinction event took place in the Kungurian.
The fossil record in the marine realm, which is much more geographically widespread and stratigraphically continuous than that of continental vertebrates, is also relevant here. Biakov (Reference Biakov2012) recognized a late Kungurian mass extinction event in northeast Asia; he documented species-level extinction rates of 78% for bivalves, 86% for brachiopods, and 60% for foraminifera. These data tip the balance in favor of a Kungurian mass extinction event in the OES grade, even though more research is needed to settle this question more conclusively.
Conclusion
Our new analyses confirm the presence of a mild extinction event near the Carboniferous/Permian boundary in synapsids, which we initially detected in varanopids (Laurin and Didier Reference Laurin and Didier2025), but which apparently also affected the OES grade. Climatic and floristic changes around that time may explain this mild crisis. Diversification of the OES grade stopped at that time, perhaps because of global warming and aridification throughout the Cisuralian, as the LPIA waned. We detect a strong mid-Kungurian crisis in the OES grade, which we did not detect (perhaps because of technical limitations and the low sample size) among varanopids. That crisis may result from rapid global warming linked partly to intensive volcanism (Fig. 4). The observed biodiversity pattern may represent a pulse–press–pulse pattern, reminiscent of the simpler press–pulse pattern that may characterize many mass extinction events, according to Arens and West (Reference Arens and West2008). Ancient mass extinction events, which provide a deep-time perspective on biodiversity loss, deserve careful study based on sophisticated analytical methods and datasets incorporating detailed stratigraphic and phylogenetic data. Such studies may improve our understanding of the current, anthropogenic biodiversity crisis, which constitutes a major societal challenge (Wake and Vredenburg Reference Wake and Vredenburg2008; Barnosky et al. Reference Barnosky, Matzke, Tomiya, Wogan, Swartz, Quental, Marshall, McGuire, Lindsey and Maguire2011; Ceballos et al. Reference Ceballos, Ehrlich, Barnosky, García, Pringle and Palmer2015, Reference Ceballos, Ehrlich and Dirzo2017).
Acknowledgments
The authors thank D. Bogdanov (Chelyabinsk, Russia), who drew Figures 1 and 3; R. Yevseyev (Moscow, Russia), who reworked part of Figure 4; and A. Khaleghparast (Tehran, Iran) for establishing the contact between one of us (M.L.) and D. Bogdanov All these individuals gave us permission to publish these images. We thank two anonymous reviewers for their suggestions to improve the article. This work was supported by the recurring grants to the CR2P from the French ministry of research, the CNRS and Sorbonne Université (M.L.). This is the Paleobiology Database (PBDB) contribution number 550.
Competing Interests
The authors declare no conflicts of interest.
Data Availability Statement
The data used in this paper are available at https://doi.org/10.5061/dryad.qfttdz0x8. The code (along with the data) is available at https://github.com/gilles-didier/TimePiecewiseFBD2.
Open access
