Geochemical evolution of the early Palaeozoic Earth system

Geochemical data covering the late Neoproterozoic and early Phanerozoic (560–420 Ma) are summarised in Fig. 1. In general, there is evidence for considerable heterogeneity in the evolution of ocean redox chemistry from the late Neoproterozoic through to the mid-Cambrian, as might be expected in a pervasively low oxygen world. However, much of this apparent heterogeneity may be due to difficulties in adequately sampling sediments from a range of water depths in individual studies. For example, the existence of an anoxic oxygen minimum zone (OMZ) along productive continental margin settings has been advocated for the early Cambrian28. Samples from within the OMZ would give a very different redox signature compared to samples from oxic shallower and deeper waters. Nevertheless, an overall progression in ocean oxygenation across this period is now emerging, aided by the application of redox indicators that provide a more global indication of ocean redox chemistry29,30,31.

Fig. 1 Compilation of geochemical data from the Neoproterozoic into the Palaeozic. a Mo abundances, as compiled by ref. 29. Higher oxygen levels lead to higher abundances of Mo. Note that other proxies providing support for the ocean oxygenation state outlined here are compiled in Supplementary Fig. 1. b Sulphate-S isotopes (δ34S SO4 ) as compiled by ref. 69. δ34S SO4 increases through the Cambrian, with a return to lighter values at the GOBE. c Pyrite fraction of sulphur burial (f pyr = pyrite burial/(pyrite burial + gypsum burial)), full range of estimates, as presented in ref. 44. d Carbonate-C isotopes (δ13C carb ) as compiled by ref. 70. Higher values indicate a higher rate of organic carbon burial. e Mixed layer depth, reproduced from ref. 7. Black lines on panels b and d are local regression (LOESS) fits. Grey shaded areas indicate the Cambrian explosion (540–521 Ma) and the Great Ordovician Biodiversification Event (GOBE; 470–450 Ma). Blue shaded line indicates the Hirnantian glaciation Full size image

There is robust evidence for widespread deep ocean oxygenation in the late Neoproterozoic30,32,33,34,35. Building upon this, selenium isotope evidence30 also suggests progressive oxygenation through the Neoproterozoic. However, multiple lines of evidence demonstrate a short-lived return to widespread ocean anoxia at the Precambrian-Cambrian boundary (Fig. 1a and Supplementary Fig. 1)33,36,37,38. This anoxic episode was likely too short-lived to be captured by our modelling approach (see below), but evidence from sedimentary Mo concentrations and isotopes, U isotopes, and rare earth element concentrations suggests that the global ocean then became progressively oxygenated through the early Cambrian up until the height of the Cambrian explosion at ~520 Ma29,31,34. These same redox proxies demonstrate a subsequent return to more widespread anoxia after ~520 Ma (for simplicity we demonstrate this with Mo concentration data in Fig. 1, but other proxies are compiled in Supplementary Fig. 1), which is also consistent with evidence from S isotopes, Fe speciation, and trace metal concentrations indicating widespread euxinia in the later Cambrian ocean39.

Evidence from sulphur isotope systematics suggest this widespread euxinia continued during the early and mid-Ordovician40. Other geochemical records through the Ordovician are relatively sparse, although ocean redox conditions appear to have been subject to temporal variability, which is consistent with an overall low-oxygen world41. However, by the end of the Ordovician and during the Silurian (<460 Ma), oxygen started to increase27, coincident with the earliest instance of fossilised charcoal, indicating near-modern levels of atmospheric O 2 42. This rise in oxygen has been attributed to the evolution of land plants, which culminated in the Devonian and Carboniferous with the development of roots and vasculature27,43.

The sedimentary record of oceanic sulphate δ34S shows significant variability (Fig. 1b), but the general trend suggests an increase in δ34S SO4 from the Ediacaran (25‰) to the mid Cambrian (39‰), followed by a recovery to pre-Palaeozoic values of 25‰, just before the Great Ordovician Biodiversification Event (GOBE). Increases in δ34S SO4 are often linked to increased rates of pyrite burial, but this may not necessarily be the case during the Early Palaeozoic44 (note the very large range in f pyr ; the fraction of ocean sulphate that is buried as pyrite; Fig. 1c). A long-term perturbation during the Cambrian and Ordovician is apparent in the sedimentary carbonate δ13C record. While this record reveals fluctuations in the carbon cycle on short timescales of a few million years, on a longer time scale the δ13C carb record shows generally lower values from the middle Cambrian to the late Ordovician (~0‰) as compared to the late Ediacaran and early Cambrian (~4‰), with a return to higher values in the Silurian (~2‰). The shift in δ13C carb at around 521 Ma is consistent with a reduction in net organic carbon burial, and a resultant decrease in the production rate of O 2 . The recovery to more oxic conditions in the Silurian was likely coupled to increased organic carbon burial, as indicated by the more positive δ13C carb values (Fig. 1d).

The analysis of sediment fabric disturbance suggests early animals only reworked the seafloor superficially, with shallow burrowing appearing during the Cambrian and continuing into the Ordovician, whereas only towards the end of the Silurian do burrow systems become deeper (mixed layer depth shown in Fig. 1e)7. Various scalings of the biogeochemical response to bioturbation have been shown to be consistent with subsets of available geochemical proxies: a rapid and non-linear response at low bioturbation intensities may have increased phosphate retention in the sediment when the first shallow-burrowing animals appeared in the Cambrian22, leading to a decrease in the oxygen production source from organic carbon burial, and driving a return to anoxic ocean conditions after 520 Ma (Fig. 1a). In contrast, a more protracted response has been shown to be consistent with low sulphate concentrations throughout the Palaeozoic7. However, the proposed scenarios for the evolution of bioturbation have not been evaluated using multiple geochemical proxies or model simulations combining feedbacks between the carbon, oxygen, phosphorus and sulphur cycles.

The COPSE model

To obtain more robust constraints on the timing and environmental consequences of the rise of bioturbation, we modified the COPSE model26,27 (Carbon Oxygen Phosphorus Sulphur Evolution), which is a synthesis of the ‘Redfield’ model45 and the GEOCARB model46. The COPSE model simulates the coupled evolution of the major biogeochemical cycles over the Phanerozoic by describing burial, weathering and degassing processes, which transport chemical species between the atmosphere, oceans and sediments over geological timescales (Fig. 2). COPSE aims to capture trends in biogeochemical cycling over the timescale of 10–100 s of millions of years, but not shorter-term fluctuations on timescales of ~100,000 years. The model produces estimates for the global abundance of oxygen, carbon dioxide, phosphate and sulphate, alongside records of whole-ocean δ13C carb and δ34S SO4 , which can be used to test hypotheses by comparison to data. For an overview of the biotic and tectonic controls covered by the current COPSE model, see ref. 47. The evolution of bioturbation and its feedback on global biogeochemistry has not been explicitly considered in COPSE: previous model versions have implicitly assumed that bioturbation is always active and that the bioturbation intensity remains independent of the oceanic oxygenation state.

Fig. 2 Diagram of key processes in the COPSE model. a Carbon cycle. Hydrospheric CO 2 is transferred to sediments as organic C or carbonate by burial (B). Sedimentary C is returned to the ocean/atmosphere via weathering and metamorphism (W). Buried organic C is isotopically lighter than the carbon it is derived from. Burial of reduced organic carbon results in a net source of O 2 , whereas oxidative weathering of sedimentary organic carbon consumes O 2 . b Sulphur cycle. Burial of reduced pyrite is a net source of O 2 , whereas oxidative weathering of sedimentary pyrite consumes O 2 . c Oceanic phosphorus (P) cycle. Dissolved, bio-available P is delivered to the ocean by chemical weathering via rivers, and is buried either as organic phosphorus, or with iron or calcium minerals. Dashed lines show burial processes that are influenced by bioturbation (but are not considered so in the baseline model) Full size image

By using COPSE we employ a forward modelling approach, which enables a comparison of model predictions of δ13C carb and δ34S SO4 trends to the independent geological record. This contrasts with inverse modelling, where geological records are used as a model forcing, leaving no potential for quantitative testing of the model results. As with all models, comparison with the geological record requires some assumptions. Foremost, the model predicts δ13C carb and δ34S SO4 trends that are representative of the global marine dissolved inorganic carbon and sulphate reservoirs, reflecting the globally averaged operation of the long-term geochemical cycles (in essence changes in organic carbon and pyrite burial). However, changes in the geological isotope record are not solely dependent on changes in the global biogeochemical cycling of carbon and sulphur, but also incorporate possible effects of diagenesis, or evolutionary changes to the fractionation factors associated with photosynthesis and microbial sulphate reduction. Moreover, it is possible that some data represent regional signals rather than global trends. Nevertheless a quantitative comparison of our model predictions to the geological isotope record provides a useful test of the assumptions underlying the COPSE model.

For the late Ediacaran to mid Ordovician, the ‘baseline’ COPSE model (i.e., the model version as presented in ref. 27) generates stable conditions with high atmospheric CO 2 (~15× Present Atmospheric Level (PAL)), a high degree of ocean anoxia (0.8, or 80% of the ocean surface resides under anoxic water) and δ13C carb around 0‰ (Fig. 3). Broadly, these results are driven by the absence of terrestrial productivity, and a suppression of silicate mineral weathering before land plant evolution (i.e., suppressed burial of both organic and inorganic carbon). The predictions of the baseline model reveal discrepancies with the available geochemical data, which suggest more dynamic ocean redox conditions in the early Palaeozoic and higher δ13C carb values in the late Ediacaran and early Cambrian (Fig. 1).

Fig. 3 COPSE baseline model simulation. Simulation as presented in ref. 27. a Atmospheric CO 2 . b Average δ34S SO4 of seawater. c Pyrite fraction of sulphur burial. d Average δ13C carb of seawater. e Degree of ocean anoxia (1 = completely anoxic, 0 = completely oxic). f Summary of the evolution of sedimentary Mo concentrations over time. Model outcomes (in blue) are compared to δ13C carb and δ13S SO4 data and the sedimentary Mo concentrations, which is reflective of the extent of ocean oxygenation and is supported by multiple independent proxies (see Supplementary Note 1). Dotted lines in panels b and d represents a local regression (LOESS) fit to the data. Grey shaded areas indicate the Cambrian explosion (540–521 Ma) and the Great Ordovician Biodiversification Event (GOBE; 470–450 Ma). Blue shaded line indicates the Hirnantian glaciation Full size image

The effect of bioturbation on sedimentary elemental cycling

Here, we update the COPSE model parameterisation for the burial of organic carbon, pyrite sulphur and phosphorus to include a response to bioturbation. We introduce a bioturbation parameter, f biot , which provides a single measure of the biogeochemical impact of bioturbation on sedimentary cycling (0 represents no bioturbation effect, 1 represents maximum bioturbation impact). This f biot parameter must be expressed as a function of the strength of bioturbation, which is traditionally represented by the bio-mixing depth (L b ) and the bio-mixing intensity (D b ) parameters, that feature in early diagenetic models48. As D b and L b are essentially linked (L b ~ √D b , see ref. 49 for a theoretical justification), the parameter f biot may be expressed as some function of either one of these parameters (Fig. 4). Currently, there is a lack of data to constrain the exact nature of the relation between f biot and bioturbation intensity.

Fig. 4 The effect of bioturbation on sediment geochemistry. a The effect of bioturbation intensity (D b ) on the sulphur cycling rate. Model results reproduced from ref. 16. b Solid line (Scenario 1): the bioturbation impact on geochemistry (f biot ) is linearly correlated with the depth and intensity of burrowing. Dashed line (Scenario 2): the bioturbation impact is maximal with shallow burrowing, but the areal expansion of bioturbation increases gradually throughout the early Palaeozoic. Dash-dotted line (Scenario 3): the bioturbation impact on geochemistry (f biot ) is already at full strength by the end of the Cambrian explosion. c Simplified conceptual model of the effect of bioturbation on the sedimentary cycles of carbon, phosphorus and sulphur. Arrow sizes denote relative changes in flux sizes Full size image

To examine the link between the emergence of bioturbation and global biogeochemistry, we describe three different formulations for the evolution of f biot over time that are designed to represent the envelope of possibilities. In ‘Scenario 1’, the biogeochemical impact of bioturbation is assumed to be weak during the Cambrian and Ordovician, and becomes important when deeper and more intensively burrowing organisms evolve in the late Silurian and Devonian, increasing the mixing depth and creating larger burrow networks that are intensely flushed with overlying water (solid line in Fig. 4b). Throughout the Cambrian and Ordovician, the mixed layer depth may have been at least 5 times lower than the observed depths in the late Silurian and Devonian (Fig. 1d and ref. 7), and therefore we choose f biot = 0.2 (thereby assuming a linear relationship between the effects of bioturbation and the mixed layer depth). Mixed layer depths likely increased to higher levels in the late Silurian and the Devonian7. Therefore, we allow f biot to increase to 0.5 by the end of the Silurian, which implies that the biogeochemical impact of bioturbation only peaked after the Silurian (as suggested by ref. 7).

In Scenarios 2 and 3, we suggest that low depths of bioturbation mediate a disproportionally large response in sediment geochemistry. Both experimental25 and theoretical16 studies have proposed that geochemical state variables and rates (e.g., elemental cycling rates, partitioning of redox acceptors, stimulation of organic carbon breakdown) respond in a highly non-linear way to increasing levels of bioturbation, with 80–90% of the maximum response attained at low bioturbation intensities (D b < 1 cm2 yr−1 and L b < 3 cm; see ref.16). For example, diagenetic modelling shows that the cycling rate of sulphur (the number of times a sulphur atom entering the sediment column is cycled between its oxidised and reduced states before it is eventually buried) rapidly increases from 5 to 10 when D b increases from 0 to 1 cm2 yr−1, to stabilise at 11 when D b values > 1 cm2 yr−1 (Fig. 4a). This increased redox cycling of sulphur potentially inhibits the rate at which reduced sulphur compounds are buried (by stimulating re-oxidation)14,17. Additionally, experimental work has shown that meiofauna (micron-scale animals that burrow ~1 cm) stimulate organic carbon breakdown as much as large animals (that burrow >10 cm depth)25.

During the Cambrian, the mixed layer depth increased from 0 to <0.5 cm 7, which may correlate with a significant biogeochemical response. In Scenario 3, the advent of shallow mixing thus invokes a large biogeochemical effect, and so the bioturbation impact (f biot ) increases exponentially from 0 to 1 during the Cambrian Explosion, and remains constant afterwards (dash-dotted line in Fig. 4b). We also define an intermediate situation (Scenario 2), in which we assume that shallow bioturbation has a large impact, but that the areal extent of bioturbation (and accordingly f biot ) gradually increases throughout the early Palaeozoic (dashed line in Fig. 4b). Furthermore, in all scenarios, two distinct options for the response of bioturbation towards anoxia have been tested; (i) no anoxia limitation, and (ii) anoxia limitation that scales with the fraction of the ocean that is anoxic (solid vs. dashed lines in model outputs, see Methods for more information).

The effect of bioturbation on the elemental cycling of C, P and S in marine sediments is summarised in Fig. 4c. Overall, the main effects of bioturbation are driven by the increase in oxygen exposure in the anoxic part of the sediment. In a sediment without bioturbation (e.g., the Ediacaran seafloor), organic matter is broken down less efficiently11,50,51 and sulphur is more efficiently sequestered as pyrite, leading to high burial rates of organic carbon and pyrite. With the introduction of burrowing fauna, organic matter mineralisation is enhanced and pyrite is more efficiently re-oxidised14, and so the burial of carbon and sulphur is reduced3,17. At the same time, bioturbation leads to an increase in polyphosphate sequestering, which then leads to an increase in organic phosphorus burial19. We implemented these relations between bioturbation and geochemical C, S and P cycling in COPSE by adapting the model equations for marine organic carbon burial, marine pyrite sulphur burial and marine organic phosphorus burial, and introduce new parameters that describe pre-bioturbation values for organic carbon, pyrite sulphur and organic phosphorus burial (see Methods).

Overall, this model description allows us to test the working hypotheses; (i) Scenario 1: the effect on global biogeochemical cycles scales with bioturbation depth and intensity, and increased only markedly in the Silurian–Devonian, (ii) Scenario 2: the effect on global biogeochemical cycles increased gradually throughout the Palaeozoic, and (iii) Scenario 3: the effect on global biogeochemical cycles was rapid and occurred in the Early Cambrian. We tested the validity of these hypotheses by comparing the quantitative model output to the geological record.

New model results

Under a protracted response to bioturbation (Scenario 1; Fig. 5a–g), δ13C carb values slightly decrease with the early development of bioturbation at the start of the Cambrian, and remain constant throughout the Cambrian and Ordovician. After rising in the late Ordovician, δ13C carb drops gradually during the Silurian, in response to decreased primary production and decreased carbon burial due to bioturbation (Fig. 5d). In the model simulation of ocean anoxia, the limited effects of bioturbation during the early Palaeozoic imply a high rate of organic carbon burial and therefore higher concentrations of atmospheric and oceanic oxygen, and a limited prevalence of anoxia (Fig. 5e). Both of these predictions are to some extent at odds with the geochemical record, which shows evidence for significant ocean anoxia after 520 Ma 31 until the middle Ordovician (470–460 Ma)40, and does not support high δ13C carb (~2‰) throughout this interval (Fig. 1).

Fig. 5 COPSE model with the addition of the evolution of bioturbation. Scenario 1 shows the effect of the sedimentary response that scales linearly with bioturbation intensity. Scenarios 2 and 3 assume that the effects of bioturbation on sediment geochemistry occur non-linearly (strong response for low levels of bioturbation), where Scenario 2 follows a gradual increase of the areal extent of bioturbation and Scenario 3 shows the maximum effect at the Ediacaran-Cambrian boundary (see panels g, n, u). a, h, o Atmospheric CO 2 . b, i, p Average sulphate δ34S of seawater. c, j, q Pyrite fraction of sulphur burial. d, k, r Average δ13C of carbonate. e, l, s Degree of ocean anoxia. Model outcomes (in red) are fitted to the δ13C carb and δ13S SO4 proxies (grey dotted line represents a LOESS fit), predictions for the relative importance of pyrite for the total sulphur burial rate (blue dotted lines represent the range of model results presented in ref. 44) and compared to a summary of the evolution of sedimentary Mo concentrations over time (f, m, t). Solid lines represent the model outcomes with anoxia feedback, dashed lines represent the model outcomes without anoxia feedback. Grey shaded areas indicate the Cambrian explosion (540–521 Ma) and the Great Ordovician Biodiversification Event (GOBE; 470–450 Ma). Blue shaded line indicates the Hirnantian glaciation Full size image

With a gradually increasing areal extent of bioturbation across the Palaeozoic (Scenario 2; Fig. 5h–n), organic matter burial gradually decreases, while sedimentary phosphate retention gradually increases, both leading to an increase in ocean anoxia (Fig. 5k, l). While the associated protracted decrease in δ13C carb values is not entirely at odds with the geochemical record, the gradual increase of ocean anoxia across the Palaeozoic does not fully agree with geochemical proxy evidence for an increase in ocean anoxia in the early Cambrian31. Furthermore, the absence of a transient increase in δ34S SO4 also disagrees with the geochemical record (Fig. 5i).

When the model simulations incorporate a strong biogeochemical response to shallow bioturbation in the early Cambrian (Scenario 3; Fig. 5o–u), the emergence of bioturbation results in significantly enhanced oxidation of marine organic carbon and pyrite, as well as benthic phosphate retention, which limits oceanic primary production. The initial decrease in marine organic carbon burial at 520 Ma (marked by the drop in δ13C carb from 2‰ to 0‰; Fig. 5r) is accompanied by an increase in ocean anoxia (Fig. 5s), which then provides a negative feedback and limits bioturbation. The resulting reduction in marine organic carbon burial results in increased atmospheric pCO 2 (by >1000 ppm; Fig. 5o), and suggests that the evolution of burrowing organisms in the ocean could have triggered significant climate warming, consistent with the ‘greenhouse’ climate invoked for the Cambrian and early Ordovician, as indicated by elevated sea levels52 and oxygen isotope systematics53. We note that the direct effect of the organic carbon cycle on CO 2 levels is not as widely discussed in the literature as the changes resulting from silicate weathering, but that it is an important part of the coupled C and O cycles54.

The increase in anoxia induced by bioturbation leads to an increase in δ34S SO4 (~30‰; Fig. 5p). We do not find, as previously suggested17, a rise in marine sulphate concentrations coincident with the evolution of bioturbation (Supplementary Fig. 2), because of our implementation of gypsum (CaSO 4 ) burial. In COPSE, gypsum burial scales with the oceanic sulphate concentration, while gypsum burial was previously considered important only at higher sulphate concentrations17. Given the uncertainty over gypsum formation, it is difficult to make any strong conclusions regarding this effect. In general, the effect of the evolution of bioturbation on ocean sulphate concentration remains uncertain. Eventually, the effects of the evolution of bioturbation are largely reversed by the rise of early land plants in the late Ordovician, which increases net carbon burial and phosphate supply, as well as drawing down atmospheric CO 2 . Assumed high Palaeozoic erosion rates contribute to low δ13C carb in the model runs54, but it is unlikely that the rapid drop in δ13C carb at the end of the Cambrian explosion and the subsequent rise at the end of the Ordovician can be driven solely by rapid erosion and sedimentation, which continues into the Devonian. Overall, Scenario 3 brings the COPSE model predictions reasonably close to the available data from the geological record.

Caution is of course required, as bioturbation is likely not the only driver for variation in the isotope data, especially for variations on timescales <1 million years, which the model cannot capture. However, these limitations do not invalidate the hypothesis testing here, which shows that Scenario 3 (large geochemical effects due to shallow burrowing) can reproduce the broad changes seen in δ13C carb , δ34S SO4 and ocean anoxia, whilst Scenarios 1 and 2 (more protracted global biogeochemical responses) produce discrepant predictions for δ13C carb or δ34S SO4 . Reversing these conclusions would require a very specific, and unlikely, set of additional forcing factors (see Supplementary Note 2 for further discussion). Nevertheless, many additional mechanisms were likely at work over the studied time interval. For example, nutrient stress and limited primary production were undoubtedly a factor in explaining the low-oxygen conditions in the early Palaeozoic. The COPSE model attempts to take into account such mechanisms: for example it includes process-based ocean nutrient cycles, including a representation of phosphorus removal with organic matter, calcium, and iron species and how these sinks may have varied with time26. However, as with all box models, these processes are greatly simplified, and apparent mismatches with geochemical data suggest that there are still processes missing or poorly represented (e.g., the P cycle). The COPSE model predictions are improved by our additional representation of bioturbation, suggesting that it was a major driver of biogeochemical cycling during the Palaeozoic.

Synthesis

By recognising that moderate levels of shallow bioturbation have a large impact on sediment geochemistry (Fig. 4a), our model results are able to broadly reconcile the bioturbation record (Fig. 1e) in relation to various proxies in the geochemical record (Fig. 5). We propose a significant bioturbation-driven step change in environmental conditions and geochemical cycling in the early Cambrian, well before benthic fauna reached their full capacity in terms of sediment reworking. This appears to have resulted in a ~100 Myr period of prevalent ocean anoxia and greenhouse climate conditions that is consistent with the available geological evidence. This transition period between the Cambrian explosion and the Great Ordovician Biodiversification Event (GOBE) was also marked by sizable fluctuations in δ13C carb values and variable ocean redox conditions (Fig. 1) on shorter timescales, alongside repetitive extinction and recovery events that sustained a radiation plateau55 which was eventually followed by the GOBE56.