Episodes of extinction of families and genera of marine organisms have been reported to show an underlying cycle of 26.4 to 27.3 My over the past 260 My (Raup and Sepkoski 1984, 1986; Rampino and Caldeira 2015) although extinction periodicity is still a contentious issue (Omerbashich 2006; Melott and Bambach 2014, 2017; Rampino and Caldeira 2015; Erlykin et al. 2017). Estimates of diversity in the well-sampled record of calcareous plankton for the past 230 My show a similar 26 to the 29-My period (Prokoph et al. 2004). In order to test for potential cycles in parallel extinctions of non-marine tetrapods (amphibians, reptiles, birds and mammals) over a similar interval, we compiled an updated record of distinct tetrapod-extinction episodes from the relevant literature. These are brief intervals of intensified extinctions that apparently occurred close to chronostratigraphic stage boundaries, which for the last 300 My are tightly dated, with error bars ≤ ± 1 My in most cases (Table 1) (Ogg et al. 2016).

We chose to expand on previous work on identifying distinct extinction episodes (Benton 1985, 1989, 1995, 2010; Maxwell and Benton 1990), rather than using data on total diversity or rates of extinction and origination from current databases (Alroy 2008, 2014). Use of the ages of discrete recognised extinction events reduces problems related to uncertainties in stage duration (that affect the calculation of extinction and origination rates), avoids questions concerning completeness of the fossil record, and minimises sampling problems (Melott 2008) that can affect estimates of total diversity (Foote 1994). This is the same strategy taken in prior searches for periodicity in marine-extinction events (Raup and Sepkoski 1984, 1986; Rampino and Caldeira 2015). There are quantitative estimates of % extinction of families for six of the non-marine extinction episodes, and we relied on semi-quantitative information as to the severity of the four other recognised events (Benton 1985, 1989, 2010; Maxwell and Benton 1990). Our search in the literature turned up no other events that could be interpreted as non-marine extinction episodes. We estimate that we could detect extinction pulses with at least ~10% loss of tetrapod families (Benton 1989).

The known stratigraphic ranges of families of non-marine tetrapods were initially compiled by Benton (1985) who recognised 6 distinct episodes of intensified extinctions since the origin of tetrapods in the Late Devonian: the early Permian (290.1 Ma–the Sakmarian-Artinskian Stage boundary–58% tetrapod family extinctions), the latest Permian (251.9 Ma–49% extinction), the latest Triassic (201.4 Ma–22% extinction), the latest Cretaceous (66 Ma–14% extinction), the early Oligocene (33.9 Ma–8% extinction) and the late Miocene (7.25 Ma – 2% extinction). In more recent studies, the latest Jurassic (145.7 ± 0.8 Ma), also shown as a distinct non-marine extinction episode by Benton (1989), (1995), (2010); Maxwell and Benton (1990), is confirmed as a significant event (Tennant et al. 2017). Benton later identified an additional non-marine extinction event at the end of the Late Triassic Carnian Stage (Benton 1986; Miller et al. 2017) dated at 228.5 ± 0.25 Ma (Ogg et al. 2016) (Table 1). A slightly older (~232 Ma) mid-Carnian age would put the extinction event at the time of the ‘Carnian Pluvial Event’ (Ogg et al. 2016).

Some corrections in Benton’s age data have been applied. Benton (1985) listed the latest Triassic tetrapod extinction as occurring at the end of the Norian Stage (209.5 or 206.5 Ma), but this is a result of problems in defining the final stage of the Triassic at the time. This non-marine extinction event clearly occurred at the end of the latest Triassic Rhaetian Stage at 201.4 ± 0.17 Ma (Hesselbo et al. 2001; Lucas and Tanner 2015; Ogg et al. 2016) coincident with the latest Triassic marine-mass extinction (Table 1) (Raup and Sepkoski 1984, 1986). Benton also placed the mid-Cenozoic tetrapod extinction event (the ‘Grande Coupure’) in the Rupelian Stage (early Oligocene, 28.1 Ma) (Benton 1985). It is clear now that in Europe and Asia, the tetrapod extinctions actually took place close to the Eocene-Oligocene boundary at 33.9 Ma (Sun et al. 2014) (Table 1).

In the Early Permian, the Sakmarian-Artinskian boundary event is now dated at 290.1 ± 0.2 Ma (Chumashov et al. 2013; Ogg et al. 2016) (Table 1). Another proposed Early Permian event, Olson’s Extinction (at ~270 Ma) is now considered to be an artefact of a geographic shift in the sampling of the tetrapod fossil record (Benton 2012; Benson and Upchurch 2013) and is not included in our extinction dataset. The other Late Palaeozoic non-marine tetrapod extinction events apparently lack any such substantial sampling bias (Sahney and Benton 2008; Benson and Upchurch 2013).

Based on subsequent work, two additional apparent non-marine tetrapod extinction episodes have come to light, one at the mid-late Norian boundary (215 ± 0.25 Ma) especially in western North America (the boundary between the Adamanian and Revueltian Land Vertebrate Faunachrons) (Parker and Martz 2011; Atchley et al. 2013), which is marked also by marine extinctions (Onoue et al. 2016). Another, now well-documented, severe tetrapod extinction event apparently occurred at the end of the Guadalupian (Capitanian) Stage (259.8 ± 0.5 Ma) (Retallack et al. 2006; Day et al. 2015; Metcalfe et al. 2015; Ogg et al. 2016) correlative with a strong marine-extinction episode (Rampino and Shen 2019; Rampino et al. 2019) (Table 1).

In related work, Benton (1989) showed six major contractions of the ecospace occupied by non-marine tetrapods at the end-Guadalupian, end-Permian, end-Triassic, end-Jurassic, end-Cretaceous and the Eocene-Oligocene boundary, matching six times of recognised tetrapod-extinction pulses, and supporting the idea that these were times of species loss and not merely related to decreases in originations.

Thus, we recognised a total of at least 10 events of intensified extinctions of non-marine tetrapods over the last 290 My, and note a one-to-one correlation of 8 out of the 10 tetrapod-extinction episodes with reported times of coeval marine-extinction events (Raup and Sepkoski 1984, 1986; Rampino et al. 2019) (Table 1), which suggests a common cause, and perhaps a common 26 to 27-My periodicity.

We tested for periodicity in the ages of the non-marine tetrapod extinction events first by applying a relevant circular spectral method developed by Stothers (1991) and Lutz (1985), which was designed to test for cycles in a time series of discrete events without the use of amplitude information. This method has been used previously to search for cycles in records of the ages of impact craters and marine-extinction events (Rampino and Caldeira 2015).

The time series consists of Dirac delta-functions, where only the times of the extinction events (which have small calculated errors of ≤ ± 0.8 My) and not the magnitudes, are important. In the time series, some of the events may represent an underlying periodic process with some superimposed noise, but some of the events could be unrelated to the periodic process (Stothers 1991).

This method effectively measures the degree of clustering of events when the dates of the events are mapped to a data modulus of the test period under consideration. It produces a normalised distance metric, which should be minimised at the preferred cycle. The method works well for time-series that lack accurate amplitude information, may be a mixture of periodic and non-periodic events, with potential missing events, and consist of unevenly spaced data (Stothers 1991). Furthermore, in our statistical methods, any small dating errors arising from random or systematic errors in the geologic timescale, or any missing events do not significantly affect the value of the periods detected, although they may affect the relative heights of spectral peaks (Lutz 1985; Raup and Sepkoski 1986; Stothers 1989, 1991).

For circular spectral analysis, a timeline can be ‘wrapped’ around a circle, the circumference of which represents a trial period (Lutz 1985; Stothers 1991). For each occurrence, we calculate a unit vector from the origin. A series that is not periodic will tend to plot randomly around the circle, regardless of the trial period selected. A periodic series, however, will tend to form a cluster at one point on the circumference when the correct trial period P is selected. The angular location of the cluster relative to 0° (the present) gives the phase (t ₀).

The event times t_i are mapped onto a circle by conversion to angles a_i and b_i : (1) $\begin{array}{c}{a}_i=sin2\pi /P\left(t_i\right)\\ b_i=cos2\pi /P\left(t_i\right)\\ S=1/N{\sum ^N}_{i=1}\left(a_i\right)\\ \begin{array}{c}C=1/N{\sum ^N}_{i=1}\left(b_i\right)\\ R={\left(S^2+C^2\right)}^{1/2}\end{array}\end{array}$(1)

where i ranges from 1 to N (the number of events), and P is the trial period. S and C are the summations, S = (∑sin  a_i)/N and C = (∑cos  b_i)/N.

Application of circular statistics leads to a mean vector magnitude, R = (S ² + C ²)^1/2 (a normalised measure of goodness of fit). The direction of the vector that maximises R, and therefore minimises the dispersion, at the trial period P indicates the phase, which can be computed: t ₀ = (P/2π) tan⁻¹ (S/C), if C > 0; or as t ₀ = (P/2) + (P/2π) tan⁻¹ (S/C), if C < 0. When R is plotted against P, then maximal values of R would correspond to periods in the series, t ₁, t ₂, t_N . If, however, t ₁, t ₂, …,t_N are randomly distributed, then (a_i, b_i ) would define a random walk, and the sum R will be small (Stothers 1991).

In this method, the normalised distance metric does not have a direct interpretation in terms of statistical significance (Lutz 1985; Stothers 1991). To assess statistical significance, we address the question: What is the likelihood that a time series drawn at random from the universe of possible similar time series would have a smaller normalised distance metric at the test period than does the actual time series?

To ensure that our results are a consequence of true periodicity, we constructed the set of test time series as follows: First, we calculated the length of time between each of the N events in the test time series, giving N – 1 intervals. We then took a random permutation of these N – 1 Intervals and generated a time series of N events by taking the cumulative sum of this random permutation. Thus, we are testing against a set of time series that contain the same set of intervals as the test time series but permuted in a different order. This makes it clear that the periodicity detected is a consequence of the ordering of the intervals and not the statistical distributions of the intervals themselves.

We compared the Stothers (1991) normalised distance metric calculated at each test period with the same metric computed for 100,000 pseudo-time series generated by the permutation approach described above, each having 10 pseudo-events over an interval of 290 My. This significance test is powerful because the same set of intervals is present in the test time series and every one of the pseudo-time series, and the total length of all of the pseudo-time series is identical to that of the test time series, eliminating artefacts associated with varying record lengths (Lutz 1985)^. The significance test thus ascertains the likelihood that a random permutation of the same set of intervals would produce as strong evidence of periodicity as the time series being tested.

The power spectrum produced by our circular spectral analysis for a range of trial periods running from 5 to 50 My is shown in Figure 1. A number of high-frequency peaks are present, which may be a result of artefacts in the extinction dataset. The highest spectral peak between 15 and 50 My, however, appears at a period of 27.3 My. The 0.95, 0.99, and 0.999 confidence lines in Figure 1 represent the 95^th, 99^th and 99.9^th percentile results of 100,000 random pseudo-data sets. The 27.3-My spectral peak is significant at the 95% confidence level. Two other relatively strong low-frequency spectral peaks occur at ~32 My and 36.9 My, which may be the 6/5 and 4/3 harmonics of the 27.3-My cycle, or an indication that the 27.3-My underlying period is somewhat irregular (Lutz 1985; Stothers 1991). The spectral analysis methods used in this study may commonly produce relatively high spectral peaks that are harmonics of the main period (Lutz 1985; Stothers 1989, 1991).

P-values in Figure 2 represent the number of cases out of 100,000 Monte Carlo simulations that had higher scores in our statistical method. For a period of 27.3 My, a P-value of 0.05 means that of 100,000 Monte Carlo simulations only 5,000 cases had higher scores in our method of time-series analysis. Again, two additional low-frequency spectral peaks occur at 32 and 36.9 My, which may be harmonics of the 27.3 My cycle.

To check our results with the circular method, we also performed more traditional Fourier transform analysis (Gasquet and Witomski 1999) of the ages of the 10 recognised non-marine extinction episodes of the last 290 My. This is similar to searches for periodicity in marine-extinction events, using parametric and non-parametric methods of time-series analysis (Raup and Sepkoski 1984). For the Fourier method, we first rounded the original data to the nearest million years, then utilised a standard Tukey window with a window size of 6 My. In the end, a Fourier transform was applied to the processed data (Gasquet and Witomski 1999). The highest spectral peak occurs at a period of 27.5 My, essentially the same as for the circular spectral method (Figure 3).

We computed the significance of the result by generating 100,000 test datasets and comparing the spectral power at 27.5 Myr (Figure 3). To simulate the test datasets, the intervals between every two consecutive events were calculated. After permuting the order of the intervals, a new time series was re-calculated based on the new set of intervals. This method ensures that the test time series are between 0 and 260 My and the number of events remains at 10. Ultimately, less than ~1% of the test datasets produced a higher spectral power at 27.5 My, indicating a confidence level of at least 99% for that period (Figure 3).

The detection of a similar strong 27.3 to 27.5 My spectral peak in the non-marine extinctions, utilising the two different spectral-analysis techniques, supports the reality of the underlying ~27.5 My periodicity. We also performed Fourier spectral analysis (Gasquet and Witomski 1999) of the combined age data for non-marine and marine-extinction events (a total of 23 events in the last 290 My) (Table 1) resulting in a spectrum with a single strong spectral peak at 27.5 My (Figure 4) suggesting a single 27.5 My periodicity in the marine and non-marine data sets.

Considering extinctions over a longer time interval, Rampino and Haggerty (1996) performed Fourier transform analysis of 21 marine-extinction peaks back to 541 Ma and found a similar significant spectral peak at 27.3 My. They also performed Fourier analyses on a series of truncated marine-extinction time series starting from 0 to 541 Ma, and subtracting one extinction event at a time. The truncations exhibited a stable spectral peak between 26.5 and 27.3 My, which remained the dominant feature in the spectra. In a number of the truncation analyses, the second-highest or third highest spectral peak in the marine extinctions occurred at ~36 My (Rampino and Haggerty 1996) a potential harmonic similar to that which appeared in our results for extinction events of non-marine tetrapods (Figures 1 and 2).

Rampino and Stothers (1998) performed linear spectral analysis (Stothers 1991) on the 260-My record of marine-extinction events and found that the highest spectral peak occurred at ~28 My, with two somewhat smaller low-frequency peaks at ~32 and 36 My, again similar to the potential harmonic spectral peaks for non-marine extinctions (Figures 1 and 2), suggesting a common structure in the marine and non-marine time series analysed.

In a related study, Stothers (1989) found a weak 28-My period in the ages of Mesozoic and Cenozoic stratigraphic stage boundaries, which might be expected, as the stage boundaries are commonly defined by faunal changes; 13 of the 39 stage boundaries (33%) are coincident with marine-extinction events, and 10 of the stage boundaries (25%) are marked by non-marine extinction pulses, both of which show the ~26 to 27-My period. To follow up on Stothers (1989) analysis, but using the latest data, we performed a circular spectral analysis (Stothers 1991) of the ages of an updated list of Mesozoic/Cenozoic stage boundaries as re-dated by Ogg et al. (2016) and found a weak 27-My period, as expected.

The correlations and similar cycles in marine and non-marine extinction episodes suggest a common cause. The eight recognised co-occurrences of non-marine and marine-extinction events are in agreement with the ages of Large Igneous Province eruptions, and five of those events (which are all continental flood basalts) are marked by stratigraphic volcanogenic Hg anomalies, further suggesting a close relationship between flood-basalt volcanism and the extinction pulses (Table 1) (Rampino et al. 2019).

Furthermore, however, three of the co-extinction episodes (at 66, 145 and 215 Ma) are also coeval with the three largest recorded impacts (craters ≥100 km in diameter) of the last 260 My (Rampino et al. 2019) (Table 1), which are apparently capable of causing significant extinction events (Rampino 2020). This suggests a connection among the three phenomena (extinctions, LIP eruptions and large impacts). For example, we note that the end of the Cretaceous (66 Ma) (Table 1) is marked by the Deccan Basalts in India, and the very large (180-km diameter) Chicxulub impact in the Yucatån. The end-Jurassic extinction event (145.7 ± 0.8 Ma) is correlative with the formation of the large submarine Tamu Massif on the Shatsky Rise oceanic plateau in the northwest Pacific (Sager et al. 2013) (Table 1), dated by the ⁴⁰Ar/³⁹Ar method at 144.4 ± 1.0 Ma (Geldmacher et al. 2014). The end-Jurassic is also marked by the large (~130-km diameter) South African Morokweng crater (145.2 ± 0.8 Ma) (Rampino et al. 2019).

A third extinction event, the mid-Norian marine/non-marine extinctions (215 ± 0.25 Ma) is close in age to the large (100-km diameter) Manicouagan impact in Quebec (214.6 ± 0.05 Ma), and may also be correlative with the less well-known Angayucham Basalts in the Brooks Range in Alaska (Rampino et al. 2019), a former oceanic plateau roughly dated at 214 ± 7 Ma (Table 1). These co-occurrences suggest multiple causes for some of the extinction episodes (Keller 2008).

LIP eruptions are capable of producing severe environmental conditions including brief periods of intensely cold climate from atmospheric aerosols (Self et al. 2006), acid rain, ozone destruction and increased UV-B radiation (Black et al. 2014; Benca et al. 2018). Longer-term volcanic effects include lethal greenhouse heating (Sun et al. 2012; Benton 2018) from the release of CO₂ and CH₄, ocean acidification and anoxia (Rampino et al. 2019), and potential release of poisonous H₂S from the anoxic oceans (Kump et al. 2005). On the other hand, large impacts (with craters ≥100 km in diameter) can apparently create widespread dark and cold conditions, wildfires, acid rain, ozone depletion and other effects that would stress non-marine as well as marine organisms (Toon et al. 2016; Rampino 2020).

The potential correlation of the periodic non-marine and marine extinctions with large LIP eruptions and large-body impacts is noteworthy, as both the eruptions and impact craters have also shown some evidence of an underlying ~30-My cycle (Rampino and Caldeira 1993, 2015). Although it is likely that LIP eruptions are the result of internal Earth dynamics (Rampino and Caldeira 1993) the pacing may be partly related to astrophysical factors, notably the ~30-My interval between passes of the solar system through the midplane of the disk-shaped Galaxy (Rampino and Stothers 1984), which might modulate encounters with disk-dark matter leading to periodic comet storms (Randall and Reece 2014). Furthermore, astrophysicists have predicted that coeval capture and annihilation of some of the disk dark-matter particles in the Earth’s interior might cause periodic thermal disturbances that could trigger concurrent mantle-plume activity and resulting flood-basalt volcanism (Abbas and Abbas 1998; Rampino 2015). It has not escaped our attention that coincidently large impacts and flood-basalt eruptions might indicate a more direct connection between the two phenomena (Rampino 1987, 2017; Richards et al. 2015). Further coupled geological (e.g., improved dating of geological events, better estimates of extinction magnitudes) and astrophysical research (e.g., searches for disk dark matter) (Shaviv et al. 2014; McKee et al. 2015; Kramer and Randall 2016) should shed additional light on these potential associations.