Thresholds of temperature change for mass extinctions
Paleotemperature estimatorsPhanerozoic paleotemperature data are derived from oxygen isotope (δ18O), clumped isotope (Δ47), and organic geochemistry (TEX86) data. Many published studies have used different conversion formulae to get from measured values to temperatures for each paleothermometry method, which makes it difficult to compare results from various studies. Because of this, we used a uniform formula to re-calculate paleotemperatures. Different estimators generally show consistent trends for a given time interval, e.g. Cretaceous planktonic foraminifera δ18O and TEX8640, and Jurassic belemnite δ18O41.Most Paleozoic and Mesozoic paleotemperature data come from oxygen isotope paleothermometry. Cenozoic paleotemperature data derive mainly from TEX86 and oxygen isotope of planktonic foraminifera. At present, there are major debates over the first-order trend of Phanerozoic surface temperatures calculated from oxygen isotopes12,42,43,44. For instance, it is uncertain whether δ18Osw evolved towards more enriched values due to exchange at mid-ocean ridges, potentially impacting interpretations of seawater temperatures in the early Paleozoic, for instance43. However, this is not an issue in our study because we focus only on changes in temperature rather than absolute values.The data we used for each time bin was selected from a large paleotemperature dataset12 and published literature by adopting the following criteria: 1) data with well-constrained ages; 2) data measured at high time resolution; 3) data from tropical/subtropical regions. For multiple sets of data in the same time bin, we use the average values. Where possible, in order to test the reliability of climate events, we use temperature data from two different proxies or two different regions. For example, δ18O of planktonic foraminifera and TEX86 were used for the late Eocene (Pg4) and Maastrichtian (K8) intervals; and conodont δ18O from South China and Armenia were used for the Lopingian (P5) interval. The results show that data from two different proxies or different regions have a similar magnitude and rate of temperature change (Supplementary Fig. 4).Oxygen isotopesOxygen isotope values from carbonate fossils (i.e., planktonic foraminifera, brachiopod, oyster, and belemnite) are converted to SSTs using the BAYFOX Bayesian model (https://github.com/jesstierney/bayfoxm)45. For phosphate fossil conodont, we used Monte Carlo simulations to propagate parameter uncertainties to temperature estimation by using the equation of Pucéat et al.46:$${SST}=118.7(pm4.9)-4.22(pm0.20)times left({delta }^{18}{O}_{{con}}-{delta }^{18}{O}_{{SW}}right)times 1000$$
(1)
where δ18Ocon is the oxygen isotope composition of conodonts, and δ18OSW is the oxygen isotope composition of seawater.Seawater oxygen isotope composition is affected by salinity and changes in continental ice volume43,47. As such, δ18O data from localities with abnormal salinities (e.g., evaporite facies, upwelling regions) were excluded (see paleotemperature estimators above). Changes in continental ice volume were also considered in our calculation of paleotemperatures. To do this, the oxygen isotope value of seawater was set to −1‰ (VSMOW) for ice-free time intervals48 and +1‰ (VSMOW) during glacial maximum intervals, e.g., the Pennsylvanian Glacial Maximum and the Pleistocene Last Glacial Maximum49. Most oxygen isotope data are from tropical and subtropical regions (Supplementary Fig. 5), therefore, no latitudinal seawater corrections were made. In addition, seawater pH has a further important influence on the δ18O of foraminifera50. During global warming events that are related to CO2 release, seawater pH can decrease, which would lead to underestimated SSTs50,51. Unfortunately, pH estimates for the Phanerozoic time bins are few, and their values are highly uncertain52. Therefore, we did not pH-correct δ18O estimates. Oxygen isotope values that were likely affected by diagenetic alteration were also removed from the database. Diagenetic screening criteria included δ18O values that are unrealistically negative or positive, and δ18O values from carbonate fossil shells with [Mn] > 250 ppm and [Sr] 0.4 were not used in this study. Other screening criteria of TEX86 were also employed. Notably, data were removed from the database if Methane Index (MI) >0.5 (ref. 58), delta-Ring Index (ΔRI) > 0.3 (ref. 59), %GDGT-0 >67% (ref. 60), and/or fCren′:Cren′ + Cren >0.25 (ref. 40).Age modelGeochronologic constraints can provide absolute age at the initiation and termination of warming/cooling events. In our database, only one warming event, at the Permian-Triassic boundary, meets this criterion. Other dates were obtained using a comprehensive approach including isotope geochronology, astrochronology, and biostratigraphy with reference to the Geologic Time Scale 201261. Age data were applied in the following order of priority: isotope geochronologic age, astrochronologic age, and biostratigraphic age. If isotope geochronologic ages were available, we gave priority to these absolute dates. In the absence of absolute age data, astrochronologic ages were preferred. If neither of these numerical data was available, the climatic events were timed based on the age of biozones in the Geologic Time Scale 201261. Relative ages within individual biozones were constrained based on the stratigraphic position. Age uncertainty was calculated for all events by using Monte Carlo simulation based on the assumption of uniform distribution of ages.The magnitude and rate of temperature changeThe maximum magnitudes (ΔT) of temperature change within each of the 45 defined time intervals were calculated as:$$Delta T={T}_{1}-{T}_{0}$$
(2)
where ({T}_{0}) and ({T}_{1}) represent the initial and terminal temperature of warming/cooling events, respectively. Supplementary Fig. 1 illustrates how the parameters ({T}_{0}) and ({T}_{1}) are derived from a climatic time series (t).The time scale Δt represents the duration of time during which the temperature increases/decreases. Δt was calculated as:$$Delta t={t}_{1}-{t}_{0}$$
(3)
where ({t}_{0}) and ({t}_{1}) represent the initial and terminal time of warming/cooling events, respectively.The rate (R) of temperature change was computed as:$$R=Delta T/Delta t$$
(4)
The ratio (R) represents the rate of a single warming/cooling event over geological time.ΔT, R and their uncertainties were calculated by using Monte Carlo simulation based on the distribution of data from two groups around t0 and t1.Temperature databaseThe temperature database is composed of the most significant warming/cooling events in 45 time intervals from the Late Ordovician (445 Ma) to early Miocene (16 Ma) (Supplementary Data 162). The time intervals are consistent with the time bins used to compute biodiversity and evolutionary rates63,64, and are defined by one or several neighboring geologic stages with roughly uniform durations (averaging 9.71 Myr). Time bins range between 4.7 and 18.9 Myr. The 45 bins are Katian and Hirnantian (Or5), Llandovery (S1), Wenlock (S2), Ludlow and Pridoli (S3), Lochkovian and Pragian (D1), Emsian (D2), Eifelian and Givetian (D3), Frasnian (D4), Famennian (D5), Toutnaisian (C1), early Visean (C2), late Visean and Serpukhovian (C3), Bashkirian (C4), Moscovian and Kasimovian (C5), Gzhelian (C6), Asselian and Sakmarian (P1), Artinskian (P2), Kungurian and Roadian (P3), Wordian and Capitanian (P4), Wuchiapingian and Changhsingian (P5), Induan and Olenekian (T1), Anisian and Ladinian (T2), Carnian (T3), Norian (T4), Rhaetian (T5), Hettangian and Sinemurian (J1), Pliensbachian (J2), Toarcian and Aalenian (J3), Bajocian-Callovian (J4), Oxfordian (J5), Kimmeridgian and Tithonian (J6), Berriasian and Valanginian (K1), Hauterivian and Barremian (K2), Aptian (K3), Albian (K4), Cenomanian (K5), Turonian-Santonian (K6), Campanian (K7), Maastrichtian (K8), Paleocene (Pg1), Ypresian (Early Eocene, Pg2), Lutetian (early Middle Eocene, Pg3), Bartonian and Priabonian (late Middle-Late Eocene, Pg4), Oligocene (Pg5), and Aquitanian and Burdigalian (Early Miocene, Ng1).Original proxy data, calculated temperature data, geologic age, time span, paleolatitudes, and relevant references are provided in Supplementary Methods and Source Data. Most paleotemperature data are sea surface temperatures (SST) from tropical and subtropical regions (between 40°N and 40°S, Supplementary Fig. 5). Only one collection is from a mid-latitude region with a paleolatitude of 42.71°N (Source Data). Paleolatitudes were reconstructed using PointTracker v7 rotation files published by the PALEOMAP Project65.For most time bins, the major climate change events (warming/cooling) were restricted to that bin. A few events occurred that cross two adjacent bins, e.g., the rapid climate warming around the Permian-Triassic and the Triassic-Jurassic boundaries17,66,67. Both these events had relatively short durations (61–400 kyr) and were initiated towards the end of the first (earliest) time bin. Therefore, these warmings would have primarily impacted organisms in the first bin too. Accordingly, we didn’t divide the warming event into two intervals belonging to the two adjacent bins, but take instead the warming event as belonging to the first bin.Extinction rateThere are a few well-established methods to calculate extinction rate in geological history, e.g., boundary-crosser, gap-filler (GF), and three-timer (3 T) rates63,68,69. Both the GF and 3 T rates are calculated from occurrence data and have higher accuracy than the boundary-crosser method64. The boundary-crosser method is more susceptible to major biases such as the Signor-Lipps effect and sampling bias because it uses full age ranges instead of investigating sampling patterns across limited temporal windows64. The 3 T method can be noisy in the case of high turnover rates and poor sampling. The GF method is more precise when sampling is very poor.Gap-filler and three-timer extinction rates of marine animals were calculated using data from the Paleobiology Database (PBDB, http://paleobiodb.org), which was downloaded on 4 January 2021. The fossil dataset includes all metazoans except for Arachnida, Insecta, Ostracoda, and Tetrapoda and consists of 850,840 fossil occurrences of 37,134 genera. These four groups that are excluded (in keeping with some previous studies64,69) were not used because many of them (Arachnida, Insecta, and Tetrapoda) are terrestrial and appear in marine strata. Ostracoda also has a record in terrestrial rocks. Similar to most studies using data from PBDB, we use genus-level occurrences because genus-level taxonomy is better standardized than the taxonomy of species. GF extinction rate (({r}_{{GF}})) is calculated by using the equation:$${r}_{{GF}}={{{{{rm{log }}}}}}left(frac{2T+{pT}}{3T+{pT}+{GF}}right)$$
(5)
Three-timer extinction rate (({r}_{3T})) is calculated by using the equation:$${r}_{3T}={{{{{rm{log }}}}}}left(frac{2{T}_{i}}{3T}right)+{{{{{rm{log }}}}}}left(frac{3T}{3T+{pT}}right)$$
(6)
where ({GF}) represents gap-fillers, i.e., the number of genera in time bins (i-1) and (i+2) and not within bin (i+1); (2{T}_{i}) represents the number of genera sampled before and within the ith time bin; (3T) represents the number of genera sampled in three consecutive time bins; and ({pT}) represents the number of genera sampled before and after but not within a time bin.Autocorrelation analysisWe tested for serial correlation in the time series of extinction rate, ∆T, and log R, because autocorrelation could impact the calculated correlation coefficients between these data. The autocorrelation functions (ACFs) of extinction rate, ∆T, and log R are very low at all lags >1 (Supplementary Fig. 6), indicating that all-time series lack significant serial correlation.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
