Possible links between extreme oxygen perturbations and the Cambrian radiation of animals.

The role of oxygen as a driver for early animal evolution is widely debated. During the Cambrian explosion, episodic radiations of major animal phyla occurred coincident with repeated carbon isotope fluctuations. However, the driver of these isotope fluctuations and potential links to environmental oxygenation are unclear. Here, we report high-resolution carbon and sulphur isotope data for marine carbonates from the southeastern Siberian Platform that document the canonical explosive phase of the Cambrian radiation from ~524 to ~514 Myr ago. These analyses demonstrate a strong positive covariation between carbonate δ13C and carbonate-associated sulphate δ34S through five isotope cycles. Biogeochemical modelling suggests that this isotopic coupling reflects periodic oscillations in atmospheric O2 and the extent of shallow ocean oxygenation. Episodic maxima in the biodiversity of animal phyla directly coincided with these extreme oxygen perturbations. Conversely, the subsequent Botoman-Toyonian animal extinction events (~514 to ~512 Myr ago) coincided with decoupled isotope records that suggest a shrinking marine sulphate reservoir and expanded shallow marine anoxia. We suggest that fluctuations in oxygen availability in the shallow marine realm exerted a primary control on the timing and tempo of biodiversity radiations at a crucial phase in the early history of animal life.

The early Cambrian witnessed a dramatic diversification of animal body plans and behaviours1, as well as between-species interactions and palaeocommunity innovations2,3, ultimately leading to modern animal ecosystems. Ocean oxygenation is a commonly invoked environmental pre-requisite4–6. However, some recent studies suggest that despite probable low-oxygen conditions, the oceans exceeded requisite oxygen thresholds for simple animals, such as sponges, well before the Cambrian Period7,8. Many of the new animal body plans and lifestyles that appeared during the early Cambrian were associated with considerably higher oxygen demands9,10. Fluctuations in the maximum dissolved oxygen content of surface waters, or the extent of shallow ocean oxygenation, could therefore have played an important role in regulating the pattern of Cambrian radiations. This brings into question the role of oxygen in early animal evolution, which is exacerbated by a lack of convincing evidence for a direct link between Earth’s oxygenation history and early Cambrian bio-radiations and extinctions11.

High-resolution records of the sulphur and carbon cycles, when considered in the context of the fossil record may, however, afford an opportunity to resolve potential environmental controls on early animal evolution. The marine biogeochemical sulphur and carbon cycles interconnect via their respective redox-sensitive reservoirs and fluxes. Both elements have a single, large oxidised oceanic reservoir (dissolved sulphate and inorganic carbon), the isotopic composition of which is governed by isotope fractionation during microbially-mediated reduction to sulphide (ultimately preserved as pyrite) and organic carbon. Burial of these reduced species represents the two main net sources of oxygen to the surface environment12–14, and also imprints on both the seawater sulphate sulphur isotope (δ34S, as recorded by carbonate-associated sulphate) and carbon isotope (δ13C, as recorded in carbonate) records, allowing redox changes in the surface environment to be traced through geologic time.

Here we present paired carbon and sulphur isotope data from lower Cambrian marine carbonates from the southeastern Siberian Platform. These data provide a continuous, high-resolution record from Cambrian Stage 2 through to Stage 4 (~524–512 Myr ago; Fig. 1), and allow a direct assessment of potential links between ocean redox variability, atmospheric oxygenation, and the major biological events of the early Cambrian.

Fig. 1

Carbonate carbon and carbonate-associated sulphate sulphur isotope records from Cambrian Stage 2 to Stage 4 of Siberian Aldan-Lena rivers sections.

Regional stage subdivisions are shown next to the global subdivision plan for comparison15 (F.: Fortunian Stage; N.–D.: Nemakit–Daldynian Stage; TST: Transgressive System Tract36; Fm.: Formation; A.: archaeocyaths; SSFs: small shelly fossils). Names for the positive δ13C peaks (III, IV, V, VI, VII) are consistent with those of previously suggested δ13C curves15. FAD: first appearance datum.

Carbon and sulphur isotope systematics

Carbonate δ13C and carbonate-associated sulphate δ34S analyses (see Methods) were primarily performed on well-preserved micritic limestone samples collected from sections along the Aldan and Lena rivers in Siberia. These sections archive a continuous and highly fossiliferous sedimentary record from a shallow, open ocean carbonate platform, and preserve over half of all fossil diversity currently known from the Cambrian radiation interval worldwide, thus providing a unique window into early Cambrian shallow marine ecosystems (see Supplementary Information for geological and palaeontological context, sample details, diagenesis evaluation and all data).

Our carbon isotope data record five cycles through Stage 2 and Stage 3 of the lower Cambrian. Positive excursions are labelled here as III to VII (Fig. 1), consistent with previous studies of the Siberian Platform15,16, but these excursions are also found elsewhere17,18. The new sulphur isotope data range from +16‰ to +36‰, demonstrating that seawater sulphate δ34S values fell from a peak (~40–45‰) during the late Ediacaran19,20 to lower values by the early Cambrian. Significantly, these data also demonstrate for the first time that oceanic sulphate δ34S values varied across five cycles that directly correlate with excursions in seawater δ13C (Fig. 1; see Supplementary Table S1 for statistical correlation parameters). In sharp contrast to the coupled δ13C-δ34S trends during the Cambrian stages 2-3, the δ34S trend across the early Cambrian Stage 4 Botoman–Toyonian extinctions (BTE; the first animal mass extinction of the Phanerozoic Eon)21,22 is characterised by rapid fluctuations of large magnitude that are decoupled from the carbon isotope record (Fig. 1).

Over long timescales the excess oxidant generated by increased organic carbon burial (as indicated by higher carbonate δ13C) may be balanced by reduced rates of pyrite burial (lower seawater sulphate δ34S), and vice-versa, which results in relatively stable atmospheric oxygen levels and an inverse relationship between the first-order global seawater δ13C and δ34S records23,24. However, the positive correlations we observe between δ13C and δ34S in Cambrian stages 2–3 likely reflect higher rates of both organic carbon and pyrite sulphur burial, which may have been associated with large distinct pulses in atmospheric oxygenation, as previously suggested for the late Cambrian SPICE (Steptoean Positive C-isotope Excursion) event14.

The rate of change of seawater sulphate sulphur isotope ratios allows us to estimate marine sulphate concentrations through this interval. Using the ‘rate method’ model25,26 (see Methods for model details), and taking the average values of the lower end of the data envelopes shown in Supplementary Fig. S2, an upper estimate can be obtained for marine sulphate of ~1.0–6.6 mM for the interval from ~524 to ~514 Myr ago, followed by ~0.4–1.4 mM for ~514 to 512 Myr ago. Estimates for the earlier interval are broadly consistent with previously modelled estimates of ~5–10 mM25 and with fluid inclusion-based estimates of ~4.5–11 mM for the early Cambrian27, but trend toward lower values. Thus, the early Cambrian ocean was characterized by a relative paucity of sulphate, when compared with the modern ocean (~28–29 mM). Our data document a significant drawdown of more than half of the sulphate pool during early Cambrian Stage 4 (~514–512 Myr ago), coincident with the BTE.

Environmental oxygenation and animal radiations

The covariant behaviour of the carbon and sulphur isotope systems during Cambrian Stage 2 to late Stage 3 can be explained by coupled burial of pyrite and organic carbon in marine sediments under highly productive, anoxic conditions23,28,29. Such conditions result in enhanced preservation and burial of organic carbon, and simultaneously enhance microbial sulphate reduction (MSR), leading to a high pyrite burial flux. Since pyrite and organic carbon are enriched in the lighter isotopes (32S and 12C respectively), elevated burial fluxes on a global scale would drive the positive excursions in seawater sulphate δ34S and inorganic δ13C.

A biogeochemical box model12,30,31 (see Methods for model details) was applied to test whether measured trends in S isotopes can be reproduced from the coupled burial of sulphur (as pyrite) and carbon (as organic carbon). The model infers the rate of organic carbon burial using the δ13C record and an isotopic mass balance, while the rate of pyrite burial is calculated by assuming a linear relationship with organic carbon burial, allowing prediction of δ34S values. Results (Fig. 2c) show that both the amplitude of positive sulphur isotope excursions and their long-term trend from ~524–514 Ma can be replicated in this way. The model assumes that the isotopic composition of carbon and sulphur inputs (δ13Cin, δ34Sin), and the background carbon and sulphur cycle input fluxes through weathering and metamorphism remained constant. Variations in these processes may help to explain the slight drift of the baseline δ34S in model average predictions when compared to the observed δ34S data. The shaded areas in Fig. 2 show the result of varying δ13Cin between -5‰ and -8‰, allowing the model to encompass most of the data. Our model requires a low concentration of sulphate in seawater (best-fit shown is 1 mM), in order to match the rate and amplitude of δ34S variations, consistent with the lower end of maximum estimates derived from the ‘rate method’ model.

Fig. 2

Carbon and sulphur cycle model output.

a. This model takes measured δ13C values as an input parameter. b. Burial rates of organic carbon (Corg) are inferred from isotope mass balance and δ13C record, and burial rates of pyrite are assumed to be controlled by modelled organic matter availability. c. Comparison between analysed δ34S data (green curve) and simulated seawater sulphate δ34S values (pink); Dashed part of the green curve shows the sampling gap. d. Variations in modelled net oxygen production. For all plots, the uncertainty window represents an alteration of the δ13C values of carbon inputs between -5‰ and -8‰.

The coupled carbon and sulphur isotope swings show repeated cycles of approximately 0.5–2 Myrs duration that reflect cyclical changes in the burial rates of organic carbon and pyrite, which may have been induced by episodic expansion of bottom-water anoxia/euxinia on the deeper portions of continental shelves and slopes. Ultimately, coupled burial of both reduced species in marine sediments results in the release of oxygen and other marine oxidants14. Each rising limb and the peak of the positive isotope swing thus represents enhanced net oxygen production and a pulse of atmospheric oxygen, which initially increased the extent of oxygenated waters and/or the maximum dissolved O2 in the shallower realm. Subsequently, increased ventilation of the deep ocean would have resulted in a reduced flux of reductant (organic carbon and pyrite) to seafloor sediments32. This acted to decrease the net oxidant flux, which ultimately buffered against further oxygenation. Furthermore, positive feedbacks between ocean ventilation and phosphorus retention in sediments33 may have driven rapid bottom-water oxygenation, and in this case the decrease in the net oxidant flux may be substantial, leading to a re-establishment of anoxia, and potentially giving rise to the repetitive isotope cycles33.

Alternatively, isotope cyclicity might be driven by orbital forcing via climatic impacts on weathering, similar to the ~1–2 Ma “third-order” eustatic sequences of the Mesozoic and Cenozoic Eras34. However, neither the timing, duration and frequency of early Cambrian third-order sea-level fluctuations18,35, nor regional sequence stratigraphy data from Siberia36 (Supplementary Table S3), appear to match the isotope cycles identified in this study. Similarly, an erosional driver37 for the observed isotope cycles is incompatible with their combined high amplitude and frequency, which would require very large (~3-5 fold37) changes in global erosion over geologically-short timescales. Furthermore, an erosional driver is not supported by contemporaneous changes in seawater 87Sr/86Sr38. Fluctuations in oxygen minimum zone depth39,40, alongside biological feedbacks such as enhanced diurnal vertical migration via increased expansion of metazoan mobility41, may also have contributed to the perturbations in shallow ocean oxygenation.

To summarize, our model indicates the potential for large variations in the net atmospheric oxygen production flux (±50% around the baseline value; Fig. 2d). We propose that periods of rising δ13C represent enhanced burial of reductants under anoxic bottom-water conditions and atmospheric oxygenations, whereas the falling limbs record the decrease of reductant burial under a more widely oxygenated deep ocean. A more direct estimate of oxygen production rates can be made within our model by treating both δ13C and δ34S as input parameters, thus inferring rates of organic carbon and pyrite burial, respectively, for the time points where we have input information for δ34S. These estimates are shown in Fig. 3 and are similar in magnitude to those of the carbon-only model, which is to be expected as the carbon-only model produced a reasonable fit to the δ34S data.

Fig. 3

Animal diversity, biological events and their correlation to the isotope records and oxygenation pattern across Cambrian stages 2-4.

Global oxygen production is inferred from isotope mass balance modelling, using inputs of δ13C only (light shade), or δ13C and δ34S (dark shade). Archaeocyathan species (blue line) and total animal species (green line) diversity records are expressed as the mean number of species per sampling unit (grey box) in Siberia; OP: oxygenation pulse; BH: biodiversity high; F.: Fortunian Stage; N.–D.: Nemakit–Daldynian Stage. FAD: first appearance datum.

One direct impact of pulses in atmospheric oxygenation during the early Cambrian was episodic oxygenation of marginal shallow marine environments. Shallow carbonate platforms, such as the Aldan-Lena rivers region, evidence relatively high animal origination rates and biodiversity42,43. Within shallow ocean ecosystems, biogenic reefs serve as critical evolutionary cradles and net sources of marine biodiversity44. Comparing the isotopic cycles and estimated oxygen production curves with species diversity curves for the Siberian Platform (see Supplementary information for full palaeontological data), oxygenation pulses (III, IV, V, VI, VII) generally coincided with regional biodiversity highs in either reef-building archaeocyathan or total animal species (Fig. 3). Although no significant total animal biodiversity high was associated with oxygenation pulse IV, the number of archaeocyathan species increased dramatically by ~60%. Moreover, the rising limb of isotope excursion IV coincided with the first emergence of trilobites, bivalved arthropods, and stenothecoids possessing relatively thick biomineralised skeletons, as well as a geographic expansion of possible burrowing filter-feeding arthropods over the Siberian Platform, as recorded by the appearance of Thalassinoides-type trace fossils36. A significant increase in the inter-habitat (beta)-diversity of reefal palaeocommunities was also restricted to the IV interval in the Aldan-Lena rivers region2, reflecting a differentiation of species between assemblages, and thus ecological diversification within the shallow marine environment.

On a global scale, positive isotope excursion V appears to coincide with major radiations of large predatory arthropods and radiodonts, increased durophagy, and the first appearance of pelagic motile deuterostomes, evidenced by the Chengjiang biota and similar faunas11,45. Similarly, excursion VII coincided with a global radiation of echinoderms and archaeocyaths. The latter is revealed by the inter-regional (gamma)-diversity peak reflecting formation of numerous isolated faunal provinces2. By contrast, minor extinction events here and elsewhere appear to be associated with the negative excursions11,18,46. In the deeper ocean setting of northern Siberia and South China, multi-proxy analyses reveal broadly similar oceanic redox fluctuations4,47–49, which coincide with the positive carbon isotope excursions in the early Cambrian16,46,50. These episodic redox oscillations, evident from the δ13C record and, in places, as δ13C and δ34S covariance16 (also see Supplementary Fig. S6 for δ13C-δ34S covariance from the Cambrian Stage 2 ZHUjiaqing Carbon isotope Excursion (ZHUCE) in the Xiaotan section, South China), suggest that these coupled isotope excursions record a global phenomenon. We therefore propose that perturbations to shallow ocean oxygen budgets were driven by fluctuations in atmospheric oxygen. High oxygen levels would have suited various newly evolved animal body plans and lifestyles, and so oxygen fluctuations likely resulted in episodic expansions/contractions of the habitable zone within shallow ocean ecosystems. This shallow ocean oxygen control is likely reflected in contemporaneous fluctuations of animal origination and speciation rates, and thus possibly regulated the global radiation patterns of early Cambrian animals.

Expanded shallow ocean anoxia and sulphate reduction across the BTE

In contrast to the coupling of carbon and sulphur isotopes during Cambrian stages 2-3, the decoupled δ13C-δ34S records and unsystematic temporal fluctuations in δ34S values observed across the BTE (Fig. 1) appear to reflect a significant and persistent decline in oceanic sulphate concentration (Supplementary Fig. S2). At reduced marine residence times, δ34S is more responsive to perturbations to the sulphur cycle. A fall in seawater sulphate concentration is generally attributed to enhanced evaporite deposition or widespread anoxia, and indeed, there are a number of thick evaporite deposits in the global rock record during this interval51,52. However, these evaporites are restricted to the innermost isolated basins of the Siberian Platform and the Australian part of Eastern Gondwana, and their stratigraphic distribution does not correlate with the interval of low sulphate inferred for the BTE. This suggests that anoxic/euxinic conditions likely prevailed in the shallow marine realm at this time. The expansion of shallow ocean anoxia is consistent with an observed accumulation of over ~750,000 km2 of black organic-rich carbonate-rich sediments (comprising bituminous limestone, chert and argillaceous calcareous sapropelic shale) in the Sinsk Formation across the Siberian Platform, as well as enrichments in pyrite, V, As, Cr, Cu and Ni, and the presence of abundant biomarkers indicative of anaerobic bacteria as a major source of organic matter21,53. Such phenomena have previously been linked to shoaling of oxygen-depleted waters during a major marine transgression21,36, which has been suggested as the cause of the major extinction pulse of the BTE (Sinsk event; Fig. 1). Thus, while bottom-water anoxia on the deeper portions of continental shelves and slopes may have contributed to the episodic burial of reductant and oxygenation of the atmosphere and shallow oceans in Cambrian stages 2-3, shoaling of anoxic waters in Cambrian Stage 4 may have driven a mass extinction, and therefore a reduction in primary productivity and overall reductant burial.

Implications for early animal diversification

Oxygenation of the early Cambrian shallow marine environment can be inferred from the coupled behaviour of the carbon and sulphur cycles. Episodic shallow ocean oxygenation corresponded to pulses of animal diversification, and so provides a plausible environmental explanation for the step-wise nature of the Cambrian radiation of animals. In the modern and ancient oceans, well-oxygenated waters are generally associated with larger body sizes, higher diversity, advanced skeletal biomineralization, and increased motility and carnivory9,10,54,55. Pulses of shallow ocean oxygenation in the early Cambrian likely expanded the global proportion of habitable marginal ocean to provide new ecological opportunities and biodiversity cradles. Similarly, the extended radiation of the Great Ordovician Biodiversification Event (~490–450 Ma) also appears to have been facilitated by pulses in atmospheric oxygenation56. A prolonged pause in biological diversification, which lasted over 20 million years and was associated with recurring extinctions (BTE, SPICE-trilobite extinctions18), occurred between these two major diversification events. Environmental stress caused by the persistent development of oxygen-deficient conditions in shallow marine realms due to low net atmospheric oxygen production57 is likely to have been a major contributing factor. Thus, the global extent of well-oxygenated shallow ocean habitats during the early Paleozoic, as well as the maximum dissolved oxygen content of surface waters, played a vital role in regulating the emergence and radiation of early animal life.

Methods

Carbonate-associated sulphate (CAS) extraction and δ34S analysis

Well preserved carbonate samples composed primarily of micrite were targeted for CAS extraction. Where that was not possible, few samples were selected with sparitic or dolomitic textures. Large blocks (>200 g) of carbonate rocks were cut and polished under running water to trim weathered surfaces prior to powdering. Blocks were then cut into small chips using a water-cooled, diamond tipped bench circular saw. Rock chips were ground to a fine powder (flour-like consistency, <10 µm) using a Retsch® Agate Mortar grinder. We applied a high-fidelity miniaturized CAS extraction protocol, which is an extension of two published approaches58,59. The protocol was established following tests involving twelve consecutive leaching steps on five carbonate samples from different stratigraphic horizons of the Aldan-Lena river sections and three samples from the Ediacaran Nama Group19. Approximately 10 g of the fine powder for each sample was leached in 40 ml of 10% NaCl solution for 24 hours to remove the non-CAS sulphur-bearing compounds and easily soluble sulphate. During leaching, samples were constantly agitated using a roller shaker at room temperature. Residues were rinsed in ultrapure water three times between each leach and five times after the final leach. After each leach, the leachate was retained, and the presence of sulphate was tested by adding saturated barium chloride solution and allowing three days to precipitate barite. As illustrated in Supplementary Fig. S3, the amount of sulphur removed during sequential NaCl leaching of test samples exhibited a sharp decline through multiple NaCl leaches and reached blank levels in the 3rd or 4th leachates, suggesting five leaches is sufficient for complete removal of all soluble sulphur-bearing constituents from ~10 g of carbonate powder. All five-times pre-leached carbonate samples were treated with 6 M HCl, which was added in calculated aliquots based on total HCl-leachable carbonate content. This step was completed within 30 minutes to minimise the potential for pyrite oxidation during dissolution. The insoluble residue was separated from the solution by centrifugation in 50 ml tubes followed by filtration through VWR® 0.2 µm Polypropylene membrane syringe filters. Saturated barium chloride solution was then added to the filtered solution and left to precipitate within the housing of a sealed tube over three days at room temperature. Where no visible precipitate was observed after 24 h, 2 mg isotopic-grade sulphur-free quartz powder was added, which served as an inert medium onto which barium sulphate could precipitate59. Each sample was centrifuged, and the supernatant replaced with ultrapure water repeatedly until the pH attained neutral values. Washed samples were then dried prior to isotope analysis. 34S/32S analysis of barium sulphate precipitates was undertaken using an Elementar® Pyrocube elemental analyzer linked to an Isoprime® 100 mass spectrometer operated in continuous flow mode at the Lancaster Environment Centre, Lancaster University. Pellets of BaSO4, resulting from sulphate extraction with or without the addition of quartz powder, were combusted in tin capsules in the presence of excess vanadium pentoxide (V2O5) at 1030°C to yield SO2 for the determination of δ34S. All samples and standards were matrix matched, and values were corrected against VCDT using within-run analyses of international standards NBS-127 and SO5 (assuming δ34S values of +20.3‰ and +0.49‰, respectively). Within-run standard replication was below 0.3‰ (1sd). Procedural standard solutions of calcium sulphate precipitated as barium sulphate were used to test the integrity of the method59. These yielded δ34S values of +2.7‰ (±0.3‰, 1sd, n=12) compared to values of +3.0‰ (±0.3‰, 1sd, n=13) for analysis of the raw calcium sulphate powder. Blank contamination associated with δ34S determination was zero.

CAS concentrations and sulphur content in NaCl leached solution

The concentration of CAS and sulphur content in each leaching step was measured in aliquots of filtered solution using a Varian ® 720 Inductively Coupled Plasma Optical Emission Spectrometer (ICP-OES) at the London Geochemistry and Isotope Centre (LOGIC), University College London. Wavelength 182.5 nm was selected to minimise interference with calcium ions, and analysis was conducted using the N2-purging polyboost function to avoid oxygen interference in the system.

Carbonate carbon and oxygen isotopes

Micritic limestone was targeted for δ13C analysis. Where that was not possible, we selected a few sparitic or dolomitised samples and fossiliferous samples with skeletal components known to secrete low-Mg calcite. About 20 mg of powder drilled from a rock chip was analyzed for stable C and O isotopes. Limestone samples were reacted with 100% H3PO4 at 25°C for more than 12 h, and dolostone samples were reacted with 100% H3PO4 at 50°C for more than 24 h. Prepared gas samples were analysed for 13C/12C and 18O/16O using the Chinese national standard, an Ordovician carbonate from a site near Beijing (reference number GBW04405: δ13C= 0.57 ± 0.03‰ VPDB; δ18O= -8.49 ± 0.13‰ VPDB). The analyses were performed using the Finnigan® MAT 253 mass spectrometers at the Nanjing Institute of Geology and Palaeontology, Chinese Academy of Sciences.

Elemental analysis

For concentrations of diagenesis-diagnostic elements, including Ca, Mg, Mn, and Sr, an aliquot of approximately 50 mg of power was micro-drilled from a rock chip and dissolved with excess 6 M hydrochloric acid at room temperature for 12 h. The concentration of acid used here is identical to the concentration used during CAS extraction. The reaction was facilitated using an ultrasonic bath and roller shaker. After centrifugation, aliquots of the supernatant were analysed for elemental concentration using a Varian® 720 ICP-OES at University College London. Solution standards of certified reference materials, SRM1c (argillaceous limestone) and SRM120b (Florida phosphate rock), were run at the start of the analyses along with a blank to monitor the accuracy of the bulk elemental analysis. Laboratory control solution standards were also run after every batch of 20 samples to monitor drift and precision. Analytical precision for elemental concentrations was generally better than 5%.

‘Rate method’ model

Maximum seawater sulphate concentrations are calculated using the modified ‘rate method’25,26. The model was constructed based on the observed rate of change in seawater sulphate (carbonate-associated sulphate) δ34S, fractionation between oxidized (sulphate) and reduced sulphur (pyrite) reservoirs and equation (1) that connects the two parameters, where Fx represents the input and output fluxes, Δ34Sx represents isotopic difference of δ34S values between fluxes (Q = total input flux of sulphur, SUL = seawater sulphate, PY = pyrite burial, SW = seawater/sulphate deposition) and MSW represents the mass of sulphate in the ocean. The maximum rates of δ34S change are attained when sulphur input flux to the ocean approaches zero (FQ = 0), and the standing oceanic sulphate reservoir is removed as pyrite. Equation (1) is then transformed to equation (2) to calculate the size of seawater sulphate reservoir. Because the observed rates of seawater sulphate δ34S change in a normal marine environment should never exceed the theoretical maximum rates of change (dδ34S/dt), the calculation of MSW using equation (2) should provide the maximum estimate of seawater sulphate concentration. The definition of FPY, Δ34SSUL-PY, and unit-conversion constants (gram to mM) are consistent with the values applied for the long-term secular variation of seawater sulphate concentration25. FPY = 4 × 1013 g yr-1 is suggested for a normal marine environment. Δ34SSUL-PY = 35‰ is suggested for the fractionation during MSR. The variation of seawater sulphate concentration ([SO42-]) between ~524 Myr ago and ~512 Myr ago is represented based on a point-to-point calculation (Supplementary Fig. S2). Because the sampling density between δ34S values is generally below 0.1 Myr (Supplementary Table S3), this study uses a 0.1 Myr gridded data smoothing curve (red line in Supplementary Fig. S2) to represent the best estimate of seawater [SO42-]. Besides, the maximum concentration for an individual point could be under or overestimated due to fluctuations and anomalies in the rate of δ34S changes. To overcome this bias, the resulting [SO42-] data are binned into 0.5 Myr bands. The lower envelope (black dotted line in Supplementary Fig. S2) of the [SO42-] red curve, which links the lowest value for each band, is expected to represent the maximum rates of δ34S change and thus the theoretical estimate of maximum seawater sulphate concentration through time.

Coupled carbon and sulphur cycle model

A simple model of the global carbon and sulphur cycles was applied to explore the proposed mechanisms for isotopic variations in the system. This follows the work of Garrels and Lerman30, Berner12 and Bergman et al.31. The model calculates the global rate of organic carbon burial using isotope mass balance, and then attempts to predict the operation of the sulphur system based on the supply of organic matter. Supplementary Fig. S4 shows the model processes as a diagram; Supplementary Table S2 shows the model flux and parameter values. The model estimates long-term fluxes between the ocean and sediments for both carbon and sulphur. Carbon is modelled as CO2 in the atmosphere and ocean (A), and will be buried either as organic carbon (G) or carbonate (C). Similarly, sulphur can exist as oceanic sulphate (S), and will be buried as pyrite (PYR) or gypsum (GYP). Weathering (and metamorphism) constitutes the return flux from the sediments to the ocean and atmosphere. We set the weathering inputs to constant values in line with previous models12,31. We allow for around half of present total organic carbon burial (and weathering) due to the absence of land plants, and an enhanced burial flux of pyrite sulphur due to anoxia. The weathering rate of gypsum is held constant, but the burial rate is adjusted so that the model maintains a constant sulphate concentration. Due to the relatively short model timeframe relative to the residence times of the vast sedimentary reservoirs, these reservoirs are assumed to have a fixed isotopic composition and are assumed not to vary in size. The ocean and atmosphere reservoirs are allowed to vary in size and isotopic composition. Organic carbon burial is calculated via isotope mass balance12,30, which uses the total carbon input fluxes and isotopic composition of seawater (δA) to calculate the required burial rate of isotopically depleted organic carbon (equation (3)): It is assumed that pyrite burial is governed by the supply rate of organic carbon to microbial sulphate reducers, and therefore scales with the burial rate of organic carbon (equation (4)). The proportionality constant (0.5) is chosen to balance pyrite weathering. Variation in the ocean and atmosphere carbon is calculated as: Variation in ocean sulphate is calculated as: Variation in the isotopic composition of ocean sulphate is calculated as: Net oxygen production flux is calculated from the burial rate of organic carbon and pyrite: The model is solved in MATLAB using the ODE (Ordinary Differential Equation) suite. The model broadly reproduces the duration and magnitude of fluctuations in δ34S (Fig. 2c). It also predicts similar fluctuations in oxygen production (Fig. 2d). The model does not calculate the concentration of oxygen in the atmosphere and ocean, and all fluxes are assumed to be oxygen-independent. More detailed modelling, which takes into account the variation in oxygen sinks, is required to analyse the overall long-term trends in atmospheric oxygen levels.

An alternative version of the model is run in Fig. 3 that estimates pyrite burial rates directly from the δ34S record. In this version, equation (4) is replaced by equation (9), and equation (7) is not required.

Total marine animal species diversity

Supplementary Table S4 shows the distribution and diversity of total and individual animal species of Cambrian stages 2-4 of the Siberian Platform. This dataset is an upgrade of a previously published version21 (see supplementary information for detailed description and source of data). Siberian biozones (archaeocyathids/trilobite) are selected as the sampling units for diversity data collection. The finalised animal diversity record is generated by plotting total species diversity against sampling units (grey boxes in Fig. 3).