A dataset was formulated containing occurrences within each formation of amniotes, diadectomorphs, microsaurs and lysorophians, occurring between their origin in the Carboniferous and the end of the Cisuralian. The data was derived primarily from the Paleobiology Database (PBDB;²) and supplemented from the published literature and observations from museum collections. The PBDB data was downloaded in July 2020, and so taxa published after this were not added to the dataset.

The interval of time under study was divided into time bins, using the Formation Binning approach (Dean et al., 2020). This approach uses the top and bottom ages of formations to assess the most suitable position for the interval boundaries, rather than using bins from the international commission of stratigraphy that are largely based on marine faunal turnover. This method was implemented in R v3.6.1 (R Core Team, 2017) using the functions provided by Dean et al. (2020).

Due to recent studies suggesting that global estimates of diversity are more representative of spatial patterns in sampling than genuine patterns of species richness (Close et al., 2020a, b), this study will focus on more local (alpha) diversity patterns. Within each time bin, the formations were grouped into bioregions. These are areas of endemicity, defined by the taxa present in each, and so their boundaries better represent barriers to dispersal instead of arbitrarily defined boundaries based on geography, geology or sampling e.g., continents, basins, localities (Brocklehurst and Fröbisch, 2018b). The bioregions were defined with the approach of Brocklehurst and Fröbisch (2018b), where formations within each time bin are grouped using two cluster analyses, one based on taxonomic distances, and one based on geographic distances. Clusters shared between the two analyses represent bioregions. This approach allows diversity to be assessed at different spatial scales, defined by the geographic cluster node height, since patterns of diversity will vary at different scales (Palmer and White, 1994; Rosenweig, 1995; Whittaker et al., 2001; Field et al., 2009; Keil et al., 2012; Brocklehurst and Fröbisch, 2018b; Brocklehurst et al., 2018b).

In this study, taxonomic distances between formation were calculated using the modified Forbes metric (Alroy, 2015), computed in R using the script provided in Brocklehurst et al. (2018a). Geographic distances were based on palaeocoordinates derived from the PBDB and calculated using functions in the R package letsR (Vilela and Villalobos, 2015). Diversity within each time bin was analyzed at two spatial scales: a local scale (formations within 100 km of each other united into their bioregions) and a regional scale (formations within 1,000 km of each other united into their bioregions). In the former, the boundaries between bioregions will be more reflective of local environmental/habitat variability, while in the latter geophysical boundaries will be of greater importance.

Within each bioregion, diversity was assessed using shareholder quorum subsampling (SQS) (Alroy, 2010; Chao and Jost, 2012). This approach accounts for sampling heterogeneity by subsampling occurrences to fixed levels of coverage, rather than sample size as in rarefaction approaches. The analysis was carried out in R using functions from the package iNEXT (Hsieh et al., 2016). A coverage quorum of 0.9 was applied, following research suggesting diversity estimates below this coverage become imprecise (Brocklehurst et al., 2018a; Close et al., 2018).

As far as possible, taxa were defined to the species level, but attempts were made to include occurrences that could not be assigned to that taxonomic resolution. If an occurrence was assigned to a higher taxon, and no occurrences of species within that higher taxon are known from that bioregion, then the higher taxon was considered a unique species within that bioregion.

The above procedures were carried out on eight different taxonomic schemes (Table 1), representing different combinations of taxonomic assignments of the disputed clades (varanopids, diadectomorphs, recumbirostrans and lysorophians. Note that while the dataset contains a broader sampling of microsaurs, only recumbirostrans and lysorophians are tested within reptiles; others are included in the dataset to aid future analysis in case of further revisions to microsaurs (see “Discussion”). Within each bioregion, all taxa were assigned to either Reptilia or Synapsida based on the taxonomic scheme in use, and the species richness of each within the bioregion was calculated.

Richards Spur is an Artinskian aged (Woodhead et al., 2010) locality in Oklahoma representing the most diverse tetrapod-bearing locality currently known from the early Permian, with more than 30 taxa currently described (MacDougall et al., 2016). The site represents a cave deposit, with clay and mudstone filling fissures in limestone (Sullivan et al., 2000). Being a cave deposit, unlike many of the lower Permian formations in North America, the mode of preservation is more favorable toward smaller taxa (Muñoz-Durán and Van Valkenburgh, 2006; Jass and George, 2010; Brown et al., 2019). Richards Spur was therefore analyzed separately, to provide a comparison of the diversity of reptiles and synapsids at small sizes.

Diversity of reptiles and synapsids within Richards Spur was again calculated using SQS. However, as a single locality is being analyzed, specimens were subsampled rather than occurrences. Abundance data was drawn from the PBDB and the published literature (Supplementary Data 3). Diversity was calculated at a range of coverage quora, at intervals of 0.05 between 0.4 and 0.95, to examine the impact of varying sampling intensity. Also, due to concerns surrounding the performance of SQS when abundance distributions are highly uneven (Close et al., 2018; Alroy, 2020; see “Discussion”), a second method more robust to these issues was used to test the results: Squares (Alroy, 2018).

Richards Spur contains numerous varanopid species, but only four recumbirostrans and some fragmentary unnamed diadectid specimens (Reisz and Sutherland, 2001). It was therefore decided not to be worth calculating diversity with all eight taxonomic schemes, but to test only a “best case for synapsids” (the varanopids and diadectid specimens included within synapsids and recumbirostrans not included), and a “best case for reptiles (both varanopids and recumbirostrans included within reptiles).

South Grandfield is another well-sampled Cisuralian locality (Kungurian age) containing a diverse array of tetrapods, including good preservation of several small taxa (Daly, 1973). Although the small number of species and large number of singletons (at least when recumbirostrans are not included) makes this locality unsuitable for examination via the species accumulation curves using SQS, squares was used to examine the diversity of reptiles and synapsids within the two taxonomic schemes, again using abundance data drawn from the PBDB and literature (Supplementary Data 4).

While altering the taxonomic scheme employed varies the relative diversities of synapsids and reptiles, the overall trends of species richness remain consistent (Figures 2, 3). The earliest interval where sampling is sufficient to assess diversity of both reptiles and synapsids is the Moscovian (Pennsylvanian). Median synapsed diversity is stable through the latter stages of the Carboniferous, before rising gradually to a late Artinskian peak, visible when examining diversity at local scales (at regional scales the late Artinskian does not have the spatial extent of sampling to examine diversity in this time bin). Median reptile diversity also remains relatively stable in the Pennsylvanian, before rising to an early Artinskian peak, visible at both local and regional scales. During the late Kungurian, synapsed diversity at both local and regional scales falls sharply. At local scales reptile diversity remains stable, but a decline is observed in their regional diversity.

When using taxonomic schemes 1 (“best case for synapsids”) and 2 (pre-2017 taxonomy), synapsids are found to be more diverse than reptiles throughout much of Pennsylvanian and Cisuralian, with higher median diversity at both local and regional scales (Figures 2A,B, 3A,B). There are only two time bins where median reptile diversity exceeds that of synapsids when using taxonomic schemes 1 and 2: the early Artinskian and, at local but not regional scales, the latest Kungurian. As more of the recent taxonomic revisions are incorporated, however, median reptile diversity is pushed closer to that of synapsids (Figures 2C–G, 3C–G). In taxonomic schemes 7 (all recent taxonomic changes applied) and 8 (“best case for reptiles”), median diversities of reptiles and synapsids are consistent throughout the Pennsylvanian and Cisuralian (Figures 2G,H, 3G,H). In fact, at local scales during the latest Kungurian, reptile diversity substantially exceeds that of synapsids when using these taxonomic schemes (Figures 2G,H).

When using the taxonomic scheme most favorable to synapsids, at low levels of coverage synapsed diversity at Richards Spur substantially exceeds that of reptiles (Figure 4A). However, at higher levels of coverage (above 0.75), reptile diversity is found to be higher (Figure 4A). When using the taxonomic scheme most favorable to reptiles, reptile diversity in Richards Spur is found to exceed that of synapsids at all levels of coverage (Figure 4B). In both taxonomic schemes, analysis using squares produced higher diversity of reptiles than synapsids at Richards Spur. When using the best case for synapsids, synapsed diversity calculated with squares is 8, reptile diversity 17.84. When using the best case for reptiles, synapsed diversity calculated with squares is 3, reptile diversity 26.31. At South Grandfield, analysis using squares also produced higher diversity of reptiles than synapsids regardless of the taxonomic scheme. When using the worst case for reptiles, reptile diversity calculated using squares is 6.39; when using the best case for reptiles it is 10.61. Synapsid diversity using both schemes is 2.71.

The recent revisions of amniote phylogeny will doubtless be subject to further debate and modification, and there may well be other clades whose position will come under scrutiny. It is therefore not worth discussing at this point in time which of the global curves of diversity through time is most likely to represent the true signal. It is more useful to examine general patterns across the different datasets, which will indicate what we are still able to say with confidence, and what work is most important in the future.

First, the relative diversity patterns during the latest Carboniferous and earliest Cisuralian are most strongly affected by the assignment of the varanopids to either synapsids or reptiles. The Asselian stage documents a substantial increase in varanopid diversity in North America, particularly west of the Hueco seaway. This shallow inland seaway separated the formations in Texas and Oklahoma from those in New Mexico, Colorado, Arizona, and Utah (Lucas et al., 2011), and there are noticeable differences in the faunas on either side (Brocklehurst and Fröbisch, 2018b; Brocklehurst et al., 2018b). The Cutler Group of New Mexico has a considerable diversity of varanopids, including Ruthiromia (Eberth and Brinkman, 1983), two species of Aerosaurus (Romer, 1937; Langston and Reisz, 1981), and potentially Nitosaurus [(Romer, 1937), although this taxon’s affinity and validity is uncertain (Reisz, 1986)]. There is also an indeterminate varanopid in the Placerville localities of Colorado (Lewis and Vaughn, 1965; Brocklehurst et al., 2016). Their assignment to reptiles pushes median reptile diversity during the Gzhelian and Asselian close to that of synapsids and alters the timing and rapidity of the reptile radiation (although note Maddin et al. (2020) discuss how ontogenetic changes in varanopids may lead to taxonomic oversplitting). When varanopids are assigned to Synapsida, reptile diversity remains low in the initial stages of the Cisurialian, and the rise to their Artinskian peak is rapid (Figures 2A,B,D,F, 3A,B,D,F). When they are assigned to Reptilia, reptile diversity rises gradually through the early stages if the Cisuralian (Figures 2C,E,G,H, 3C,E,G,H).

The global diversity estimates also provide an interesting perspective on the latest Cisuralian amniote diversity. The late Kungurian and early Roadian were intervals of substantial change in tetrapod faunas, with a decline in the pelycosaurs and amphibians that had thrived during the early Permian, and a radiation of therapsid synapsids and procolophonian parareptiles (Kemp, 2006; Benson and Upchurch, 2013; Brocklehurst et al., 2013, 2017). Substantial ecological shifts also occurred, with the establishment of herbivore-dominated ecosystems and increasingly complex food webs with more trophic levels (Olson, 1966), as well as increases in faunal provinciality (Brocklehurst and Fröbisch, 2018b). This transition has been suggested to have been accompanied by a mass extinction, dubbed Olson’s Extinction (Sahney and Benton, 2008). Substantial debate does surround this event, it having been argued to be an artifact of spatial patterns of sampling (Benson and Upchurch, 2013), incorrect dating of formations (Lucas, 2017, 2018), and inappropriate time-binning methods (Lucas, 2017). However, thorough study of these issues still supports the extinction as a genuine event (Brocklehurst et al., 2017; Brocklehurst, 2018, 2020).

All diversity curves presented here indicate a substantial decrease in synapsed diversity during the latest Kungurian (Figures 2, 3). However, the pattern in reptiles is more complex and the signal varies depending on the spatial scale examined. At the local scale, reptile diversity remains stable between the middle and late Kungurian (Figure 2), whereas at regional scales there is a decrease (Figure 3). This pattern would indicate a more cosmopolitan reptile fauna developing during the late Kungurian relative to that of synapsids. When amniote-bearing localities are combined into larger bioregions, the number of synapsed species in each region increases, indicating different faunas are being combined. This is not the case in reptiles, potentially indicating a broader range of environmental tolerances; different local-scale bioregions have similar taxa present. A similar conclusion has been drawn from event-based biogeography analysis (Brocklehurst et al., 2018b).

While ichnotaxa are not considered in this analysis, the tetrapod footprint record provides support for decreasing diversity of synapsids during the Kungurian, and their being superseded in abundance by the reptiles. The Kungurian stage is roughly contemporary with the Erpetopus ichno-biochron (Voigt and Lucas, 2018). This biochron is characterized by a diverse array of trackways usually assigned to reptilian taxa, in particular Erpetopus, Varanopus and Hyloidichnus, while the Dimetropus tracks abundant in lower stratigraphic levels and usually associated with pelycosaurian grade-synapsids are considerably rarer (Haubold and Lucas, 2003; Voigt and Lucas, 2018; Marchetti et al., 2019). The footprint record provides a less spatially restricted sample that the body-fossil record, providing data on this reptilian diversification in North America, Europe and North Africa, in which reptile tracks make up between 50 and 100% of ichnoassociations (Marchetti et al., 2019).

Apparent reptile diversity during Olson’s Extinction is particularly influenced by the exclusion or inclusion of recumbirostrans and lysorophians. At local scales it is this which determines whether late Kungurian reptile diversity is similar to that of synapsids or substantially higher. At regional scales it determines whether reptile diversity is substantially lower than that of synapsids or is more similar. These two lineages are among the most diverse and abundant taxa in the Clear Fork group and to a lesser extent the Hennessey Formation, the set of formations that most clearly documents Olson’s extinction (see Supplementary Data 1). In fact, all lineages previously assigned to microsaurs show an increase in abundance and richness at this time, underscoring the necessity for further taxonomic revision of this group. Phylogenetic analysis and anatomical revisions of taxa are increasingly indicating that “Microsauria” is a polyphyletic assemblage of taxa, with many showing similarities in morphology with amniotes (e.g., Szostakiwskyj et al., 2015; Pardo et al., 2017; Mann et al., 2019, 2020; Marjanović and Laurin, 2019). A recent example is Asaphestera from Joggins (Bashkirian age) previously assigned to tuditanid microsaurs, but recently found to be a chimera of which part was reassigned to Synapsida, making it the earliest known member of that lineage (Mann et al., 2020).

While there will doubtless be further debate regarding the taxonomy of amniotes, and which lineages should be assigned to reptiles or synapsids, what is apparent from these analyses is that many of the recent taxonomic revisions are pushing the diversity of reptiles closer to that of synapsids during the Pennsylvanian and Cisuralian. Even adding diadectomorphs to synapsed diversity does not substantially increase the median diversity in most time bins. When using a taxonomy incorporating all recent revisions, including the addition of diadectomorphs to synapsids, there is little difference in median diversity of reptiles and synapsids throughout the Pennsylvanian and Cisuralian (Figures 2G, 3G). If these recent revisions are further supported by future research, the apparent dominance of pelycosaurs during the Cisuralian may have to be reconsidered.

The majority of the data on Cisuralian tetrapods comes from North America, in particular the Texas Red Beds (Romer and Price, 1940; Kemp, 2006; Benson and Upchurch, 2013; Brocklehurst et al., 2013). These formations are heavily biased in their mode of preservation toward larger taxa (Behrensmeyer, 1988), a bias that will favor the preservation of synapsids over reptiles. While synapsids reached large sizes early in their evolution, with multiple lineages independently evolving sizes of greater than 40 kg during the Carboniferous (Reisz and Fröbisch, 2014; Brocklehurst and Brink, 2017; Brocklehurst and Fröbisch, 2018a), reptiles did not reach such sizes until the latest Cisuralian with the evolution of Moradisaurinae (Brocklehurst, 2016).

The mode of preservation in the cave deposits at Richards Spur provides a unique taphonomic window into the small vertebrates present during the Cisuralian (MacDougall et al., 2016). Richards Spur contains a considerably more diverse reptile fauna than elsewhere in the Cisuralian and is partially responsible for the early Artinskian peak in reptile diversity observed in the global analysis. There is an unparalleled diversity and abundance of captorhinids (deBraga et al., 2019), with the genus Captorhinus being represented by hundreds of specimens (Richards, 2016). There is also a substantial parareptile diversity (MacDougall et al., 2016), including lineages not found elsewhere until the middle Permian (MacDougall and Reisz, 2014).

Richards Spur’s unusual combination of taxa and different relative abundances of species compared to contemporary localities has been explained by it being an upland locality, while most early Permian formations represent lowland fluival/lacustrine environments (Sullivan et al., 2000). However, there are two key issues that make Richards Spur an extremely important data point. First, lagerstätten (areas of exceptional preservation) frequently preserve assemblages of taxa atypical of non-lagerstätten formations (Walker et al., 2020). This is likely due to the unique modes of preservation in such localities, that often selectively preserve small delicate taxa that are usually easily destroyed by taphonomic processes, but in the environments that characterize many lagerstätten are easier to bury rapidly than larger taxa (Brocklehurst et al., 2012; Dean et al., 2016; Brown et al., 2019). Selective preservation of smaller taxa is characteristic of cave deposits (Muñoz-Durán and Van Valkenburgh, 2006; Jass and George, 2010; Brown et al., 2019) like Richards Spur.

The second point that must be remembered is that, in modern ecosystems, most species are small. The greater species richness of small animals is observed across a wide range of taxa, both vertebrate and invertebrate (e.g., May, 1978; Brown and Nicoletto, 1991; McClain, 2004), and is underpinned by theoretical considerations (Hutchinson and MacArthur, 1959), so there is no reason to expect that different patterns were prevalent in Paleozoic ecosystems. Therefore, the differences observed between the faunas preserved in lagerstätten and those in other formations should not be interpreted as the lagerstätten being aberrant. Instead, formations that selectively preserve smaller taxa should be considered more representative of the true diversity patterns in an ecosystem.

Thus, Richards Spur, instead of representing an aberrant fauna, could represent a crucial insight into the organisms that would be expected to make up the majority of the early Permian fauna, but are less well-sampled in the majority of available localities due to the size-based preservation biases. South Grandfield (Hennessey Formation, Kungurian in age), which represents a lowland fluvial system (Daly, 1973) is another locality that contains unusually good preservation quality of small taxa, and so is also useful in this context.

Analyses of reptile and synapsed diversity within Richards Spur support reptiles as being more diverse, whether using a taxonomic scheme favorable to reptiles or synapsids. It is true that, when using the scheme favorable to synapsids, synapsed diversity is found to be higher when subsampling to lower coverage levels (Figure 4A), but this is likely an artifact of a bias affecting the SQS method. SQS is heavily influenced by the shape of the relative abundance distribution (Close et al., 2018; Alroy, 2020): when it is uneven, with a small number of hyper-abundant taxa easily sampled, and most others being rare, it is easier to reach low coverage quora after sampling a small number of taxa. The extreme abundance of Captorhinus aguti at Richards Spur is likely confounding the estimates of reptile diversity in this way. When subsampling to higher coverage levels that produce more robust results (Brocklehurst et al., 2018a; Close et al., 2018) and when using the Squares method that is more resilient to this bias (Alroy, 2018, 2020), reptiles are found to be more diverse at Richards Spur no matter which taxonomic scheme is employed (Figure 4).

South Grandfield, which represents a lowland fluvial system but also shows unusual preservation quality of small taxa, preserves a similar fauna to that of Richards Spur: a great abundance of captorhinid reptiles, with large synapsids forming a relatively minor component of the fauna (Daly, 1973). Reptile diversity within South Grandfield is here found to be higher than that of synapsids whether or not the numerous recumbirostran species known from the locality are included within Reptilia.

Pennsylvanian Largestätte, also with preservation modes favoring of small taxa, could also be argued to show a greater emphasis on reptile diversity, although most of these contain better sampling of non-amniote tetrapods and did not contain the coverage of both synapsids and reptiles necessary to perform the species-accumulation analyses applied to Richards Spur. The Linton and Nýřany tetrapod assemblages of Ohio and Czechia, respectively (both of Moscovian age) contain only a single synapsed (Archaeothyris) but multiple reptiles and numerous microsaurs (Carroll and Baird, 1972; Reisz, 1975; Hook and Baird, 1986, 1988). The Garnett Quarry (of Kasimovian age) contains a greater observed species richness of synapsids, but by far the most abundant taxon is the diapsid Petrolacosaurus (Peabody, 1952; Wehrbein, 2017).

While this study focusses on the Pennsylvanian and Cisuralian, the interval previously thought dominated by pelycosaurian-grade synapsids, it is worth discussing the middle and late Permian (Guadalupian and Changhsingian) in the context of this lagerstätten effect. Reptile and synapsed diversity estimates have been compared during these intervals within the most productive tetrapod-bearing basin of that time: the Karoo of South Africa. Synapsids (then mostly represented by the therapsids rather than pelycosaurs) were found to be consistently more diverse than the reptiles (Fröbisch, 2013). However, this basin, and many other contemporary tetrapod-bearing basins, again show preservation heavily selective toward larger, more robust specimens (Kammerer, 2016), perhaps again biasing diversity estimates in favor of synapsids. A formation which provides a better record of small-bodied amniotes is Mezen in Russia, of early Guadalupian age (Golubev, 2015). Mezen’s depositional environment is thought to have been repeated, low-current shallow flooding events, which overtook small amniotes but allowed larger taxa to escape (Efremov, 1940). This fauna, where small taxa have been preferentially preserved, is again overwhelmingly dominated by reptiles. More than 70% of specimens from Mezen are parareptiles (see Supplementary Text and Data in Brocklehurst et al., 2017), while therapsid synapsids are rare (Lozovsky, 2005; Brocklehurst et al., 2017). There are also abundant varanopids at Mezen (Reisz and Berman, 2001; Anderson and Reisz, 2004), previously considered pelycosaurian-grade synapsids, but now potentially assigned to reptiles (Ford and Benson, 2019, 2020). While it must be noted that Mezen is geographically separated from other contemporary Russian localities, potentially providing an alternative explanation for the unusual fauna (Ivakhnenko, 2001), it is also possible that Mezen has allowed the preservation of the abundant small-bodied taxa that were prevalent elsewhere but less easily preserved (see Supplementary Text in Brocklehurst et al., 2017). In Mezen further corroboration is found of the signal found in Richards Spur: the selective preservation of smaller taxa, which we expect to make up the bulk of diversity, produces a reptile-dominated fauna.