Atmosphere-ocean oxygen and productivity dynamics during early animal radiations.

The proliferation of large, motile animals 540 to 520 Ma has been linked to both rising and declining O2 levels on Earth. To explore this conundrum, we reconstruct the global extent of seafloor oxygenation at approximately submillion-year resolution based on uranium isotope compositions of 187 marine carbonates samples from China, Siberia, and Morocco, and simulate O2 levels in the atmosphere and surface oceans using a mass balance model constrained by carbon, sulfur, and strontium isotopes in the same sedimentary successions. Our results point to a dynamically viable and highly variable state of atmosphere-ocean oxygenation with 2 massive expansions of seafloor anoxia in the aftermath of a prolonged interval of declining atmospheric pO2 levels. Although animals began diversifying beforehand, there were relatively few new appearances during these dramatic fluctuations in seafloor oxygenation. When O2 levels again rose, it occurred in concert with predicted high rates of photosynthetic production, both of which may have fueled more energy to predators and their armored prey in the evolving marine ecosystem.

The “Cambrian explosion” refers to a relatively short time interval, ca. 540 to 520 Ma, during which nearly every major animal group appeared in the fossil record, mineralized skeletons became widespread, and some large, motile animals began to plow through the sediment surface (1–4). The Cambrian explosion is one of a series of radiations that began in the late Ediacaran and continued through the early Paleozoic, during a time of significant fluctuations in oceanic redox conditions and nutrient supply (4–6). How these fluctuations relate to animal radiations is not well understood (4, 7–11). One hypothesis is that rising O2 availability in the marine realm allowed the evolution of carnivory, which in turn drove predator–prey “arms races” and concomitant diversification (3, 8). This radiation could result from either higher molar concentrations of O2 in the surface waters or at the seafloor, which would, in turn, expand oxygenated areas for the benthic animal fauna. With higher molar concentrations of O2 and/or greater habitats for benthic animals, oxygen could promote animal-driven mixing of sediments. However, as bioturbation ventilates the sediments, less phosphate is released from the sediments back into the water column resulting in decreased nutrients for primary producers and, thus, decreased atmospheric O2 production (7). In addition to O2 availability, food availability for prey is another factor that can drive predator–prey evolution toward greater diversity (12), but periods with high nutrient levels also can drive expansive ocean anoxia and limit the animal habitable zone. Therefore, we can envision animal proliferation to be correlated to either rises in surface O2 availability, increasing productivity, or declining seafloor anoxia. Alternatively, if food and O2 availability in the marine realm played only a subordinate role in the radiation, we may discover that radiation was not driven by environmental change and could simply result from new biological innovations and styles of interactions (13).

Attempts to test the food and oxygenation hypotheses have given ambiguous results (7, 9, 14, 15), in part due to insufficient spatial and/or temporal resolution of the available paleo-O2 and paleoproductivity proxies. To fill this gap, we compiled a record of global seafloor oxygenation constrained by high-resolution uranium isotope data in marine carbonates from the Terreneuvian Series of the Cambrian (∼539 to 521 Ma). Furthermore, we model an atmospheric pO2 curve from carbon, sulfur, and strontium isotope records in the same sedimentary sequences to constrain variations in O2 availability in the surface oceans through time. We propose to use these records in combination to indirectly constrain O2 consumption within the ocean interior, which, to a first order, is a function of the size of the global marine nutrient pool that controls new marine biomass production. With global estimates on how net productivity and O2 availability in the early Cambrian ocean varied through time, we explore potential causal relationships between the availability of food, O2, and animal radiations by comparing geochemical and paleobiological records from the same sedimentary successions.

The Uranium Isotope Proxy for Global Seafloor Oxygenation

The uranium isotope composition (δ238U, the permil deviation of the 238U/235U ratio relative to the CRM 145 standard) of seawater is a promising global oxygenation proxy because of its sensitivity to ocean oxygenation at a globally integrated scale. This attribute results from the uniform δ238U distribution in the well-mixed oceans and a strong U isotope fractionation that occurs when sediments are buried in anoxic environments (16–18). Model calculations suggest that oceans still would have been homogeneous in U isotope composition during the Cambrian, despite potentially greater U burial in anoxic settings that were generally more common at the time (19). The δ238U proxy has been applied to better constrain oceanic anoxic events (OAEs) in the past and reliably recorded negative δ238U excursions at times when anoxic water masses expanded in the oceans (19–23). The OAEs studied to date are the Cambrian SPICE event (−0.2‰), Late Ordovician extinction (−0.3‰), Frasnian–Famennian extinction (−0.3‰), End Permian extinction (−0.5‰), End Triassic extinction (−0.7‰), and the Cretaceous OAE2 (−0.1‰) (19, 24–28).

The translation of seawater δ238U into the carbonate rock record has been studied both experimentally and in modern environments (29–32). Uranium exists in the form of uranyl carbonate anions in seawater that are incorporated into primary biogenic carbonate minerals with δ238U indistinguishable from seawater (−0.39 ± 0.04‰, 2 SD) (16). As carbonate cements form inside the sediments, they tend to adopt a slightly higher and variable δ238U value. Even in the very oxic oceans today, many of the pore waters in sediment columns are anoxic. In modern Bahamas bank sediments, uranium reduction occurs near the sediment–water interface allowing for nonquantitative U removal that leaves the bulk sediments on average +0.24 ± 0.14‰ (1 SD) heavier than overlying seawater (29–31). The scatter in this offset is not well understood but should be linked to the ability of dissolved U to escape from the sediments. The precursor sediment in the Terreneuvian oceans included both aragonite and calcite phases (33). Today, metastable aragonite and high-Mg calcite are preserved even at several hundred meters of burial depths (34), so U release from these phases would correspondingly occur at these depths. The dissolution process itself causes no U isotope fractionation (32), but the subsequent U(IV) incorporation into secondary and more stable calcite phases does cause U isotope fractionation (35), as may calcite recrystallization under pore water conditions distinct from seawater (36). In all cases, the U isotope fractionation is only expressed in the remaining sediments if U escapes from the system (e.g., in a “fluids-buffered” system as opposed to “rock-buffered” with respect to U). Apart from oxidative weathering in the meteoric vadose zone, where dissolved U(VI) exists, diagenetic alteration at depth appears to induce no systematic changes in δ238U, presumably because of the low solubility of U(IV) phases that form already at very shallow burial depths (29–31). In this study, we use petrography and various diagenetic proxies (including Mn/Sr, Sr/Ca, Mg/Ca, and δ18O) that also are affected by the fluid vs. rock-buffered style of diagenesis. This enabled us to identify samples with the smaller diagenetic δ238U offset from contemporaneous seawater. Furthermore, we also adopt a high-stratigraphic resolution sampling strategy (∼10 per Ma) recommended to correct for stochastic effects in the reconstruction of secular changes in seawater δ238U from bulk carbonate sedimentary rocks (30).

The Carbon, Sulfur, and Strontium Isotope Proxy for Atmospheric Oxygenation

The combination of carbon, sulfur, and strontium isotope records from shallow carbonate sediments potentially record global changes in the chemical composition of the oceanic reservoir, and if that is the case, these proxies allow us to reconstruct variations in the atmospheric O2 budget during the Terreneuvian (∼539 to 521 Ma) using a simple mass balance approach for these marine element cycles (37, 38). At this time, 2 of the largest positive carbon isotope excursions of the Phanerozoic are recorded in marine carbonate sediments worldwide (33, 39, 40). We recognize that part of the geochemical signals from shallow carbonates reflects local and/or diagenetic processes, and a major task is to further develop ways to tease apart these local signals from the global. However, in some cases where the inferred seawater isotope trajectories covary over approximately million-year timescales consistent with the respective oceanic residence timescales, we ask what we expect to see, if these signals record global change. The principal atmospheric O2 sources, organic carbon and pyrite burial in marine sediments, are constrained from the seawater δ13C and δ34S records, whereas atmospheric O2 sinks (mainly reducing gasses emanating from Earth’s interior) and oxidative weathering on land are assumed proportional to changes in weathering fluxes on land, constrained by the 87Sr/86Sr record (41).

The changes in atmospheric O2 sourced from the marine carbon cycle are calculated from organic carbon burial flux (JORG) in excess of oxidative carbon weathering:Furthermore, we assume that total oxidative carbon weathering flux is proportional to continental weathering of carbon bearing rocks [JW(C)] and the organic proportion of carbon undergoing weathering (RORG). A similar equation describes the O2 equivalent released from the marine sulfur cycle (see details in ).

If the carbonate δ13C record is related to the global dissolved inorganic carbon pool, then the global organic carbon burial flux is calculated from the shallow marine water δ13C record (abbreviated δC) and its derivative (dδC/dt) as a function of average C isotope offset between seawater and buried organic matter (ΔC), δ13C of weathering input [δW(C)], and the oceanic C inventory (MC) (37):A similar relationship for the marine sulfur cycle (Eq. 14) constrains total pyrite burial from the seawater δ34S record. The total change in atmospheric pO2 is modeled by summation of the integrated O2 flux from the marine cycles of carbon (Eq. ) and sulfur (Eqs. 15 and 16). In this way, we let the isotope curves drive changes in atmospheric O2 from the excess O2 production via sediment burial relative to O2 destruction via oxidative weathering. Subsequently, we explore the governing feedbacks at play.

The consumption of O2 through oxidative weathering of organic carbon and pyrite (Eqs. and 16) is allowed to vary in proportion to the global continental weathering rate, which is proportional to the continental Sr runoff [JW(Sr-R)] obtained from simultaneously fitting the strontium isotope evolution of seawater, RSr = (87Sr/86Sr)seawater:(see derivation in ). Although, the Sr isotope constraint depends on boundary conditions, including average isotope compositions (Rx), marine Sr inventory (MSr), and fluxes of major sources (JW,x), including continental rock weathering (Sr-R), diagenetic carbonates (Sr-DC), and hydrothermal Sr input (Sr-HT), this equation tells us that the long-term changes in continental weathering fluxes compared to other oceanic Sr sources were actually quite modest during the Terreneuvian.

Samples

Samples and geological settings are described in . Briefly, we analyzed δ238U in 187 (and δ34S in 93) marine carbonate rock samples deposited at 6 different localities across 3 paleocontinents during the Terreneuvian when small animals with skeletons and shells diversified in the oceans (Fig. 1). The sections from Xiaotan and Laolin in China, and from the Kotuikan River and the upper part of Sukharikha River in Siberia, mainly consist of limestone with low Mn/Sr, whereas the Kuonamka River section in Siberia displays atypically low U contents, and the Oued Sdas section consists mostly of dolostones. Together, these sections cover most of the Terreneuvian time interval, and specifically record 2 massive perturbations in the marine carbon cycle of hitherto-unknown biogeochemical affinity, manifested as large (>5‰) positive carbon isotope excursions recognized globally in carbonate sediments and in our studied sections (33, 39, 40). The ages of the sediments were approximated using a global age model (33, 40) with minor modifications that combines worldwide chronostratigraphic, biostratigraphic, and chemostratigraphic data worldwide ().

Fig. 1.

Sample localities and sections in (A) Morocco, (B) Siberia, and (C) South China and (D) paleogeographic reconstruction for ∼525 Ma, modified after ref. 33 (but see alternative reconstruction in ref. 74).

Results

Uranium Isotopes.

Our U isotope data are plotted as a function of time in Fig. 2, and as a function of stratigraphy at each locality in . There is a clear systematic coupling between the marine carbon and uranium isotope records, expressed as large, negative δ238U excursions (of magnitude −0.4‰) correlated with positive carbon isotope excursions (δ13C). This coupling occurs in several sections worldwide and is not expected for any known diagenetic processes. Furthermore, the large isotope oscillations do not correspond to geochemical and petrographic signals of diagenetic alteration. Instead, we suggest these events were associated with global expansions of anoxic water masses in the oceans.

Fig. 2.

(A) Carbon and (B) uranium isotope data of carbonate associated uranium (δ238UCAU) during the Terreneuvian (∼539 to 521 Ma). Stratigraphic correlations are obtained at estimated ±1-My relative age resolution using carbon isotope chemostratigraphy, consistent with biostratigraphic records, and tuned to the age-model obtained from chronostratigraphic constraints in Morocco and China (33) (see for details). An interpolated curve through minima (δ238UMIN, black curve) represents the best estimate of the composition of Terreneuvian seawater, and it parallels the smoothing spline fit obtained through least altered carbonates (LAC) (δ238ULAC, red curve). The δ238ULAC and δ238UMIN curve carry smaller average offsets from coeval seawater than average marine carbonates (∼0.1 to 0.2‰ and ∼0.0 to 0.1‰, respectively). The composition of the oceanic input is assumed to be similar to modern rivers and continental rocks at −0.30‰ (18).

That said, the δ238U record shows both systematic oscillations and single outliers reflecting a combination of secular changes in seawater chemistry and isotope variability likely generated locally during diagenetic modifications within the sediments. Although the vast majority of our samples experienced a low degree of diagenetic alteration, recorded as Mn/Sr < 1 and no dolomite (Mg/Ca < 0.1), samples with a more fluid-buffered (open-system) diagenetic alteration signature also are associated with distinctly higher average δ238U (+0.14‰) compared to the least altered samples (P < 0.01 for both measures; see for details). Accordingly, we also find that samples with Sr/Ca below the median indicative of Sr (and U loss) through open-system diagenesis on average display +0.11‰ higher mean δ238U value (P < 0.01 for all 3 measures; see for details). This offset is within the range occurring during syndepositional or very early diagenesis on the modern Bahamas banktop (+0.24 ± 0.14‰, 1 SD), suggesting the least altered samples are, on average, ∼0.1‰ closer to contemporaneous seawater than modern carbonates are to modern seawater (29, 31, 36). We use a smoothing spline fit through the least altered carbonate samples and obtain a coherent curve, δ238ULAC, and compare this signal to a curve obtained by interpolation between the lowest δ238U values within 0.4-My intervals, δ238UMIN, representing samples with the smallest offset from contemporaneous seawater. The offset between these 2 curves is on average ∼0.1‰, which means that δ238UMIN, on average, is even closer to contemporaneous seawater and resembles δ238U of Terreneuvian seawater with only a small average offset of 0.09‰ (). In fact, the scatter in Terreneuvian samples matches well with that observed in modern carbonates. Importantly, the δ238UMIN and δ238ULAC curves are parallel and both capture the major oscillations. Therefore, we infer that δ238UMIN best resembles secular variability in Terreneuvian seawater with values between −0.9‰ and −0.4‰, and we adopt a conservative uncertainty estimate of ±0.14‰ at any given time step (1 SD).

Carbon, Sulfur, and Strontium Isotopes.

Published carbon isotope records from Mongolia, Siberia, China, and Morocco provide the chemostratigraphic framework for the age model () (33, 39, 40, 42). Our uranium isotope data are negatively correlated to the 2 largest positive δ13C excursions (>4‰) characteristic of expansions of marine anoxia at the global scale. This correlation is robust, even if the age model is altered.

To constrain changes in the Earth’s atmospheric O2 budget, we use the carbon isotope data (δ13CCARB) and compiled sulfur isotope data of carbonate-associated sulfur (δ34SCAS) from the same carbonate rocks where 87Sr/86Sr, and δ238UCAU also are recorded. The δ34SCAS data (Fig. 3) shows persistently high values (30 to 45‰) in Terreneuvian seawater compared to modern seawater and carbonate banktop sediments (21‰) (43) that are weakly correlated to the δ13CCARB excursions. The δ34S data fills major gaps in the existing sulfur isotope record (44, 45). There is no covariation between δ34SCAS and diagenetic indicators. If the δ34SCAS record resembles the isotope composition of overlying seawater sulfate, then the high δ34SCAS values indicate that pyrite burial accounted for the majority of total marine sulfur burial (∼60 to 100%). This calculation assumes the marine sulfur cycle was at steady state and that the δ34S of oceanic input (4 to 8‰) and average isotope offset of buried pyrite relative to contemporaneous seawater (35‰) were similar to modern day (Dataset S3).

Fig. 3.

Isotope composition of carbonate-associated (A) carbon (δ13CCARB) (39, 46, 75), (B) sulfur (δ34SCAS; this study), and (C) strontium (87Sr/86Sr) (33, 46, 76). Smoothing splines were fitted through the data in A and B to represent the chemical evolution of seawater. Diagenetic effects are known to shift 87Sr/86Sr to higher values than contemporaneous seawater, so the Sr isotope composition of latest Ediacaran to Terreneuvian seawater would evolve through minimum values of available records. The red dotted curve represents the seawater 87Sr/86Sr calculated from the weathering forcings in Fig. 5 and the ordinary differential equation “ODE model” described in Eq. .

Fig. 5.

Atmospheric O2 levels in the Terreneuvian is simulated from (A) organic carbon burial (Eqs. and ) and (B) pyrite sulfur burial (Eqs. 11 and 12) using typical Cambrian parameter values (Dataset S3) obtained from seawater carbon and sulfur isotope curves and forced by a shift toward more continental weathering 524 to 523 Ma. (C) The weathering fluxes forcings applied to fit the Sr isotope data in Fig. 3 is shown together with the expected oceanic Sr residence time. (D) O2 release to the atmosphere from organic C burial (light blue) and pyrite S burial (blue) are the principal sources of atmospheric O2. In combination with global O2 removal through oxidative weathering (blue curves in A and B), we estimate atmospheric pO2 change. (E) The absolute atmospheric pO2 curve (black) is fit to a model-dependent constraint at 521 Ma 4 to 7 atm% (black, dashed curves). A shift from preferential aragonite to calcite precipitation in the ocean ∼525 Ma (here shown a few million years too early to distinguish the effect from enhanced weathering) would have changed the Sr residence time (33), but has no effect on the atmospheric pO2 curve. Note that enhanced weathering both affects atmospheric O2 supply through sediment burial and atmospheric O2 removal through oxidative weathering, and therefore has little effect on atmospheric pO2 in this model.

Published strontium isotope data from sedimentary carbonates through the Neoproterozoic–Cambrian interval show a long-term 87Sr/86Sr decline from ∼550 Ma indicative of enhanced basaltic weathering, perhaps by increased seafloor spreading or warming abyssal waters (33, 46) that was finally broken at the inflection point around ∼524 Ma, suggesting the relative amount of continental versus seafloor weathering increased (Fig. 3). A simple model for the marine Sr budget shows that a 20% increase in continental weathering relative to basaltic weathering alone is sufficient to cause the inflection at ∼524 Ma (see details in ).

Discussion

Massive Expansions in the Extent of Seafloor Anoxia.

The uranium isotope record indicates that the modeled oxygenation state of the early Cambrian oceans was neither steady nor monotonous (7–9, 47, 48), but rather oscillated dramatically in concert with changes to the global marine carbon cycle as indicated by the δ13C curve. Fig. 4 shows the calculated proportion of total marine U burial occurring in anoxic marine settings (fU) from the estimated seawater δ238U trends. The calculation assumes a constant U input from rivers and a constant isotope fractionation between seawater and the U sinks. The change in weathering regime favoring more granitic than basaltic rock weathering should not affect the concentration or isotope composition oceanic input (18, 49). The modern oceans have a homogenous δ238U composition of −0.39 ± 0.04‰ (all errors are 2 SE) slightly below riverine input and crustal rocks (−0.30 ± 0.03‰) due to isotope fractionation occurring when U is buried in reducing sediments (18). Given that 12 to 25% of marine U burial occurs in anoxic sediments with a large positive +0.5 ± 0.1‰ offset from seawater today, all other U sinks (mainly carbonates and reducing sediments below oxic waters) carry an average isotopic offset of −0.01 ± 0.06‰ (18, 45, 50). Applying these constants to the Terreneuvian implies that the early Fortunian (∼539 to 532 Ma) oceans were initially less oxygenated than today (δ238U = −0.6 ± 0.1‰, fU = 59 ± 27%). Over millions of years, the modeled oxygenation state of the oceans oscillated between extreme anoxia and extreme oxia with respect to marine U cycling. At peak anoxia, essentially all marine U burial occurred in anoxic settings (δ238U = −0.80 ± 0.05, fU = 100 ± 28%) and oceans were as oxygenated as today during peak oxia in the immediate aftermath of the events (fU ∼ 18 ± 16%, δ238U = −0.40 ± 0.05‰). In the latest part of the Terreneuvian (“Tommotian,” ∼525 to 521 Ma, fU = 69 ± 22%, δ238U = −0.65 ± 0.05‰), oceans returned to an intermediately anoxic state before the major radiation of arthropods, when a brief and dynamically distinct ocean oxygenation episode is inferred from Siberian sections at 521 to 520 Ma, as seen in Figs. 2 and 4 (45).

Fig. 4.

(A) The fraction of total anoxic sediment burial obtained from the δ238U record (fU) and (B) predicted δ98Mo evolution of Terreneuvian seawater. The best estimated fU (A, black curve) is obtained from δ238UMIN and a conservative underestimate obtained from δ238ULAC is shown for comparison in red. (B) Seawater δ98Mo curves predicted from the δ238UMIN (via fU) according to Eqs. 8 and 9 for various γ, where γ is the Mn-oxide proportion of all oxic Mo burial. A high γ-value (0.75 to 1.0) is inferred before the OAEs and lower γ (0.25 to 0.5) is acceptable after ∼525 Ma. Since euxinic shales are fractionated to lower δ98Mo values than overlying seawater during nonquantitative Mo removal, this U–Mo isotope comparison suggest that δ98Mo in ∼90% of the samples must have been fractionated by −2 to −1‰ relative to contemporaneous open oceans. The local fractionation occurs in settings where total hydrogen sulfide levels are either low or only intermittently present (e.g., <100 µM at pH = 8). A few samples (1 phosphorite and 3 shales of unknown redox affinity) plot in the forbidden zone above the predicted curves. However, the uncertainty of the age is large enough (±1 My or more) that this is an artifact of imprecise correlation between the siliciclastics and carbonate sedimentary records, when biostratigraphy is poor.

A first-order estimate for the areal proportion of anoxic seafloor (AANOX) can be derived from fU and the ratio of average U burial rates in anoxic and oxic continental shelf settings (rANOX/rOX = 15 ± 10; ):The estimated seafloor area only includes areas where sedimentation occurs at rates relevant for continental shelves, since anoxia in areas with little or no sediment accumulation will go undetected. The U burial fluxes in oxic and anoxic settings are assumed to scale with seawater uranium concentration in a similar fashion and therefore cancel out in Eq. . The resulting anoxic seafloor coverage on a sedimentation-weighted basis is shown in Fig. 6. During peak anoxia, we find that anoxic water masses might have expanded across 70 to 80% of the continental shelf areas.

Fig. 6.

Comparison of O2 levels on Earth, inferred new productivity, and animal biodiversification. (A) Simulated atmospheric pO2 (from Fig. 5) with absolute uncertainty shown in red and green. (B) The proportion of anoxic seafloor (AANOX) is calculated from fU (from Fig. 4) using a linear relationship between oxic and anoxic burial fluxes and their respective areal coverage, and anoxic U burial rates that are 15 ± 10 times higher than in oxygenated settings (see for details). (C) Estimated global marine new production normalized to modern day using the average and high pO2 estimate (black and blue curves in A) and the high AANOX estimate (dashed curve in B). Other combinations of pO2 and AANOX curves leads to erroneous low productivity at times and have been omitted. (D) Average bioturbation index (scale: 1 to 6) and maximum bioturbation index (59). (E) Number of new animal genera appearing in 2-My time bins during the Terreneuvian (33, 40). After 520 Ma, animal diversity further explodes (>10-fold) (77, 78).

These OAEs in the Terreneuvian would have affected all redox-sensitive metal marine cycles and their isotopes. For example, the oceanic U inventory is predicted to decrease ∼2- to 4-fold during some Phanerozoic OAEs, due to higher U removal rate in anoxic settings relative to oxic settings (19, 23, 51). Interestingly, uranium enrichments in Moroccan carbonates systematically oscillate with about the predicted amplitude in concert with the putative Terreneuvian OAEs (see details in ) (33). However, most other sections do not show clear stratigraphic patterns in U contents and U/Th ratios. A dramatic increase in U abundances in shales has been interpreted as evidence for Terreneuvian ocean oxygenation (8). However, the U content in both carbonates and shales predominantly reflect processes in the local depositional environment that affect the U content of sediments more than the subtle variations in the size of the marine U inventory. For example in carbonates, the U content can vary 100-fold depending the relative abundances of primary carbonate minerals and the extent of early diagenetic U loss. U enrichments found in shales are dependent on their total organic carbon content (52). A careful inspection of the U shale record on a total organic carbon-normalized basis and statistical analyses () allows us to reject a previous conjecture that the oceanic U inventory steadily increased from the Neoproterozoic into the Cambrian (8).

The changes in global ocean oxygenation state can be independently verified in the molybdenum isotope record, especially if hydrogen sulfide was abundant in the expanding anoxic water masses. This connection appears because both molybdenum and uranium are well-mixed elements in the oceans, and the isotope compositions of seawater (δ98Mo, δ238U) are both positively related to the proportion of global sediment burial in oxygenated settings. Local processes tend to drive the 2 isotope systems in opposite directions, such that the carbonate δ238U record is equal to or higher than seawater δ238U, therefore overestimating global ocean oxygenation, whereas the δ98Mo record in shales adopt seawater δ98Mo or lower values, therefore underestimating global ocean oxygenation. In principle, this relationship means that the combination of the two isotope systems brackets the global ocean oxygenation state. In Fig. 4, we show the predicted seawater δ98Mo curves derived from the δ238U record if changes in seawater composition were dictated primarily by sediment burial in anoxic marine settings (details in ). The comparison of model and data suggests that ∼90% of all shale samples carry a −1 to −2‰ offset of the δ98Mo values relative to the predicted seawater δ98Mo composition, in good agreement with observations from modern open marine euxinic settings (45). A few samples plot in the forbidden zone above the modeled seawater δ98Mo curves. However, this offset can be ascribed to inaccurate age assignment of the shales (up to ∼1 to 2 My), since there is no chemostratigraphic or chronostratigraphic control that can be used to correlate the siliciclastic sediments to the carbonate record at this temporal resolution. Therefore, the δ238U record from carbonates can be viewed as superior to the existing shale-based Mo isotope record in terms of reliably recording secular changes in global extent of seafloor oxygenation.

Nevertheless, the δ98Mo record and model independently verify the large oscillations in marine oxygenation and have low values, similar to the Proterozoic average (∼1.0‰), during peak anoxia and values as high as the modern (up to 2.3‰) during peak oxia, corresponding to the proportion of oxic Mo removal that occurs with marine Mn-oxide burial. The Mo isotope data demand a substantial Mn-oxide sink for Mo before and between the OAEs (25 to 75%) comparable to today (41 ± 3%) and somewhat smaller (0 to 25%) after the OAEs (53). As a large part of this Mo sink occurs in the deep abyssal oxygenated settings in the ocean today, the high Mo isotope values (>1.5‰) before and in between OAEs suggest that oxygenation was a global phenomenon that expanded well beyond the continental shelves and into the deep ocean where Mn-oxide deposition is common in oxygenated areas (14).

Atmospheric pO2 Decline, Oscillations, and Modest Rise in the Terreneuvian.

We now turn to the question of whether seafloor oxygenation was linked to changes in atmospheric O2 levels. Although, the C, S, and Sr isotope records do not uniquely determine absolute atmospheric pO2, and shelfal processes and diagenesis might push these signals in different directions, our point here is to emphasize that a straightforward solution to the observed dynamics exists that consistently fits with all observed isotope records and with our current understanding of the governing feedbacks between atmospheric and oceanic oxygen levels.

Atmospheric O2 production (Fig. 5) results from global organic carbon and pyrite burial fluxes that can be calculated from marine carbon and sulfur isotope mass balance (37). Solving for the burial fluxes requires both the isotope records and their first time derivative. To avoid amplification of noise in the derivative curves, we first smoothed the δ13CCARB and δ34SCAS records using cubic smoothing splines. The splines represent a visual best fit to the data (Fig. 3) and preserve the important timing and magnitudes of observed excursions.

Continental weathering of phosphorous (P) modulates global marine productivity, organic carbon burial, and atmospheric O2 levels on our planet because P is the biolimiting nutrient in the oceans over million-year timescales, where atmospheric O2 levels change. In the absence of mycorrhizal symbionts that amplify P weathering by plants on the continents today, the Terreneuvian P weathering rate would have been lower than today (54). We fixed P weathering at one-quarter of the modern-day value and explored the O2 dynamics required by the isotope records. The relative proportion of continental vs. basaltic weathering varied in a fashion that affected the strontium isotope composition of seawater (Fig. 3, dotted red curve). High rates of seafloor spreading and ocean crust alteration during the early Terreneuvian is suggested to account for the 87Sr/86Sr decline and arguably lowered seawater Mg/Ca to the extent that calcite rather than aragonite predominantly precipitated in the marine realm from ∼524 Ma and well into the Paleozoic (33, 55). The inflection in the strontium isotope record at this time fits with a modest ∼20% relative increase in continental weathering, and this estimate is unaffected by the higher Sr residence times in calcitic seas (Figs. 3 and 5). The concentrations of dissolved inorganic carbon and sulfate in early Cambrian seawater have been estimated to ∼9 and ∼1 mM, respectively (33, 45, 56), and the organic carbon and pyrite proportions of weathered material are kept constant at characteristic values for marine sediments at this time (RORG,W = 0.20, RPY,W = 0.90). A sensitivity analysis shows that the atmospheric pO2 trajectory depends strongly on the absolute P weathering rate, the average proportions of organic carbon to total carbon burial and pyrite sulfur to total sulfur burial in weathering products and, to a lesser degree, the oceanic sulfate pool (Dataset S3). Although these parameter values affect the baseline atmospheric pO2 level and the amplitude of pO2 change, they cannot drive the observed dynamics. The age model is uncertain, and that also affects the predicted atmospheric pO2 trajectory, but we suggest that the general features discussed below are robust. The change in continental weathering at ∼524 Ma only had little influence on simulated atmospheric pO2 because weathering promotes both O2 consumption by oxidative weathering and O2 production by organic carbon sediment burial in our model (Fig. 5). The calculated pyrite burial flux is rather stable and acts as an O2 source throughout the Terreneuvian. The calculated organic C burial flux varies dramatically through the early Terreneuvian and occasionally yields erroneous, negative values for brief periods when seawater δ13C dropped rapidly to values below volcanic input of −5‰. These intervals are attributed to preferential oxidative weathering with excess organic carbon delivering carbon with lower δ13C signature into the oceans and acting as an O2 sink at these times.

The changes in atmospheric pO2 are reconstructed from the isotope records, but the absolute pO2 levels are difficult to constrain. We anchored the atmospheric O2 curve on a recent model-dependent constraint of 4 to 7 atm% at ∼521 Ma (45). Remarkably, in light of the simplicity and uncertainty of the model, this approach readily leads to a physically reasonable (nonzero) atmospheric pO2 curve. The resulting atmospheric pO2 trajectory shows a notable decline from ∼7 atm% at the Precambrian–Cambrian boundary to a minimum of ∼4 atm% ∼533 Ma followed by rising O2 levels during the OAEs before ultimately settling at ∼6 atm% after the second OAE.

Global New Production.

To distinguish between reduced O2 supply from the atmosphere or O2 removal within the oceans as the key driver for expanding ocean anoxia, we provide an estimate of global marine photosynthetic production, which consumes O2 as it sinks into the deeper parts of the oceans. The first-order relationship between global marine new production (newp′), global extent of anoxia (AANOX) and atmospheric pO2 (o′) has been obtained from Earth system ocean–atmosphere models (57):The constants kanox = 12, ku = 0.5 are determined so that 0.25% of the seafloor is overlain by anoxic water masses similar to modern productivity and oxygen levels (newp′ = o′ = 1) and the formula gives a fairly sharp transition of AANOX from 23 to 77% as atmospheric pO2 declines from 13 to 8 atm%, consistent with biogeochemical models for the spatial distribution of oxygen in the oceans (57, 58). Clearly, this relationship would have changed in concert with changes in organic matter export rates and how ocean circulation ventilated the deeper parts of the oceans. Until values for the Cambrian ocean have been established, we suggest that the marine productivity curve obtained this way is qualitatively correct. This approach consistently shows that the Cambrian oceans were distinctly less oxygenated than today due to lower atmospheric pO2, and that the overall availability of nutrients oscillated in concert with the OAEs.

We observe 3 key qualitative features of the simulated global net biological productivity during the Terreneuvian that help explain the driving mechanism for early Cambrian OAEs (Fig. 6). First, new photosynthetic production and nutrient availability in the early Fortunian slowly declined from the Precambrian–Cambrian boundary until ∼533 Ma in concert with atmospheric pO2 levels. Second, both OAEs were associated with high nutrient levels promoting massive blooming of phytoplankton worldwide. Third, a significant, albeit small, increase in marine new production after the OAEs did not drive a third expansion of anoxia. Instead, atmospheric pO2 rise kept pace with the enhanced biological production reaching slightly higher O2 levels than the Fortunian and resulting in no further expansion of oceanic anoxia.

Driving Factors for Terreneuvian Atmosphere–Ocean Oxygenation Dynamics.

The combination of the simulated global seafloor oxygenation, atmospheric pO2, and marine productivity allows us to assess plausible first-order causal feedbacks at play in the Earth system during the Terreneuvian.

First, the mutual decline in new photosynthetic production and atmospheric pO2 are the key features predicted for nutrient burial in ventilated sediments by bioturbating animals (7). Theory predicts that, as animal activities stir marine sediments, they cause the biolimiting nutrient phosphorous to be retained within organic matter and iron oxyhydroxides in marine sediments. Interestingly, the first widespread motile benthic animal communities, plowing horizontally on the seafloor, emerged near the Ediacaran–Cambrian boundary (59, 60). Therefore, we conclude that benthic animal communities may have limited global ocean phosphorous availability, new marine photosynthetic production, and initiated atmospheric pO2 decline (7).

Second, an increase in marine productivity hit when atmospheric pO2 levels were at their lowest in the Terreneuvian (e.g., ∼4 atm% at ∼533 Ma). Marine productivity continued to rise and peaked with the first OAE. The prolonged interval leading into the first OAE appears much longer than most OAEs in the Phanerozoic, and in contrast to younger OAEs triggered by volcanism, there is no evidence of extensive volcanism to trigger expansive anoxia at this time (33, 61). Hence, we suggest that limited O2 supply from the atmosphere could have initiated the expansion of anoxia in some areas of the continental shelves, and that expanding anoxia ultimately continued via an accelerating feedback loop acting to increase productivity by efficient recycling of P released from expanding seafloor areas covered with anoxic bottom waters (62).

The atmospheric pO2 decline would have slowly made different parts of the oceans more prone to anoxia depending on the productivity, organic carbon export, water circulation, and O2 supply into the deeper parts of each basin. When sustained for more than ∼105 to 106 y, excess organic carbon burial acts to increase atmospheric pO2 to levels where P recycling ultimately collapses and ocean basins again experience higher oceanic O2 supply from the atmosphere (63). Apparently, our modeled atmospheric pO2 levels show that O2 increased from 533 to 529 Ma while the ocean anoxia also continued to expand at the globally integrated scale. Thus, atmospheric O2 had to increase significantly above the levels that initiated expansive anoxia in the first place, before runaway anoxia was ultimately halted.

The second OAE indicates that the system did not immediately stabilize itself. Rather, atmospheric pO2 again crashed to low levels (∼5 atm%) and triggered a second expansion of marine anoxia lasting for another ∼1 to 2 million years. At least 3 potential triggers for the second event exist. First, the Earth system could have entered a state with natural oscillations, as suggested for the Cretaceous OAEs (63). Second, a change in ocean circulation patterns perhaps linked to slowly rising sea levels may have stratified the oceans and made some parts of the oceans more prone to expanding anoxia even without atmospheric pO2 decline, once again initiating the positive-feedback loop involving P regeneration in anoxic basins (64). Third, benthic animal communities evolved new infaunal feeding strategies in low-energy offshore deposits and extended the maximum burrowing depth and average bioturbation index at this time [Cambrian stage 2, ∼529 to 521 Ma (59); but see ref. 65 for a different view]. Vertical sediment mixing and enhanced irrigation could have accelerated nutrient removal from the oceans into oxidized sediments and set the system out of balance once again, triggering renewed atmospheric O2 decline.

Notably, the atmosphere–ocean O2 oscillations did not produce a third event of similar size. Instead, from ∼525 Ma onward, atmospheric pO2 levels were higher (∼6 atm%) and above the critical levels where O2 supply had previously initiated OAEs (∼4 to 5 atm%). From ∼524 to 523 Ma onward, the modest expansion of anoxia, unrelated to atmospheric pO2 decline, we see evidence for increased marine O2 consumption either due to external nutrient injection (e.g., from the continents) or enhanced nutrient recycling (e.g., by ocean stratification).

In view of a comparatively modest relative increase in continental weathering at this time consistent with the strontium isotope record, we suggest enhanced physical erosion on the continents is a plausible driver that also would have promoted nutrient loading into the oceans as well as sediment burial in deltaic sediments. Deltaic settings are oxygenated and preferentially affect the marine C and S cycle through organic carbon and pyrite sulfur burial (to a greater degree than deltaic burial affects the marine U cycle) (19). As such, sediment burial in deltas preferentially releases O2 into the atmosphere and surface oceans that may not have balanced the effect that enhanced weathering had on oceanic nutrient loading and seafloor anoxia. In this view, the change in weathering regime interrupted the anoxic oscillations at ∼524 Ma and both increased O2 and food availability for marine animals from this time onward.

Implications for Animal Evolution.

Atmospheric pO2 levels dictate the O2 availability in the surface oceans. For example, air-saturated water at 25 °C with a salinity of 35‰ would, at equilibrium, have an O2 concentration of ∼50 μmol·kg−1 at 4 atm% pO2. Our results show that the dissolved O2 concentration in the Terreneuvian surface ocean would have varied between ∼40 and ∼110 μmol·kg−1, a range well above the level of severe hypoxia ∼9 μmol·kg−1 where carnivorous taxa disappear (8). Hence, the Terreneuvian surface ocean was sufficiently oxygenated to support both predators and prey living in seawater at all times.

To test hypotheses for the coevolution of animal life and the environmental change, we compare our oxygenation and productivity records to records of fossil first appearances during this time (33, 40). The fossil record show 3 or 4 pulses of animal appearances in the Terreneuvian recorded in distinct parts of the world. A first pulse of skeletal animal appearances is recorded in Mongolia 536 to 534 Ma and was soon followed by 1 or 2 pulses in China and Siberia 532 to 530 Ma. This was followed by an interval ca. 530 to 525 Ma with relatively few first appearances, before a final “Tommotian” pulse of appearances began in the late Terreneuvian ca. 525 Ma (Fig. 6) (33, 40). Looking for links between biological evolution and environmental change must be considered tentative, since it is still not known to what extent the appearance pulses reflect stratigraphic and preservational biases, or originations of new animal species (40).

The first pulse of animal appearances recorded in Mongolia 536 to 534 Ma occurs at a time when inferred atmospheric O2 levels declined and marine productivity were rather stable. Although anoxia likely expanded diachronously, enhancing nutrient recycling and productivity in various ocean basins at various times, our records cannot causally link this pulse to a global environmental trigger (66, 67).

The following fossil appearance pulses in China and Siberia at ∼532 and ∼530 Ma occurred during both rising productivity and surface ocean O2 levels that would have elevated food and O2 supply to the marine food webs. Experiments show that these factors can promote animal diversification (12). Further investigations on the diachroneity of these events are needed and local paleoredox and/or productivity proxies with high spatial resolution should address whether fossil appearances were correlated to enhanced biological productivity at these times.

After the 2 expansions of seafloor anoxia, animal fossil appearances rose around 524 ± 2 Ma (Tommotian) in both Siberia and China when predicted O2 levels in the surface oceans were higher and the seafloor was also rather oxygenated, despite slightly higher productivity. A change in weathering regime at this time could have settled atmospheric pO2 at higher levels than in the earlier Terreneuvian, despite increasing levels of bioturbation. With more food and O2 available, the carrying capacity for prey and carnivores in the marine ecosystem would have been higher, and this could have promoted further animal radiations (12).

We note that the rate of fossil first appearances was generally lower throughout the intervals with expanded anoxia, although new animal fossils begin to emerge in Siberia during the second OAE, and there are few new fossil taxa recorded in the oxygenated period between the OAEs. We suggest that the Tommotian origination pulse was linked to seafloor oxygenation and the change in weathering regime, similar to later Cambrian bioradiation events (68). We also acknowledge that that the rate of new appearances during the intermediate oxygenated interval has possibly been underestimated either because of sampling or preservational biases or because this particular time interval was shorter than predicted by the age model (40).

Finally, in the greater context of what triggered the emergence of complex animal ecosystems on Earth, our isotope data and models suggest that there were ample amounts of O2 for animals in the surface ocean at the onset of the Cambrian (69). Our model also shows that atmospheric pO2 levels declined in concert with global new marine production in the earliest Terreneuvian, which is the predicted result of evolving benthic animal communities that promoted phosphorous burial in bioturbated sediments (7). The first pulse of fossil appearance rates occurred before peak eutrophication and may not have been environmentally triggered. The late Terreneuvian (Tommotian) pulse occurred at a time with both elevated atmospheric O2 levels and marine productivity likely associated with a change in weathering regime. Ultimately, the trigger for the developing networks of ecological interaction lies in biological innovations (e.g., motility, eyes, armor), but our results suggest that the rate at which animal ecosystems evolved was modulated by O2 availability as well as chemical energy flow into the marine food web—both of which, in turn, evolving animal ecosystems impacted themselves.

Methods

Carbonate-Associated Elements.

Elemental analyses including U, Th, Mn/Sr, and Mg/Ca were carried out on the mild acetic leachates using Thermo Scientific X-Series 2 quadrupole mass spectrometer at the University of Copenhagen. The leaching procedure utilize a 10% acetic acid with pH >3.6 for <16 h calibrated to maximize the carbonate-associated uranium (CAU) yield, while minimizing the extraction of uranium associated with carbonate fluorapatite, ferric oxyhydroxide-associated P, silicates, and organic-bound U (70). The extraction yield was typically ∼50 to 70%.

Mineralogy.

The bulk mineralogy was characterized by powder X-ray diffraction using a Bruker-AXS D8 Advance diffractometer equipped with the primary beam Ge(111) monochromator, Cu tube (Kα1 wavelength, 1.54059 Å), and silicon strip detector (LynxEye) at the University of Copenhagen. The minerals documented in each sample were identified using the Powder Diffraction File (PDF) database of the International Centre for Diffraction Data.

δ238UCAU.

CAU was obtained using a calibrated soft leaching procedure that mainly extracts uranium in calcite and avoids U associated with authigenic carbonate fluorapatite and ferric oxyhydroxide-associated phosphate minerals or any detrital or organic-bound U. The procedure is described in ref. 45.

For the isotope analyses, the concentrations of uranium and thorium were measured on ∼1% of the leachate, and the remaining fraction (99%) or up to 600 ng U was used for isotope analysis. A calibrated 233U–236U double spike (IRMM-3636) was added to the analyte and allowed to equilibrate with the sample to generate a homogeneous mixture with 235U/236U ∼ 1 and 235U/233U ∼ 1 (71). Chemical purification was conducted using 4 sequential ion chromatography columns to ensure high U purity of the analyzed solution similar to the reference standard solution (CRM-145) (72). The spike/sample ratio was kept constant for all samples and standards during each analytical session. Each sample solution was measured in 5 replicates at signal intensities matched to the bracketing standard solutions. The long-term reproducibility of a matrix-matched reference material, SRM-1d is better than 0.03‰ (2 SE, n = 5). Therefore, we estimate the long-term reproducibility and accuracy of our 238U/235U data to be ±0.03‰ (2 SE). Other reference materials include BCR-2 (whole rock) and IAPSO seawater at δ238U = −0.271 ± 0.017‰ and −0.372 ± 0.022‰ (2 SE, n = 5), respectively, which are indistinguishable from results obtained in other laboratories (18, 45).

δ34SCAS.

The isotope composition of carbonate-associated sulfur was obtained using a traditional approach at Virginia Polytechnic Institute and State University (73). Briefly, 3 to 8 g of fine sample powder were treated with a NaCl solution, followed by several distilled water rinses to remove any soluble sulfates that might have been present in the samples. The samples were then dissolved using 4N HCl for less than 1 h. The leachates were centrifuged and vacuum filtered (45 µm) to remove any insoluble portion of the sample, before ∼100 mL of saturated BaCl2 solution was added to allow sulfate to precipitate as BaSO4 from the solution. Samples were left for at least 3 d to ensure complete precipitation. Barium sulfate was separated from the remaining solution via filtration. The 34S/32S isotope ratios were analyzed through on-line combustion using an Elementar Isotope Cube elemental analyzer coupled to Isoprime 100 continuous-flow stable isotope ratio mass spectrometer. Data are reported in standard delta notation (δ34S) as per mil deviations from Vienna Canon Diablo Trolite (VCDT). International isotope standards NBS-127, IAEA SO-5, and IAEA SO-6 were used to convert measured isotope ratios to the VCDT scale. Long-term reproducibility of the international isotope standards was <0.2‰ (1 SD). The external reproducibility of our δ34S data are estimated to be better than 1‰ based on replicate dissolution of samples (45).

Data and Materials Availability.

All data needed to evaluate the conclusions in the paper are present in the paper and/or .