## No precedent

Human activities are creating a mass extinction event. The intensity of this event is unprecedented during human times, but there have been several comparable events during Earth's history. Roopnarine and Angielczyk examined one of the largest, the Permian-Triassic Extinction (see the Perspective by Marshall). The structure and diversity of communities were key predictors of stability through the event. Furthermore, extinctions were not random, with smaller-bodied species being more prone to extinction. This pattern is in direct contradiction to the patterns seen in our current extinction. Thus, the current anthropogenically driven extinction is fundamentally different from previous catastrophic extinctions.

## Abstract

The fossil record contains exemplars of extreme biodiversity crises. Here, we examined the stability of terrestrial paleocommunities from South Africa during Earth's most severe mass extinction, the Permian-Triassic. We show that stability depended critically on functional diversity and patterns of guild interaction, regardless of species richness. Paleocommunities exhibited less transient instability—relative to model communities with alternative community organization—and significantly greater probabilities of being locally stable during the mass extinction. Functional patterns that have evolved during an ecosystem's history support significantly more stable communities than hypothetical alternatives.

Accelerating rates of extinction will have negative impacts on Earth’s ecosystems and human well-being (*1*), but they have no modern precedents, and ecosystem responses are uncertain (*2*). The Earth's history, however, has recorded episodes of extreme biodiversity loss. Here, we use the most severe of those events, the Permian-Triassic mass extinction (PTME), to examine how the short-term stability of ecological communities may vary as their functional structures are altered by extinction. We modeled paleocommunity stability as local stability (*3*), and transient dynamics after perturbation (*4*). Both measures describe paleocommunity dynamics when species populations are perturbed. Paleocommunity structures were determined by species richnesses, species partitioning among functional groups or trophic guilds, and the distribution of interspecific trophic interactions (Fig. 1A). We examined the relationships between these parameters and stability by comparing paleocommunity food webs to hypothetical alternatives, wherein species richness, the number of trophic interactions, and average interaction strength were held constant as we varied patterns of functional partitioning and trophic interaction distributions (Fig. 1, B to F). We thus characterized paleocommunity stability before, during, and after the mass extinction.

We assembled paleocommunity data for the Middle Permian to Middle Triassic terrestrial ecosystems of the Beaufort Group in the Karoo Basin of South Africa (*5*) (Fig. 2). Only four tetrapod genera of over 50 are known to have crossed the Permian-Triassic boundary (PTB) (*6*–*8*), and richness remained low in the Early Triassic (*9*). Our paleocommunities include the Permian *Pristerognathus*, *Tropidostoma*, *Cistecephalus*, and *Dicynodon* assemblage zones, and the Triassic *Lystrosaurus* and *Cynognathus* assemblage zones (table S1). We used a two-pulsed extinction (*10*) to reconstruct three communities from the *Dicynodon* assemblage zone (DAZ): (i) a Phase 0 (Ph0) community comprising all taxa that occurred within the DAZ; (ii) a Phase 1 (Ph1) survivor community that spanned the lower part of the extinction interval; and (iii) a Phase 2 (Ph2) survivor community comprising taxa that became extinct just before the PTB or survived into the Early Triassic. We refer to Ph1 and Ph2 throughout as mass extinction intervals, and all other paleocommunities, including Ph0, as intervals of background extinction.

Species were partitioned among guilds according to trophic ecology and body size (*11*) (Fig. 1A). Amniotes were divided into ten guilds, with herbivores feeding on a single producer guild, while faunivores preyed on amniote herbivores and faunivores up to two size classes larger and smaller than themselves, with the smallest two faunivore guilds also preying on arthropods (tables S2 and S3). Temnospondyl guilds preyed on temnospondyl and amniote guilds up to two size classes larger and smaller than themselves, as well as fish, insects, and aquatic invertebrates. Fish fed on aquatic producers, insects, aquatic invertebrates, and temnospondyls. Insects were divided into herbivorous, omnivorous, and predatory guilds (*5*). Few noninsect invertebrates are known from the Beaufort Group, but our data include molluscs, myriapods, and conchostracans. Plant species were not treated individually because taxonomic resolution is currently insufficient. We instead assume that producer species behaved neutrally with respect to competition and herbivory, and so we aggregated them into four guilds: aquatic microphytes, aquatic macrophytes, terrestrial production accessible to insects and amniotes, and terrestrial production available to insects only (arboreal amniotes were absent from the Beaufort Group).

Although paleontological data lack details of biotic interactions that are available for modern species, paleotrophic interactions can be inferred from predatory traces, preserved gut contents (*12*), functional morphological interpretations, body size relationships, and habitat (*13*). These define partitions of paleocommunity species richness (*S*) into a fixed number of trophic guilds, as well as a pattern of functional organization or interactions among guilds. Thus, of all food web topologies possible given *S* (*14*), only a subset would be consistent with paleocommunity data (*15*) (Fig. 1B). We constrained this subset further by assuming a hyperbolic decay distribution for the number of prey species per consumer species (*5*, *16*), with more species being relatively specialized, and the distribution’s heavy tail representing fewer, relatively generalized consumers. Finally, we sampled food webs from the constrained subset, describing species interactions with generalized Lotka-Volterra (LV) functions, and each food web as a *S × S* matrix of interspecific interactions. Time-averaging and incomplete preservation limit measurement of the number of biotic interactions and their strengths, so we parameterized interaction number with stochastic draws of the number of prey per consumer from the hyperbolic distribution, as well as corresponding interaction coefficients from uniform distributions (*5*). Positive off-diagonal elements represent predation, and negative elements represent the impacts of predatory species. Diagonal terms describe population growth in the absence of interspecific interactions and are always negative in our models. Each paleocommunity's matrices are therefore of statistically similar connectance and average interaction strength.

We formulated a series of alternative food web models, each one relaxing a constraint on topology, to understand how changing functional organization would have affected paleocommunity stability during the PTME (Fig. 1, C to F). Each model expands the subset of food webs to include webs that are inconsistent with some aspect of the paleocommunity. For example, the first model retains species richness, as well as guild structure, but randomizes species assignments to consumer guilds (Fig. 1C). Subsequent models randomize interactions and species assignments to consumer guilds (holding guild number constant) (Fig. 1D); or remove guild structure altogether but retain the interaction distribution (Fig. 1E); or assign trophic interactions randomly, yielding short-tailed Poisson distributions (Fig. 1F) (*17*).

Paleocommunity and model stabilities were assessed as the proportion of stochastically generated webs that are locally stable and the transient instability of each web. Local stability is a mathematically simplified but tractable measure of stability given typical community complexity, with a community in equilibrium returning to equilibrium after a minor perturbation. The simplification restricts local stability analysis to minor perturbations but is, nonetheless, a useful first approximation (*3*). We measured local stability as the real (noncomplex) part of the dominant eigenvalue of a LV matrix. Negative eigenvalues indicate that any perturbation will decline asymptotically to zero, whereas webs with positive eigenvalues do not return to equilibrium. At least 100 LV matrices were generated stochastically for each model and paleocommunity.

No matrices of the background extinction paleocommunities, or their alternative models, were locally stable. This is consistent with the negative relationship between the probability of local stability, and species richness, density of interspecific interactions, and average interaction strength (*17*). *S* ranged from 71 to 140 for those paleocommunities. The fractions of locally stable webs do increase for Ph1 and Ph2, to 25% and 84%, respectively, as *S* declined to 43 and 19, respectively. The increase is not a function of declining richness but, instead, of functional organization, because, whereas holding *S* and functional organization constant but partitioning species randomly (Fig. 1C) yielded frequencies of 56% and 88% in Ph1 and Ph2, respectively, frequencies were significantly lower when functional organization was also randomized (Fig. 1D); 2% and 24%, respectively [Ph1: χ^{2}(df = 2, *n* = 400) = 74.78, *P* < 0.0001; Ph2: χ^{2}(df = 2, *n* = 400) = 134.86, *P* < 0.0001]. Removing guild compartmentalization (Fig. 1E) also reduces the frequencies significantly, to 0% and 41% for Ph1 and Ph2, respectively [χ^{2}(df = 1, *n* = 300) = 30, *P* < 0.0001; χ^{2}(df = 1, *n* = 300) = 58.33, *P* < 0.0001, respectively]. Poisson-distributed models (Fig. 1F) yielded no locally stable webs. Real patterns of functional organization therefore promoted significantly greater probabilities of yielding locally stable food webs during the PTME, regardless of how species were partitioned.

Multispecies models generally amplify perturbations and return to equilibrium nonmonotonically. Such transient short-lived instability is common in real communities (*18*) and contributes to dynamics when perturbations are frequent, as would be expected during the PTME. We therefore measured three transient responses to perturbation: reactivity, maximum amplification, and the timing of maximum amplification (*4*, *5*). Reactivity is the instantaneous rate of displacement in response to perturbation. The dynamics of the displacement yield the maximum amplification during the return to equilibrium. We measured locally stable matrices only, using a Monte Carlo optimization approach to ensure local stability (*5*) even if a stochastically generated matrix was not initially stable.

Paleocommunities were significantly more reactive than alternative models (e.g., Ph0, *F*_{4,497} = 1119.07, *P* < 0.0001) (table S4), and reactivity decreases from the most structured model to the least (Fig. 1, C to F). Almost all background extinction paleocommunities amplify perturbations to a significantly lesser degree, and earlier in time, than do models with alternative guild structure or hyperbolically distributed interactions, although Poisson-distributed webs are least amplified (e.g., Ph0, Kruskal-Wallis *H* = 97.734, df = 4, *P* = 0.0001) (Fig. 3, A and B). Species responses to perturbation would thus be significantly dampened in paleocommunities relative to alternatively structured communities. The mass extinction Ph1 paleocommunity is no more amplified than alternative models (*H* = 136.385, df = 4, *P* = 0.0001) (Fig. 3, C and D), except the Poisson model, which remains the most stable. During Ph2, paleocommunity perturbations are amplified significantly less than those of the functionally altered models (*H* = 32.613, df = 2, *P* = 0.0001) (Fig. 3F), but there is no longer any distinction based on interaction distributions (*H* = 2.176, df = 2, *P* = 0.1402) (Fig. 3E). Paleocommunities during the PTME were therefore as stable or significantly more stable than models with alternative functional organization. The Early Triassic *Lystrosaurus* assemblage zone (LAZ) paleocommunity, in the aftermath of the PTME, is indistinguishable from the guild-randomized and hyperbolic interaction distribution models (*H* = 3.759, df = 3, *P* = 0.2887) (Fig. 3, F and G) but differs significantly from the Poisson-distributed model (*H* = 53.455, df = 1, *P* = 0.0001) (Fig. 2F). This paleocommunity was therefore no more stable than the randomized models, an observation consistent with LAZ’s inferred hypersensitivity to simulated disruptions of primary productivity (*11*). By the Middle Triassic, though, the *Cynognathus* assemblage zone community exhibited dynamics similar to those of the premass extinction communities (table S4).

The increased frequencies of locally stable paleocommunities and greater transient stability of the Ph1 and Ph2 paleocommunities, compared with alternative models, suggests that the extinctions were not random. We tested this by comparing the Ph1 paleocommunity to simulated communities produced by applying the Ph1 extinction randomly to the Ph0 paleocommunity. The simulated communities differed significantly in composition from the Ph1 community [nonmetric multidimensional scaling (NMDS) Bray-Curtis dissimilarity, *n* = 100, Hotelling’s *T*^{2} = 10.355, *P* = 0.0076] (Fig. 4A). The average probability of producing any single simulated community by random extinction, although exceedingly small (*P* = 1.19 × 10^{−29}, *n* = 100), is significantly larger than the probability of producing the Ph1 paleocommunity itself (*P* = 1.61 × 10^{−34}) (*z* test, *P* < 0.001) (*5*). The bias stems from a significantly elevated frequency of extinction of small-bodied amniotes, in contrast to medium and very large amniotes (Fisher’s exact test, *P* = 0.016) (Fig. 2). Furthermore, although the proportion of locally stable simulated communities does not differ from that of the Ph1 community (19%, Fisher’s exact test, *P* = 0.308), they are significantly more reactive (*F*_{1,199} = 869.9, *P* < 0.0001) and have greater values of maximum amplification (*H* = 10.426, *P* = 0.0012). Biased extinction therefore produced a Ph1 community more stable than expected had extinctions been unbiased. There is no evidence of bias, however, in the Ph2 extinction.

Communities before the PTME therefore exhibited less transient instability than alternative model communities because of their functional organization. During the PTME, communities biased by the selective extinction of small-bodied amniotes during Ph1 exhibited no more, and in some cases less, transient instability than alternative models and were more likely to be locally stable. Community stability and its maintenance were, therefore, defining features of this Permian-Triassic ecosystem.

## Supplementary Materials

## REFERENCES AND NOTES

- ↵
- ↵
- ↵
- ↵
- ↵Materials and methods are available as supplementary materials on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**ACKNOWLEDGMENTS:**Data are available in the supplementary materials. We thank two anonymous reviewers for comments. This work was supported by NSF Earth Life Transitions Grant 1336986 (to P.D.R., K.D.A., and C. Sidor).