Research Article

A high-resolution summary of Cambrian to Early Triassic marine invertebrate biodiversity

See allHide authors and affiliations

Science  17 Jan 2020:
Vol. 367, Issue 6475, pp. 272-277
DOI: 10.1126/science.aax4953

A finer record of biodiversity

We have pressing, human-generated reasons to explore the influence of environmental change on biodiversity. Looking into the past can not only inform our understanding of this relationship but also help us to understand current change. Paleontological records depend on fossil availability and predictive modeling, however, and thus tend to give us a picture with large temporal jumps, millions of years wide. Such a scale makes it difficult to truly understand the action of environmental forces on ecological processes. Enabled by a supercomputer, Fan et al. used machine learning to analyze a large marine Paleozoic dataset, creating a record with time intervals of only ∼26,000 years (see the Perspective by Wagner). This fine-scale resolution revealed new events and important details of previously described patterns.

Science, this issue p. 272; see also p. 249


One great challenge in understanding the history of life is resolving the influence of environmental change on biodiversity. Simulated annealing and genetic algorithms were used to synthesize data from 11,000 marine fossil species, collected from more than 3000 stratigraphic sections, to generate a new Cambrian to Triassic biodiversity curve with an imputed temporal resolution of 26 ± 14.9 thousand years. This increased resolution clarifies the timing of known diversification and extinction events. Comparative analysis suggests that partial pressure of carbon dioxide (Pco2) is the only environmental factor that seems to display a secular pattern similar to that of biodiversity, but this similarity was not confirmed when autocorrelation within that time series was analyzed by detrending. These results demonstrate that fossil data can provide the temporal and taxonomic resolutions necessary to test (paleo)biological hypotheses at a level of detail approaching those of long-term ecological analyses.

Understanding patterns of global diversity can reveal the history of the biosphere and relations between environmental changes and diversity fluctuations, and can provide insights into how the fossil record might inform current biodiversity concerns. Early global-scale quantitative analysis identified what have come to be known as the “big five” mass extinctions. However, such efforts depend on the quality and temporal resolution of paleontological data, which have improved substantially since the 1990s, most recently through the intensive data compilation of the Paleobiology Database (16). Analyses of those data have increased our understanding of paleobiodiversity (4, 710).

Previous deep-time paleobiodiversity reconstructions (1, 11) were limited by coarse age determinations of taxon occurrences. The relatively long and uneven duration of age bins (stage or series level) used in these studies imposed complexly structured limits on resolving power across different intervals. Resolutions were generally no better than 8 to 11 million years (Myr) with standard deviations of 2.4 to 3.2 Myr, although some trials have been made to achieve better resolution for the early Paleozoic (6). Taxon age assignments were subject to error, not equally applicable to all clades, and quickly became outdated by new correlations or updated age estimates. Previous analyses have also been performed at taxonomically broad and phylogenetically suspect family or genus levels. Such resolutions are often too crude and imprecise to assess diversification rates or patterns associated with various global events (gradual, stepwise, or abrupt) and may mask multiple events as well as finer-scale fluctuations (7, 12, 13).

Here, we used a new parallel computing implementation of the constrained optimization method (CONOP.SAGA) run on the Tianhe II supercomputer. This approach uses inferred stratigraphic correlations to construct composite biodiversity curves for Cambrian to Triassic marine invertebrate genera and species (Fig. 1) and has demonstrated the capacity to establish finely resolved, traceable time zones over wide geographic areas (14).

Fig. 1 General trajectories of Paleozoic genus and species diversity and species diversity for 10 major fossil groups.

(A) Genus and species diversity. (B) Species diversity. The light and dark green shading in (A) represent 1σ and 2σ standard deviations, respectively, which are based on 500 bootstrap runs. 2σ approximately equals the 95% confidence interval. Gray bars on either side show the buffer zones with edge effects. 1, GOBE; 2, end-Ordovican mass extinction; 3, early Silurian radiation; 4, Middle to Late Devonian diversity decline; 5, late Carboniferous–early Permian biodiversification event; 6, end-Permian mass extinction.

Data and methods

Data compilation and standardization were conducted through the Geobiodiversity Database (15). This database is particularly suitable for biodiversity studies because, unlike the Paleobiology Database (Fig. 2), it is based on section data and provides quality control at a bed-by-bed level using an online, interactive system for recording expert taxonomic opinions (15). Taxonomic and age assignments used in this investigation were vetted by a team of 11 paleontologists, who checked and updated each taxonomic record. We also cross-checked these species names for synonyms.

Fig. 2 Comparisons between previous diversity curves and the present study for the Paleozoic.

(A) Diversity curves of previous studies. (B) ~10-Myr-binned species and genus diversity curves based on the present dataset. The high-resolution unbinned species diversity in Fig. 1 has been rescaled to the 26 uneven time bins adopted by Alroy et al. (25). The two y-axes on the right only refer to the red curves.

Because the Geobiodiversity Database records local taxon occurrences and their positions in stratigraphic sections, we were able to construct a composite sequence of assemblages and calibrate this sequence to a current estimate of the geological time scale using the best available chronostratigraphic data (16). Our study focused on marine invertebrates and used data from 3766 published stratigraphic sections, including 266,110 local records of the stratigraphic ranges of 45,318 taxonomic units, covering all Chinese Cambrian to Lower Triassic tectonic blocks (fig. S1). Although our data were largely derived from Chinese sections, the tectonic blocks on which they reside were situated in paleolatitudes stretching from southern Gondwanan to northern Boreal realms (17). Accordingly, these data reflect global biodiversity patterns (figs. S2 and S3).

Our initial analyses revealed that the rarity of Silurian–Devonian data in China (due to worldwide regression) hampered regional and global correlations. Consequently, we added a small amount of European Silurian–Devonian data to improve the correlations in this interval. These additional data did not alter the generality of our results because different Chinese tectonic blocks were located in different regions during the Paleozoic, with some residing close to Europe (fig. S3). Our study interval terminated at the late Middle Triassic marine regression.

Taxonomic names in open nomenclature, questionable taxa, and taxa unidentifiable to the species level were not included. Species recorded from only one locality were also removed to avoid the “monograph effect” (18). The resulting final dataset contained 116,060 local records of total stratigraphic ranges of 11,268 species from 3112 published stratigraphic sections.

To avoid the need to use coarse time bins, we used constrained optimization (CONOP) (19) stratigraphic correlation to reconstruct the Paleozoic biodiversity history of marine invertebrates. The CONOP correlation method, which applies a simulated annealing algorithm to infer a globally optimized sequence of stratigraphic datums, has been used previously for local high-resolution biochronostratigraphic studies (14, 20, 21). However, the original CONOP algorithm (22, 23) did not support parallel or high-performance computing and it would have required dozens of years to calculate one CONOP composite for this dataset. To overcome this “big data” problem, we modified the original CONOP algorithms to parallelize the sequencing problem. We also designed a special hybrid strategy of simulated annealing and genetic algorithm for the parallel computing application, CONOP.SAGA (16).

CONOP.SAGA iteratively compares species ranges from many local range charts to assemble the global first and last occurrence datums into a single, global, best-fit sequence, thereby reducing the effect of local-section incompleteness substantially (22). As is often the case in complex optimal modeling problems, CONOP.SAGA does not yield a single solution, but rather a set of equally optimal solutions. To obtain robust results, we created a medium-sized dataset and made >10 complete CONOP.SAGA calculations. If the results were similar for each run, then this suggested that the true global optimal solution had been found. Basic settings of the CONOP.SAGA parameters for the full dataset were estimated from these results. We then uploaded the full dataset and ran it while varying these parameters. Each calculation took 20 to 60 hours to complete and produced one composite sequence. Structured increase in these parameters allowed us to determine the effect of their variation on the diversity curve. Despite >1.2 billion iterations per run, all runs produced consistent results (>98.6%) with nearly indistinguishable variations. We then used the best (i.e., lowest misfit) (19) of the final three complete-dataset solutions to construct the final diversity curve estimate (Fig. 1).

Imprecision in estimating age and duration of biozones or other time intervals was prevented because local, isotopically dated species’ first and last appearance levels in the composite sequence were calibrated directly by regression against the high-precision geochronologic dates used to calibrate the Cambrian to Early Triassic geologic time scale (24, 25) with estimated ages for the bases of major international standard biozones (fig. S5 and table S2) (24). This calibrated composite is composed of 11,326 discrete temporal levels for the interval from 538.85 to 244.41 million years ago (Ma), yielding a mean temporal resolution of ~26.0 ± 14.9 thousand years (kyr).

Both species- and genus-level diversity summaries were assembled using an unbinned method (16, 20) to avoid issues associated with different counting strategies (26). A simple diversity curve was created by counting the number of species and genera occurring at each temporal level (Fig. 1A). Fossil data are inevitably biased by incomplete preservation and sampling across time, geographic range, and environmental setting. We used a bootstrap technique to estimate the range of variation of the species-diversity curve (Fig. 1A) (20). The diversity at each point was also standardized by average diversity and the number of sections to reduce effects of sampling effort (figs. S6 and S7) (16).


Our results (Fig. 1A and fig. S6) revealed a sharp increase in diversity associated with the Cambrian explosion, a pause through the late Cambrian Steptoean positive carbon isotope excursion event (27), followed by a nearly threefold increase in species diversity during the Early Ordovician. Species/genus ratio data suggest that the Great Ordovician Biodiversification Event (GOBE) reflects species-level diversification, which resulted in a facultative expansion of marine ecosystems. By contrast, the Cambrian radiation was associated mainly with increases at the genus level or among higher taxa (Fig. 1A). The GOBE (6, 28) is evident as a sustained 29.72-Myr-long diversity increase (from 497.05 to 467.33 Ma) until the Middle Ordovician. By contrast, the Alroy et al. diversity curve [(3); Fig. 1] showed a steady radiation from the early Cambrian until the Emsian (Early Devonian) with a minor drop across the Ordovician–Silurian transition (Fig. 2A).

The end-Ordovician mass extinction is evident as a rapid diversity decline from the late Katian to the Hirnantian; two phases corresponding to the development and disappearance of the Hirnantia fauna were previously identified (18, 29), but that distinction is less evident in our Chinese species-level summary (Fig. 1A and fig. S8). This interval may represent either a single event (30) or a community-level turnover, presumably in response to Hirnantian glaciation. An immediate recovery and radiation occurred near the Ordovician–Silurian boundary (444.66 Ma) and persisted until the early Silurian (Telychian: 437.08 Ma). Previous studies have not recognized this earliest Silurian radiation (4).

The Late Devonian Frasnian–Famennian extinction was broadly evident in the Alroy et al. curve (4) with a rebound in the Famennian. By contrast, our higher-resolution summary shows the Frasnian–Famennian event to be embedded within a protracted diversity decline that began in the Eifelian (392.72 Ma) and continued to ~368.24 Ma with no subsequent major recovery. No discrete Frasnian–Famennian event (31, 32) was evident in our results (Fig. 1A and fig. S6).

The Late Paleozoic ice age was accompanied by a previously unrecognized great biodiversification event (the Carboniferous–Permian Biodiversification Event) comparable to the GOBE on the basis of both the increase of species richness and species/genus ratio (Fig. 1A and fig. S9). The Alroy et al. result has this biodiversification event extending from the early Cisuralian into the middle Guadalupian (Figs. 1A and 2A and fig. S8) (4).

We found little evidence for an end-Guadalupian mass extinction, confirming previous qualitative results (33). The temporal resolution of earlier global compilations was too coarse to reveal the abrupt nature of the end-Permian mass extinction, but this is resolved clearly in our results. Both species and genus diversities declined rapidly from 252.73 to 251.95 Ma, followed by a sudden drop over ~63 kyr (Fig. 1A), which is consistent with the duration for this event calibrated with high-precision geochronology based on the Meishan section (table S1 and fig. S6) (34). This extinction was followed immediately by a minor Late Early Triassic radiation at 249.58 Ma, but marine diversity remained low through the Early and Early Middle Triassic (fig. S6).

To compare these results with diversity curves from the Sepkoski and Paleobiology databases, we generated generic and species-diversity patterns using ~10-Myr time bins with durations similar to those used by Alroy et al. (4). This coarse graining of the data yielded generally comparable patterns for the GOBE, the Early Silurian Radiation, and the Carboniferous–Permian Biodiversification Event, along with three major diversity crises (Late Ordovician, Middle to Late Devonian, and end-Permian) (Fig. 2). More importantly, this comparison shows the dependence of paleobiodiversity estimation on temporal resolution. If our high-resolution diversity pattern (Fig. 1A and fig. S6) is compared with the patterns of the same data using the ~10 Myr-bins, binning spreads out the end-Ordovician mass extinction and the Early Triassic radiation largely disappears (Fig. 2). These differences reflect the effects of low temporal resolution and uneven time-bin durations, varying from 4.97 to 32 Myr (fig. S11).

We used different time-bin durations (from 0.1 to 10 Myr) to calculate a series of diversity curves (Fig. 3) in addition to parsing our data into the unequal time bins used by Alroy et al. (4). The diversity pattern produced using 10-Myr bins is very close to these (4) but the influence of coarse and uneven sampling remained evident. Bins of up to 0.5-Myr duration retained a largely distortion-free representation of fine-scale structure in the present paleobiodiversity analysis (~300 Myr; Fig. 3). Previous studies suggest that resolutions as high as 0.1 to 0.5 Myr are needed to test causal-association hypotheses because this is the estimated duration of many physical perturbations to environmental states (e.g., large igneous province eruptions, bolide impacts, ocean anoxia events).

Fig. 3 Composite diversity patterns produced in different even time bins (0.1 to 10 Myr) to show the biases and effects of different temporal resolutions on paleobiodiversity estimation.

Dashed blue line represents the result fitting in the unequal time bins adopted by Alroy et al. (4).

In terms of Sepkoski’s three marine evolutionary faunas (35) (fig. S8) and different fossil groups (Fig. 1B and fig. S10), the Cambrian fauna represent the bulk of Cambrian biodiversity and remained substantial during the Ordovician, although they were later overwhelmed by diversification of the Paleozoic fauna. Trilobites, graptolites, and hyolithids lost their dominance after the Silurian.

The Paleozoic evolutionary fauna, consisting largely of brachiopods and other epifaunal filter feeders, dominated marine habitats from the Ordovician through the Devonian diversity decline. Conodont species diversity reveals two distinct radiations, one in the Early and Middle Ordovician and the other from the Middle Devonian to the late Carboniferous. Conodonts declined gradually during the early Carboniferous and experienced considerable losses during the Late Paleozoic Ice Age. Brachiopods experienced a rapid latest Ordovician to early Silurian diversification. This radiation was interrupted by a minor fluctuation in response to waxing and waning of the Hirnantian unipolar ice sheet. Brachiopod diversity declined during the Middle and Late Devonian before a steady increase began from the late early Carboniferous. The Lopingian had the most diverse brachiopod faunas of the Paleozoic, but this clade experienced catastrophic losses during the end-Permian mass extinction (36). Fusulinids exhibit a substantial diversification from the mid-Carboniferous to early Cisuralian and constitute an important component of the Carboniferous–Permian Biodiversification Event along with brachiopods and crustaceans. Fusulinid diversity declined slightly until the Lopingian, with the clade disappearing during the end-Permian mass extinction (Fig. 1B and fig. S10).

Mechanisms of diversity change

Many factors have been invoked as potential driving factors of biodiversity changes, including changes in paleoclimate, sea level, nutrient flux, ocean-atmospheric circulation, total habitable area, and intercontinental connectivity (8, 10, 3739). There are two challenges to establishing the relationships between environmental drivers and diversity: the need for high-resolution proxy data for appropriate environmental indicators (Fig. 4) and the confounding problems of potential cross-correlation of autocorrelated time series. Data of appropriate resolution are mostly lacking, particularly for changes in global and regional climate (38, 40). Currently available data suggest that the relationship between diversity and climate during the Paleozoic was complex. The GOBE and the Carboniferous–Permian Biodiversification Event were associated with climatic cooling (41, 42), which stands in contrast to the modern association between warm climates and diversity (Fig. 4). The early Silurian radiation was associated with overall warming until the middle Llandovery (43). The late Asselian (294.19 Ma) diversity peak coincided with the acme of Late Paleozoic Ice Age (41, 44). Species diversity declined from the Sakmarian to middle Guadalupian as climate ameliorated (45). The end-Ordovician extinction was associated with a short glaciation (46, 47).

Fig. 4 Correlation between the species-diversity trajectory and the trends of multiple environmental proxies.

(A) 87Sr/86Sr ratio (60). (B) δ13C (53). (C) δ18O (62). (D) Continental fragmentation index (10). (E) Sedimentary material (63). (F) Estimated Paleozoic Pco2 (58, 59). (G) Rescaled species-diversity trajectory of the present study.

The Middle to Late Devonian diversity decline also had a complex temperature history. Paleotemperature increased from the Eifelian–Givetian transition to the Frasnian–Famennian boundary, then decreased after the Frasnian–Famennian boundary (48). On the basis of the latest δ18Oapatite data, the Lower and Upper Kellwasser events were each associated with distinct cooling events (49). The end-Permian mass extinction was followed immediately by a rapid warming of 8 to 10°C (50, 51), but low diversity in the Early Triassic coincided with a lethal “hothouse” (52).

The carbon isotope record reflects changes in diversity and abundance that affect the global carbon cycle (53). Therefore, large-scale biotic events are often associated with large carbon isotope excursions (54). The end-Ordovician mass extinction is associated with an ~4‰ positive excursion of δ13Ccarb (55), and the end-Permian mass extinction includes an ~5‰ negative excursion (34). Carbon isotope excursions have been reported from the Late Devonian Frasnian–Famennian event (56) but do not appear to coincide with distinct diversity changes (54).

Changes in the atmospheric concentrations of Pco2 have had a major impact on earth system dynamics, but there have been major discrepancies between previous reconstructions of secular trends (57). A recent long-term Pco2 reconstruction (58) and the stable carbon isotopic fractionation associated with photosynthesis (59) show a secular trend similar to that documented by our late Silurian to Early Triassic diversity pattern (Fig. 4, F and G). The increasing trend in Pco2 is associated positively with increasing long-term diversity before the early Visean (58). After the post-Visean appearance of a more modern fauna, diversity increased more rapidly than the increase of Pco2. The similarity of these secular trends (Fig. 5, A and B) might imply that both Pco2 and diversity were responding to a common set of driving factors. However, the results obtained from further correlation analysis (Fig. 5, C and D, and fig. S13) do not show any significant correlation or causal association between Pco2 and diversity. Similarly, other environmental proxies, (e.g., 87Sr/86Sr ratio, δ13C, δ18O, continental fragmentation index) do not exhibit compelling trends with diversity changes (fig. S12). However, such comparisons are hampered by the absence of long-term high-resolution environmental proxy data (Fig. 4 and fig. S12).

Fig. 5 Correlation analysis between species diversity and Pco2 [data from (58, 59)].

Results in (A) show strong secular trends of both diversity curves and Paleozoic Pco2. These secular trends generally decrease from 253 to 342 Ma and then increase from 343 to 419 Ma. The Pearson correlation analysis without considering the effects of the trends gives the results, coefficient of determination R2 = 0.70, Student’s t-value = 13.2, and n = 77 for the curve from 419 to 343 Ma and R2 = 0.52, t-value = 9.72, and n = 90 for the other range between 342 and 253 Ma (B). To test for the effects of trends or autocorrelations in the time-series data, further correlation analysis was applied to the changes in species diversity and the changes in Pco2 according to Yule’s interpretation. The changes of values are calculated by first difference (Diff) method (C and D) in which each point is subtracted from the point that came before it. The results indicate that there is no significant correlation between the two curves of changes.

Evidence for a decline in continental block fragmentation from the Cambrian to the Triassic (10) is negatively related to the general increase in diversity (Fig. 4E). Changes in the 87Sr/86Sr ratio are related to continental breakup, oceanic midridge spreading, and more intense continental weathering. Comparison between the 87Sr/86Sr ratio profile and the diversity pattern shows that each of the Paleozoic diversity crises (end-Ordovician, Middle to Late Devonian, and end-Permian) coincides with transitions from low 87Sr/86Sr ratios to increasing values, suggesting a reduction in seafloor spreading, an increase of the continental weathering, or both. However, the GOBE is associated with a steady decline in the 87Sr/86Sr ratio, whereas during the Carboniferous–Permian Biodiversification Event, the 87Sr/86Sr ratios were high (Fig. 4A), indicating increased continental weathering (44, 60).


The combination of our new Chinese data compilation, a new parallel computing implementation of CONOP.SAGA stratigraphic correlation algorithms, and the parallel processing power of the Tianhe II supercomputer have allowed the construction of a high-resolution composite species-diversity history with an average resolving power of 26.0 ± 14.9 kyr. Our results indicate that the coarse and uneven temporal resolutions used by previous summaries artificially influenced paleobiodiversity estimations. This analysis confirms the existence of end-Ordovician and end-Permian mass extinctions, a long-term Middle to Late Devonian diversity decline, and a markedly subdued Frasnian–Famennian event. Three biodiversification events are also evident in our results: (i) from late Cambrian to Middle Ordovician, (ii) in early Silurian, and (iii) from late Carboniferous to Cisuralian. The proposed mid-Carboniferous and end-Guadalupian “mass extinctions” are only evident as minor species-diversity fluctuations. A recent long-term Pco2 reconstruction (59) shows similar secular trends with the diversity data between the Silurian and Early Triassic. However, correlation between these two types of time-series curves shows strong autocorrelation and definite cause-and-effect links between specific environmental factors and the diversity changes require more study because of the lack of long-term, high-resolution, environmental proxy data and/or discrepancies between the dates assigned to our biodiversity curve and the environmental proxy curves.


Materials and Methods

Figs. S1 to S13

Tables S1 to S3

References (6476)

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: We thank G. Hunt, M. Foote, and an anonymous reviewer for their helpful reviews, which greatly improved the manuscript. We also thank Z. Shen for improving Fig. 5. Funding: This research was supported by the Strategic Priority Research Programs of the Chinese Academy of Sciences (grant nos. XDB26000000, XDB18000000, and XDA19050100), the National Natural Science Foundation of China (grant nos. 41725007 and 41420104003), and the National Key R&D Program of China (grant no. SQ2018YFE020456). The data from Tibet are supported by the Second Tibetan Plateau Scientific Expedition and Research (grant no. 2019QZKK0706). This paper is a contribution to the IUGS “Deep-Time Digital Earth” Big Science Program. Author contributions: J.F. and S.S. designed the research. S.S. and J.F. wrote the main text and supplementary materials with inputs from D.H.E., N.M., P.M.S., and Q.C. and data contributions from all authors. J.F., X.H., J.Y., and P.M.S. designed the CONOP.SAGA program for parallel computing. J.F. and S.S. performed most of the analyses. Competing interests: The authors declare no competing interests. Data and materials availability: All data and code are available from the OneStratigraphy ( and Dryad (61) repositories.

Stay Connected to Science

Navigate This Article