Nonequilibrium evolution of volatility in origination and extinction explains fat-tailed fluctuations in Phanerozoic biodiversity.

Fluctuations in biodiversity, large and small, pervade the fossil record, yet we do not understand the processes generating them. Here, we extend theory from nonequilibrium statistical physics to describe the fat-tailed form of fluctuations in Phanerozoic marine invertebrate richness. Using this theory, known as superstatistics, we show that heterogeneous rates of origination and extinction between clades and conserved rates within clades account for this fat-tailed form. We identify orders and families as the taxonomic levels at which clades experience interclade heterogeneity and within-clade homogeneity of rates, indicating that families are subsystems in local statistical equilibrium, while the entire system is not. The separation of timescales between within-clade background rates and the origin of major innovations producing new orders and families allows within-clade dynamics to reach equilibrium, while between-clade dynamics do not. The distribution of different dynamics across clades is consistent with niche conservatism and pulsed exploration of adaptive landscapes.

INTRODUCTION

Biodiversity has not remained constant nor followed a simple trajectory through geologic time (–). Instead, it has been marked by fluctuations in the richness of taxa, both positive in the case of net origination and negative in the case of net extinction. Major events, such as adaptive radiations and mass extinctions, have received special attention (, ), but fluctuations of all sizes are ubiquitous (, ) and follow a fat-tailed distribution, where large events are more probable compared to, e.g., a Gaussian distribution. Understanding the fat-tailed nature of these fluctuations continues to elude paleobiologists and biodiversity theoreticians.

The fat-tailed distribution of fluctuations in taxon richness inspired earlier researchers to invoke ideas from complex systems with similar distributions. These ideas include the hypotheses that biological systems self-organize to the brink of critical phase transitions (, ) and that environmental perturbations are highly nonlinear (). More recent data and analyses have not, however, supported these hypotheses at the scale of the entire Phanerozoic marine invertebrate fauna (, ). Other studies have modeled the mean trend in taxon richness as tracking a potentially evolving equilibrium (, , ) and yet ignore the role of stochasticity and nonequilibrium dynamics in producing observed patterns (, , ). Individual-, population-, and local ecosystem–scale processes that could produce complex dynamics, such as escalatory coevolutionary interactions (), have not been documented to scale up to global patterns () and should not be expected to do so (). Thus, we still lack a theory to describe the striking fat-tailed nature of fluctuations throughout the Phanerozoic.

Despite the heterogeneity of explanations of Phanerozoic biodiversity, consensus has emerged on one property of macroevolution: Clades experience different rates of morphological evolution, origination, and extinction (, , , ). Here, we show that the simple fact of conserved rates within clades and variable rates across clades is sufficient to describe pervasive, fat-tailed fluctuations in taxonomic richness throughout the marine Phanerozoic. This biological mechanism has a precise correspondence to the nonequilibrial theory from statistical physics known as “superstatistics” (), which has been applied across the physical and social sciences (, ). We leverage this correspondence to explain the distribution of fluctuations in the standing richness of marine invertebrates preserved in the Phanerozoic fossil record. We further show that the specific mathematical form of this superstatistical distribution is consistent with niche conservatism (, ) and pulsed exploration on an adaptive landscape by higher taxa (, –). We operationally define “adaptive landscape” to mean a clade’s set of characteristics, and the fitness they impart to the clade, that influences its macroevolution. Those characteristics could be ecological [e.g., substrate preference (, , )], morphological [e.g., body plan ()], or macroecological [e.g., range size (, )].

Superstatistics of fossil biodiversity

Superstatistics () proposes that nonequilibrial systems can be decomposed into many local subsystems, each of which attains a unique dynamic equilibrium. The evolution of these dynamic equilibria across subsystems occurs more slowly. This separation in time scales allows local systems to reach equilibrium, while the system as a whole is not (). In the context of macroevolution, we propose that a clade with conserved macroevolutionary rates corresponds to a subsystem in dynamic equilibrium.

In statistical mechanics, local subsystems can be defined by a simple statistical parameter β often corresponding to inverse temperature. In macroevolutionary “mechanics,” we define β of clade k as the inverse variance of fluctuations x in the number of genera within that clade, i.e., fluctuations in the genus richness. β values thus represent the inverse variances, what we term volatilities, of the origination-extinction processes of genera with clades. The details of this origination-extinction process, e.g., whether it is linear or subject to a carrying capacity, are not of central importance to our analysis; so long as fluctuations can reach a stationary distribution and are observed over time-averaged intervals in a temporally coarse-grained fossil record, they will be approximately Gaussian [see section S1; ()].

We make the hypothesis of dynamic equilibrium within a clade following MacArthur and Wilson () in recognition that while the identity and exact number of taxa will fluctuate stochastically from random origination and extinction (taking the place of local immigration and extinction in island biogeography), the overall process determining the number of taxa, and, by extension, fluctuations in that number, is in equilibrium. The different regions of adaptive space occupied by different clades can be conceptualized as islands with unique dynamic equilibria, albeit with macroevolutionary processes determining the “colonization” of adaptive peaks, as opposed to short time scale biogeographic processes.

The volatility of richness fluctuations will vary across these islands in adaptive space as an emergent trait of a clade resulting from the macroevolutionary fitness of the clade and the shape of the surrounding adaptive landscape. Ultimately, volatility emerges from the life histories, ecologies, and evolutionary histories that drive each clade’s macroevolutionary fitness and characterize its occupancy of different regions of an adaptive landscape. We do not attempt to diagnose which characteristics of different regions account for volatility differences, but others have found rates of origination and extinction to depend on larval type (), body plan (), body size (), range size (, ), and substrate preference (). Not all of these traits would be considered dimensions of an ecological niche or characteristics of a guild (, , ), but they all point to different strategies that influence a clade’s macroevolutionary success. These characteristics result from interactions between heritable traits and environments, which themselves may be thought of as semi-heritable (). Thus, different regions of adaptive space, and the clades occupying them, will experience different magnitudes of stochastic fluctuations in taxonomic richness. As clades occasionally split to fill new regions of adaptive space, their pulsed diversification determines the nonequilibrium nature of the entire biota.

Real paleontological data to test superstatistics

To uncover the superstatistical nature of the marine invertebrate Phanerozoic fauna, we analyzed the distribution of fluctuations in genus richness (the lowest reliably recorded taxonomic resolution) using the Paleobiology Database (PBDB; paleobiodb.org). We corrected these raw data for incomplete sampling and bias using a new method described in Materials and Methods. Occurrences from the PBDB were matched to 49 standard time bins, all of approximately 11–million year (Ma) duration following previous publications (, ). Fluctuations in genus richness were calculated as the simple difference between bias-corrected richnesses in adjacent time bins.

To focus attention on the variance of fluctuations, we zero-centered each clade’s fluctuation distribution. In this way, we focus on fluctuations about any possible trend toward net diversification or extinction. Because “equilibrium” in the statistical mechanical sense means that a system undergoes coherent, concerted responses to perturbation, the mean trend line (positive or negative) is of less interest than deviations from it. We also note that the distributions of fluctuations for most clades are already very close to a mean of 0 (mean at the family level: 0.038 ± 0.176 SD), and so centering has little influence on clade-specific fluctuation distributions, consistent with the observation that origination is often roughly equal to extinction (). Following (), we also ignore all instances of no change (i.e., zero fluctuation).

We define potentially equilibrial subsystems based on taxonomic hierarchies, as a full phylogenetic hypothesis for all marine invertebrates is lacking. Taxa ideally represent groups of organisms that descend from a common ancestor and share similar ecologically and evolutionary relevant traits (, ). Thus, our model assumes that at a given higher taxonomic level, within-taxon fluctuations in richness are driven by equilibrial processes characterized by Gaussian distributions. We further assume that new higher taxa arise due to the emergence of sufficiently novel traits (be they ecological, morphological, life history, or macroecological) so that those new taxa occupy a new region of an adaptive landscape. We lastly assume that different regions of adaptive space are characterized by different volatilities in origination and extinction.

To evaluate the optimal taxonomic level for subsystem designation, we test our superstatistical theory using taxonomic levels from family to phylum. In addition, we compare our results to randomized taxonomies and confirm that the observed fit of superstatistics is not an artifact of arbitrary classification but instead represents real, biologically relevant diversification processes within and between clades. We find that families and orders conform to the assumptions of our superstatistical model, while classes and phyla do not.

RESULTS

We first evaluate the local equilibria of clades from family level to phylum. We find that family-level fluctuation distributions are well approximated by Gaussians (Fig. 1 and fig. S3). Three exemplar family-level dynamics are highlighted in Fig. 1 to illustrate how different volatility equilibria express themselves as actual richness time series. This Gaussian approximation also largely holds for orders, but classes and phyla increasingly show deviations from Gaussian with greater kurtosis, corresponding to more frequent outliers at these taxonomic levels (fig. S3).

Fig. 1

The distributions of within-family fluctuations in genus richness and across family volatilities.

(A) Richness fluctuation trajectories of three exemplar families. (B) Empirical cumulative density function of richness fluctuations aggregated across all families, highlighting the exemplar families. To display all families simultaneously, we collapse their fluctuation distributions by dividing by their SDs. If families conform to the Gaussian hypothesis, then their scaled fluctuations should fall along the cumulative density line of a normal N(0, 1) distribution, as shown in (B). We further confirm this normal distribution in fig. S3. (C) The distribution of inverse variances β across all families matches very closely to a gamma distribution (black line); exemplar families are again highlighted.

To predict the superstatistical behavior of the entire marine invertebrate Phanerozoic fauna, we must integrate over all possible local equilibria that each clade could experience. The stationary distribution of β values describes these possible equilibria, specifying the probability that a given clade, chosen at random, will occupy a region of adaptive space characterized by β.

We estimate the distribution of β values as the maximum likelihood distribution describing the set of volatilities for all families, orders, classes, or phyla. Phanerozoic marine invertebrate families follow a gamma distribution in their β values (Fig. 1). The gamma distribution also holds for orders but shows increasing deviations again for classes and especially for phyla (fig. S4).

Using the observation of within-family statistical equilibrium and gamma-distributed β parameters, we can calculate, without further adjusting free parameters, the distributions of family-level fluctuations for the entire marine Phanerozoic, P(x), aswhere is the distribution of fluctuations within a family and is the stationary distribution of volatilities in richness fluctuations. The integral in Eq. 1 leads to

This corresponds to a non-Gaussian, fat-tailed prediction for P(x), which closely matches aggregated family-level fluctuations in the bias-corrected PBDB (Fig. 2).

Fig. 2

Distribution of fluctuations in genus richness within different taxonomic groupings of marine invertebrates in the PBDB () after sampling correction.

(A to D) The distribution is fat-tailed as compared to the maximum likelihood estimate of the normal distribution (blue line). At the family and order levels, the empirical distribution of fluctuations are well described by our superstatistical approach, both when computed from integrating over the distribution of observed variances (red line) and when fit via maximum likelihood [95% confidence interval (CI); red shading in (A) and (B)].

To quantitatively evaluate how well the superstatistical prediction matches the family- and order-level data, we constructed a 95% confidence envelope from bootstrapped maximum likelihood estimates of P(x). Observed fluctuations for both taxonomic levels fall within these 95% confidence envelopes (Fig. 2), indicating that the data do not reject the superstatistical prediction. For further comparison, we fit a Gaussian distribution to the observed fluctuations, which corresponds to the equilibrium hypothesis that all families conform to the same dynamic. Using Akaike Information Criterion (AIC), we find that observed fluctuations are considerably better explained by the superstatistical prediction than by the Gaussian hypothesis (ΔAIC = 1895.622). Thus, as expected under the superstatistical hypothesis, the fat-tailed distribution of fluctuations arises from the superposition of independent Gaussian statistics of fluctuations within families. Computing the distribution of aggregated fluctuations using orders also closely matches the observed data (Fig. 2), but as we further coarsen the taxonomy to classes and phyla, we see increasingly poorer correspondence between data and theory (Fig. 2).

We quantify this change in the goodness of fit with the Kolmogorov-Smirnov statistic (Fig. 3). We can see that both families and orders have low Kolmogorov-Smirnov statistics, and order-level designation of equilibrial subsystems performs slightly better than the family level. Classes are substantially worse, and phyla are worse yet, with the Kolmogorov-Smirnov statistic of phyla being no different than the null randomized taxonomies described below.

Fig. 3

Distribution of Kolmogorov-Smirnov statistics from randomly permuting genera within families (gray shading represents 95% confidence interval).

Solid colored lines are observed Kolmogorov-Smirnov statistics at different taxonomic levels, as indicated.

However, if superstatistical theory explains the data, then this worsening fit with increasing taxonomic scale is expected, as the different classes and phyla should not represent dynamically equilibrial subsystems in their fluctuation dynamics. Instead, classes and phyla aggregate increasingly disparate groups of organisms and thus effectively mix their associated Gaussian fluctuations, meaning that one statistic should no longer be sufficient to describe class- and phylum-level dynamics. We see this confirmed by the increasing frequency of outlier fluctuations in within-class and phylum-level fluctuation distributions (fig. S3). We can also see that families and orders represent, on average, one to two ecospace hypercubes [defined by taxon environment, motility, life habit, vision, diet, reproduction, and ontogeny (, , )], respectively. In contrast, classes and phyla represent, on average, 8 to 30 hypercubes, respectively (fig. S5).

Our analysis indicates that orders and families are evolutionarily coherent units with all subsumed taxa sharing key ecological and evolutionary attributes, allowing them to reach steady-state diversification independently from other clades at the global scale. The fact that both orders and families conform to theoretical predictions is consistent with superstatistics. If superstatistics operates at the order level, then the families subsumed by these orders should represent random realizations of their order’s stationary volatility. The sum of gamma random variables is still gamma, but with new parameters, and so the family-level distribution of is still gamma.

To further test the evolutionary coherence of families, we conducted a permutation experiment in which genera were randomly reassigned to families while maintaining the number of genera in each family. For each permutation, we calculated the superstatistical prediction and its Kolmogorov-Smirnov statistic. The permutation simulates a null model in which common evolutionary history is stripped away (genera are placed in random families), but the total number of observed genera per family is held constant. Because we ignore all instances of no change (i.e., zero fluctuation), we remove any possible large and artificial gaps in the genus occurrences of these permuted clades. Controlling for the total number of genera per family is key because this could be purely an artifact of an arbitrary taxonomic process (, ) and genus richness alone could be solely responsible for differences in β across clades.

We test the possibility that richness is responsible for variation in β in two ways. First, we find that the distribution of genus richnesses within families is not itself distributed gamma (fig. S7), indicating that there is not a simple equivalence between β and the richness of family k. Second, we find that the number of genera in a family and that family’s β value are negatively correlated (fig. S6). A negative correlation between clade richness and β is not unexpected because fluctuations are the sums of the random variables representing genus origination and extinction events; the more of these random variables in the summation (i.e., the more genus richness in a clade), the higher the variance of the summation. Because , increasing richness should lead to decreasing β values. Thus, we want to know if this correlation accounts for all downstream superstatistical results. The permutation test is specifically designed to determine if the β values imposed by this correlation with richness are sufficient to explain the observed superstatistical fit.

Repeating the null permutation of genera in families 500 times yields a null distribution of Kolmogorov-Smirnov statistics that is far separated from the observed values at the family and order levels (Fig. 3), suggesting that the good fit at these levels is not merely a statistical artifact of classification or the richness of clades but carries important biological information. Classes approach the null and phyla are no different. It should also be noted that the width of 95% confidence interval of this null distribution is not far from the distance between the Kolmogorov-Smirnov statistics of orders versus families, suggesting that differences of fit between these taxonomic levels are at least partially accounted for by the randomness of the sampling distribution of Kolmogorov-Smirnov statistics.

DISCUSSION

Our analysis makes no assumption that orders and families should correspond to superstatistical subsystems, but identifies them as the appropriate level for marine invertebrates. Our study is the first to demonstrate that complex patterns in the fluctuation of taxon richness in the fossil record are the result of a simple underlying process analogous to the statistical mechanisms by which complexity emerges in large, nonequilibrium physical () and social systems (). We do so by identifying the biological scale at which clades conform to locally independent dynamic equilibria in fluctuations. Equilibrium could result from many processes, including neutrality (, ), diversity dependence (, ), and processes that dampen, rather than exacerbate, fluctuations in complex ecological networks (). These candidate processes are directly opposed to the presumption of instability underlying the self-organized criticality hypothesis for paleobiodiversity (, ).

We show that the distribution describing the evolution to different equilibria between orders and families is gamma (Fig. 1). A gamma distribution, while consistent with multiple processes, could result from evolution of diversification rates across an adaptive landscape that promotes niche conservatism and pulsed exploration of niche space (). Specifically, if β values are associated with a clade’s macroevolutionarily relevant traits, and those traits evolve via Ornstein-Uhlenbeck–like exploration of an adaptive landscape, then the resulting stationary distribution of β will be gamma (). For macroevolutionary rates to vary in a way consistent with the observed superstatistical description of fluctuations, this landscape cannot be flat (i.e., equal fitness everywhere) but instead must be rugged. Thus, niche conservatism around local fitness optima in adaptive space interrupted by adaptive exploration is likely (, ). The specifics of how this adaptive landscape is shaped and is traversed by evolving clades determine the exact form of the distribution of β volatilities, in the case of the marine Phanerozoic resulting in a gamma distribution. Our work thus motivates further study of the trait spaces and evolutionary shifts consistent with gamma-distributed equilibria in richness fluctuation volatilities.

We show that the pulsed shift to different equilibria between orders and the families they subsume is sufficient to explain the characteristically fat-tailed distribution of richness fluctuations when the marine Phanerozoic invertebrate fauna is viewed as a whole macrosystem. Armed with an understanding of the statistical origin of this diversification pattern, we can explore which models of niche conservatism and pulsed adaptive radiation are consistent with the statistical behavior of the Phanerozoic. Our statistical theory provides new motivation for identifying the eco-evolutionary causes of innovations between lineages and how those innovations are eventually conserved within lineages. Using the superstatistical prediction as a theoretical baseline, we can also go on to identify and robustly examine the mechanisms underlying deviations from statistical theory. For example, some clades wax and wane systematically, and possibly nonsymmetrically, through time (, , ), a pattern that we cannot explain with superstatistics alone.

Superstatistics could also be applied to other areas of evolution and macroecology. For example, recently proposed phylogenetic models already consider heterogeneous rates of diversification [e.g., ()], as expected between different subsystems. The superstatistics of clades in adaptive landscapes could motivate models that jointly predict changes in traits and diversification, a research area currently struggling with model inadequacy (). This framework could also provide a previously unexplored paradigm in modeling the distributions of richness, abundance, and resource use in non-neutral communities, which can be viewed as emerging from the combination of locally equilibrium subsystems. Non-neutral models in ecology are criticized for their overparameterization (), yet a persistent counterargument to neutral theory () is the unrealistic assumption of ecological equivalency and poor prediction of real dynamics (). If ecosystems are viewed as the superposition of many individualistically evolving clades, each exploiting the environment differently and thus obeying a different set of statistics, then diversity dynamics could be parsimoniously predicted with superstatistics while incorporating real biological information on ecological differences between taxa.

Superstatistics is a powerful tool to derive macroscale predictions from locally fluctuating subsystems whose evolution is driven by interesting, but complex and difficult to model, biological mechanisms. Hence, applications of superstatistics to a wide variety of patterns in ecological and evolutionary systems are ripe for exploration.

MATERIALS AND METHODS

All data processing and analyses were performed in R (), and all codes needed to reproduce our study are provided, with added explanation, in Supplementary Appendix A.

Paleobiology database data download and filtering

Data on individual fossil occurrences and the ecospace characteristics of Phanerozoic marine invertebrates were downloaded from PBDB (https://paleobiodb.org) on 16 November 2018 via the database’s application programming interface (API) (data retrieval and processing script available in the Supplementary Materials). Collections were filtered using the same approach as Alroy et al. () to ensure that only well-preserved marine invertebrate occurrences were used in subsequent analyses. This filtering resulted in 815,222 unique genus-level occurrences. These were further filtered to exclude those occurrences without family-level taxonomy and those collections with age estimate resolutions outside the 11-Ma time bins proposed by Alroy et al. (), resulting in 454,033 occurrences. Time bins were compiled from http://fossilworks.org with a custom script reproduced in the Supplementary Materials. The first and last of these time bins, corresponding to the earliest Cambrian and the latest Cenozoic, were excluded from analysis because their sampling completeness (see below) could not be assessed.

Correcting for imperfect and potentially biased sampling

We developed and used a new and flexible method to correct for known sampling incompleteness and biases in publication-based specimen databases (, ). Incompleteness is inherent in all biodiversity samples, the fossil record being no exception (–). In addition to incompleteness, bias may result from preferential publication of novel taxa (), which exacerbates the difference between poorly sampled and well-sampled time periods. We therefore developed a simple two-step method: We first corrected genus richness for incomplete sampling using the “three-timer” correction () and then further corrected this three-timer richness estimate by accounting for any correlation between the number of genera and the number of publications in a time period.

The three-timer correction estimates the probability of failure to observe a genus in a given time period p as the number of times any genus is recorded before and after that period but not during, divided by the number of genera whose occurrence histories span the period t. To calculate the sampling-corrected richness of a clade k in the time period in question, the observed genera within that clade and time period were divided by 1 − p and their occurrences were summedwhere j ∈ k designates genera in clade k and I is an indicator equal to 1 if genus j occurs in time period t.

is the maximum likelihood estimator of richness in a simple occupancy through time type model assuming binomial sampling () and, in that way, mimics other proposed methods for the fossil record (, ). We avoid parametrically modeling the sampling process through time by instead taking a sliding window of time bins from the Cambrian to the Cenozoic. It should be noted that the three-timer correction compares favorably to other similar methods to account for imperfect detection ().

To eliminate further bias due to preferential publication of novel taxa (), we divided the three-timer–corrected number of genera per family per time period by the expected number of genera, given publications in that time period. The expected number was calculated by regressing the log-transformed three-timer–corrected number of genera on log-transformed number of publications. There is only a weak trend toward higher richness with more publications (fig. S1), meaning that the most important correction comes from the three-timer correction.

Our new method rescales each genus occurrence from 0 or 1 (absent or present) to a weighted number continuously ranging between 0 and 1. Because these weighted numbers represent sampling and bias-corrected occurrences, we can add them arbitrarily, corresponding to the membership of any given genus in any given higher taxonomic group. We must, however, choose a taxonomic level at which to evaluate the relationship between richness and publications; we chose the level of family because this is the most finely resolved option.

We opted not to use subsampling methods (, , ) because these approaches would not be advisable for clades with few genera. However, our new method achieves similar results to subsampling procedures at the global scale across all clades. We directly compared our predicted time series of global fluctuations in genus richness, with results derived from rarefaction and shareholder quorum subsampling (SQS) in fig. S2. Our method shows very minor differences with these subsampling-based predictions, and any discrepancies do not affect the statistical distribution of fluctuations (fig. S2).

Superstatistical methods

We first derived the superstatistical distribution P(x) by fitting Gaussian distributions to clade-level distributions of fluctuations p(x), extracting the inverse variances β of those p(x), testing the best function to describe the distribution of β, and then integrating P(x) = ∫βp(x∣β)f(β). This process allows no free parameters to hone the fit of P(x) to the data. However, each inverse variance must, of course, be estimated for each clade, making its good fit to data all the more surprising. To do so, we used least squares instead of maximum likelihood because the asymmetric fluctuation distributions of small clades were more reliably fit with curve fitting than with maximum likelihood.

We also estimated P(x) directly from the family-level data using maximum likelihood to compare the fit of our superstatistical prediction and that of a simple Gaussian distribution using AIC. To calculate a likelihood-based confidence interval on our prediction, we bootstrapped the data, subsampled fluctuations with replacement from all families, and fit the superstatistics using maximum likelihood to the aggregated fluctuation distribution of each bootstrap replicate.