Persistent global marine euxinia in the early Silurian

The second pulse of the Late Ordovician mass extinction occurred around the Hirnantian-Rhuddanian boundary (~444 Ma) and has been correlated with expanded marine anoxia lasting into the earliest Silurian. Characterization of the Hirnantian ocean anoxic event has focused on the onset of anoxia, with global reconstructions based on carbonate δ238U modeling. However, there have been limited attempts to quantify uncertainty in metal isotope mass balance approaches. Here, we probabilistically evaluate coupled metal isotopes and sedimentary archives to increase constraint. We present iron speciation, metal concentration, δ98Mo and δ238U measurements of Rhuddanian black shales from the Murzuq Basin, Libya. We evaluate these data (and published carbonate δ238U data) with a coupled stochastic mass balance model. Combined statistical analysis of metal isotopes and sedimentary sinks provides uncertainty-bounded constraints on the intensity of Hirnantian-Rhuddanian euxinia. This work extends the duration of anoxia to >3 Myrs – notably longer than well-studied Mesozoic ocean anoxic events. The Late Ordovician mass extinction has been attributed to extended marine anoxia. Here, the authors use a metal isotope mass balance model and find the marine anoxic event lasted over 3 million years, notably longer than the anoxic event associated with the Permian-Triassic extinction and Cretaceous ocean anoxic events.

T he early Silurian Rhuddanian Age directly follows the Hirnantian glaciation and Late Ordovician mass extinction 1,2 . Although biotic turnover may have begun in the Katian 3 , two extinction pulses are commonly considered to define the Late Ordovician mass extinction 2,4 : the first, around the Katian-Hirnantian boundary, is broadly attributed to global cooling and predominantly affected deep-water fauna 5 . The second, around the Hirnantian-Rhuddanian (Ordovician-Silurian) boundary, impacted animals throughout the water column 2 and has been linked to a sulfide-dominated Hirnantian ocean anoxic event [6][7][8] . The Rhuddanian Age was an interval of sustained low marine diversity 3 , with the fossil plankton (graptoloid) record indicating a dramatic increase in average extinction rates from the Late Ordovician 9,10 and benthic brachiopods exhibiting protracted post-extinction recovery 11,12 .
Geochemical data supporting the expansion of euxinic (anoxic, sulfide-rich) conditions are increasingly linked to extinction intensity around the Hirnantian-Rhuddanian boundary. A shift in the global marine redox landscape has been inferred from the extrapolation of regional iron and sulfur records 6,7 , and mass balance modeling of carbonate uranium isotope data 8 . It has also long been noted that there is an abundance of black shales in the early Silurian 13 and the Rhuddanian stands out as an especially productive interval in the distribution of organic-rich mudrocks throughout the Phanerozoic 14 . Although this lithological transition has been linked to rapid post-glacial warming and marine anoxia in the early Silurian 13,15,16 , the negative shift in carbonate uranium isotopes (indicating the expansion of bottom-water euxinia) precedes sequence stratigraphic and geochemical evidence for a major glacial pulse in the Late Ordovician 1,8,17 . Whether the low-diversity interval and extended organic-rich shale deposition through the Rhuddanian Age is linked to an extended Hirnantian-Rhuddanian ocean anoxic event, persisting through glacial extremes and Rhuddanian deglaciation, therefore remains an open question.
Existing quantitative constraints on Hirnantian-Rhuddanian global marine anoxia have been based on a commonly used twosink uranium isotope mass balance modeling framework with fixed flux and fractionation parameters 8 . There is increasing consensus that ferruginous (anoxic, sulfide-poor) conditions were common in the early Paleozoic 18 , although understanding of the implications of these depositional environments for trace metal isotope cycling remains limited. Separating non-euxinic and euxinic anoxic sinks in the δ 238 U mass balance significantly impacts model reconstructions of the global seafloor redox landscape 19 , as do varying fluxes and fractionations within reasonable limits based on observations from modern marine environments 20 . The combined effects of small deviations from commonly used model assumptions has the potential to dramatically alter implications for ancient marine environments. Therefore, moving forward, there is an obvious need for both improved statistical treatment of uncertainty and the use of multiple globally representative redox proxies and/or sedimentary archives in trace metal isotope mass balance models. With respect to the Hirnantian-Rhuddanian ocean anoxic event, both the complete duration and uncertainty-bounded estimates of the intensity and biogeochemical nature of the event are yet to be fully constrained.
In this study, we investigate the extent, nature, and duration of Hirnantian-Rhuddanian euxinia, and associated correlations with marine biodiversity and organic-rich shale deposition. Iron speciation measurements are used to establish the local water column conditions under which the Rhuddanian black shale succession from the Tanezzuft Shale Formation of the E1-NC174 core, Murzuq Basin, Libya, was deposited. Shale molybdenum and uranium concentrations are evaluated in the context of organic carbon correlations to infer the global applicability of measured trace metal isotope geochemistry. Molybdenum (δ 98 Mo) and uranium (δ 238 U) stable isotope data from the E1-NC174 core are then evaluated using a coupled Monte Carlo global mass balance model, in combination with previously generated carbonate δ 238 U data 8 . Our improved statistical treatment of uncertainties in metal mass balances, and incorporation of multiple trace metal isotopes and sedimentary archives, allows us to improve quantitative constraint on the spatial and temporal dynamics of Rhuddanian ocean euxinia (and associated uncertainties) during the Hirnantian-Rhuddanian ocean anoxic event. We confidently show that euxinic bottomwaters were likely two orders of magnitude more widespread than today for a period of >3 Myrs following the second pulse of the Late Ordovician mass extinction, significantly longer than well-studied Mesozoic ocean anoxic events.

Results
Geologic setting. The E1-NC174 core records approximately three million years of Rhuddanian black shale deposition through the Tanezzuft Shale Formation of the Murzuq Basin, Libya [21][22][23] . The Murzuq Basin is an intracratonic depression that is believed to have originally extended northward across the passive margin of Gondwana 24 ( Supplementary Fig. 1), recording Cambrian through Silurian siliciclastic marine sedimentation. The economic importance of the mid-Rhuddanian hot shale as a hydrocarbon source rock has led to the detailed biostratigraphic characterization of the E1-NC174 core 21,25,26 , indicating that the stratigraphic record presented here spans from at least the lowermost Rhuddanian (Normalograptis tilokensis Biozone in Gondwana, Akidograptus ascensus Biozone in global biostratigraphy) to the uppermost Rhuddanian (Neodiplograptus fezzanensis Biozone in Gondwana, Coronograptus cyphus Biozone in global biostratigraphy; Fig. 1; Supplementary Note 1). The mid-Rhuddanian hot shale varies in thickness across northern Gondwana, a trend that has been interpreted as the product of sea level change relative to post-glacial paleorelief 24 . Molybdenum and uranium enrichment factors covary in proportion with the modern seawater Mo-U ratio ( Supplementary Fig. 2), and the concentrations of both elements correlate linearly with organic carbon concentrations ( Supplementary Fig. 3). These trends suggest that, despite this inferred depositional relationship with paleobathymetry, changes in sea level did not markedly impact the resupply of molybdenum and uranium from the open ocean to the Murzuq Basin.
Molybdenum and uranium isotopes as tracers of marine redox. Molybdenum and uranium stable isotope ratios (δ 98 Mo and δ 238 U) have been developed as tracers of global marine redox conditions and widely applied to the geologic record [27][28][29] . There are distinct isotopic fractionation factors associated with the incorporation of both elements into oxic vs. euxinic sedimentary sinks, meaning that the extent of euxinic seafloor exerts strong controls on the δ 98 Mo and δ 238 U compositions of seawater 30,31 . Molybdenum and uranium have modern residence times of >400 kyr in seawater 30,32 and are therefore expected to be well-mixed on the timescales of ocean circulation, even during times of more widespread anoxia 33 . The heavy isotopic compositions of euxinic sedimentary sinks in both the molybdenum and uranium systems mean that seawater δ 98 Mo and δ 238 U values are predicted to decrease during periods of expanded euxinia. Euxinic shales track the seawater molybdenum and uranium isotopic composition, albeit with a predictable fractionation factor in most cases 30,31 . Increases in the areal extent of bottom-water anoxia (and, by extension, authigenic enrichments in marine sediments) further lead to drawdown of dissolved seawater molybdenum and uranium concentrations 34 , highlighting the value of mass balance models that fully couple seawater metal concentrations and isotope values in reconstructions of global redox (see Methods-Mo-U mass balance model). The redox-sensitivity of predictable authigenic enrichment differs between molybdenum and uranium 34 , specifically with respect to the high H 2 S concentrations required for quantitative molybdenum reduction 35 . We therefore propose that interpreting sedimentary records of both trace metal isotope systems in combination improves constraint on the biogeochemical nature of anoxic intervals, as does the combination of δ 238 U data from multiple sedimentary archives.
One of the appeals of the molybdenum and uranium proxies is that they can be easily translated into a quantitative framework to estimate the areal extent of euxinic conditions 8,27,29,36 . The global molybdenum and uranium mass balances are described by a series of ordinary differential equations relating the concentrations and isotopic compositions of molybdenum and uranium in seawater to relevant sources and sinks (Methods-Mo-U mass balance model). In contrast to the modern ocean, ferruginous (anoxic, sulfide-free) bottom-waters are thought to have been prevalent in the early Paleozoic 18 . Consequently, both trace metal cycles are modeled here with three redox-sensitive sinks: oxic, euxinic, and a broadly defined 'reducing' sink intended to capture the behavior of ferruginous, as well as traditionally suboxic, settings (following Reinhard et al. 36 ). The fluxes and isotopic fractionations associated with the trace metal sources and sinks in mass balance models are based on observations of trace metal behavior in modern environments [36][37][38] . We propose that the fixed parameterization employed in conventional mass balance approaches, although useful in quantifying general trends, inherently underrepresents uncertainty in the implications of trace metal cycling for global marine redox. To accommodate variations in trace metal behavior between modern depositional environments, and corresponding uncertainty about trace metal cycling in deep time, we employ Monte Carlo simulations, coupling the global molybdenum and uranium isotope mass balances. We randomly subsampled all model parameters that are known to vary across modern environments between minimum and maximum observed values (see Table 1 for parameter space and justification of parameter ranges). For all parameters, we assumed a uniform distribution, rather than forcing a most likely scenario with normal distributions. Broadly reducing, specifically ferruginous, environments are proposed to create some of the biggest uncertainties in this mass balance approach. We therefore widely varied the extent of broadly reducing environments for each logarithmically scaled scenario of ocean euxinia, between 0% of the global ocean floor and the remaining (non-euxinic) ocean floor available. Owing to the scarcity of modern ferruginous analogs 39 , we applied widely bounded fractionation factors for the broadly reducing sinks to incorporate plausible fractionation factors not observed in modern suboxic environments but expected based on laboratory experiments 40 . Supplementary  Fig. 6 illustrates that alternative parameterizations of these fractionation factors do not significantly alter the results of the fully coupled δ 238 U eux + δ 238 U carb + δ 98 Mo eux model.
The effects of expanding anoxic regions on local trace metal burial are also incorporated into our model using a pseudospatial scaling algorithm 36 . Estimates of trace metal burial in euxinic settings are conventionally based on observations from relatively shallow modern environments 37,38 . As widespread euxinia develops at a global oceanographic scale, authigenic enrichments in reducing environments occur at increasing depth as a physical consequence of expanded low-oxygen conditions. Trace metal complexes therefore travel greater distances through the water column, and are increasingly susceptible to remineralization 36 . The scaling factor used here links trace metal accumulation to organic   Fig. 1 Geochemistry of the E1-NC174 core. Core depth is presented alongside the global and regional Rhuddanian graptolite biozones based on Loydell 21,22 . a Total organic carbon concentrations (TOC) are presented from Loydell et al. 25 . b Iron speciation highly reactive to total iron ratios; samples with Fe HR / Fe T > 0.38 (dashed red line) are interpreted as anoxic. c Iron speciation pyrite to highly reactive iron ratios; samples with Fe Py /Fe HR > 0.7-0.8 are interpreted as euxinic (red box, lighter red=possibly euxinic). d Mo concentrations normalized to total organic carbon concentrations (ppm/wt%). e U concentrations normalized to total organic carbon concentrations (ppm/wt%). f δ 98 Mo illustrated relative to modern euxinic shale range 47 . g δ 238 U illustrated relative to observed modern euxinic/anoxic shale range 44 . In Mo and U concentration and stable isotope plots, white points represent bulk measurements, and gray points represent measurements corrected for detrital input (for molybdenum, the difference between these values is often negligible). Error bars associated with trace metal concentrations illustrate 2 SD uncertainty on crustal concentrations used in detrital corrections. Error bars associated with trace metal isotopes illustrate combined analytical error (2 SE) and 2 SD uncertainty on crustal concentrations used in detrital corrections. Computed error bars for δ 98 Mo (Supplementary Data 1) are generally smaller than the data points illustrated here.  Fig. 2). Mo/TOC and U/TOC ratios ( Fig. 1, Supplementary Fig. 3) exhibit no directional stratigraphic trend through the E1-NC174 core, indicating that increased enrichments in molybdenum and uranium through the hot shale interval ( Supplementary Fig. 4) are likely a product of organic carbon loading 41,42,46 . Similar trends are observed in vanadium enrichments, but not in metals such as chromium that are more sensitive to detrital input 34 ( Supplementary Fig. 4).
Modeling global marine redox conditions. Uranium and molybdenum isotope data exhibit low variance around isotopically light mean values (relative to modern euxinic shales) throughout the E1-NC174 core ( Fig. 1). Molybdenum isotope fractionations in low-sulfide conditions are expected to be highly variable 35,47 , whereas in high-sulfide conditions (with extensive thiomolybdate formation) fractionations are expected to exhibit low variance with near quantitative capture of seawater δ 98 Mo. We therefore make a case that low variances in δ 98 Mo (mean = 0.69‰, SD = 0.13‰) and δ 238 U (mean = −0.02‰, SD = 0.07‰) values are further indications that molybdenum and uranium were reduced with consistent fractionations (or, in the case of molybdenum, likely no fractionation) from seawater under the locally euxinic conditions overlying the Rhuddanian black shales of the Murzuq Basin. The absence of a directional stratigraphic response in δ 98 Mo and δ 238 U to organic carbon loading (and likely associated sea level change 25 ) is further evidence that the trace metal isotope compositions of the E1-NC174 core shales were not controlled by variation in local redox conditions and therefore reliably record global seawater signals. Furthermore, the offset (~0.5‰) between δ 238 U eux data presented here and published δ 238 U carb data 8 (reported mean value = −0.45‰) for the early Rhuddanian is consistent with modern observations of the fractionation between shallow-water carbonates and euxinic sedimentary sinks 44,48 ( Supplementary Fig. 5), adding support that the Murzuq Basin shales and Anticosti Island carbonates both record global trace metal isotope signals. Long-term variations in the marine molybdenum and uranium budgets are expected to be primarily controlled by the relative extent of strongly reducing sedimentary sinks due to high metal accumulation rates 30,31,37 . Isotopically light δ 98 Mo and δ 238 U values (relative to modern euxinic shales) are therefore broadly interpreted as recording globally expansive anoxic, likely euxinic, bottom-water conditions throughout the Rhuddanian. Molybdenum and uranium mass balance model results are illustrated as frequency distributions of predicted euxinic shale δ 98 Mo and δ 238 U values corresponding to a range of logarithmically scaled scenarios for global marine euxinia (Fig. 2, f eux describes the fraction of global seafloor with euxinic bottomwaters). The two-dimensional density plots presented here are applicable to all intervals of geologic time with modern-style tectonic and weathering processes (see Supplementary Fig. 6 for the impact of elevated weathering rates), assuming that the molybdenum and uranium cycles are in steady-state and the trace metals are well-mixed in seawater. Smoothed frequency distributions in Fig. 2a illustrate the distributions of modeled f eux values compatible with the range of measured euxinic shale values from the E1-NC174 core, and Rhuddanian shallow-water carbonates from Anticosti Island 8 . The corresponding f eux ranges illustrate the implications of δ 238 U eux , δ 238 U carb (from Bartlett et al. 8 ), δ 98 Mo eux , and δ 238 U eux + δ 238 U carb + δ 98 Mo eux values for the global extent of marine euxinia through the Rhuddanian Age. The combined δ 238 U eux + δ 238 U carb + δ 98 Mo eux distribution is proposed as the bestconstrained model of f eux for this time interval, indicating that euxinic bottom-waters were most likely two orders of magnitude more widespread than modern (median f eux estimate = 15.8%, mean = 29.2%), with a significant spread of feasible f eux values (5th percentile=2.0%, 95th percentile = 100%). Despite these uncertainties, fully coupled Rhuddanian model results are not compatible with modern levels of ocean oxygenation (Fig. 2), and indicate modally different, dominantly euxinic ocean biogeochemistry throughout the Rhuddanian.
The implications of the coupled Mo-U Monte Carlo mass balance model across the Ordovician-Silurian boundary are illustrated in Fig. 3a, based on Bartlett et al. 8 and the new geochemical data presented in Fig. 1 (see Methods-Stratigraphic age models). The best-fit distribution is defined by cross-validated LOESS models fitted to percentiles of the distribution of f eux estimates at each timestep. The time-dependent distribution is bounded by cross-validated LOESS models fitted to the 5th and 95th percentile f eux estimates. These are considered the most representative (and conservative) estimate of uncertainty across the logarithmic range of marine redox scenarios investigated here. The onset of globally widespread euxinia modeled here is well correlated with the second pulse of the Late Ordovician mass extinction 2,4,5 (Fig. 3), as indicated by the δ 238 U carb data of Bartlett et al. 8 . Our model constrains uncertainty on the spatial extent of euxinic seafloor, demonstrating that under all parameterizations of the uranium isotope cycle this geochemical data requires a 1-2 order-of-magnitude change in seafloor redox landscape. Furthermore, the incorporation of a second redoxsensitive δ 238 U archive and δ 98 Mo data allow us to increase quantitative constraint on the extent of Rhuddanian euxinia (relative to δ 238 U carb data alone, Fig. 2), and to extend the record of euxinia throughout the Rhuddanian Age.

Discussion
The black shale trace metal isotope data, and model evaluation of new and published 8 data, presented here, support previous inferences that the earliest Rhuddanian experienced widespread ocean anoxia [6][7][8] and extends the duration of the Hirnantian ocean anoxic event by >2 million years to include the entire Rhuddanian Age (Fig. 3). Our coupled approach, using euxinic shale molybdenum and uranium isotope measurements in conjunction with existing carbonate δ 238 U data from the Hirnantian and lower Rhuddanian 8 , allows us to dramatically improve quantitative constraints on the intensity of global marine euxinia. Namely, previously published Rhuddanian δ 238 U carb data 8 , when analyzed with a three sink redox model and full accounting of uncertainty in trace metal cycles, results in similar mean (17.9%) and median (6.3%) f eux results to our fully coupled model, but with a 5th percentile at 0.3% (Fig. 2). This broad distribution of compatible f eux estimates means that modeling of carbonate δ 238 U data alone cannot confidently reject f eux estimates compatible with, or at least similar in magnitude to, modern anoxic seafloor 36,37 for the Rhuddanian Age. A similar resulting f eux distribution is obtained for our δ 238 U eux data alone (Fig. 2), and these analyses echo other recent modeling results, demonstrating that the uranium isotope system in isolation has relatively limited resolution when all applicable redox sinks and sources of uncertainty are taken into account 19 . However, global sensitivity analyses demonstrate that for the combined δ 238 U eux + δ 238 U carb + δ 98 Mo eux model, even with considerable deviations from modern molybdenum and uranium cycling, it is very unlikely that the isotopic record through the Rhuddanian could have been generated without strongly reducing depositional environments likely two orders of magnitude more extensive than today (Fig. 2). The multiproxy approach applied here therefore indicates that early Silurian marine redox dynamics were modally different to the modern ocean, and adds statistically robust support to previous inferences of Rhuddanian euxinia 8 . Furthermore, these extensively euxinic conditions were relatively stable for at least the duration of the Rhuddanian Age, extending the duration of the Hirnantian-Rhuddanian OAE by approximately three times (to at least three million years). The Hirnantian-Rhuddanian OAE therefore lasted significantly longer than well-studied Mesozoic ocean anoxic events [49][50][51] , and was more similar in duration to other Paleozoic anoxic intervals in the Cambrian 33 and Devonian 52 .
The emerging picture of the Rhuddanian Age is one of poorly oxygenated oceans, with relatively low benthic marine biodiversity 3 and enhanced organic carbon burial 13,14 . Recent models of extinction timing based on global database analyses 3,4 (full curve shown in Fig. 3d, schematic representation of extinction pulses shown in Fig. 3b) add support to inferences that the second pulse of the Late Ordovician mass extinction coincided with the onset of the Hirnantian-Rhuddanian ocean anoxic event 8 (Fig. 3). The capture-recapture analysis of early Paleozoic marine biodiversity 3 illustrated in Fig. 3d is taken here as the most reliable reconstruction of both the extinction pulses and subsequent low-diversity interval owing to the explicit consideration of taphonomy and sampling bias, and the marked improvement in temporal resolution relative to previous reconstructions. High-resolution extinction rate analyses of brachiopod 5 and graptoloid 9,10 fossil occurrence data, although more vulnerable to sampling biases, also support a major extinction event that is well correlated with the onset of global euxinia (Fig. 3b). Post-extinction recovery appears to have been limited in the Rhuddanian, with the Hirnantian-Rhuddanian OAE as described in this study (Fig. 3a) constituting a relative low-point in Middle Ordovician through Silurian biodiversity 3 (Fig. 3d). Global biodiversity estimates suggest that marine biodiversity did not reach pre-extinction levels for more than five million years 3 (Fig. 3d). Specifically, Silurian brachiopod biodiversity shows a protracted post-extinction recovery, lasting until at least the Telychian 11,12,53 (starting~436 Ma), and high rates of planktonic (graptoloid) turnover persisted throughout the Silurian 9,10 (Fig. 3b). The correlation between the second pulse of the Late Ordovician extinction 4 and the prolonged low-diversity interval 3 with the onset and persistence of widespread global marine euxinia indicates a major potential role for global ocean deoxygenation (and/or its external drivers) in governing global marine biodiversity.
The Rhuddanian trace metal isotope record presented in this study quantitatively supports existing observations of extended euxinic shale deposition from sedimentological compilations [13][14][15] . This euxinic interval corresponds not only to the North African hot shale source rocks, but to established gas shales such as the Longmaxi Formation in China 54 and potential hydrocarbon targets such as the Bardo Formation in the Baltic Basin 55 .  The inferred expansion of global marine euxinia is thus correlated with the global appearance of economically important black shales in the early Silurian 14 , which our work indicates were probably deposited over large swaths of the global seafloor. Although regional basin dynamics likely control specific differences in the hydrocarbon potential of Rhuddanian black shales, the analyses presented here demonstrate the governing role that global-scale ocean biogeochemistry can play in the distribution of organic-rich mudrocks. The Hirnantian-Rhuddanian ocean anoxic event appears to be sustained across both peak Hirnantian glacial cooling 1,8 and postglacial Rhuddanian warming 8,13,15,16 . Therefore, although the onset of the anoxic event may coincide with a major pulse of the Hirnantian glaciation (cooccurring with the inception of Late Ordovician Glacial Cycle 3 8,17 , and predating the major late Hirnantian cooling episode indicated in oxygen isotope 8 and clumped isotope 1 data), the extended global record of Rhuddanian euxinia presented here demonstrates that the modally different ocean biogeochemistry we describe is non-specific to global climate state. The lack of correlation between model reconstructions of euxinia and both stratigraphic and geochemical reconstructions of climatic change indicate that this episode of ocean euxinia defies conventional, Cenozoic-Mesozoic models of temperature-driven marine euxinia 56 . The reducing postextinction oceans described here also likely persisted much longer than more frequently discussed Cenozoic and Mesozoic deep-time analogs for ocean deoxygenation. In other words, the Hirnantian-Rhuddanian ocean anoxic 'event' was an extended state change in marine redox, lasting a minimum of approximately six times longer than Cretaceous OAE2, and approximately two times the duration of early Triassic anoxia, based on current geochronological models and proxy data 23,50,51 . This adds to the existing examples of extended-duration Paleozoic 'OAEs' in the Cambrian 33 and Devonian 52 . Secular changes in atmospheric oxygen 57,58 and the biological pump 57,59 , as well as feedbacks between global climate, tectonics, and ocean circulation 16,60 , are all plausible driving mechanisms for this difference in the nature of ocean anoxic events through the Phanerozoic. Statistically robust models of global marine redox based on multiple geochemical proxies and/or sedimentary archives, as presented here, will provide valuable constraint for approaches that seek to address this problem through combined ocean circulation and biogeochemical modeling.
In summary, the Tanezzuft Shale Formation of the Murzuq Basin, Libya, hosts a continuous record of persistently euxinic black shale deposition, constraining the succession's exceptional suitability for the reconstruction of globally redox-sensitive trace metal isotope records. Combined Monte Carlo mass balance modeling of δ 238 U eux , δ 238 U carb , and δ 98 Mo eux measurements support reconstructions of a Hirnantian ocean anoxic event continuing into the earliest Silurian ascensus Zone, and further demonstrate that extended marine euxinia then lasted throughout the Rhuddanian Age (likely >3 million years). Euxinic seafloor area can be confidently reconstructed as at least an order of a Modeled distribution of f eux values across the Ordovician-Silurian boundary based on the geochemical data presented here and in Bartlett et al. 8 . Cross-validated LOESS models are fitted to percentile-bracketed estimates of global seafloor euxinia for each timestep (Methods-Mo-U mass balance model). The temporal extent of carbonate 8 and euxinic shale records analyzed in this model framework are shown above. b Per-capita brachiopod extinction rates 5 and per lineage million year (LMY) graptoloid extinction rates 9 with extinction pulses from capture-recapture analyses 3,4 schematically illustrated above. Envelope of graptoloid extinction rates illustrates mean values±68% quantiles from bootstrap means (1000 iterations) of median values for each 0.25 Myr time bin. We note that the relatively broad time bins of Kröger et al. 4 leave the timing and duration of the first extinction pulse illustrated loosely constrained. c Sea-surface temperature estimates based on clumped isotope paleothermometry 1 -with a line fitted (as in Finnegan et al. 1 ) to the lowest temperature estimates for each time bin. Error bars illustrate 1 SD analytical error reported in original study. d Genus-level diversity calculated by capture-recapture methods from the Paleobiology Database 3 , error envelope illustrates 95% confidence intervals for capture-recapture analyses. Vertical gray box illustrates time interval depicted in a-c, with the extended Cambrian-Silurian diversity trends illustrated as context for the end-Ordovician mass extinction and subsequent low-diversity interval. Age models are updated from original publications to reflect the Geologic Timescale 2012 23   ARTICLE magnitude, and more likely two orders of magnitude, more widespread in the early Silurian than the modern. Prolonged recovery from the Late Ordovician mass extinction supports the notion that ocean deoxygenation has repeatedly had prolonged impacts on marine ecosystems at the timescales of geological stages/ages. The early Silurian marine system appears to be defined by low benthic marine biodiversity, enhanced organic carbon burial, and dramatically expanded euxinic bottom-waters.

Methods
Mo-U mass balance model. The concentration and isotope values of molybdenum and uranium in seawater are described by Eqs. 1-4.
Fluxes for redox-sensitive sinks are defined as Where b i represents per area burial flux, A i is the seafloor area of the redoxsensitive sink, α i is a pseudospatial scaling coefficient relating burial rate to the area of the sink 36 , [Me] sw is the mean modeled concentration of the metal in seawater, and [Me] M.sw is the mean modern concentration of the metal in seawater.
The area of the redox-sensitive sink is defined as Where A is the total modern seafloor area, and f i is the proportion of the global seafloor that the redox-sensitive sink covers.
In the Monte Carlo simulations, f eux is modeled in 31 equally spaced, logarithmically scaled steps, with 1000 model runs per f eux scenario. The areal extent of the broadly reducing sink, f red , is included in the global sensitivity analysis, and allowed to vary between 0 and 100-f eux %. The range of f ox.lim , the maximum extent of oxic deposition, is defined in Table 1 based on Reinhard et al. 36 and physical oceanographic limits.
The burial scaling factor, α, introduces a pseudospatial burial coefficient to the model, such that burial in redox-sensitive sinks is scaled to the anticipated effects of organic carbon remineralization, depending on the extent of euxinic and broadly reducing conditions across the continental shelf. Where z = water depth, b oxic scales linearly such that: Burial rates in euxinic settings scale with the expansion of euxinia, following the remineralization model 61 used by Reinhard et al. 36 . In this scaling algorithm: Where α eux:min ¼ 1:58 À 0:16 ln min z eux ð Þ ð Þ and N(Z eux ) describes the number of depths modeled, and following Reinhard et al. 36 "we assume first that ∼5% of the shallow seafloor remains essentially authigenically neutral unless it becomes absolutely necessary to encroach upon this area (i.e., above 95% seafloor anoxia or euxinia)". min(Z eux ) and max(Z eux ) are defined by a crossvalidated LOESS model fitted to the global ocean depth-area relationships of Menard & Smith 62 based on the approach of Reinhard et al. 36 Burial rates in broadly reducing settings are also expected to scale with the extent of euxinic and broadly reducing conditions. Following previous application of Menard & Smith 62 to metal mass balances 36 , broadly reducing (ferruginous +   c References indicate that negative fractionations are possible under reducing conditions. As quantitative constraint on the magnitude of these fractionations is poor, the magnitude of the negative fractionation is parameterized to match the positive fractionation. See Supplementary Fig. 6 for alternative parameterizations. suboxic) settings are therefore modeled as expanding below the euxinic depositional environments along continental shelves in the pseudospatial scaling algorithm. The effect of α red = 1 (no pseudospatial scaling in broadly reducing settings) is illustrated in Supplementary Fig. 6.
The purpose of this mass balance approach is to generate frequency distributions of feasible synthetic sedimentary Mo-U isotope values for a range of paleoredox scenarios. This is achieved by calculating euxinic shale and carbonate values as described in Eqs. 10-12, with fractionation factors in local settings allowed to vary independently of the global means.
Mass balance equations are solved using a variable coefficient ordinary differential equation solver to fully couple marine trace metal concentrations and isotope values (rather than assuming instantaneous steady-state), and 1000 time-dependent Monte Carlo simulations are run to steady-state for 31 logarithmically scaled scenarios of ocean euxinia. Monte Carlo simulations are performed by global sensitivity analyses using the R packages deSolve 63 and FME 64 to subsample uniform prior distributions between prescribed minimum and maximum values. Minimum and maximum values for Monte Carlo variables are listed, with associated references, in Table 1. Fixed model parameters are listed with associated references in Table 2.
Time-dependent f eux reconstructions in Fig. 3 are generated by fitting crossvalidated LOESS models to percentiles of the f eux distribution modeled for each timestep. In this case, cross-validated LOESS models are fitted to the δ 238 U carb data from Bartlett et al. 8 and the measured δ 98 Mo eux and δ 238 U eux data presented here, based upon their associated age models (see Stratigraphic age models). At 20 kyr timesteps between 444.7 and 440.8 Ma, 'measured' geochemical values are generated from these local regression models for each metal isotope proxy archive that (based on stratigraphic age models-see below) has recorded data for that time interval. The next step in this time-dependent modeling exercise is to describe modeled distributions of f eux for each geological timestep. In order to preserve the resolution of distributions presented in Fig. 2, we binned both modeled and measured sedimentary geochemical data before filtering the modeled data set based on measured sedimentary values. Measured geochemical data for each proxy were subdivided into 10 proxy-specific bins (spaced uniformly between the minimum and maximum measured values). Modeled data were then binned at the same resolution for each of the three metal isotope proxy archives investigated. The binned f eux distribution from global sensitivity analyses was then filtered to only include results compatible with measured geochemical data at each timestep in order to calculate the distributions described in Fig. 3a. Cross-validated LOESS models are used throughout this modeling workflow to reduce modeling subjectivity.
Sample digestion and geochemical analyses. Shale samples were obtained from drill core and crushed to fine powder in a tungsten carbide shatterbox. Iron speciation measurements were performed at Stanford University following the methods of Poulton & Canfield 65 for sequential extraction and Canfield et al. 66 for chromium reduction of sulfur (CRS). Pyrite sulfur isotope values were analyzed on silver sulfide from the CRS extraction, measured on an Isoprime 100 isotope ratio mass spectrometer interfaced with an Elementar vario ISOTOPE cube elemental analyzer at Virginia Tech. Total organic carbon concentrations and carbon isotope values were first presented in Loydell et al. 25 , and were measured on a Finnegan Mat 251 mass spectrometer coupled to a Fisons 1108 elemental analyzer. Major and trace metal concentrations were measured at Bureau Veritas/ACME Labs using standard four acid digestion techniques and Q-ICP-MS/ICP-OES.
Uranium and molybdenum isotopes were measured at the Yale Metal Geochemistry Center using previously described acid digestion, column chromatography, and MC-ICP-MS methods 44,67,68 : All powdered samples were first ashed at 600°C to remove organic matter, then dissolved at 105°C in a concentrated HF-HNO 3 mixture followed by an aqua regia digestion of the dried residues. Final residues were dissolved in 5 ml 6 N HCl as stock solutions. Splits from these solutions were then used for Mo and U column chromatography and subsequent analyses via MC-ICP-MS. Bulk trace metal concentrations and isotope measurements were corrected for detrital input with associated error estimates for authigenic values 69 (Methods-Detrital corrections).
Uranium isotope analyses. Uranium was separated using UTEVA ion exchange resin and the 233 U-236 U double spike method (see methods of Weyer et al. 44 and Wang et al. 67 ). Double spike was added to samples prior to column chemistry to achieve a 236 U/ 238 U ratio of~30 based on previously measured U concentrations. Samples were then evaporated to dryness, brought up in 3 N HNO 3, and then purified on UTEVA resin columns using the method of Wang et al. 67 with 100 ng loaded on the columns. After purification, U samples were dissolved in 0.3 N HNO 3 to achieve 50 ppb solutions for analysis. Samples were measured on a Thermo-Finnigan Neptune Plus multi-collector ICP-MS. Samples were introduced using an Apex IR sample introduction system and measured at low resolution. Isotopes 232 Th, 233 U, 235 U, 236 U, and 238 U were measured simultaneously on Faraday collectors connected to 10 11 Ω amplifiers achieving a signal of~40-45 volts on 238 U.
The CRM 112a (New Brunswick Laboratory, US Department of Energy) standard was analyzed every three samples to monitor drift in instrumental mass bias, and samples were normalized to the average of bracketing standards. The drift of the standard was typically <0.08‰. Blank levels were always <0.1% of sample voltage, and no blank correction was made (the correction would be smaller than the instrumental counting error). External reproducibility was 0.1‰ based on duplicate sample preparations and analyses of both the USGS geostandard NOD-A-1 (0.60‰, 2 SD = 0.04‰, n = 4) and sample duplicates (n = 5). Values for U isotope measurements and associated error can be found in Supplementary Data 1.
Molybdenum isotope analyses. Molybdenum isotope measurements were also performed at the Yale Metal Geochemistry Center following methods in Planavsky et al. 68 . An aliquot of the bulk digests was doped with a 97 Mo-100 Mo double spike, prepared gravimetrically from Oak Ridge Laboratory metal powders 68 . The spiked sample was dried to ensure sample-spike equilibration. We used a two-stage column separation 68 . We first ran the sample through an anion resin (AG-MP-1M) and then a cation resin (AG50W-X8) to remove any remaining Fe. The Mo isotopic ratios were analyzed using a Neptune Plus MC-ICP-MS. Molybdenum isotope compositions are reported using the δ notation, where δ Detrital corrections. Detrital corrections are made for bulk molybdenum and uranium isotope and concentration measurements to infer authigenic isotope and concentration values. Authigenic molybdenum estimates are generally very similar to bulk measured values, as expected based on understanding of detrital trace metal shuttling 45 . However, both metals are standardized for consistency. Uranium values are standardized to thorium based on the conservative behavior demonstrated in Cole et al. 69 ; molybdenum values are standardized to aluminum based on the methods of Noordmann et al. 71 , using crustal average data from Rudnick & Gao 72 . Standardizing uranium to aluminum results in lighter δ 238 U than standardizing to thorium (corresponding to more widespread euxinia; Supplementary Data 1), so thorium standardization is presented as the most conservative estimate of authigenic δ 238 U.  Stratigraphic age models. Extrapolating the time-independent models of Fig. 2 to the time-dependent biogeochemical change modeled in Fig. 3a is dependent on stratigraphic age models. For the E1-NC174 core (data measured in this study), we use a continuous stratigraphic age model beginning at the Hirnantian-Rhuddanian boundary (443.8 Ma) and ending at the Rhuddanian-Aeronian boundary (440.8 Ma). Although the core does not definitively contain these boundaries based on biostratigraphy, the core does contain all graptolite biozones of the Rhuddanian 21 and this is therefore proposed as the most appropriate age model. For the published δ 238 U carb data, we use a continuous stratigraphic age model based on the description in Bartlett et al. 8 and references therein (including Wickson 74 ). The age models used in both Finnegan et al. studies 1,5 presented in Fig. 3 are updated to reflect the Geologic Timescale 2012 23 and illustrate stratigraphic relationships relative to δ 238 U carb data 8 where applicable. Notes justifying the placement of placement of stratigraphic boundaries in these stratigraphic age models are included in Supplementary Data 2.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data presented in this manuscript are included in Supplementary Data File 1 (all new geochemical data presented in Fig. 1 and Supplementary Figs. [2][3][4][5] and Supplementary Data File 2 (data from published studies used in Fig. 3 with revised age models and associated justifications).

Code availability
The model code required for the analyses presented in this manuscript can be freely accessed on Github [https://github.com/richardstockey/Mo.U.Model].