Multiple sulfur isotope evidence for massive oceanic sulfate depletion in the aftermath of Snowball Earth

The terminal Neoproterozoic Era (850–542 Ma) is characterized by the most pronounced positive sulfur isotope (34S/32S) excursions in Earth's history, with strong variability and maximum values averaging δ34S∼+38‰. These excursions have been mostly interpreted in the framework of steady-state models, in which ocean sulfate concentrations do not fluctuate (that is, sulfate input equals sulfate output). Such models imply a large pyrite burial increase together with a dramatic fluctuation in the isotope composition of marine sulfate inputs, and/or a change in microbial sulfur metabolisms. Here, using multiple sulfur isotopes (33S/32S, 34S/32S and 36S/32S ratios) of carbonate-associated sulfate, we demonstrate that the steady-state assumption does not hold in the aftermath of the Marinoan Snowball Earth glaciation. The data attest instead to the most impressive event of oceanic sulfate drawdown in Earth's history, driven by an increased pyrite burial, which may have contributed to the Neoproterozoic oxygenation of the oceans and atmosphere.

S ulfate represents one of the most important metabolic electron acceptors on the planet, and constraining its biogeochemical cycle is crucial for understanding the longterm redox evolution of the oceans and atmosphere [1][2][3][4] . Classically, the sedimentary record of sulfur cycling is probed through the sulfur isotopic composition of sedimentary pyrite, d 34 S pyr , and carbonate-associated sulfate (CAS), d 34 S CAS (where d 34 S ¼ 10 3 Â ( 34 S/ 32 S sample / 34 S/ 32 S Canyon Diablo Troilite À 1), both of which may record fluctuations in past marine sulfate isotope composition 5 . Variations in oceanic sulfate isotope composition (d 34 S SO4 ) are usually interpreted under steady-state assumptions [6][7][8][9][10] , whereby the oceanic sulfate content (M SO4 ), here written for 32 S, remains constant: F in being the input flux related to weathering and F out the output flux related to pyrite and evaporite burial. Long-term variations in d 34 S SO4 therefore depend on the fraction of sulfur buried as pyrite relative to evaporite (f pyr, ranging from 0 to 1) as well as on the isotopic composition of sulfate delivered to the ocean (d 34 S in ), according to the conservative isotopic mass equation: where D 34 S SO4-pyr ¼ d 34 S SO4 À d 34 S pyr is the difference between the average d 34 S SO4 of evaporites and/or CAS and the average d 34 S pyr of sedimentary pyrite at a given time. Considering a modern d 34 S in of þ 3%, a D 34 S SO4-pyr of 40%, and a modern seawater d 34 S SO4 of þ 21%, the resulting present-day f pyr is close to 0. 45 (ref. 11).
In a steady-state framework (equation (2)) and assuming modern d 34 S in and D 34 S SO4-pyr values, the strong positive d 34 S CAS recorded in Neoproterozoic 6,7,9,10,12,13 , sediments would only hold for f pyr values above unity which are therefore inconsistent. A concomitant increase of both f pyr and D 34 S SO4-pyr for such excursions is thus necessary. Similarly, assuming that D 34 S SO4-pyr is constant, the high d 34 S CAS values (accompanied by high d 34 S pyr ) need to be accounted by a concomitant increase in both f pyr and d 34 S in . For example, the d 34 S CAS positive Ara anomaly (B545 Ma) 10 in Oman requires anomalously high f pyr (0.9) and d 34 S in (15%) values for which a geological driver remains to be identified. Others have instead observed an increase in d 34 S CAS values not accompanied by a change in d 34 S pyr values, resulting in a D 34 S SO4-pyr increase that may reflect the advent of sulfur disproportionation (SD) metabolism 7,8 . An alternative possibility is that the steady-state model itself does not hold for the Neoproterozoic. This hypothesis has already been envisaged 13,14 but has never really been tested so far because it increases further the number of unconstrained parameters compared with the steady-state model (for example, d 34 S in , D 34 S SO4-pyr , F in , F out and M SO4 ).
Here we use multiple sulfur isotopes ( 33 S/ 32 S, 34 S/ 32 S and 36 S/ 32 S ratios) of CAS and pyrite to investigate dynamic models for the Neoproterozoic sulfate reservoir evolution in the aftermath of the Marinoan glaciation. The results show that the steady-state assumption does not hold in the aftermath of the Marinoan Snowball Earth glaciation and attest to an impressive event of oceanic sulfate drawdown.

Results
Stratigraphy and age constraints. The studied sedimentary sequence corresponds to the Mirassol d'Oeste and Guia formations (Araras group, Central Brazil) and starts with a typical cap dolostone (B635 Ma, refs 15-17) directly overlying Marinoan glacial deposits. We focus primarily on CAS as its isotopic composition is generally considered to reflect a seawater signal 5 , provided that diagenesis and recrystallization are limited. Samples were selected based on previous petrographic and geochemical results 17 to avoid as much as possible diagenetic overprints. Selected samples were then carefully washed to avoid contamination (Methods section). It has been shown that secondary pyrite oxidation to sulfate and atmospheric sulfur contamination may lead to lower d 34 S CAS values 18,19 . The consistent isotopic pattern displayed by the generated data set (see below) suggests that contamination was not significant in our samples.
Three slightly overlapping sections were sampled in freshly exposed quarries (Terconi, Camil and Carmelo quarries) along a basinward profile, from the innershelf to the outershelf of the Araras carbonate platform 17 . The base of the composite section ( Fig. 1) corresponds to the post-Marinoan glaciation diamictitedolostone contact that is globally dated at 635. 5 Ma (refs 20,21). This is further constrained by Sr and C isotopes correlations and a Pb-Pb age of 627 ± 32 Myr (see Supplementary Information of ref. 22). The age of the top of the section is constrained by the acritarch assemblage of the overlying Nobres Formation (Cavaspina acuminata, Chlorogloeaopsis sp., Obruchevella sp., Ericiasphaera sp., Appendisphaera barbata, Tanarium irregulare, Tanarium conoideum and Micrhystridium pisinnum 23 ), which corresponds to the ECAP biozone of Grey (2005), ref. 24, that is, to an age interval of 570-580 Myr. This estimate is coherent with recently obtained detrital zircon and mica ages at 544 ± 7 Myr on the overlying Diamantino Formation 25 . We can therefore reasonably estimate a maximum depositional time (from 635 to 575 Ma) of 60 Myr for Mirassol do Oeste and Guia formations (300 m of sediments). Outcrop-based facies analysis, complemented by petrographic description of representative samples, reveals a transgressive systems tract, with the deepest part of the platform corresponding to a ramp depositional environment 17 .
Sulfur isotopes. We extracted sulfate from 16 micritic carbonate samples. Most CAS analyses (12 out of 16) were performed on samples from the Carmelo quarry, which contains the thickest post-glacial sedimentary sequence ( Fig. 1). At the base of the section (that is, in the direct aftermath of glaciation), d 34 S CAS is þ 14% (lower than that of the modern ocean, B þ 21%) and increases steadily up to B þ 38%, locally reaching þ 50.6% (n ¼ 2, Fig. 1 and Supplementary Table 1). Our range of values is consistent with data reported for other Marinoan post-glaciation deposits in Namibia 6 and North China 26 .
The isotope composition of pyrite was also analysed (n ¼ 35, from Terconi, Carmelo and Camil sections). Pyrite is absent from the first 50 metres of the composite section. Above the base, its content varies between 0.01% and 3.06%. Scanning electron microscopy investigations show that pyrite is present as a mixture of framboidal aggregates (that is, early diagenetic) and euhedral crystals 17 (that is, later-diagenetic). d 34 S pyr increases upsection from À 9.9 to þ 26.2%, with strong variations in the upper part of the section that are associated with lithological variations (Fig. 1). d 34 S pyr and d 34 S CAS are broadly correlated (R 2 ¼ 0.79), yielding a mean D 34 S CAS-pyr of 33.5 ± 5.3% (1s).

Discussion
Deviations of D 33 S CAS and D 36 S CAS from zero together with correlations of d 34 S CAS À D 33 S CAS (Fig. 2), d 34 S CAS À D 36 S CAS and D 36 S CAS À D 33 S CAS are observed for the first time in Neoproterozoic sections and result from mass-conservation effects 27 . Our results can only be produced in a non-steady-state system where F in aF out , the non-zero D 33 S CAS values being a consequence of the subtle interplay of sulfate input and removal from the ocean (Methods section). Sulfate input lowers the oceanic sulfate D 33 S CAS through mixing processes ( Supplementary Fig. 1), whereas sulfate removal from the ocean by hydrothermal or biological activity, the latter including bacterial sulfate-reduction (BSR) coupled to pyrite burial, increases the D 33 S CAS of the residual oceanic sulfate (Supplementary Figs 2  The analytical results and observed correlations can be quantitatively modelled by combining the dynamic equations of mass balance (equation (3)) and isotopic mass balance (for example, equation (4) for 34 S and 32 S), which govern seawater sulfate concentrations and its isotopic compositions: and where 34 a sulfide-sulfate is the fractionation factor between sulfur buried as pyrite and sulfur buried as evaporite at the global scale, which is usually inferred from the average sedimentary D 34 S SO4-pyr of equation (2), (ref. 10). Similar equations can be written for 33 S and 36 S isotope ratios. Because sulfur isotope fractionation factors 33 a and 36 a can be related to 34 a by the 33 b and 36 b exponents, respectively (for example, 33 b ¼ ln( 33 a)/ln( 34 a), ref. 27), these equations can all be written as a function of 34 a. Here for 33 S: 33 b is typically close to 0.515 for abiotic processes and varies from 0.509 to 0.516 for microbial sulfate-reduction 28,29 . As illustrated by equations (4 and 5), the modelled D 33 S-d 34 S trend of oceanic sulfate (blue line in Fig. 2) only depends on three parameters, namely F in /F out -ratio, 33 b and 34 a sulfide-sulfate values. Most importantly, as opposed to previous approaches, in this model, 34 a sulfide-sulfate is a free parameter that we can explore to determine the best-fit scenario, that is, it is not deduced from the D 34 S SO4-pyr values. The model does not depend on strong temporal constraints; therefore no a priori assumption was made on the initial sulfate residence time. Equally, no attempt was made to fit the isotope trend in Fig. 1, which would require a well-constrained deposition rate. However, we emphasize that the observed increase in ARTICLE F in /F out , 33 b and 34 a sulfide-sulfate solving simultaneously equations (3)(4)(5). All parameters, namely 33 b À and 36 b À factors, F in /F out , 34 a sulfide-sulfate and S-isotope compositions of inputs were considered constant with time. Initial oceanic sulfate-sulfur isotope values were chosen from the data at the base of the sequence, which correspond to the onset of the observed excursion with d 34 S initial B þ 12%, D 33 S initial ¼ À 0.016% and D 36 S initial ¼ À 0.1% (Fig. 2). The initial values may in fact have been lower and, in this respect, our approach is conservative. Note that they are consistent with the few other available Neoproterozoic sulfate multi-isotopic values obtained on Amadeus Basin (Australia), MacKenzie Fold Belt (Canada) and Oman (refs 14,30). We explored the space of solutions for different combinations of F in /F out -ratio, 33 b À factor and 34 a sulfide-sulfate by calculating the F in /F out -ratio that best fits the d 34 S À D 33 S slope defined by our data for a given set of 34 a and 33 b values. For that, we used 51 values of 33 b (from 0.511 to 0.516) and 231 values of 34 a (expressed as 1,000ln( 34 a) from À 17% to À 60%). A total of 11,781 combinations of 33 b, 34 a and F in /F out values compatible with the observed d 34 S À D 33 S slope were produced (Fig. 3). We can restrict step by step the space of solution of our model using constraints available for d 34  Valid combinations of 33 b, 34 a and F in /F out are restricted to the coloured, lower half-right of each panel of Fig. 3, which is delimited by curve #1. There are no viable solutions above curve #1 because the corresponding 33 b and 34 a-values would produce D 33 S too close to zero compared with our data (Fig. 2). The first constraint used here is the d 34 S CASfinal À d 34 S CASinitial parameter. Combinations of 33 b, and 34 a compatible with our data for different d 34 S CASfinal À d 34 S CASinitial values are shown in Fig. 3a. Curve #2 delimitates the field below which the difference between d 34 S CASfinal and d 34 S CASinitial values ( þ 38% and þ 12%, respectively) is too low (that is, d 34 S CASfinal À d 34 S CASinitial o26%) to account for the high d 34 S CAS -values measured in our section. Therefore, only results above curve #2 (that is, d 34 S CASfinal À d 34 S CASinitial 426%) are valid. This further constrains the space of solution to between curves #1 and #2, with 33 b-exponent o0.514 and 1,000ln( 34 a)o À 30%, as defined by the intercept between curves #1 and #2 (Fig. 3a). Figure 3b represents the range of F in /F out -ratios compatible with our model. It illustrates the calculated F in /F out -ratios that fit the observed d 34 S CAS versus D 33 S CAS slope for each pair of 33 b and 34 a. The white curve #3 highlights steady-state conditions, where F in ¼ F out . This figure clearly shows that for the field defined by curves #1 and #2, F in /F out -ratios are always below unity (F in oF out ), pointing to a decrease in oceanic sulfate concentration through time. This demonstrates that whatever the input parameters, the observed trends between d 34 S CAS versus D 33 S CAS cannot be reproduced in a steady-state model (Supplementary Fig. 4). Taken together, Fig. 3a and Fig. 3b constrain the solution space to F in oF out and 33 bo0.514 (between curves #1 and #2). Another interesting outcome of the model is that b 33 -values are distinct from those expected for sulfate-reduction to hydrogen sulfide under thermodynamic equilibrium (green curve in Fig. 3a), a result consistent with previous studies [31][32][33] .
The space of solution can be further restricted taking into account the fact that 36 b must also fit the observed D 33 S versus D 36 S correlation (Fig. 4a). Figure 3c represents the combinations of 33 b and 34 a compatible with our data for different 36 b-values. Our space of solutions (between curves #1 and #2) is compared with 34 a sulfide-sulfate data and their respective 33 b and 36 b exponents estimated from batch culture experiments for the two main sulfur metabolisms, namely SD and BSR (refs 28,29). Our space of solutions intersects only the field delimited by the BSR data set while SD data (red dots in Fig. 3c) plot outside. This allows us to assume that BSR is the main mechanism leading to sulfate drawdown and to further restrict our space of solutions to its intersection with the BSR cultures data set. It is worth noting, however, that available data from culture experiments are still limited and present significant variability. In Fig. 5, all available 34 a sulfide-sulfate , 33  (full grey circle). Although available data most clearly relate to BSR, one cannot rule out that part of the signal may reflect a mixed signature between SD and BSR organisms (Fig. 5).
The final space of solutions in Fig. 3 is represented by a grey polygon which corresponds to the following combination of values: 34 a sulfide-sulfate B0.960±0.005, 33 bB0.5125±0.0005 and F in /F out ¼ 0.30 ± 0.25. The high sensitivity of the model allow to have precise values for the best-fit scenario, indeed slight modifications of each modelled parameter ( 34 a sulfide-sulfate , 33 b and F in /F out ) shows significant effects on the d 34 S CAS versus D 33 S CAS slope (Fig. 4b).
This multi-isotopic approach allows us, using the best-fit values deduced above for 33 b, 34 a and F in /F out in equations (4 and 5), to quantify the contraction of the sulfate reservoir responsible for the observed increase in both d 34 S CAS and D 33 S CAS (blue line in Fig. 2). For the strong increase in d 34 S-value from þ 12% at the base of the section to þ 38% at the top, the model indicates that water column sulfate concentrations decreased dramatically by nearly 50%. In other words, at the end of deposition of the Guia Formation, sulfate concentration would be only half its initial value (Fig. 1). A single extreme d 34 S CAS value of þ 50% is observed in a typical event bed at the top of the section 17 ,   Fig. 3a,b, and the field defined by culture experiments- Fig. 3c). characterized by hummocky cross stratification. Given its sedimentary characteristics, it is unclear whether this single extreme value should be considered; if representative, it would indicate a 60% decrease compared with the sulfate concentration observed at the base of the Araras composite section. The high d 34 S values reported in post-Marinoan glacial deposits from Namibia 7 and North China 26 can also be accounted for by a drawdown of oceanic sulfate in a dynamic non steady-state sulfur cycle, without invoking the extreme modifications in both f pyr and d 34 S in needed in a steady-state model approach.
The model also provides fundamental constraints on the lower limit of marine sulfate concentration during the Neoproterozoic. For sulfate concentrations below 1 mM, the sulfur isotope fractionation associated with BSR decreases, reaching 0% below 200 mM (ref. 38) or lower 39 . The results showing 34 a sulfide-sulfate close to 0.960 (D 34 S ¼ À 40%) without significant changes throughout the sequence, suggests that sulfate concentrations remained well above 200 mM even by the end of the sulfate reservoir drawdown. The fractionation factor would otherwise have decreased significantly 38 . This constrains the lower limit of marine sulfate concentration in the immediate aftermath of the glaciation to well above 400 mM. These estimates are consistent with the work of Kah et al. 13 , who constrained the upper limits of Neoproterozoic marine sulfate concentrations to between 7 and 10 mM using a completely independent approach.
The post-glacial marine environment recorded in the Snowball aftermath deposits is thought to be characterized by an enhanced delivery of phosphate 40 and bioavailable iron 41 . The resulting planktonic bloom suggested by Elie et al. 42 accompanying the Snowball deglaciation may have increased the organic matter flux to the sediment exhausting its dissolved O 2 content and enhancing anaerobic respiration of organic matter. We thus suggest here that post-glacial conditions were adequate for anaerobic sulfatereduction metabolism triggering significant sulfate removal from the water column. If widespread, such a sulfate drawdown by BSR and pyrite burial would have important consequences for other biogeochemical cycles including the global oxygen budget during the end of the Neoproterozoic Era 2,43 . Pyrite burial plays a key role in the long-term sources of atmospheric O 2 and its transient increase may have contributed to the accumulation of O 2 in the oceans and atmosphere 44,45 . While burial of organic matter is generally considered to be the most important driver of O 2 production on geological timescales (present-day value 10 Â 10 12 mol O 2 per year), refs 44,45, today pyrite burial contribution to this long-term redox balance is of the same order of magnitude 44,45 .
We therefore estimated a minimum range of O 2 fluxes produced by the post-glacial intense BSR activity and pyrite burial implied by our model. Because the flux of sulfate delivery to the ocean (F in ) is unknown for this period, we only calculate the net O 2 flux due to a 50% decrease of the ocean sulfate concentration without continuous sulfate inputs to the ocean. In that sense the fluxes in the following calculation has to be taken as minimum values. From a global perspective and assuming the modern ocean volume 46 , the combined results of our study and others 13,46-48 constrain the total marine sulfate reservoir of the late Neoproterozoic to between 5 Â 10 17 mol and 13 Â 10 18 mol (SO 4 2 À concentrations of 0.4-10 mM, respectively). Considering the stoichiometry of BSR coupled to pyrite precipitation, where reduction of 1 mol SO 4 2 À and burial as pyrite equates the release of 2 mol O 2 , ref. 49, a 50% decrease in the late Neoproterozoic sulfate reservoir by BSR coupled to pyrite burial is equivalent to the net total production of 0.5-13 Â 10 18 mol of O 2 .
To express this number as a flux of O 2 , one needs tight temporal constraints that are lacking for the Araras Group (see above). Figure 6 presents the flux of O 2 produced from pyrite buried per year as a function of the duration of the d 34 S CAS increase episode.
Using the conservative deposition time of 60 Myr, the decrease in sulfate concentration observed would lead to an O 2 production comprised between 0.1 and 1.2 Â 10 12 mol per year (Fig. 6). These values correspond respectively to 0.1 and 3.2% of the modern production of O 2 attributed to pyrite burial per year (Fig. 6) 49,50 . If the increase in d 34 S-values has occurred over a shorter interval of time, which is supported by the fact that the observed increase in d 34 S-values occurs within the first 175 m of the section, the O 2 flux is then higher. For a depositional time period of 10 Myr, the O 2 production flux would rise to between 1.5 and 20% of the present O 2 production by BSR activity (Fig. 6). We conclude, therefore, that in the aftermath of the Marinoan glaciation enhanced BSR and pyrite burial represents a viable mechanism contributing to the Neoproterozoic oxygenation event of the ocean-atmosphere system.

Methods
Sample preparation. Five to 100 g of carbonate samples (with carbonate contents typical 470 wt% of the total rock) were powdered, leached of soluble sulfates in a 5% NaCl solution, followed by four rinses in deionized (DI) water. This step was repeated three times, then the powder was dissolved in 4 N HCl (12 h). The acidified samples were filtered, on a 0.45-mm nitrocellulose paper and an excess of 250 g l À 1 of BaCl 2 was added to the filtrate to precipitate BaSO 4 .
Multiple sulfur analyses. The samples were prepared and analysed for their multiple sulfur isotope compositions at the Stable Isotope Laboratory of the Institut de Physique du Globe de Paris.The barium sulfate was subsequently reacted with Thode reagent 51 in a helium atmosphere extraction line. The released H 2 S was converted to silver sulfide (Ag 2 S) by reaction with a silver nitrate solution and silver sulfide was fluorinated to produce SF 6 . The d 34 S values are presented in the standard delta notation against V-CDT with an analytical reproducibility ofr0.1%. We report these values against an assumed D 33 S and D 36 S for Vienna Cañon Diablo troilite (V-CDT) that yields d 34 S, D 33 S and D 36 S values for IAEA S-2 (n ¼ 23) of 5.224%, 0.030 and À 0.203%, respectively. Based on duplicate and triplicate analyses, uncertainties of D 33 S and D 36 S values by the SF 6 technique are estimated at 0.01 and 0.2% in 2s, respectively.
Model and concepts associated with non-zero D 33 S and D 36 S. The plot of 33 S/ 32 S and 34 S/ 32 S fractionations displays a slight curvature expressed by: The b-exponent is not arbitrary and can be deduced from the high-temperature approximation of the reduced partition function, from the mass (in atomic mass unit) of the considered isotopes, for example, ref. 52. Thus, for 33 S/ 32 S and 34 S/ 32 S  A mixing between two isotopically different pools (A and B) will fall along a mixing line (equations (9); ref. 55) that deviates from the theoretical equilibrium curve: where X denotes for the fraction of A. In other words, two sulfur reservoirs at isotope equilibrium will lie on the same fractionation line, with the same D 33 S. In contrast mixing will be expressed along a secondary fractionation line with negative D 33 S. The resulting D 33 S-anomaly is maximum for 50% mixing being approximately À 0.05% for two pools differing by 40% in d 34 S. Given that the mixing of two pools leads to negative D 33 S, the formation of two sulfur pools at isotope equilibrium will move them along a secondary fractionation line above that of the starting composition. These effects can be enhanced by Rayleigh distillation (that is, open system fractionation) and can be demonstrated starting from the well-known equation 56 : With d A the isotopic composition of the residual component of A and and d A(0) being the starting isotope composition, f is the residual sulfate concentration in our case, and a is the isotopic fractionation corresponding to the given system studied (bacterial sulfato-reduction in our study).
Data availability. All results that support the findings of this study are available in Supplementary Table 1.