Impact statement
What are the biggest consequences of climate change for marine ecosystems? Is it deoxygenation, thermal stress, ocean acidification or any combination thereof? The Permian–Triassic climate crisis was an episode of severe and rapid climate warming with similarities to the worst-case projected scenarios for the near future. To better understand which consequences of this climate event led to one of the most severe biodiversity crisis ever, we implemented a novel approach of statistically integrating high-resolution fossil data with high-resolution geochemical data. Our results demonstrate that for equatorial, marine ecosystems, oxygen isotope (temperature proxy) and cadmium isotope (primary productivity proxy) dynamics best explain the marine extinction. This suggests that the biggest threats to past and modern biodiversity in these settings are the impacts of thermal and nutrient stress, as well as associated trophic knock-on effects.
Introduction
The most distinct and widely acknowledged causes of extinction in marine ecosystems today are pollution, habitat loss, overexploitation, introduction of invasive species and climate change (Bonebrake et al., Reference Bonebrake, Guo, Dingle, Baker, Kitching and Ashton2019; IPBES, 2019). Understanding how these threats will reduce populations or drive species to extinction is a core component of modern-day conservation and policymaking. One issue is that horizon threats, like climate change, occur on global and centennial scales that are much broader in scope than knowledge that can be acquired using modern-day datasets alone (Bonebrake et al., Reference Bonebrake, Guo, Dingle, Baker, Kitching and Ashton2019). The rock record, however, provides the only record of long-term biotic responses from disturbances and information on ecosystem re-establishment, which is now a priority in the Intergovernmental Panel on Climate Change (IPCC, 2021; Finnegan et al., Reference Finnegan, Harnik, Lockwood, Lotze, McClenachan and Kahanamoku2023; Kiessling et al., Reference Kiessling, Smith and Raja2023). We can, therefore, use different hyperthermal events of the past to provide key information on how horizon threats operate at community, ecosystem, and even biome levels.
The Permian–Triassic climate crisis is an exceptionally rapid warming event (around 8–12 °C rise in 60 ± 48 ka at low latitudes) from the latest Permian into the Early Triassic (late Griesbachian) (Joachimski et al., Reference Joachimski, Lai, Shen, Jiang, Luo, Chen, Chen and Sun2012, Reference Joachimski, Alekseev, Grigoryan and Gatovsky2020; Sun et al., Reference Sun, Joachimski, Wignall, Yan, Chen, Jiang, Wang and Lai2012; Chen et al., Reference Chen, Shen, Li, Xu, Joachimski, Bowring, Erwin, Yuan, Chen, Zhang, Wang, Cao, Zheng and Mu2016; Gliwa et al., Reference Gliwa, Wiedenbeck, Schobben, Ullmann, Kiessling, Ghaderi, Struck and Korn2022). This climate crisis is thought to have been caused by the simultaneous eruptions of the Siberian Traps Large Igneous Province and the combustion of organic-rich sedimentary rocks (Burgess and Bowring, Reference Burgess and Bowring2015), leading to a large and rapid injection of CO2 and volatiles into the atmosphere (Svensen et al., Reference Svensen, Planke, Polozov, Schmidbauer, Corfu, Podladchikov and Jamtveit2009; Joachimski et al., Reference Joachimski, Müller, Gallagher, Mathes, Chu, Mouraviev, Silantiev, Sun and Tong2022). This event is also associated with the Permian–Triassic mass extinction, the most catastrophic mass extinction on Earth, which was highly selective against taxonomic groups that dominated pre-extinction marine communities (Foster et al., Reference Foster, Ayzel, Münchmeyer, Rettelbach, Kitzmann, Isson, Mutti and Aberhan2022a, Reference Foster, Allen, Kitzmann, Münchmeyer, Rettelbach, Witts, Whittle, Larina, Clapham and Dunhill2023a), with an estimated loss of 81–96% of species (Erwin, Reference Erwin1993; Stanley, Reference Stanley2016). Multiple environmental perturbations occurred simultaneously during the climate crisis, making it difficult to disentangle which specific environmental changes were most significant in causing the extinctions. In addition, environmental stressors can interact in antagonistic or synergistic ways, where one stressor could reduce the impact of another or where a multiple of stressors can lead to an additive response (Benton, Reference Benton2018; Penn et al., Reference Penn, Deutsch, Payne and Sperling2019). Furthermore, the drivers of extinction are expected to be spatially heterogeneous, as factors such as carbonate saturation state and the polar amplification of climate warming lead to heterogeneous patterns (Feldl and Merlis, Reference Feldl and Merlis2021). Therefore, it is not unequivocally known exactly which factors played a major role in causing the biodiversity crisis. This lack of understanding is also due to the poor geographical coverage of continuous Permian–Triassic successions and the small number of sections that have been investigated at a high-resolution with multiple proxies for environmental and biodiversity changes.
Along the Meishan Hill, Zhejiang, China, a Permian–Triassic succession extends 2 km laterally and has been the subject of many paleontological and geochemical studies (Chen et al., Reference Chen, Yang, Luo, Benton, Kaiho, Zhao, Huang, Zhang, Fang, Jiang, Qiu, Li, Tu, Shi, Zhang, Feng and Chen2015), which combined make the Meishan composite section the only place that can currently be quantitatively investigated to better understand which environmental proxies relate to biodiversity loss. In addition, Meishan’s Permian–Triassic succession has a well-defined stratigraphic framework with each bed and sub-bed numbered allowing accurate correlations between studies performed over the last three decades. Furthermore, the Meishan D section has now been digitized into an open access, 3D interactive model that can be accessed via https://outcrop3d.deep-time.org/?model=194a20d7-958e-d871-11c2-ab18dfb62e16. In contrast, other regions with a rich paleontological and geochemical record for the Permian–Triassic transition, such as the Dolomites in Italy, do not yet have the same clear stratigraphic scheme or diversity in analyses that make studies like this one possible. During the Permian–Triassic transition, the Meishan section represents an outer slope setting in an equatorial (ca. 20°N) epicontinental sea (Yin et al., Reference Yin, Zhang, Tong, Yang and Wu2001). This means that the Meishan section can provide an analog into the causes of extinction during an extreme climate crisis for equatorial, shallow marine ecosystems (i.e., for the worst-case Shared Socioeconomic Pathways (SSP5–8.5) scenario, which predicts a total temperature increase of 3.3 to 5.7 °C by 2100 (IPCC, 2021). Here, we have conducted a data analysis by (a) creating a database of the fossil record to define the timing of extinction among different marine taxa, (b) assembled a database of 18 geochemical proxies (Figure 1 and Supplementary Table S1) for different environmental changes from the Meishan section that have been hypothesized to have had a critical role in the marine extinctions, and (c) quantitatively investigated which environmental changes associated with the climate crisis best explain the marine extinctions.
An outcrop and depositional model for the Meishan section. (A) The Deep-time Digital Earth 3D Outcrop model of the Meishan D section. The GSSP for the Wuchiapingian/Changhsingian and Permian/Triassic boundaries are indicated by stars. The interactive model can be accessed via https://outcrop3d.deep-time.org/?model=194a20d7-958e-d871-11c2-ab18dfb62e16. (B) Schematic of the paleoenvironmental setting, indicating the inorganic geochemical proxies that were selected to investigate the role of different environmental changes at Meishan. Data come from δ7Li (Sun et al., Reference Sun, Joachimski, Wignall, Yan, Chen, Jiang, Wang and Lai2012), δ13Ccarb (Shen et al., Reference Shen, Cao, Zhang, Bowring, Henderson, Payne, Davydov, Chen, Yuan, Zhang, Wang and Zheng2013), δ13Corg (Cao et al., Reference Cao, Love, Hays, Wang, Shen and Summons2009) supplemented at beds 23–38 with (Huang et al., Reference Huang, Luo, Bai and Tang2007) and (Sial et al., Reference Sial, Chen, Korte, Pandit, Spangenberg, Silva-Tamayo, Drude de Lacerda, Pinto Ferreira, Barbosa, Gaucher, Pereira and Riedel2021), δ15N (Cao et al., Reference Cao, Love, Hays, Wang, Shen and Summons2009), δ18Oapatite (J. Chen et al., Reference Chen, Shen, Li, Xu, Joachimski, Bowring, Erwin, Yuan, Chen, Zhang, Wang, Cao, Zheng and Mu2016), ∆33S and δ 34S (Shen et al., Reference Shen, Crowley, Wang, Bowring, Erwin, Sadler, Cao, Rothman, Henderson, Ramezani and Zhang2011a), δ44/40Caapatite (Hinojosa et al., Reference Hinojosa, Brown, Chun, DePaolo, Paytan, Shen and Payne2012), δ66Zncarb (Liu et al., Reference Liu, Wu, Shen, Jiang, Zhang, Lv, Zhang and Li2017), 87Sr/86Srapatite (Song et al., Reference Song, Wignall, Tong, Song, Chen, Chu, Tian, Luo, Zong, Chen and Lai2015), Th/Uapatite, ΩCeapatite, (Song et al., Reference Song, Wignall, Tong, Bond, Song, Lai, Zhang, Wang and Chen2012), δ114/110Cd (Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018), 187Os/188Os (Liu et al., Reference Liu, Selby, Zhang and Shen2020), Hg/TOC (Sial et al., Reference Sial, Chen, Korte, Pandit, Spangenberg, Silva-Tamayo, Drude de Lacerda, Pinto Ferreira, Barbosa, Gaucher, Pereira and Riedel2021), FeHR/Fetot and Fepy/FeHR (Xiang et al., Reference Xiang, Zhang, Schoepfer, Cao, Zheng, Yuan, Cai and Shen2020). The bathymetry follows Zhang et al. (Reference Zhang, Tong, Yin and Wu1997) at the time of the Permian/Triassic boundary.
Materials and methods
Fossil data
Using the Geobiodiversity Database (http://www.geobiodiversity.com), Paleobiology Database (https://paleobiodb.org) and a literature search, we constructed a database of all known fossil occurrences from the Meishan section that spans from the Longtan Formation (Wuchiapingian) to the Nanlinghu Formation (Dienerian). The wide breadths of these time intervals were chosen to reduce the impact of edge effects. The clades included in the dataset were the Arthropoda, Brachiopoda, Bryozoa, Chlorophyta, Chordata, Cnidaria, Foraminifera, Mollusca, Radiolaria, Rhodophyta, and Problematica. The occurrences were manually vetted to ensure that typographic errors were corrected, so species did not appear with multiple spellings, and to ensure that individual species were not represented within multiple genera in the database due to taxonomic synonymy, in which case the most up-to-date species identification was followed. Freshwater and palynomorph fossils were removed.
To calculate the stratigraphic range of each species, occurrences of species with open nomenclature (“”,?, aff., cf., informal) were taken into account. In older references from the Meishan section, the beds do not have the same stratigraphic subdivision as today and the occurrence of a specimen is considered present in all respective sub-beds. These occurrences that were not specified to a single bed/sub-bed as recognized in this study were subsequently excluded, as for them the timing of extinction is poorly constrained.
The resulting database included 603 species from 6,457 occurrences.
Inorganic geochemical data
To investigate changes in environmental conditions, we downloaded the raw datasets of inorganic geochemical proxies for the Meishan section. We initially obtained all the articles for each proxy investigated for the Meishan section. Where multiple records of a single proxy were collected, we selected the most robust record, that is, we avoided mixing datasets collected from the same beds by different studies and selected the most extensive record. The resulting dataset included 18 proxies (Supplementary Table S1; Supplementary Figures S1–S4). The sample heights were standardized according to the Permian/Triassic boundary, with 0 cm marking the base of bed 27c. Data from the Meishan core, which is located 550 m west of the Meishan section D (GSSP section), were scaled to correlate with the section from Meishan D as the beds demonstrate considerable thickness variations.
Statistical analysis
Determining the patterns of extinction can be confounded by subjective interpretations and by the Signor-Lipps effect. Therefore, to quantify the nature of extinction, here we used a modified version of the two-step extinction pulse algorithm of Wang and Zhong (Reference Wang and Zhong2018) (see also extended materials and methods).
To quantitatively determine the number of breakpoints in the segmented regression analysis, the selgmented() function from the segmented package was used (Muggeo et al., Reference Muggeo, Atkins, Gallop and Dimidjian2014). The segmented() function was then used to statistically determine where these breakpoints occur for each geochemical proxy. For data imputation of the geochemical data, a segmented regression was used because (a) it is less affected by anomalous data points, (b) estimates are based upon overall trends in the data, (c) it does not assume that a relationship between different proxies exists, and (d) it recognizes significant shifts in data trends and is more dynamic than a single regression model (see also Supplementary Figures S5–S6).
We applied General Linear Models (GLM) with a Poisson distribution to test the effects of multiple geochemical proxies on changes of species richness through the study interval. Only proxies that showed significant correlations with species richness were included in the model, and model selection was carried out by exploring the value-inflated factors and factors that are highly correlated were successively dropped from the model. This resulted in two GLMs, one where δ114/110Cd was dropped because it highly correlated with δ18Oapatite, and vice versa.
To investigate relationships between fossil incidence data and geochemical proxies, we carried out a partial-distance-based redundancy analysis using the Jaccard distance measure. Model selection was carried out using value-inflated factors, and factors that are highly correlated were successively dropped from the model. Variables that were insignificant in explaining incidence data dynamics using a permutation test for partial-dbrda were also dropped from the final model.
All analyses were carried in R v.3.4.3. Data and relevant code for this research work are stored and publicly available in GitHub: https://github.com/wjf433/QMEI
Results
Nature of the mass extinction event
The nature of the Permian–Triassic mass extinction is hotly debated and has been interpreted either as a single pulse, interval, or a two-pulsed extinction (Jin et al., Reference Jin, Wang, Wang, Shang, Cao and Erwin2000; Shen et al., Reference Shen, Crowley, Wang, Bowring, Erwin, Sadler, Cao, Rothman, Henderson, Ramezani and Zhang2011a; Song et al., Reference Song, Wignall, Chu, Tong, Sun, Song, He and Tian2014; Wang et al., Reference Wang, Sadler, Shen, Erwin, Zhang, Wang, Wang, Crowley and Henderson2014). Quantifying the number of pulses of extinction, which considers confidence intervals of stratigraphic ranges (Figure 2 and Supplementary Figures S7–S11), demonstrates that the nature of the mass extinction is complex and varies between different phyla. Near the Permian/Triassic boundary, the traditional extinction horizon (bed 25) (Wang et al., Reference Wang, Sadler, Shen, Erwin, Zhang, Wang, Wang, Crowley and Henderson2014) marks numerous last occurrences (LAD), leading to composition shifts (Figure 2 and Supplementary Figures S7–S11). Our analysis shows, however, that a single pulse of extinction (the final LADs) occurs at bed 28 for mollusks and bed 29a for foraminifera, brachiopods, and conodonts (Figure 2). Ostracods, instead, record two earlier pulses of extinction at beds 23a and 24d (Figure 2), the latter coinciding with a sequence boundary. The species richness of the remaining groups is too low to detect the timing of extinction, but the highest occurrences do not occur above bed 27c. This suggests that the nature of the mass extinction event at Meishan, except for ostracods, is best characterized as an extinction interval (51 cm), from beds 25 to 29a (C. meishanensis – I. isarcica conodont zones). Such an interpretation is also relatively consistent with a stark reduction in bioturbation and tiering depth at the base of bed 25 (Zhao and Tong, Reference Zhao and Tong2010).
Stratigraphic ranges of fossil species (vertical lines) from the Meishan section. Stratigraphic ranges of (A) foraminifera, (B) arthropods (all ostracods, except one trilobite species), (C) brachiopods, (D) mollusks, (E) conodonts, and (F) other (includes: bryozoans, corals, calcareous algae, and Tubiphytes). Quantitatively determined extinction pulses for each phylum indicated (horizontal red line). Singletons are excluded from the figure and from determining the number of extinction pulses. Bed numbers and sedimentology follow Zhang et al. (Reference Zhang, Tong, Yin and Wu1997) and Yin et al. (Reference Yin, Yu, Ding, Zhang, Tong and Yang1995). 0 meters is taken as the base of bed 27c, which is the biostratigraphic position of the Permian/Triassic boundary that is defined by the first appearance of Hindeodus parvus (Yin et al., Reference Yin, Zhang, Tong, Yang and Wu2001). W. = Wuchiapingian, Ind. = Induan, Yink. = Yinkeng Formation.
On the contrary, ostracods record an earlier major extinction interval from beds 22-23a and a subsequent minor pulse at bed 24d (Crasquin et al., Reference Crasquin, Forel, Feng, Yuan, Baudin and Collin2010), suggesting that this group of organisms was more sensitive to initial environmental changes or responded to different environmental changes before the main extinction interval. When the timing of extinction is investigated with all the species included, the mass extinction event is consistent between the different data splits, with a 2-pulse event at beds 23a and 29a, with bed 23a reflecting the selective extinction of ostracods (Supplementary Table S2). In addition to the extinction interval that spans the Permian/Triassic boundary, ostracods, mollusks, and brachiopods also record a minor extinction pulse earlier in the Changhsingian (beds 9 and 12, Figure 2).
These changes are also reflected by the breakpoints and a rapid decline in species richness at beds 23, 25, and 29a (Figure 3). Lithology changes and sequence stratigraphic boundaries play a role in determining the stratigraphic position of the LAD (Holland and Patzkowsky, Reference Holland and Patzkowsky2015; Nawrot et al., Reference Nawrot, Scarponi, Azzarone, Dexter, Kusnerik, Wittmer, Amorosi and Kowalewski2018; Zimmt et al., Reference Zimmt, Holland, Finnegan and Marshall2021). At Meishan, the extinction interval includes a lithostratigraphic boundary at the base of bed 25, and a transgressive surface at bed 27a, suggesting that these sedimentological changes affect our interpretations of the nature and timing of the extinction. Despite that, radiometric dating proposes that beds 25 to 28 only represent 60 ± 48 ka (Burgess et al., Reference Burgess, Bowring and Shen2014), and any hiatuses associated with sequence stratigraphic surfaces during the extinction interval are of relatively short duration.
Stratigraphic correlation of selected paleoenvironmental proxies with species diversity at the Meishan section, South China, with segmented regression lines overlain. δ13Ccarb (Shen et al., Reference Shen, Cao, Zhang, Bowring, Henderson, Payne, Davydov, Chen, Yuan, Zhang, Wang and Zheng2013), δ18Oapatite (VSMOW) (Chen et al., Reference Chen, Shen, Li, Xu, Joachimski, Bowring, Erwin, Yuan, Chen, Zhang, Wang, Cao, Zheng and Mu2016), δ114/110Cd (Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018), FeHR/Fetot (Xiang et al., Reference Xiang, Zhang, Schoepfer, Cao, Zheng, Yuan, Cai and Shen2020), ΩCeapatite (Song et al., Reference Song, Wignall, Tong, Bond, Song, Lai, Zhang, Wang and Chen2012), δ44/40Caapatite (Hinojosa et al., Reference Hinojosa, Brown, Chun, DePaolo, Paytan, Shen and Payne2012). The main extinction interval (beds 25–29a) is highlighted in orange and with two horizontal dashed lines. Note: only paleoenvironmental proxies that showed significant relationships with diversity are included, for a full figure with all the paleoenvironmental proxies see Supplementary Figures S1–S4.
Changes in species richness, in particular the earlier onset of ostracod extinctions, are problematic when trying to compare extinctions with geochemical proxies. This is because many of the proxies that have been investigated at the Meishan section only span a short interval, for example, δ114/110Cd only spans beds 22–33 (Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018), after species diversity has already started to decline (Figure 3). Analyses linking geochemical and fossil data were, therefore, restricted to beds 22–29a. Investigations did not extend beyond bed 29a because the protracted low diversity after the extinction interval can be attributed to a delayed recovery rather than environmental conditions.
Quantifying the causes of extinction
Quantifying the causes of extinction is complex, as environmental changes will manifest with different patterns and may be reflected as either state-shifts, anomalies, or correlations that can be associated with biodiversity dynamics. For example, the investigated proxies associated with volcanism, for example, Hg/TOC, δ66Zn and 187Os/188Os, are expected to appear as anomalies or spikes. The proxies Hg/TOC, δ66Zn and 187Os/188Os show anomalies that coincide with the onset of the mass extinction interval (Supplementary Figure S3), with the Hg/TOC, δ66Zn, and 187Os/188Os anomalies from beds 24b-24e being interpreted to reflect volcanism associated with the Siberian Traps coming along with input of volcanic ashes (Liu et al., Reference Liu, Wu, Shen, Jiang, Zhang, Lv, Zhang and Li2017; Liu et al., Reference Liu, Selby, Zhang and Shen2020; Sial et al., Reference Sial, Chen, Korte, Pandit, Spangenberg, Silva-Tamayo, Drude de Lacerda, Pinto Ferreira, Barbosa, Gaucher, Pereira and Riedel2021).
A segmented regression analysis, which can be used to quantify significant temporal shifts in proxies (i.e., state-shifts), recognizes significant changes for δ13Ccarb, δ18Oapatite, δ114/110Cd, and δ15N at the onset of the extinction interval (bed 25, Figure 3 and Supplementary Figures S1–S4). In addition, TOC shows a breakpoint at bed 22 (Supplementary Figure S1), corresponding with the main extinction pulse of ostracods (Figure 3). However, Th/Uapatite ratios show a state shift at bed 29 (Supplementary Figure S4), with an interpreted state shift from oxic to anoxic conditions (Song et al., Reference Song, Wignall, Tong, Bond, Song, Lai, Zhang, Wang and Chen2012), and corresponding with a plateau of low richness.
One issue with comparing the different proxies and changes in species richness or incidence data is the difference in resolution between the different datasets. To allow for statistical exploration of the data, the data were aggregated to the same resolution as the species dataset, that is, bed and sub-bed level resolution. In addition, not all the beds record proxy values, and, therefore, data were interpolated using the segmented regression curves for each proxy (Supplementary Figures S5 and S6). Another issue is that many of the different geochemical proxies correlate with one another (Supplementary Figures S12–S14), and these correlations are not necessarily direct causal effects, but could rather be associated with autocorrelation effects within each time series, indirect links, or common drivers (Runge et al., Reference Runge, Nowack, Kretschmer, Flaxman and Sejdinovic2019). This makes it difficult to disentangle whether the proxy is robust enough to interpret environmental changes, if an environmental change is causing a decline in diversity, or both diversity and proxy dynamics have a common cause. A correlation plot shows that δ13Ccarb, δ18Oapatite, δ114/110Cd, and δ15N are significantly correlated, which suggests these proxies share a common cause.
Many of the proxies show a correlation with changes in species richness (Supplementary Table S3). A Poisson regression model was performed to identify which proxies best explain the diversity dynamics. Value-inflated factors show that the correlation between δ13Ccarb, δ18Oapatite, δ114/110Cd, and δ15N significantly affects the quality of the model. For this reason and because the δ15N data is at a low resolution, δ15N was dropped from the model, whereas for δ114/110Cd and δ18Oapatite two separate models were run. The generalized linear models show that δ13Ccarb, δ18Oapatite, δ114/110Cd, and δ44/40Caapatite have significant relationships with changes in species diversity at Meishan (Table 1). In addition, no proxies showed a significant relationship between proxy variance and extinction rate.
Generalized linear model of significant environmental variables (geochemical proxies) and changes in diversity
Note: Model selection was based on proxies that showed consistent and significant linear relationships with diversity (Supplementary Table S3). δ 114/110Cd and δ 15N were dropped from the first model because they showed a significant correlation with δ 18O apatite that negatively impacted the model (Supplementary Material). A second model swapping δ 18O apatite and δ 114/110Cd was done to investigate the best model.
A partial-distance-based redundancy analysis (partial-dRDA) was undertaken to investigate the changes in fossil incidence data for beds 22 to 29a and changes in geochemical proxies (Figure 4). Once more, value-inflated factors show that the correlation between δ13Ccarb, δ18Oapatite, δ114/110Cd, and δ15N significantly affects the quality of the model. δ114/110Cd and δ15N were, therefore, dropped from the model. The partial-dRDA showed that δ13Ccarb, δ18Oapatite, and 87Sr/86Sr best-explained changes in the incidence data. Swapping δ114/110Cd with δ18Oapatite shows that δ18Oapatite is a more significant proxy for explaining incidence data dynamics. When only the significant factors are included in the partial-dRDA model, only δ13Ccarb and δ18Oapatite record significant relationships (Figure 4). It is also evident that the fossil incidence data cluster according to lithology (Figure 4), highlighting how lithological changes reflect changes in the environment affecting species loss.
Partial-Distance-based Redundancy Analysis (capscale) for fossil assemblages and geochemical proxies from the Meishan section. Included vectors are the geochemical proxies that were determined as having a significant relationship with the fossil assemblages. Sample point shapes relate to bed lithology: filled circles = limestone, open squares = silty limestone, and filled triangles = clay. The bed numbers for each assemblage are indicated, and only beds 22-29a are included due to limited coverage of geochemical proxies at the Meishan section. Smooth contours of the oxygen isotope values underlie the ordination plot to demonstrate the relationship with the fossil assemblages.
Discussion
Due to the large suite of geochemical proxies investigated for the Meishan section, a number of different environmental changes have been proposed as possible causes of the Permian–Triassic mass extinction. Our quantitative analysis, combining diversity and proxy data, demonstrates that δ13Ccarb, δ18Oapatite, and δ114/110Cd are key for understanding the cause(s) of the mass extinction event. The geochemical signature of these proxies is mostly generated in the euphotic zone with a high chance of transfer into sediments without much alteration, suggesting that the relationships between these proxies and the fossil record reflect past environmental-life interactions. In addition, the lack of relationship between the fossil incidence data and the geochemical proxies that are affected by lithology changes supports that our interpretations are robust. δ13Ccarb is often interpreted as reflecting the release of large quantities of isotopically light carbon into the atmosphere, changes in primary productivity, and changes in carbon burial rates (Cao et al., Reference Cao, Wang and Jin2002). δ13Ccarb can, therefore, signify an environmental disturbance and even the trigger of the mass extinction (Cui et al., Reference Cui, Kump, Ridgwell, Schmidt, Fristand and Elkins-Tanton2015), but it cannot be inferred to identify the underlying environmental changes that drove species to extinction. In addition, δ13Ccarb, δ18Oapatite, and δ114/110Cd are significantly correlated, which we infer as being impacted by a common cause.
A negative excursion in δ18Oapatite is interpreted to reflect a rapid, 8–12 °C, warming associated with the mass extinction event at Meishan (Joachimski et al., Reference Joachimski, Lai, Shen, Jiang, Luo, Chen, Chen and Sun2012; Sun et al., Reference Sun, Joachimski, Wignall, Yan, Chen, Jiang, Wang and Lai2012; Chen et al., Reference Chen, Shen, Li, Xu, Joachimski, Bowring, Erwin, Yuan, Chen, Zhang, Wang, Cao, Zheng and Mu2016) and is consistently the best explanatory factor for diversity dynamics at Meishan (Table 1 and Figure 4). Thermal stress is understood to limit the performance of aerobic marine organisms, because the pejus temperature is close to the temperature optimum on the upper thermal limit, and increasing temperatures beyond a marine organism’s optimum temperature range rapidly lead to a reduction in the aerobic scope of marine organisms (excess capacity supporting activity, growth, and reproduction) (Pörtner, Reference Pörtner2012; Pörtner et al., Reference Pörtner, Bcok and Mark2017). An expectation from this mechanism would be an observed decrease in body size as temperatures increase and primary productivity declines. The only body size data from Meishan with enough measurements comes from two species of foraminifera, Diplosphaerina inaequalis and Frondina permica (Song et al., Reference Song, Tong and Chen2011), that record decreasing size with more negative δ18Oapatite values (Figure 5), suggesting a decrease in aerobic scope before their LADs in beds 27c and 29a, respectively. The paleoequatorial setting of the Meishan section also means that this locality would also have experienced some of the highest climate velocities (sensu Burrows et al., Reference Burrows, Schoeman, Buckley, Moore, Poloczanska, Brander, Brown, Bruno, Duarte, Hakpern, Holding, Kappel, Kiessling, O’Connor, Pandolfi, Parmesan, Schwing, Sydeman and Richardson2011) with Earth system models for the Permian–Triassic mass extinction demonstrating the consequent loss of aerobic habitats able to support the metabolism of marine ectotherms (Penn et al., Reference Penn, Deutsch, Payne and Sperling2019). This idea is supported, by the poleward migrations of Permian holdover radiolarians, sponges, and conodonts away from equatorial settings associated with the warming (Foster et al., Reference Foster, Asatryan, Rauzi, Botting, Buchwald, Lazarus and Kiessling2023b). The significant relationship between δ18Oapatite with diversity, compositional changes, body size and poleward migrations supports the interpretation that temperature-driven hypoxia was fundamental in causing equatorial extinctions/extirpations during the Permian–Triassic climate crisis.
A scatter plot showing the relationship between foraminifera test size and δ18Oapatite (a temperature proxy) from the Meishan section. (A) Measurements of Frondina permica from beds 13a-27c (B) Measurements of Diplosphaerina inaequalis from beds 13a-29a. An order 2 polynomial trend line is underlain to illustrate the relationship between δ18Oapatite values and test size. As δ18Oapatite values decrease they can be used to infer a warming of the climate and vice versa, which is indicated on the axis. The test sizes of the foraminifera were converted to geometric sizes and log-transformed. Body size data are from Song et al. (Reference Song, Tong and Chen2011) and oxygen isotope data from Chen et al. (Reference Chen, Shen, Li, Xu, Joachimski, Bowring, Erwin, Yuan, Chen, Zhang, Wang, Cao, Zheng and Mu2016). Other species could not be included due to either the lack of species-level identification or size measurements.
δ18Oapatite is also significantly correlated with δ114/110Cd and δ15N, where both δ114/110Cd and δ15N record negative excursions that are interpreted to reflect primary productivity dynamics (Cao et al., Reference Cao, Love, Hays, Wang, Shen and Summons2009; Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018). δ114/110Cd reflects nutrient utilization by phytoplankton and is an indirect proxy for primary productivity. At Meishan a negative excursion in δ114/110Cd coincides with the mass extinction interval reflecting a collapse in primary productivity (Figure 2; Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018), which is also associated with the major extinction of radiolarians (the primary fossil record of planktic biodiversity) (O’Dogherty et al., Reference O’Dogherty, Carter, Goričan and Dumitrica2010). A stark reduction in primary productivity would be catastrophic for marine ecosystems because the cascading effect of extinction causes knock-on effects on species populations in successive layers of the marine food web (Huang et al., Reference Huang, Chen, Roopnarine, Benton, Zhao, Feng and Li2023). Huang et al. (Reference Huang, Chen, Roopnarine, Benton, Zhao, Feng and Li2023) inferred that the loss of radiolarians at the onset of the extinction interval could have been a cascading effect of a collapse in primary productivity. The migration of radiolarians to thermal refugia at higher latitudes and deeper waters (Foster et al., Reference Foster, Asatryan, Rauzi, Botting, Buchwald, Lazarus and Kiessling2023b) could also suggest additive effects of different environmental stressors, that is, nutrient and thermal stress. In addition, aerobic metabolism is not only affected by thermal stress but also nutritional stress, which can exacerbate the effects of climate change in marine ectotherms (Saulsbury et al., Reference Saulsbury, Moss, Ivany, Kowalewski, Lindberg, Gillooly, Heim, McClain, Payne, Roopnarine, Schöne, Goodwin and Finnegan2019). The impacts on the aerobic metabolism of marine organisms can also be inferred from the observed decrease in the body size of surviving taxa (He et al., Reference He, Twitchett, Zhang, Shi, Feng, Yu, Xu and Peng2010).
The hypothesis that primary productivity collapsed during the Permian–Triassic mass extinction is controversial, where some proxies suggest a collapse (e.g., Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018) while others suggest enhanced primary productivity (e.g., Qiu et al., Reference Qiu, Song, Hu, Wignall and Song2019). In part, this can be explained by spatial heterogeneity in primary productivity rates as evidenced by spatial variations in productivity proxies (e.g., Shen et al., Reference Shen, Schoepfer, Feng, Zhou, Yu, Song, Wei and Algeo2015), but conflicting results are also known from the Meishan section (e.g., Zhang et al., Reference Zhang, Wen, Zhu, Fan and Cloquet2018; Qiu et al., Reference Qiu, Song, Hu, Wignall and Song2019). δ114/110Cd is not a redox proxy, but in sulfide-bearing anoxic sediments there are more negative values than in oxic surface waters (Hohl et al., Reference Hohl, Galer, Gamper and Becker2017). Zhang et al. (Reference Zhang, Wen, Zhu, Fan and Cloquet2018), however, noted that there is no correlation between sulfur concentrations and δ114/10Cd values, inferring that the inferred trends reflect changes in primary productivity. Other proxies for primary productivity come from lipid biomarkers and δ15N, which show changes in the archaeal and bacterial communities as a consequence of the environmental changes (Xie et al., Reference Xie, Pancost, Yin, Wang and Evershed2005; Cao et al., Reference Cao, Love, Hays, Wang, Shen and Summons2009). Taken together, this suggests that not only did primary productivity appear to collapse in South China, but there were also phytoplankton community shifts which would have led to nutrient stress.
In the equatorial paleosetting of the Meishan section, both thermal and nutrient stress are interpreted to best explain the extinctions. It has been shown that despite inter-specific differences, there are clear differences in hypoxia tolerance among higher taxa (Song et al., Reference Song, Wignall, Tong and Yin2013). Ostracods and crustaceans have the least tolerance to hypoxia compared to other invertebrate groups (Song et al., Reference Song, Wignall, Tong and Yin2013), and their earlier onset of extinction at the Meishan section also corresponds to the initial changes in the δ18Oapatite negative excursion (Figure 2). δ18Oapatite dynamics have been divided into two phases associated with different rates and magnitudes of warming (Wu et al., Reference Wu, Cui, Chu, Song, Tong, Del Corso and Ridgwell2023), with the first phase, coinciding with the major extinction of ostracods, being slower and of a smaller magnitude. Pre-extinction changes in δ18Oapatite in equatorial settings have also been related to body size changes in ammonoids (Gliwa et al., Reference Gliwa, Wiedenbeck, Schobben, Ullmann, Kiessling, Ghaderi, Struck and Korn2022) and correspond to pre-extinction changes in brachiopod assemblages (Zhang et al., Reference Zhang, Shi, Wu, Yang, He, Yuan and Lei2017). This suggests that pre-extinction slower warming and the following rapid warming led to different timings of extinction for different marine organisms, depending on their sensitivity to temperature and oxygen-concentration changes.
Widespread anoxic conditions throughout shallow and deep marine basins have long been associated with marine extinctions. The Meishan section has been the subject of several paleoredox studies, utilizing: pyrite framboids (Wei et al., Reference Wei, Zhang and Qiu2020; Chen et al., Reference Chen, Yang, Luo, Benton, Kaiho, Zhao, Huang, Zhang, Fang, Jiang, Qiu, Li, Tu, Shi, Zhang, Feng and Chen2015), sulfur isotopes (Shen et al., Reference Shen, Farquhar, Zhang, Masterson, Zhang and Wing2011b) as well as iron speciation, and redox-sensitive metals (Xiang et al., Reference Xiang, Zhang, Schoepfer, Cao, Zheng, Yuan, Cai and Shen2020). Iron speciation has been established as a proxy for local water column redox in clastic successions and some carbonate successions, yet the FeHR/Fetot data for Meishan records almost persistently anoxic conditions for bed 21 to 34 (Xiang et al., Reference Xiang, Zhang, Schoepfer, Cao, Zheng, Yuan, Cai and Shen2020). Considering the abundant and deeply penetrated trace fossils (>20 cm depth) found throughout most of this interval (Zhao and Tong, Reference Zhao and Tong2010), the interpretation of consistently anoxic conditions is equivocal. Hence, the anoxic FeHR/Fetot signals at best represent dynamic redox conditions with intermittent anoxic intervals or an alternative explanation is the redeposition of dissolved Fe at Meishan that was released from nearby oxygen minimum zones (similar to the Guaymas Basin in the Gulf of California; Scholz et al., Reference Scholz, Schmidt, Hensen, Eroglu, Geilert, Gutjahr and Liebetrau2019). Multiple sulfur isotope signals from pyrite and pyrite aggregate sizes are suggested to support the development of episodic anoxic water column conditions during the deposition of beds 22–24 (Shen et al., Reference Shen, Farquhar, Zhang, Masterson, Zhang and Wing2011b). Depending on the extent and duration of these anoxic episodes, they could have contributed to the loss of ostracods before the main extinction interval, as deposit-feeding ostracods are negatively impacted by falling oxygen levels (Lethiers and Whatley, Reference Lethiers and Whatley1994). Conversely, ΩCeapatite anomalies from the same beds mean deposition was in an oxygenated setting (Song et al., Reference Song, Wignall, Tong, Bond, Song, Lai, Zhang, Wang and Chen2012). Above bed 24d, redox-sensitive metal enrichment factors record decreasing trends (Xiang et al., Reference Xiang, Zhang, Schoepfer, Cao, Zheng, Yuan, Cai and Shen2020), reflecting increasing oxygenation of the water column. This is also supported by pyrite δ34S and ∆33S signals that do not support an anoxic water column interpretation (Shen et al., Reference Shen, Farquhar, Zhang, Masterson, Zhang and Wing2011b). Overall, the water column redox conditions at Meishan were likely dynamic during the Permian–Triassic transition, but the timing and extent of anoxic episodes are too poorly constrained to unequivocally conclude that anoxia or dysoxia played a role in the relatively shallow setting of Meishan.
Despite the intense geochemical and paleontological research on the Meishan section, this study highlights some limitations that must be addressed in future research. The restriction of investigations of geochemical proxies over a short interval at the Permian/Triassic boundary hinders our ability to understand how environmental conditions evolved over the Changhsingian and how that relates to the climate crisis (e.g., δ44/40Caapatite; Hinojosa et al., Reference Hinojosa, Brown, Chun, DePaolo, Paytan, Shen and Payne2012). Therefore, even though δ44/40Caapatite, a potential proxy for ocean pH (Hinojosa et al., Reference Hinojosa, Brown, Chun, DePaolo, Paytan, Shen and Payne2012), is recorded as having a significant relationship with changes in species richness (Table 1), the short record still makes this interpretation equivocal. Even if the δ44/40Caapatite trends are seen as robust, concerns of using this proxy to determine ocean acidification have been raised (Komar and Zeebe, Reference Komar and Zeebe2016; Foster et al., Reference Foster, Hirtz, Farrell, Reistroffer, Twitchett and Martindale2022b), and, therefore, the role of ocean acidification in the extinctions for these settings is still unknown. In addition, the lack of abundance data and other ecological data from paleontological studies (e.g., Song et al., Reference Song, Tong and Chen2009; Crasquin et al., Reference Crasquin, Forel, Feng, Yuan, Baudin and Collin2010) means it is not yet possible to investigate the ecological impacts of the Permian–Triassic climate crisis beyond the timing of extinction. Therefore, several ecological changes, such as changes in relative abundance, dominance, or body size and how they relate to environmental changes, cannot yet be explored. Finally, the cause(s) of extinction is expected to vary on various spatial scales, and more high-resolution studies from other sections and regions are, therefore, required. Despite these shortcomings, our statistical analysis demonstrates that the extreme impact of environmental changes on the aerobic metabolism of marine ectotherms and the cascading effects of extinction best explain the cause of extinction in epicontinental, equatorial settings. This means that for the worst-case RCP scenario, the biggest climate threats to modern-day shallow marine, equatorial biodiversity are thermal and nutrient stress.
Open peer review
To view the open peer review materials for this article, please visit http://doi.org/10.1017/ext.2024.9.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/ext.2024.9.
Data availability statement
All analyses were carried in R v.3.4.3. Data and relevant code for this research work are stored and publicly available in GitHub: https://github.com/wjf433/QMEI.
Acknowledgements
We would like to thank Steve Wang (Swathmore College) for his work and help in developing the parallelization of his extinction pulse algorithm. We would like to thank the CEN-IT at Universität Hamburg for access to the computer cluster required for our analyses. We would like to thank all the scientists who have entered fossil occurrence data into both the Geobiodiversity Database and Paleobiology Database. This is Paleobiology Database official publication No. 487 and Geobiodiversity Database official publication No. 2. W.J.F. and X.W would like to thank Prof. S-Z. Shen and Dr. H. Xu (Nanjing University) for taking and introducing them to Meishan. We would like to thank Prof. H. Song for providing the foraminifera body size data. We would also like to acknowledge the entire Deep-time Digital Earth Outcrop3D team for digitizing the Meishan section.
Author contribution
W.J.F. conceived the project. W.J.F., A.B.F., and Q.L. collected the data. W.J.F. wrote the code and ran the analysis. W.J.F. wrote the initial manuscript draft and all authors contributed thereafter.
Financial support
This research was funded by the DFG Grant FO 1297/1 awarded to W.J.F.
Competing interest
The authors declare no competing interests exist.
Open access
Comments
Dear Editor,
We hereby submit our paper entitled “Thermal and nutrient stress drove Permian-Triassic shallow marine extinctions” for consideration as a research article.
The research presented in this manuscript represents a comprehensive meta-analysis that investigates the causes of marine extinctions during the Permo-Triassic climate crisis. Our study combines the examination of geochemical proxies, paleontological and sedimentological changes, and statistical analyses to provide valuable insights into this significant event. One novel aspect of this study is that our results propose that a combination of thermal and nutrient stress best account for the extinctions at equatorial latitudes. Essentially these results provide an analog for the worst-case RCP scenario, which predict a total temperature increase of more than 4.5°C by 2100.
We anticipate that the results of this work will stimulate future research and will become a highly cited piece of research, especially as we provide a template for future studies looking to integrate high-resolution data and identify a number of different avenues that future research will need to undertake to better understand how climate crises relate to mass extinction events.
This discovery is not only interesting to geochemists, earth-system modelers, sedimentologists and paleontologists as we propose a novel combination for the cause of equatorial marine extinctions (e.g., Penn et al. 2013 Science 362, eaat1327), but a general audience interested in the impact of climate crises and extinction, other scientists who are working on “hot-topic” climate-feedback mechanisms, and policy makers that use geological results to make key decisions (e.g., IPCC Working Group II).
We are currently not producing any other articles and are not aware of any other articles reporting on the same topic as the submitted manuscript. However, this project is the first from a series of projects from the €1.6 million DFG-funded Ocean life on the brink project (Project No. FO1297/1-1) and will serve as a template for future research and research directions that will stem from this work.
Corresponding author: William Foster, william.foster@uni-hamburg.de
If you have any further questions or require any more information, please do not hesitate to contact me. I look forward to hearing from you.