Non-technical Summary
The chances that existing taxa such as genera and species are replaced by new taxa decreases over geological time, and older taxa persist for longer and longer. This long-term trend has been interpreted as resulting from repeated culling and non-replacement of groups with high turnover or, in other words, as resulting from long-term selection for extinction resistance. It is therefore assumed that taxon-specific properties are critical for this trend. However, taxa are part of eco-genealogical units, that is, they share a history of ecological interaction and a genealogy. To date, structural aspects, such as the changing composition of eco-genealogical units, have received little attention as a possible explanation.
Here we analyzed the temporal changes within large-scale fossil assemblages of bottom-dwelling organisms (called “Phanerozoic benthic marine mega-assemblages,” or PMAs). PMAs are groups of genera of marine animals delineated through hierarchical modular partitioning in a spatiotemporal multilayer network representation of fossil occurrences. PMAs are interpreted as representing eco-genealogical units. For each higher-level PMA, we calculated how its median genus age, its diversity, and the evolutionary rates of its genera changed over time. The resulting estimates support the assumption that the assemblages are meaningful eco-evolutionary units, because their temporal patterns differ from those of randomly chosen Phanerozoic time intervals.
PMAs differ from randomly chosen intervals in that the increase of genus age is associated with a decline in origination proportions and at the same time with increasing between-site diversities. This contradicts the hypothesis of a simple Phanerozoic sorting toward extinction resistance and a taxon’s ability to tolerate a wide range of environmental conditions. The observed pattern can be explained as a process of eco-evolutionary entrenchment resulting from conditioned assembly and weak within-lineage directionality toward specialization and geographic restriction.
Introduction
One of the major features of the Phanerozoic evolution of ecological communities is the long-term decline in taxon extinction rates (Raup and Sepkoski Reference Raup and Sepkoski1982; van Valen Reference Van Valen1984; Foote Reference Foote2000). This decline is linked to a decrease in origination rates (Van Valen and Maiorana Reference Van Valen and Maiorana1985; Sepkoski Reference Sepkoski1998; Gilinsky Reference Gilinsky1994; Foote Reference Foote2000; Bambach et al. Reference Bambach, Knoll and Wang2004; Alroy Reference Alroy2008), a pattern that has been termed the “extinction/origination ratchet” (Stanley Reference Stanley2007).
There is convincing evidence that extinction and origination rates are taxon specific, that is, they are seen as an intrinsic property of a given taxon (e.g., van Valen Reference Van Valen1973, Reference Van Valen1985; Krug and Jablonski Reference Krug and Jablonski2012), and that they indeed are causally coupled (Stanley Reference Stanley1985, Reference Stanley, Allmon and Ross1990, Reference Stanley2007). This led to the conclusion that the observed rate decline is the result of a sorting of taxa by culling and non-replacement of groups with high turnover (Van Valen Reference Van Valen1987; Gilinsky Reference Gilinsky1994). Accordingly, and in line with these findings, it has been argued that during the Phanerozoic, the capacity to resist the fate of extinction became the most important determinant of taxon diversity (Finnegan et al. Reference Finnegan, Payne and Wang2008; Knope et al. Reference Knope, Bush, Frishkoff, Heim and Payne2020).
A logical consequence of this Phanerozoic pattern of declining rates and increasing diversity is that the biota increasingly consists of long-lived taxa during the Phanerozoic, as fewer taxa go extinct while fewer new ones originate. Long-lived taxa have low evolutionary rates (i.e., they are less susceptible to extinction), and they are characterized by low diversity volatilities, which is probably the reason why they accumulate in the biota over time (Gilinsky Reference Gilinsky1994).
In this perspective, the properties of taxonomic groups are at the center of consideration, that is, their resistance to extinction (Van Valen Reference Van Valen1987; Knope et al. Reference Knope, Bush, Frishkoff, Heim and Payne2020) and diversity volatility (Gilinsky Reference Gilinsky1994) are seen as the decisive factors in the dynamics. However, individuals of species and higher taxonomic groups live and lived in evolving local communities (i.e., “sets of species governed by species interaction, resource acquisition, and microhabitat” [Jarzyna et al. Reference Jarzyna, Ohyama, Economo, Gill, Mascarenhas, Okie and Qin2025: p. 2]). Ecological communities, in turn, are formed by regional species pools—sets of species that can potentially influence the community composition and dynamics (Jarzyna et al. Reference Jarzyna, Ohyama, Economo, Gill, Mascarenhas, Okie and Qin2025). Furthermore, local communities and regional species pools can be viewed as part of the nested eco-genealogical hierarchy of the evolutionary units of life, known as the “Bretskyan hierarchy” (Spiridonov and Eldredge Reference Spiridonov and Eldredge2024). This combines the genealogic hierarchy of Linnean taxa; the community hierarchy of, for example, molecules, cells and organisms; and the ecosystem hierarchy of biotic and abiotic chemical reaction networks (Spiridonov and Eldredge Reference Spiridonov and Eldredge2024). The units through which the Bretskyan hierarchy of life is formed, such as holobionts and geobiomes, are relatively stable in terms of both time and space; they are known as “Bretskyan units” (Spiridonov and Eldredge Reference Spiridonov and Eldredge2024: p. 195).
Members of a species interact with each other and with other species within local communities, which are shaped by these interactions and in turn shape the members of the species. They are also part of broader, evolving entities that influence them and which they also influence. The evolving structure of the communities, ecosystems, and geobiomes themselves, therefore, must be seen as a factor in the Phanerozoic patterns of diversity and evolutionary rates. In other words, resistance to extinction and the volatility of diversity could be seen not as properties of the taxon, but of the ecological community and of the broader, evolving entities to which it belongs (see also Dussault and Bouchard Reference Dussault and Bouchard2017; Roopnarine Reference Roopnarine2025).
The distribution of ages of taxa of ecological communities, ecosystems, and regional faunas and floras is one structural aspect of eco-genealogical systems that has received little attention. A few studies that explicitly considered the evolutionary age of taxa found differences indicating that past environmental changes had influenced the taxonomic age structure of faunas (Boyajian Reference Boyajian1992; Van Der Meulen et al. Reference Van Der Meulen, Peláez‐Campomanes and Levin2005; Jablonski et al. Reference Jablonski, Roy and Valentine2006; Harnik et al. Reference Harnik, Jablonski, Krug and Valentine2010). It can therefore be hypothesized that a community or regional species pool composed of taxa with a long shared evolutionary history will respond differently to external change compared with one composed mainly of evolutionary newcomers.
The concept that perhaps best describes such a shared evolutionary history is the concept of “ecological memory,” which is defined as the “capacity of past states or experiences to influence present or future responses of the community” by Padisák (Reference Padisák1992: p. 225). If one considers species and higher taxa to be units of an ecological community that share certain combinations of developmental, functional, habitual, and morphological traits that emerge and disappear alongside the taxon, then the members of a species or genus can be seen as forming an ecological memory of the period during which they emerged.
The original, relatively narrow neontological meaning of the concept refers to the ability of limnic microorganisms to persist within local refugia over long time intervals. When conceptualized more broadly as the capacity of ecological systems to drive their own conditions based on antecedent conditions (e.g., see definition in Benito et al. Reference Benito, Gil-Romera and Birks2020: p. 1), this meaning can be transformed toward geological time frames. One example from paleobiology is the potential recovery of taxa from nearshore or deep-sea refugia after mass extinctions (Erwin Reference Erwin and Valentin1996, Reference Erwin2001) or more generally after disruptive eco-evolutionary transitions (Zambito et al. Reference Zambito, Brett, Baird, Kolbe and Miller2012). From this perspective, the ecological memory is the capacity of the postextinction system to draw on the pre-extinction interval.
Ecological communities, ecosystems, and geobiomes in paleobiology can be reconstructed via fossil assemblages. In paleobiology, it has become established practice to distinguish between community groups sensu Boucot (Reference Boucot1981: p. 350), that is, recurrent fossil associations of similar species in a particular biogeographic region through geological time, or ecological evolutionary subunits (EESUs), that is, recurrent fossil associations exhibiting relative spatial and temporal stability over geological time (Brett et al. Reference Brett, Ivany, Zambito, Welych-Flanagan and Baird2025). At larger scales, empirical evidence supports the distinction of 12 global ecological evolutionary units (EEUs), which characterize tens of millions of years of relatively stable mega-assemblages of taxa punctuated by short intervals of rapid and drastic reorganization (Boucot Reference Boucot1983; Brett et al. Reference Brett, Ivany, Zambito, Welych-Flanagan and Baird2025). These units, in turn, can be understood as subdivisions of Sepkoski’s (Reference Sepkoski1981) three global scale evolutionary faunas (Sheehan Reference Sheehan1996), which found subsequent support from multivariate statistics (Alroy Reference Alroy2010) and network analysis (Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021; therein termed “Phanerozoic marine benthic mega-assemblages” [PMAs]). Community groups, EESUs, EEUs, Sepkoski’s Evolutionary Faunas, and PMAs can be seen as representing units at specific levels of the Bretskyan hierarchy.
The Phanerozoic increase in taxon longevity (Boyajian Reference Boyajian1992; Gilinsky Reference Gilinsky1994) suggests a general trend toward increased ecological memory within ecological communities and more broadly within the units of the eco-genealogical Bretskyan hierarchy over the past half billion years, with significant setbacks during periods of mass extinction. Long-lived species exhibit combinations of traits that developed earlier in their place of origin or in the past. These traits could therefore increase ecological memory. Here, the PMAs of Rojas et al. (Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021) are used to test whether and how the age distribution and associated paleobiological indicators, such as α-, β-, and γ-diversities and evolutionary rates changes within the units of the eco-genealogical Bretskyan hierarchy.
PMAs are network representations obtained from a multilayer network analysis of fossil occurrences at the genus level and at per-stage and per-paleogeographic grid cell resolution. These are retrieved from the Paleobiology Database (see Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021). A multilayer network analysis is a mathematical framework that models the interdependent connections between entities across different layers. In this case, the layers are chronostratigraphic stages. The multilayer network analysis results in a hierarchical partition into divisions called “modules” that represent spatiotemporally constrained sets of genera. These sets are interpreted as PMAs (Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021). Here we used the three highest levels (PMA levels 1–3) of Rojas et al. (Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021; Fig. 1). The first level contains four modules with PMAs starting at the following boundaries: earliest Cambrian (base Fortunian), Late Cambrian (base Jiangshanian), late Permian (base Changhsingian), and Cretaceous (base Hauterivian), with an average temporal length of 161 Myr. The second level contains 16 modules with an average temporal length of 62 Myr; and the third level contains 62 modules with an average temporal length of 32 Myr, of which 50 are longer than three stages.
Temporal range of benthic marine Phanerozoic mega-assemblages (PMAs). A, First-level PMAs. PMA number (nr.) 1 is equivalent to the Cretaceous–Quaternary (mK–Q) mega-assemblage; PMA nr. 2, to the Paleozoic (Pz) mega-assemblage; PMA nr. 3, to the Triassic–Late Cretaceous (Tr–lK) mega-assemblage; and PMA nr. 4, to the Cambrian (Cm) mega-assemblage of Rojas et al. (Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021). B, Second-level PMAs. C, Third-level PMAs. The y-axes indicate the PMA numbers (nr.) that function as identifiers at the respective levels. In B and C, the numbered PMAs are sorted along the y-axes of the diagrams based on their higher-level PMA affiliation. Assemblages with a range of fewer than three stages are excluded from analysis. Color codes in B and C refer to those of the respective first-level PMAs (A). Dashed red vertical lines are five mass extinction events at Ordovician/Silurian, Devonian/Carboniferous, Permian/Triassic, Triassic/Jurassic, and Cretaceous/Paleogene boundaries. Abbreviations: Cm, Cambrian; O, Ordovician; D, Devonian; C, Carboniferous; P, Permian; Tr, Triassic; J, Jurassic. K, Cretaceous; Pg, Paleogene; Ng, Neogene.
Figure 1. Long description
A three-panel multi-plot displays P M A ranges against a geological timescale in M a. The x-axis at the bottom shows time from 500 M a to 0, with colored blocks for periods: C m, O, S, D, C, P, T r, J, K, P g, N g, and Q. Five dashed red vertical lines mark mass extinction events.
* Panel A (1st level): Four broad horizontal bars. P M A n r 4 is light green in the C m. P M A n r 2 is teal, spanning from O to the end of P. P M A n r 3 is purple, spanning from T r to mid-K. P M A n r 1 is orange, spanning from mid-K to Q.
* Panel B (2nd level): Shows nested sub-assemblages. Light green bars are in C m. Teal bars (n r 6 to 10) span O through P. Purple bars (n r 11 to 14) span T r through J. Orange bars (n r 1 to 5) span K through Q.
* Panel C (3rd level): Shows the highest granularity with numerous short horizontal bars. Teal bars (n r 20 to 37) are concentrated between 480 and 250 M a. Purple bars (n r 38 to 50) are concentrated between 250 and 145 M a. Orange bars (n r 1 to 19) are concentrated from 145 M a to the present.
Vertical grey bands in the background indicate specific geological stages across all three panels.
Using these PMAs as approximations of units in the eco-genealogical Bretskyan hierarchy enables us to gain a better understanding of the impact of Phanerozoic assemblage dynamics on taxon properties such as evolutionary rates. This should also be seen as an approach to overcoming perspectives that, strictly speaking, regard Linnaeus’s genealogical hierarchy (and not the entire spectrum of the Bretskyan hierarchy in the sense of Spiridonov and Eldredge [Reference Spiridonov and Eldredge2024]) as the decisive explanatory basis for the paleobiological trends of the Phanerozoic.
Methods
Multilayer Network Representation
We used a multilayer network analysis, using the Map Equation Framework, based on fossil occurrences retrieved from the Paleobiology Database (see Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021). The input network represents chronostratigraphic stages as ordered layers and includes two sets of physical nodes, taxa (i.e., genus-level occurrences), and spatiotemporal units (i.e., spatially explicit and stage-restricted grid cells created by aggregating fossil occurrences). These sets of nodes are connected through weighted links indicating their connectivity strength, calculated using a common neighbors index (Rojas et al. Reference Rojas, Patarroyo, Mao, Bengtson and Kowalewski2017) that accounts for the collection-based structure of the underlying data (Peters and McClennen Reference Peters and McClennen2016). This multilayer network representation captures both spatial (intra-stage) and temporal (inter-stages) relationships between grid cells and benthic marine genera through the Phanerozoic. No original data, network representation, or aggregation approaches were used in our study. Instead, we used the inputs and procedures described in Rojas et al. (Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021) to obtain a multilevel network partition of the Phanerozoic marine biodiversity. We chose this approach because it provides a robust framework for reconstructing taxon age distributions across different hierarchical levels in the organization of marine biodiversity. We used the genus-level taxonomic resolution because it is generally seen as a meaningful proxy for species and provides a good compromise between data availability (i.e., occurrence data in the Paleobiology Database), spatial grain resolution, the underlying species-level evolutionary dynamics (e.g., see Hendricks et al. [Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014], Sigwart et al. [Reference Sigwart, Sutton and Bennett2018], and Graham et al. [Reference Graham, Storch and Machac2018] for an in-depth discussion on what a genus is and the justification and limitations on their uses in analyses like the one presented here).
The resulting hierarchical modular solution reveals PMAs at different spatiotemporal scales (Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021; Eriksson et al. Reference Eriksson, Carletti, Lambiotte, Rojas, Rosvall, Battiston and Petri2022). At the highest scale, four mega-assemblages can be distinguished: (1) Cambrian (Cm), (2) Paleozoic (Pz), (3) Triassic–Late Cretaceous (Tr–lK), and (4), mid-Cretaceous–Quaternary (mK–Q) assemblages following Rojas et al. (Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021). Only the three highest levels of the network partition are considered herein, because their modules are taxonomically overlapping. The upper two PMA levels are global in scale. In PMA level 3, the mega-assemblages form regionally distinct, geographically coherent units (see Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021: fig. 6). Although the stage-restricted grid cells can belong only to one module at each hierarchical level in the network partition, modules (= PMAs) exhibit chronostratigraphic overlap when grid cells from the same layer (or geological stage) are clustered into different modules (Fig. 1). In contrast, the same genera can occur in multiple PMAs at each hierarchical level. This difference between taxa and grid cells arises because state nodes representing the same physical node—and carrying information about geological stages—can be assigned to different modules. Each grid cell has only one state node, whereas a given taxon can have multiple state nodes, one for each geological stage in which it occurs (see Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021: fig. 1). The PMAs have on average a per-stage size of 44 cells at the first level, 31 cells at the second level, and 23 cells at the third level. The PMA sizes change through time, but the sizes are generally not (positively or negatively) associated with the PMA duration (Supplementary Table 3). The PMAs of the first level contain on average 3874 genera; those of the second level, 939 genera; and those of the third level, 292 genera. The percentage of genera assigned to PMAs from the total pool of genera occurring in the respective stage is on average 70% for the PMAs of the first level, 55% for the second level, and 38% for the third level.
Additionally, 100 randomized time intervals were created for each PMA level. These time intervals are positioned randomly at the start of stage boundaries within the Phanerozoic Eon. Their temporal length is random, using the frequency distributions of the time lengths of the PMAs at each level as a model. The cohort of genera occurring in these intervals are treated as random control assemblage. However, because the mega-assemblages at PMA level 3 are not global in scale, the randomized samples are less comparable with the PMAs at level 3.
The absolute time from the beginning of a PMA to its end (its geological duration) is measured in millions of years. To enable comparison of PMAs of different durations, their durations were scaled to one. The scaled time measures the time that has elapsed since the start of the respective PMA interval (time 0) relative to the total duration (time 1) of that PMA.
Paleobiological Estimates
For each PMA and randomized time interval longer than at least three stages (three stages is our deliberately chosen minimum time-series length to estimate temporal trends), the per-stage median genus ages, the richness estimates, turnover rates, and β-diversities were calculated from the cohort of genera occurring in the respective interval (the data and code to produce the estimates are available at https://doi.org/10.5281/zenodo.19564332.)
The per-stage median genus age is the median age of the cohort of genera of the PMA or a random interval occurring in one stage, measured as the time length from its global first occurrence to the upper boundary of the respective stage. The genus richness (γ3timers-diversity) is estimated as “corrected sampled-in-bin diversity” of Kocsis et al. (Reference Kocsis, Reddin, Alroy and Kiessling2019), which is based on “three-timers,” genera that occur in three consecutive time intervals (Alroy Reference Alroy2008). The turnover rates are based on the originations and extinctions of genera; therefore, a genus may have originated before its first occurrence in a PMA and/or it may have occurred after its last occurrence in a PMA. The originations and origination rates are calculated as “second-for-third” corrected rates of Alroy (Reference Alroy2015) using the R package divDyn v. 0.8.2 (Kocsis et al. Reference Kocsis, Reddin, Alroy and Kiessling2019). Because this application of the three-timers rule in the calculations can come at the expense of boundary effects at the beginning and end of the PMA durations, we additionally calculate the γ-diversity (herein: γ3D) with the per-bin rarefaction–extrapolation method of the iNEXT framework developed by Chao et al. (Reference Chao, Chiu and Jost2014) using the R package iNext.beta3D v. 1.0.2 (Chao and Hu Reference Chao and Hu2023; Chao et al. Reference Chao, Thorn, Chiu, Moyes, Hu, Chazdon and Wu2023) with a fixed coverage level of 0.6.
Three different β-diversity measures are calculated to evaluate effects of incomplete sampling and changes in γ-diversity. These three measures were selected because they each correct for these effects in fundamentally different ways. All three measures are based on per-stage genus occurrences in grid cells. Before calculation, cells with fewer than five genera were filtered out. Additionally, all stages with fewer than 25 occurrences, 5 genera, and 5 cells were culled. With this procedure, the β-diversity estimation focuses on the dominance fraction of the sampled taxa and reduces the impact of poor sampling. The three β-diversity measures are β3D, β tot, and βCNESSa.
The β3D is based on the iNEXT.beta3D standardization of Chao et al. (Reference Chao, Thorn, Chiu, Moyes, Hu, Chazdon and Wu2023) under the Hill number framework. Hill numbers are an effective and intuitive way to quantify abundance-based taxon diversities (for a review, see Chao et al. Reference Chao, Chiu and Jost2014), and they allow for a unified partitioning scheme into α-, β-, and γ-diversity components (see Jost Reference Jost2007). The iNext.beta3D standardization approach employs rarefaction and extrapolation methods to provide β-diversity values that statistically remove their dependence from α and γ and thus “purely reflects among-assemblage differentiation (Chao et al. Reference Chao, Thorn, Chiu, Moyes, Hu, Chazdon and Wu2023: p. 6). Here, using the R package iNext.beta3D v. 1.0.2 (Chao and Hu Reference Chao and Hu2023), the β-diversity is calculated as a coverage based (coverage level = 0.6), rarified, and extrapolated Hill number under the parameter q = 0, which is equivalent to the per-stage generic change among cells.
The β tot is calculated based on the framework suggested in Cardoso et al. (Reference Cardoso, Rigal, Carvalho, Fortelius, Borges, Podani and Schmera2014); it reflects the species replacement and the spatial loss or gain between pairs of sites, with dissimilarity fractions here derived from the Sørensen index, using function beta() of the R package BAT v. 2.9.6 (Cardoso et al. Reference Cardoso, Rigal and Carvalho2015), here resulting as an arithmetic mean from 100 rarefactions by the minimum abundance among all sites. The β tot is based on presence–absence data.
The βCNESSa is based on randomly drawn individuals from pairs of samples (the “CN” in the measure’s name stands for “cord-normalized”; Trueblood et al. Reference Trueblood, Gallagher and Gould1994), assuming a hypergeometric abundance distribution in each sample; it measures the Horn-Morisita index based on the expected shared species (ESS) and is here calculated using the equation given in Zou et al. (Reference Zou, Zhao, Wu, Lai, Peres‐Neto and Axmacher2025: p. 3) with the R package rarestR v. 1.1.1 (Zou et al. Reference Zou, Zhao, Wu, Lai, Peres‐Neto and Axmacher2025). The respective per-stage estimates of βtot and βCNESSa are the arithmetic means of all pairwise β-diversities.
The methods employed herein for estimating evolutionary rates and diversities comprise different standardization approaches that should minimize sampling and preservation biases. The problems arising from heterogenous sampling coverage over geological time are complex and include various mutual dependencies (e.g., see Peters and Foote Reference Peters and Foote2001; Na et al. Reference Na, Li, Krause, Zhu and Kiessling2023; Antell et al. Reference Antell, Benson and Saupe2024). These problems are not discussed in detail here, because our analyses aim to provide initial insights into potential patterns.
Additionally, the per-stage number of genus occurrences (G obs), the per-stage number of singletons (G s) (i.e., genera with one occurrence per stage), the number (n) of occupied cells (n cells), the number of cells containing only a single genus (n singleton-cells), and the number of one-cell genera (i.e., genera that occur only in one geographic cell) are counted for each PMA. These give information on how the estimates are associated with raw abundances.
We checked for autocorrelation in all variables before calculating correlations between them using the acf function in R. Because all variables showed autocorrelation, we detrended the data by calculating the first differences and used the first differences for Kendall’s rank correlation, except for correlations with scaled time.
Results
Genus Age
The median ages of the genera that occur in a stage are not associated with the durations of the respective stage (τ = −0.073, p = 0.295). The differing lengths of stratigraphic stages, therefore, have no systematic effect on the estimated median ages.
The plot of the median ages of genera of all PMAs that occur in each stage shows a highly volatile, asymptotic growth trend (Fig. 2A). There is a rapid increase of the values during the Paleozoic, from values well below 10 Myr in the Cambrian toward a leveling off at values of around 30 Myr in the Carboniferous to Permian. Notably, the lowest median ages occur in the recovery intervals just after the end-Ordovician (Aeronian Stage: median genus age = 6.7 Myr, first quartile genus age = 4.9 Myr, third quartile genus age = 22.4 Ma), the end-Permian (Olenekian Stage: median genus age = 4.97 Myr, first quartile genus age = 4.1 Myr, third quartile genus age = 32.1 Ma), and the end-Cretaceous (Danian Stage: median genus age = 9 Myr, first quartile genus age = 4.4 Myr, third quartile genus age = 63.4 Ma) mass extinctions. Generally, high median ages occur in intervals just before these mass extinctions. This results in sharp median age declines during mass extinctions and a relatively slow trend of general increase in the times between mass extinction intervals (τ for the respective intervals: Cambrian–Ordovician = 0.78; Silurian–Devonian = 0.28; Carboniferous–Permian = 0.62; Triassic = 0.37; Jurassic–Cretaceous = 0.64; Cretaceous–Pleistocene =0.50; p-values in all cases < 0.05). This pattern roughly resembles the median age curves for families published by Boyajian (Reference Boyajian1992: fig. 2), in that it has a steep Paleozoic increase punctuated by drastic falls during mass extinction intervals.
Per-stage median genus age of: A, the total set of benthic marine invertebrates of the Phanerozoic Eon assigned to Phanerozoic mega-assemblages (PMAs); and B, the cohorts of genera belonging to first-level PMAs. Colors of trajectories are the same as in Fig. 1A. See also Fig. 1.
Figure 2. Long description
Two vertically stacked line graphs share a horizontal axis representing time in M a from 541 to 0, labeled with geological periods: C m, O, S, D, C, P, T r, J, K, P g, and N g. Vertical dashed red lines mark major extinction boundaries.
Panel A: The vertical axis is median age in M y r ranging from 0 to 50. A single black line with circular data points fluctuates significantly. It starts near 10 M y r in the C m, shows a general upward trend with high volatility, peaks near 48 M y r at the end of the P, drops sharply to 5 M y r at the start of the T r, and ends near 45 M y r in the N g.
Panel B: The vertical axis is median age in M y r ranging from 0 to 200. Four distinct colored trajectories represent different P M A cohorts:
- A light green line starts in the C m and ends in the O.
- A teal line begins in the O, showing a steady linear increase from 10 M y r to approximately 160 M y r by the end of the P.
- A purple line begins at the start of the T r at 0 M y r and increases to about 80 M y r by the end of the J.
- An orange line begins in the K at 40 M y r and rises to 55 M y r by the N g.
Both panels feature vertical gray bars of varying widths in the background.
This zig-zag pattern is even more pronounced when plotting the per-stage median age of the cohort of genera of each PMA at PMA level 1 individually (Fig. 2B). This temporal increase of the per-stage median genus-age, however, is not a common feature of PMAs at all three levels (Fig. 3). The increase is statistically significant only at the highest PMA level (Table 1). Hence, at the beginning, the first-level PMA is dominated by young genera, and by the end of the PMAs, the genera are, on average, older. In the randomly selected intervals, such an increase is apparent at all levels, although with decreasing steepness toward lower-level intervals (Fig. 4, Table 2).
Per-stage estimates of median genus age and γ- and β-diversities of Phanerozoic mega-assemblages (PMAs) of three hierarchical levels along their respective scaled durations, where 0 is the beginning and 1 the end of the duration of a PMA. Each dot represents one per-stage estimate of one PMA at the respective level (PMA1st level n = 4; PMA2nd level n = 16; PMA3rd level n = 50). Black lines are the estimates, averaged in 10 time bins. Shaded areas are the 95% confidence intervals of the binned averages. See also Fig. 1. Point color identifies respective PMA level 1 membership (see Fig. 1).
Figure 3. Long description
The grid is organized by columns representing P M A Level 1, P M A Level 2, and P M A Level 3, and rows representing different metrics. All panels share an x-axis of duration scaled from 0.00 to 1.00.
Top Row: Median genus age in M y r.
- P M A Level 1 shows a steady linear increase from approximately 20 to 100 M y r.
- P M A Level 2 shows a fluctuating but generally stable trend around 40 M y r.
- P M A Level 3 shows a stable trend around 35 M y r with high scatter.
Middle Row: Gamma 3timers in n genera.
- P M A Level 1 shows a jagged increase peaking near 600.
- P M A Level 2 shows a hump-shaped curve peaking at 0.75 duration.
- P M A Level 3 shows a relatively flat trend with high density between 300 and 600.
Bottom Row: Beta C N E S S a.
- All three levels show an initial sharp increase followed by an asymptotic plateau near 0.95.
In all panels, individual data points are colored purple, green, and orange. A thick black line represents the binned average, and a light gray shaded area represents the 95 percent confidence interval. The scatter of points is most dense in the P M A Level 3 column.
Results of Kendall’s rank correlation of first differences of Phanerozoic mega-assemblage (PMA) estimates per PMA level (except correlations with scaled time, which are based on original values). Age: median genus age; duration: scaled time; extinctions: extinction proportions; originations: origination proportions. βCNESSa: β-diversity, normalized after Zou et al. (Reference Zou, Zhao, Wu, Lai, Peres‐Neto and Axmacher2025); γ3timer: γ-diversity, estimated as “corrected sampled-in-bin diversity” of Kocsis et al. (Reference Kocsis, Reddin, Alroy and Kiessling2019). All correlations refer to the respective per-stage estimates. Bold: significance level ≤ 0.05. Green shading: positive association; red shading: negative association
Table 1. Long description
The table consists of four columns: the variable pair being correlated, P M A level 1, P M A level 2, and P M A level 3. Data points include tau values and p-values. Significant correlations (p less than or equal to 0.05) are in bold.
* Duration versus age: P M A 1 tau 0.475 (p less than 0.05); P M A 2 tau 0.058 (p less than 0.05); P M A 3 tau 0.070 (p 0.12).
* Duration versus gamma sub 3timer: P M A 1 tau 0.096 (p 0.13); P M A 2 tau 0.288 (p less than 0.05); P M A 3 tau 0.053 (p 0.19).
* Duration versus beta sub C N E S S a: P M A 1 tau 0.282 (p less than 0.05); P M A 2 tau 0.264 (p less than 0.05); P M A 3 tau 0.203 (p less than 0.05).
* Duration versus originations: P M A 1 tau negative 0.205 (p less than 0.05); P M A 2 tau negative 0.248 (p less than 0.05); P M A 3 tau negative 0.278 (p less than 0.05).
* Duration versus extinctions: P M A 1 tau negative 0.079 (p 0.27); P M A 2 tau 0.106 (p less than 0.05); P M A 3 tau 0.231 (p less than 0.05).
* Gamma sub 3timer versus age: P M A 1 tau 0.011 (p 0.88); P M A 2 tau negative 0.242 (p 0.31); P M A 3 tau negative 0.138 (p 0.28).
* Gamma sub 3timer versus beta sub C N E S S a: P M A 1 tau negative 0.093 (p 0.20); P M A 2 tau 0.055 (p 0.83); P M A 3 tau 0.238 (p 0.06).
* Gamma sub 3timer versus originations: P M A 1 tau 0.357 (p less than 0.05); P M A 2 tau 0.242 (p 0.31); P M A 3 tau 0.348 (p 0.01).
* Gamma sub 3timer versus extinctions: P M A 1 tau 0.294 (p less than 0.05); P M A 2 tau 0.333 (p 0.16); P M A 3 tau 0.378 (p 0.01).
* Beta sub C N E S S a versus age: P M A 1 tau 0.085 (p 0.24); P M A 2 tau negative 0.091 (p 0.76); P M A 3 tau 0.043 (p 0.73).
* Beta sub C N E S S a versus originations: P M A 1 tau 0.029 (p 0.72); P M A 2 tau negative 0.030 (p 0.94); P M A 3 tau negative 0.065 (p 0.66).
* Beta sub C N E S S a versus extinctions: P M A 1 tau negative 0.075 (p 0.35); P M A 2 tau 0.333 (p 0.15); P M A 3 tau negative 0.065 (p 0.66).
* Age versus originations: P M A 1 tau 0.046 (p 0.57); P M A 2 tau 0.273 (p 0.28); P M A 3 tau negative 0.188 (p 0.19).
* Age versus extinctions: P M A 1 tau 0.014 (p 0.86); P M A 2 tau 0.151 (p 0.54); P M A 3 tau negative 0.071 (p 0.63).
* Originations versus extinctions: P M A 1 tau 0.244 (p less than 0.05); P M A 2 tau 0.485 (p less than 0.05); P M A 3 tau 0.440 (p less than 0.05).
Per-stage estimates of median genus age and γ- and β-diversities of randomly selected time intervals equivalent to the three hierarchical levels of Phanerozoic mega-assemblages (PMAs). The y-axis shows the position of the estimate within the scaled duration of the respective interval, where 0 is the beginning and 1 the end of the duration of an interval. Gray dots represent individual per-stage measurements. Shaded areas (contours) represent five point densities, with darkest areas representing highest densities. Black lines are the estimates, averaged in 10 time bins. The 95% confidence intervals of the binned averages are too narrow to be visible.
Figure 4. Long description
A grid of nine scatter plots arranged in three rows and three columns. All plots share a common x-axis labeled duration scaled, ranging from 0.00 to 1.00.
Top row: Median genus age in m y.
- Level 1 equivalent: A strong linear increase from 25 to 100 m y.
- Level 2 equivalent: A moderate linear increase from 20 to 50 m y.
- Level 3 equivalent: A shallow, nearly flat trend hovering around 25 to 35 m y.
Middle row: gamma 3timers in n genera.
- Level 1 equivalent: The black line remains stable near 500 until a duration of 0.75, then increases toward 800.
- Level 2 equivalent: A gradual increase from 500 to 800 across the duration.
- Level 3 equivalent: A slight rise and plateau between 500 and 750.
Bottom row: beta C N E S S a.
- All three levels (1, 2, and 3) show nearly identical horizontal trends. The data is highly concentrated in a dense band at the top of the y-axis, with the black average line remaining flat at approximately 0.80.
In all panels, gray dots represent individual measurements, while shaded contours indicate point density, with the darkest regions representing the highest concentration of data points. The black lines represent binned averages.
Results of Kendall’s rank correlation of first differences of estimates from randomly selected time intervals (except correlations with scaled time, which are based on original values). Age: median genus age; duration: scaled time; extinctions: extinction proportions; originations: origination proportions. βCNESSa: β-diversity, normalized after Zou et al. (Reference Zou, Zhao, Wu, Lai, Peres‐Neto and Axmacher2025); γ3timer: γ-diversity, estimated as “corrected sampled-in-bin diversity” of Kocsis et al. (Reference Kocsis, Reddin, Alroy and Kiessling2019). All correlations refer to the respective per-stage estimates. Bold: significance level ≤ 0.05. Green shading: positive association; red shading: negative association. Compare with Table 1
Table 2. Long description
The table consists of 15 rows of variable correlations across three columns: random P M A level 1, random P M A level 2, and random P M A level 3. Each cell contains a tau value and a p-value.
* Duration versus age: tau values are 0.519, 0.365, and 0.182 respectively, all with p much less than 0.05.
* Duration versus gamma sub 3timer: tau values are 0.104, 0.115, and 0.076 respectively, all with p much less than 0.05.
* Duration versus beta sub C N E S S a: tau values are 0.011 (p equals 0.11), negative 0.007 (p equals 0.19), and negative 0.016 (p much less than 0.05).
* Duration versus originations: tau values are negative 0.027, negative 0.036, and negative 0.023 respectively, all with p much less than 0.05.
* Duration versus extinctions: tau values are negative 0.012 (p equals 0.09), 0.015 (p much less than 0.05), and 0.039 (p much less than 0.05).
* Gamma sub 3timer versus age: tau values are negative 0.046, negative 0.080, and negative 0.148 respectively, all with p much less than 0.05.
* Gamma sub 3timer versus beta sub C N E S S a: tau values are 0.210, 0.219, and 0.218 respectively, all with p much less than 0.05.
* Gamma sub 3timer versus originations: tau values are 0.317, 0.321, and 0.318 respectively, all with p much less than 0.05.
* Gamma sub 3timer versus extinctions: tau values are 0.227, 0.273, and 0.283 respectively, all with p much less than 0.05.
* Beta sub C N E S S a versus age: tau values are 0.007 (p equals 0.33), 0.009 (p equals 0.22), and 0.013 (p much less than 0.05).
* Beta sub C N E S S a versus originations: tau values are 0.211, 0.249, and 0.263 respectively, all with p much less than 0.05.
* Beta sub C N E S S a versus extinctions: tau values are 0.009 (p equals 0.20), 0.064 (p much less than 0.05), and 0.088 (p much less than 0.05).
* Age versus originations: tau values are negative 0.010 (p equals 0.16), negative 0.030 (p much less than 0.05), and negative 0.064 (p much less than 0.05).
* Age versus extinctions: tau values are 0.042, 0.015, and negative 0.008 respectively, all with p much less than 0.05.
* Originations versus extinctions: tau values are 0.196, 0.217, and 0.236 respectively, all with p much less than 0.05.
Differences exist also among the individual PMAs at each level. At the highest level, the Cm-assemblage and the mK–Q assemblage have the lowest relative age increase (15 Myr and 13.72 Myr, respectively). The relative age increase is highest within the Pz assemblage (143.53 Ma). Initially high median ages (Induan Stage: median genus age = 72 Myr, first quartile median genus age = 8.7 Myr, third quartile media genus age = 156.4 Ma) in the Tr–lK assemblage are remarkable and reflect the mass extinction and postextinction interval across the end-Permian mass extinction.
γ-diversity
The per-stage genus richness values of the PMAs are estimated with two independent standardization methods as γ3timer-diversity and as γ3D-diversity. Both estimates resulted in roughly similar age trends (Fig. 3, Supplementary Fig. 1); the γ-diversity is relatively low at the start of the PMAs, leveling off quickly within the first half of their duration. At the first and second PMA levels, strong differences between the individual PMAs can be detected with a general unimodal hat-shaped pattern among the PMAs of the Cm, Pz, and Tr–lK assemblages, and with a general γ-diversity increase in the PMAs of the mK–Q assemblage (Supplementary Fig. 2). The combined per–PMA level results give an inconsistent picture with no clear association between PMA duration and γ-diversity (Table 1).
This pattern strongly differs from the consistent γ-diversity increase among the randomly selected control intervals (Fig. 4, Table 2). This increase is an effect of the well-known overall trend of Phanerozoic diversity increase (e.g., Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fursich, Harries, Hendy and Holland2008).
The γ3timer-diversity of the PMAs is strongly positively associated with the respective raw genus occurrences (τ = 0.661, p << 0.05), with the number of cells occupied by the PMA in the respective interval (τ = 0.546, p << 0.05), and it is weakly negatively associated with the number of singleton cells (i.e., cells occupied by a single genus) (τ = −0.134, p << 0.05). However, it is not associated with the respective number of singletons (i.e., genera that occur in stage in a PMA only once) (τ = −0.095, p = 0.01).
β-diversity
The β-diversity of PMAs increases with increasing PMA duration up to an asymptote of 1, which corresponds to the upper limit specified by the method. This trend occurs consistently in all three PMA levels (Fig. 3, Table 1) when using βCNESSa as a measure. Consistency checks with two other measures (β3D, β tot; see “Methods”) confirm this general trend for the two upper PMA levels (Supplementary Fig. 3, Supplementary Table 1). Generally, the β-diversity is very low at the beginning of a PMA, leveling off within the first half of its duration, with only a slight increase afterward. In the randomly selected intervals, no such trend is evident (Fig. 4, Table 2). The three β-diversity measures are positively associated with the raw genus occurrences (G obs) (except βCNESSa), the number of occupied cells, and the γ-diversity of the respective PMA stage interval (Supplementary Table 2). Furthermore, a positive association is evident between βCNESSa and the origination proportion, and between βCNESSa and γ3timer in the random samples, but not in the PMAs (Tables 1, 2). The βCNESSa is not associated with the proportion of one-cell-genera (PMA level 1: τ = 0.05, p = 0.49; PMA level 2: τ = −0.36, p = 0.79; PMA level 3: τ = −0.05, p = 0.72). This is consistent with the finding that the proportion of one-cell genera does not change systematically with scaled PMA time (PMA level 1: τ = −0.04, p = 0.51; PMA level 2: τ = −0.02, p = 0.60; PMA level 3: τ = −0.02, p = 0.63). This is important to note, because this excludes decreasing geographic genus range as a potential factor for the changes in β. Likewise, changing PMA size can be ruled out as a decisive factor for the trends in β, because there is no consistent trend in PMA size over time (see Table 1).
Evolutionary Rates
A general difference between PMA estimates and estimates from random samples is also evident in the trends of the origination proportions: they tend to fall during PMA duration, but not so in the randomly selected intervals (compare Figs. 5 and 6, Tables 1 and 2). The sig nal from the extinction proportions is more ambiguous; they tend to increase in the lower two PMA levels, but no trends are evident from PMA level 1. There is also no extinction trend within the randomly sampled intervals (Table 2).
Per-stage estimates of origination and extinction proportions of Phanerozoic mega-assemblages (PMAs) of three hierarchical levels along their respective scaled duration, where 0 is the beginning and 1 the end of the duration of a PMA. Each dot represents one per-stage estimate of one PMA at the respective level. Black lines are the estimates, averaged in 10 time bins. Shaded areas are the 95% confidence interval of the binned averages. See also Fig. 1. Point color identifies respective PMA level 1 membership (see Fig. 1).
Figure 5. Long description
A multi-panel grid of six scatter plots arranged in two rows and three columns.
Row 1: Proportion Originations.
* Panel 1 (P M A Level 1): The Y-axis ranges from 0.0 to 0.8. Data points in purple, green, and orange are scattered widely. A black trend line starts at 0.3, fluctuates with a peak near 0.65 duration, and ends near 0.1. A grey shaded confidence interval follows the line.
* Panel 2 (P M A Level 2): The trend line shows a sharp initial peak at 0.15 duration, followed by a gradual decline toward 0.05 at the end of the duration.
* Panel 3 (P M A Level 3): The trend line remains relatively flat between 0.15 and 0.20 until 0.6 duration, then drops sharply toward zero.
Row 2: Proportion Extinctions.
* Panel 4 (P M A Level 1): The trend line shows high volatility with two distinct peaks at 0.15 and 0.45 duration, ending around 0.1.
* Panel 5 (P M A Level 2): The trend line is more stable, hovering between 0.1 and 0.2 throughout the entire scaled duration.
* Panel 6 (P M A Level 3): The trend line starts very low near zero, rises to a plateau of 0.15 at 0.25 duration, and maintains that level with minor fluctuations until the end.
Across all panels, the X-axis is labeled duration scaled from 0.00 to 1.00. Individual data points are colored by P M A Level 1 membership, showing high variance around the binned averages.
Per-stage estimates of origination and extinction proportions of randomly selected time intervals equivalent to the three hierarchical levels of Phanerozoic mega-assemblages (PMAs). The y-axis shows the position of the estimate within the scaled duration of the respective interval, where 0 is the beginning and 1 the end of the duration of an interval. Gray dots represent individual per-stage measurements. Shaded areas (contours) represent five point densities, with darkest areas representing highest densities. Black lines are the estimates, averaged in 10 time-bins. The 95% confidence intervals of the binned averages are too narrow to be visible.
Figure 6. Long description
The figure contains six panels arranged in a two-by-three grid. The top row is labeled proportion originations on the vertical axis, and the bottom row is labeled proportion extinctions. The horizontal axis for all panels is duration scaled from 0.00 to 1.00. The columns are labeled Level 1 equivalent, Level 2 equivalent, and Level 3 equivalent from left to right.
Each panel features a background of gray dots representing individual per-stage measurements. Overlaid on these dots are shaded contour areas representing five levels of point density, with the darkest shading concentrated near the bottom of the graphs between 0.0 and 0.2 on the vertical axis. A thick black line runs horizontally across each panel representing the average of the data in 10 time-bins.
In the top row (originations):
- Level 1 equivalent: The black line remains relatively stable near 0.2, with a slight dip toward 0.15 at the end of the duration.
- Level 2 equivalent: The black line shows a subtle hump, peaking around 0.2 at the midpoint of the duration.
- Level 3 equivalent: The black line shows more volatility, peaking at 0.23 early in the duration before leveling off near 0.2 and dropping at the very end.
In the bottom row (extinctions):
- Level 1 equivalent: The black line is stable near 0.2, with a slight downward trend after the 0.75 mark.
- Level 2 equivalent: The black line is nearly flat, hovering just below 0.2 across the entire duration.
- Level 3 equivalent: The black line shows a slight upward trend, starting near 0.18 and ending near 0.22.
The origination and extinction proportions are highly positively associated with the γ-diversity in the randomly sampled intervals but not consistently so in the PMAs. However, a positive association of origination and extinction proportions is a consistent feature in all estimations (Tables 1, 2). This consistent association between origination and extinction is mirrored in consistent patterns of association of the two rates estimates with the raw genus occurrences (G obs) and with the relative number of cells per PMA stage (positive association) and with the relative number of singletons (i.e., our measure of rarity) (negative association) (Table 3).
Results of Kendall’s rank correlation of first differences of selected per-stage PMA estimates with the raw genus count (G obs), the relative number of one-cell genera (rarity), the raw cell count (n cells), and the relative number of cells containing only a single genus (n singleton-cells) of the respective stage. γ3timer: γ-diversity, estimated as “corrected sampled-in-bin diversity” of Kocsis et al. (Reference Kocsis, Reddin, Alroy and Kiessling2019). Bold: significance level ≤ 0.05. Green shading: positive association; red shading: negative association
Table 3. Long description
A table with four columns of metrics and four rows of variables. The columns are gamma sub 3timer, Extinctions proportion, Originations proportion, and Genus age. The rows are G sub obs, Rarity, n sub cells, and n sub singleton-cells.
* G sub obs row: tau equals 0.661 p much less than 0.05 for gamma sub 3timer. tau equals 0.317 p much less than 0.05 for Extinctions. tau equals 0.399 p much less than 0.05 for Originations. tau equals minus 0.076 p equals 0.06 for Genus age.
* Rarity row: tau equals minus 0.095 p equals 0.01 for gamma sub 3timer. tau equals minus 0.191 p much less than 0.05 for Extinctions. tau equals minus 0.149 p much less than 0.05 for Originations. tau equals 0.029 p equals 0.49 for Genus age.
* n sub cells row: tau equals 0.546 p much less than 0.05 for gamma sub 3timer. tau equals 0.314 p much less than 0.05 for Extinctions. tau equals 0.432 p much less than 0.05 for Originations. tau equals minus 0.069 p equals 0.09 for Genus age.
* n sub singleton-cells row: tau equals minus 0.134 p much less than 0.05 for gamma sub 3timer. tau equals 0.081 p equals 0.08 for Extinctions. tau equals 0.185 p much less than 0.05 for Originations. tau equals minus 0.041 p equals 0.34 for Genus age.
Positive associations are shaded green and negative associations are shaded red. Significant values are in bold.
Discussion
PMAs are Meaningful Eco-evolutionary Units
The hierarchically nested modules that result from a multilevel clustering of the multilayer network representation of a global dataset of marine invertebrate fossils have been interpreted as global PMAs corresponding to Sepkoski’s (Reference Sepkoski1981) evolutionary faunas and to Boucot’s (Reference Boucot1983) EEUs (Rojas et al. Reference Rojas, Calatayud, Kowalewski, Neuman and Rosvall2021). This interpretation assumes that PMAs are meaningful eco-evolutionary units, that is, that they are built through accumulation and disintegration of taxon co-occurrences, with richness peaks randomly distributed at any time during each module’s (= PMA) duration (Raulo et al. Reference Raulo, Rojas, Kröger, Laaksonen, Orta, Nurmio, Peltoniemi, Lahti and Žliobaitė2023). The averaged richness curves, therefore, should describe a hat-shaped waxing and waning pattern along the scaled time of the completed PMAs (however, cf. Hohmann and Jarochowska Reference Hohmann and Jarochowska2019). This is the case for most assemblages (Supplementary Fig. 2). The lower values near the boundaries are not an artifact of the method (such as the three-timer method), because consistent trends result from γ3timers and γ3D, which are methodologically different.
Importantly, and in contrast to the γ trends of the PMAs, no such waxing and waning pattern results from the randomly selected intervals. Here, a continuous richness increase occurs through interval time. This pattern is scale independent and reflects the general Phanerozoic trend of an increasing richness (e.g., cf. Sepkoski et al. Reference Sepkoski, Bambach, Raup and Valentine1981; Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fursich, Harries, Hendy and Holland2008). The scale independency in the waxing and waning pattern of the PMAs, in turn, probably reflects the fractal symmetry of the nested structure of the global-scale modules.
The increase in β-diversity and the decrease in origination proportions over the respective PMA durations are also scale-independent patterns that do not occur in randomly selected intervals. Decreasing evolutionary rates are well known from total Phanerozoic trends (Foote Reference Foote2000; Alroy Reference Alroy2010; cf. Newman and Eble Reference Newman and Eble1999; Plotnick and Sepkoski Reference Plotnick and Sepkoski2001). The volatility of these published estimates across time intervals is however relatively high and the total Phanerozoic decrease is mostly due to the very high rates during the earliest Paleozoic (cf. Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fursich, Harries, Hendy and Holland2008: fig. 1). Randomizing the time intervals, as was done herein, will therefore almost certainly destroy the known global trend of the Phanerozoic curves. It is therefore remarkable that statistically significant trends occur in PMAs at all three levels.
The increasing extinction rates in the low-level PMAs are somewhat unexpected. It is unlikely that they result from sampling effects, because the potential causes of such sampling effects, such as differences in G obs or in the proportion of singletons, have similar associations with origination and extinction proportions at all three PMA levels, respectively (see Table 3). Therefore, it is assumed that they either reflect a real trend or are the result of a different, unknown artifact resulting, for example, from the method (such as a boundary effect).
The consistently relatively high origination and extinction proportions at the beginnings of the first-level PMAs (Fig. 5), probably reflect the exceptionally high turnover pulses from rebounds after mass extinctions (e.g., cf. Erwin Reference Erwin1998; Krug and Jablonski Reference Krug and Jablonski2012; Hull Reference Hull2015). This difference between scales could be another line of evidence for the hypothesis that mass extinctions represent a macroevolutionary regime that differs from times dominated by “background” extinctions (Jablonski Reference Jablonski1986).
Increase of Genus Age
The increase in median age through time is a consistent feature in the randomly selected time intervals. In the PMAs, in contrast, such an increasing trend occurs consistently only at the highest level (Fig. 3). This feature has long been appreciated from Phanerozoic-scale analyses (Boyajian Reference Boyajian1992; Gilinsky Reference Gilinsky1994). Here, it is apparent that it results from a Phanerozoic zigzag trend, with massive setbacks (rapid decline in median genus ages) during PMA transitions and slow rebounds afterward (Fig. 2). This means, at global scale, new PMAs start with significantly lower median genus ages than they end with.
It is, however, important to note that the age trajectories differ significantly between individual PMAs. Some of the high PMA levels show a pattern that can best be explained as an effect of mass extinctions: the Pz and some of the mK–Q assemblages start with relatively high genus median ages, which then decline rapidly in the early part of the PMA, only to recover and to rise above previous levels toward the end. These early intervals with relatively high median ages can be interpreted here as reflecting the postextinction interval after mass extinctions, with a prevalence of generalized genera (e.g., Erwin Reference Erwin2001; Jablonski Reference Jablonski2005) that is followed by a rebound with many newcomers.
In the lower two PMA levels examined, in which the PMAs have shorter absolute durations (mean = 32 Myr in PMA level 3, compared with a mean of 161 Myr for the first-level PMAs), no statistically relevant age trends can be detected. This is despite increasing median ages being the average trend in the level 2 and level 3 equivalent random intervals, respectively (see Fig. 4, Table 2). Hence, in the randomizations, the massive declines that occur within short intervals (at mass extinction events) can be regarded as isolated cases that have no significant influence on the overall trend, which is the globally slowly increasing median age trend. The fact that this global trend is not reflected in individual low-level PMAs, where the median age may rise, fall, or remain constant, indicates that at lower-level eco-genealogical units, either different dynamics are at play than at the global level and at the longest timescales or that at this level (genus, paleogeographic cell) the dynamics are not resolved (see discussion in “Assembly of Global-Scale Eco-genealogical Units”).
Contradictory Model Predictions
The Phanerozoic increase of genus age is linked to a general decline in evolutionary rates (e.g., Gilinsky Reference Gilinsky1994; Finnegan et al. Reference Finnegan, Payne and Wang2008). This is consistent with the neutral theory of biodiversity (NT). NT predicts that in a stochastically evolving community of fixed size and under zero-sum assumption, early-appearing members have higher extinction rates than latecomers and that extinction probability declines with time (Saulsbury et al. Reference Saulsbury, Parins-Fukuchi, Wilson, Reitan and Liow2024). It is also consistent with findings that long-lived taxa of an assemblage tend to be more widespread, less specialized, and tend to have lower evolutionary rates than younger taxa (e.g., Saupe et al. Reference Saupe, Qiao, Hendricks, Portell, Hunter, Soberón and Lieberman2015; Raia et al. Reference Raia, Carotenuto, Mondanaro, Castiglione, Passaro, Saggese and Melchionna2016; Plotnick and Wagner Reference Plotnick and Wagner2018; Tamre and Parsons Reference Tamre and Parsons2024).
Notably, against expectation, extinction rates do not decline with scaled time in the PMAs (Fig. 5, Table 1). This suggests a more complex dynamic than assumed under NT and not a simple sorting toward more widespread and generalist genera. This is also evidenced by the increase in β-diversity through scaled PMA time, which suggests the opposite, namely, that geographically restricted genera and specialists accumulate in the assemblages with time. Such a tendency is otherwise not directly apparent in our data, but the fact that the change in β-diversity is positively associated with changes in γ, with the origination proportions, and with occupancy (Table 3) suggests that it reflects diversification and specialization intensity.
Such a tendency toward specialization, in turn, has been suggested as a clade-inherent feature (not to be confused with the assemblage-inherent features discussed earlier) under the model of “weak directionality” (Raia and Fortelius Reference Raia and Fortelius2013; Raia et al. Reference Raia, Carotenuto, Mondanaro, Castiglione, Passaro, Saggese and Melchionna2016). Weak directionality is the tendency of lineages to produce more and more specialized types as time progresses (Raia et al. Reference Raia, Carotenuto, Mondanaro, Castiglione, Passaro, Saggese and Melchionna2016: p. 1). It is unlikely that such a lineage-inherent trend will prevail in a randomly assembled set of taxa. But when lineages of an assemblage share a common history, with an aligned succession of appearance, aging, and replacement, lineage-inherent trends become linked with assemblage trends. This is exactly the difference seen between the randomly selected intervals and the PMAs. PMAs are collections of genera that appear in a particular order, where the chances of genus originations are higher at the beginning of a PMA than at the end. Hence, originations are concentrated at the beginning of the PMA. This large group of genera that emerges at the beginning, theoretically (and empirically; see Supplementary Fig. 4) also has the longest life expectancy when the end of the assemblage is set as the absolute limit.
Assembly of Global-Scale Eco-genealogical Units
Ecological communities and community groups (sensu Boucot Reference Boucot1983) must always be understood as a succession of assemblages. The concept of ecological succession has a long history and continues to be a central concept in ecological research (e.g., Clements Reference Clements1916; Prach and Walker Reference Prach and Walker2011; Chang and Turner Reference Chang and Turner2019). From ecology, it is well known that arrival histories can have various effects on the later succession of a community via facilitative or inhibitory niche modification or other processes that are collectively known as priority effects (e.g., Connell and Slatyer Reference Connell and Slatyer1977; Almany Reference Almany2004; Fukami Reference Fukami2015; Debray et al. Reference Debray, Herbert, Jaffe, Crits-Christoph, Power and Koskella2022).
Directional change of community assembly at geological timescales or the assembly of the higher levels of the eco-genealogical Bretskyan hierarchy is less well studied, but it can be assumed that it is a ubiquitous feature from smallest to largest spatiotemporal scales. For instance, Walker and Alberstadt (Reference Walker and Alberstadt1975) described directional change in fossil assemblages at small local scale, such as in fossil reefs, in short regional stratigraphic time intervals, such as in storm-dominated sediments, and in large regional spatiotemporal mosaics of sets of strata. Allmon and Martin (Reference Allmon and Martin2014) describe directional change at the Phanerozoic scale in eco-evolutionary successions within the marine trophic chain. Antell and Saupe (Reference Antell and Saupe2021) review various Phanerozoic-scale bottom-up pathways of eco-evolutionary change within marine ecosystems. And Jablonski (Reference Jablonski2007: pp. 96–97) provides additional examples of long-term directional change occurring at all scales in the fossil record.
Thus, it is important to recognize that the assembly of Bretskyan units is not a purely random process, but it is conditioned on its prior history. It seems, that for most of the times during the Phanerozoic and at the global scale, conditioned assembly was the dominant scenario. Except for times of mass extinctions. During times of mass extinctions, contingency effects were likely to soar (see Gould [Reference Gould1989: p. 306] for the original argument; see also Jablonski [Reference Jablonski2005]; see Monarrez et al. [Reference Monarrez, Heim and Payne2023] for evidence of increased variability of extinction selectivity during mass extinctions). Subsequently, effects of incumbency (for the role of incumbency effects in macroevolution, see, e.g., Rosenzweig and McCord [Reference Rosenzweig and McCord1991]) and evolving conditional dependencies played an important role. The large-scale transitions between first-level PMAs in this sense can be seen as representing the boundary intervals of a long periods of conditional assembly of global-scale eco-genealogical units.
The increasing median ages during these long intervals could be the result of the interplay of two factors: The first is a kind of priority effect, because the assembly begins preferably with generalist taxa (especially after mass extinctions), which form the basis for the successive addition of novel taxa, and which have lower probabilities to go extinct. Assuming the model of weak directionality, the novel taxa become increasingly specialized as time progresses. This second factor leads to decreasing chances of new taxa reaching occupancy levels that would allow them to establish themselves (decline in origination rates) and increasing chances that the cohort of latecomer (i.e., young) taxa go extinct (increase in extinction rates). This model is in accordance with global-level empirical data presented in Boyajian (Reference Boyajian1992), where the median age of the cohort of extinct families of a stage is consistently lower than the total, except during some mass extinctions.
Notably, an increasing median genus age is not a feature of the lower PMA levels. There, age trajectories vary from PMA to PMA without a significant general trend. Consequently, in the second and third PMA levels, the median age of the extinction cohort can be higher than or equal to the median age of the total assemblage at any time during the PMA duration. In contrast to the first-level PMAs, the lower-level PMAs contain only a fraction of the total of the genus occurrences of a given stage. Moreover, the third-level PMAs are not global in scale. Consequently, the probability that genera that are part of these low-level PMAs occur elsewhere and independently before or after the actual PMA duration is much higher than in the first-level PMA. The lack of median age trends in these low-level PMAs, therefore, can be interpreted as reflecting the temporal peculiarities and complexities of regional faunal shifts and migrations during the PMA durations. Additionally, one major difference that sets the lower-level PMAs apart from the first-level PMAs and the random intervals is that in the former only a minor number of PMAs begin at or include mass extinction intervals with their effects of a radical setback of genus ages (see Fig. 2).
This model offers a hierarchical, multi-causal view of a conditional assembly of eco-genealogical units. It explains that, at a global level, a consistent pattern of increasing median age emerges, resulting from the interplay between short periods of mass extinction and long periods of conditional assembly.
Eco-evolutionary Entrenchment and Memory
Thus, the assemblage history can be considered as important a factor in the composition of eco-genealogical units as the intrinsic properties of the taxa involved. At global scale, the assemblage can be seen as a kind of storage process, where species and higher taxa function as entities in which a certain moment or snapshot in time is inscribed and preserved via genes, traits, behavior, developmental pathways, ecological, and paleogeographic relations of the taxa. The clades preserve certain conditions that prevail during their subsequent persistence. In this sense, they represent a memory of the eco-genealogical units (memory here seen as the capacity of past states to influence the presence; see “Introduction”).
The increasing median genus age within first-level PMAs hence suggests that with time progressing, the margin between memory length (in the sense of Ogle et al. [Reference Ogle, Barber, Barron-Gafford, Bentley, Young, Huxman, Loik and Tissue2015]) and the present increases. The decrease in origination rates also suggests that this increase is linked with decreasing opportunities for new species and higher taxa to evolve. And finally, the positive association of β-diversity with age (Table 1) suggests that as time progresses and genera become older, the spatial differences between assemblages increase. In other words, the assemblages become increasingly eco-evolutionary entrenched and constrained when time progresses and memory increases.
The process can be compared with a community succession in a spatial mosaic. Hendry and McGlade (Reference Hendry and McGlade1995) demonstrated with cellular automaton models that memory length (here measured as the longevity of the individuals of a taxon) affects the structure of the spatial mosaic of an ecological community; in a community with long memory, the mosaic is highly aggregated and structured, it has random distributions when memory is absent (see also Peterson Reference Peterson2002).
This, in turn, is a process that can be compared with the hypothesis that during the evolution of clades, possible developmental pathways become increasingly nested and interdependent, resulting in an increase of developmental constraints (e.g., Galis et al. Reference Galis, Metz and Van Alphen2018) and in a “generative entrenchment” (Wimsatt Reference Wimsatt and Bechtel1986). Such effects of memory on temporal dynamics in networks also have been shown in random walk models where memory (defined here as the influence of previous steps in the random walk on the probability of subsequent step) increases the constraints on flow and structures the pathways taken (Rosvall et al. Reference Rosvall, Esquivel, Lancichinetti, West and Lambiotte2014).
The eco-evolutionary entrenchment in this scenario, and almost in a dialectic of remembering and forgetting, would only be interrupted when the eco-genealogical unit collapses and when memory is reset. This is what has been termed by van Valen (Reference Van Valen1984: p. 51) as a resetting of “the clock of community evolution.”
Here, such a scenario is apparent only at the global scale. In lower-level PMAs, the reset during PMA transitions (i.e., eco-genealogical unit transition) has less impact on median age trajectories, because here the genera more often also occur outside the spatiotemporal limits of the respective PMAs. It is therefore possible that within lower-level eco-genealogical units, similar pattern of conditional assembly and entrenchment are apparent when analyzed at the species level. The existence of a shared history of taxa in long-lived eco-evolutionary units has long been recognized and is probably best expressed in the concept of “coordinated stasis” within EEUs (Brett et al. Reference Brett, Ivany and Schopf1996). A promising area of future research would be to search for lower-level eco-genealogical units, such as EEUs or EESUs, in which species’ spatial and temporal limits are constrained. Within these units, potential increasing median age trends could be investigated.
Conclusions
The chances that existing taxa are replaced by new taxa decreases over Phanerozoic time, and older taxa persist for longer and longer. This long-term trend has been interpreted as resulting from repeated culling and non-replacement of groups with high turnover or, in other words, from long-term selection for extinction resistance. It has therefore been assumed that the respective characteristics of the taxa are decisive for this trend.
However, taxa are part of evolving ecological communities, regional species pools, biomes, and more broadly, of eco-genealogical Bretskyan units of life. Structural aspects of change of these evolving units, such as the composition of the evolving eco-genealogical units, has received little attention as a possible explanation. Ecological communities and the higher-level eco-genealogical units of the geological past are accessible via fossil assemblages. Here we combined a multilayer network representation with multilevel clustering to reconstruct the complex organization of the Phanerozoic benthic marine biota. The resulting modular structures are treated as PMAs and interpreted as representing global marine eco-genealogical or Bretskyan units. For each PMA, we calculated how its median genus age, its richness, the evolutionary rates, and the spatial β-diversity changed over time. The resulting estimates support the assumption that the assemblages represent meaningful eco-genealogical units, as the temporal pattern differs from randomly chosen Phanerozoic time intervals. Remarkably, in the PMAs, a waxing and waning pattern of the genus richness is apparent, which can be expected from ecological communities with successions of emergence, growth, decline, and collapse.
The PMAs also differ from the randomly chosen intervals in their combination of trends of decreasing proportions of originations and increasing β-diversities. This contradicts the hypothesis of a simple Phanerozoic sorting toward extinction resistance and generalization. Contrary to this expectation, assemblages become more specialized, sympatric, and geographically restricted over time. This is, in turn, a pattern that would be predicted from theory and empirics to occur as a within-lineage evolutionary trend. The within-lineage trends are apparent in the PMAs because they are not random sets of taxa, but they are ecologically meaningful successions in which clade histories are aligned. On average, the assemblages of the highest-level PMAs start with a high proportion of young genera, with high origination rates and low spatial differentiation. As time progresses, the genera and their descendants age and fewer and fewer new genera originate, which also have fewer chances to reach old ages. This process can be described as a conditioned assembly, where genera that appear early in an assemblage have the longest life expectancy. Hence, it is not necessarily the particular property of the genus that determines its later fate, but the timing of its appearance. At the highest hierarchical level—the global scale—the PMAs can therefore be seen as eco-genealogical units that were built during long periods of conditional assembly. They were bounded by short transition intervals, during which the assembly process was reset and median genus ages decreased dramatically. Declining origination proportions and the increasing β-diversities indicate that the assembly was probably accompanied by an eco-evolutionary entrenchment and by a weak within-lineage directionality toward specialization and geographic restriction.
Acknowledgments
This work is a result of ongoing discussions within the team of the Project “Comparing Evolutionary Processes in Nature and Society” funded by the Kone Foundation (Finland). We are grateful to L. Lahti (Turku, Finland), S. Nurmio (Helsinki, Finland), M. Peltoniemi (Tampere, Finland), A. Raulo (Oxford, UK), and I. Žliobaitė (Helsinki, Finland) for inspiring discussions on change in systems. The quality of the paper profited much from the constructive and supportive reviews of earlier versions by A. Spiridonov (Vilnius, Lithuania) and one unnamed reviewer. A.R. is supported by the Johannes Kepler University Linz, Linz Institute of Transformative Change (LIFT_C).
Data Availability Statement
Data available from the Zenodo Digital Repository: https://doi.org/10.5281/zenodo.19564332.
Competing Interests
The authors declare none.
Open access
