Post-extinction recovery of the Phanerozoic oceans and the rise of biodiversity hotspots


 The fossil record of marine invertebrates has long fueled the debate on whether or not there are limits to global diversity in the sea1–4⁠. Ecological theory states that as diversity grows and ecological niches are filled, the strengthening of biological interactions imposes limits on diversity5–7⁠. However, the extent to which biological interactions have constrained the growth of diversity over evolutionary time remains an open question1–4,8–12⁠, largely because of the incompleteness and spatial heterogeneity of the fossil record13–15⁠. Here we present a regional diversification model that reproduces surprisingly well the Phanerozoic trends in the global diversity of marine invertebrates after imposing mass extinctions. We find that the dynamics of global diversity is best described by a diversification model that operates broadly within the exponential growth regime of a logistic function. A spatially resolved analysis of the diversity-to-carrying capacity ratio reveals that only < 2% of the global flooded continental area exhibits diversity levels approaching ecological saturation. We attribute the overall increase in global diversity during the Late Mesozoic and Cenozoic to the development of diversity hotspots under prolonged conditions of Earth system stability and maximum continental fragmentation. We call this the "diversity hotspots hypothesis", which is proposed as a non-mutually exclusive alternative to the hypothesis that the Mesozoic marine revolution led this macroevolutionary trend16,17.


Abstract
The fossil record of marine invertebrates has long fueled the debate on whether or not there are limits to global diversity in the sea 1-4 . Ecological theory states that as diversity grows and ecological niches are filled, the strengthening of biological interactions imposes limits on diversity [5][6][7] . However, the extent to which biological interactions have constrained the growth of diversity over evolutionary time remains an open question 1-4, [8][9][10][11][12] , largely because of the incompleteness and spatial heterogeneity of the fossil record [13][14][15] . Here we present a regional diversification model that reproduces surprisingly well the Phanerozoic trends in the global diversity of marine invertebrates after imposing mass extinctions. We find that the dynamics of global diversity is best described by a diversification model that operates broadly within the exponential growth regime of a logistic function. A spatially resolved analysis of the diversityto-carrying capacity ratio reveals that only < 2% of the global flooded continental area exhibits diversity levels approaching ecological saturation. We attribute the overall increase in global diversity during the Late Mesozoic and Cenozoic to the development of diversity hotspots under prolonged conditions of Earth system stability and maximum continental fragmentation.
We call this the "diversity hotspots hypothesis", which is proposed as a non-mutually exclusive alternative to the hypothesis that the Mesozoic marine revolution led this macroevolutionary trend 16,17 .

Main text
The question of whether or not there is an equilibrium diversity that the biota, or portions of the biota, cannot exceed has led to decades of debate between those who think that there is a limit to the global diversity that the Earth can carry 2,3,11,18 (i.e., a carrying capacity or saturation level) and those who think that diversity can increase in an unlimited fashion over time or, alternatively, that the biosphere is so far from the equilibrium diversity (i.e., its carrying capacity) that we can ignore the existence of any limit 8,9,12 . This question has traditionally been addressed by examining the shape of global fossil diversity curves 3,19 . For example, the Paleozoic plateau in marine invertebrate diversity is generally taken as strong evidence for the existence of ecological limits to further diversification 3,20 . However, because diversity varies dramatically among geographic regions, and each geographic region has its own geological and environmental history, addressing this question requires simultaneously reconstructing the dynamics of regional diversity in both space and time 14,21 . If diversity dynamics were governed by diversity-dependent feedbacks on speciation and extinction rates, then regional diversity should remain stable regardless of time once carrying capacity had been reached (i.e., the logistic model). The reasoning is the same as that used to explain the logistic growth model in population dynamics in which the per capita rate of increase decreases as the population approaches its maximum size or carrying capacity. Conversely, if evolutionary rates were independent of standing diversities, then we should observe positive relationships between evolutionary time-within-regions (or time-for-speciation) and diversity; the older the habitat the longer the lineages have had to diversify and fill empty niches or explore new ones (i.e., the exponential model). The reasoning in this case is the same as that used to explain the exponential growth model in population dynamics in which the per capita rate of increase does not depend on the population size but only on the modulating effects of environmental conditions. Determining which diversification model best 3  describes the dynamics of regional diversity over time is key to understanding the mechanisms underlying biogeographic patterns and macroevolutionary trends. However, the fossil record is biased by the inequality of the geographic and stratigraphic sampling effort 13,14 , and the inequality in the rock record available for sampling 22 , hindering our ability to investigate the effect of geographic variability in evolutionary time and diversification rate.
In order to overcome this limitation, we couple two alternative models of diversification, logistic and exponential, to a global plate motion model that constrains evolutionary timewithin-regions (i.e. the age of the seafloor for the deep ocean and the time underwater for the flooded continental regions), Then, we reconstruct the spatial distributions and time trajectories of marine benthic animal divesity throughout the Phanerozoic. In both diversification models, the net diversification rate varies within a fixed range of values as a function of seawater temperature and food supply, which are reconstructed using a spatiallyexplicit Paleo-Earth system model (see Methods). The effects of temperature and food supply are parameterized, respectively, using a Q10 coefficient and a non-linear food limitation factor, under the premise that increasing temperature and food supply increase the rate of genus origination by shortening individuals' generation times and increasing population sizes, respectively. In the logistic model, the spatially resolved effective carrying capacities (Keff) are allowed to vary within a fixed range of values (Kmin and Kmax) as a positive linear function of the food supply in each ocean region and time. That is, those regions with the lowest and highest food supply are assigned Kmin and Kmax, respectively, while those regions with intermediate levels of food supply are assigned Keff values within the range Kmin to Kmax, accordingly. Finally, mass extinctions are imposed by imputing negative net diversification rates to regional communities and assuming non-selective extinction. The percentage of diversity loss as well as the starting time and duration of mass extinctions are extracted from three fossil diversity curves of reference, namely Sepkoski 23 , Alroy 24  of these curves provide alternative insights into the Phanerozoic history of marine animal diversity based on uncorrected range-through genus richness estimates 23,25 and sampling standardized estimates 24 .

Seeking diversity hotspots in Phanerozoic oceans
There are clear differences in the spatial distributions of diversity generated by the logistic model and the exponential model as illustrated for 4 representative time-slices in the Phanerozoic (Figure 1 and Extended Data Figures 1 and 2; see Supplementary Video s 1 and 2 (password: video2021) for the full Phanerozoic sequences). Regardless of the diversification model, most of the diversity is concentrated in shallow marine environments, where high temperatures and abundant food supplies increase the rates of diversification compared to deep-sea benthic habitats (Fig. 1, Extended Data Fig. 1 and 2). However, while regional diversity increases unconstrained in the exponential model, carrying capacities limit the growth of regional diversity in the logistic model, preventing the development of ocean regions with exceptionally high levels of diversity, hereinafter diversity hotspots.
Diversity hotspots occur in tropical shelf seas of the Early Devonian, Permian, Late Cretaceous and Cenozoic ( Fig. 1e- Pacific provinces (Fig. 1g, h, Extended Data Fig. 1g, h and 2g, h). This temporal trend in the prominence of diversity hotspots cannot be explained by a secular increase in the maximum lifetime of shelf seas, a proxy for the maximum potential evolutionary time-withinregions. Sedimentary data (i.e., magnetic anomalies in ancient continental margins trapped within orogenic belts) 26 and global tectonic reconstructions 27 , including our reconstruction ( Supplementary Fig. 1), show no evidence of an increase in the lifespan of passive continental margins or in the maximum ages of the seafloor over the Phanerozoic. Rather, we argue that the temporal proximity between the Ordovician-Silurian (Hirnantian), Late Devonian (Frasnian-Famennian), and Permian-Triassic mass extinctions, coinciding with a long-lived phase of continental coalescence and coastline destruction, interrupted the full development of diversity hotspots during the Paleozoic. By contrast, the comparatively long expanse of time that separated the mass extinctions of the end-Triassic and end-Cretaceous extended the time-for-speciation under conditions of increasing continental fragmentation, giving rise to exceptionally high diversity regions before the Cretaceous-Paleogene mass extinction. The extraordinary diversity of Late Cretaceous hotspots ensured the continuity of relatively high diversity levels in the aftermath of the end-Cretaceous mass extinction, facilitating the subsequent development of diversity hotspots during the Cenozoic.
Consistent with the patterns emerging from the fossil record of animal-like protists, such as large benthic foraminifera 28 , the exponential model is able to reproduce the biogeography of diversity hotspots across the Late Cretaceous and Cenozoic (Fig. 1g, h, Extended Data Fig.   1g, h and 2g, h, and Supplementary Video 2 ) (password: video2021). Our analysis suggests that these diversity hotspots were a consequence of the long residence time of shallow water seas within the tropical belt, leading to the establishment of three main tropical high diversity loci: the Paratethys-Mediterranean Sea, the Atlantic Caribbean-East Pacific, and the Indo-West Pacific. These spatial distribution patterns are remarkably consistent with

Reconstructing global diversity dynamics
Each of the two diversification models tested here produces a total of 82 spatially-explicit reconstructions of diversity spanning from the Cambrian to the present. On each of the diversity distribution maps, we trace hundreds of line transects from diversity peaks to their nearest diversity troughs and integrate the total diversity in each transect by assuming a decay function in taxonomic similarity with geographic distance (see Material and Methods and Supplementary Fig. 2). Then, for each of the 82 time intervals, all integrated diversities along transects are re-integrated step-wise, from the transect with the greatest diversity to the transect with the lowest one, assuming the same distance-decay function applied to individual transects. The resulting global diversity estimates are plotted against the mid point value of the corresponding time interval to generate a synthetic global diversity curve. Both the logistic model and the exponential model produce relatively similar global diversity dynamics over time (Fig. 2, Extended Data Fig. 3). This was to be expected since the global diversity curves produced by both models are equally influenced by long-term variations in the global area of tropical shallow shelf seas 31,32 (Extended Data Fig. 4), which harbour the vast majority of the diversity of marine benthic animals. However, while both models show similar diversity dynamics, the amplitude of global diversity variations differ markedly between models depending on whether or not regional-scale diversities self-limit their increase over time. The exponential model gives rise to conspicuous increases in global diversity from the Cambrian to Late Ordovician, Silurian to Early Devonian, Carboniferous (Lower to Upper Pennsylvanian), Early to Late Cretaceous, and Paleocene to present. The Permian-Triassic mass extinction event lowered global diversity to Early Paleozoic levels, but later diversification led Late Cretaceous and Neogene faunas to exceed the Mid-Paleozoic global diversity peak. These trends emerge consistently regardless of the mass extinctions pattern imposed, be it Sepkoski 23 , Alroy 24 , or Zaffos et al. 25 (Fig. 2a, b, c, respectively). Nevertheless, it is worth noting how well the exponential model reproduces the 'uncorrected' Sepkoski fossil diversity curve. The logistic model also reproduces the initial increase in diversity, from the Cambrian to Upper Ordovician and from the Silurian to Early Devonian (Fig. 2). However, unlike the exponential model, in the logistic model this initial upward trend is followed by a convex diversity pattern interrupted by a modest increase during the Cretaceous, which rarely exceeds the mid-Paleozoic global diversity peak in our set of simulations (Fig. 2).
Our logistic model allows the spatially-resolved effective carrying capacities (Keff) to vary within a fixed range of values (from Kmin to Kmax) as a positive linear function of food availability in each ocean region and time. All other things being equal, the higher the Kmin and Kmax values, the longer the evolutionary time required to reach diversity saturation.
Consequently, the choice of Kmin and Kmax critically influences the extent to which regional biotas reach saturation. In order to calibrate the Kmin and Kmax parameters, we run simulations of pair-wise Kmin and Kmax combinations in a geometric sequence of base 2, from 2 to 256 genera, and test the effect of changing the Kmin and Kmax values on the concordance between the normalized diversities generated by the model and those estimated from the fossil record ( Fig. 3). Unlike other correlation coefficients, the Lin's concordance correlation gold standard, that is, the 1:1 line. We focus the analysis on the time series data between the end of one mass extinction and the beginning of the next, that is, considering those time intervals dominated by rising diversity trajectories, which are influenced by the mode of regional diversification. The Lin's CCC increases with increasing Kmin and Kmax until reaching a plateau except for the mass extinction pattern of Sepkoski for which it continues to increase even at the highest Kmin and Kmax values (Fig. 3, Extended Data Fig. 5, 6, and 7). These results are consistently replicated using alternative values for the parameters of the model that define the temperature-and food-dependence of the net diversification rate (Fig. 3, grey lines in insets, Extended Data Table 1).
High Kmin and Kmax values (and thus high Keff), imply the need for longer evolutionary times and/or higher diversification rates to reach diversity saturation. Therefore, the results of the calibration analysis suggest that, in broad regions of the ocean, diversity could have been systematically far from saturation. In order to corroborate it, we re-run the logistic model using and analyse the spatial and temporal variability of the diversity-to-carrying capacity (Keff) ratio.
This ratio represents a quantitative index of how far (ratios close to zero) or how close (ratios close to one) are the regional faunas from noticing the effect of diversity-dependent ecological factors. The diversity-to-Keff ratio falls below 0.25 in most of the ocean and throughout the Phanerozoic ( Fig. 4a- Table 1), and represent its frequency distributions ( Fig. 4m-o, Extended Data Fig. 8). Most of the estimates fall within the exponential growth regime of the logistic function (i.e. diversity-to-Keff ratio < 0.25). On average, less than 10% of the estimates exceed the threshold of 0.25, and only < 2% of the estimates, those associated with well-developed diversity hotspots, exceed the threshold of 0.5.
It is unlikely that our calibration analysis leads to an overestimation of Keff since the CCC threshold of 0.7 imposed for the calculation of Kmin and Kmax would have led, if anything, to an understimation of Kmin and Kmax and thus, Keff. Nonetheless, a deliberate decrease of 25% in the Kmin and Kmax values of the model does not alter significantly the shape of the diversity-to-Keff ratio frequency distributions (Extended Data Fig. 9), which indicates that the resulting patterns are indeed robust. Furthermore, the paucity of diversity hotspots during the Paleozoic points to a scenario in which the relatively short time elapsed between successive mass extinctions interrupted their development. In fact, by deactivating the Late Devonian mass extinction in the model, we find that the full development of diversity hotspots before the end of the Permian leads to global diversities two to four times greater than those generated by the same 'calibrated' logistic model but with all mass extinctions enabled (Extended Data   Fig. 10). We cannot reject the hypothesis that diversity saturation slowed down diversification processes in ocean regions where diversity hotspots i) had long development times, ii) evolved fast or iii) evolved from relatively high initial diversities (Fig. 4). Nevertheless, the slowdown and eventual halt of diversity growth at hotspots would not have prevented global diversity from continuing to grow as new diversity hotspots emerged elsewhere. We have shown that the development of biodiversity hotspots accounts for much of the increase in the global diversity of marine benthic animals throughout the Phanerozoic. Our analysis reveals that the temporal proximity between successive mass extinction events, along with the long-term reduction in the length of the coastline during the assembly of Pangea, interrupted the development of diversity hotspots during the Paleozoic. We find evidence of regional biota approaching diversity saturation at post-Paleozoic diversity hotspots, the development of which helps explain the increase in global diversity during the Late Mesozoic and Cenozoic. We call this the 'diversity hotspots hypothesis', which is proposed as an non-mutually exclusive alternative or supplement to the hypothesis that the Mesozoic marine revolution, that is, the evolutionary emergence of the durophagous predators and the ensuing cascade diversification, led this macroevolutionary trend 16,17 .
With the possible exception of well-developed diversity hotspots, our results indicate that the diversity of marine benthic animals has remained well below saturation levels throughout their evolutionary history, shedding light on one of the most controversial questions in evolutionary ecology 2,3,8,9,11,12,18,34 . A taxonomic diversification model operating widely within the exponential growth regime of the logistic function implies a concave-upward relationship between the magnitude of diversity loss (x-axis) and the subsequent rebuilding time. This mode of diversification provides the most plausible explanation for the observed decoupling between mass extinctions and evolutionary radiations 35 over the Phanerozoic. We envision that our spatially-explicit reconstructions of diversity may shed light on other long-standing questions in (paleo)biogeography and macroevolution.   times. In these cases, we assume that the seafloor has a zero-age at the time the intraoceanic arc first develops, then remains predominantly underwater for the rest of its lifetime.
Supplementary figure 1a-d shows the estimated age of the seafloor for open ocean and flooded continental shelves using the approach described in this section.
Therefore, for each of the 82 paleogeographic reconstructions, we annotate 0.5º by 0.5º grids as continental, flooded continental shelf, or oceanic for later use in model coupling and production of regional diversity maps.

Paleo-environmental conditions: cGenie Earth System model
We use cGENIE 41  for the analysis.

Regional diversification model
We test two models of diversification, the logistic model and the exponential model, describing the dynamics of regional diversity over time. In both models, the net diversification rate (ρ), with units of inverse time (Myr -1 ), is the difference between the rates of origination and extinction and varies within a pre-fixed range of values as a function of seawater temperature and food availability. The implicit mechanism is that high temperatures and abundant food supplies increase the genus origination rates by shortening individual's generation times (i.e. higher metabolic rates) and increasing population sizes (i.e. higher where ρmin and ρmax set the lower and upper net diversification rate limits within which ρ is allowed to vary, and Qtemp and Qfood are non-dimensional limitation terms with values between 0 and 1 that define the dependence of ρ on temperature and food, respectively (see Supplementary Table 1). These temperature and food supply limitation terms vary in space and time as a result of changes in seawater temperature and particulate organic carbon export rate, respectively, thereby controlling the spatial and temporal variability of ρ.
The temperature-dependence of ρ is calculated using the following equation: where the Q10 coefficient measures the temperature sensitivity of the origination rate.
Assuming a constant background extinction rate, a Q10 coefficient of 2 would correspond to a doubling of the net diversification rate for every 10 ºC increase in temperature. In the The food limitation term is parameterized using a Michaelis-Menten formulation as follows: where POC flux (mol m -2 year -1 ) is the particulate organic carbon export flux, which is used as a surrogate for food availability, at a given location and time of the simulated seafloor. The Zaffos et al 25 (Supplementary Fig. 4). Each of these fossil diversity curves provides different insights into the Phanerozoic history of marine animal diversity based on uncorrected rangethrough genus richness estimates 23,25 and sampling standardized estimates 24 . Regional-scale processes, such as sea level fall during marine regressions and/or seafloor destruction at plate boundaries (either by subduction or uplift), are simulated by the combined plate tectonic/paleo-elevation model, and constrain the time span that seafloor habitats have to accumulate diversity.
Letting D represent regional diversity (number of genera within a given seafloor point) and t represent time, the logistic model is formalized by the following differential equation: where D(t) is the number of genera at time t and Keff is the effective carrying capacity or maximum number of genera that a given seafloor point (i.e., grid cell area after gridding) can carry at that time, t. In our logistic model, Keff is allowed to vary within a fixed range of values (from Kmin to Kmax) as a positive linear function of the POC flux at a given location and time as follows: where In the logistic model, the net diversification rate decreases as regional diversity approaches its Keff. The exponential model is a particular case of the logistic model when Keff approaches infinity and, therefore, neither the origination rate nor the extinction rate depend on the standing diversities. In this scenario, diversity grows in an unlimited fashion over time only truncated by externally imposed mass extinctions and/or by the dynamics of the seafloor (creation versus destruction). The exponential model is thus as follows: where the rate of change of diversity (the time derivative) is proportional to the standing diversity D such that the regional diversity will follow an exponential increase in time at a speed controlled by the temperature-and food-dependent net diversification rate. Even if analytical solutions exist for the steady-state equilibrium of the logistic and exponential functions, we solved the ordinary differential equations (4) and (6)  Estimation of global diversity from regional diversity Our regional diversity maps are generated by separately interpolating ocean point diversity and flooded continental point diversity into the 0.5º by 0.5º annotated grids provided by the paleogeographic model. We calculate global diversity at each time step from each of the regional diversity maps following a series of steps to integrate diversity along line transects from diversity peaks (maxima) to diversity troughs (minima) (Supplementary Fig. 2). To select the transects, first, we identify on each of the regional diversity maps the geographic position of the diversity peaks. The peaks are defined as the grid cells identified as local maxima (i.e., with diversity greater than their neighbour cells). In the case of grid cells with equal neighbour diversity, the peak is assigned to the grid cell in the middle point of the grids with equal diversity. We subsequently identify the geographic position of the diversity troughs, which are defined as newly formed ocean grid cells (age = 0 Myr) and, therefore, with diversities equal to one. The troughs are mostly located at mid-ocean ridges.
On each of the 82 spatial diversity maps, we trace a line transect from each diversity peak to its closest trough, provided that the transect does not cross land in more than 20 % of the grid cells along the linear path. On average, for each spatial diversity map, we trace 400 (σ = ±75) linear transects. This sampling design gives rise to transects of different lengths, which may bias the estimates of global diversity. To minimize this bias, we cut the tail of the transects to have a length of 555 km equivalent to 5º at the equator. We test an alternative cutoff threshold; 1110 km, and the results do not alter the study's conclusions.
We apply Bresenham's line algorithm 50 to detect the grid cells crossed by the transects and annotate their diversity. To integrate regional diversity along the transects, we develop a method to simplify the scenario of peaks and troughs heterogeneously distributed in the 2D diversity maps. The method requires i) a vector (the transect) of genus richness (αn) at n different locations (grids) arranged in a line (1D) of L grids, and ii) a coefficient of similarity (V , +1) between each two neighbouring locations, n and n+1. Vn,n+1, the coefficient of similarity, follows a decreasing exponential function with distance between locations. The number of shared genera between n and n+1 is V , +1 * min(α ; α +1). We integrate diversity from peaks to troughs and assume that, along the transect, α +1 is lower than α . We further 25  assume that the genera present in n and n+2 cannot be absent from n+1. Using this method, we integrate the transect's diversity (γi) using the following equation: To integrate the diversity of all transects (γi) on each 2D diversity map (or time slice), we apply the same procedure as described above (Supplementary Fig. 2). We first sort the transects in descending order from the highest to the lowest diversity. Then, we assume that the number of shared genera between transect i and the rest of transects with greater diversity {1, 2,…, i -1} is given by the distance of its peak to the nearest neighbour peak [NN(i)] of those already integrated {1,2,… i -1}. Thus, we perform a zigzag integration of transects' diversities down gradient, from the greatest to the poorest, weighted by the nearest neighbour distance among the peaks already integrated. As a result, the contribution of each transect to global diversity will depend on its diversity and its distance to the closest transect from all those transects already integrated. With this method, we linearize the problem to simplify the cumbersome procedure of passing from a 2D regional diversity map to a global diversity estimate without knowing the identity (taxonomic affiliation) of the genera. Being γtotal the global diversity at time t: Finally Lin's CCC combines measures of both precision and accuracy to determine how far the observed data deviate from the line of perfect concordance or gold standard (that is, the 1:1 line). Lin's CCC increases in value as a function of the nearness of the data's reduced major axis to the line of perfect concordance (the accuracy of the data) and of the tightness of the data around its reduced major axis (the precision of the data).

Model parameterization and calibration
The diversification models are parameterized assuming a range of values that constrain the lower and upper limit of the genus-level net diversification rate (ρmin and ρmax, respectively)  (Figure 3, Extended Data Fig. 5, 6,   and 7). We perform these simulations independently for each of the fifteen parameter 27  settings selected previously (Extended Data Fig. 1). Each combination of Kmin and Kmax produces a global diversity curve, which is evaluated as described above using Lin's CCC.
Calculating estimates of global diversity from regional diversity maps in the absence of information on genus-level taxonomic identities requires assuming a spatial turnover of taxa with geographic distance (Supplementary Fig. 2). Distance-decay curves are routinely fitted by calculating the ecological similarity (e.g. Jaccard similarity index) between each pair of sampling sites, and fitting an exponential decay function to the points on a scatter plot of similarity (y-axis) versus distance (x-axis). Following this method, we fit an exponential decay function to the distance-decay curves reported in Miller et al 51  where Joff = 0.06 (n.d.) is a small offset, Jmax = 1.0 (n.d.) is the maximum value of the genusbased Jaccard similarity index, and λ= 0.0024 (Km -1 ) is the distance-decay rate.
The Jaccard similarity index (J) between consecutive points n and n+1 is bounded between 0 and min(α ; α +1)/max(α ; α +1). A larger value for J would mean that there are more shared genera between the two communities than there are genera within the least diverse community, which is ecologically absurd. However, using a single similarity decay function can lead the computed value of J to be locally larger than min(α ; α +1)/max(α ; α +1). To prevent this artifact, we use the Simpson similarity index or "overlap coefficient" (V) instead of J. V corresponds to the percentage of shared genera with respect to the least diverse community (min(α ; α +1)). V is bounded between 0 and 1, whatever the ratio of diversities. As the pre-existing estimates of similarity are expressed using J 51 , we make the conversion from J to V using the algebraic expression V = (1 + R) * J / (1 + J) where R = max(α ; α +1)/min(α ; α +1) (see Annex 1). In the cases in which J exceeds the min(α ; α +1)/max(α ; α +1), V becomes > 1 and, in those cases, we force V to be <1 by assuming R = 1, that is α = α +1.

Fossil data
We digitized three fossil diversity curves of reference, namely, Sepkoski 23 , as depicted in    Table 1). This range of carrying capacity values is arbitrarily selected to emphasize the differences between the logistic model and the exponential model.
The same plots but after imposing the mass extinction patterns extracted from the fossil diversity curves of Alroy 24  *The CCCs are for the relationship between the normalized diversities estimated from the fossil record and those generated by the exponential (Exp) and the logistic (Log) models. The fifteen combinations of model parameters that gave the highest CCC for each mass extinction pattern were selected. Of these, the combination that gave the highest CCC for the relationship between the fossil diversities and the diversities generated by the calibrated logistic (Cal. Log) model was selected as the best (Extended Data Table 1