Estimating spatial variation in origination and extinction in deep time: a case study using the Permian–Triassic marine invertebrate fossil record

Origination and extinction are two of the most fundamental processes in evolution, structuring taxonomic diversity across space and time (Allen and Gillooly 2006; Jablonski 2008; Schluter and Pennell 2017). Understanding spatial variation in origination and extinction can help to unravel the mechanisms underlying macroevolutionary patterns, such as latitudinal diversity gradients (Mittelbach et al. 2007; Mannion et al. 2014; Powell et al. 2015; Schluter and Pennell 2017; Saupe et al. 2019; Meseguer and Condamine 2020), and provide deeper insight into biotic crises and radiations (Jablonski 2008; Saupe et al. 2014; Kocsis et al. 2018; Reddin et al. 2019; Flannery-Sutherland et al. 2022).

Although methods have been developed for estimating global origination and extinction rates from the fossil record, no framework exists for applying these methods to restricted spatial regions. Raw evolutionary rates have been calculated previously within latitude bands (e.g., Powell et al. 2015; Song et al. 2020), but this approach does not take into account variation in the sampling completeness of the fossil record (Benson and Upchurch 2013; Vilhena and Smith 2013; Mannion et al. 2014; Fraser 2017; Close et al. 2020; Benson et al. 2021; Jones et al. 2021; Shaw et al. 2021), and the relationship between these estimates and the true rates is uncertain. Other studies have attempted to circumvent the issue by categorizing taxa based on the latitude at which they are most abundant (e.g., Clapham et al. 2009; Reddin et al. 2019), but such an approach can involve compromises when taxon ranges cross latitudinal bins. Regional origination and extinction have also been explored in some studies (e.g., Dunhill et al. 2018; Flannery-Sutherland et al. 2022), but this approach can conflate immigration with true speciation and extirpation with true extinction.

Here, we examine how well three commonly used origination and extinction metrics perform when applied regionally, using simulated data to investigate whether spatial differences in origination and extinction can be estimated reliably from fossil data. These three metrics were then applied to empirical datasets of Permian and Triassic marine invertebrate occurrences, to examine the evidence for contrasting origination and extinction rates among clades and between high and low latitudes.

Extinction and origination rates varied substantially during the Permian and Triassic (~300–200 Ma), an interval marked by several major biotic crises coincident with the eruption of large igneous provinces (LIPs). Heightened extinction at these times was likely driven by associated global warming and ocean acidification and anoxia (Wignall et al. 2010; Payne and Clapham 2012; Sun et al. 2012; Wignall 2015; Penn et al. 2018; Dal Corso et al. 2022). Extreme fluctuations in evolutionary rates, combined with the relative spatial completeness of the fossil record at this time (Jones et al. 2021), makes the Permian–Triassic an ideal interval for testing how origination and extinction may vary spatially (e.g., Flannery-Sutherland et al. 2022).

The Capitanian biotic crisis (CBC, also known as the end-Guadalupian extinction; ~260 Ma) occurred near the end of the middle Permian and coincided with the eruption of the Emeishan LIP (Clapham et al. 2009; Wignall et al. 2009; Bond et al. 2010; Wignall 2015; Rampino and Shen 2019). A coincident fall in global sea level may have further affected marine ecosystems and their preservation potential (Shen and Shi 2004; Clapham et al. 2009; Clapham 2015). Estimates of the taxonomic severity of the CBC have varied (McGhee et al. 2013; Rampino and Shen 2019), and diversity changes have been attributed previously to back-smearing of the later Permian–Triassic extinction (Foote 2007) or reduced origination in the Capitanian and Wuchiapingian (Clapham et al. 2009). Extinctions associated with the CBC appear to have been highly ecologically selective, resulting in the loss of sponge–microbial reefs, high extinction levels in foraminifera and calcareous algae, and rapid turnover in ammonoids (Wignall et al. 2009; Bond et al. 2010; McGhee et al. 2013; Clapham 2015; Rampino and Shen 2019).

The Permian–Triassic mass extinction (PTME; ~252 Ma), the most severe in Earth history, saw an estimated loss of 83% of marine genera (McGhee et al. 2013). Losses were particularly profound for marine reefs (Wignall 2015; Martindale et al. 2019). The PTME was caused by the eruption of the Siberian Traps LIP, which resulted in extreme warming and anoxia and acidification in the oceans (Kiehl and Shields 2005; Payne and Clapham 2012; Sun et al. 2012; Wignall 2015; Penn et al. 2018; Martindale et al. 2019; Dal Corso et al. 2022). Considering the relative stability of the environmental niches of marine invertebrates over geological timescales (Saupe et al. 2014), previous discussion of kill mechanisms associated with the PTME offers two hypotheses concerning the spatial distribution of extinctions:

1. Increasing temperatures in lower latitudes rendered these regions inhospitable for most animals, driving high extinction rates at the equator and poleward migration (Sun et al. 2012; Bernardi et al. 2018; Allen et al. 2020; Song et al. 2020; Flannery-Sutherland et al. 2022).

2. Increasing temperatures (and anoxia) in the polar regions left cool-adapted organisms with no temperature-suitable habitat, leading to high extinction rates at high latitudes (Penn et al. 2018).

These hypotheses are not mutually exclusive and can be tested independently. Evidence supporting the first hypothesis has been reported previously for marine environments during other intervals of global warming in Earth history, including the Triassic–Jurassic mass extinction (Kiessling and Aberhan 2007; Dunhill et al. 2018; Reddin et al. 2019). We use the origination and extinction metrics tested in our simulation framework to evaluate support for these two hypotheses.

We used two simulation frameworks to test the efficacy of several rate estimation methods when applied regionally. The first assumed no extinction selectivity, while the second investigated the effect of extinction selectivity with respect to abundance. We then applied the rate estimation metrics to Permian–Triassic marine invertebrate occurrence data.

All analyses were conducted in R (R Core Team 2021) using the tidyverse (Wickham et al. 2019), pspearman (Savicky 2014), and lsr (Navarro 2015) packages. All R code demonstrating how to apply these metrics in the manner outlined here is provided within the Supplementary Data.

We constructed a simulation to produce fossil occurrence data (Fig. 1) using the protocol of Barido-Sottani et al. (2020). Our approach allowed “true” origination and extinction to be measured as a benchmark for method comparison, which cannot be achieved using empirical fossil data. Simulated datasets were designed to fit the standard format of fossil occurrences available in the Paleobiology Database (PBDB), that is, a list of occurrences, each representing the presence of a particular species within a “collection” or locality, agglomerated across a specified geographic area. Simulation input parameters were initially based on values calculated from Permian–Triassic marine invertebrate occurrences (see “Applying Metrics to Permian–Triassic Fossil Data”) and subsequently amended to increase the range of values included within the simulation outputs (Supplementary Fig. S1).

Initial starting conditions (t₀) consisted of a “world” split into six spatial bins, each containing 1000 species occurrences (here considered akin to 30° latitudinal bands, but they could equally represent any six spatial subdivisions, such as bioregions, marine basins, or continents; Fig. 1). The size of the global species pool was drawn at random, containing between 100 and 800 species, and each spatial bin was generated independently by drawing species’ identities at random from the global species pool.

The simulated occurrences were then subjected to three iterations of “origination” and “extinction” to produce a four-slice time series (t₀, t₁, t₂, t₃). To produce each subsequent time slice, a random proportion of occurrences from a given spatial bin, between 0% and 20%, were selected to survive, with the others going extinct.

Origination was simulated by adding a random number of occurrences to the spatial bin, between 0 and 300, with their identities selected from a pool consisting of species present in any of the six spatial bins in the previous time slice, plus a random number of new species, selected between 0 and 400. These processes were carried out independently for each spatial bin, operating at the level of a local population. This procedure allowed migration of species between spatial bins across the different time slices: for example, a particular species could suffer local extinction(s) in t₁ but be selected for local origination(s) in t₂, a process that would not be counted as “true” extinction or origination events, regardless of the spatial bins within which these phenomena took place (Fig. 1B). Our simulation design permitted migration between any pair of spatial bins, irrespective of their identities (or apparent positioning relative to one another), in adjacent time slices.

Once each bin had been populated with occurrences, we replicated the sampling filters known to exist in the fossil record by subsampling each spatial bin in every time slice once. The number of occurrences to subsample was drawn from a uniform distribution between 0 and the total number of occurrences, with all occurrences having an equal probability of being selected.

The described simulation protocol was repeated across 10,000 iterations to incorporate variation in origination, extinction, and sampling completeness (Supplementary Fig. S1). We also ran simulations of 5000 iterations using different combinations of low (0 to ~0.2), medium (~0.1–0.3), and high (~0.3–0.5) origination and extinction levels, to investigate the impact of turnover rate on the accuracy of estimates (Supplementary Figs. S6–S11).

Our simulation was purposefully simplistic and therefore may not replicate the complexity of empirical fossil datasets. For example, the six latitudal bins were initially allocated equal numbers of occurrences, which does not reflect the differences in habitat area available across latitudinal bands. Sampling bias was modeled as random, which is an oversimplification of the interaction of the various facets of sampling bias, many of which are at least somewhat systematic (e.g., Darroch et al. 2020; Benson et al. 2021). The ability of species to migrate between any two latitudinal bands is also unrealistic. However, the ways in which species’ ranges shift on geological timescales is poorly understood, and the approach used here avoids applying potentially false assumptions to this facet of the simulation. In this simulation framework, we assumed that common and rare species were equally likely to experience local extinction, but we investigated the impact of this assumption on metric accuracy with further simulations (see “Simulation with Abundance Selectivity”).

Three metrics designed to estimate global origination and extinction were applied to the simulated data: (1) raw, (2) boundary-crosser (BC) (Alroy 1996; Foote 1999), and (3) three-timer (3T) (Alroy 2008). Origination and extinction were calculated as proportions, indicating the fraction of species within the focal time slice, using these metrics (Fig. 1B).

• Raw values were calculated as the proportion of taxa in a given spatial bin that were not represented in the previous (origination) or subsequent (extinction) time interval. Raw values calculated using the complete (unsampled) datasets were deemed the “true” values to which all other estimates were compared.

• The BC method is based on cohort analysis (Raup 1978) and evaluates the proportion of taxa that do not cross the “bottom” (originations) or “top” (extinctions) of a given time interval. While the method was developed by Foote (1999), the proportions calculated here are based on the equations of Alroy (1996), which exclude singletons (species known only from a single time slice).

• The 3T method, described by Alroy (2008), is also based on cohort analysis, but incorporates an estimate of sampling completeness. The method uses the proportion of “part-timers,” taxa that are present in the first and third slices of a time series but not the second, relative to “three-timers,” which are taxa present in all three time slices, to adjust the calculated proportions of origination and extinction.

The 3T approach has been developed further, including the “gap-filler” (Alroy 2014) and “second-for-third” (Alroy 2015) modifications. These metrics are more complex to calculate, which gives them improved performance when interrogating changes in global biodiversity through time compared with the above three metrics (Alroy 2014, 2015). However, their parameterization makes them difficult to apply to spatially restricted datasets, and thus they are not tested here.

These origination and extinction metrics are usually applied to global datasets, but their application here was altered to estimate origination and extinction within individual spatial bins. The implementation of raw, BC, and 3T methods to assess origination and extinction for a single spatial bin is demonstrated in Figure 1B. Only the focal time slice (t₁ for extinction, t₂ for origination) was filtered to the relevant spatial bin, and all comparisons were made with global datasets in the other time slices. This meant that species that migrated between spatial bins, experiencing local but not global origination and extinction, were not included in estimates. In this approach, sampling bias was particularly important in determining whether species were correctly identified as present within the spatiotemporal bins they occupied, and our subsequent discussion of sampling bias refers to this as the definition of “completeness.” We also calculated metric estimates for each simulation iteration globally, to provide a benchmark for their performance as they are usually applied.

The performance of the three metrics was assessed using two different approaches. Differencing, that is, subtraction, was used to evaluate the absolute accuracy of estimates for individual spatial bins. Student's t-tests and Cohen's D were used to evaluate whether the distributions of error for each metric had mean values different from 0 and from one another's mean values. The ability to recreate the gradient of origination and extinction proportions across the six spatial bins in a single iteration was also assessed by: (1) evaluating whether the maximum and minimum values were attributed to the same spatial bins; (2) calculating Spearman's rank correlation coefficients between the estimated and “true” values; and (3) quantifying the number of spatial bins with overestimated values, rather than underestimated or identical values, to indicate the consistency of the direction (or sign) of numerical difference within individual iterations or “worlds.”

The first simulation framework assumed no relationship between the number of occurrences of a species and its probability of extinction at a local level. However, this assumption may not be valid in real biological systems (e.g., Brown and Kodric-Brown 1977; Fagan et al. 2002; Kiessling and Aberhan 2007; Hooftman et al. 2016; Casey et al. 2021). We therefore ran additional simulations in which this assumption was relaxed and investigated the possibility of an interaction between sampling completeness and abundance selectivity when estimating evolutionary rates.

To generate the initial dataset, occurrences for each of 400 taxa were drawn from those of Wordian brachiopod genera in the PBDB, to generate a realistic relative abundance distribution. We then simulated an extinction event in this population, with a “true” taxon extinction proportion of around 50%. Extinction selectivity in each trial was determined using a random value from a Poisson distribution, treating zero values as extinction and nonzero values as survival (Fig. 2). Sequences of 400 evenly spaced numbers were generated, with each sequence centered around a value of λ = 0.69 to result in approximately 50% zero and 50% nonzero values. We varied selectivity by altering the range of the sequence around the center, with broader ranges leading to greater differences in extinction risk for rare and common taxa. For example, a sequence of 400 λ values between 1.29 and 0.09 would correspond to extinction risk of 27.5% for the rarest taxon and 91.4% for the commonest taxon. A sequence of 400 λ values between 0.39 and 0.99 would correspond to extinction risk of 67.7% for the rarest taxon and 37.1% for the commonest taxon. A sequence of 400 λ values with a range of 0 (i.e., all λ values are equal to 0.69) would correspond to an extinction risk of approximately 50% for all taxa. We varied the range of the sequences from −1.2 (1.29 to 0.09) to 1.2 (0.09 to 1.29) at increments of 0.002 units, resulting in 1201 different selectivity scenarios.

After assigning survival or extinction status to each taxon, we calculated “true” selectivity using a logistic regression describing the probability of survival as a function of the number of occurrences. The log-odds ratio describes how the odds of survival change for each one-unit increase in a taxon's occurrence count. A positive log-odds ratio corresponds to increased survival among common taxa, and a negative log-odds ratio corresponds to increased survival among rare taxa.

Five different sampling levels were then applied, drawing 10% (“smallest”), 32.5% (“small”), 55% (“medium”), 77.5% (“large”), and 100% (“perfect”) of occurrences. The simulations used all combinations of sampling level and selectivity ranges, resulting in 6005 trials. To evaluate the impact of extinction selectivity based on abundance on the accuracy of extinction estimates, we subtracted the “true” extinction proportion for each simulation from the estimated raw extinction proportion after subsampling.

We also investigated the performance of raw (range through), BC, and 3T extinction metrics under conditions of both varying extinction selectivity and incomplete sampling in specific time bins. We used the same five sampling levels and extinction selectivity method described above, varying sampling in time bin t₀ (the interval immediately preceding the focal time interval) and sampling in time bin t₂ (immediately following the focal time interval). We then examined the combined effects of incomplete sampling in the focal interval t₁, the preceding interval t₀, and the succeeding interval t₂. To keep the number of permutations manageable, for this analysis we used two sampling levels for t₁, 30% (“small”) and 60% (“moderate”) of occurrences, and three sampling levels for t₀ and t₂, 30% (“small”), 60% (“moderate”), and 100% (“perfect”) of occurrences.

Having investigated the efficacy of the evolutionary rate estimation metrics, we applied these metrics to the Permian and Triassic marine invertebrate fossil record. Occurrences identified to genus or species level from four major invertebrate clades (Ammonoidea, Bivalvia, Brachiopoda, Gastropoda) were downloaded from the PBDB in March 2021. The occurrences date from the Artinskian (middle early Permian) to the Norian (middle Late Triassic), enabling origination and extinction proportions to be estimated for the Roadian (early middle Permian) to the Ladinian (late Middle Triassic). Uncertain generic identifications were excluded, as were any collections from nonmarine environments or those dated less precisely than to a single stage. The cleaned dataset consisted of a total of 90,209 occurrences.

The paleolatitudes occupied by these fossils were determined by rotating their modern-day locality coordinates back to their time of deposition, stage-by-stage, using the PALEOMAP Global Plate Model (Scotese 2016) implemented in GPlates (Müller et al. 2018). Following this, occurrences were allocated to 30° latitudinal bands, approximately representing tropical (0°–30°), temperate (30°–60°), and polar (60°–90°) regions in each hemisphere. Raw, BC, and 3T proportions of origination and extinction were calculated for each latitudinal bin containing more than five genera in every stage, as described in Figure 1B. To increase the completeness of the available fossil record, origination and extinction were evaluated at the genus level, but estimates were also calculated for species (Supplementary Fig. S13).

To test the accuracy of the three rate estimation methods when applied regionally, we calculated the numerical difference between estimated proportions of origination and extinction and their corresponding true values for individual spatial bins (n = 6 spatial bins × 10,000 iterations = 60,000 bins). Randomly sampling the simulated occurrences, to emulate the patchy nature of the fossil record, reduced both the precision and accuracy of origination and extinction estimates for all three metrics (Fig. 3, Supplementary Table S1). Following sampling, raw and BC metrics tended to overestimate both origination and extinction (Table 1). The 3T estimates were more accurate, with mean and median differences from the true values much closer to zero, and narrower interquartile ranges (IQRs) (Fig. 3, Table 1). For all metrics, extinction estimates were slightly more accurate than origination estimates. These trends were also seen when applying the metrics to the global datasets (n = 10,000 iterations), but with slightly smaller overall ranges and IQRs (Fig. 3).

Examining estimates based on sampling completeness of the relevant spatial bin revealed that, for all three metrics, the ranges and IQRs of difference reduced as sampling completeness increased (Supplementary Fig. S2b). However, although 3T estimates approached the true values with increased sampling, the raw and BC estimates did not, instead tending toward overestimation of ~0.1–0.2. Additional analyses examining the effect of low versus high origination and extinction levels (0 to ~0.5) on estimate accuracy indicated that the range of differences increased with higher turnover (Supplementary Fig. S8).

The accuracy of rate estimation metrics was also assessed by considering the gradient in origination and extinction proportions across spatial bins within a single iteration (n = 10,000 iterations). We compared true and estimated proportions of origination and extinction to test (1) whether maximum and minimum values were attributed to the same spatial bins, (2) whether true and estimated values were correlated with one another, and (3) whether estimated values were consistently overestimated (or underestimated) within each iteration. The gradients of BC and 3T estimates for any given iteration were identical, because the 3T method is effectively the BC method with an added sampling correction that can only be applied at a global level.

After sampling, spatial bins with the highest and lowest proportions of origination and extinction were identified correctly in only a third of the iterations, regardless of the metric used (Supplementary Fig. S3). Raw estimates were slightly more successful than BC and 3T estimates. Across metrics, the spatial bin with the highest proportional origination was identified correctly in more iterations than the bin with the lowest proportional origination. For extinction, bins with the highest and lowest proportional extinction were approximately equally likely to be identified correctly. Additional analyses examining low versus high origination and extinction levels showed that higher turnover resulted in a slight reduction in the ability to correctly identify bins with the lowest and highest proportional origination and extinction (Supplementary Fig. S9).

Rank correlation between post-sampling estimates and true origination and extinction values was assessed using Spearman's ρ (Supplementary Fig. S4). In most iterations, the correlation between post-sampling estimates and true origination and extinction values was positive and of intermediate strength (~0.3–0.7), but with corresponding p-values indicating nonsignificant relationships. These results were also seen in the simulations with low, medium, and high origination and extinction levels (Supplementary Fig. S10).

The proportion of spatial bins within each iteration that overestimated proportional origination and extinction (Supplementary Fig. S5) had a similar distribution to that of individual spatial bins (Fig. 2). Following sampling, raw and BC estimates generally overestimated proportional origination and extinction in five, or all six, of the spatial bins. In contrast, 3T estimates tended to overestimate proportional origination and extinction in two to four of the six spatial bins, suggesting a relatively even distribution of differences above and below the true values within a single iteration.

Our simulations indicated that estimated extinction proportions were increasingly biased as the strength of the relationship between abundance and extinction increased and as sampling completeness decreased (Fig. 4). If common taxa were more likely to survive, as is often assumed to be the case, smaller subsamples resulted in a lower extinction estimate. This bias exists because a smaller subsample is less likely to capture rare taxa, and common taxa will be relatively overrepresented, leading to lower estimated extinction levels in the subsample. Conversely, if rare taxa are more likely to survive, smaller subsamples will tend to overestimate extinction.

Incomplete sampling of the t₀ interval produced a similar bias in the BC and 3T methods and had the same outcome as incomplete sampling of the focal interval t₁ (Fig. 5, Supplementary Fig. S12). This is because removing a taxon from t₀ meant that it was no longer a boundary-crosser or “two-timer,” and was therefore not counted in the extinction estimate calculation. Rare taxa are preferentially removed as sampling becomes poorer, resulting in a relative overrepresentation of common taxa in the BC cohort in t₁ (the focal interval). The range through method was unaffected when the t₀ bin was undersampled, because this method calculates extinction using data from only the t₁ and t₂ bins.

Incomplete sampling of the t₂ interval had the largest effect on extinction estimate bias for the BC method and also strongly affected the accuracy of the range through method (Fig. 5, Supplementary Fig. S12). This occurs because failure to sample a taxon in t₂ means it will be erroneously marked as extinct in t₁, rather than marked as surviving. The extent of the bias also differed as a function of selectivity; if rare taxa were more likely to survive the extinction, they were also more likely to be missed when t₂ was incompletely sampled, inflating extinction estimates. The 3T extinction metric was unaffected by incomplete sampling of t₂ in this simulation, because the calculation is designed to use part-timer taxa (those present in t₁ and t₃, but missing from t₂) to correct for sampling incompleteness in t₂.

None of the three methods tested here are immune from bias due to abundance selectivity. However, the 3T method was generally less biased than the raw (range through) and BC methods, unless the t₂ interval was perfectly sampled. When the t₀ interval was poorly sampled but t₂ sampling was perfect, the range through method performed best.

The marine invertebrate fossil record for the Permian and Triassic varies considerably in sampling completeness across space, through time, and among clades (Supplementary Tables S2, S3). As a result, origination and extinction metrics could not be calculated for some of the spatiotemporal (latitude by stage) bins, particularly using the BC and 3T metrics, which have a higher sample size threshold (Fig. 6). Species-level analyses produced similar trends to those at the generic level but with generally higher levels of origination and extinction (Supplementary Fig. S13).

Origination and extinction appear to have been uniformly low across latitudes for marine invertebrates in the Permian, predominately below 50% of genera (Fig. 6). There is little evidence of heightened extinction in the Capitanian, although ammonoids underwent considerable diversification in the Wuchiapingian, particularly at low latitudes (around 80% of genera).

The effect of the PTME is clearly visible, with high extinction globally in the Changhsingian. Ammonoids experienced highest extinction levels at low latitudes (around 90%) and reduced extinction in the southern midlatitudes (around 40%). Brachiopods and bivalves appear to have exhibited the opposite gradient, with slightly lower extinction proportions at low latitudes (around 70% at low latitudes vs. 90% at high latitudes for brachiopods, and 40% at low latitudes vs. 65% at high latitudes for bivalves).

The brachiopod, bivalve, and gastropod records are sparse during the Early Triassic, particularly in the Southern Hemisphere, hindering reliable rate estimation. Origination appears not to have been unusually high in the Induan, in the aftermath of the mass extinction, except perhaps for the ammonoids (around 85%). Although extinction was generally reduced and spatially uniform in the Induan, brachiopods experienced a high proportion of extinction in the low northern latitudes. Ammonoids underwent both high origination and extinction globally in the Olenekian (both around 85%), indicating rapid turnover in the clade, whereas bivalves experienced reduced extinction rates (around 20%). In the Middle Triassic, origination and extinction gradients became more spatially uniform and stable for all four clades, although ammonoids appear to have undergone turnover in the high northern latitudes during the Ladinian.

Our results indicate that in the absence of abundance-based extinction selectivity, raw, BC, and 3T origination and extinction metrics can be applied to spatially subdivided datasets of fossil occurrences, producing estimates with only slightly less accuracy than when used globally (Fig. 3, Supplementary Fig. S2). This inaccuracy is likely partly due to the reduction in sample size (number of occurrences) resulting from creating spatial subdivisions within a global dataset. However, sampling incompleteness reduces the accuracy of origination and extinction rate estimation considerably (Supplementary Fig. S2b), as does the potential presence of extinction selectivity based on taxon abundance (Fig. 5).

The 3T estimates were generally more accurate than those calculated using raw or BC methods (Figs. 3, 5, Supplementary Fig. S2, Table 1). However, the 3T method requires occurrences to be present across multiple time slices (Fig. 1), a condition that may not always be met. BC estimates produced the same gradient of origination and extinction in latitude bins within a given iteration, and so can be used as a reasonable compromise when comparing estimates within a time slice. Raw estimates were slightly more successful at identifying spatial bins with the highest and lowest levels of origination and extinction (Supplementary Fig. S3). As a result, choice of metric should depend on the analysis required; we recommend calculating more than one metric, when possible, to allow for comparison.

The most accurate estimates were produced using the 3T method on spatial bins with a high sampling completeness (Fig. 5, Supplementary Fig. S2b), and variation in sampling completeness between spatial bins appears to hinder efforts to identify the bins that experienced highest and lowest origination and extinction (Supplementary Figs. S3, S9). Using higher taxonomic levels may increase the number of usable occurrences, but the patterns observed should not be assumed to be analogous to trends at lower taxonomic levels (Hendricks et al. 2014).

Estimates produced for low levels of origination and extinction were slightly more accurate than those for high levels of origination and extinction (Supplementary Fig. S8). This was expected for BC and 3T estimates; these methods exclude singletons, so when turnover rates increase, a larger proportion of species are represented as singletons in the fossil record, and the effective sample size is reduced. However, the difference in accuracy observed in the simulated data is relatively small, and therefore the uncertainty around high origination and extinction estimates should not be considered greater than that around smaller estimates.

Additional investigation using more complex simulations might reveal how the three rate metrics used here could be modified to produce more accurate results in a regional context, or whether a correction could be used to account for any regular overestimation in rate estimates (Supplementary Fig. S2b). Potential also exists for the development of metrics specific to estimating spatial variation in rates, including adaptation of the gap-filler and second-for-third methods, which perform better than raw, BC and 3T metrics when applied to global datasets (Alroy 2014, 2015). The resilience of these metrics to abundance-based extinction selectivity will also need to be tested. Simulation approaches might also hold potential for investigating how best to estimate migration rates from the fossil record.

Analyses of Permian–Triassic marine invertebrates provide little evidence for high extinction in any of the studied clades during the CBC (Fig. 6). Brachiopods experienced slightly higher extinction in the low and high northern latitudes (around 40% of genera), which supports the results of Clapham (2015), who identified a generic turnover in Capitanian–Wuchiapingian faunas from Iran and South China. Extinction in the clade was higher at the species level, particularly in the Northern Hemisphere (around 75%; Supplementary Fig. S13), which is somewhat lower than the 87% reported by Shen and Shi (1996) for South China.

Globally, low origination and extinction levels (both around 20%) indicate that bivalves appear to have been largely unaffected by the CBC, in agreement with Bond et al. (2010). By contrast, ammonoids and gastropods experienced moderate extinction levels at low latitudes (around 50%) during the CBC. Our results support a generic turnover of ammonoids during the Capitanian (Villier and Korn 2004; Rampino and Shen 2019), but the clade's fossil record is highly spatially restricted at this time, so it is unclear whether this turnover took place globally or only in the southern low latitudes. We found high ammonoid origination levels across low to midlatitudes in the Wuchiapingian, in agreement with Clapham et al. (2009) and McGhee et al. (2013). Collectively, the results presented here support the hypothesis that the CBC was highly selective (McGhee et al. 2013), with some marine clades undergoing turnover, particularly at low latitudes (ammonoids and gastropods), while others were seemingly unaffected (bivalves and brachiopods).

High extinction levels during the Changhsingian, which mark the PTME, are seen globally in all four marine invertebrate clades. For bivalves, gastropods, and brachiopods, the highest extinction levels were observed at mid- to high latitudes. These results support those of Reddin et al. (2019), who reported that benthic marine invertebrates experienced higher extinction rates away from the equator during the PTME. However, slightly reduced extinction levels at low latitudes in these clades may be an artifact caused by the occurrence of “mixed faunas.” These are communities of benthic invertebrates, particularly well known from sections in low-latitude South China, that survived the Changhsingian only to become extinct in the earliest Induan (Song et al. 2013). The presence of a low-latitude Induan extinction peak for brachiopods marks the loss of these mixed faunas. In contrast, ammonoids experienced the highest Changhsingian extinction levels at low latitudes. The late Permian ammonoid fossil record is relatively spatially restricted in comparison with the Early Triassic record, which may contribute to this result (Dai and Song 2020); that is, low extinction in the southern midlatitudes corresponds to a handful of localities in northern India and may therefore reflect relatively incomplete sampling.

Bivalves, gastropods, and ammonoids appear to have experienced highest origination at low latitudes during the Induan. The discrepancy between raw and BC rates for some ammonoid and brachiopod latitudinal bins during this interval is likely due to “disaster faunas,” which contributed a high number of temporal singletons. Ammonoids underwent rapid turnover across latitudes in the Olenekian, which coincides with an observed increase in their endemism during this time, and the influence of a possible extinction event around the Smithian/Spathian boundary (~250 Ma) (Brayard et al. 2006; Song et al. 2011; Wignall 2015; Dai and Song 2020). For all four clades, origination and extinction reduced and became more globally homogenous in the Anisian and Ladinian, corresponding to the amelioration of global temperatures and the re-establishment of diverse benthic communities (Sun et al. 2012; Martindale et al. 2019).

In summary, both origination and extinction were relatively consistent across latitudes for all four clades throughout the middle Permian to Middle Triassic (Fig. 6). Minor variation in rates was observed through time and between latitudinal bins, but caution should be exercised in interpreting any observed differences as indicative of contrasting evolutionary dynamics across latitudes rather than reflecting noise or sampling bias. The taxonomic identity of a genus appears to be a more important factor in determining extinction vulnerability than the paleolatitudes occupied, indicating a stronger phylogenetic signal in selectivity than that for latitude.

Although the scale of estimated variation between latitudinal bands is limited, the Changhsingian latitudinal gradients of extinction seen in benthic clades (bivalves, brachiopods, gastropods) align with the hypothesis that extinction during the PTME was most severe at high latitudes due to loss of suitable habitat for cool-adapted taxa, as advocated for by Penn et al. (2018). By contrast, the extinction gradient of ammonoids, the only pelagic clade examined here, better fits the hypothesis that extremely high equatorial temperatures resulted in tropical extinction and poleward migration, a mechanism proposed by Sun et al. (2012) and supported by Song et al. (2020) and Flannery-Sutherland et al. (2022). Our results therefore suggest a possible contrast in spatial patterns of PTME extinction based on paleoecology. Variation in the rate of post-PTME recovery has previously been identified between these two ecologies, with ammonoids and conodonts recovering more rapidly than benthic animals (Brayard et al. 2006; Song et al. 2011, 2018; Martindale et al. 2019), a trend also observed after the Triassic–Jurassic mass extinction (Dunhill et al. 2018). Penn et al. (2018) hypothesized that these two contrasting extinction gradients may reflect differences in primary kill mechanism, with extinctions at higher latitudes linked to oxygen depletion, and extinctions at low latitudes linked to temperatures beyond thermal tolerance levels. Oceanic oxygen depletion tends to be more severe in bottom waters (e.g., Liao et al. 2010), leaving benthic animals more vulnerable to anoxia than pelagic animals. In addition, the elevated extinction levels experienced by ammonoids during the late Smithian thermal maximum may support a higher sensitivity to extremes of temperature within the clade (Dai and Song 2020).

Although it is tempting to accept the above narrative, the results derive from an incomplete and biased fossil record, imperfect rate estimation metrics, and small differences between rate estimates across time and latitudes. Extinction selectivity based on abundance may have also reduced the accuracy of our estimates. As such, we advocate for additional assessment of latitudinal variation in Permian–Triassic evolutionary rates in the future, with larger datasets and more robust metrics, ideally designed specifically to estimate variation in evolutionary rates across space. Investigating spatial variation in origination and extinction in other pelagic clades, such as conodonts and fishes, and directly testing the timing and direction of migration in marine invertebrates, may also provide further insight into the nature of biotic responses to the PTME in the marine realm.

Using simulations, we demonstrate that it is possible to estimate origination and extinction patterns across space from the fossil record with reasonable accuracy, albeit under conditions of relatively complete sampling and no extinction selectivity based on taxon abundance. We applied rate estimation metrics to marine invertebrates from the Permian and Triassic to test previously posited hypotheses on the latitudes at which the highest extinction rates occurred during the PTME. Both origination and extinction were found to be generally uniform across latitudes throughout the study interval. Larger fossil datasets with better geographic spread and more robust metrics would be necessary to confirm whether any deviations in origination and extinction are not simply a result of sampling bias, and of a sufficient magnitude to be deemed evolutionarily important. Our results indicate that origination and extinction levels during the PTME were more dependent on taxonomic identity than latitudinal occupation.

We thank members of the Palaeo@Leeds, Ecosystem Resilience and Recovery from the Permo-Triassic Crisis (Eco-PT), and Computational Evolution (cEvo) research groups for discussion, and P. Mannion and T. Aze for their constructive criticism of earlier versions of this article. We are grateful to M. Powell and an anonymous reviewer, whose comments greatly improved this article. We thank all those who have contributed toward the Paleobiology Database collections used in this study. This is Paleobiology Database publication no. 440. B.J.A. commenced this work during her Ph.D. at the University of Leeds, funded by a U.K. Natural Environment Research Council studentship (NE/L002574/1) and has since been supported by an ETH+ grant (BECCY) funded by ETH Zurich. P.B.W and A.M.D. were both supported by the Eco-PT grant (NE/P0137224/1), part of the Biosphere Evolutionary Transitions and Resilience (BETR) project, jointly funded by the U.K. Natural Environment Research Council and National Natural Science Foundation of China. E.E.S. was supported by a Leverhulme Prize and Leverhulme grant (DGR01020) and NERC grant (NE/V011405/1).

The authors declare no competing interests.

Data available from the Zenodo Digital Repository: https://doi.org/10.5281/zenodo.7355273.