Respiratory protein-driven selectivity during the Permian-Triassic mass extinction

Extinction selectivity determines the direction of macroevolution, especially during mass extinction; however, its driving mechanisms remain poorly understood. By investigating the physiological selectivity of marine animals during the Permian-Triassic mass extinction, we found that marine clades with lower O2-carrying capacity hemerythrin proteins and those relying on O2 diffusion experienced significantly greater extinction intensity and body-size reduction than those with higher O2-carrying capacity hemoglobin or hemocyanin proteins. Our findings suggest that animals with high O2-carrying capacity obtained the necessary O2 even under hypoxia and compensated for the increased energy requirements caused by ocean acidification, which enabled their survival during the Permian-Triassic mass extinction. Thus, high O2-carrying capacity may have been crucial for the transition from the Paleozoic to the Modern Evolutionary Fauna.


Fossil occurrence data
Fossil data used to calculate diversity variation were obtained from a previously published database of Permian-Triassic marine fossils 1,2 .The database contains 52,322 occurrences at the generic level from 1,768 published papers and the Paleobiology Database, spanning the late Permian Changhsingian to the Late Triassic Rhaetian (Data S1).Our analysis is based on the occurrences of genera, as specieslevel identifications is often inaccurate and inconsistent.For example, a specimen may be species A in Europe and species B in Asia, even though they are the same species.Within the considered interval, a total of 1,097 genera belong to 13 major clades, including two clades of protozoa (foraminifera and radiolarians), nine clades of invertebrates (corals, sponges, brachiopods, bryozoans, ostracods, cephalopods, gastropods, bivalves, and echinoderms), and two clades of vertebrates (conodonts and fishes).For marine arthropods, we used only ostracod data because ostracods are abundant in the fossil record during the late Permian.Other marine arthropods are very rare in this time interval.For example, only two genera of trilobite, one genus of chelicera, and one genus of decapod are recorded in the Changhsingian bin compared to > 100 genera of ostracods in the Paleobiology Database.We did not consider background extinction in the late Permian because many studies have shown that the background extinction rate of marine taxa in the Changhsingian was negligible compared to the mass extinction interval around the Permian-Triassic boundary [3][4][5][6] .
Therefore, the results using the Changhsingian and Induan fossil data reflect a selectivity pattern of the Permian-Triassic mass extinction rather than background extinction.

Body size data
We updated the comprehensive database of Schaal et al. 7 to assign body size expressed as the maximum length for each species (see Data S2).Using the maximum size for each taxon is a common approach for body size studies, as the effects of juvenile specimens in the database can be avoided [7][8][9][10][11][12] .We followed the same methods to compile additional data for taxa not included in this database.A number of recently published databases were used to compile the size data, including references [13][14][15][16][17][18] .Other size data were mainly obtained from the published taxonomic literature.Only common taxa from both Changhsingian and Induan are included because these taxa have abundant fossil data to study their size change during the Permian-Triassic interval, i.e., foraminifera, brachiopods, ostracods, gastropods, cephalopods, bivalves, conodonts, and fishes.Other taxa, including corals, sponges, radiolarians, bryozoans, and echinoderms, are absent/very rare in the Induan bin and accordingly are not included in this study.The Changhsingian and Induan body size dataset is composed of 1495 species in 635 genera belonging to eight common clades (Data S2).

Extinction
The environmental crisis spanning the Permian-Triassic boundary, e.g., extreme high temperatures 19 , ocean anoxia [20][21][22] , and potential ocean acidification 23 , was just the start of a series of perturbations that extended throughout the Early Triassic 19,24,25 .Prolonged environmental stress not only exacerbated species extinctions but also negatively affected Early Triassic recovery, which further hampered the rebound of biodiversity 2,[26][27][28] .Here, we focus on the selectivity in extinction of animals with different respiratory proteins.
To investigate selectivity, we used multiple methods to calculate the extinction, i.e., the proportion 29 , per capita 29 , three timer 30 , and gap filler 31 extinction estimators.These estimators can decrease the Signor-Lipps effect and edge effects to some extent 31 .Error bars for the proportion of extinction represent binomial 95% confidence intervals 32 .Considering that the extinction episodes straddled the Permian-Triassic boundary [3][4][5]33 and lasted approximately 60 kyr from the latest Changhsingian to earliest Induan 34 , species that occurred in the extinction interval from the Clarkina meishanensis to Isarcicella staeschei zones are included in the Changhsingian bin (Figure S2).

Body size reduction
Body size change during the Permian-Triassic mass extinction was analyzed by the Mann-Whitney U test and resampling method.We performed a two-tailed paired Mann-Whitney U test using SPSS (version 25).We used a bootstrap method with 1000 replicates to estimate the size reduction of major clades and oxygencarrying types.Size reduction between the two adjacent time bins, i.e., Changhsingian and Induan, is represented by the following equation: Size reduction = log10 (SizeCh) -log10 (SizeIn), where SizeCh = median value of Changhsingian body size; SizeIn = median value of Induan body size.The standard deviation (1 SD) was calculated from 1000 bootstrap replicates of size reduction.Resampling experiments were performed in R. The results of the Mann-Whitney U test and the resampling method matched well (Figure 1).Size data (Data S4-6) were analyzed at both the genus level (Figure 1, Figure S3) and species level (Figure S4).The genus-level analysis of size change parallels the genus-level analysis of mass extinction.The results of the species-level analysis are consistent with the genus-level results, both supporting the significant selectivity among taxa with different oxygen carrying capacities.

Logistic regression
We applied multiple logistic regression models to evaluate the association of physiological and ecological traits with extinction.This analytical approach has been successfully applied to assess extinction selectivity in past and future biodiversity crises 12,[35][36][37][38] .We performed logistic regression analysis using SPSS (version 25).
Logistic regression was performed to evaluate the relationship between physiological and ecological variables and selectivity in extinction.The analysis for evaluating extinction used pre-extinction genera and their extinction (0) and survival (1) as the outcome variables (Data S3).
The logistic regression analysis used mineralogy, oxygen-carrying capacity, physiological buffering capacity, geographic range, motility, species richness and circulatory system as predictors.A multivariate framework for oxygen-carrying capacity, physiological buffering capacity, and motility was used, as it provides a more detailed approach 39,40 (Table S3).Oxygen carrying capacity is a continuous variable, expressed as the maximum O2-carrying capacity of each oxygen carrier (Table S1).Considering the variation in oxygen carrying capacity for each respiratory protein, we also used the median O2-carrying capacity to perform a logistic regression analysis.The results showed a consistent selectivity between extinction and O2carrying capacity (Tables S4-S8).There are also uncertainties in diffusion.For example, temperature can affect the rate of diffusion.The physiological buffering capacity category includes no carbonate load and less vulnerable to hypercapnia, moderate carbonate load and some buffering capacity, and heavy carbonate load with little buffering capacity following references [41][42][43][44][45] .The mineralogy category consists of carbonate and non-carbonate following references 36,46 .Silica and phosphate skeletons are not ranked because they are unordered categorical variables and were combined with soft-bodied animals into one variable, i.e., non-carbonate.Motility is composed of fast mobile, slow mobile, facultatively mobile, stationary unattached and stationary attached 40,43 .The circulatory system includes no circulatory system, open circulatory system, and closed circulatory system (see Table S2).The occurrence is the number of occurrences in the fossil database.Geographic range and occurrence are continuous variables and were log-transformed (loge).The categories of physiological buffering capacity, mineralogy, motility, and circulatory system are ordered factor variables.A correlation matrix of the variables shows that a few variables are correlated, e.g., geographic range and number of occurrences (Table S9).
In general, the correlation values are low.Variance inflating factors (VIFs) and tolerances were used to analyze the severity of multicollinearity in the regression analysis (Tables S6, S8).The results of the analysis including all seven variables show a maximum VIF of 10.6 in the physiological buffering capacity category, suggesting significant multicollinearity.This is likely due to the association between variables of the circulatory system and physiological buffering capacity.Therefore, we performed logistic regression after removing the circulatory system.We chose to keep physiological buffering capacity because it is a more common variable for the analysis of extinction selectivity.The results using six variables with variance VIF < 7 and tolerances > 0.1 (see Table S6) suggest that the estimated coefficient is not inflated by other factors.

Meta-analysis
We used a mixed-effects meta-regression model to perform linear regression between maximum oxygen capacity and proportion of extinction and size reduction 47 .
Data that were used to perform the meta-analysis are shown in Data S8-S11.
Different from standard linear regression, this method could consider the sampling variances or uncertainty of data, as our extinction magnitude and size reduction data have uncertainty.This analysis was performed by the function rma in meta for the R package 47 .Regression lines and their 95% confidence intervals, regression coefficients, and p values in Figure 1 and Figure S3 are derived from the analysis results of function rma in meta.

Earth system model simulations
The cGENIE, an Earth system model of intermediate complexity, is used to simulate ocean O2, H2S, and pH changes during the Permian-Triassic mass extinction.
We ran the model on a 36x36 grid with 16 vertical levels in the ocean, with specific modules and boundary conditions taken from previous work 48 .In particular, a temperature-dependent POM remineralization module is applied in this simulation, which describes the temperature-dependent microbial metabolism that has a significant effect on ocean anoxia and euxinia 48,49 .The ocean phosphate concentrate is set as 1.0 and 2.0 times the present oceanic levels (2.159 μmol kg −1 ) in the late Permian and earliest Triassic, respectively, following reference 48 .We also test other levels of ocean phosphate concentrations in the earliest Triassic (1.0, 1.5, 2.5 and 3.0 × modern levels, see Figs.S5-S7).We only modified the boundary conditions of atmospheric CO2 levels at the Permian-Triassic boundary, which are different from those set in reference 48 .The pCO2 in our simulations comes from previous CO2 reconstruction based on high-resolution δ 13 C of C3 plant remains from terrestrial sections in southwestern China 50 , which is set as 426 ppmv in the late Permian and 2507 ppmv in the earliest Triassic.All experiments are run for 10 kyr to reach the steady state.Table S3.Ecological traits used to analyze the selectivity of extinction.Physiological buffering capacity after references [41][42][43][44][45] .Shell mineralogy after references 36,46 .Silica and phosphate skeletons are not ranked because they are unordered categorical variables, and were combined into one variable, i.e., non-carbonate.Motility after references 40,43 .The values of O2carrying capacity and circulatory system types are from references in Tables S1 and S2.

Occurrence
Natural logarithmic (number of occurrences)

Fig. S1 .
Fig. S1.The relationships between extinction rate and maximum oxygen capacity of marine animals with carbonate shells.(A) Per capita extinction.(B) Three timer extinction.(C) Gap filler extinction.

Fig. S3 .
Fig. S3.Changes in body size of marine animals at the genus level during the Permian-Triassic mass extinction.For each clade, four dotplots and boxes represent the size distribution of, from left to right, victims, survivors (pre, Changhsingian), survivors (post, Induan), and new comers (originators).Diff, Hr, Hc, and Hb represent different types of animals that use diffusion, hemoglobin, hemocyanin, and hemoglobin, respectively, to transport oxygen from their surroundings to their bodies.Bivalve-non p represents hemoglobin non-protobranchia bivalve.

Fig. S4 .
Fig. S4.Changes in the median size of marine animals at the species level during the Permian-Triassic mass extinction.(A) Body size reduction of major marine clades.(B) Correlation between size reduction and oxygen capacity for marine clades.(C) Correlation between size reduction and oxygen capacity for animals with carbonate shells.Vertical bars represent the standard deviation (1 SD) of 1000 bootstrap replicates.Diff (orange color), Hr (purple color), Hc (blue color), and Hb (magenta color) represent different types of animals that use diffusion, hemerythrin, hemocyanin, and hemoglobin, respectively, to transport oxygen from the surroundings to their bodies.Bivalve-non p represents hemoglobin non-protobranchia bivalve.Dashed lines and shades represent linear regression lines and their 95% confidence intervals.For the methods used to calculate size reduction, see Materials and Methods.For fossil occurrence and size data, see Data 1 and 2.

Table S1 . Types and characteristics of oxygen carriers.
In this work, we assumed that 250 million-year-old organisms had the same protein type and O2-carrying capacity as modern organisms of the same clade.

Table S2 . Types of oxygen carriers and circulatory systems in animals
. (A) Types of oxygen carriers.(B) Types of circulatory systems.

Table S4 . Results of logistic regression for extinction. (A)
Results of logistic regression including all six variables.(B)Results of logistic regression including four variables that have significant logistic regression relationships in panel A. (C)Results of logistic regression with O2carrying capacity as a categorical covariate.OCC, maximum value of oxygen carrying capacity; SE, standard error; Exp, exponential; CI, confidence interval.

Table S5 . Results of logistic regression for extinction using median values of oxygen carrying capacity
. OCC, median value of oxygen carrying capacity; SE, standard error; Exp, exponential; CI, confidence interval.

Table S7 . Results of logistic regression including all seven variables
. OCC, oxygen carrying capacity; SE, standard error; Exp, exponential; CI, confidence interval.

Table S8 . Results of multicollinearity diagnostics including all seven variables
. OCC, oxygen carrying capacity; SE, standard error; VIF, variance inflating factor.