Sulfidic anoxia in the oceans during the Late Ordovician mass extinctions – insights from molybdenum and uranium isotopic global redox proxies

The Late Ordovician Mass Extinction wiped out 85% of animal species in two phases (LOME1 and LOME2). The kill mechanisms for the extinction phases are debated, but deteriorating climate and the expansion of marine anoxia appear to have been important factors. Nevertheless, the spatial extent and intensity of marine anoxia and its temporal relationship with the extinctions are not well understood. Here, we review existing global paleoredox proxy data based on molybdenum (Mo) and uranium (U) isotopes from four paleocontinents combined with new Mo isotope data from Dob's Linn, Scotland. Individually, these sedimentary records demonstrate significant redox fluctuations, but our coupled dynamic oceanic mass balance model for the evolution of the marine Mo and U cycles reveals that globally expansive ocean anoxia is best constrained by δ 238 U in carbonates from Anticosti Island that record expansive anoxia during LOME2. In addition, we consider periodic sulfidic anoxia developing in well-ventilated parts of the shallow oceans (e.g. during warmer periods with greater solar insolation) to have produced temporarily high seawater δ 98 Mo values during LOME1 in accordance with trends to high values observed in the sedimentary records. In this view, oceanic oxygen loss had a causal role during both extinction phases in the Late Ordovician.


Introduction
The Late Ordovician Mass Extinction (LOME) ~444 million years ago (Ma) was one of the largest extinctions in animal history, wiping out about a quarter of all animal families, more than half of the marine invertebrate genera, and ~ 85% of all animal species (Harper et al., 2014). The classic view of the event denotes two extinction phases with the first interval (LOME1) occurring near the Katian-Hirnantian Global Stage boundary and the second interval (LOME2) in the later Hirnantian and start of the Metabolograptus persculptus Graptolite Biozone. Both of these intervals have been inferred to coincide with dramatic environmental changes recorded by sedimentary evidence for extensive glacial deposits in high-latitude Gondwana with associated faunal response across a number of clades (e.g., Colmenar et al., 2018;Ghienne et al., 2014), a sudden drop in seawater surface temperatures (Finnegan et al., 2011), positive carbon isotope excursions in marine carbonates and sedimentary organic carbon near the start of the Hirnantian (Brenchley et al., 2003), as well as associated eustatic sea level falls in the order of 100 m (Brenchley et al., 2006;Haq and Schutter, 2008;Nielsen, 2004). All of the above lines of evidence point to advancing icehouse conditions during the early Hirnantian first extinction episode (Finnegan et al., 2011). That is, climatic cooling is inferred to drive glacial advances in Gondwana and falling sea level globally. The second extinction phase is linked to global warming and the expansion of marine anoxia, deglaciation and continental flooding 1 (Brenchley, 1984;Brenchley et al., 1994;Hammarlund et al., 2012;Harper et al., 2014). There has been increased focus on geochemical evidence for sulfidic water column anoxia in local basins (euxinia) as well as for intensified volcanism both during LOME1 and LOME2 (Bond and Grasby, 2020;Hammarlund et al., 2012;Zou et al., 2018). In the current study we place global constraints on expansive sulfidic marine anoxia (euxinia) in the context of global cooling. Still, the spatial extent and intensity of euxinia and its temporal evolution across the first extinction episode is not well understood.
Here, we aim to reconstruct the globally integrated scale of ocean oxygenation across both phases of the Late Ordovician Mass Extinction using molybdenum and uranium isotope records from shales and limestones in South China, Libya, Scotland and Canada. The comparison of global redox tracers from sedimentary sequences worldwide provides an overlapping set of (although imperfect) oxygenation records permitting an evaluation of the relationship between the expansion of anoxic or euxinic 'dead zones' in the oceans and the two phases of the LOME.

Theorytracking ocean oxygenation with Mo and U isotopes
Molecular oxygen (O 2 ) is a reactive gas not reliably preserved in the geological record, but its presence and abundance in the past can be indirectly traced. This tracing is possible since atmospheric and oceanic O 2 levels control oxidation and reduction (redox) reactions of for example iron, molybdenum, and uranium at the Earth's surface. Molybdenum (Mo) and uranium (U) isotopes in ancient sediments have become powerful redox-sensitive proxies that can track global marine anoxia including the proportion of sediment buried in anoxic versus oxic settings in the oceans. Both systems build on the same two fundamental attributes: 1) seawater has a homogeneous Mo and U isotope composition, and 2) the heavier isotopes ( 98 Mo and 238 U) accumulate in the oceans when a larger proportion of marine sediment burial occurs in oxygenated settings. Extensive reviews of the modern marine Mo and U isotope cycles are found elsewhere (Andersen et al., 2017;Kendall et al., 2017;Zhang et al., 2020b), and we summarize the important aspects regarding the Mo and U removal pathways from seawater to sediments as necessary for this study.

The Mo isotope proxy
Nearly all molybdenum in the ocean is sourced via rivers and ground water (90-99.6%) with a small contribution from hydrothermal systems (Kendall et al., 2017;McManus et al., 2002;Miller et al., 2011). Today, the marine molybdenum input is thought to be balanced by sediment burial that occurs through three major sinks: 1) fully oxygenated sedimentary settings with Mn-oxides present (OX, ~35% of total burial); 2) oxic sediments with sulfide accumulation in the sediments (SAD, ~55% of total burial); and 3) euxinic sediments deposited under an anoxic and sulfidic water column (EUX, ~10% of total burial) (Kendall et al., 2017). The Mo residence time in the modern oceans is ~440 kyr today, which greatly exceeds the average ocean water mixing time of 1.5 kyr (Miller et al., 2011). Therefore, the modern oceans have a uniform Mo concentration of 105 nM and a uniform Mo isotope composition, δ 98 Mo = 2.3‰, where the δ notation describes the deviation in ‰ relative to a universal scale defined by NIST-SRM-3136 Nägler et al., 2013;Nakagawa et al., 2012). Mass balance models suggest that the oceans could have been homogeneous with respect to the Mo isotope composition since at least the Neoproterozoic (Dahl et al., 2011;Reinhard et al., 2013).
Due to the preferential incorporation of lighter Mo isotopes in the three sedimentary sinks (more on this below), oceans have an approximately 2‰ heavier δ 98 Mo than oceanic input with an average δ 98 Mo of 0.2 to 0.8‰, and also higher than average crustal rocks at ~0‰ (Archer and Vance, 2008;Horan et al., 2020;King et al., 2016;King and Pett-Ridge, 2018;Revels et al., 2021;Voegelin et al., 2014;Willbold and Elliott, 2017). Traditionally, the magnitude of Mo isotope fractionation expressed in sediments is thought to increase in more oxygenated settings (e.g., Poulson Brucker et al., 2009), so that seawater will become enriched in the heavier Mo isotopes at times in Earth history when a larger proportion of marine sediments were buried under oxygenated waters. However, a relatively large Mo isotope fractionation is also expressed in mildly and intermittently euxinic settings (Brüske et al., 2020;Neubert et al., 2008;Noordmann et al., 2015). Although these settings are probably not a major marine Mo sink today, we explore the possibility that intermittently and/or mildly sulfidic bottom water conditions prevailed over a larger portion of the open continental shelves in the Late Ordovician.
We assess how intermittent euxinia on the global scale would affect the Mo isotope composition of seawater by introducing the major marine Mo sources and sinks, and by paying special attention to the isotope fractionation processes associated with sediment burial. In oxygenated seawater and river water, molybdenum exists as a hexavalent oxyanion molybdate, MoO 4 2− that is soluble and behaves conservatively with only relatively slow adsorption onto Fe-and Mn-oxyhydroxides. A large equilibrium isotope fractionation (Δ MnOx-SW = − 2.9‰) is expressed for Mo adsorbed onto Mn-oxides relative to ambient seawater solution, and this offset is expressed in Mn-oxide rich sediments (Barling et al., 2001;Barling and Anbar, 2004;Siebert et al., 2003;Wasylenki et al., 2008). The burial of Mn-oxides is notoriously slow, but can be observed in ferromanganese nodules in the deep abyssal oceans, where there is little or no terrigenous sedimentation (e.g. Calvert and Piper, 1984). In contrast to the slow deposition of Mo in deep-sea oxic settings, Mo burial rates are high in anoxic and sulfidic basins where the aqueous hydrogen sulfide (H 2 S) concentrations exceed a threshold of 11 μM in the bottom waters and molybdate is converted to thiomolybdates (Ericksson and Helz, 2000). The isotope fractionation preserved depends on the completeness of the conversion, and this is controlled by both time and access to sulfide. In short, the conversion of MoO 4 2− to MoO 4-x S x 2− occurs in a sequence of sulfidation reactions (in four steps), where each subsequent step produces a chemical form with increased affinity for reaction with particulate matter (Ericksson and Helz, 2000;Vorlicek et al., 2004;Vorlicek and Helz, 2002). Each sulfidation reaction, x = 0 → 1 → 2 → 3 → 4, involves a substantial equilibrium Mo isotope fractionation of ~1‰ between reactant and product (Kerl et al., 2017;Tossell, 2005). The kinetics of each sulfidation reaction is ~10-times slower than the previous step (Ericksson and Helz, 2000), and the increased particle reactivity for the more sulfidized species likely arises because thiomolybdates with two sulfur ligands (MoO 2 S 2 2− , MoOS 3 2− , and MoS 4 2− ) can be reduced to Mo(IV)-polysulfides by S 0 or other reductants (Dahl et al., 2013a;Dahl et al., 2010a;Helz et al., 1996;Vorlicek et al., 2004). Presumably, Mo reduction lowers the coulombic resistance to forming cuboidal structures on pyrite (Vorlicek et al., 2004). Thus, one hypothesized mechanism for Mo scavenging is that Mo(IV)-polysulfides react with particulate sulfide minerals and/or in association with organic matter Dahl et al., 2010a;Vorlicek et al., 2004). Alternatively, thiomolybdates can be removed by conversion to a Fe(II)-Mo(VI)-sulfide mineral, and Mo reduction occurs at a later stage (Helz and Vorlicek, 2019). In any case, the sulfidation reactions are rate-limiting for Mo scavenging and, isotope fractionation will be expressed if the reactive Mo species are being scavenged before the sulfidation reactions to the most reactive species are complete. The completeness of Mo removal is controlled mainly by three environmental factors: total H 2 S availability, pH, and the time available for sulfidation in the zone where thiomolybdate formation occurs (Dahl et al., 2010a;Helz and Vorlicek, 2019;Nägler et al., 2011). In addition, Mn-oxides settling across the chemocline can also deliver isotopically fractionated molybdenum into the euxinic zone. The oxides will dissolve and release isotopically light Mo that also undergoes sulfidation, contributing further to a negative δ 98 Mo offset in the sediments relative to the open ocean composition (Noordmann et al., 2015;Reitz et al., 2007;Scholz et al., 2013).
Euxinic sediments deposited near the chemocline, where O 2 and H 2 S meet, typically record a lower δ 98 Mo value than sediments in (usually) deeper, more sulfidic parts of the basin (Dahl et al., 2010a;Neubert et al., 2008;Noordmann et al., 2015;Scholz et al., 2013). Empirically, modern euxinic sediments display a wide range of negative δ 98 Mo offsets from seawater with the largest offset expressed in more open basins with lower H 2 S availability and faster water renewal rate into the sulfidic zone. The smallest offsets, of 0-0.8‰, are expressed in highly euxinic restricted settings (Kyllaren Fjord in Norway, deep part of the Black Sea; modern Lake Cadagno in Switzerland; and the seawater-fed Lake Rocogniza in Croatia), intermediate isotope fractionations of − 1.0 to − 1.6‰, are observed in more open euxinic settings such as Saanich Inlet in Canada and the Cariaco Basin in Venezuela. Large isotope fractionations of up to 3‰ are found near the chemocline of the Black Sea (Amini et al., 2016;Arnold, 2004;Arnold et al., 2012;Brüske et al., 2020;Bura-Nakić et al., 2018;Dahl et al., 2010a;Neubert et al., 2008;Noordmann et al., 2015;Scholz et al., 2018). Today, the modern euxinic sink comprises predominantly sediments deposited in silled basins (e.g. the Black Sea, Baltic sea) that, on average, are associated with a small isotope fractionation relative to overlying oxic seawater (− 0.5‰). If anoxic waters with low H 2 S concentrations, even intermittently, invaded more open-settings on the continental shelves (e.g. modern Peru and Namibian upwelling systems; (Schunck et al., 2013;Weeks et al., 2004), we would predict less quantitative conversion of Mo to particle-reactive phases and, perhaps also more Mn-oxide shuttling across the chemocline, therefore a larger magnitude of Mo isotope fractionation (− 3.0 to − 1.0‰) associated with the intermittently euxinic Mo sink.
In oxygenated settings where sulfide is confined to the sediment pore fluids, Mo is delivered to the sediments both via diffusion in aqueous form and with particulate Fe-and Mn-oxyhydroxides, with subsequent release as the oxyhydroxides dissolve. The Mo is subsequently sequestered into sulfide minerals (Zheng et al., 2000). In such settings, Mo removal from the sulfidic pore fluids can be more complete than from sulfidic seawater in intermittently euxinic settings, resulting in a smaller isotope offset in the sediments relative to contemporaneous seawater likely because diffusive escape out of the sulfidic zone is less efficient than is dissolved Mo escape from an intermittently sulfidic water column (Clark and Johnson, 2008). Thus, sediments with sulfide at depth (SAD) have intermediate δ 98 Mo values offset on average by − 0.7‰ (− 0.3 to − 1.8‰) relative to overlying seawater (Kendall et al., 2017;Poulson Brucker et al., 2009;Siebert et al., 2003).

The U isotope proxy
The marine U cycle is also sensitive to the global extent of ocean oxygenation. Like Mo, uranium is delivered into the oceans in its most oxidized form via rivers and ground waters with little or no isotope fractionation relative to continental crustal rocks (Noordmann et al., 2016;Tissot and Dauphas, 2015). Uranium has a long residence time (~400 kyr) compared to ocean mixing time and, thus, a homogeneous composition in the modern oceans (Dunk et al., 2002;Tissot and Dauphas, 2015). The marine removal pathways for U are sediment burial into: 1) euxinic sediments (9 ± 6% of total removal) and the rest in other sediments. The other sediments include 1) reducing sediments under oxygenated waters (40 ± 10% of total removal); 2) biogenic carbonates (30 ± 10% of total removal); 3) oxic sediments (12 ± 6% of total removal); and 4) hydrothermal alteration of crust at high temperature (3 ± 3% of total removal) and low temperature (6 ± 6% of total removal) (Lu et al., 2020).
In oxygenated seawater, uranium exists in its hexavalent oxidation state in various uranyl carbonate complexes (UO 2 (CO 3 ) 3 4− , CaUO 2 (CO 3 ) 3 3− , MgUO 2 (CO 3 ) 2 3− , and Ca 2 UO 2 (CO 3 ) 3 0 ) (Chen et al., 2016;Dong and Brooks, 2006). Uranium can be incorporated into aragonite with seawater U/Ca ratio (~3.1 ppm / 40 wt%) and essentially no U isotope fractionation (Livermore et al., 2020), presumably because the entire uranyl triscarbonate unit, UO 2 (CO 3 ) 3 4− , is passively incorporated into the crystal structure (Reeder et al., 2000). In calcite, U is hosted in a distinct molecular coordination environment, and U is 'excluded' from the mineral lattice with lower U/Ca than seawater (typically, U < 0.5 ppm / 40 wt%). Abiotic calcite precipitation is, however, also not associated with U isotope fractionation according to controlled laboratory experiments (Chen et al., 2016), and biogenic primary calcites also carry the same isotope composition as the seawater in which they grew (Livermore et al., 2020). Although primary marine carbonates seem to record seawater composition, marine carbonate sediments with diagenetic carbonate that precipitated from the pore fluids contain U with distinctly higher 238 U/ 235 U ratio than contemporaneous seawater (Zhang et al., 2020b and references therein). The U isotopic offset in marine carbonate sediments varies around an average value of 0.27 ± 0.14‰ (1SD; Chen et al., 2018;Romaniello et al., 2013;Tissot et al., 2018). In anoxic settings, U is preferentially removed from seawater to sediments with large U isotope fractionation expressed in the sediments as would be the case during both abiotic and biotic reduction of U(VI) to U(IV) (Brüske et al., 2020;Bura-Nakić et al., 2018;Holmden et al., 2015;Noordmann et al., 2015;Rolison et al., 2017;Stirling et al., 2007;Weyer et al., 2008). Reductive uranium removal in anoxic settings predominantly occurs below the sediment water interface. For example, the majority of U in the anoxic parts of the Black Sea is retained as U(VI) and little U is removed on particulate matter sinking through the water column (Anderson et al., 1989). Most U isotope fractionation occurs during reduction in fully anoxic settings with euxinic sediment burial producing a positive offset from overlying seawater of Δ U,EUX = 0.60 ± 0.20‰ (Andersen et al., 2017;Brüske et al., 2020;Zhang et al., 2020b). Other marine U sinks include reducing settings with anoxic sediments deposited under oxic waters (~26% of oceanic U input), deltaic sediments in the coastal zone (~19% of oceanic U input), basaltic alteration (10% of oceanic U input) and ferromanganese crust (~2% of oceanic U input) (Dahl et al., 2014;Dunk et al., 2002;Tissot and Dauphas, 2015). All of these U sinks are associated with smaller isotope fractionations than the anoxic sink. For the purpose of tracking the anoxic sink, these sinks are grouped as one sink with an average isotope offset relative (Δ U, OTHER ≈ 0.05 ± 0.10‰; discussed further below). As a result, seawater δ 238 U is particularly sensitive to the anoxic fraction of global U burial.

Predicting seawater δ 98 Mo from seawater δ 238 U
Sediment burial in euxinic sediments is a significant sink for both marine Mo and U that, in turn, largely determines the Mo and U isotope composition of the oceans. On this basis, Dahl et al. (2019) proposed a general first-order relationship between the euxinic fraction of global marine Mo burial (f Mo,eux ) and the euxinic fraction of global marine U burial (f U,eux ): The three parameters a = 0, b = 1, λ = 1.24± 0.38 are constants determined from two extreme end member criteria: 1) if f Mo,eux = 0, then f U,eux = 0, and 2) if f Mo,eux = 1, then f U,eux = 1, and a third constraint from the euxinic burial fractions in the modern ocean (f U,eux , f Mo,eux ) = (5 ± 3%, 9 ± 6%) (Dahl et al., 2019;Lu et al., 2020). With this relationship between euxinic Mo and U burial, seawater δ 98 Mo can be predicted from the seawater δ 238 U record if the following eight parameters are constrained: The oxic proportion of non-euxinic Mo sinks, denoted γ = f Mo,OX / (1 -f Mo,EUX ); the isotope compositions of oceanic input (δ 98 Mo IN, δ 238 U IN ); and the average isotopic offsets from seawater in each sink (Δ Mo,OX , Δ Mo,SAD, Δ Mo,EUX , Δ U,EUX , Δ U,OTHER ). At steady state, the Mo isotope composition of seawater is determined by the balance between oceanic input and outputs via burial in three distinct redox settings (Kendall et al., 2017): Similarly, the seawater U isotope composition is at steady state simplified to (Andersen et al., 2017): Previous comparisons of U and Mo isotope records suggest the euxinic sediment burial fraction is better constrained by U isotopes in marine carbonates than from the Mo isotope record in marine shales (Dahl et al., 2019). The present study evaluates this assertion for the Late Ordovician records. To illustrate large scatter in the ~445-443 Ma Mo isotope records of euxinic shales, we compare the steady state seawater δ 98 Mo composition predicted from the δ 238 U record to the observed δ 98 Mo records from Dob's Linn and other euxinic sedimentary sequences that span the LOME interval.

Main geological locality and sampling
We report new molybdenum isotopic and geochemical data from a deep-water sequence of graptolitic shales from the Upper Hartfell and Birkhill Shales (55 • 25 ′′ 46.799 N, 3 • 16 ′′ 16.977 W) at Dob's Linn, Southern Uplands, Scotland ( Fig. 1), which is the Global Stratotype Section and Point (GSSP) for the base of the Silurian (~443.8 ± 1.5 Ma) (Cohen et al., 2013). The shales were deposited in a tropical, continental margin setting facing the Iapetus Ocean (Scotese and McKerrow, 1990;Fig. 2).
A total of 37 hand-size samples were collected from outcrop across ~18 m of stratigraphy spanning the upper Katian (part), the Hirnantian and some lower Rhuddanian (lower Silurian) strata. The two main extinction pulses LOME1 and LOME2 both occur within ~1.5 m of stratigraphic section, near the Katian-Hirnantian boundary and ~ 1 m below the period boundary at the GSSP, respectively. The Dob's Linn section is rather condensed with only 3.5 m of Hirnantian deposits. The Hirnantian spans 1.4 Myr according to the Geological Timescale 2020, but see Ling et al., 2019 for a shorter duration (0.47 ± 0.34 My) corresponding to a sedimentation rate of 7.4 mm/kyr (4-27 mm/kyr) at Dob's Linn.

Analytical methods
The bulk rock elemental abundances, including Mo, U Al, Fe, and the Mo isotopes were analyzed on a Thermo X-Series ICPMS and Neptune MC-ICPMS at Arizona State University, respectively.
Elemental abundances were measured in the samples along with an international standard SDO-1 (Devonian Ohio Shale) and blanks to verify instrument accuracy and monitor for any contamination in the laboratory. The Al-normalized enrichment factors (EF) of Mo and U were calculated relative to the Upper Continental Crust (UCC) (EF = [metal/Al] sample / [metal/Al] UCC ). The UCC values for Al and Mo are 8.04 wt% and 1.5 μg/g, respectively (Taylor and McLennan, 1995).
Iron speciation data, including the pyrite and other highly reactive iron phases (Fe HR = Fe-carbonates + Fe-oxides + magnetite + pyrite) were extracted using calibrated protocols at NordCEE, University of Southern Denmark (Poulson Brucker et al., 2009). The long-term reproducibility (2SD) of standard material PACS-2 for the highly reactive to total iron ratio, Fe HR /Fe T , and the ratio of pyrite to highly reactive iron (Fe P /Fe HR ) are 0.44 ± 0.06, and 0.40 ± 0.09‰, respectively.
The molybdenum isotope composition of bulk samples was analyzed in 30 shales of which 14 samples were previously reported by Dahl et al., 2010b. Additional samples were collected to fill gaps in the stratigraphic record with isotope data carried out as described in our previous study. Briefly, approximately 10-50 g of shale was crushed and powdered, and a small fraction (~100 mg) of rock powder was processed. First, the samples were ashed at >700 • C and completely digested using HF and aqua regia before Mo was separated from matrix elements using a 2-column ion chromatographic purification scheme. Mo blanks were below <1% and international reference materials (SDO-1, Nod-A1, Nod-P1, Seawater) were analyzed in the laboratory for comparison of data with other laboratories. The Mo isotope composition was measured as 98 Mo/ 95 Mo deviation in per mil from the inhouse RochMo2 Mo standard solution (Johnson Matthey Company (Alfa Aesar), Specpure, ICP Mo standard solution, lot 802309E) with δ 98/95 Mo value of − 0.09 ± 0.05‰ relative to NIST 3134 defined as 0.25‰. We report all Mo isotope data on the international scale, adopting the abbreviation δ 98 Mo = δ 98/ 95 Mo NIST3134 . Two aliquots of SDO-1 were processed with the samples and gave δ 98 Mo = 1.01 ± 0.06 (n = 2) and δ 98 Mo = 1.15 ± 0.15 (n = 2) identical to the recommended value of 1.05 ± 0.14‰ (2SD) (Goldberg Terrigenous material carries minerals with small amounts of lithogenic Mo to sediments. This lithogenic Mo is unreactive in the marine Mo cycle, but included in our isotope measurements, and therefore, this needs to be excluded from the assessment of seawater composition. To correct for this detrital contamination, we use Al as a tracer for the detrital material and calculate the isotope composition of authigenic Mo, denoted δ 98 Mo auth , relative to the UCC as follows: The detrital δ 98 Mo composition is assumed to be 0.3‰ representing average upper crust (Kendall et al., 2017;Voegelin et al., 2014).

Local redox conditions at Dob's Linn, Scotland
A summary of geochemical data from the GSSP is shown in Fig. 3.
The local redox conditions on the continental margin of Avalonia at Dob's Linn are constrained from iron speciation and trace metal enrichments in the studied shales (Table 1), and the inferred redox conditions are shown in the figure as black (euxinic) or white (non-euxinic) bars. Apart from a brief interval in the upper Katian Dicellograptus T. typicus -D. mirus Graptolite Biozone (locally known as anceps zones), most samples display Fe HR /Fe T < 0.38 indicative of deposition in what were potentially oxic bottom waters. These samples with Fe HR /Fe T < 0.38 are also gray shales characterized by low Mo enrichments (Mo EF < 2) consistent with sediment deposition from O 2 -bearing bottom waters that contrast with the black Birkhill Shale above, suggesting that these sediments likely deposited under oxic conditions. Anoxic conditions prevailed from the middle of the M. persculptus Biozone and into the Silurian ( Fig. 3; Table 1). Half of the samples where the ratio of Fe HR /Fe T indicate water column anoxia have Fe P /Fe HR ratios between 0.62 and 0.74, which falls just at and below the threshold used   (Brenchley et al., 2003). The local redox conditions are inferred from the Fe speciation, Mo and U data (euxinic = black, non-euxinic = white). The red curves are smoothing spline fits through the data.  to characterize deposition in an euxinic setting (Poulton and Canfield, 2011). Strictly, these samples are better defined as ferruginous, but may represent periods of (mild) euxinic conditions. We note, however, that in these outcrop samples, the weathering of pyrite might have increased the pool of Fe-oxides, leading to an underestimation of the Fe P /Fe HR ratios and a bias towards ferruginous anoxia. Each sample represents ~2-4 cm of stratigraphic thickness, so these signatures represent the prevailing time-averaged redox conditions over ~1-10 kyr (albeit the signal is biased towards periods with higher sedimentation rate). Redox fluctuations might well have occurred over shorter time scales, however, this is beyond the resolution of our data.

Global redox constraints from δ 98 Mo at Dob's Linn, Scotland
The δ 98 Mo auth curves from Dob's Linn show no systematic trends, but rather large point-to-point variations from − 0.38‰ to 1.43‰ (− 0.72‰ to 1.42‰) around a median value of 0.65‰ that is close to the value of average oceanic input today, and 0.6-0.8‰ below that of modern open marine euxinic settings (Brüske et al., 2020;Fig. 6). The highest δ 98 Mo auth value is 1.4‰. The magnitude of the observed Mo isotopic variations at Dob's Linn (1.7‰) compares in magnitude to that of weakly euxinic settings today, but part of the observed variation may also represent temporal variation in the composition of Late Ordovician seawater (Kendall et al., 2017).

Discussion
The Dob's Linn section records sedimentation at the margin of the southern edge of Laurentia, where the local redox conditions, as inferred from iron speciation, molybdenum enrichments and facies, changed from mainly oxic in the mid-Hirnantian to anoxic ferruginous and euxinic conditions from the upper Hirnantian (middle of the M. persculptus Biozone) and into the Rhuddanian.
The local redox variability involving ferruginous and euxinic conditions may explain part of the large δ 98 Mo variations (1.7‰), but it is difficult to distinguish variable isotope fractionation in the local basin from the global changes in the open ocean. Further, fluctuating bottom water redox conditions in the Hirnantian are also observed at Nanbazi and Wangjiawan, South China, where high δ 98 Mo values occur in association with 100-fold Mo enrichments in the sediments at this time, suggestive of water column euxinia at ~100-150 m water depths on the Yangtze shelf (Zhou et al., 2015). Therefore, we consider periodic and mild euxinia in shallower shelf settings as plausible disturbances during LOME1, and compare our new data with global redox proxies from other basins in an attempt to reconstruct the molybdenum and uranium isotope record of the oceans in detail.

Anticosti Island, Canada
The Laframboise Section at Anticosti Island, Canada, shows high δ 238 U values (− 0.1 ± 0.1‰) in the middle Hirnantian carbonates, with values similar to modern marine carbonates, and a dramatic transition to lower δ 238 U values (− 0.6‰, down to − 0.8‰) beginning in the M. persculptus Biozone and extending into the Silurian (Bartlett et al., 2018) (Fig. 4). This negative shift (− 0.5‰) in the sedimentary record has been previously ascribed to a shift in seawater δ 238 U composition associated with an expansion of the global anoxic U sink (Bartlett et al., 2018).
We use this δ 238 U record to calculate theoretical steady-state seawater δ 98 Mo curves following eqs. 1-3, allowing us to compare this record to Mo isotope records from (euxinic) shales at Dob's Linn and other localities.
We start by estimating the δ 238 U composition of Late Ordovician seawater from the marine carbonate record at Laframboise Point by assuming a constant offset between seawater and sediments similar to today, Δ SW-Sed = 0.27 ± 0.14‰ (1 SD) Romaniello et al., 2013;Tissot et al., 2018). The magnitude of this offset might have been larger if the carbonates were deposited under anoxic waters, but the sedimentology and paleobiology indicate moderately oxygenated conditions. Importantly, the Laframboise Point section records a large negative δ 238 U shift (0.5-0.7‰) greater than the maximum local variation observed in modern marine carbonate settings. This trend is not associated with diagenetic alteration and is also decoupled from systematic Ca isotope variations higher in the section that have been suggested diagenetic in origin and the cause of a positive δ 13 C excursion (Fig. 4, Jones et al., 2020). Local redox variations could at most account for approximately half the magnitude of the negative δ 238 U trend, and since the sedimentology and paleobiology in the Laframboise section suggest no major shifts in local redox conditions, we follow the original interpretation and also assume that the dramatic negative shift in the M. persculptus Biozone mirrors a change in open ocean δ 238 U composition with a constant isotopic offset.
Next, we use eq. 3 to calculate the euxinic proportion of global U burial, f U,EUX , from the seawater δ 238 U curve assuming the marine U cycle is close to steady state at all times. From this, we predict theoretical curves for f Mo,EUX (eq. 1) and calculate a theoretical seawater curve (δ 98 Mo δ238U ) using eq. 2. The latter also depends on the relative proportion of sediment burial in oxygenated setting with manganese oxides preserved, γ = f OX /(1-f EUX ). Fig. 5 shows three curves (γ = 0.2, 0.5 and 0.8) that mainly differ during times when a more oxygenated state was predicted for the oceans. Given that the modern ocean has γ = 0.4, we expect lower γ values if Mn-oxides contributed less to Mo removal in the oxygenated part of the Late Ordovician oceans than today. The predicted δ 98 Mo values of early-mid Hirnantian seawater are between 1.5‰ (γ = 0.2) and 2.3‰ (γ = 0.5).

Comparison of global redox constraints
The theoretical seawater Mo isotope curves (e.g. γ = 0.5) show a negative shift in the M. persculptus Biozone (e.g. from 2.0 to 0.8‰). The available records of Late Ordovician euxinic shales at Dob's Linn (Scotland), Nanbazi and Wangjiawan (South China) and Libya do not independently confirm this shift, but because euxinic sediments carry δ 98 Mo values that are similar to or lower than the contemporaneous seawater composition, we find that the Mo isotope curves are consistent with the theoretically predicted trajectory from Anticosti Island (denoted δ 98 Mo δ238U ) overlain by variably negative isotope offsets similar to that observed in open euxinic settings today. We conclude that most of the δ 98 Mo variability in the shale records should be ascribed to variations in the local offset between seawater and sediment, and therefore the best constraints on the redox evolution of the global ocean come from the U isotope record at Anticosti Island.
Interestingly, local redox proxies at Dob's Linn signal a persistent shift from oxic to anoxic deposition approximately in the middle M. perscupltus zone, when U isotopes signal a global expansion of marine anoxia (Figs. 3 and 4). The fluctuations in the δ 98 Mo record at Dob's Linn are more dramatic at this time perhaps reflecting a combination of variable local Mo isotope offset (up to 2‰) and a negative shift in global seawater composition (up to − 1.2‰), or perhaps representing a stochastic effect associated with slightly higher sampling resolution in this interval. Thus, the available paleo-redox proxy data are consistent with a global expansion of marine anoxia during LOME2 from a more oxygenated state.

Was there a global oceanic anoxic event during LOME1?
It has been previously proposed that an expansion of euxinia in the oceans in combination with a contemporaneous sea level drop compressed the zone with intermediate oxygen concentrations in the water column and that contributed to the killing during LOME1 (Hammarlund et al., 2012). Yet, the sporadically high δ 98 Mo values (up to 2.0‰) observed in D. mirus and M. extraordinarius zones in South China are at odds with expansive euxinia at this time, since sediments always carry lower δ 98 Mo than overlying seawater. In essence, high seawater δ 98 Mo values require that (at least) one of the major Mo sinks was associated with large Mo isotope fractionation relative to contemporaneous seawater.

Shallow-water euxinia during LOME1
In the context of the extensive marine anoxia in the early Paleozoic (Dahl et al., 2010b;Sperling et al., 2015b;Stolper and Keller, 2018;Wallace et al., 2017) at a time of potentially rising atmospheric O 2 levels (Edwards et al., 2017;Lenton et al., 2016) and evidence for fluctuating local redox conditions both at Dob's Linn and on the Yangtze shelf, we explored expansion of mild and intermittent euxinia with large isotope fractionation in shallow oceans settings as an alternative solution to a large oxic Mo sink during LOME 1. As a driver, we envision marine anoxia to prevail during periods of warmer climate and/or at time of highest productivity that might be paced by orbital variations in solar insolation.
To do this, we used a dynamic mass balance model for the coupled evolution of the marine Mo and U cycles (Dahl et al., 2019). In order to stabilize the oceanic Mo and U inventories when burial fluxes change, the model assumes that the burial fluxes in all settings (F i ) scale with oceanic inventory (M); i.e. F i ~ M i α where α = 0.33-1.00 has been suggested for the Mo cycle (Dahl et al., 2011;Reinhard et al., 2013). Once euxinia expands into shallower parts on the continental shelves, Fig. 5. Comparison of global marine redox constraints inferred from Mo and U isotope data in all studied sections. The model-predicted δ 98 Mo curve, derived from the δ 238 U data at the Laframboise section, is shown in blue. Geochemical data from the Wangjiawan and Nanbazi sections in South China (Zhou et al., 2012, Zhou et al., 2015. The local redox conditions are inferred from Fe speciation data and Mo enrichments, where euxinic/borderline ferruginous conditions (black) are distinguished from oxic conditions (white) by > 25 ppm Mo and/or Fe HR /Fe T > 0.38 and Fe P /Fe HR > 0.7 (Scott and Lyons, 2012;Dahl et al., 2013b).
we predict that the euxinic Mo and U burial fractions (f Mo,EUX, f U,EUX ) increase and produce larger Mo and U isotope offsets from overlying seawater than in restricted euxinic settings (Brüske et al., 2020;Bura-Nakić et al., 2018;Lu et al., 2020). Empirically, we know from modern settings (Fig. 6) that the seawater-sediment isotope offset for both Mo and U (Δ Mo,EUX, and Δ U,EUX ) is greater in more open settings due to more incomplete Mo and U removal from the water column during sediment burial. The coupling of Δ Mo,EUX, and Δ U,EUX is expressed in a linear negative correlation because sediments are preferentially enriched in the heavier U isotopes and lighter Mo isotopes due to fundamentally distinct fractionation mechanisms (mass dependent fractionation for Mo vs. nuclear volume fractionation for U).
To fit the high δ 98 Mo values near the Katian-Hirnantian boundary in South China and large spread of δ 98 Mo values in the various euxinic settings worldwide, we explored a scenario where the euxinic sink expanded periodically; for example in response to Milankovitch-driven warmer climate states. Warmer ocean temperatures may have contributed to subsurface anoxia due to lower oxygen solubility, increased marine O 2 consumption by higher rates of aerobic respiration and enhanced water column stratification (Sperling et al., 2015a;Ulloa et al., 2012). A short period of warming climate during LOME1 has been suggested based on an inflection in the oxygen isotope record (Bond and Grasby, 2020). Fig. 7 exemplifies a scenario where mildly euxinic sediments are a substantial marine Mo sink oscillating with a period of 400 kyrs with an amplitude of 300% relative to an initial Late Ordovician ocean state with an anoxic deep ocean (f Mo,eux = 0.25; f Mo,SAD = 0.75, f Mo,ox = 0) (e.g., Dahl et al., 2010b;Lu et al., 2020). To simulate the effect of periodic shallow euxinia, we forced the model with a variable Mo isotopic offset in the euxinic sink, Δ Mo,EUX , starting at model time t = 1 Myr with oscillations of 405 kyr as an example of insolation driven climate change. The oscillations varied from 0.5‰ to 2.0‰ consistent with the full range observed in shallow euxinic settings today, and the U isotope fractionation (Δ U,EUX ) to respond accordingly as predicted from modern open ocean settings (Fig. 6). The latter mimics 'shallow water euxinia' in areas where water renewal from the oxygenated zone is high and H 2 S concentrations remain low. The best modern analogue settings are the upwelling zones off the west coast of South America and South Africa where hydrogen sulfide occasionally breaches into the water column during productive periods (Schunck et al., 2013;Weeks et al., 2004). Yet, we envision shallow euxinia was far more extensive in the  with an amplitude factor of 3 and a period of 400 kyr relative to a more anoxic ocean state than today without Mn-oxide deposition in the deep ocean (f Mo,eux = 0.20, f Mo,SAD = 0.75). The magnitude of isotope fractionation in the euxinic sink (for Mo, green curve) followed the open-basin trend (see Fig. 6) until t = 2.0 Myr (LOME2), when a transition to a more persistent deep-water euxinic state occurs (25-fold greater euxinic Mo sink) with isotope fractionations following the restricted euxinic trend. In this model scenario, the marine Mo and U residence time are ~100 ky and 340 kyr, respectively, as predicted if the sensitivity of sedimentary Mo and U burial fluxes to their respective marine inventories scale with exponent β = 0.33 (see text for further details). This situation leads to a heterogeneous Mo distribution in the ocean after LOME2 with shorter residence time than the oceanic mixing time (~1.5 ky) (Sarmiento and Gruber, 2006). The vertical dotted lines represent modern oceanic input (red) and modern seawater (yellow) for comparison, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) pastboth in space and time.
The studied shelf settings at Dobs Linn, Nanbazi and Wangjiawan all preserve geochemical evidence for occasional local water column euxinia near the Katian-Hirnantian boundary (Hammarlund et al., 2012;Zhou et al., 2015) and therefore we consider these sites as representative of regional expressions of periodic euxinia that were common in the oceans at this time. Our model demonstrates that periodically high seawater δ 98 Mo values reaching values similar to today's ocean can occur in a far less oxygenated ocean, consistent with the extremely high δ 98 Mo values found in the late Katian and early Hirnantian part of the sections from the Yangtze Shelf Sea in South China. High amplitude oscillations of seawater δ 98 Mo occur, when the Mo residence time is shorter than the period of the forcing, and it also requires that the largest isotope fractionation occurs in phase with peak euxinia, as expected if sulfidic waters invaded shallower parts of the ocean where sedimentation rates are naturally higher. The largest δ 98 Mo swings are observed when the Mo residence time in the ocean is substantially shorter than the period of shallow-water euxinic oscillations; for example in the chosen parameterization (Table 2) δ 98 Mo oscillates from 1.6 to 2.1‰ for τ Mo = 100 kyr and from 1.6 to 1.9‰ for τ Mo = 250 kyr. The shorter Mo residence time in this simulation results from a lower Mo burial flux sensitivity to oceanic Mo inventory (β = 0.33 for τ Mo = 100 kyr, and β = 1.0 for τ Mo = 250 kyr). The peaks in seawater δ 98 Mo are lagging behind peak shallow euxinia by ~50 kyr (vs. ~70 kyr). Thus, there is a tight relationship between periods of extensive shallow euxinia, the marine Mo residence time, and the actual time represented by the sedimentary records. We find that periodic shallow-euxinia is only a viable solution to the Hirnantian Mo and U isotope records if the driving force for the shallow water euxinia (e.g. Milankovitch climate cycles of ≤400 kyr) is sufficiently short compared to the length of the Hirnantian stage, e.g. 1.4 Myr or 0.4 Myr, (Cohen et al., 2013;Ling et al., 2019), and the Mo residence time in the oceans (< 250 kyr).
The model in Fig. 7 also shows a transition into a more stable euxinic state, at model time t = 2.0 Myr, simulating plausible conditions during LOME2. Here, we find that a further ~25-fold expansion of the euxinic Mo sink would be consistent with the U and Mo isotope constraints across the Ordovician-Silurian boundary (Fig. 5). Even if the size of euxinic sink continued to oscillate, muted isotope fractionations representative of euxinia extending into deeper and more restricted parts of the oceans would lead to more stable low δ 98 Mo and δ 238 U values in post-LOME2 seawater.

LOME1
The first extinction interval (LOME1) occurs near the Katian-Hirnantian boundary at a time when marine species richness had been declining for several million years (Rasmussen et al., 2019). LOME1 affected most major animal clades across ecological niches, such as graptolites (pelagic) and trilobites that swam in the pelagic zone and brachiopods and corals in the benthic zone (Harper et al., 2014). LOME1 appears to have been selective in targeting deep-water, higher latitude faunas (Finnegan et al., 2016;Rasmussen and Harper, 2011a). Specifically, the main victims were species inhabiting shallow cratonic seaways that presumably drained as Gondwanan ice sheets grew, as well as species adapted to relatively warm temperatures (Finnegan et al., 2016;Sheehan, 1973). There is geochemical evidence, based on Fe-speciation and Mo concentrations, for weak or intermittent euxinia developed during LOME1 in the Yangtze Shelf Sea, South China (Zou et al., 2018), and a similar tendency is also observed at Dob's Linn, Scotland (Bond and Grasby, 2020;Hammarlund et al., 2012). Due to the large sea level drop at this time, the marine euxinia may have expanded far offshore into depositional settings that are not preserved in the geological record (Hammarlund et al., 2012), or perhaps as suggested here, intermittent or mild euxinia invaded shallower waters if the chemocline moved upwards and closer to the well-mixed and oxygenated part of the surface ocean (for reasons that we discuss further below).
Our new model suggests episodic instability of marine euxinia that might have occasionally also affected deeper parts of the shelf-slope settings. This finding would explain the preferential removal of deepwater species, such as the Foliomena brachiopod fauna that disappeared during the first pulse of the extinctions (Finnegan et al., 2016;Rasmussen and Harper, 2011b). The combination of expanding euxinia and sea level fall would further reduce the zones of habitable niches, increasing competition as shelf volume narrowed. Habitat destruction as a cause for the late Katian extinction pulse has been suggested (Brenchley, 1984;Rasmussen et al., 2019), and the most species-rich regions, such as island arc settings near continental margins, disappeared as continents started to amalgamate as a consequence of the changing Late Ordovician paleogeographic configuration (Rasmussen and Harper, 2011a). Habitat destruction through slow plate tectonic processes must have been ongoing for millions of years, gradually destroying habitable niches. This trend is mirrored by large global richness datasets, showing that extinctions commenced much earlier than the Hirnantian (Fan et al., 2020;Rasmussen et al., 2019;Sepkoski, 1997) and that conditions got worse as glaciations expanded and the sea level dropped in the late Katian (Brenchley et al., 2006;Haq and Schutter, 2008;Nielsen, 2004) displaced shallow-water faunas farther offshore while euxinia at the same time removed the Foliomena-type faunas that previously had occupied the deepest-water settings (Finnegan et al., 2016;Rasmussen and Harper, 2011b). As habitable space reduced, increased competition favored the incumbent lower shelf-upper slope faunas over the now displaced shallow-water faunas. This is the classic view on this first strike of the extinctions (Sheehan, 1975). Expanded deep-water euxinia that even invaded shallower settings is in line with this hypothesis as it increases the area of inhabitable regions considerably.
The trigger for LOME1 has been debated. Paleotemperature records suggest the cooling trend into an icehouse state started in the Middle Ordovician (Finnegan et al., 2011;Rasmussen et al., 2016;  2008; Vandenbroucke et al., 2010;Young et al., 2009). Traditionally this cooling has been linked with atmospheric CO 2 decline that can be explained by high rates of atmospheric CO 2 removal via silicate weathering in the aftermath of the Middle Ordovician emplacement of the Taconics volcanic terrain (e.g. Swanson-Hysell and Macdonald, 2017). The drawdown of CO 2 might have also been assisted by evolution of non-vascular plants (Lenton et al., 2012) and a cessation in volcanic outgassing rates (Young et al., 2009), or both. A rapid increase in atmospheric CO 2 removal via silicate weathering in the late Katian leading into the Hirnantian glaciation has been suggested based on a positive shift to more radiogenic initial 187 Os/ 188 Os values at Dob's Linn (Fig. 3) (Finlay et al., 2010). Osmium isotopes indicate more weathering of radiogenic continental crustal rocks, but this proxy is primarily sensitive to the provenance of the local rock type delivering solid weathering products to the basin and much less sensitive to the absolute weathering flux (e.g. riverine Si, Mg or Ca fluxes, Fig. 3). Thus, the Os isotope signal reflects changes in the local weathering source and may not have a global significance.
Lithium isotopes in marine sediments can constrain changes in the global weathering regime further, and specifically the global weathering intensity (Pogge von Strandmann et al., 2017), defined as the ratio of chemical weathering to total denudation rate (W/D) (Dellinger et al., 2015).. There is evidence that climatic cooling during LOME1 was driven by a decline in volcanic CO 2 outgassing rather than enhanced CO 2 removal (Pogge von Strandmann et al., 2017), and we suggest that the Li isotope records also points to enhanced weathering during LOME2. This would fit well with our modeling results, since oceanic P input from land could have fertilized the oceans and promoted the transition to more stably anoxic ocean conditions.
The lithium isotopic evidence comes from sediments deposited on different paleocontinents in the Late Ordovician oceans (Pogge von Strandmann et al., 2017). There is a positive δ 7 Li shale excursion at Dobs Linn rising below LOME1 and declining into LOME2. Previously, this excursion has been correlated to a positive δ 7 Li carb excursion in carbonates on Anticosti Island based on carbon isotopic chemostratigraphic correlation (Pogge von Strandmann et al., 2017), but the biostratigraphic correlation in Figs. 3 and 4 shows that the positive δ 7 Li carb excursion is confined to the M. persculptus Zone and is not positively correlated to the positive excursion in Dobs Linn shales. By comparison to modern river systems, more positive δ 7 Li values in detrital material (and thus in δ 7 Li shale ) occurs at lower weathering intensity (Dellinger et al., 2017). There is no δ 7 Li carb record across LOME1 to test whether lower weathering rates was a global phenomenon that changed seawater composition occurred at this time.
A negative correlation between δ 7 Li carb and δ 7 Li shale appears at the onset of LOME2, suggestive of intensified weathering when oceans entered a more stable anoxic state. Comparing to modern river systems, the δ 7 Li value of the dissolved riverine discharge increases with weathering intensity until an intermediate level (W/D ~ 0.05) (Dellinger et al., 2015). Therefore, an increase from low to intermediate weathering intensity can explain the positive δ 7 Li excursion in seawater and marine carbonates (δ 7 Li carb ) observed during the onset of LOME2, and simultaneously produce a δ 7 Li decline of detrital Li in marine shales (δ 7 Li shale ) in a manner consistent with the inferred negative δ 7 Li carb -δ 7 Li shale correlation at this time. In summary, enhanced weathering causing ocean fertilization and marine anoxia is compatible with existing records for LOME2, but the Early Hirnantian and the Silurian lithium isotope records are still too sparse to evaluate whether the negative correlation (fingerprinting global weathering change) extends over longer periods of time.
Apart from the inorganic CO 2 sources and sinks, the organic sinks (increasing organic carbon burial) and loss of organic CO 2 sources (declining oxidative weathering) might also contribute to CO 2 draw down and cooling in the Late Ordovician. Increased organic carbon burial fluxes should also simultaneously increase atmospheric O 2 , as suggested at this time (Adiatma et al., 2019;Edwards et al., 2017;Lenton et al., 2016). Thus, a wide range of potential triggers and kill mechanism(s) for the first phase of the Late Ordovician Mass Extinction still exist and can be constrained by more U and Li isotope data.

LOME2
Whereas the first phase of the Late Ordovician extinctions exhibit biogeographical and bathymetric selectivity (Finnegan et al., 2016), the second phase was more taxonomically selective than the first phase. In particular, the shelly benthos, such as brachiopods, corals and trilobites suffered preferential losses. Notably, brachiopods and trilobites had adapted to the colder climate, evolving into the globally occurring successful Hirnantia-Mucronaspis Fauna. After the extinction came the Edgewood-Cathay fauna which albeit less diverse thrived in the latest Ordovician shallow-water zones globally (Rong et al., 2020). The mobile benthos (including some conodonts) and plankton (including graptolites) on the other hand, quickly radiated as temperatures rose even though these clades had been among those hardest hit by the first wave of the extinctions (Harper et al., 2014).
Expansive marine anoxia is commonly evoked as a kill mechanism during Phanerozoic extinction events, and paleo-redox proxies also suggest that animals disappeared during LOME2 in part due to oxygendeficiency, sulfide-poisoning in the water masses, or both (Bond and Grasby, 2020;Hammarlund et al., 2012;Stockey et al., 2020;this study). Dark graptolitic mudstones rich in pyrite and organic matter are widely deposited during the transgression (in shallow settings) after the melting of the Hirnantian ice sheets developed on continental shelf-like areas across North Africa, Southern Europe and the Middle East (Berry and Wilde, 1978). Local redox proxies also suggest anoxic water masses returned in shelf settings on Baltica (Billegrav-2, Denmark), on the margin of Northern Gondwana (Carnic Alps, Austria), in outer shelf settings of Avalonia (Dob's Linn, Scotland), and in the deeper parts of the Yangtze Shelf Sea (Wangjiawan, South China) (Hammarlund et al., 2012;Zhou et al., 2015;Zou et al., 2018). Further, global redox proxies based on molybdenum and uranium isotope analyses support that the global fraction of sediment burial in anoxic parts of the oceans increased in the late Hirnantian and across the Ordovician-Silurian boundary and that anoxic ocean state and a cold climate prevailed for several million years into the Silurian (Bartlett et al., 2018;Dahl et al., 2010b;Finnegan et al., 2011;Stockey et al., 2020;Trotter et al., 2016).
An abrupt intensification of chemical weathering during LOME2 might have been forced by the emerging terrestrial ecosystems (Dahl and Arens, 2020). The earliest evidence of vascular plants comes from trilete spores in Late Ordovician deposits (~450 Ma) and the oldest body fossils of vascular plants (e.g. Cooksonia) are from the Wenlock Series (~430 Ma). Vascular plants with shallow roots have the capacity to retain finer soil grains on land and therefore expose more mineral surface area to weathering fluids (McMahon and Davies, 2018). Also, water transport over the continents in the Early Paleozoic would have been faster in the absence of root plants and absence of meandering rivers compared to that subsequent evolving terrestrial ecosystems (Gibling et al., 2014). Therefore, we expect Li uptake into secondary minerals to have increased as weathering fluids interacted with more secondary clays forming in the terrestrial environment (i.e., lithium weathering would have been highly incongruent before vascular ecosystems with shallow roots evolved when inferred W/D was low).
Enhanced chemical weathering during LOME2 could have fertilized the ocean with phosphorous and triggered euxinia in shallow parts of the oceans, in a fashion also suggested for anoxic oceanic events in the Late Devonian when seed plants evolved and invaded drier upland areas (Algeo and Scheckler, 1998;Joachimski and Buggisch, 1993;Zhang et al., 2020a). Yet, the positive δ 7 Li carb excursion (intensified weathering) at Anticosti Island is short in duration and extensive marine anoxia prevailed for a long (~3 Myr) period long after weathering intensity declined (e.g. Stockey et al., 2020). Thus, expansive anoxia into the Silurian must have crossed a tipping point that reinforced marine anoxia via a runaway feedback; for example via a global scale anoxia-phosphorous feedback (Cappellen and Ingall, 1996).

Conclusion
In conclusion, results from our utilization of both Mo and U isotope data from four paleocontinents in a coupled dynamic oceanic mass balance model suggests that the expansion of sulfidic anoxia (intermittent or through low sulfide concentrations) into well-ventilated parts of the shallow oceans can produce δ 98 Mo variations in accordance with the fluctuations observed in the sedimentary records near the level of LOME1. Yet, Mo isotopes in euxinic sediments are a less reliable record of global seawater composition than is the U isotope composition in marine carbonates. For the LOME2, the model demonstrates that globally expansive ocean anoxia is best constrained by δ 238 U in carbonates from Anticosti Island, and that the Mo isotope record is consistent with the U isotope record. Thus, oceanic oxygen loss might have played a role during the first extinction phase and more certainly did during the second extinction phase in the Late Ordovician.

Declaration of Competing Interest
None.