# Decreasing Phanerozoic extinction intensity as a consequence of Earth surface oxygenation and metazoan ecophysiology

Richard Stockey, Alexandre Pohl, Andy Ridgwell, Seth Finnegan, Erik Sperling

► To cite this version:

Richard Stockey, Alexandre Pohl, Andy Ridgwell, Seth Finnegan, Erik Sperling. Decreasing Phanerozoic extinction intensity as a consequence of Earth surface oxygenation and metazoan ecophysiology. Proceedings of the National Academy of Sciences of the United States of America, 2021, 118 (41), pp.e2101900118. 10.1073/pnas.2101900118. hal-03364856

HAL Id: hal-03364856

https://hal.science/hal-03364856

Submitted on 6 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
## Main Manuscript for

# Decreasing Phanerozoic extinction intensity as a consequence of Earth surface oxygenation and metazoan ecophysiology

Richard G. Stockey<fnref n="1" />*, Alexandre Pohl<fnref n="2" />,<fnref n="3" />, Andy Ridgwell<fnref n="2" />, Seth Finnegan<fnref n="4" />, and Erik A. Sperling<fnref n="1" />

**Email:** Richard Stockey (rstockey@stanford.edu)

**Author Contributions:** R.G.S. and E.A.S. designed the research; R.G.S. conducted Earth system modeling with input from A.R.; A.P. and A.R. generated the Ordovician model configuration; R.G.S. conducted ecophysiological modeling; R.G.S. and S.F. conducted quantitative paleobiological analyses; R.G.S wrote the paper with contributions from all authors.

**Classification:** Physical Sciences/Earth, Atmospheric, and Planetary Sciences

**Keywords:** Extinction, oxygen, ecophysiology, temperature-dependent hypoxia, Earth system evolution

**This PDF file includes:**

- Main Text

- Figure Captions

<footnote n="1">¹Department of Geological Sciences, Stanford University, Stanford, CA, USA</footnote>

<footnote n="2">²Department of Earth and Planetary Sciences, University of California, Riverside, CA, USA</footnote>

<footnote n="3">³Biogéosciences, UMR 6282, UBFC/CNRS, Université Bourgogne Franche-Comté, 6 boulevard Gabriel, F-21000 Dijon, France</footnote>

<footnote n="4">⁴Department of Integrative Biology, University of California, Berkeley, Berkeley, CA, USA</footnote>

<footnote>*To whom correspondence should be addressed.</footnote>
# 1 Abstract

The decline in background extinction rates of marine animals through geologic time is an established but unexplained feature of the Phanerozoic fossil record. There is also growing consensus that the ocean and atmosphere did not become oxygenated to near-modern levels until the mid-Paleozoic, coinciding with the onset of generally lower extinction rates. Physiological theory provides us with a possible causal link between these two observations – predicting that the synergistic impacts of oxygen and temperature on aerobic respiration would have made marine animals more vulnerable to ocean warming events during periods of limited surface oxygenation. Here, we evaluate the hypothesis that changes in surface oxygenation exerted a first-order control on extinction rates through the Phanerozoic using a combined Earth system and ecophysiological modeling approach. We find that although continental configuration, the efficiency of the biological carbon pump in the ocean, and initial climate state all impact the magnitude of modeled biodiversity loss across simulated warming events, atmospheric oxygen is the dominant predictor of extinction vulnerability, with metabolic habitat viability and global ecophysiology type extinction exhibiting inflection points around 40% of present atmospheric oxygen. Given this is the broad upper limit for estimates of early Paleozoic oxygen levels, our results are consistent with the relative frequency of high-magnitude extinction events (particularly those not included in the canonical Big Five mass extinctions) early in the Phanerozoic being a direct consequence of limited early Paleozoic oxygenation and temperature-dependent hypoxia responses.

22

# 23 Significance Statement

The decline in extinction rates through geologic time is a well-established but enigmatic feature of the marine animal fossil record. We hypothesize that this trend is driven largely by secular changes in the oxygenation of the atmosphere and oceans, as physiological principles predict that marine animals would have been more vulnerable to ocean warming during intervals of geological time with limited atmospheric oxygenation. We
test this at a global oceanographic scale by combining models of ocean biogeochemistry and animal physiology. We show that atmospheric oxygen exerts a first-order control on the simulated extinction vulnerability of marine animals, highlighting its likely importance in controlling extinction trends through geologic time.

33

## Introduction

One of the most striking and poorly understood features of the marine fossil record is a long-term decline in apparent extinction rates through the Phanerozoic (1–7). First revealed by analyses of stratigraphic range compilations ~40 years ago (1), this decline has been a persistent feature in subsequent analyses using expanded databases (8) and improved extinction rate and richness metrics since then (e.g. (2–7); Figs. 1, S1). Declining genus-level extinction rates are a consistent feature of Phanerozoic analyses using a broad range of different richness and extinction rate metrics ((7); Fig. S2, Table S1), and are not substantially altered by varying temporal binning or controlling for species/genus ratio (Figs. S1-2, Table S1; Materials and Methods). In particular, the Cambrian and Ordovician Periods stand out as intervals of especially high faunal turnover ((6, 7, 9); Figs. S1-2), to the point that the early Paleozoic is often treated separately from the rest of the Phanerozoic in quantitative analyses of extinction (e.g. (7, 10, 11)). Changes in the distribution of lineage ages and geographic range sizes, area of continental shelf environments, intensity of species interactions, frequency of geological triggers, stabilization of Earth's climate, and rock-area sampling biases have all been proposed as drivers of this trend (10–14). However, despite the long history of study and numerous proposed explanations, no consensus has been reached regarding the drivers of much higher early Paleozoic extinction rates vs. those of the later Paleozoic, Mesozoic and Cenozoic.

Secular trends in the oxygenation of the ocean and atmosphere are an underexplored factor that may help explain changing extinction intensity through time. A growing body of geochemical proxy evidence suggests that Earth's surface did not become oxygenated to near-modern levels until the early-mid Paleozoic (Fig. 1) (15–20)
and the most recent generation of long-term carbon cycle models support only limited early Paleozoic oxygenation (21–23). While early Paleozoic marine oxygen concentrations were self-evidently above minimum thresholds for aerobic metazoan metabolism (at least, in the shallow ocean), modern observations clearly show that variations in dissolved oxygen exert a strong control on marine biodiversity (24–26).

Moreover, experimental respirometry, metabolic theory, and ecophysiological modeling have demonstrated that the synergistic effects of oxygen and temperature on capacity for aerobic respiration significantly influence the metabolic viability of marine habitat (27–32). Following the principle of oxygen and capacity limited thermal tolerance (OCLTT – see Pörtner (28) for review), the critical partial pressure of oxygen (pO₂) required to sustain resting metabolism increases exponentially with temperature, and the increase in oxygen supply required for ecologically sustainable populations results in a further steepening of this exponential relationship (Fig. 2; (33)). Meta-analyses investigating the sensitivity of modern marine ectotherms to climate-related stressors and the vulnerability of their fossil representatives to extinction further suggest that physiological traits of modern marine animals can be generally scaled to extinction sensitivity to ancient climate perturbations (34).

Seawater pO₂ would have been dramatically lower in an early Paleozoic Earth system with limited surface oxygenation, leading to the prediction that marine ectotherms would have been substantially more vulnerable to temperature-dependent hypoxia, and therefore climate-driven extinction, due to reduced thermal safety margins in poorly oxygenated ancient oceans (Fig. 2; (35)). This logic further leads to the implication that dissolved oxygen was likely important to the structure of marine ecosystems at concentrations well above critical hypoxic thresholds that have previously been invoked for the evolution of early animals (36, 37), as is observed at regional scales in modern marine environments (29–31, 38). Recent advances in Earth system modeling provide an opportunity to test this hypothesis at the planetary-scale and to explicitly distinguish the role of changes in surface oxygenation through the Phanerozoic from variations in continental configuration, biological carbon cycling and climate state.
Our approach is two-fold – incorporating both Earth system and ecophysiological modeling. First, using an Earth system model of intermediate complexity (cGENIE (39) – see *Materials and Methods* and *Supplementary Information*) and an ensemble experimental approach, we generate a suite of three-dimensional realizations of potential marine environmental conditions. These Earth system model experiments span a range of scenarios for atmospheric O₂ and CO₂ concentrations (Fig. 3a) in order to address uncertainty associated with reconstructions of these variables through geologic time (e.g. Fig 1b). Atmospheric CO₂ concentrations are varied to force global climate in our experiments and solar forcing is kept constant for each ensemble (Fig. 3). Simulating the potential three-dimensional structure of past oceans in these experiments is critical: in the upper ocean, oxygen and temperature both decrease with depth, but have opposite gradients with respect to latitude (33, 38). Temperature gradients decrease towards the poles, while subsurface oxygen gradients increase. Because aerobic habitability for marine animals involves the synergistic combination of oxygen and temperature (29, 34, 35, 38), elucidating the spatial complexity of the environmental landscape of ancient oceans is required to appropriately understand biological responses to environmental perturbations.

We additionally identify key boundary conditions that have significantly changed through Earth history and may offer alternative explanations for the trends in extinction rates observed in the fossil record. 1) Continental configuration, which affects latitudinal habitat (i.e. continental shelf) distribution, as well as ocean circulation (and resulting temperature and oxygen gradients), and could therefore have modulated the susceptibility of animals to extinction (33, 40, 41). 2) The efficiency of the biological carbon pump in the ocean, which exerts a control on the distribution of oxygen in the upper water column (42) and has been previously identified as important in differentiating earlier Paleozoic vs. Mesozoic-Cenozoic ocean redox conditions (17). 3) Climate state, which directly impacts ocean surface temperatures (and by extension dissolved marine oxygen) and therefore influences how close animals at a given latitude are living to their physiological limits (29, 38).
Utilizing model ocean temperatures and dissolved oxygen concentrations from this suite of 3D realizations, we then employ the recently developed Metabolic Index (29, 33, 38) in order to evaluate the proportion of modern marine ectotherms (based upon available experimental respirometry data) that could survive in these simulated environments. While we do not attempt to directly simulate taxonomic diversity, the Metabolic Index provides us with a means to mechanistically assess the physiological ecotypes (herein, ecophysiotypes) that can contribute to biodiversity. We apply a probabilistic representation of temperature-dependent hypoxia responses (Fig. 2) to define ‘metabolic habitat viability’, a global metric describing the percentage of feasible ecophysiotypes that can live anywhere in the shelf environments of our 3D ocean models. Following other paleoecological modeling approaches (e.g. (41)), we focus our analyses on shelf environments because they constitute the majority of the Phanerozoic fossil record and host the majority of marine biodiversity in the modern, although we provide complimentary simulations for the entire global ocean.

Finally, we estimate global extinction sensitivity across a range of atmospheric O₂ concentrations to link our ecophysiological models more directly to the trends we observe in the geologic record. We model ‘global ecophysiology extinction’ as the predicted loss of ecophysiotypes across a nominal warming event of ~5°C (as defined by mean temperature change in the equatorial surface ocean; Fig. S3, Table S2). These warming events are represented by Earth system models separated by two doublings of atmospheric CO₂, and are designed to mimic relatively extreme (but not infrequent) examples of climate variability known to have occurred throughout Earth history (43, 44). We performed this extinction sensitivity estimate for all scenarios of atmospheric O₂ and for each ensemble (continental configuration, biological pump efficiency; Figs. 3, S3) across multiple baseline climate states. We further explore the consistency of our results across different thresholds for minimum habitat area (number of equal area ocean model cells) and tolerances for proportional habitat loss, as well as across deviations from the distributions of physiological tolerances observed in modern marine ectotherms. Thus, by combining ecophysiology and ensemble Earth system modeling, we evaluate the vulnerability of marine ectotherms to warming events at a range of atmospheric O₂ levels
146 in a framework that is neutral to the selective advantage of different ecophysiotypes and
147 explores the relative impacts of other aspects of Phanerozoic Earth system evolution on
148 marine animal ecophysiology.

149

150 **Results**

151 *Metabolic habitat viability*

152 Metabolic habitat viability – the total number of ecophysiotypes that can live in
153 the shelf environments of each ancient ocean realization – is shown in Figure 4 for our
154 ensemble of simulations employing an Ordovician continental configuration and Figure
155 S4 for all other Earth system model ensembles. In all ensembles, metabolic habitat
156 viability exhibits a positive relationship with atmospheric O₂ and a negative relationship
157 with global temperature, resulting in a sigmoidal relationship between atmospheric O₂
158 and metabolic habitat viability for a given climate state. While our primary analyses
159 focus on shelf environments, very similar trends are observed when the entire ocean is
160 considered as potential habitat for marine animals (Figs. S5-6). Using entire 3D ocean
161 models often results in slightly higher estimates of metabolic habitat viability because
162 cooler and (at higher atmospheric O₂ levels) well-oxygenated deep ocean environments
163 may act as refugia, despite being unlikely diversity hotspots for other reasons such as
164 food supply (45, 46). The latitudinal distribution of viable ecophysiotypes is similar
165 between continental configurations (Figs. 3, S3), with more ecophysiotypes generally
166 able to inhabit cooler polar regions of the surface ocean. Identifying oceanographic and
167 tectonic drivers of regional differences in viable ecophysiotype biogeography between
168 continental configurations is an interesting avenue for future research, but this would
169 require tuning model ensembles to precisely equivalent climate states and is therefore
170 beyond the scope of this study.

171 Changes in metabolic habitat viability between cooler and warmer climate states
172 allow us to begin linking ecophysiology to extinction sensitivity under different scenarios
173 of atmospheric oxygenation (Fig. 4). The percentage change in metabolic habitat viability
174 between colder and warmer ocean simulations increases as atmospheric oxygen decreases
175 from modern levels (100% preindustrial atmospheric levels, PAL), reaching maximum
176 values at atmospheric O₂ concentrations reconstructed for the early Paleozoic (<40-50%
177 PAL - grey envelope in Figs. 4, S4). Furthermore, because extinction is expected to track
178 the proportional rather than absolute loss of ecophysiotypes (33, 40), differences in
179 metabolic habitat viability would therefore be amplified to result in high levels of
180 warming-driven extinction at lower levels of atmospheric oxygenation reconstructed for
181 the early Paleozoic (e.g. 20% PAL). In such cases, the number of viable ecophysiotypes
182 in any low-O₂ scenario is much lower, leading to greater proportional change despite a
183 similar absolute change in metabolic habitat viability to that predicted at higher O₂ levels
184 (e.g. 80% PAL). Our analyses thus demonstrate that the expected proportional loss of
185 ecophysiotypes during a global warming event would be substantially higher under the
186 limited surface oxygenation scenarios hypothesized for the early Paleozoic than for
187 modern (or, more broadly, post-mid-Paleozoic) levels of oxygenation.

188

189 *Sensitivity of global extinction simulations*

190 When we simulate the impacts of warming events on the extinction of marine
191 animals, the general trend of greater predicted extinction at lower atmospheric oxygen
192 levels is a common feature of all ensemble experiments (Figs. 5, S7-12). There is
193 variation in the absolute magnitude of simulated global ecophysiology extinction between
194 the different continental configurations tested in our ensembles. However, we do not
195 observe any time-dependent trends in extinction sensitivity between progressive
196 continental reconstructions and the most striking feature of the results remains the non-
197 linear increase in global ecophysiology extinction with declining oxygenation (Figs. 5,
198 S7-12). Below 10% PAL O₂, the relationship between modeled extinction and
199 atmospheric oxygen is less predictable. However, the initial pool of viable pre-extinction
200 ecophysiotypes in these simulations is very small at low O₂ levels (Figs. 3, 4, S4-6),
201 suggesting that this inconsistency may partially be a product of sample size (as supported
202 by wider envelopes in Figures 5, S7-8, S11) and the limited experimental respirometry
203 data for animals that live in low-oxygen environments (47). We find that modeled
204 extinction is high at low atmospheric oxygen levels for the Ordovician configuration
205 regardless of the efficiency of the biological pump (Fig. S7), suggesting that despite its
206 role in controlling the distribution of oxygen in the upper water column this factor is of
207 secondary importance compared to atmospheric O₂. A shallower remineralization depth
208 does slightly reduce simulated extinction sensitivity at low atmospheric O₂ (~5-20%
209 PAL), although the same caveat of small initial ecophysiology population size applies.

210 Increasing the baseline temperature in our extinction simulations (i.e. testing for
211 any initial climate state dependence of extinction susceptibility), results in increased
212 predicted global ecophysiology extinction for all three continental configurations (Fig.
213 S8). A corresponding decrease in modeled ecophysiology extinction is generally observed
214 at cooler baseline temperatures. Cool, low O₂ (<20% PAL) Ordovician simulations are a
215 minor exception to this general rule, suggesting that continental configuration may exert a
216 second-order control on the climate state dependence of extinction vulnerability. The
217 impact of initial climate state on simulated ecophysiology extinction appears less
218 dramatic than the impact of atmospheric oxygen, particularly in light of the apparent
219 inflection in modeled extinction observed around upper estimates for early Paleozoic O₂
220 levels. However, we note that pre-extinction climate state impacts modeled
221 ecophysiology extinction more than any other Earth system boundary condition explored
222 in this study.

223 Altering the assumptions of our ecophysiological modeling approach does not
224 dramatically impact the observed relationship between predicted extinction and
225 atmospheric oxygen. In our primary analyses, we define ecophysiology extinction as only
226 occurring when no cells in our simulated shelf environments are viable for an
227 ecophysiology, or (equivalently) when 100% of metabolically viable habitat area is lost.
228 Increasing the minimum number of inhabited ocean cells required to survive a warming
229 event or decreasing the proportional habitat loss required for extinction moderately
230 increases modeled global ecophysiology extinction (Figs. S9-10). The differences in the
231 shape of these responses highlight the importance of these parameters and the second-
order effects of continental configuration. However, the overall trend of decreasing extinction intensity with increased atmospheric O₂, with an apparent inflection point generally around 40% PAL, is consistent across all of these simulations. Similar trends are observed when investigating the impacts of varying the distributions used to parameterize physiological responses (A₀ – the inverse of the hypoxic threshold; and E₀ – the temperature-dependency of the hypoxic threshold). While varying the distributions of A₀ and E₀ does impact the magnitude of simulated global ecophysiotype extinction and the steepness of the inflection around 40% PAL O₂ (Fig. S11), the overarching trends of increasing predicted extinction with decreasing atmospheric oxygen (down to at least 20% PAL O₂) are not impacted. The size of the synthetic population of ecophysiotypes (*Materials and Methods*) has no noticeable effect on the analyses illustrated in Figure 5 (Fig. S12).

244

## Discussion

Using a coupled Earth system and ecophysiological modeling approach, our analyses suggest that oxygenation of the shallow ocean (as modulated by atmospheric oxygen levels) is the dominant factor in governing the vulnerability of marine animals to extinction during warming events. These results remain consistent when complex oceanographic relationships between temperature and seawater oxygenation and uncertainties in ocean circulation, biological pump efficiency, climate state, and habitat area thresholds are taken into account. Critically, the relationship we observe between simulated global ecophysiotype extinction and atmospheric oxygen is non-linear, with an apparent inflection point around 40% PAL O₂ (Fig. 5). This broadly coincides with the upper limit of estimates for early Paleozoic atmospheric oxygen (Fig. 1; (15, 17)), supporting the hypothesis that relatively low atmospheric oxygen levels created boundary conditions under which high-magnitude extinction events were more likely to occur as a product of climate variability.
While atmospheric oxygen is the primary control on extinction sensitivity in our simulations, our analyses show that other factors modulate the response to varying degrees. This is unsurprising considering how these factors affect the three-dimensional oxygen-temperature landscape of the ocean. In particular, our analyses indicate that continental configuration may play an important secondary role in governing ecophysiology of extinction vulnerability at the global scale (Fig. 5), although this is difficult to clearly separate from minor climatic differences between ensembles. Notably, we do not observe time-correlative trends between the three continental configurations that could help explain the secular extinction rate observations from the fossil record. In other words, we do not find that extinction sensitivity in the polar-biased and relatively compact continental configuration of the Ordovician is particularly different from that of the more evenly latitudinally distributed and dispersed configuration of the Paleocene when compared to the differences associated with changing atmospheric O₂. The efficiency of the biological pump (simulated here as a difference in prescribed remineralization depths in the ocean) is also thought to have changed through the Phanerozoic (17). However, we find that the differences in simulated extinction sensitivity are minimal and (if anything) in opposition to the influence of changes in atmospheric oxygenation (Fig. S7). This is likely because both continental configuration and carbon remineralization depth mainly change oxygen distributions in the deeper ocean, whereas most animals are living (and dying) in the surface ocean, which is closely coupled to the O₂ concentration of the atmosphere. Finally, while the habitat area thresholds we used to define simulated extinction and specific parameterization of physiological traits do impact the slope and magnitude of the observed relationship between atmospheric oxygen and extinction (in expected directions; Figs. S9-11), the first order trends remain consistent regardless of how we parametrize our ecophysiological extinction model.

We additionally find that an important climate state-dependence exists for ecophysiological extinction sensitivity. Marine animals are expected to be living closer to their thermal limits as defined by oxygen and capacity limited thermal tolerance in warmer environments based on observations from modern tropical marine habitats (Fig. 5).
289 2; (29, 38)), with implications for broader areas of the ocean in ancient greenhouse climates. We provide further model support for this hypothesis by demonstrating increased simulated extinction magnitudes in warmer climate states (Fig. S8). Very warm 290 early Paleozoic sea surface temperatures have been reconstructed using geochemical 291 paleothermometry (48–50), suggesting that this phenomenon could also provide some 292 explanatory power for the frequency of high magnitude extinctions. However, while 293 recent advances in paleothermometry have shown considerable potential in constraining 294 Ordovician climate dynamics (50), there is still currently limited consensus on directional 295 trends in early Phanerozoic sea surface temperatures, especially in the earliest Paleozoic 296 (51–53). In contrast, reconstructions of atmospheric oxygen are increasingly consistent, 297 even if the absolute concentrations are poorly constrained (Fig. 1). This discussion should 298 not be framed as a question of oxygen *versus* temperature; inherent to the concept of 299 temperature-dependent hypoxia is that both interact to determine aerobic safety margins 300 (Fig. 2; (29, 34, 38)). We therefore suggest that atmospheric oxygen can confidently be 301 established as a first-order control on early Phanerozoic extinction rates, and that 302 potentially high sea surface temperatures through the same time interval may have 303 amplified the effect of limited surface oxygenation. Changes in the latitudinal focus of 304 fossil sampling through the Phanerozoic may also bias extinction rate estimates slightly 305 towards environments with warmer baseline temperatures in the early Paleozoic when 306 well-sampled paleocontinents were mainly situated at low latitudes (54). However, the 307 lack of similarly high long-term extinction rates across other intervals with low-latitude 308 sampling biases indicates that this cannot fully explain the trends we seek to explain here 309 and that atmospheric oxygen (potentially compounded by warm global climate) is a more 310 probable first-order driver. 311

312
The plausibility of temperature-dependent hypoxia as a driver of high early 313 Paleozoic extinction rates is supported by links between sea-level changes, ocean 314 deoxygenation and extinctions during this interval of geologic time. The Cambrian- 315 Ordovician biomere-style extinction events are a major component of the 316 characteristically high early Paleozoic extinction rates (9). These events, predominantly 317 observed in the trilobite fossil record (e.g. (55)), have long been recognized as associated 318
with sedimentological evidence of sea-level change, possible thermocline migration and ocean anoxia (56, 57). Stable isotope geochemistry provides further evidence that the biomere events were linked to major marine carbon and sulfur cycle perturbations, consistent with ocean warming and deoxygenation (58–61). Similar inferences of ocean anoxia and rapid sea-level change have been made for the early-middle Cambrian Botomian extinction event (62). In contrast, multiple Mesozoic and Cenozoic hyperthermals (including the Paleocene-Eocene Thermal Maximum and Cretaceous Oceanic Anoxic Event 2) that exhibit similar broad patterns of environmental change (43, 63–65) do not correlate with global marine extinctions of comparably high magnitude (Fig. 1; (9)). Our analyses demonstrate that the effects of changing atmospheric oxygenation on the thermal safety margins of marine ectotherms (Fig. 2) may account for this discrepancy in extinction magnitude between events with similar signatures of environmental change. Our simulations of metabolic habitat viability also offer a potential mechanism for the high origination rates observed through the early Paleozoic (6, 7) – particularly if simple logistic models of extinction and subsequent recovery (i.e. origination; e.g. (4, 66)) can be broadly applied to the Phanerozoic fossil record. However, we focus upon warming driven extinction events because fewer assumptions regarding the evolutionary advantages of different ecophysiotypes (and how they relate to the rate at which physiological diversity increases in cool climate states) are required for these models, which simply simulate the instantaneous loss of viable ecophysiotypes.

Our combined Earth system and ecophysiological modeling approach provides a methodological advance in moving from simple correlation to more mechanistic understanding of how environmental change impacted marine ecosystems through Earth history. The use of physiology as a conceptual bridge between the fossil and geochemical records is of clear appeal for establishing a mechanistic understanding of animal-environment interactions from the geologic record (67). However, the power of Earth system models in simulating the ecophysiological responses of marine animals to deep-time environmental change is only recently being explored (e.g. (33, 41)). We apply simple assumptions about the diversification of ecophysiotypes (i.e. that ecophysiotypes proliferate in times of stability and become inviable during warming events) to test the
hypothesis that limited atmospheric oxygenation would increase physiological vulnerability to ocean warming events on geologic timescales. For increasingly accurate models of animal-environment interactions through the Phanerozoic (for example, ones that directly simulate latitudinal gradients in taxonomic diversity), other key variables involved in these interactions will need to be incorporated. At the organism-scale, physiological traits including vulnerability to ocean acidification, food supply limitation and possible variations in physiological responses are all likely important variables. For example, the different buffering capacities of marine animal groups have previously been invoked to explain differential extinction selectivity (68, 69). Meta-analytical approaches also suggest that sensitivity to climate related stressors correlates with genus survival across ancient extinction events in higher-level taxonomic groups (34), further indicating that differences in physiological tolerances between marine clades may be an important consideration when applying ecophysiological principles to Earth history. While we demonstrate that our results are consistent across a range of deviations from mean modern physiological responses, it will hopefully become feasible to begin linking physiological traits more directly to extinct animal groups based on phylogenetic and adaptive principles as the availability of experimental respirometry data characterizing temperature-dependent hypoxic responses of marine ectotherms increases. At the macroecological scale, migration (as applied in Saupe et al. (41)), species interactions, and resource competition also play key roles in governing population and species responses to environmental change. As model frameworks for deep-time animal-environment interactions develop in terms of both ecophysiological and biogeochemical realism, the temporal-spatial resolution of paleontological and geochemical databases will also be critical to the explanatory power of these methodologies, as will techniques addressing the limitations and biases inherent to the geologic record (4–6, 13, 54, 70–72).

374

## 375 Conclusions

Using a novel ensemble Earth system and ecophysiological modeling framework, we demonstrate that atmospheric oxygen levels likely exerted a first-order control on
extinction vulnerability through the Phanerozoic. Our model analyses illustrate that the
theoretical predictions of oxygen and capacity limited thermal tolerance are expected to
have globally significant implications, resulting in dramatically increased extinction of
marine animals under varying climatic conditions during geological periods with limited
surface oxygenation. While warmer initial climate states increase the magnitude of our
simulated extinctions, and continental configuration and the strength of the biological
pump also have minor impacts on our predictions, the trend of increased extinction
vulnerability with early Paleozoic levels of atmospheric oxygen is dominant across all of
our simulations. We therefore argue that the exceptional frequency of high-magnitude
extinction events in the early Paleozoic was primarily a consequence of limited surface
oxygenation and temperature-dependent hypoxia responses in marine animals.

389

390 **Materials and Methods**

391 *Quantitative paleobiology*

We performed a series of quantitative paleobiological analyses to evaluate evidence for declining extinction rates through the Phanerozoic. We adapted methods from Kocsis et al. (7) to generate new reconstructions and statistical analyses of extinction rates through the Phanerozoic specifically for marine ectotherms, using an updated download from the Paleobiology Database (PBDB) and omitting endothermic and/or exclusively air-breathing higher-level marine taxa (marine mammals, marine reptiles and turtles). We generated reconstructed extinction rates using three richness metrics (classical rarefaction, raw occurrences and shareholder quorum subsampling) and four extinction rate metrics (second-for-third substitution, corrected three-timer, gap-filler and per-capita). We applied these metrics to generate reconstructions for all genera both at stage resolution and using the PBDB ~10 million year bins. To investigate the potential impact of changing species to genus ratios through the Phanerozoic (e.g. (73)), we further conducted the same analyses on filtered datasets, generating reconstructions all genera, genera with a single accepted named species, genera with 2 species, genera with 3
species, and genera with 4 or more accepted named species (Fig. S1). For shareholder
quorum subsampling analyses we used a quorum level of 0.5 to accommodate lower
sampling intensity in filtered datasets (as supported by Boag et al. (32)) and remain
within the >0.4 quorum range suggested by previous authors (7). Figure 1 is a summary
of the top-right panel of Figure S1, including all genera at stage resolution. We further fit
cross-validated LOESS models to our reconstructed extinction rates to evaluate whether
early Paleozoic decreases in extinction rates were a major feature of the full Phanerozoic
record (Fig. S2). Finally, we generated Spearman's correlation coefficients to investigate
the statistical significance of temporal declines in extinction rates across all of our
treatments (Table S1). Adapted code from Kocsis et al. (7) and filtered Paleobiology
Database download are available at https://github.com/richardstockey/cGENIE-
metabolic_index.extinction and are assigned a DOI: 10.5281/zenodo.5143583

418

419 *Earth system modeling*

420 cGENIE is an Earth system model of intermediate complexity, designed for the spatially-
421 explicit biogeochemical characterization of global-scale paleoceanographic problems
422 (39). We create a range of 3D realizations of ocean environmental conditions by running
423 ensembles of model experiments. The experiments in each ensemble are run for 10,000 or
424 20,000 years (depending on the time required for dissolved marine O2 and seawater
425 temperature to reach steady state), with a range of assumptions regarding atmospheric
426 pCO2 specific to each continental configuration based upon changing solar luminosity
427 and albedo (Fig. 3), and a range of atmospheric pO2 scenarios based upon the range of
428 estimates from geochemical proxy data and biogeochemical box modeling (Fig. 1). We
429 apply this ensemble treatment to three ancient continental configurations, chosen for their
430 relatively even temporal spread across the Phanerozoic and history of use in previous
431 paleoclimate studies – a Paleocene configuration (~55 Ma) (43), a Permian configuration
432 (~251 Ma) (74), and an Ordovician configuration (~450 Ma) (based on the 8 PAL CO2
433 FOAM simulation of Pohl et al. (75) – see Supplementary Information). For each of these
434 ancient configurations, we apply age-appropriate solar luminosity (76). Finally, in the
435 Ordovician continental configuration, we explore the possible influence of postulated

436 Phanerozoic changes in organic matter remineralization depth in the ocean (17, 42) –

437 testing ~1/3 modern, or specifically 200m, e-folding depth (42).

438 We extract the mean annual ocean temperature and ocean oxygen concentration

439 during the last simulated year for each ocean grid cell and generate corresponding

440 seawater oxygen partial pressure ($pO_2$) values (following the methods of Hofmann et al.

441 (77)). Changing solar luminosity and continental configuration results in different mean

442 sea surface temperatures for the same atmospheric $CO_2$ assumption. We therefore

443 identify the $pCO_2$ values for each continental configuration that create comparable

444 climate simulations (Table S2 and *Supplementary Information*). These climate scenarios

445 are compared using the mean equatorial (3.2°S to 3.2°N, 0-80m depth) surface ocean

446 temperature for each model (Figs. 4, S4-6, Table S2).

447 The code for the version of the ‘muffin’ release of the cGENIE Earth system
448 model used in this paper is tagged as v0.9.19 and is assigned a DOI:

449 10.5281/zenodo.4473048. Configuration files for the specific experiments presented in

450 the paper can be found in the directory: genie-userconfigs/MS/stockeyetal.PNAS.2021.

451 Details on the experiments, plus the command line needed to run each one, are given in

452 the readme.txt file in that directory. All other configuration files and boundary conditions

453 are provided as part of the code release.

454 A manual detailing code installation, basic model configuration, tutorials covering

455 various aspects of model configuration and experimental design, plus results output and

456 processing, is assigned a DOI: 10.5281/zenodo.4469678.

457

458 *Ecophysiological simulations of animal-environment interactions*

459 We combine our ensemble Earth system modeling experiments with a probabilistic

460 ecophysiological model of extinction vulnerability. The Metabolic Index, as defined by

461 Deutsch et al. (29) and Penn et al. (33), is configured here to calculate the proportion of
modern marine ectotherms (based on described temperature-dependent hypoxia responses) that can inhabit each ocean cell in the cGENIE Earth system model.

$$ \text{Metabolic habitat viability} = \max_{\substack{A_o, E_o, \phi_{crit} \\ \text{min } A_o, E_o, \phi_{crit}}} \sum_{A_o} \frac{pO_2}{\exp\left[-\frac{E_o}{k_B}\left(\frac{1}{T} - \frac{1}{T_{\text{ref}}}\right)\right]} > \phi_{crit} \quad \text{(Eqn. 1)} $$

At the individual organism scale, $A_o$ is the inverse of the hypoxic threshold of the organism (the minimum required seawater $pO_2$ ($P_{crit}$) to sustain resting aerobic metabolism), $E_o$ is the temperature-dependency of the hypoxic threshold and $\phi_{crit}$ is the multiplicative increase in oxygen supply that is required to support ecologically sustainable populations (calculated based on biogeographic distribution, relative to experimental respirometry measurements of resting metabolism). In our modeling approach, these variables are parametrized by probability density functions based on observations from laboratory experiments and species distributions (Fig. 2; following Penn et al. (33)). Seawater temperature and $pO_2$ are environmental variables generated through the Earth system modeling approach described above (Fig. 3). All other variables are constant ($k_B$ is the Boltzmann constant, $T_{ref}$ is a reference temperature of 15°C). In our analyses, we simulate 1000 ecophysiotypes, sampled from probability distributions of $A_o$, $E_o$ and $\phi_{crit}$, and map which cells (if any) in each cGENIE ocean model are habitable for each ecophysiology type (see Figure S12 for exploration of the number of ecotypes used).

Metabolic habitat viability, as defined here, quantifies the proportion of these modeled ecophysiotypes that can live in the non-polar shelf environments (defined as any cell adjacent to continental land mass in the top 3 layers of the cGENIE ocean, or <283.8m water depth) of each global ocean model. See Figures S5-6 for sensitivity analyses using the entire global ocean model, rather than shelf environments alone. Polar environments (>70°N and <70°S) are excluded from all analyses due to the coarse model resolution of polar ocean environments (e.g. in a Paleocene northern polar ocean) in cGENIE and to reduce local non-linear responses to sea ice extent.
We further simulate extinction as the global loss of ecophysiotypes between two CO₂ scenarios at the same O₂ level. Specifically, we evaluate the proportional loss of ecophysiotypes between baseline climate states (relatively cool, approximately pre-industrial mean equatorial surface ocean temperatures – *Supplementary Information*) and equivalent model scenarios that have undergone ~5°C equatorial warming (Table S2; modeling a nominal hyperthermal event of roughly the same magnitude as the Paleocene-Eocene Thermal Maximum (43)) at each atmospheric O₂ level. Global ecophysiotype extinction is thus defined as the proportion of ecophysiotypes that were viable in shelf environments of the colder baseline simulation but are no longer viable in any continental shelf cell of the warmer ocean simulation. Sensitivity analyses exploring the number of viable cells required to survive a simulated warming event, and the proportional habitat loss required to drive extinction, are shown in Figures S9-10. Additional sensitivity analyses exploring deviations from the distributions of A₀ and E₀ defined in Penn et al. (33) are shown in Figure S11. In these analyses, “low” and “high” distributions have mean values defined by the 25th and 75th percentiles of the distributions used in our primary analyses, and half of the variance of those distributions. In plots of global ecophysiotype extinction in shelf environments, we summarize the distribution of results from 100 different simulated populations of 1000 randomly sampled ecophysiotypes to illustrate the range of sampling-related uncertainty in our probabilistic approach to ecophysiological modeling. All plots of single ecotype populations (i.e. those with points rather than envelopes – Figs. 4, S4-6) use the same randomly sampled population.

R scripts to reproduce the ecophysiological analyses and associated figures are available at https://github.com/richardstockey/cGENIE-metabolic_index.extinction and are assigned a DOI: 10.5281/zenodo.5143583

## Acknowledgements

We thank Thomas Boag and Justin Penn for helpful discussion. R.G.S. and E.A.S. acknowledge support from National Science Foundation grant EAR1922966. R.G.S. acknowledges support from the NASA Astrobiology Institute Early Career Collaboration Award. A.R acknowledges support from the Heising-Simons Foundation. This project
517 has received funding from the European Union's Horizon 2020 research and innovation

518 programme under the Marie Sklodowska-Curie grant agreement No. 838373. cGENIE

519 simulations were performed on Earth system modelling cluster facility at University of

520 California Riverside. Ecophysiological analyses of cGENIE models were performed on

521 the Sherlock cluster at Stanford University. We thank Stanford University and the

522 Stanford Research Computing Center for providing computational resources and support

523 that contributed to this research.

524

525

526

527
528 References

529 1. D. M. Raup, J. J. Sepkoski, Mass extinctions in the marine fossil record. *Science* **215**, 1501–1503
530 (1982).

531 2. M. Foote, Origination and extinction components of taxonomic diversity: general problems. *Paleobiology* **26**, 74–102 (2000).

532 3. M. Foote, Origination and extinction through the Phanerozoic: A new approach. *J. Geol.* **111**, 125–148 (2003).

533 4. J. Alroy, Dynamics of origination and extinction in the marine fossil record. *Proc. Natl. Acad. Sci. U. S. A.* (2008) https://doi.org/10.1073/pnas.0802597105.

534 5. J. Alroy, A more precise speciation and extinction rate estimator. *Paleobiology* **41**, 633–639 (2015).

535 6. B. Kröger, F. Franeck, C. M. Ø. Rasmussen, The evolutionary dynamics of the early Palaeozoic marine biodiversity accumulation. *Proc. R. Soc. B Biol. Sci.* **286** (2019).

541 7. Á. T. Kocsis, C. J. Reddin, J. Alroy, W. Kiessling, The r package divDyn for quantifying diversity dynamics using fossil sampling data. *Methods Ecol. Evol.* **10**, 735–743 (2019).

543 8. J. Alroy, et al., Effects of sampling standardization on estimates of Phanerozoic marine diversification. *Proc. Natl. Acad. Sci.* **98**, 6261–6266 (2001).

545 9. R. K. Bambach, A. H. Knoll, S. C. Wang, Origination, extinction, and mass depletions of marine diversity. *Paleobiology* **30**, 522–542 (2004).

547 10. S. Finnegan, J. L. Payne, S. C. Wang, The Red Queen revisited: reevaluating the age selectivity of Phanerozoic marine genus extinctions. *Paleobiology* **34**, 318–341 (2008).

549 11. J. L. Payne, S. Finnegan, The effect of geographic range on extinction risk during background and mass extinction. *Proc. Natl. Acad. Sci. U. S. A.* **104**, 10506–10511 (2007).

551 12. S. E. Peters, M. Foote, Determinants of extinction in the fossil record. *Nature* **416**, 420–424 (2002).

552 13. S. E. Peters, Environmental determinants of extinction selectivity in the fossil record. *Nature* **454**, 626–629 (2008).

554 14. J. L. Payne, A. Bachan, N. A. Heim, P. M. Hull, M. L. Knope, The evolution of complex life and the stabilization of the Earth system. *Interface Focus* (2020) https://doi.org/10.1098/RSFS.2019.0106.

557 15. E. A. Sperling, et al., Statistical analysis of iron geochemical data suggests limited late Proterozoic oxygenation. *Nature* **523**, 451–454 (2015).

559 16. T. W. Dahl, et al., Devonian rise in atmospheric oxygen correlated to the radiations of terrestrial plants and large predatory fish. *Proc. Natl. Acad. Sci. U. S. A.* **107**, 17911–17915 (2010).

561 17. W. Lu, et al., Late inception of a resiliently oxygenated upper ocean. *Science* **361**, 174–177 (2018).

562 18. M. W. Wallace, et al., Oxygenation history of the Neoproterozoic to early Phanerozoic and the rise of land plants. *Earth Planet. Sci. Lett.* **466**, 12–19 (2017).

564 19. D. A. Stolper, C. B. Keller, A record of deep-ocean dissolved O2 from the oxidation state of iron in submarine basalts. *Nature* **553**, 323–327 (2018).

566 20. E. A. Sperling, et al., A long-term record of early to mid-Paleozoic marine redox change. *Sci. Adv.* **7**, eabf4382 (2021).

568 21. A. J. Krause, et al., Stepwise oxygenation of the Paleozoic atmosphere. *Nat. Commun.* **9**, 1–10 (2018).

570 22. T. M. Lenton, S. J. Daines, B. J. W. Mills, COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time. *Earth-Science Rev.* **178**, 1–28 (2018).

572 23. R. Tostevin, B. J. W. W. Mills, Reconciling proxy records and models of Earth's oxygenation during the Neoproterozoic and Palaeozoic. *Interface Focus* (2020) https://doi.org/10.1098/rsfs.2019.0137.

575 24. E. A. Sperling, C. A. Frieder, L. A. Levin, Biodiversity response to natural gradients of multiple stressors on continental margins. *Proc. R. Soc. B Biol. Sci.* **283** (2016).

577 25. E. Sampaio, et al., Impacts of hypoxic events surpass those of future ocean warming and acidification. *Nat. Ecol. Evol.* (2021) https://doi.org/10.1038/s41559-020-01370-3.

579 26. L. A. Levin, "Oxygen Minimum Zone Benthos: Adaptation and Community Response to Hypoxia"
580 in *Oceanography and Marine Biology, An Annual Review, Volume 41*, (CRC Press, 2003), pp. 1–45.

581 27. H. O. Pörtner, A. P. Farrell, Physiology and climate change. *Science* **322**, 690–692 (2008).

583 28. H. O. Pörtner, Oxygen- And capacity-limitation of thermal tolerance: A matrix for integrating climate-related stressor effects in marine ecosystems. *J. Exp. Biol.* **213**, 881–893 (2010).

585 29. C. Deutsch, A. Ferrel, B. Seibel, H. O. Pörtner, R. B. Huey, Climate change tightens a metabolic constraint on marine habitats. *Science* **348**, 1132–1135 (2015).

587 30. E. M. Howard, et al., Climate-driven aerobic habitat loss in the California Current System. *Sci. Adv.* **6**, eaay3188 (2020).

589 31. M. I. Duncan, N. C. James, W. M. Potts, A. E. Bates, Different drivers, common mechanism; the distribution of a reef fish is restricted by local-scale oxygen and temperature constraints on aerobic metabolism. *Conserv. Physiol.* **8** (2020).

592 32. T. H. Boag, W. Gearty, R. G. Stockey, Metabolic tradeoffs control biodiversity gradients through geological time. *Curr. Biol.* (2021) https://doi.org/10.1016/j.cub.2021.04.021.

593 33. J. L. Penn, C. Deutsch, J. L. Payne, E. A. Sperling, Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction. *Science* **362**, eaat1327 (2018).

594 34. C. J. Reddin, P. S. Nätscher, Á. T. Kocsis, H. O. Pörtner, W. Kiessling, Marine clade sensitivities to climate change conform across timescales. *Nat. Clim. Chang.* **10**, 249–253 (2020).

598 35. T. H. Boag, R. G. Stockey, L. E. Elder, P. M. Hull, E. A. Sperling, Oxygen, temperature and the deep-marine stenothermal cradle of Ediacaran evolution. *Proc. R. Soc. B Biol. Sci.* **285**, 20181724 (2018).

601 36. C. T. Reinhard, N. J. Planavsky, S. L. Olson, T. W. Lyons, D. H. Erwin, Earth’s oxygen cycle and the evolution of animal life. *Proc. Natl. Acad. Sci.* **113**, 8933–8938 (2016).

603 37. E. A. Sperling, A. H. Knoll, P. R. Girguis, The Ecological Physiology of Earth’s Second Oxygen Revolution. *Annu. Rev. Ecol. Evol. Syst.* **46**, 215–235 (2015).

605 38. C. Deutsch, J. L. Penn, B. Seibel, Metabolic trait diversity shapes marine biogeography. *Nature* **585**, 557–562 (2020).

607 39. A. Ridgwell, et al., Marine geochemical data assimilation in an efficient Earth system model of global biogeochemical cycling. *Biogeosciences* **4**, 87–104 (2007).

609 40. C. J. Reddin, Á. T. Kocsis, W. Kiessling, Climate change and the latitudinal selectivity of ancient marine extinctions. *Paleobiology* **45**, 70–84 (2019).

611 41. E. E. Saupe, et al., Extinction intensity during Ordovician and Cenozoic glaciations explained by cooling and palaeogeography. *Nat. Geosci.* **13**, 65–70 (2020).

613 42. K. M. Meyer, A. Ridgwell, J. L. Payne, The influence of the biological pump on ocean chemistry: implications for long-term trends in marine redox chemistry, the global carbon cycle, and marine animal ecosystems. *Geobiology* **14**, 207–219 (2016).

616 43. M. Gutjahr, et al., Very large release of mostly volcanic carbon during the Palaeocene–Eocene Thermal Maximum. *Nature* **548**, 573–577 (2017).

618 44. G. L. Foster, D. L. Royer, D. J. Lunt, Future climate forcing potentially without precedent in the last 420 million years. *Nat. Commun.* **8** (2017).

620 45. S. N. C. Woolley, et al., Deep-sea diversity patterns are shaped by energy availability. *Nature* **533**, 393–396 (2016).

622 46. E. A. Sperling, R. G. Stockey, The temporal and environmental context of early animal evolution: Considering all the ingredients of an “explosion.” *Integr. Comp. Biol.* **58**, 605–622 (2018).

624 47. B. A. Seibel, C. Deutsch, Oxygen supply capacity in animals evolves to meet maximum demand at the current oxygen partial pressure regardless of size or temperature. *J. Exp. Biol.* **223**, jeb.210492 (2020).

627 48. T. W. Hearing, et al., An early Cambrian greenhouse climate. *Sci. Adv.* **4**, eaar5690 (2018).

628 49. J. A. Trotter, I. S. Williams, C. R. Barnes, C. Lecuyer, R. S. Nicoll, Did Cooling Oceans Trigger Ordovician Biodiversification? Evidence from Conodont Thermometry. *Science* **321**, 550–554 (2008).

631 50. S. L. Goldberg, T. M. Present, S. Finnegan, K. D. Bergmann, A high-resolution record of early Paleozoic climate. *Proc. Natl. Acad. Sci.* **118**, e2013083118 (2021).

633 51. C. Vérard, J. Veizer, On plate tectonics and ocean temperatures. *Geology* **47**, 882–885 (2019).
634 52. B. J. W. Mills, et al., Modelling the long-term carbon cycle, atmospheric CO2, and Earth surface temperature from late Neoproterozoic to present day. *Gondwana Res.* **67**, 172–186 (2019).

636 53. K. D. Bergmann, et al., A paired apatite and calcite clumped isotope thermometry approach to estimating Cambro-Ordovician seawater temperatures and isotopic composition. *Geochim. Cosmochim. Acta* **224**, 18–41 (2018).

639 54. R. A. Close, R. B. J. Benson, E. E. Saupe, M. E. Clapham, R. J. Butler, The spatial structure of Phanerozoic marine animal diversity. *Science* **368**, 420–424 (2020).

641 55. A. R. Palmer, Trilobites of the Late Cambrian pterocephaliid biomere in the Great Basin, United States. *U. S. Geol. Surv. Prof. Pap.* **493**, 1–104 (1965).

643 56. A. R. Palmer, The biomere problem: evolution of an idea. *J. Paleontol.* **58**, 599–611 (1984).

644 57. R. A. Fortey, There are Extinctions and Extinctions: Examples from the Lower Palaeozoic. *Philos. Trans. R. Soc. B Biol. Sci.* **325**, 327–355 (1989).

646 58. M. R. Saltzman, C. T. Edwards, J. M. Adrain, S. R. Westrop, Persistent oceanic anoxia and elevated extinction rates separate the Cambrian and Ordovician radiations. *Geology* **43**, 807–810 (2015).

649 59. B. C. Gill, et al., Geochemical evidence for widespread euxinia in the Later Cambrian ocean. *Nature* **469**, 80–83 (2011).

651 60. M. Elrick, S. Rieboldt, M. Saltzman, R. M. McKay, Oxygen-isotope trends and seawater temperature changes across the Late Cambrian Steptoean positive carbon-isotope excursion (SPICE event). *Geology* **39**, 987–990 (2011).

654 61. M. A. LeRoy, B. C. Gill, E. A. Sperling, N. R. McKenzie, T.-Y. S. Park, Variable redox conditions as an evolutionary driver? A multi-basin comparison of redox in the middle and later Cambrian oceans (Drumian-Paibian). *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 110209 (2021).

657 62. A. Y. Zhuravlev, R. A. Wood, Anoxia as the cause of the mid-Early Cambrian (Botomian) extinction event. *Geology* **24**, 311–314 (1996).

659 63. W. Yao, A. Paytan, U. G. Wortmann, Large-scale ocean deoxygenation during the Paleocene-Eocene thermal maximum. *Science* **361**, 804–806 (2018).

661 64. H. C. Jenkyns, Geochemistry of oceanic anoxic events. *Geochemistry, Geophys. Geosystems* **11** (2010).

663 65. M. O. Clarkson, et al., Uranium isotope evidence for two episodes of deoxygenation during Oceanic Anoxic Event 2. *Proc. Natl. Acad. Sci.* **115**, 201715278 (2018).

665 66. D. H. Erwin, Lessons from the past: Biotic recoveries from mass extinctions. *Proc. Natl. Acad. Sci. U. S. A.* **98**, 5399–5403 (2001).

667 67. A. H. Knoll, Systems paleobiology. *Bull. Geol. Soc. Am.* **125**, 3–13 (2013).

668 68. A. H. Knoll, R. K. Bambach, J. L. Payne, S. Pruss, W. W. Fischer, Paleophysiology and end-Permian mass extinction. *Earth Planet. Sci. Lett.* **256**, 295–313 (2007).

670 69. J. L. Payne, M. E. Clapham, End-permian mass extinction in the oceans: An ancient analog for the twenty-first century? *Annu. Rev. Earth Planet. Sci.* **40**, 89–111 (2012).

672 70. J. X. Fan, et al., A high-resolution summary of Cambrian to early Triassic marine invertebrate biodiversity. *Science* **367**, 272–277 (2020).

674 71. A. Mehra, et al., Curation and Analysis of Global Sedimentary Geochemical Data to Inform Earth History. *GSA Today* **31**, 4–9 (2021).

676 72. Ú. C. Farrell, et al., The Sedimentary Geochemistry and Paleoenvironments Project. *Geobiology*, gbi.12462 (2021).

678 73. K. W. Flessa, D. Jablonski, Declining Phanerozoic background extinction rates: Effect of taxonomic structure? *Nature* **313**, 216–218 (1985).

680 74. K. M. Meyer, L. R. Kump, A. Ridgwell, Biogeochemical controls on photic-zone euxinia during the end-Permian mass extinction. *Geology* **36**, 747–750 (2008).

682 75. A. Pohl, et al., Glacial onset predated Late Ordovician climate cooling. *Paleoceanography* **31**, 800–821 (2016).

684 76. D. O. Gough, “Solar Interior Structure and Luminosity Variations” in *Physics of Solar Variations*, (Springer Netherlands, 1981), pp. 21–34.

686 77. A. F. Hofmann, E. T. Peltzer, P. M. Walz, P. G. Brewer, Hypoxia by degrees: Establishing definitions for a changing ocean. *Deep. Res. Part I Oceanogr. Res. Pap.* **58**, 1212–1226 (2011).
688 78. C. T. Edwards, M. R. Saltzman, D. L. Royer, D. A. Fike, Oxygenation as a driver of the Great

689 Ordovician Biodiversification Event. *Nat. Geosci.* **10**, 925–929 (2017).

690

691
692 **Figures**

693
<figure><img src="image_2.png" /><figcaption>Figure 1: Trends in atmospheric oxygen and extinction rates through the Phanerozoic. A) Genus-level extinction rates for marine ectotherms (adapted from Kocsis et al. (7), *Materials and Methods*) – line and points represent mean values, envelope represents maximum and minimum estimates from all 12 methodological approaches (3 treatments: raw, classical rarefaction and Shareholder Quorum). Subsampling-standardized, for each of 4 metrics: average per capita rates, corrected three-timer rates, gap-filler rates and second-for-third substitution rates). B) Reconstructions of atmospheric oxygen from biogeochemical models (21, 22, 78) and geochemical proxy records (15, 17).</figcaption></figure>

694

695

696

697

698

699

700

701

702
Metabolic Index parameters:

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Example ecophysiology type from $A_o$, $E_o$, $\Phi_{crit}$ distributions:

<figure><img src="image_6.png" /></figure>

Figure 2: Schematic of ecophysiological model based on the Metabolic Index (Equation 1). The frequency distributions of three key ecophysiological parameters – $A_o$, $E_o$ and $\Phi_{crit}$ – are illustrated to the left. The implications of these parameters for the ecophysiological modeling framework employed in this study are illustrated for an example ecophysiology on the right, based on the median values of these distributions. The minimum seawater $pO_2$ required for resting metabolism increases at an exponential rate, primarily due to temperature effects on metabolic rate. The minimum seawater $pO_2$ required for sustainable populations on ecological timescales also increases exponentially with temperature (faster than the minimum $pO_2$ required for resting metabolism), as the product of an exponential function defined by $A_o$, $E_o$ and $\Phi_{crit}$ (see *Materials and Methods*, Equation 1 for full details). At a given temperature, the thermal safety margin of a hypothetical ecophysiology is predicted to decrease with decreasing seawater $pO_2$ (either due to a change in equilibrium oxygen levels resulting from decreasing atmospheric oxygen, or a decrease in saturation resulting from biogeochemical processes such as organic carbon remineralization). Surface seawater $pO_2$ distributions (far right) are taken from Earth system model analyses presented in this study (8 PAL CO$_2$ Ordovician configuration; points represent mean values, and error bars 2 SD).

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719
<figure><img src="image_7.png" /><figcaption>Figure 3: Ensemble Earth system modeling approach using cGENIE. A) Example of ensemble O₂-CO₂ Earth system modeling experiment using Ordovician continental configuration with modern remineralization depth. Color scale illustrates the proportion of the modeled ecophysiotypes that can live in each cell of the surface ocean (top 0-80m layer of the 16 layer cGENIE ocean is illustrated here). B) Sensitivity analyses performed in this study. Variations in continental configuration, remineralization depth and CO₂ forcing are detailed here. The atmospheric O₂ levels investigated in each ensemble experiment are the same. We present full surface ocean models here for ease of visual comparison, although main text analyses use only non-polar shelf environments (as illustrated in Fig. S3).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4: Metabolic habitat viability for the ensemble O₂-CO₂ Earth system modeling experiment with an Ordovician continental configuration and modern remineralization depth. Metabolic habitat viability is calculated as the proportion of the ecophysiotypes modeled that can live in shelf environments of each three dimensional ocean model produced using cGENIE (*Materials and Methods*). Circles represent individual experiments and are color coded by mean equatorial sea surface temperature (-3.2°S to 3.2°N, 0-80m depth), as a product of atmospheric CO₂ forcing.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 5: Global ecophysiotype extinction predicted to result from a ~5°C warming event from an approximately pre-industrial baseline sea surface temperature at different atmospheric O₂ levels, illustrating sensitivity analyses evaluating the influence of continental configuration on predicted extinction. Global ecophysiotype extinction is defined as the proportion of ecophysiotypes that were viable in the shelf environments of the colder baseline simulation but are no longer viable in any cell of the shelf environments of the warmer ocean simulation (*Materials and Methods*). Baseline CO₂ levels used: 1 PAL – Paleocene; 8 PAL – Permian; 8 PAL – Ordovician. Darker envelopes represent 25ᵗʰ–75ᵗʰ percentiles of 100 simulation results, lighter envelopes represent 5ᵗʰ–95ᵗʰ percentiles.</figcaption></figure>