Richard G Stockey
email: rstockey@stanford.edu

Alexandre Pohl

Andy Ridgwell

Seth Finnegan

Erik A Sperling

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

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

Captions

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.

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 nearmodern 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 ecophysiotype 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.

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

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 [START_REF] Raup | Mass extinctions in the marine fossil record[END_REF][START_REF] Foote | Origination and extinction components of taxonomic diversity: general problems[END_REF][START_REF] Foote | Origination and extinction through the Phanerozoic: A new approach[END_REF][START_REF] Alroy | Dynamics of origination and extinction in the marine fossil record[END_REF][START_REF] Alroy | A more precise speciation and extinction rate estimator[END_REF][START_REF] Kröger | The evolutionary dynamics of the early Palaeozoic marine biodiversity accumulation[END_REF][START_REF] Kocsis | The r package divDyn for quantifying diversity dynamics using fossil sampling data[END_REF]. First revealed by analyses of stratigraphic range compilations ~40 years ago [START_REF] Raup | Mass extinctions in the marine fossil record[END_REF], this decline has been a persistent feature in subsequent analyses using expanded databases [START_REF] Alroy | Effects of sampling standardization on estimates of Phanerozoic marine diversification[END_REF] 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. [START_REF] Kocsis | The r package divDyn for quantifying diversity dynamics using fossil sampling data[END_REF][START_REF] Finnegan | The Red Queen revisited: reevaluating the age selectivity of Phanerozoic marine genus extinctions[END_REF][START_REF] Payne | The effect of geographic range on extinction risk during background and mass extinction[END_REF]). 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 [START_REF] Finnegan | The Red Queen revisited: reevaluating the age selectivity of Phanerozoic marine genus extinctions[END_REF][START_REF] Payne | The effect of geographic range on extinction risk during background and mass extinction[END_REF][START_REF] Peters | Determinants of extinction in the fossil record[END_REF][START_REF] Peters | Environmental determinants of extinction selectivity in the fossil record[END_REF][START_REF] Payne | The evolution of complex life and the stabilization of the Earth system[END_REF]. 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) [START_REF] Sperling | Statistical analysis of iron geochemical data suggests limited late Proterozoic oxygenation[END_REF][START_REF] Dahl | Devonian rise in atmospheric oxygen correlated to the radiations of terrestrial plants and large predatory fish[END_REF][START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF][START_REF] Wallace | Oxygenation history of the Neoproterozoic to early Phanerozoic and the rise of land plants[END_REF][START_REF] Stolper | A record of deep-ocean dissolved O2 from the oxidation state of iron in submarine basalts[END_REF][START_REF] Sperling | A long-term record of early to mid-Paleozoic marine redox change[END_REF] and the most recent generation of long-term carbon cycle models support only limited early Paleozoic oxygenation [START_REF] Krause | Stepwise oxygenation of the Paleozoic atmosphere[END_REF][START_REF] Lenton | COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time[END_REF][START_REF] Tostevin | Reconciling proxy records and models of Earth's oxygenation during the Neoproterozoic and Palaeozoic[END_REF]. 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 [START_REF] Sperling | Biodiversity response to natural gradients of multiple stressors on continental margins[END_REF][START_REF] Sampaio | Impacts of hypoxic events surpass those of future ocean warming and acidification[END_REF][START_REF] Levin | Oxygen Minimum Zone Benthos: Adaptation and Community Response to Hypoxia[END_REF]. 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 [START_REF] Pörtner | Physiology and climate change[END_REF][START_REF] Pörtner | Oxygen-And capacity-limitation of thermal tolerance: A matrix for integrating climate-related stressor effects in marine ecosystems[END_REF][START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Howard | Climate-driven aerobic habitat loss in the California Current System[END_REF][START_REF] Duncan | Different drivers, common mechanism; the distribution of a reef fish is restricted by local-scale oxygen and temperature constraints on aerobic metabolism[END_REF][START_REF] Boag | Metabolic tradeoffs control biodiversity gradients through geological time[END_REF]. Following the principle of oxygen and capacity limited thermal tolerance (OCLTTsee Pörtner [START_REF] Pörtner | Oxygen-And capacity-limitation of thermal tolerance: A matrix for integrating climate-related stressor effects in marine ecosystems[END_REF] for review), the critical partial pressure of oxygen (pO2) 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; [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF]). 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 [START_REF] Reddin | Marine clade sensitivities to climate change conform across timescales[END_REF].

Seawater pO2 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; [START_REF] Boag | Oxygen, temperature and the deep-marine stenothermal cradle of Ediacaran evolution[END_REF]). 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 [START_REF] Reinhard | Earth's oxygen cycle and the evolution of animal life[END_REF][START_REF] Sperling | The Ecological Physiology of Earth's Second Oxygen Revolution[END_REF], as is observed at regional scales in modern marine environments [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Howard | Climate-driven aerobic habitat loss in the California Current System[END_REF][START_REF] Duncan | Different drivers, common mechanism; the distribution of a reef fish is restricted by local-scale oxygen and temperature constraints on aerobic metabolism[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF]. 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 O2 and CO2 concentrations (Fig. 3a) in order to address uncertainty associated with reconstructions of these variables through geologic time (e.g.

Fig 1b

). Atmospheric CO2 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 [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF]. 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 [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Reddin | Marine clade sensitivities to climate change conform across timescales[END_REF][START_REF] Boag | Oxygen, temperature and the deep-marine stenothermal cradle of Ediacaran evolution[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF], 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 [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF][START_REF] Reddin | Climate change and the latitudinal selectivity of ancient marine extinctions[END_REF][START_REF] Saupe | Extinction intensity during Ordovician and Cenozoic glaciations explained by cooling and palaeogeography[END_REF]. 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 [START_REF] Meyer | 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[END_REF] and has been previously identified as important in differentiating earlier Paleozoic vs. Mesozoic-Cenozoic ocean redox conditions [START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF]. 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 [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF].

Utilizing model ocean temperatures and dissolved oxygen concentrations from this suite of 3D realizations, we then employ the recently developed Metabolic Index [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF] 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 O2 concentrations to link our ecophysiological models more directly to the trends we observe in the geologic record. We model 'global ecophysiotype 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 CO2, and are designed to mimic relatively extreme (but not infrequent) examples of climate variability known to have occurred throughout Earth history [START_REF] Gutjahr | Very large release of mostly volcanic carbon during the Palaeocene-Eocene Thermal Maximum[END_REF][START_REF] Foster | Future climate forcing potentially without precedent in the last 420 million years[END_REF]. We performed this extinction sensitivity estimate for all scenarios of atmospheric O2 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 O2 levels in a framework that is neutral to the selective advantage of different ecophysiotypes and explores the relative impacts of other aspects of Phanerozoic Earth system evolution on marine animal ecophysiology.

Results

Metabolic habitat viability

Metabolic habitat viability -the total number of ecophysiotypes that can live in the shelf environments of each ancient ocean realization -is shown in Figure 4 for our ensemble of simulations employing an Ordovician continental configuration and Figure S4 for all other Earth system model ensembles. In all ensembles, metabolic habitat viability exhibits a positive relationship with atmospheric O2 and a negative relationship with global temperature, resulting in a sigmoidal relationship between atmospheric O2 and metabolic habitat viability for a given climate state. While our primary analyses focus on shelf environments, very similar trends are observed when the entire ocean is considered as potential habitat for marine animals (Figs. S5-6). Using entire 3D ocean models often results in slightly higher estimates of metabolic habitat viability because cooler and (at higher atmospheric O2 levels) well-oxygenated deep ocean environments may act as refugia, despite being unlikely diversity hotspots for other reasons such as food supply [START_REF] Woolley | Deep-sea diversity patterns are shaped by energy availability[END_REF][START_REF] Sperling | The temporal and environmental context of early animal evolution: Considering all the ingredients of an "explosion[END_REF]. The latitudinal distribution of viable ecophysiotypes is similar between continental configurations (Figs. 3,S3), with more ecophysiotypes generally able to inhabit cooler polar regions of the surface ocean. Identifying oceanographic and tectonic drivers of regional differences in viable ecophysiotype biogeography between continental configurations is an interesting avenue for future research, but this would require tuning model ensembles to precisely equivalent climate states and is therefore beyond the scope of this study.

Changes in metabolic habitat viability between cooler and warmer climate states allow us to begin linking ecophysiology to extinction sensitivity under different scenarios of atmospheric oxygenation (Fig. 4). The percentage change in metabolic habitat viability between colder and warmer ocean simulations increases as atmospheric oxygen decreases from modern levels (100% preindustrial atmospheric levels, PAL), reaching maximum values at atmospheric O2 concentrations reconstructed for the early Paleozoic (<40-50% PAL -grey envelope in Figs. 4,S4). Furthermore, because extinction is expected to track the proportional rather than absolute loss of ecophysiotypes [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF][START_REF] Reddin | Climate change and the latitudinal selectivity of ancient marine extinctions[END_REF], differences in metabolic habitat viability would therefore be amplified to result in high levels of warming-driven extinction at lower levels of atmospheric oxygenation reconstructed for the early Paleozoic (e.g. 20% PAL). In such cases, the number of viable ecophysiotypes in any low-O2 scenario is much lower, leading to greater proportional change despite a similar absolute change in metabolic habitat viability to that predicted at higher O2 levels (e.g. 80% PAL). Our analyses thus demonstrate that the expected proportional loss of ecophysiotypes during a global warming event would be substantially higher under the limited surface oxygenation scenarios hypothesized for the early Paleozoic than for modern (or, more broadly, post-mid-Paleozoic) levels of oxygenation.

Sensitivity of global extinction simulations

When we simulate the impacts of warming events on the extinction of marine animals, the general trend of greater predicted extinction at lower atmospheric oxygen levels is a common feature of all ensemble experiments (Figs. 5, S7-12). There is variation in the absolute magnitude of simulated global ecophysiotype extinction between the different continental configurations tested in our ensembles. However, we do not observe any time-dependent trends in extinction sensitivity between progressive continental reconstructions and the most striking feature of the results remains the nonlinear increase in global ecophysiotype extinction with declining oxygenation (Figs. 5, S7-12). Below 10% PAL O2, the relationship between modeled extinction and atmospheric oxygen is less predictable. However, the initial pool of viable pre-extinction ecophysiotypes in these simulations is very small at low O2 levels (Figs. 3,4, S4-6), suggesting that this inconsistency may partially be a product of sample size (as supported by wider envelopes in Figures 5, S7-8, S11) and the limited experimental respirometry data for animals that live in low-oxygen environments [START_REF] Seibel | Oxygen supply capacity in animals evolves to meet maximum demand at the current oxygen partial pressure regardless of size or temperature[END_REF]. We find that modeled extinction is high at low atmospheric oxygen levels for the Ordovician configuration regardless of the efficiency of the biological pump (Fig. S7), suggesting that despite its role in controlling the distribution of oxygen in the upper water column this factor is of secondary importance compared to atmospheric O2. A shallower remineralization depth does slightly reduce simulated extinction sensitivity at low atmospheric O2 (~5-20% PAL), although the same caveat of small initial ecophysiotype population size applies.

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

Altering the assumptions of our ecophysiological modeling approach does not dramatically impact the observed relationship between predicted extinction and atmospheric oxygen. In our primary analyses, we define ecophysiotype extinction as only occurring when no cells in our simulated shelf environments are viable for an ecophysiotype, or (equivalently) when 100% of metabolically viable habitat area is lost.

Increasing the minimum number of inhabited ocean cells required to survive a warming event or decreasing the proportional habitat loss required for extinction moderately increases modeled global ecophysiotype extinction (Figs. S9-10). The differences in the 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 O2, 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 (Ao -the inverse of the hypoxic threshold; and Eothe temperature-dependency of the hypoxic threshold). While varying the distributions of Ao and Eo does impact the magnitude of simulated global ecophysiotype extinction and the steepness of the inflection around 40% PAL O2 (Fig. S11), the overarching trends of increasing predicted extinction with decreasing atmospheric oxygen (down to at least 20% PAL O2) 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).

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 O2 (Fig. 5). This broadly coincides with the upper limit of estimates for early Paleozoic atmospheric oxygen (Fig. 1; [START_REF] Sperling | Statistical analysis of iron geochemical data suggests limited late Proterozoic oxygenation[END_REF][START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF]), 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 ecophysiotype 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 O2. 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 [START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF]. 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 O2 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. 2; [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF]), 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 early Paleozoic sea surface temperatures have been reconstructed using geochemical paleothermometry [START_REF] Hearing | An early Cambrian greenhouse climate[END_REF][START_REF] Trotter | Did Cooling Oceans Trigger Ordovician Biodiversification? Evidence from Conodont Thermometry[END_REF][START_REF] Goldberg | A high-resolution record of early Paleozoic climate[END_REF], suggesting that this phenomenon could also provide some explanatory power for the frequency of high magnitude extinctions. However, while recent advances in paleothermometry have shown considerable potential in constraining Ordovician climate dynamics [START_REF] Goldberg | A high-resolution record of early Paleozoic climate[END_REF], there is still currently limited consensus on directional trends in early Phanerozoic sea surface temperatures, especially in the earliest Paleozoic [START_REF] Vérard | On plate tectonics and ocean temperatures[END_REF][START_REF] Mills | Modelling the long-term carbon cycle, atmospheric CO2, and Earth surface temperature from late Neoproterozoic to present day[END_REF][START_REF] Bergmann | A paired apatite and calcite clumped isotope thermometry approach to estimating Cambro-Ordovician seawater temperatures and isotopic composition[END_REF]. In contrast, reconstructions of atmospheric oxygen are increasingly consistent, even if the absolute concentrations are poorly constrained (Fig. 1). This discussion should not be framed as a question of oxygen versus temperature; inherent to the concept of temperature-dependent hypoxia is that both interact to determine aerobic safety margins (Fig. 2; [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF][START_REF] Reddin | Marine clade sensitivities to climate change conform across timescales[END_REF][START_REF] Deutsch | Metabolic trait diversity shapes marine biogeography[END_REF]). We therefore suggest that atmospheric oxygen can confidently be established as a first-order control on early Phanerozoic extinction rates, and that potentially high sea surface temperatures through the same time interval may have amplified the effect of limited surface oxygenation. Changes in the latitudinal focus of fossil sampling through the Phanerozoic may also bias extinction rate estimates slightly towards environments with warmer baseline temperatures in the early Paleozoic when well-sampled paleocontinents were mainly situated at low latitudes [START_REF] Close | The spatial structure of Phanerozoic marine animal diversity[END_REF]. However, the lack of similarly high long-term extinction rates across other intervals with low-latitude sampling biases indicates that this cannot fully explain the trends we seek to explain here and that atmospheric oxygen (potentially compounded by warm global climate) is a more probable first-order driver.

The plausibility of temperature-dependent hypoxia as a driver of high early Paleozoic extinction rates is supported by links between sea-level changes, ocean deoxygenation and extinctions during this interval of geologic time. The Cambrian-Ordovician biomere-style extinction events are a major component of the characteristically high early Paleozoic extinction rates [START_REF] Bambach | Origination, extinction, and mass depletions of marine diversity[END_REF]. These events, predominantly observed in the trilobite fossil record (e.g. ( 55)), have long been recognized as associated with sedimentological evidence of sea-level change, possible thermocline migration and ocean anoxia [START_REF] Palmer | The biomere problem: evolution of an idea[END_REF][START_REF] Fortey | There are Extinctions and Extinctions: Examples from the Lower Palaeozoic[END_REF]. 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 [START_REF] Saltzman | Persistent oceanic anoxia and elevated extinction rates separate the Cambrian and Ordovician radiations[END_REF][START_REF] Gill | Geochemical evidence for widespread euxinia in the Later Cambrian ocean[END_REF][START_REF] Elrick | Oxygen-isotope trends and seawater temperature changes across the Late Cambrian Steptoean positive carbon-isotope excursion (SPICE event)[END_REF][START_REF] Leroy | Variable redox conditions as an evolutionary driver? A multi-basin comparison of redox in the middle and later Cambrian oceans (Drumian-Paibian)[END_REF]. Similar inferences of ocean anoxia and rapid sea-level change have been made for the early-middle Cambrian Botomian extinction event [START_REF] Zhuravlev | Anoxia as the cause of the mid-Early Cambrian (Botomian) extinction event[END_REF]. 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 [START_REF] Gutjahr | Very large release of mostly volcanic carbon during the Palaeocene-Eocene Thermal Maximum[END_REF][START_REF] Yao | Large-scale ocean deoxygenation during the Paleocene-Eocene thermal maximum[END_REF][START_REF] Jenkyns | Geochemistry of oceanic anoxic events[END_REF][START_REF] Clarkson | Uranium isotope evidence for two episodes of deoxygenation during Oceanic Anoxic Event 2[END_REF] 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. [START_REF] Alroy | Dynamics of origination and extinction in the marine fossil record[END_REF][START_REF] Erwin | Lessons from the past: Biotic recoveries from mass extinctions[END_REF]) 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 animalenvironment interactions from the geologic record [START_REF] Knoll | Systems paleobiology[END_REF]. However, the power of Earth system models in simulating the ecophysiological responses of marine animals to deeptime environmental change is only recently being explored (e.g. [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF][START_REF] Saupe | Extinction intensity during Ordovician and Cenozoic glaciations explained by cooling and palaeogeography[END_REF]). 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 [START_REF] Knoll | Paleophysiology and end-Permian mass extinction[END_REF][START_REF] Payne | End-permian mass extinction in the oceans: An ancient analog for the twenty-first century?[END_REF]. Meta-analytical approaches also suggest that sensitivity to climate related stressors correlates with genus survival across ancient extinction events in higher-level taxonomic groups [START_REF] Reddin | Marine clade sensitivities to climate change conform across timescales[END_REF], 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 animalenvironment 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).

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.

Materials and Methods

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. [START_REF] Kocsis | The r package divDyn for quantifying diversity dynamics using fossil sampling data[END_REF] 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 airbreathing 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 percapita). 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. [START_REF] Boag | Metabolic tradeoffs control biodiversity gradients through geological time[END_REF]) and remain within the >0.4 quorum range suggested by previous authors [START_REF] Kocsis | The r package divDyn for quantifying diversity dynamics using fossil sampling data[END_REF]. 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. [START_REF] Kocsis | The r package divDyn for quantifying diversity dynamics using fossil sampling data[END_REF] 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

Earth system modeling cGENIE is an Earth system model of intermediate complexity, designed for the spatiallyexplicit biogeochemical characterization of global-scale paleoceanographic problems [START_REF] Ridgwell | Marine geochemical data assimilation in an efficient Earth system model of global biogeochemical cycling[END_REF]. We create a range of 3D realizations of ocean environmental conditions by running ensembles of model experiments. The experiments in each ensemble are run for 10,000 or 20,000 years (depending on the time required for dissolved marine O2 and seawater temperature to reach steady state), with a range of assumptions regarding atmospheric pCO2 specific to each continental configuration based upon changing solar luminosity and albedo (Fig. 3), and a range of atmospheric pO2 scenarios based upon the range of estimates from geochemical proxy data and biogeochemical box modeling (Fig. 1). We 75) -see Supplementary Information). For each of these ancient configurations, we apply age-appropriate solar luminosity [START_REF] Gough | Solar Interior Structure and Luminosity Variations[END_REF]. Finally, in the Ordovician continental configuration, we explore the possible influence of postulated Phanerozoic changes in organic matter remineralization depth in the ocean [START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF][START_REF] Meyer | 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[END_REF] testing ~1/3 modern, or specifically 200m, e-folding depth [START_REF] Meyer | 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[END_REF].

We extract the mean annual ocean temperature and ocean oxygen concentration during the last simulated year for each ocean grid cell and generate corresponding seawater oxygen partial pressure (pO2) values (following the methods of Hofmann et al. [START_REF] Hofmann | Hypoxia by degrees: Establishing definitions for a changing ocean[END_REF]). Changing solar luminosity and continental configuration results in different mean sea surface temperatures for the same atmospheric CO2 assumption. We therefore identify the pCO2 values for each continental configuration that create comparable climate simulations (Table S2 and Supplementary Information). These climate scenarios are compared using the mean equatorial (3.2°S to 3.2°N, 0-80m depth) surface ocean temperature for each model (Figs. 4, S4-6, Table S2).

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

Ecophysiological simulations of animal-environment interactions

We combine our ensemble Earth system modeling experiments with a probabilistic ecophysiological model of extinction vulnerability. The Metabolic Index, as defined by Deutsch et al. [START_REF] Deutsch | Climate change tightens a metabolic constraint on marine habitats[END_REF] and Penn et al. [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF], 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.

Metabolic habitat viability

= / A ! 𝑝O " exp 5 -E ! k # 9 1 T - 1 T $%& <= > ϕ '$() *+, (/ ! , 1 ! , 2 "#$% ) *(4 (/ ! , 1 ! , 2 "#$% ) (Eqn. 1)
At the individual organism scale, Ao is the inverse of the hypoxic threshold of the organism (the minimum required seawater pO2 (Pcrit) to sustain resting aerobic metabolism), Eo is the temperature-dependency of the hypoxic threshold and ϕ '$() 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. [START_REF] Penn | Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction[END_REF]). Seawater temperature and pO2 are environmental variables generated through the Earth system modeling approach described above (Fig. 3). All other variables are constant (kB is the Boltzmann constant, Tref is a reference temperature of 15°C). In our analyses, we simulate 1000 ecophysiotypes, sampled from probability distributions of Ao, Eo and ϕ '$() , and map which cells (if any) in each cGENIE ocean model are habitable for each ecophysiotype (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 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 [START_REF] Krause | Stepwise oxygenation of the Paleozoic atmosphere[END_REF][START_REF] Lenton | COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time[END_REF][START_REF] Edwards | Oxygenation as a driver of the Great Ordovician Biodiversification Event[END_REF] and geochemical proxy records [START_REF] Sperling | Statistical analysis of iron geochemical data suggests limited late Proterozoic oxygenation[END_REF][START_REF] Lu | Late inception of a resiliently oxygenated upper ocean[END_REF]. 1). The frequency distributions of three key ecophysiological parameters -Ao, Eo and ϕ !"#$ -are illustrated to the left. The implications of these parameters for the ecophysiological modeling framework employed in this study are illustrated for an example ecophysiotype on the right, based on the median values of these distributions. The minimum seawater pO2 required for resting metabolism increases at an exponential rate, primarily due to temperature effects on metabolic rate. The minimum seawater pO2 required for sustainable populations on ecological timescales also increases exponentially with temperature (faster than the minimum pO2 required for resting metabolism), as the product of an exponential function defined by Ao, Eo and ϕ !"#$ (see Materials and Methods, Equation 1 for full details). At a given temperature, the thermal safety margin of a hypothetical ecophysiotype is predicted to decrease with decreasing seawater pO2 (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 pO2 distributions (far right) are taken from Earth system model analyses presented in this study (8 PAL CO2 Ordovician configuration; points represent mean values, and error bars 2 SD). Figure 3: Ensemble Earth system modeling approach using cGENIE. A) Example of ensemble O2-CO2 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 CO2 forcing are detailed here. The atmospheric O2 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).

  apply this ensemble treatment to three ancient continental configurations, chosen for their relatively even temporal spread across the Phanerozoic and history of use in previous paleoclimate studies -a Paleocene configuration (~55 Ma) (43), a Permian configuration (~251 Ma) (74), and an Ordovician configuration (~450 Ma) (based on the 8 PAL CO2 FOAM simulation of Pohl et al. (

  10.5281/zenodo.4473048. Configuration files for the specific experiments presented in the paper can be found in the directory: genie-userconfigs/MS/stockeyetal.PNAS.2021. Details on the experiments, plus the command line needed to run each one, are given in the readme.txt file in that directory. All other configuration files and boundary conditions are provided as part of the code release. A manual detailing code installation, basic model configuration, tutorials covering various aspects of model configuration and experimental design, plus results output and processing, is assigned a DOI: 10.5281/zenodo.4469678.

( 33 )Figures 692 Figure 1 :

 336921 Figures 692

Figure 2 :

 2 Photosynthetic fractionation COPSE GEOCARBSULFOR

29 Figure 4 : 9 Figure 5 :

 29495 Figure4: Metabolic habitat viability for the ensemble O2-CO2 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 CO2 forcing.

  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 Ao and Eo defined in Penn et al.

	We further simulate extinction as the global loss of ecophysiotypes between two
	CO2 scenarios at the same O2 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 O2 level. Global ecophysiotype
	extinction is thus
	<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.

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

has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 838373. cGENIE simulations were performed on Earth system modelling cluster facility at University of California Riverside. Ecophysiological analyses of cGENIE models were performed on the Sherlock cluster at Stanford University. We thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to this research.