Atmospheric oxygenation at the onset of Earth’s Great Oxidation forced enhanced marine anoxia

Treated as a biosignature in the search for life, free molecular oxygen is typically considered a prerequisite to sustain widespread eukaryotic life (Nursall, 1959). Nevertheless, despite oxygen comprising 21% of Earth’s contemporary atmosphere, this life-sustaining gas was conspicuously absent for half of our planet’s existence. In its broadest sense, the oxygenation of Earth’s atmosphere occurred via two key steps around the transition from the Archaean to the Palaeoproterozoic, commonly referred to as Great Oxidation Event (GOE) (e.g., Holland, 2002, 2006) and again at the end of the Neoproterozoic, known as Neoproterozoic Oxidation Event (NOE) (e.g., Och and Shields-Zhou, 2012). The oxygenation of the Earth most probably originated as a consequence of the evolution of oxygenic phototrophs and perturbations in the oxygen source-sink balance of the planet (Brocks et al., 1999; Catling et al., 2001; Holland, 2002; Kump and Barley, 2007). Interestingly, signs of marine oxygenation occur much earlier than atmospheric oxygenation with heavy stable isotope systems (e.g., Fe, Cr, Mo, U, Tl) and redox-sensitive trace element systematics disclosing small to medium scale marine oxygenation in so-called ‘oxygen oases’ up to several hundred million years prior to the GOE (Brocks et al., 1999; Rosing and Frei, 2004; Wille et al., 2007; Duan et al., 2010; Falcón et al., 2010; Planavsky et al., 2014; Kurzweil et al., 2015; Eroglu et al., 2017; Eickmann et al., 2018; Ossa Ossa et al., 2018; Albut et al., 2019; Ostrander et al., 2019a; Roué et al., 2021). A key indicator for atmospheric O₂ evolution is the transition from mass-independent fractionation of sulphur isotopes (MIF-S) to mass-dependent fractionation of sulphur isotopes (MDF-S), which occurs once pO₂ concentrations reach a threshold value of 10^-5 to 10^-6 of the present atmospheric levels (PAL) (Farquhar et al., 2000; Zahnle et al., 2006). This is currently best constrained to have occurred for the first time between ~2.45 Ga and 2.32 Ga in the global geological record (Philippot et al., 2018; Warke et al., 2020), but a deeper understanding of improved stratigraphic correlations and more absolute age constraints within sedimentary basins and between cratons are still needed. Two fundamental formations in this respect are the stratigraphically correlative Duitschland and Rooihoogte formations of the Transvaal Supergroup (South Africa), which both record some of the most prominent and abrupt MIF-S to MDF-S transitions in the global geological record (Guo et al., 2009; Luo et al., 2016; Poulton et al., 2021; Izon et al., 2022). Thus, these two formations have been studied intensively in order to understand atmospheric oxygen dynamics during the early stages of the GOE. However, the contemporary marine realm has received comparatively less attention during these two formations’ deposition. Hence, we measured a series of geochemical environmental tracers, including total sulphur (TS) and total organic matter (TOC) concentrations as well as trace element data, and stable isotope redox proxies (Mo, U, V, δ³⁴S, δ^98/95Mo_3134+0.25, Fe³⁺/Fe_tot) on four drill cores intersecting the Duitschland and Rooihoogte formations. Hereby, we aim to improve the current knowledge of the temporal redox evolution in the marine ecosystem at this critical junction in Earth’s history and understand the interaction between atmospheric and marine redox change during the onset of the GOE.

The Transvaal Supergroup, (South Africa), is preserved within two main structural basins – the Griqualand West Area (GWA) and the Transvaal Area (TA) – covering an area of about 500 000 km² (Beukes, 1987; Coetzee, 2001; Beukes et al., 2002). Representing remnants of a single marine basin that once covered the entire Kaapvaal Craton, the GWA and TA capture principally slope- and platform-type sedimentation, respectively. The TA features three unconformity-bound sedimentary packages (Figure 1) that start with the predominantly siliciclastic Wolkberg Group (Beukes, 1987; Sumner and Beukes, 2006), overlain by transgressive – regressive chemical sedimentation of the Chuniespoort Group (Beukes, 1984; Beukes, 1987; Sumner and Beukes, 2006) that, itself, yields to the largely siliciclastic Pretoria Group (Eriksson and Botha, 1988; Coetzee, 2001; Schröder et al., 2016). In the GWA, the Schmidtsdrift Subgroup and Ghaap Group represent respective equivalents of the Wolkberg and Chuniespoort groups of the TA (Beukes, 1984; Beukes, 1987). The Postmasburg Group GWA has no preserved correlative in the TA and consists of various lithologies, ranging from glacial diamictite, volcanic rocks and chemical precipitates in ascending order. It is unconformably deposited on the Ghaap Group GWA and chronostratigraphically appears before the Pretoria Group TA (Moore et al., 2012; Gumsley et al., 2017) (Figure 1).

The Duitschland and Rooihoogte formations mark the onset of the Pretoria Group within the TA, consisting of two lobes that are referred to as the Eastern Transvaal Area (ETA) and Western Transvaal Area (WTA) Transvaal areas that are separated by the Bushveld Igneous Complex (Figure 1). The Duitschland formation is geographically limited to the northernmost part of the eastern lobe. In contrast, the Rooihoogte Formation spans the remainder of the eastern and western lobes (Figure 1). Despite having different lithostratigraphic names, we consider the Duitschland and Rooihoogte formations to be different environmental expressions of the same sedimentary formation deposited at different distances to the alaeo-shoreline, based on arguments lined out in Havsteen et al. (2023). The strongest argument for stratigraphic equivalence between the Duitschland and Rooihoogte formations is their remarkably consistent internal stratigraphy, which is generally traceable throughout the entire Transvaal Basin. Specifically, both formations entail a basal conglomerate and diamictite, a mid-formational unconformity and a capping chert breccia that is locally karstified. Additionally, there are no discernable differences in the major and trace element compositions, U-Pb distribution patterns of detrital zircon, or radiogenic Sr-Nd-Hf isotope compositions of materials from the two formations (Coetzee, 2001; Schröder et al., 2016; Zeh et al., 2020; Havsteen et al., 2023). However, we note that other researchers regard each formation as a distinct sedimentary package, placing the Duitschland Formation depositionally before the Rooihoogte Formation (Moore et al., 2012; Bekker, 2015; Gumsley et al., 2017; Poulton et al., 2021; Senger et al., 2023). This view especially arises from not regarding the uppermost chert contact as a locally karstified surface, but rather as an erosional contact between the Duitschland and the Timeball Hill formations and a sharp conformable contact between the Rooihoogte and the Timeball Hill formations (Gumsley et al., 2017). Furthermore, an ongoing debate surrounds the interpretation of the maximum depositional age of the Duitschland Formation inferred from U-Pb dating of detrital zircons, as opposed to the well-constrained and undisputed depositional age at the transition between the uppermost Rooihoogte Formation and lowermost Timeball Hill Formation of 2 316 ± 7 Ma by Re-Os dating of black shales (Hannah et al., 2004). Zeh et al. (2020) originally interpreted the maximum depositional age of the Duitschland Formation to be 2 342 ± 18 Ma, based on the U-Pb date of the youngest, singular detrital zircon of their data set, and thus similar to the age of the Rooihoogte Formation. However, Senger et al. (2023) reinterpreted the maximum depositional age of the Duitschland Formation to be 2 427 ± 18 Ma based on the youngest reproducible zircon cluster from the data set of Zeh et al. (2020), and thus placing the Duitschland Formation to be significantly older than the Rooihoogte Formation.

Both formations are composed of argillite, arenite and interbedded conglomerate with subordinate chert and carbonate beds, sharply overlain by the pyritic black shales of the lower Timeball Hill Formation (Coetzee, 2001). The Duitschland and Rooihoogte formations can be divided into three intervals (Coetzee, 2001; Warke and Schröder, 2018; Havsteen et al., 2023):

• the lowerD/R interval, represented by a basal glacial diamictite and a chert pebble conglomerate (sometimes breccia) known as the Bevets conglomerate;

• the middleD/R interval consisting of two upward coarsening deltaic sequences terminated by the mid-Duitschland unconformity (MDU), and

• the upperD/R interval, comprising interbedded argillites and carbonates that overall represents a shallowing of the system.

Above lies the deep water pyrite-rich black shales of the lower Timeball Hill Formation (lower TBH Fm.), commonly with a chert layer marking the transition.

This study leverages samples from four scientific drill cores obtained during the CIMERA-Agouron GOE and Biomarker Drilling Project. All major and trace element analyses, detailed core descriptions, and sample positions are reported in Havsteen et al. (2023), while this work only provides a condensed summary of the core descriptions. The Rooihoogte Formation is captured in three cores, two obtained near Gopane (AGP-1 and AGP-2; 25°17.359S, 25°44.493E) and the remainder near Ngodwana (ANW; 25°36.239S, 30°37.002E). Except for the transition to the Timeball Hill Formation, AGP-1 captures almost the entire Rooihoogte Formation, while AGP-2 transects the upper Rooihoogte Formation and the lower Timeball Hill Formation. The ANW core transects the upper Rooihoogte Formation and lower Timeball Hill Formation, while the Duitschland Formation is represented within the ADL core, obtained near the Duitschland Farm KS95 type-area (24°17.268S, 29°07.491E) in the northeastern ETA. The ADL drill core intersects the upper part of the middle Duitschland and most of the upper interval of the Duitschland Formation. Fine-grained siliciclastic lithologies (shale–silt) were targeted for sampling, taking care to avoid veining and secondary mineralisation.

Depending on the samples Mo contents, between 100 and 300 mg of powdered sample were weighed into Savillex PFA beakers along with an appropriate amount of a ⁹⁷Mo-¹⁰⁰Mo double-spike tracer to achieve the desired 1:1 sample–spike ratio. A 2:5 mixture of 28 M HF and 14 M HNO₃ was added to the sample powders, and sample digestion and spike equilibration were achieved via closed vessel heating atop a hot plate at 100°C for three days with intermittent treatments in an ultrasonic bath. After overnight drying at 80°C, the fluoride complexes that formed during the initial digestion step were dissolved by a three-day reflux with 2 mL of 6 M HCl at 120°C, leaving a clear, residue-free digestate. Finally, after dissolution in 2 mL 3 M HCl, molybdenum was purified via the one-step anion exchange chromatographic protocol described in Willbold et al. (2016), using 7 mL columns filled with 2 mL of Eichrom AG1-X8 (100 to 200 mesh) resin.

Molybdenum isotope measurements were made using the ThermoFisher^© Scientific NeptunePlus multi-collector inductively coupled plasma mass spectrometer (MC-ICP-MS) housed within the Isotope Geochemistry Group at the University of Tuebingen, Germany. Purified Mo was introduced to the MC-ICP-MS as 25 ppb solutions in 0.3 M HNO₃ via a CETAC Aridus II™ desolvating nebuliser system. Static measurements comprising 90 cycles with 4.2 second cycle integration times were performed in low-resolution mode. Collectors were arranged so that potential isobaric interferences from Ru could be monitored via the ⁹⁹Ru ion beam, allowing hypothetical correction to the ⁹⁸Mo and ¹⁰⁰Mo signals. All samples, standards and procedural blanks were bracketed by measurements of the carrier solution 0.3 M HNO₃, allowing on-peak-zero background subtraction. Double spike deconvolution was performed assuming compliance with exponential mass fractionation law, effectively correcting for any mass-dependent fractionation induced through purification and analysis (Rudge et al., 2009).

All data are reported in δ-notation relative to the international reference material, NIST-SRM-3134 (NIST 3134; Equation 1), and expressed in per mille by multiplication with a factor of 1000.

To facilitate comparison to earlier studies that leveraged a Johnson Matthey (JM) Bern standard, however, a 0.25‰ offset was applied to compensate for the isotopic differences between the JM-Bern standard and NIST SRM3134 reference material (Nägler et al., 2014) as employed for example by Craig (1957) for O isotopes and, and also recently applied for Zn isotopes by Rosca et al. (2021).

Inclusive of unknowns, duplicates, reference materials and blanks, a total of 96 Mo isotope measurements were made over four analytical sessions. Throughout the study, the NIST-3134 and the JM-Bern standard solution yielded respective average δ^98/95Mo₃₁₃₄ values of 0.000 ± 0.044‰ (2SD, n = 58) and -0.277 ± 0.051‰ (2SDσ, n = 43), in accordance with published values (Greber et al., 2012; Goldberg et al., 2013; Kurzweil et al., 2016; Ossa Ossa et al., 2018) and, indeed, the Tuebingen lab long-term offset (Δ⁹⁸Mo_{(JM-NIST 3134)} = -0.274 ± 0.056‰ (2SD, n = 223). In Tables 1 to 4, we report both δ^98/95Mo₃₁₃₄ and δ^98/95Mo_3134+0.25 values, but we use δ^98/95Mo_3134+0.25 values in all figures and the text, as suggested by Nägler et al. (2014). Repeated processing of the reference material, JB-2 (Geological Survey of Japan), yielded an average δ^98/95Mo_3134+0.25 value of 0.279 ± 0.059‰ (2SD; n = 6, digests = 4), which is indistinguishable from the compiled reported average (0.298 ± 0.033‰ (2SD, n = 51) (Freymuth et al., 2015; Willbold et al., 2016; Zhao et al., 2016; Casalini et al., 2019; Chen et al., 2019; Villalobos-Orchard et al., 2020) and the long-term in-house average (0.284 ± 0.046‰, 2SD, n = 34). Total procedural Mo blanks were below 0.5 ng, with one outlier at 2.7 ng. However, when compared to δ^98/95Mo_3134+0.25 of duplicate digests of the same samples measured in different batches, the values were identical (Table 4), and thus no isotopic blank correction was performed.

Following the protocol outlined in Izon et al., (2022), major sulphur isotope data were generated by fluorination within the Geobiology Laboratory at the Massachusetts Institute of Technology, USA. Briefly, aliquots of powdered sample were boiled with acidified ethanoic chromous chloride, allowing their sulphide-sulphur inventories to be captured as zinc sulphide. Dropwise addition of silver sulphide then displaced zinc, leaving silver sulphide precipitates that were washed and captured. Thereafter, 3 mg aliquots of dry silver sulphide were converted to sulphur hexafluoride (SF₆) via overnight reaction at 300°C with fluorine gas. Finally, the analyte was cryogenically separated from residual F₂ and impurities before preparative-gas-chromatography, which ultimately yielded pure SF₆ for mass spectrometric analysis. Monitoring SF⁵⁺ ion beams at mass/charge ratios (m/z) of 127 and 129, ³⁴S/³²S isotope ratios of samples were expressed in conventional delta-notation (i.e., δ³⁴S) relative to the international reference standard, Vienna Canyon Diablo Troilite (VCDT).

Following the Pratt Method (Maxwell, 1968), the proportion of ferrous to ferric iron in selected samples from the AGP-1 drill core was determined via rapid acid digestion and subsequent colorimetric redox-titration. Depending on iron content, two doublets of sample powder of either 150 to 200 or 200 to 250 mg were digested for 12 min at 170°C in a mixture of 27.5 M HF (5 mL), 18.0 M H₂SO₄ (5 mL) and water (5 mL) in a teflon liner. The samples were then immediately transferred to an Erlenmeyer flask containing 50 mL of previously boiled and still warm milli-Q18.2 MΩ · cm H₂O, 30 mL of saturated H₃BO₃ solution and 6 mL of 14.6 M H₃PO₄ and diluted with a further 140 mL of warm 18.2 MΩ · cm H₂O, and titrated directly against a KMnO₄ solution.

Calibration of 0.1 M KMnO₄ against 0.1 M oxalic acid was performed daily and used for the titration. Each measurement session was initiated with two blank measurements, followed by four rock reference materials (BIR-1a and OU-6) before sample measurements were instigated. In the rare instance where the FeO content of the duplicate digestions deviated by more than 0.5 wt.%, an additional measurement was performed. Processed alongside our unknowns, replicate analyses of USGS BIR-1a yielded a mean FeO content of 8.75 ± 0.61 wt.% (1SD; n = 16), in agreement with published estimates (Saikkonen and Rautiainen, 1993; Schuessler et al., 2008; Babechuk et al., 2019); Likewise, OU-6 (IAG) was processed 24 times, returning a mean FeO content of 1.87 ± 0.23 (1SD), inseparable from its assigned value (1.65 ± 0.019 wt.%; Kane (2004)).

Total carbon, nitrogen, and sulphur concentrations were determined by flash combustion using the Elementar VARIO EL III analyzer within the Soil Science and Geomorphology Group, University of Tuebingen, Germany. Approximately 40 mg of sample was wrapped into tin foil capsules and analysed using oxidative heat combustion at 1 150°C with tungsten trioxide as a catalyst. Triplicates of the organic-rich reference material, OAS IVA 33802150, were measured every 30 samples to test accuracy and precision of the method. The long-term (n >2000) internal relative standard deviations (1rsd) are 1.3%, 2.5% and 2.4% for total carbon, nitrogen and sulphur determinations, respectively. These numbers are adapted as an estimate of method precision.

Carbonate-carbon was removed from 20 mg aliquots of powdered samples via multiple treatments with 5% HCl. The residues were washed with 18.2 MΩ · cm H₂O and dried at 60 °C for subsequent mass spectrometric analyses. Dependent on their carbon contents, 0.05 to 1.3 mg of decarbonated residue was wrapped in Sn-foil and loaded into the autosampler of a Carlo-Erba (NC 2500) elemental analyzer interfaced with a Finnigan Delta Plus XL isotope ratio mass spectrometer. Samples were combusted at 1 050°C in an oxygen stream and the combustible component was swept through the oxidation–reduction furnace, and the resultant CO₂ was ultimately cleaned by gas chromatography before admittance into the isotope ratio mass spectrometer (IRMS). The sample CO₂ is measured and evaluated relative to an internal laboratory tank gas standard calibrated against internal and international standards (Acetanilide). All δ¹³C_org values are given in permille relative to Vienna Pee Dee Belemnite (VPDB), and the external precision calculated over 10 to 15 standards is typically in the range of 0.1‰.

The authigenic components of the samples’ redox-sensitive trace elements (RSTE) molybdenum, uranium and vanadium were determined by correcting their detrital background with two different methods. The first method (model 1; Algeo and Tribovillard, 2009) corrects for any detrital background components by normalising to aluminium (Al), under the assumption that all Al is detrital and the siliciclastic portion of the sediments resembles the composition of the Post Archaean Average Shale (PAAS) (Taylor and McLennan, 1985). By using modern restricted basins as analogues, covariations in U and Mo enrichment factors were subsequently utilised to differentiate marine redox environments in the geological record (Figure 2). Enrichment factors (EF) for U and Mo were calculated according to equation 3.

Where Al and X are the weight concentration of Al in wt.% and element ‘X’ in μg/g in the sample and the selected normalisation material. Most samples from all drill cores plot within the suboxic field, with a few samples extending towards an anoxic or particulate shuttling environment (Figure 2). It is important to note that some samples fall below 1 in Mo and/or U enrichment factors and can thus not be distinguished from detrital background levels.

The second method (model 2; Bennett and Canfield, 2020) used to determine the authigenic enrichment of RSTE is avoiding normalisation to a crustal average (e.g., PAAS or upper continental crust (UCC)), which may not be a good representative for the particular terrigenous input source and speed of the Transvaal sediments. Instead, this model uses receiver operating characteristic curve analysis to quantitively determine enrichment factor thresholds, allowing differentiation between various marine settings with equal weight on minimising false positives and negatives (Bennett and Canfield, 2020). Enrichments are calculated using equation 4.

Where X_E is the enrichment of element X (V, Mo, Re, U) and Al and X are the concentration of Al and element X in the sample. The threshold values between different redox settings are defined in Bennett and Canfield (2020). Most samples from all cores plot in the oxic field regardless of tracer system used (V/Al versus U/Al or V/Al versus Mo/Al; Figure 3). However, some samples plot in the ‘beneath perennial-oxygen minimum zone’ (beneath P-OMZ), indicative of oxic bottom waters. This is especially true for samples from AGP-1 and ANW drill cores.

The entire Duitschland and Rooihoogte succession extending into the lowermost Timeball Hill Formation is only captured by combination of drill cores AGP-1 and -2 (Figure 4A). The middleD/R and upperD/R intervals display decreasing trends in δ³⁴S and δ¹³C from approximately +15‰ to -20‰ and approximately -30‰ to -35‰, respectively (AGP-1, 151 to 48 m; Figure 4A, columns 2 and 6). In concert with this isotopic decrease, sulphur contents rise from levels typically below the limit of detection (LOD) to concentrations as high as 0.6 wt.% in the upperD/R interval. Likewise, TOC shows stable values with an average of 0.59 wt.% in the lowerD/R interval (AGP-1, 151 to 105 m) that increases to as much as 5.9 wt.% in the upperD/R interval (AGP-1, 92 to 70 m; Figure 4A, column 5). Over the same stratigraphic section, δ^98/95Mo_3134+0.25 values also increase from - 0.657‰ to as much as 1.072‰, while Fe/Mn ratios overall gradually increase throughout the upperD/R interval (AGP-1, 92 to 47 m; Figure 4A, columns 1 and 3). The transition from the upperD/R interval into the lower TBH Formation (AGP-2, 102 to 82 m) records low and stable δ³⁴S and δ¹³C values around -23‰ and -37‰, respectively (Figure 4A, columns 2 and 6). Over the same stratigraphic section, TS and TOC values stabilise, with average concentrations of 0.370 wt.% and 1.507 wt.%, respectively, while Fe/Mn ratios steadily increase (AGP-2, 95 to 83 m) before stabilising around 77 m and upwards (Figure 4A, columns 3, 4 and 5). From 108 to 83 m in AGP-2, δ^98/95Mo_3134+0.25 values display a generally negative trend from 0.788 to approximately -0.057‰ before embarking on an oscillatory trend between 83 and 59 m.

The trends observed in cores ANW (Figure 4B) and ADL (Figure 4C) show good agreement with those derived from the complete succession captured by the AGP cores. This is most evident in ANW, which records the transition from the upperD/R interval into the lower TBH Formation. Like the AGP drill cores, the ANW captures a decrease in δ³⁴S and δ¹³C values that manifest as respective approximately 30‰ and 35‰ negative excursions within the upperD/R interval (ANW, 49 to 37 m) before stabilising in the lower TBH Formation (Figure 4B, columns 2 and 6). This isotopic decrease corresponds with an increase in TS and TOC, with values rising abruptly at 41 m in the upperD/R interval, accompanied by a steady increase in Fe/Mn ratios and increasing δ^98/95Mo_3134+0.25 values from -0.727 to 0.672‰ at the stratigraphic level from approximately 43 to 38 (Figure 4B, columns 1 and 3). In the upper part of the upperD/R interval (ANW, 38 to 32 m), TS, Fe/Mn ratios and δ^98/95Mo_3134+0.25 values steadily decrease to 0.156 wt.%, 20 and 0.089‰, respectively, while δ³⁴S, δ¹³C and TOC stay relatively constant with means of -34.9‰, -23.3‰ and 5.3 wt.% (Figure 4B, columns 1 through 6). In the lower TBH Formation, δ^98/95Mo_3134+0.25 values gradually become more positive, peaking at 0.740‰ at 30 m core depth before decreasing to -0.352‰ in the topmost sample (Figure 4B, column 1). δ³⁴S and δ¹³C values are low and stable throughout the lower TBH Formation with average values of -22.2‰ and -35.4‰, respectively, while Fe/Mn and TS increase from 35 and 0.434 wt.% to 67 and 0.677 wt.%, (Figure 4B, columns 2, 3, 4 and 6). TOC shows an oscillatory trajectory, starting at 5 wt.% in the upperD/R interval and lowermost TBH Formation (ANW, 36 to 33 m), whereafter values drop to 1.5 wt.% (31 m) and shortly after peaking at 7.7 wt.% at 29 m core depth, (Figure 4B, column 5).

The ADL drill core almost solely records the upperD/R interval till below the Timeball Hill Formation. δ³⁴S values are relatively stable between 0 and 10‰ with a slightly increasing trend in a parabolic fashion between 270 to 109 m core depth, while δ¹³C values oscillate between -31 and -23‰ over the same interval (Figure 4C, columns 2 and 6). Both TS and TOC show low and stable values in the upperD/R interval (ADL, 267 to 78 m) with average concentrations of 0.1708 (LOD values omitted, see Table 4) and 0.085 wt.%, respectively, while Fe/Mn have an average of 85 and show no particular trend, (Figure 4C, column 3, 4, 5 and 6). The δ^98/95Mo_3134+0.25 compositions gradually decrease from 0.393 to -0.088‰ between 270 to 109 m, with the lowest value (-0.794 ‰) at 80 m core depth (Figure 4C, column 1).

An important feature between the drill cores is the trends observed in Fe/Mn ratios, which increase with stratigraphic height in the AGP-1, AGP-2 and ANW drill cores but with subtle differences in the relative increase (Figure 4A and B, column 3). The Fe/Mn ratios range from 30 to 200 (AGP-1), 25 to 230 (AGP-2) and 20 to 60 (ANW), while no particular trend is observed in the ADL drill, which ranges in values from 80 to 240 (ADL), (Figure 4C, column 3). Additionally, the AGP-1 drill core has two outliers of 300 (AGP-1 48.14 m) and 578 (AGP-1 47.66 m) (Figure 4A, column 3). It is important to note that the same trend is observed in the drill cores throughout the basin in the upperD/R interval and that the increasing Mn/Fe ratios seem to reflect an increasing iron content with stratigraphic height while Mn displays consistently low and stable concentrations. The exceptions are the lower part of the upperD/R interval in the ANW drill core (45 to 43 m) and the upper part of the middleD/R and lowermost upperD/R intervals in the AGP-1 drill core (105 to 86 m) (Figure 4A and B, column 3). Concerning the interval in the ANW drill core (back bulge), MnO concentrations display a mean of 0.22 wt.% (n = 4), which is significantly higher than background levels (PAAS = 0.06 wt.%; UCC = 0.08 wt.%), (Figure 4B, column 3). For the AGP-1 drill core, the lower part of the middleD/R interval (151 to 122 m) MnO concentrations display a low average of 0.08 (n = 4), which is similar to the UCC and PAAS, whereas the upper part of the middleD/R interval (105 to 92 m) display MnO content above average crustal concentrations with a mean MnO content of 0. 0.15 wt.% (n = 5), (Figure 4A, column 3).

Starting with the most distal core, ANW, the transition to the Timeball Hill Formation is not marked by a change in grain size; rather, the transition is recorded as a change from grey to black shale. This might suggest that the observed redox change is not predominantly sea-level controlled, although sea-level shifts in mudstones can, at times, manifest as solely compositional changes without significant alteration in grain size (Dinelli et al., 2007). Similarly, the ANW core captures a change in redox at 43 to 40 m core depth, which is also not accompanied by a change in grain size (Figure 4B). In the AGP-1 (middle - and upperD/R intervals), a coarsening-upwards regressive system tract (RST) is succeeded by several minor TSS. However, the change in redox parameters registered in the lower to middle part of the upperD/R interval (AGP-1, 75 to 100 m) post-dates the MDU by several meters of strata (Figure 4A). Thus, the observed changes in geochemical redox proxies seem to predominantly reflect redox changes in the basin that are not controlled by relative sea level change.

Given that trace element cycling was likely significantly different during the Palaeoproterozoic (Da Silva and Williams, 2001; Anbar, 2008; Robbins et al., 2016), direct transfer of proxies from modern ocean settings should be applied carefully to the ancient sedimentary rock record. Additionally, knowing that post-depositional diagenetic effects, microbial activity and hydrothermal activity all have the potential to perturb the elemental and isotopic budgets of ancient records, these processes must also be considered carefully (Morford et al., 2005; Johnson et al., 2013). To avoid these pitfalls, herein, we cross-calibrated and validated our RSTE-based inferences with additional well-established environmental and redox-sensitive proxies (δ^98/95Mo_3134+0.25, Fe/Mn, δ³⁴S, TS, TOC, δ¹³C_org). While we concede that these proxies are also affected by post-depositional effects, combining several independent redox proxies with differing sensitivities to alteration helps to identify compromised data and to fortify our inferences.

Both RSTE models are in broad agreement, although subtle differences are apparent. Principally suboxic conditions were identified by model 1 (Figure 2) with transient shifts toward anoxic and oxic deposition, as well as intervals influenced by Mn-shuttling. Model 2 (Figure 3), however, fingerprints a predominantly oxygenated depositional setting with the episodic development of euxinic conditions or intervals reminiscent of those that resemble deposition beneath P-OMZ. Comparatively, model 1 exposes relatively few instances of fully oxygenated conditions (e.g., ANW, 45 to 43 m; AGP-1, oscillatory above 70 m). Against a modern oxygenated environment, deposition ‘beneath P-OMZ’ (model 2) identifies a setting with oxygenated bottom waters. In a Palaeoproterozoic context, with lower elemental seawater inventories, we adopt a more conservative ‘suboxic’ category for these samples. Some sections classified as ‘oxic’ in model 2 are considered ‘particulate influenced’ in model 1. In reality, both categories indicate significant amounts of O₂ in the water column, with particulate shuttling inferring a stratified water column where Mo is shuttled across the chemocline in association with soluble Fe-Mn (oxyhydrogen) oxides (Algeo and Tribovillard, 2009; Smrzka et al., 2019; Kurzweil et al., 2021; Kurzweil et al., 2022). Thus, areas identified as those influenced by ‘shuttle processes’ in model 1 and ‘oxic’ in model 2 both point to the presence of significant amounts of oxygen in the water column.

The most significant discrepancies in Palaeoenvironmental reconstruction are between the layers classified as ‘particulate shuttle influenced’ in model 1 and ‘euxinic’ in model 2 (i.e., ANW, 29.77 to 30.55 m). Here, δ^98/95Mo_3134+0.25 values evolve from -0.335 to 0.741‰, which is inconsistent with particulate shuttling because of the lack of the significant negative fractionation that accompanies this process between MoO₄^2- in seawater and precipitated Fe-Mn-(hydro)oxides (Barling and Anbar, 2004; Wasylenki et al., 2008). A similar conclusion is reached concerning the discrepancy seen at ~41 to 44 m (ANW), where model 1 indicates particulate shuttling while model 2 identifies deposition within a suboxic/anoxic oxygen minimum zone (OMZ). In this interval, δ^98/95Mo_3134+0.25 values increase from -0.728 to 0.754‰, again providing strong evidence against the operation of Mn-shuttling. Instead, the positive shifts in δ^98/95Mo_3134+0.25 compositions rather fit with an environmental transition from oxic to suboxic/anoxic conditions as indicated by model 2. In sediments deposited in anoxic to mildly oxygenated bottom waters (O₂< 10 μM), capable of reducing both Mn-oxides and sulphate, there are two main drivers controlling the δ^98/95Mo_3134+0.25 composition of the sediments. Firstly, the relative amount of Fe versus Mn drawdown with (oxyhydr)oxides and secondly, the aqueous H₂S concentration of the sediment pore and bottom waters (Goldberg et al., 2009; Poulson Brucker et al., 2009). Under ferruginous conditions, a larger proportion of Mo adsorbs onto Fe- rather than Mn-(oxyhydr)oxides, reducing the isotopic offset between Mo in seawater and the authigenic particles in the sediments. That is due to the smaller isotopic fractionation between the solid and solution (Δ⁹⁸Mo) for Fe-(oxyhydr)oxides and seawater (e.g., Δ⁹⁸Mo_{solution-ferrihydrite} = ~1.11 ± 0.15‰; Goldberg et al. (2009)), relative to Mn-(oxyhydr)oxides and seawater of Δ⁹⁸Mo_{solution–MnOx} = ~2.4 to 2.9 ± 0.1‰ (Barling and Anbar, 2004; Wasylenki et al., 2008). Similarly, in sediments deposited under high H₂S porewater conditions (but not sulfidic; H₂S > 11 μM), the Mo isotopic offset between seawater molybdate (MoO₄^2-) and sediment tetrathiomolybdate (MoS₄^2-) is controlled by incomplete transformation of intermediate thiomolybdate species (e.g., MoO₃S^2-, MoO₂S₂^2-, MoOS₃^2-). Accumulatively viewed incomplete thiolation induces a moderate offset from seawater of approximately -0.6 to -1.00‰ in δ^98/95Mo_3134+0.25 (Kendall et al., 2017), although the fractionation associated with a specific thiolation step may be significant (Nägler et al., 2011; Kerl et al., 2017). The development of ferruginous water conditions and/or H₂S-rich pore or bottom waters could hence independently or collectively, lower the isotopic offset between authigenic Mo and its isotopically heavier seawater precursor. Evidence for the development of such conditions at ~41 to 44 and ~29 to 30 m core depth in the ANW is indicated by increasing Fe/Mn ratios and decreasing δ³⁴S compositions coupled with high TS concentrations, respectively. Consequently, the observed shift in δ^98/95Mo_3134+0.25 is likely a combination of these environmental conditions. Furthermore, genuinely euxinic conditions (anoxic and sulfidic) probably only developed for short periods, as indicated by the intervals with high TS concentrations and increasing δ^98/95Mo_3134+0.25 compositions (e.g., ANW (40 to 43 m)). Assimilated, the water column within the Transvaal Basin during the deposition of the upper Duitschland and Rooihoogte formations and lower Timeball Hill Formation was, to a large extent, oxygen-poor (sub- to anoxic). and Fe-rich (Mn-poor). This interpretation is reinforced by literature δ^98/95Mo_3134+0.25 values and Fe_HR/Fe_T data from core EBA-2 that reveals an oscillating and largely anoxic background with isolated instances of euxinia (Asael et al., 2018).

Within the lower part of the middleD/R interval (122 to 151 m), the relatively low δ^98/95Mo_3134+0.25 values seen in core AGP-1 coupled with moderate Mo enrichments and MnO contents that approximate those seen in UCC and PAAS, fingerprint particulate shuttling where reductive Mn dissolution supplies Mo oxygen-poor bottom waters where it can be fixed within the sediments. In contrast, the lower part of the upper Rooihoogte Formation in cores AGP-1 (86 to 105 m) ANW (42 to 45 m) features up to two-fold MnO enrichments when compared to PAAS, implying the preservation of Mn oxides to at least some extent, and thus molecular oxygen within the water column (Figure 5). This inference is supported by low δ^98/95Mo_3134+0.25 values (-0.738 to -0.123 ‰), generally non-detectable sulphur contents and low TOC (approximately 0.4 wt.%; Figure 5, Table. 1).

In summary, the marine oxygenation state likely evolved from a stratified ocean in the lower reaches of the middleD/R interval into a fully oxygenated ocean in the upper part of the interval and the overlying lowermost segment of the upperD/R interval. In the middle and upper reaches of the upperD/R interval, however, the marine oxygenation state transitioned into a suboxic to anoxic state, with the episodic development of euxinia. This marine oxygenation state, predominantly in the suboxic to anoxic realm, stabilised and continued during the deposition of the lower Timeball Hill Formation.

Comparison of the CIMERA-Agouron GOE drill cores (Havsteen et al., 2023) with several drill cores that capture the loss of MIF-S (EBA-2, EBA-4 and KEA-4 drill cores) within the Rooihoogte Formation (Luo et al., 2016; Izon et al., 2022) and outcrop samples from the Duitschland Formation (Guo et al., 2009), suggest that the shift towards less oxygenated marine conditions is causally linked to atmospheric oxygenation. Importantly, this argument is valid regardless of whether our preferred correlation of the Duitschland and Rooihoogte formations is adapted or not. The pronounced shift toward negative δ³⁴S values seen in the Carletonville area (EBA-2, EBA-4 and KEA drill cores) in the immediate aftermath of the initial demise of MIF-S provides a convenient chemostratigraphic marker to examine our data.

While it is important to note that the upperD/R interval captures a basinal deepening (Coetzee, 2001; Havsteen et al., 2023), it is clear that the redox shift post-dates the mid-formational unconformity several meters of strata above the sedimentological transition from the middleD/R to the upperD/R interval. Moreover, sea level high stands are usually associated with negative δ^98/95Mo_3134+0.25 values as such episodes typically expose more sediments in the basin to H₂S poor conditions (Ostrander et al., 2019b). Thus, this rather counterintuitive trend of marine deoxygenation in response to atmospheric oxygenation, is best explained by geobiological feedbacks. A slightly oxygenated atmosphere, as expressed by the loss of MIF-S, drives oxidative weathering on the continents and elevates nutrient input to the oceans, resulting in increased primary productivity and increased oxygen demand within the deep waters (Figure 6). The conceptualisation of this idea has been suggested by Canfield (1998), but is geochemically rarely demonstrated as clearly and rapidly after the initial loss of MIF-S, as documented in this study. Especially the initiation of microbial sulphate reduction, as illustrated by the low δ³⁴S of approximately -25‰, highlights the connection between the build-up of sulphate levels (expressed as an increase in TS) in the marine realm and terrestrial pyrite oxidation. The negative δ³⁴S values (approximately -25‰) succeeding the loss of MIF-S strengthens this interpretation of increasing marine sulphate concentrations far above 200 μM, since marine sulphate reducers usually produce a low isotopic offset of less than 5‰ when sulphate concentrations are less than 200 μM (Habicht and Canfield, 2001). This threshold value for sulphate concentrations is much higher than the inferred sulphate concentration of approximately 2.5 μM in the Archaean ocean (Crowe et al., 2014), and thus only a minute offset between sulphate and sulphide is expected before a substantial marine sulphate pool is established. Higher marine sulphate concentrations would further dampen atmospheric methane fluxes to the atmosphere, since the accumulation of a marine sulphate pool perturbs the balance between sulphate reducers and methanogens, leaving less organic matter for methanogenesis (Canfield, 1993). Furthermore, an ocean rich in sulphate would increase anaerobic methane oxidation, allowing less methane to escape the sediments and reach the atmosphere, thereby perturbing the atmospheric source-sink balance in favour of O₂ accumulation (Catling et al., 2007). Likewise, the observed increase in TOC and decrease in δ¹³C values imply enhanced photosynthetic activity and carbon burial, leading to a net O₂ gain to the atmosphere. The integrated framework of presented data thus exposes a geobiological feedback-driven causality between the earliest atmospheric oxygenation and marine deoxygenation, controlled by the interplay between O₂ production by cyanobacteria in the shallow marine realm and perturbations of the sulphur cycle.

This study integrates several environmental and redox-sensitive proxies (δ^98/95Mo_3134+0.25, δ³⁴S, δ¹³C_org, TS, TOC, Fe/Mn and RSTE) from four drill cores, intersecting the correlative Duitschland and Rooihoogte formations, both recording the loss of MIF-S, in order to understand the temporal marine redox evolution in response to the earliest atmospheric oxygenation.

• Two independent RSTE models suggest an oxic-suboxic marine environment preceding the loss of MIF-S, which transitions into suboxic/anoxic conditions in the immediate aftermath of atmospheric oxygenation. This conclusion is further supported by increasing trends in δ^98/95Mo_3134+0.25, TS, and TOC and decreasing trends in δ³⁴S, δ¹³C_org and Fe/Mn. The initiation of oxidative continental weathering increased nutrient delivery to the oceans, triggered cyanobacteria blooms in the surface oceans and ultimately marine deoxygenation. This is expressed biogeochemically in the decrease of δ³⁴S and δ¹³C_org values in concert with elevated TS and TOC contents, signalling non-limiting microbial sulphate reduction and the ingrowth of the seawater sulphate reservoir. The positive shift in δ^98/95Mo_3134+0.25 values can thus be linked to decreased fractionation from ambient seawater driven by increased sulphur and ferrous iron in the bottom waters.

• Combined Fe/Mn, MnO, and δ^98/95Mo_3134+0.25 systematics disclose oxygenated shallow oceans predating the loss of MIF-S signals. However, on a basin-wide scale, the Eh potential of the water column was high enough to oxidise Mo, S and Fe but not high enough to significantly oxidise Mn in most cases, except in the uppermost part of the middleD/R and lowermost upperD/R intervals.

• Based on the combination of multiple geochemical redox-sensitive proxies, the shallow marine realm likely evolved from a stratified state into oxic-suboxic conditions before the loss of MIF-S, whereafter, the marine redox state transitioned to anoxic-suboxic conditions. Our data implies a geobiological feedback-driven relationship linking early atmospheric oxygenation to a decline in the marine oxygenation state. Yet, it remains to be determined whether this trend is constrained to the Transvaal Basin or represents a global phenomenon.

Sampling assistance from Lucile Roué and laboratory support from Elmar Reitter and Bernd Steinhilber are gratefully recognised. The Agouron Institute funded the drilling campaign that yielded the material analysed herein, while the University of Johannesburg provided invaluable logistical support and access to the cores. Financial support, thereafter, was provided by the German Research Foundation (DFG) grant SCHO1071/12-1 to Ronny Schoenberg as part of the DFG priority program SPP-1833 ‘Building a Habitable Earth’. Gareth Izon gratefully recognises his MISTI Seed Award as well as long-standing financial support from the Simons Collaboration (#290361FY18) and the National Science Foundation (#2148916) via respective awards to Roger Summons and Edward Boyle at MIT. Gareth is indebted to Shadrack Tshivhiahuvhi and Lucig Khoza, who aided his initial sampling campaign. We thank two anonymous reviewers for their careful and constructive feedback on this manuscript, which significantly improved its quality.

For a general overview of scatter plots illustrating the relationship between Mo isotope systematics ^98/95Mo_NIST+0.25 and selected isotopic and elemental proxies within the lower Pretoria Group refer to Figure S1. (A supplementary data file is archived in the South African Jounal of Geology repository (https://doi.org/10.25131/sajg.127.0002.sup-mat)).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Editorial handling: A.J.B. Smith.