Earth’s interior dynamics drive marine fossil diversity cycles of tens of millions of years

We analyze short-term cyclicities in generic biodiversity from the Paleobiology Database (PBDB) and in macrostratigraphy from the Macrostrat database (Materials and Methods). We analyze various PBDB biodiversity datasets sampled at different timescales in order to evaluate the reproducibility of potential cyclicities. We also analyze another marine genus biodiversity data compilation, previously used to explore its potential cyclic makeup (12). We attempt the analyses first over the entire Phanerozoic Eon (past 542 Myr) and then over shorter intervals to check the continuity of the cyclic components. Separately, the Paleozoic and Ceno-Mesozoic intervals show distinct cyclicities (SI Appendix, Figs. S1–S10). In this study, we focus on the latter interval to assess potential correlations with a well-constrained plate tectonic model spanning the last 250 Myr of Earth history (33).

Spectral analysis of various PBDB biodiversity datasets of the Ceno-Mesozoic interval shows two common cycle bands of ∼28 to 29 and ∼36 Myr (Fig. 1, see also SI Appendix, Figs. S1–S10). Spectral analysis of the macrostratigraphic data of total marine sediment areas in North America (Materials and Methods) reveals two cycle bands of ∼29 and ∼36 to 37 Myr (Fig. 1). Finally, spectral analysis of the macrostratigraphic data of carbonate sediment area in North America (Materials and Methods) shows a stronger ∼36 Myr signal but a weaker 27.9 Myr signal (Fig. 3 and SI Appendix, Fig. S11). Results based on Rohde and Muller’s (12) biodiversity dataset point to the reproducibility of the 28 and 36 ±1 Myr cycle bands (SI Appendix, Figs. S12–S17). We have further tested the existence of the 28 and 36 ±1 Myr cycle bands using an exploratory search for periodicities, e.g., the phase distance correlation periodogram (PDCP, see Materials and Methods). In particular, PDCP of the Ceno-Mesozoic shows significantly stronger peaks of ~27 and 37 Myr (at lower iteration orders; see SI Appendix, Figs. S18). Additional statistical frequency decomposition and analysis, adaptive to geological variations, based on Gaussian processes and stochastic harmonic oscillator (SHO) kernels (Materials and Methods) further support the results obtained by the PDCP (SI Appendix, Fig. S19). We also find that the 28 and 36 ±1 Myr cycle bands are the main periodic components in the variations of the duration of the global geological stages (Fig. 1). This implies that the diversity of preserved fossils and the amount of sedimentary rocks covary with the duration of stratigraphic intervals (stages). Sedimentary rock quantities (areas), global occurrences, and stage-duration variations share similar long-term modulations at the 28 and 36 ±1 Myr cycle bands (Figs. 1 and 4 and SI Appendix, Figs. S20 and S21). This close cyclic relationship between them hints that they are linked via preservational/sampling processes. Additionally, the 62 Myr biodiversity and sediment rock cycle (12, 34) tracks the mean duration of geological periods (Fig. 1 and SI Appendix, Fig. S8). At the same time, we show that variations in the duration of geological stages are modulated by epochs, and epochs are modulated by periods (SI Appendix, Figs. S20 and S21). Epochs and Periods, and especially their boundaries, have been traditionally erected based on major paleobiological turnovers (SI Appendix, Fig. S21). The stages are time-rock units, chosen as the fundamental building blocks of chronostratigraphy (e.g., ref. 35). Divisions of the dynamic international stratigraphic scale are principally defined by the fossil content of the rocks and most primary markers for Stages are biostratigraphic events (first and last occurrences of fossil taxa). This implies that important paleoenvironmental changes preferentially coincide with stage boundaries: stages are, in other words, not arbitrary subdivisions of geologic time but instead coincide with inflections in the state of the Earth system that likely affect biodiversity and environmental turnover. The tectono-eustatic Wilson megacycle, which unambiguously mirrors major paleogeographic and paleoenvironmental changes, is detected in both epoch- and period-duration variations (SI Appendix, Figs. S20 and S21), defining roughly the Ceno-Mesozoic and Paleozoic intervals. We found cyclicities of ~28 and 36 ±1 Myr in Earth’s interior-related processes, as well as in various paleoenvironmental proxy datasets resolved independently of geological stage durations. The deep-sea stable oxygen isotopes (δ¹⁸O) show pronounced 36 Myr trends over the past 110 Ma (30) but do not correlate with stage-duration variations (SI Appendix, Fig. S22). The 36 Myr δ¹⁸O cycle has also been recorded over the past ~200 Ma along with a strong correlation with the sea-level record (36), see also SI Appendix, Fig. S14). This weakens the sampling-bias hypothesis and further supports the common-cause mechanisms. Finally, we test the effect of the sampling process on the record of cyclicities of tens of Myr and show the robustness of ~28 and ~36 ±1 Myr cyclicities (SI Appendix, Figs. S23 and S24).

The two close ~28 and 36 ±1 Myr geological cyclicities have long been considered within the same frequency band, the mean period of which evolves through time (14), although a recent study has argued for their individuality (38). Additionally, the 62 Myr cycle is another feature in biodiversity data (12).

We perform synthetic tests in order to check whether the two spectral peaks at ~28 and ~36 ±1 Myr represent distinct cycles and to test their eventual link to the 62 Myr cycle. We analyze the macrostratigraphic datasets because their spectra significantly highlight the two ~28 and ~36 ±1 Myr peaks, as well as the 62 Myr peak (Fig. 1). Results show that the ~28 and ~36 ±1 Myr cycles are recorded throughout the last 250 Myr (SI Appendix, Figs. S25 and S26), and their AM produces a long-period cyclicity that we cannot precisely determine because the analyzed 0 to 250 Ma time interval is too short for this purpose (Fig. 4). The existence of a long-term AM cycle resulting from the interference of ~28 and ~36 ±1 Myr cycles over the 250 Myr time interval reinforces the idea that the ~28 and ~36 ±1 Myr cycles are distinct (4). Additionally, we demonstrate that these two cycles are independent of the 62 Myr cycle by analyzing fossil data that contain strong ~28 and ~36 ±1 Myr signals, but devoid of the 62 Myr signal (4). All details on these analyses are provided in SI Appendix, section 6. In summary, the most continuous, common ~28 and 36 ±1 Myr cyclicities among biodiversity and macrostratigraphic datasets in the Ceno-Mesozoic interval are unlikely to solely be statistical artifacts of sampling bias but are instead most likely the result of common-cause mechanisms.

While the ~28 and 36 ±1 Myr cyclicities are common in fossil and sediment records, the 62 Myr cycle is seemingly not significant in sampled-in-bin generic diversity and in global occurrences (Fig. 1). It is likely that the studied 0 to 250 Ma time interval is not long enough to capture the longer cyclicity for its periodicity to become discernable since it would require at least five cycle repetitions. Therefore, we focus on the record of the ~28 and 36 ±1 Myr cycles and their potential causes.

The global subduction rate and variations in the areal extent of oceanic crustal production also show ~26 and ~35 Myr cyclicities, with a relative dominance of the ~35 Myr cycle (Fig. 2 and SI Appendix, Fig. S27). This ~35 Myr tectonic cycle is correlated to the longer-term global sea-level change (Fig. 3). However, the global subduction flux (plate volume) does not show significant cyclicities (SI Appendix, Fig. S28). This suggests that the evolution of the mean age and thickness of the subducting plate through time is unlikely to be periodic (SI Appendix, Fig. S29).

We additionally show a prominent ~36 Myr cyclicity in the global median trench migration dataset (Fig. 2 and SI Appendix, Figs. S30 and S31), which complements the ~35 Myr subduction and oceanic crust production cycle. These results have been interpreted as suggesting that oceanic plate subduction drives much of the large-scale global tectonic system on Earth and that the ~35 Myr cyclic changes in the global sea level are the result of both accretion at mid-ocean ridges (MORs) and slab subduction processes (38). Finally, we find a striking ~37 Myr cyclicity in the length of the MOR system (Fig. 2). Such cyclicity may reflect subduction-paced changes in the elongation of the MOR system, further augmenting the hypothesis of a subduction-paced global tectonic system. There is another 26 to 28 Myr spectral peak, with a lower amplitude in trench migration and ridge length data and a stronger amplitude in spreading rate data (Fig. 2 and SI Appendix, Fig. S27). The ~28 and 36 ±1 Myr tectonic cycle bands seem coupled and in tune with spreading and subduction processes.

The 36 ± 1 Myr signal is the common, dominant cycle in most of the tectonic datasets (Fig. 2) and in the sea level record as well (36); hence, we proceeded to its extraction for correlation among the studied geological proxy data (Fig. 3 and SI Appendix, Fig. S35). The prevalent 36 ± 1 Myr cyclicity found in biodiversity and in tectonic and eustatic data reinforces the idea that the large-scale tectonic regime-driven sea-level change is the major driving force of the longer-term evolution of the marine biosphere (5, 29). Most of the 36 ± 1 Myr biodiversity cycle minima correlate with their equivalent eustatic sea-level minima (Fig. 3 and SI Appendix, Figs. S14 and S35, further discussion below). Furthermore, the timing of the two major biotic extinctions, at the end of the Ordovician and at the Permo-Triassic transition, tracks the 36 ± 1 Myr cycle extremes in sea-level data, suggesting that lowered sea levels would have contributed, along with other Earth’s interior-related processes (e.g., volcanism), to lowered diversity of marine biota.

We should also note that these biotic cycles are not always periodic, which is a common feature in quasi-cyclic geological archives (36). For instance, the mean ~36.5 Myr periodicity over the 66 to 331 Ma interval (SI Appendix, Fig. S14E) has a longer ~39 Myr biodiversity oscillation within the Mesozoic-Paleozoic transition, while its equivalent eustatic cycle is shorter (SI Appendix, Fig. S14). Expansion of the ~36.5 Myr mean biodiversity cycle in this specific interval is due to the major biotic crisis at the Permian-Triassic transition, where the prolonged interval was due to a longer recovery period. This suggests that the five major biotic crises were paced by a complex of contributions from major glaciations, extensive volcanism from large igneous provinces, sea-level, and (at least one) bolide impact. For instance, the Cretaceous–Paleogene extinction event has two potential triggers from the Chicxulub impact (39) and the Deccan Traps (40). It has also been suggested that during some specific Earth’s climate states, such as the Late Ordovician cooling event, combined effects of paleogeographic configuration, cooling, and draining of epicontinental seaways might explain drops in biodiversity (41). These combined effects of environmental dynamics paced by climate and tectonics, together with ecological and biological evolutionary processes (31) may disrupt the cyclic nature of biodiversity.

The correlation of the 36 ± 1 Myr biodiversity cycle with sea level, subduction, oceanic crustal production, median trench migration, and the length of the MOR suggests a connection between tectono-eustatic processes and the quasi-cyclic nature of biodiversity. To further investigate these linkages, we perform cross-phase and coherence spectral correlations (Materials and Methods) at the common 36 ± 1 Myr cycle band for the multiple geological datasets (Fig. 3 and SI Appendix, Figs. S27 and S35). Sea level leads biodiversity by a few Myr (1.9 Myr mean phase) as far back as ~130 Ma. From ~130 Ma to the present, a significant phase lead in the order of 6 to 10 Myr (9.1 Myr mean phase) can be noticed. At around 120 to 130 Ma, a jump in sea level and biodiversity is seen, approaching the Wilson megacycle maximum expressed in the eustatic data (SI Appendix, Fig. S14). There is an apparent negative correlation between sea level and biodiversity at 120 to 130 Ma sea-level maximum (Fig. 3), which is most likely the result of the natural interference between the 28 ±1 and 36 ± 1 Myr cyclicities in the fossil and sediment data (Fig. 4). At the 36 Myr cycle band, there is still an effective positive correlation between sea level and biodiversity, but with a change in phase amount as mentioned above (SI Appendix, Fig. S35).

Subduction and crustal production are in-phase (Fig. 2 and SI Appendix, Fig. S27), reflecting that subduction processes have a direct effect on crustal production and destruction (33, 38). This is further supported by the in-phase relationship between subduction and median trench migration (Fig. 2 and SI Appendix, Fig. S30), which is a measure of interaction between subducting slabs and the upper mantle/transition zone. Biodiversity and crustal production show an antiphase relationship, with crustal production in the lag phase with respect to biodiversity. It has a similar response as the above-mentioned phase lag between biodiversity and sea level, i.e., there are two distinct intervals of phase lag. In the interval ~250 to 130 Ma, crustal production lags biodiversity by only 0.9 Myr, and then from ~130 Ma to the present, there is a significant lag of about 6 to 8 Myr (6.5 Myr mean phase). In summary, greater crustal production corresponds to lower biodiversity. The presence of a dominant ~36 Myr signal in carbonate macrostratigraphy to the detriment of the ~28 Myr signal (SI Appendix, Fig. S11) is also a feature in sea-level data (SI Appendix, Fig. S14), suggesting a potential link between carbonate and sea-level variations. Although the coherence between the two ~36 Myr carbonate and sea-level signals is moderate to weak (0.49, Fig. 3) possibly because of the unstable phase between the two proxies and/or the moderate quality of Meso-Cenozoic carbonate data (Materials and Methods), the shared ~36 Myr cyclicity provides a good argument for a potential link between carbonate sedimentation and sea-level processes. The Paleozoic carbonate macrostratigraphic data that are of high quality (Materials and Methods) show a 36 Myr cycle strongly correlated to sea-level data (SI Appendix, Fig. S36). This supports the hypothesis of sea-level driver of shelfal carbonate sedimentation. (Fig. 5).

In our results, sea level and oceanic crust production show an antiphase relationship, implying that the greater the rate of oceanic crustal production, the lower the sea level. For the same reason, sea level and the length of the MOR are also in antiphase. This antiphase relationship suggests that changes in ocean basin volume alone do not control long-term sea level. In addition, reconstructed basin volumes from the same tectonic model show an acyclic nature (Materials and Methods and SI Appendix, Fig. S32). This leads us to suggest that long-term sea level may be driven by a mechanism that was first hypothesized by Cloetingh and Haq (42), i.e., water entrainment at downwelling subduction zones and expulsion at the ridges associated with changes in seafloor production rates (see also refs. 43–45). Here, we provide support for this mechanism as implied by the 36 ± 1 Myr cyclicity, that is common to both sea level and subduction processes.

We modeled deep-water flux along subduction zones and deep-water reflux through the MOR system following Karlsen et al.’s (44) equations (Materials and Methods). Subduction should drive the 36 Myr deep-water cycle since the water percolating through subduction zones is much greater than water discharged via MORs (44, 46, 47). This model can explain the above antiphase correlation between sea level and oceanic crustal production data.

Time-series analysis of the modeled deep-water data shows a significant 35 Myr cyclicity in MOR deep-water reflux (Fig. 2). On the other hand, the subducting water flux in this model does not exhibit a significant ∼36 ± 1 Myr cycle (SI Appendix, Fig. S33), which could be because the model underpinning this calculation is built on the assumption that the volume of subducting water is dominated by mantle lithospheric water. Karlsen et al. (44) assumed that the older and thicker the down-going lithosphere is, the more deeply it cracks along the subduction hinge before entering the trench, allowing water to infiltrate and hydrate the mantle more pervasively, implemented as a simple linear relationship. Because changes in the mean thickness of subducting lithosphere over time are acyclic (SI Appendix, Fig. S29), our computed deep-water subduction flux is also acyclic since it is constructed on the basis of the subducting plate thickness (Fig. 2). Although relatively little is known about the hydration of oceanic mantle (see refs. 48–51), our analyses suggest that deep-water subduction may be less strongly dependent on mantle lithospheric water than assumed and instead is more likely to be driven by subducting crustal and sedimentary water. In that case, the rate of subducting plate area, which exhibits a pronounced 35 Myr cyclicity, would drive the long-term eustatic sea level change. During fast rates of subduction, deep water is excessively transported from the oceans into the mantle through subduction zones, while during that time, a deficit of deep water is returned to the oceans via the MORs due to a time lag (42). This scenario reconciles the potential link between sea-level lows and crustal production highs and longer MORs (Fig. 3 and SI Appendix, Fig. S35). Subduction-driven deep-water and sea-level data share a prominent, correlatable ∼36 Myr cyclicity (Fig. 3).

The cyclic nature of the biodiversity variations, in particular, those at tens of Myr time scales, has been debated since the 1980s (1, 8–12, 15). The 28 and 36 ±1 Myr cycles are recorded with almost equal amplitudes in biodiversity and macrostratigraphic data (Fig. 1). The 28 ± 1 Myr cycle is characterized by a weak spectral peak (low amplitudes) in tectonic and sea-level data, while the 36 ± 1 Myr cycle is prominent in these data (Fig. 2). We therefore focus on the correlation of the 36 ± 1 Myr cycle among the studied geological archives looking for a potential causal mechanism (Results).

Our finding of the 36 ± 1 Myr common, correlatable cyclicity in subduction and sea level datasets and the anti-phase relationship between sea level and crustal production suggest together that the cyclic deep-water exchanges paced by subduction processes govern cyclicities in long-term sea level (see Section above). At this same cycle band, biodiversity data correlate, in turn, to the tectonic and sea-level data, suggesting that tectono-eustatsy is the main driver of the 36 ± 1 Myr biodiversity cycle via changes to shallow-marine environments.

The 36 ± 1 Myr tectono-eustatic cyclicity is modulated by the supercontinental Wilson megacycles (36) (SI Appendix, Fig. S14). The raw diversity data faithfully capture the Paleozoic Wilson megacycle, while the Ceno-Mesozoic interval shows a strong, linear trend with a sharp rise toward the present (SI Appendix, Fig. S14). The record of the 36 ± 1 Myr equivalent biodiversity cycle superimposed on the Wilson megacycles further suggests that the evolution of life and geodynamic processes are coupled at these two timescales, in tune with deep-water recycling (44) within the lithospheric upper mantle (48).

The potential link between Earth’s interior dynamics and biodiversity has long been postulated at the Wilson megacycle timescale (5, 52–55), and the evidence supporting the tectono-eustatic driver of biodiversity at the Wilson cycle scale continues to grow (56, 57). The aggregation of continental blocks reduces provinciality (by reducing niche areas) and hence biodiversity, while their desegregation increases provinciality and biodiversity (52). At tens of Myr timescales, where there is less evidence for continental fragmentation (33), there is support for the idea that continental flooding is the major driving force through augmentation of ecological niches over epeiric shelves and interior seaways and their expansion/contraction led in turn to changes in biodiversity (26, 29). Models of recent and deep-time spatial biotic diversification successfully reproduce biodiversity patterns in tune with dispersal ability and geographic positions. For example, diversification models using deep-time reconstructions of habitat suitability have shown that marine shallow-water dynamics, including plate tectonics, influence reef and fish diversity (58). This hypothesis is reinforced by the correlation of periodicities between sea level and shallow-marine carbonate sedimentation, as this latter process is vulnerable to changes in bathymetry. Carbonate stratigraphic data (26) (Materials and Methods) contain a ~36 Myr signal (SI Appendix, Fig. S11), which correlates to the global sea level signal (Fig. 3 and SI Appendix, Fig. S36). The record of the ~36 Myr cycle in shelfal carbonates constitutes a salient argument for the connection between tectono-eustatic forcing and biodiversity via the precarious nature of shallow-marine environments that are extremely sensitive to depth variations. We conclude that greater magnitude sea-level changes induced by tectonically recycled deep water in the lithospheric upper mantle have the potential to drive the diversity of the life via cyclic continental flooding process throughout geologic time. Although we suggest that tectono-eustasy as the dominant driver of biodiversity patterns via expanding and contracting ecological niches on shelves and on epeiric seas, multiple processes can act and interact with different strengths through time and across large spatial extents, such as ecological and biological evolutionary forces (31). Speciation, extinction, dispersal, and adaptation all act and interact simultaneously within an environment (59, 60) to produce observed biodiversity patterns (61). Also, the dominant long-term forcing processes of biodiversity can be different over various time scales. While we argue for a dominant tectono-eustatic driver of biodiversity at the 36 Myr cycle band, other driving mechanisms may interfere at various cycle bands (see Section below).

Our analyses of geological cyclicities of tens of Myr reveal two potential inflection phases during the Ceno-Mesozoic, which may mirror potential turnover in biodiversity and the forcing processes. The first inflection phase in the early Cretaceous, at around 120 to 130 Ma, corresponds to a node (minimum) in a very long-period biodiversity cycle, originating from the interference (AM) of 28 and 36 ±1 Myr cycles (Fig. 4 and SI Appendix, Fig. S25). This inflection phase also matches the highest sea level, reflecting the expression of the Wilson tectono-eustatic megacycle maximum (36, 38, 43), and corresponds to the highest tectonic activity during the studied Ceno-Mesozoic times (33), matching the Cretaceous slab superflux spike (33 and SI Appendix, Fig. S28). The important phase lag between biodiversity and tectonics at 120 to 130 Ma may reflect the interference of other processes in driving biodiversity, i.e., tectono-eustasy is not the only dominant driver of biodiversity at this time interval. We suggest additional contributing factors such as the atmospheric carbon dioxide (CO₂) concentration. Models and proxies document higher levels of CO₂ during the Early Cretaceous (62, 63). The 120 to 130 Ma inflection phase separates two distinct time intervals of significant changes in the amplitudes of ~28 and 36 ±1 Myr biodiversity cycles (Fig. 1 and SI Appendix, Figs. S1 and S2). Amplitudes of these cycles are weaker before 120–130 Ma, and stronger after, suggesting a major biotic turnover likely paced by multiple driving forces including tectono-eustasy, volcanism, CO₂, and weathering (31, 41).

The second inflection phase in the Late Triassic, at around 210 Ma, shows an apparent switch from a dominant 62 Myr cycle in the Paleozoic to dominant 28 and 36 ±1 Myr cycles in the Ceno-Mesozoic. The 28, 36, and 62 Myr biodiversity cycles are recorded throughout the Phanerozoic, but the 62 Myr biodiversity cycle is seemingly more pronounced in the Paleozoic (Fig. 1 A, F, and G and SI Appendix, Figs. S1–S3). The 210 Ma inflection point also tracks the beginning of a strong linear trend toward the present in the biodiversity variation (SI Appendix, Figs. S1 and S2). This switch in the dominant cycle periodicities together with changes in biodiversity trends may be the result of changes in the very-long-term forcing processes.

Although a combination of multiple forcing on biodiversity variations is likely, we suggest that the dominant driver of the shorter 36 ± 1 Myr cycle is tectono-eustasy, while the longer 62 Myr biodiversity cycle may be the result of dominant atmospheric (CO₂) forcing, which was likely further amplified during the Paleozoic. Models show that long-term global atmospheric CO₂ was potentially higher during the Paleozoic than during the Meso-Cenozoic (63). Excess CO₂ would accentuate weathering processes, which were likely more intense during the Paleozoic, as attested by the prominent 62 Myr cyclicity in marine strontium isotope ratio ⁸⁷Sr/⁸⁶Sr (SI Appendix, Fig. S37). Significant changes in weathering processes during the Paleozoic can also be supported by the substantial drop of CO₂ at the Carboniferous-Permian, corresponding to the late Paleozoic Ice Age. This glacial event has been postulated to be controlled by the uplift of the Hercynian mountains in the tropics, as a potential trigger for intensified rainfall and erosional conditions, driving high rates of silicate weathering (64). Amplitudes of variation in subduction processes during the Paleozoic were likely more important than during the Ceno-Mesozoic, which is a controlling factor for CO₂ degassing increase (65).

The two potential inflection phases in the biotic turnover, as well as the driving mechanism of the 62 Myr biodiversity cycle, should be investigated in future studies by integrating more proxies and datasets. Specifically, the 62 Myr cycle may represent a major Earth system function since it tracks boundaries of the global geological periods (Fig. 1A and SI Appendix, Fig. S8). Finally, although the 28 and 36 Myr cycles versus the 62 Myr cycle are genetically in tune with the major geological subdivisions (stages, periods, etc.), their dominant underlying mechanisms are likely different possibly because of changes in very long-term Earth system states.

Marine and marginal marine biodiversity datasets are extracted from the PBDB using several approaches. Range-through diversity, which tabulates diversity based on the ages of first and last occurrences of each genus, was computed for 42,991 genera (42,952 of which are Phanerozoic) using 924,538 occurrences from 170,070 collections that are designated as marine and/or marginal marine in their environment of deposition. Occurrences from collections with no environmental designation in the PBDB were not used in this tabulation. Sampled-in-bin diversity was tabulated using the ages provided for each marine/marginal marine collection in the PBDB. In all cases, the taxonomy resolved by the PBDB was used (e.g., all synonyms were resolved prior to tabulation), and no attempt was made to exclude genera from marine/marginal marine collections on the basis of their classification, biology, and/or environmental preferences. Thus, terrestrial genera found in marine and marginal marine sediments are included here, although these likely compose a very small fraction of the total number of occurrences. All of the data in the PBDB are accessible via the API (66), updated continuously as new data are acquired (see https://paleobiodb.org/classic/displayDownloadGenerator). For convenience, and to allow precise reproduction of these results, concise raw data tables providing collection and taxon numbers, as well as PBDB ages used in this analysis, are provided in Datasets S1 and S2.

The macrostratigraphic datasets for North American marine sediments were extracted using the Macrostrat API (67) (see macrostrat.org/api). A total of 18,207 marine sedimentary rock units from 940 columns in North America were used to tabulate the marine sediment coverage area. Of these, 17,229 units are Phanerozoic in age. For convenience, and to allow precise reproduction of these results, a concise raw data table is provided in Dataset S3. North American marine macrostratigraphy has been demonstrated as representative of a global signal (68), although it also contains the signature of regional tectonics unique to North America. The macrostratigraphic dataset also includes carbonate sediments, which are more abundant and potentially more globally representative in the Paleozoic than in the Meso-Cenozoic (26). Accordingly, analyses and results of Meso-Cenozoic carbonate stratigraphy should be considered with caution.

Another biodiversity dataset compilation (12) has also been analyzed. This dataset is a range-through diversity estimate based on the times of first and last occurrences estimated for each genus. It is composed of 17,797 significant genera with robust ages, and excludes the uncertain taxa with single appearances to reduce noise, and is called the highly resolved biodiversity data (12) or “cleaned” data. The raw range-through biodiversity data composed of 36,380 genera in Sepkoski’s Compendium (2) are also analyzed and compared to the “cleaned” data (12). The overall shape of the raw and highly resolved biodiversity data is similar to those previously published for fossil families and genera. Additionally, shorter-lived genera (those that endured for 45 Myr or less), consisting of 13,682 genera (12), were also analyzed (SI Appendix, Figs. S18–S28).

The eustatic sea-level curve for the Phanerozoic was constructed based on sequence-stratigraphic studies of the marine sedimentary sections around the world (69). Correlations are based on the available marine biostratigraphic markers, and if a sea-level change can be demonstrated to be present in a number of non-contiguous sections, it is considered widespread and therefore eustatic. The details of these methodologies were previously described for the Meso-Cenozoic (69, 70) and for the Paleozoic (71). The original Ceno-Mesozoic eustatic curve (69) was revised through a comparison with the Arabian Platform records (72), and the ages were updated according to the latest versions of the time scales (35).

The tectonic data were inferred from the most recent global plate model (33) that is updated from an earlier model (73). This model additionally includes plate deformation along major rifts and orogens since rifting of Pangea started in the Early Triassic (240 Ma). It also considers the evolution of the age–area distribution of the ocean floor. This model provides plate subduction data and the length of the MOR, as well as the total variations in the crustal production rate through time. It also provides the global median trench migration data through time for all subduction zones. Changes in trench motion are related to changes in the net rotation of the lithosphere relative to the mantle and in the interaction between subducting slabs and the upper mantle/transition zone. Comparison of the earlier and recent tectonic model data (33, 73) as well as testing their cyclic character are provided in SI Appendix, Fig. S34.

Long-term irregular trends in the geological datasets were measured and removed using the loess weighted average (74) and the singular spectrum analysis (SSA) (75). The goal of the use of different methods of detrending is to test the effect of the detrend on the spectral detection of the targeted cycles.

Then, short-term periodicities were quantified via the multi-taper spectral method (MTM) associated with the harmonic F-test (76). The F-test was used to seek evidence for individual lines within frequency bands with elevated amplitude that could be related to specific significant cyclicities. The geological datasets within the studied 250 Ma interval contain irregular long-term trends with greater amplitudes (Fig. 2), inhibiting in some cases the statistical F-test significance, as the F-test is sensitive to the amplitude of the signal. The SSA method separates signals from noise in a sequence of signal components (reconstructed components, RCs) that are statistically independent, at zero lag, and based on signal strength (variance). The SSA allows to extract both cyclic and acyclic components in a signal. Therefore, we first applied the SSA in order to extract the most important (cyclic and acyclic) variabilities at tens of Myr timescales. This would significantly reduce noise. Then, we performed the F-test on the main SSA components (RC2+RC3) depicting variabilities at tens of Myr to see how regular they are. This method allows a robust identification of the periodic components in the analyzed geological datasets.

Results of the MTM and the associated F-test were then compared to the basic Fourier transform squared Periodogram (unsmoothed periodogram, Matlab’s periodogram.m) along with the classical autoregressive AR (1) red noise model (77). We have shown that application of the F-test to the studied geological datasets is more conservative than the classical autoregressive red-noise model (SI Appendix, Figs. S6 and S7). Thus, the red-noise model is supplemented when the F-test fails in the detection of targeted spectral peaks. We also use the Lomb–Scargle Periodogram (78, 79) for comparison with the MTM outputs. Such comparison was focused on cyclicities in the biodiversity data because they are widely debated (SI Appendix, Fig. S16).

Filtering is used to extract the significant cycles, based on the Gaussian filter (80). For cross-spectral analysis and correlations, we used the cross-MTM method provided in the Matlab routine (81). Cross-spectral analysis is first applied to the detrended data (Fig. 3) and then to the filtered data (SI Appendix, Fig. S35). The existence of unequal cyclic and acyclic terms among the studied geological proxies affects the coherence and phase of the shared cycles between them. Also, instable phase along the studied 0 to 250 Myr time interval should reduce the coherence between the cyclic signals. Thus, to minimize such effect, we isolate the targeted cycle, i.e., the 36 ±1 Myr cycle band, and subdivide the 0 to 250 Myr interval when applicable (SI Appendix, Fig. S35). We perform amplitude modulation (AM) analysis to extract low-frequency AM envelopes of the 28 and 36 Myr cycle bands using first bandpass Gaussian filtering and then the Hilbert transform. Gaussian bandpass filtering is paired with the SSA since one of the important features of SSA is that the SSA extracted signals can be amplitude and phase modulated (75, 82). Thus, the SSA method was used in a complementary way to AM analysis (i.e., Gauss–Hilbert).

We have also developed programs for frequency analysis and decomposition and their significance tests. We compute the PDCP (83), which searches for periodicities without assumptions on the shape of the periodic signal, through a measurement of statistical independence between the data and an angle varying with a certain frequency. This allows us to search for periodic signals that are not necessarily sinusoidal. Once the PDCP is computed, the data are phase-folded at the detected period, and the shape of the periodic signals is approximated with a fitted low-order polynomial with a constraint guaranteeing that it is derivable. The estimated periodic signal is then subtracted from the data, and the PDCP is applied to the residuals; this procedure is repeated in an iterative manner. Results are summarized in SI Appendix, Fig. S18, which show significant 26.3 to 27.5 Myr and 37.5 Myr periodicities over the Ceno-Mesozoic interval (Table 1).

We further explore the fitting of non–strictly periodic signals through orthonormal polynomials and Gaussian processes (84). Gaussian processes are stochastic such that any set of n values of the process at any n times follow a Gaussian distribution. Such processes are characterized by their autocovariance function or kernel: a function of two times t and t’ giving the covariance of the value of the process at t and t’. This kernel represents the self-similarity of the process. We represent the data as a sum of a polynomial of time of order three and Gaussian processes: one white noise with a free amplitude, one with an exponential decay (EXP) with free amplitude and timescale, and several SHO kernels, characterized by their amplitude, fundamental frequencies, and quality factor (85). The quality factor of the oscillator is a measure of the loss of its energy over time. SHO Gaussian processes have the same auto covariance as a harmonic oscillator excited by a random Gaussian noise. As such, it is very well suited to represent phenomena that have a certain cyclicity but not a strict periodicity: phase and amplitude might locally vary but on average the power spectrum density will have a concentration of power close to a certain frequency. The higher the quality factor is, the closer is the signal to a pure sinusoid. The rationale of using such processes here is to give the opportunity to the model to choose a non–strictly periodic behavior (pure uncorrelated noise, correlated noise with exponential time scale that is not periodic, or non–strictly periodic signals). This is appropriate to geological cycles, which are not purely periodic but exhibit generally a quasi-periodic character, due to either the forcing process or the nonlinear response of the geological process (36, 38).

All our models include a noncorrelated homoskedastic Gaussian noise, EXP Gaussian process, and we consider models with 1, 2, 3, 4, and 5 SHO Gaussian processes. We model them using the S+LEAF public software (86). We adopt log-uniform prior distributions on the parameters and compute the Bayesian evidence of each model with the PolyChord software (87, 88). We choose a uniform distribution on the value of the SD of all the processes involved (white, EXP, and SHO) from 0 to 186, the latter value would correspond to the data being simply the sum of a third-order polynomial and a single Gaussian process. The prior on the EXP timescale is log uniform from 0.5 to 500 Myr; the prior on the quality factor and period of the SHO kernels are, respectively, log uniform from 1 to 1,000 and log uniform from 1.5 to 1,000 Myr.

For 1, 2, 3, 4, and 5 SHO Gaussian processes, we find log evidence of −2968.13 ± 0.24, −2966.45 ± 0.22, −2966.45 ± 0.23, −2965.59 ± 0.24, and −2966.51 ± 0.22, respectively. The highest evidence is therefore the four SHO models, but the different models do not have very different Bayesian results. We want to detect signals with a certain accuracy value on the period of their cyclicity. It has been shown (89) that in that case, the optimal detection criterion consists in computing the posterior probability to have a signal within a certain period range. For the 1, 2, 3, 4, and 5 SHO Gaussian processes, we find that the probability to have a signal with an SHO period between 25 and 40 Myr is 18, 59, 99.5, 99.4%, and 99.6%, respectively. Weighted by the posterior probability to have 1, 2, 3, 4, and 5 SHO Gaussian processes, we find an overall probability of 90%.

To illustrate our results, for different models, we look at the vector of parameters with the highest likelihood and represent the corresponding predictive distributions of the different Gaussian processes. We find that in all the models, except the 1 SHO kernel one, the parameter vector with the highest likelihood always includes an SHO process with a 28 to 30 Myr period (for 1 SHO kernel, we find 45 Myr, but this model is less favored). This process may have power at different frequencies, and when we compute the PDC of the corresponding predictive mean, we find a 36.6 Myr periodicity in each case when we consider the Ceno-Mesozoic data only (last 250 Myr) and 37.89 Myr when considering all of the Phanerozoic data. As an example, we show in SI Appendix, Fig. S19, the values of the parameters of the processes with the highest likelihood for the three and four SHO kernels.

S.B. designed research; S.B. performed research; S.B., S.E.P., R.D.M., B.U.H., and N.H. analyzed data; and S.B., S.E.P., R.D.M., and B.U.H. wrote the paper.

The authors declare no competing interest.