Research ArticleGEOCHEMISTRY

A long-term record of early to mid-Paleozoic marine redox change

See allHide authors and affiliations

Science Advances  07 Jul 2021:
Vol. 7, no. 28, eabf4382
DOI: 10.1126/sciadv.abf4382


The extent to which Paleozoic oceans differed from Neoproterozoic oceans and the causal relationship between biological evolution and changing environmental conditions are heavily debated. Here, we report a nearly continuous record of seafloor redox change from the deep-water upper Cambrian to Middle Devonian Road River Group of Yukon, Canada. Bottom waters were largely anoxic in the Richardson trough during the entirety of Road River Group deposition, while independent evidence from iron speciation and Mo/U ratios show that the biogeochemical nature of anoxia changed through time. Both in Yukon and globally, Ordovician through Early Devonian anoxic waters were broadly ferruginous (nonsulfidic), with a transition toward more euxinic (sulfidic) conditions in the mid–Early Devonian (Pragian), coincident with the early diversification of vascular plants and disappearance of graptolites. This ~80-million-year interval of the Paleozoic characterized by widespread ferruginous bottom waters represents a persistence of Neoproterozoic-like marine redox conditions well into the Phanerozoic.


Earth’s oceanic and atmospheric oxygen levels are commonly assumed to have markedly increased in two steps, one in the Paleoproterozoic [~2450 to 2220 Ma (million years)] and one in the Neoproterozoic (~800 to 540 Ma) (1). During this second increase, researchers have historically suggested that atmospheric oxygen levels reached approximately modern values, followed by relatively small fluctuations through the Phanerozoic. This conceptual model has been challenged by studies that instead proposed either transitory oxygen fluctuations or limited Neoproterozoic-Paleozoic oxygen changes, with a delayed transition to more permanently modern levels of oxygenation later in the Paleozoic, likely in the Devonian (26). Such hypotheses fit well with the geological observation that marine black shale deposits are common throughout the early and middle Paleozoic (7, 8).

Set against this backdrop of changing atmospheric oxygen levels, it is also clear that the nature of anoxic water masses has changed through time. Specifically, the balance between oxic, euxinic (sulfide present; iron-limited), and ferruginous (ferrous iron present; sulfide-limited) bottom waters has shifted (3, 911). This balance is intimately tied to nutrient cycling and primary productivity (1214) and to biological responses to upwelling anoxic water masses. Neoproterozoic anoxic water masses have been shown to be dominantly ferruginous over long time intervals (>10 Ma) (3, 9, 11). Note that although intervals such as the Neoproterozoic may have been relatively more ferruginous, the transition between redox states within each basin is driven by a combination of local factors and global seawater chemistry, and both anoxic states have coexisted in time throughout Earth history. There is evidence for both euxinic and ferruginous conditions at points in the early Paleozoic (3, 9, 1520)—with evidence for euxinia perhaps more pronounced in the early and late Cambrian—but there remains no long-term record of the nature of Paleozoic anoxia. These ambiguities in the timing of global oxygenation and our lack of data regarding background state redox conditions have important implications for understanding the role of environmental change in both the Cambrian radiation (2123) and other early Paleozoic radiations and extinctions (24, 25).

Tracking redox change from the Neoproterozoic through the Paleozoic becomes complicated because of the general design of geochemical studies across this broad swath of time. In the former, researchers have commonly interrogated sedimentary successions that cover entire periods or even most of the era (9, 11, 2628). Because such studies sampled at relatively consistent spacing through long stratigraphic sections, there is an increasingly robust record of the Neoproterozoic “background state.” In contrast, Paleozoic redox studies have commonly focused on specific geochemical perturbations, such as the Steptoean Positive Carbon Isotope Excursion (SPICE) (15), or mass extinctions, such as the end-Ordovician (29). This provides detailed insight into these events but also results in a series of isolated snapshots without a firm understanding of the Paleozoic background state. Here, we conduct the first Paleozoic shale geochemical study at a temporal scale and sampling style equivalent to Neoproterozoic studies, focusing on the Road River Group exposed along the Peel River, Yukon, Canada (Fig. 1). This section is an almost continuous record of Furongian (late Cambrian) to Middle Devonian deep-marine environments (slope and basin floor) that provides a unique long-term window into early Paleozoic redox evolution (30).

Fig. 1 The Peel River study section.

(A) Paleogeographic map of Laurentia (ca. 480 Ma) with the study area depicted in the box [from (30)]. (B) Early Paleozoic paleoenvironmental reconstruction of Yukon, Northwest Territories, and British Columbia [modified after (30)]. The northwestern Laurentian margin was divided into a series of shallow-water carbonate platforms (blue) and deep-water shale basins (gray). The Road River Group was sampled at the Peel River locality in the Richardson trough. (C) The majority of the Road River Group consists of unbioturbated interbedded organic-rich mudstone/shale and lime mudstone with local diagenetic chert replacement. Shown here are Tremadocian dark gray to black shale and lime mudstone of the Mount Hare Formation. Seated geologist is at ~38 m in section J1727 [see (30) for detailed sedimentological discussion of the Peel River section]. (D) Trace fossils occur sporadically throughout the Road River Group, but the primary interval of extensive, continuous bioturbation occurs in the lower to middle Katian (Upper Ordovician). Bedding plane view of heavily bioturbated Katian strata from 81 to 93 m in J1518; 30-cm-long hammer for scale. MCE, Misty Creek embayment; MRE, Meilleur River embayment. Photo credit: Erik Sperling, Stanford University.

Geological background

During the early Paleozoic, northwestern Canada formed a promontory along the northern margin of Laurentia. The region was dominated by shallow-water carbonate platforms and deep-water intra-platformal troughs and basins, such as the Richardson trough and Selwyn basin (Fig. 1A) (31). These troughs accumulated deep-water strata of the Road River Group and were analogous to those dissecting the modern Bahamas carbonate platform, such as Tongue of the Ocean [albeit wider; more generally, this represents a narrow basin situated within a broader carbonate platform built on thick, rifted continental crust (30, 32)]. The Richardson trough was connected to the global ocean through its northern opening, via seawater transiting over the Yukon block carbonate platform, and from the Selwyn basin to the south (31), although the extent of this southern connection is debated. The biostratigraphy, sedimentology, and depositional interpretation of the Peel River type section (Fig. 1B) were recently updated following a comprehensive investigation (30). The studied section, covering the Cronin, Mount Hare, Tetlit, and Vittrekwa formations in ascending stratigraphic order, consists of ~2.3 km of interbedded graptolitic shale, chert, lime mudstone, and rare rudstone, all deposited below storm wave base in slope or basin-floor settings (Figs. 1 and 2 and the Supplementary Materials). Field evidence for macroscopic bioturbation in the Road River Group is rare, but sporadic bioturbation occurs straddling the Cambrian-Ordovician boundary, at points between the upper Darriwilian and lower Telychian, and in the lower Pragian (fig. S5). Moderate bioturbation is present in more resistant strata of the Tetlit Formation (fig. S5), and continuous bioturbation is found only in the lower to middle Katian (Upper Ordovician; caudatus through lower complanatus biozones; Figs. 1D and 2 and fig. S5). Confirmed benthic fossils are also rare in the Peel River section, and the only candidates occur near the base of the Road River Group. The first graptoloid graptolite appears 64 m above the base of the Mount Hare Formation, and the last occurs ~60 m below the top of the Vittrekwa. Samples from the Road River Group were combined with new data generated from other localities in Baltica, Avalonia, Gondwana, and Laurentia (>1100 new samples total; see Materials and Methods), as well as the published literature, to examine redox change from the late Cambrian through Devonian.

Fig. 2 Geochemical data from the Road River Group and overlying Canol Formation of the Earn Group.

Stratigraphic heights in meters refers to the full composite section on the Peel River detailed in (30) and continue into the Canol Formation in the RI-07-07A core. Detailed stratigraphic sections and biostratigraphic information can be found in fig. S5. From left to right, the first two geochemical columns represent highly reactive to total iron (FeHR/FeT) and pyrite to highly reactive iron (FeP/FeHR) iron speciation analyses, interpreted using accepted empirical baselines for oxic, anoxic, euxinic, and ferruginous conditions (33). The minority of samples with FeT < 0.5 wt % or FeHR/FeT > 1 (commonly in low-iron samples with higher error) are denoted as open circles. Iron speciation columns are followed by TOC weight percent and authigenic enrichments of redox-sensitive metals Moauth, Uauth, and the Moauth/Uauth ratio. Organic carbon isotopes (δ13Corg), trace metal analyses normalized to TOC, other RSMs, and detailed stratigraphic columns can be found in the Supplementary Materials. Darriwil., Darriwilian; S., Sandbian; Kat., Katian; Rhu., Rhuddanian; Aer., Aeronian; Wen., Wenlock; Shein., Sheinwoodian; Hom., Homerian; G.?, Gorstian?; Ludford., Ludfordian; Prid., Pridoli; L., Lochkovian; Pra., Pragian; Em., Emsian; E.-G., Eifelian-Givetian; Givet.-Frasn., Givetian-Frasnian; Ferr., Ferruginous; Eux., Euxinic. Sporadic bioturbation, rare and small (diameter < 1 cm) isolated burrows; moderate bioturbation, some but not all beds burrowed, with rare burrows up to 1 cm in diameter; intense bioturbation, all beds extensively burrowed, burrow diameter often ≥1 cm.

Tracking Paleozoic redox changes

Iron speciation can fingerprint ancient anoxic bottom waters based on enrichments in highly reactive iron relative to total iron above baseline values (FeHR/FeT). For anoxic samples, it can further distinguish between euxinic and ferruginous water columns based on the proportion of highly reactive iron that has been pyritized (FeP/FeHR) (33). Note that the iron speciation proxy estimates the relationship between available sulfide and reactive iron; “ferruginous” conditions represent sufficient reactive iron present to titrate available sulfide and do not necessarily denote extremely iron-rich water columns, such as those that permeated the Archean. Likewise, “euxinic” gives no indication of absolute sulfide concentrations. Next, we examined the abundance of redox-sensitive trace metals (RSMs), particularly molybdenum (Mo) and uranium (U). These metals are soluble in oxic seawater but become insoluble through complexation with organic matter, sulfides, and/or other mineral phases under reducing conditions (34). High authigenic Mo enrichments (Moauth) require the presence of appreciable free sulfide, while authigenic U enrichments (Uauth) can also form under strongly reducing but noneuxinic conditions. Therefore, higher Mo/U ratios identify likely euxinic water-column conditions (35) in a manner independent from iron speciation. While each of these proxy results represents an individual observation of an ancient local water column, we use a reweighted bootstrap technique (36) that minimizes spatial and temporal sampling biases to gain the best understanding of global redox trends in Paleozoic oceans.


With few exceptions, geochemical redox proxies show evidence for bottom water anoxia throughout deposition of the Road River Group and overlying Canol Formation (Fig. 2, fig. S5, and table S1): >98% of samples from the Peel River section have FeHR/FeT ratios >0.38, the value commonly taken to indicate an anoxic water column (33). Moreover, 82% have >1 part per million (ppm) Moauth and 90% have >1 ppm Uauth, demonstrating that most samples show evidence for authigenic RSM enrichment under a reducing water column (Fig. 2 and fig. S5). The few samples with low to absent enrichments in Moauth and Uauth are clustered at specific intervals in the lower to middle Katian, Sheinwoodian/Homerian/Ludlow, and lower Pragian, all of which record bioturbation (fig. S5).

With respect to the nature of the anoxic water column, the Devonian portion of the Peel River section has substantially more euxinic samples based on iron speciation compared to either the Ordovician or Silurian, and samples from above the Pragian falcarius graptolite Biozone were significantly more likely to be euxinic than samples below (Fig. 2, Table 1, and table S2). Moauth/Uauth (hereafter Mo/U) ratios also show a clear and significant change from generally <4 to >4 at and above this level, providing independent support for a transition to euxinic conditions in the Richardson trough at this time. Classification tree analyses confirm that the best splits in the Peel River section are at 409 Ma (iron speciation) and 401 Ma (Mo/U).

Table 1 Summary table of redox geochemical data from shales deposited beneath anoxic water columns through the Paleozoic.

Analyses are separated as Peel River [main Peel River composite section as described (30) plus data from the RI-07-07A core] and our global dataset (incorporates new data from the Road River Group, new worldwide sampling of Ordovician-Devonian shale, and available published data). All samples analyzed here have iron geochemical data for anoxic deposition. Samples were analyzed for the proportion of euxinic samples, Corg/P, and proportion with authigenic P enrichments (P > 1000 ppm) across three time bins (>467, 467 to 408, <408 Ma; chosen based on visual inspection of redox data at the Peel River locality). Proportions of euxinic samples and samples with authigenic P enrichments in each time bin are represented as percentages and were analyzed for significance with chi-square tests. Corg/P ratios represent medians and were analyzed for significance with Mann-Whitney U tests. Significance column denotes differences that are significant at a Bonferroni-corrected P value of <0.002. Letters correspond to the three age bins in ascending stratigraphic order; bins that share the same letter are not significantly different (full P values and test statistics are provided in table S2). Best splits column shows the first three splits in classification or regression tree analyses, which provides an unbiased prediction of the age of the biggest changes in the dataset.

View this table:

There is potential for some of the signal for early Paleozoic ferruginous conditions to be driven by minor oxidative weathering of outcrop samples, which will oxidize pyrite to iron oxides (9) and preferentially remobilize U relative to Mo (37). However, these strata are exposed in a sub-Arctic river canyon, with relatively minor oxidative weathering and very similar weathering/exposure throughout the stratigraphic section. Further, we only selected visually unweathered samples. Many laminated early Paleozoic samples also show substantial iron in the acetate extraction [which represents ferrous iron (9)] and low total Moauth enrichments irrespective of Mo/U ratios. Last, we find that these patterns hold even if all iron oxides are considered to result from oxidized pyrite [fig. S3, Table 1, and table S2; see (19) for further discussion of this test]. Thus, the Devonian redox shift at the Peel River is unlikely to be an artifact of weathering.

The observation of increasing occurrence of Devonian euxinic samples is supported by a new globally distributed dataset. Classification tree analyses of this dataset also indicate the biggest changes at 387/409 Ma (iron speciation) and 402 Ma (Mo/U) (Fig. 3, Table 1, and table S2). Although many time bins are dominated by data from the Road River Group, dominantly ferruginous conditions in Ordovician and Silurian samples are documented in other regions (fig. S6), and our binned analysis uses a reweighting algorithm that accounts for spatial and temporal sampling intensity (36). Reweighted bootstrap means for the global dataset indicate that ~5 to 40% of sampled Ordovician–Early Devonian (Lochkovian) anoxic shales record a euxinic iron speciation signature. The exception to this dominance of ferruginous samples is in the early Silurian Rhuddanian Stage. This interval has a higher proportion of euxinic samples both globally (Fig. 3) and in our regional analyses (fig. S6), consistent with previously published local and global redox proxy data for this interval (17, 18, 29, 38). The proportion of sampled shales deposited under euxinic bottom water conditions then rises in the Pragian; bootstrap means for Emsian and younger time bins show that ~50 to 80% of anoxic shales have a euxinic iron speciation signature (Fig. 3). The proportion of samples estimated as euxinic from Mo/U ratios also rises at this time (Fig. 3). These changes in the biogeochemical nature of anoxia in the global dataset after the Pragian (~408 Ma) are substantial and significant in analyses of iron speciation, iron speciation with all oxides considered to represent weathered pyrite, and Mo/U ratios (Table 1 and fig. S2).

Fig. 3 Geochemical patterns in Paleozoic anoxic shales and the relationship to graptolite diversity and the rise of land plants.

Analyses incorporate new data from the Road River Group, new worldwide sampling of Ordovician-Devonian shale, and available published data. Samples are plotted in 10-million-year time bins, with circles representing the mean of 1000 bootstrap replicates using a reweighted bootstrap algorithm (36) accounting for spatial and temporal sampling density. Error bars represent 2 SD of bootstrap means. Number of samples in each time bin shown as histogram above plots. (A) Proportion of sampled anoxic shales (FeHR/FeT > 0.38) where bottom waters were euxinic (FeP/FeHR > 0.7) based on iron speciation data. (B) Proportion of anoxic samples where Moauth/Uauth > 4 (35), with higher ratios indicating likely euxinic conditions. Note that these two upper plots specifically represent the proportion of shelf, upper slope, and epeiric seaway samples commonly analyzed by shale geochemists and not a representation of the global extent of ferruginous versus euxinic conditions. (C) Proportion of samples where the total phosphorus pool is likely dominated by authigenic phosphorus enrichments (rather than detrital P), determined as P > 1000 ppm. (D) Graptolite species richness follows (56) and the fossil record of land plants follows (62).

While these analyses demonstrate a clear change through time that likely had strong impacts on other biogeochemical cycles (12, 13, 39), it is important to note that the results shown in Fig. 3 demonstrating that ~50% of sampled mid-Devonian anoxic water masses were euxinic cannot be read as a literal record of the global ocean. Previous modeling work has suggested that even during times of “more pronounced” euxinia, sulfidic waters still likely covered a relatively small portion of the global seafloor (concentrated in productive outer shelf and upper slope settings or shallow epeiric seas) (40). Thus, while these results do not speak to conditions in the abyssal plain (e.g., >3 km depth), these data do document redox shifts in the environments that represent the habitat for the majority of animal life and the largest sinks for reduced components in the carbon, sulfur, and other elemental systems (41).

Both the Peel River section and the global dataset also show conspicuous changes in phosphorus contents through time. Classification tree analyses suggest that the largest changes in P contents are at ~455 Ma (Sandbian) for both datasets. In the global dataset, anoxic samples older than 455 Ma have an average of 2093 ppm P (SD = 3616 ppm), whereas samples younger than 455 Ma have an average of 667 ppm (SD = 1534 ppm). However, as indicated by the large SDs, averages do not capture the full story. Lower and lower Middle Ordovician samples have a bimodal distribution of P contents (e.g., fig. S4), and the change is best characterized as a shift in the proportion of samples with high/authigenic P enrichments (>1000 ppm) versus samples with low/detrital P contents. Forty-two percent of samples older than 455 Ma have substantial authigenic enrichment, whereas only 7% of samples after this point do. Total organic carbon to phosphorus ratios (Corg/P) rise substantially in the Early Devonian, from 77 to 238 (global dataset), with best splits in regression tree analyses at 409 Ma for the Peel River section and 367 Ma (followed by 409 Ma) in the global dataset (Table 1).


A long-lived early Paleozoic anoxic basin

Overall, geochemical data agree with field observations of low/absent bioturbation and a lack of benthic fauna and confirm that the Richardson trough bottom waters were generally anoxic [but with brief intervals of oxygenation (42)] for the vast majority of its history. In terms of connecting this marine redox record with estimates of atmospheric O2, most model estimates for the early Paleozoic suggest that atmospheric oxygen levels were lower than modern (4346). Some recent studies, though, have alternatively indicated higher-than-modern levels for the Ordovician and Silurian (24, 47). The documentation of this long-lived early Paleozoic anoxic basin could be reconciled with high atmospheric oxygen by invoking shallower Paleozoic remineralization depths (6) and/or local water mass restriction because of the paleohydrography of the Richardson trough intra-platformal setting. However, such persistent anoxia in this and other globally distributed basins (7, 8) is more consistent with modeling results suggesting lower-than-modern atmospheric O2 during this interval (4346).

Base-level fall in intra-platformal basins such as the Richardson trough may result in increased hydrographic restriction and the development of reducing water column conditions because of longer seawater residence times [e.g., (48)]. Eustatic sea level changes undoubtedly influenced the redox conditions of this basin to some extent, but a striking feature of the data is the relative homogeneity of redox state (at least within the resolution of the proxies) across multiple early Paleozoic sea level changes (25, 49). This includes the profound episode of base-level fall and associated shelf-edge failure in the Darriwilian [Sauk-Tippecanoe megasequence boundary throughout North America, marked by the delivery of decameter-scale debrites of the Aberdeen Member into the basin (30)]. While intermittent euxinia may have developed at this time (see conflicting iron speciation and Mo/U proxy data in fig. S5 through this interval), bottom waters were generally ferruginous before, during, and after this event (Fig. 2 and fig. S5). Most likely, redox conditions on the Richardson trough basin floor were more controlled by fundamental aspects of the Paleozoic Earth system (e.g., low atmospheric oxygen and iron-carbon-sulfur relationships) than by sea level change. An exception to this pattern may be in the Silurian Tetlit Formation, which was deposited during a time of broader shoaling in the Richardson trough, likely related to the Tippecanoe-Kaskaskia megasequence boundary (30, 31). This unit is marked by alternating packages of (i) unbioturbated black shales characterized by geochemical evidence for deposition under a ferruginous water column, with packages of (ii) slightly coarser-grained and commonly turbiditic strata that show ichnological and RSM evidence for moderate oxygenation (Fig. 2 and fig. S5). These likely represent a series of global oxygenation/deoxygenation events (50), modulated locally in the Richardson trough by regional and eustatic base-level changes.

Oceanographic change and the Ordovician radiation

There is no local redox change in the Road River Group through the main pulses of the Ordovician Radiation, an event that is increasingly being linked to changes in marine and atmospheric oxygenation [(23, 24), although this latter study notes that there might be a slight temporal offset between oxygen increase and diversification]. In comparing local marine redox proxy data to hypotheses of atmospheric oxygenation, the inferences drawn from observations at one site can have multiple interpretations. For example, the Peel River site may have been deposited below the chemocline and substantial chemocline deepening resulting from global oxygenation may have occurred above the depths of the Richardson trough basin floor. In this view, the appearance of continuous macroscopic bioturbation in the lower to middle Katian (that is unique within the Peel River section) may represent transient complete oxygenation to the basin floor. However, this oxygenation spans only part of the lower to middle Katian before local bottom waters returned to anoxic conditions. Thus, the Peel River section does not support a substantial redox change through the main pulses of the Ordovician Radiation. That said, it is clear that there were likely considerable oceanographic changes between the early/middle and middle/late Ordovician whose relationship to biotic radiation remain unresolved. These include broad lithological shifts toward chert-rich strata in Laurentia (51), including in the Road River Group, and changes to the sulfur cycle (52).

There is also the start of a protracted change in phosphorus contents of anoxic sediments in the Ordovician. Specifically, a slight rise in Corg/P and a clear drop in the proportion of samples with high authigenic enrichments occurs in the Sandbian (Fig. 3 and Table 1). There are myriad controls on phosphorus contents and Corg/P ratios, including the amount of detrital (unreactive) P input, the C:P ratio of organic matter delivered to the seabed [with terrestrial organic matter having high C:P ratios compared to algal material (44)], the amount of reactive phosphorus remineralized to the water column [with more remineralization in euxinic settings, e.g., (13)], the P inventory of the ocean (12), and the Ca inventory, which controls authigenic carbonate fluorapatite precipitation (CFA) (53). Notably, the Ordovician P shift is largely decoupled from changes in the redox state of sampled shales (Figs. 2 and 3, fig. S4, and Table 1), suggesting that these changes are not due to differential P remineralization under changing ferruginous versus euxinic bottom waters. The common occurrence of moderate to strong authigenic P enrichments in the Early and Middle Ordovician (Fig. 3) could be linked to high dissolved marine Ca concentrations, given that higher Ca levels foster enhanced CFA formation (53). However, the subsequent drop in P enrichments is not perfectly in step with predicted decreases in marine Ca concentrations (54). Changes in P contents during this time could also be linked to greater input of terrestrial derived, phosphorus-poor organic matter or increased solubilization of detrital phosphorus by early land plants [e.g., (44)]. Testing between these and other possible controls will require more detailed studies specifically on phosphorus cycling.

Early Devonian redox change

The Pragian change toward more sulfidic conditions in the Peel River section represents the best direct constraint on this overall trend from the early to middle Paleozoic (Fig. 3). It is possible that the development of euxinia in the Richardson trough was related to a postulated drop in global sea level at this time (49) and associated basinal restriction. However, more prevalent euxinia was maintained in the Richardson trough through subsequent sea level rise in the Early to Middle Devonian (Fig. 2 and fig. S5). Although the precise timing of this redox change must be tested in other Lower Devonian stratigraphic successions worldwide, it is likely that the development of sulfidic conditions would have had direct effects on organisms living in adjacent waters. This Pragian redox change coincides with the taxonomic decline of the dominant zooplankton group of the early Paleozoic, graptolites, and their final extinction in the early Emsian (Fig. 3). Because of this long decline (55, 56), it is unlikely that a single cause was responsible, but it is a markedly different trajectory from once-dominant groups of marine organisms like trilobites or conodonts that experienced long-term taxonomic declines but then met their ultimate demise at era- or period-bounding mass extinctions. Building on previous oceanographic hypotheses for the extinction of the graptolites (55, 57), we hypothesize that the Early Devonian transition toward euxinic conditions in deeper anoxic water masses severely restricted their available habitat space. It has long been recognized that graptolites preferred noneuxinic oxygen minimum zone–type habitats [(57) refers to denitrifying habitats, but broadly, the redox space between fully oxygenated and euxinic was more likely]. The ideal morphological adaptations of graptolites to inhabit such low-oxygen conditions—small, thin body plans with large respiratory surface area—are the converse of those ideally required to inhabit areas with intermittent or chronic sulfide stress (21). This is because both oxygen and sulfide are taken up through diffusive processes, so increasing oxygen uptake increases risk of sulfide toxicity. There is evidence from the Silurian Kozlowskii/Lau extinction event that graptolites were relatively unaffected by initial global deoxygenation but instead specifically affected by the spread of euxinic water masses (50). We posit that similar effects played out on a longer temporal scale over the Silurian-Devonian transition. Specifically, the marked reduction in their preferred dysoxic/suboxic habitat—squeezed by more fully oxygenated conditions at shallower depths and more sulfidic conditions at deeper depths—was likely a major contributor to the final extinction of graptolites.

More generally, upwelling anoxic water masses are often invoked as a kill mechanism in mass extinctions. The effects of ferruginous versus euxinic conditions, however, are not equal. Hydrogen sulfide is a respiratory toxin for animals and is easily able to diffuse across respiratory membranes (58), whereas ferrous iron, while toxic, does not easily diffuse because of its charge and larger size. The toxicity of iron would be even lower under “ferruginous” conditions if ferrous iron only exists at very low concentrations in anoxic but nonsulfidic bottom waters (59). Both clade-level resistance to sulfide and the changing nature of anoxia are underexplored topics in Phanerozoic biodiversity studies.

The documentation of this Devonian redox change begs the question of what may have caused the overall increase in euxinic conditions. One hypothesized driver is the spread of land plants and specifically more deeply rooted vascular plants and forests (60). This was proposed to have led to an increase in weathering and phosphorus flux to the oceans and the development of widespread euxinia. Given that local water column transitions from ferruginous to euxinic conditions can commonly be driven by increased primary productivity and organic carbon loading (26, 61), and primary productivity is limited by nutrient delivery, this scenario broadly fits with the combined terrestrial fossil record and our marine redox record (Fig. 3). This narrative has recently been challenged by modeling studies suggesting that terrestrial plants cannot drive sustained increases in chemical weathering (62). However, the bioavailable phosphorus flux and overall chemical weathering are not strictly linked, and even if the overall weathering flux remained constant, terrestrial plants likely decreased soil pH with increased root depth (63), which would have had a direct effect on P mobilization (64). Considering these complicated relationships between plants, weathering, nutrients, and marine redox (65), any full accounting of Paleozoic Earth system evolution must explain the pronounced change from broadly ferruginous anoxic water masses to more prevalent euxinic conditions that began—at least based on current data—in the Pragian.

Beyond the question of marine-terrestrial “teleconnections” in the Devonian, geochemical and modeling studies have increasingly suggested that changes from dominantly ferruginous to increasingly euxinic water columns likely represent a change in marine nutrient regimes irrespective of shifts in terrestrial phosphorus input. In euxinic water masses, phosphorus is generally returned to the water column via remineralization reactions, thus fueling continued productivity and leading to higher Corg/P ratios. In contrast, under widespread ferruginous water columns, bioavailable phosphorus is more likely precipitated and trapped in the sediments, lowering oceanic phosphorus levels and limiting marine primary productivity (1214, 66). Consistent with this relationship, Corg/P ratios rise substantially between the 467- to 408-Ma and post–408-Ma time bins (Table 1). The proportion of samples with authigenic P enrichments also decreases around this time, albeit potentially slightly later depending on the redox proxy (Fig. 3). While these Early Devonian phosphorus changes appear more tightly coupled to marine redox changes than in the Ordovician, the effects of changing marine Ca contents (53), plant matter input, and changes in terrestrial solubilization likely played a role here as well.

Linkages in the Earth system between primary productivity, organic carbon burial, and atmospheric oxygen indicate that lower productivity oceans are likely responsible for maintaining lower oxygenation of the ocean/atmosphere system (11, 12, 6769). Although there is considerable work to be done in understanding similarities and differences between the Paleozoic and Proterozoic Earth systems, the documentation of widespread ferruginous conditions in the Ordovician, Silurian, and Early Devonian would suggest—by analogy to Proterozoic models—that those time periods were less productive and less oxygenated than the Middle/Late Devonian. This may seem counterintuitive, because it is tempting to associate the development of the most reducing conditions (euxinia) with a lower oxygen world. However, given the biogeochemical relationships described above, euxinia most likely developed as a localized phenomenon in deep-water Devonian basins because of higher primary productivity/export production, despite higher levels of atmospheric oxygen and a more oxygenated surface ocean. Changes to both oxygen and productivity likely affected Paleozoic metazoan ecology. Low levels of atmospheric oxygen have long been discussed as a potential throttle on early animal evolution (22), but low food supply (resulting from low productivity) is also correlated with important ecosystem-level changes in animal body size, abundance, diversity, bioturbation depth, and ecological modes in the modern ocean and could have been a major driver of ancient biological transitions (70, 71). It is likely that many of the overarching environmental controls on animal evolution discussed in reference to the Proterozoic biosphere remained important well into the Paleozoic.

Paleozoic redox change in the context of Earth history

The long-term Paleozoic redox record from the Road River Group paints a picture of a relatively poorly oxygenated world. For much of its history, bottom waters in the Richardson trough were dominated by anoxic and nonsulfidic (ferruginous) conditions. Geochemical analyses of Ordovician and Silurian shale from other localities in Laurentia, Baltica, Gondwana, and other regions tell a similar story, albeit with regional differences. This establishes that widespread and temporally long-lived ferruginous conditions are not only a Proterozoic phenomenon. These data also add to the growing body of evidence that the transition to modern well-oxygenated oceans (with sulfidic conditions where anoxia exists) was not a single state change around the Proterozoic-Phanerozoic boundary. It is unclear from available data how far back “modern” biogeochemistry persisted, but likely there were multiple oscillations between periods with relatively more euxinic versus ferruginous conditions through the Paleozoic and Mesozoic, with a probable ferruginous background state (33). These redox shifts would have had a fundamental effect on other biogeochemical cycles, including nutrient regimes (1214), concentrations of bioessential metals (39), the uranium isotope composition of the oceans (72), and early diagenetic reactions affecting soft-tissue preservation (73). Changing redox states would also have constrained habitable space for animals directly via sulfide toxicity and indirectly via implications for food supply and atmospheric oxygen. Resolving the full Phanerozoic redox history and the relationship to other geochemical cycles and biotic evolution will require the generation of long-term records that comprehensively cover background intervals instead of focusing on stratigraphic “events.”


Experimental design

Our study sampled at 2- to 4-m resolution throughout the composite stratigraphic section at the Peel River locality (30). This corresponds to samples roughly every 190, 260, and 390 ka (thousand years) for the Mount Hare, Tetlit, and Vittrekwa formations, respectively, reflecting the regional sedimentation history. All samples represent fine-grained sedimentary rocks appropriate for redox geochemical studies, even when this is not readily apparent at the scale of our stratigraphic columns. These samples were combined with additional regional and global sampling, including (i) sampling from other Road River Group outcrop localities in northwest Canada (Tetlit Creek, Rock River, and Clearwater Creek) and equivalents in the Cape Phillips Formation, Canadian Arctic Islands; (ii) analysis of the Devonian RI-07-07A drillcore in Yukon, Canada; and (iii) spot sampling from researchers’ personal and museum collections (Yale Peabody Museum and Museum of West Bohemia). In total, 837 new samples were analyzed from the Peel River section, 106 new samples from other Canadian sections/cores, and 178 new spot samples from areas worldwide.

All samples were first ground to a flour in a tungsten carbide shatterbox or ball mill. Samples were then analyzed for iron in pyrite (FeP) using the chromium reduction of sulfur method (74) and for other highly reactive iron pools using standard acetate, dithionite, and oxalate sequential extractions (75). Highly reactive iron (FeHR) was considered the sum of these three sequential iron extracts plus the iron in pyrite (33). Total organic carbon (TOC) contents and organic carbon isotopic data (δ13Corg) for Peel River samples were previously published (30). For TOC contents of museum and spot samples, powder was acidified before combustion in a Carlo-Erba NA 1500 Elemental Analyzer at Stanford University. All Peel River samples shown in Fig. 2 were analyzed for major- and trace-element composition using the LF204 package at Bureau Veritas/AcmeLabs. Samples were digested in a lithium meta/tetraborate flux, with major elements measured by inductively coupled plasma optical emission spectrometry (ICP-OES) (Spectro Arctos) and trace elements measured by ICP mass spectrometry (ICP-MS) (Elan 9000). Chalcophile trace elements were separately measured via ICP-MS following a standard four-acid digestion. Total carbon and sulfur were measured with a LECO combustion analyzer. Major and trace elements for all other samples (including Cronin Formation samples from sections J1728 and T1701 at Peel River; not shown in Fig. 2) were analyzed either by four-acid digestion and measurement via ICP-OES/MS at Bureau Veritas using the MA200 package (instruments as above) or at Yale University on a Thermo-Finnegan Element XR ICP-MS. Estimates of precision and accuracy, laboratory and analysis type for each sample, relationship to previously published studies, and any exceptions to these protocols are described in the Supplementary Materials and tables S1 and S3.

Statistical analyses

Iron geochemical data were interpreted according to accepted baselines, with elevated ratios of highly reactive iron to total iron or elevated iron/aluminum ratios taken to represent deposition under an anoxic water column (FeHR/FeT > 0.38 and Fe/Al > 0.53) (33, 76). High ratios of pyrite to highly reactive iron (FeP/FeHR > 0.7) in anoxic samples were assumed to indicate euxinic water column conditions, with lower ratios representing ferruginous conditions. A minority of samples (e.g., siliceous shales) have relatively low quantities of total iron [<0.5 weight % (wt %)] and may have erroneous iron speciation values (77); these samples are marked by open circles in stratigraphic plots. Authigenic Mo (Moauth) and U (Uauth) were calculated on the basis of Al concentrations and estimates of upper continental crust (78). For Moauth/Uauth, no accepted baselines exist, and a ratio >4 was considered to represent likely euxinic conditions (35), recognizing that while theoretically sound, there is imprecision in this estimate. Sensitivity analyses indicate that results are robust to different thresholds (fig. S7).

Our new iron, carbon, and trace metal data from the Road River Group and other areas worldwide were then combined with upper Cambrian–Devonian shale geochemical data from the literature (4499 samples) to generate an updated global dataset (5619 samples total). The visually apparent geochemical changes in the Road River Group at Peel River (broadly in the lower Darriwilian, estimated at 467 Ma, and in the Pragian, estimated at 408 Ma) were used in statistical hypothesis testing to test whether there were any significant differences in the geochemistry of samples deposited before versus after these points. These analyses only considered samples with iron geochemical data for deposition under an anoxic water column. Chi-square tests were conducted on the proportion of euxinic samples from the pre–467-Ma, 467- to 408-Ma, and post–408-Ma intervals and on the proportion of samples dominated by authigenic (> 1000 ppm) versus detrital (<1000 ppm) phosphorus. Mann-Whitney U tests were conducted on ratios of Corg/P from the same time bins; a Bonferroni-corrected P value of <0.002 in all analyses was considered significant. Classification and regression tree analyses were conducted using the rpart package in R to interrogate both the Road River Group and the global database and test whether the visually observed changes match the best split in an unbiased analysis.

For analyzing geochemical changes in the global dataset, a reweighted bootstrap method (36, 79) was used to account for differential spatial-temporal sampling density. This analysis weights each sample based on its spatial (modern coordinates) and temporal proximity to other samples in the dataset and then generates 1000 bootstrapped replicate mean values, with the likelihood that a sample is included in each bootstrap proportional to its weight. Given that most samples in our dataset are from well-constrained graptolitic shale successions, the additional layer of age uncertainty (36) was not implemented. Bootstrap mean proportions were plotted in 10-Ma bins along with 2 SD error bars. All codes used to produce geochemical figures and statistical analyses can be found at


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to the Na Cho Nyak Dun and Tetlit Gwitch’in communities for granting us permission to work at the Peel River locality. We thank Fireweed and Canadian Helicopters for safe transport, and P. Mamrol, V. Wala, T. Laakso, T. Boag, T. Allen, L. Van Drecht, and P. Sack for help in conducting fieldwork. We thank P. Linarde, C. Beck, T. Boag, H. Deres, S. Dobbs, P. Singh, M. McCormick, D. Turner, D. Mucciarone, and S. Ritzer for help in analyzing geochemical samples; B. Erdtmann, P. Wilde, and S. Butts for assistance in sampling museum collections; G. Somero for discussion; J. Sclafani for help in compiling geochemical data; and P. Sadler for providing the published graptolite diversity curve from (56). We thank C. Rasmussen, an anonymous reviewer, and B. Schoene for helping improve this manuscript during review. This is Yukon Geological Survey Contribution no. 052. Funding: This work was supported by an Ocean Sciences Research Fellowship from the Alfred P. Sloan Foundation, the affiliates of the Stanford Program on Deepwater Depositional Systems, the McGee-Levorsen fund at Stanford, and National Science Foundation (NSF) grant EAR-1922966 awarded to E.A.S. J.V.S. was funded by a postdoctoral fellowship from the Agouron Institute and NSF Tectonics Grant EAR-1624131. M.J.M. was funded by a Discovery Grant from the Natural Sciences and Engineering Research Council (Canada). Author contributions: E.A.S., J.V.S., and T.F. designed the study. M.J.M. analyzed graptolite biostratigraphy. E.A.S., J.V.S., T.F., M.J.M., R.G.S., and J.M. conducted primary fieldwork on the Peel River, Yukon. A.J.M., B.C.G., B.B., D.B.C., E.A.S., J.M., J.M.V., L.B., N.J.P., R.G.S., S.A.T., S.P.-T., T.N.B., and U.C.F. conducted geochemical analyses. A.D.R., A.J.M., A.L., D.K.L., E.A.S., J.V.S., and N.J.P. provided Paleozoic shale samples. E.A.S., R.G.S., U.C.F., and N.J.P. analyzed the data. E.A.S. wrote the paper with contributions from all authors. The authors who conducted geochemical analyses or provided samples are listed in alphabetical order. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article