L Guibourdenche

P Cartigny

F Delapierre

M Natalicchio

G Aloisi

Cryptic sulfur cycling during the formation of giant gypsum deposits

Keywords: cryptic sulfur cycling, multiple sulfur isotopes, Messinian Salinity Crisis, Primary Lower Gypsum, microbial sulfate reduction, sulfide oxidizing bacteria Highlights

Orbitally-driven hydrological changes controlled sulfur cycling at the onset of Mediterranean Messinian Salinity Crisis.

-In marginal basins, hypoxic to euxinic conditions prevailed during humid period and marls deposition. -During gypsum formation, microbial sulfate reduction (MSR) and near complete reoxidation of sulfide dominated sulfur cycle. -More than 80% of the sulfate was cycled through this cryptic sulfur cycling.

-Benthic organic carbon oxidation rates by MSR were comparable to those of the most productive areas of the modern ocean.

Introduction

At the end of the Miocene, starting from about 7.2 Ma, the progressive reduction of water exchanges between the Atlantic Ocean and the Mediterranean Sea dramatically affected Mediterranean environments, increasing their sensitivity to orbitally-driven climate perturbations [START_REF] Bulian | Messinian West Alboran Sea record in the proximity of Gibraltar: Early signs of Atlantic-Mediterranean gateway restriction[END_REF]Hsü et al., 1973;[START_REF] Kouwenhoven | Paleoenvironmental evolution of the eastern Mediterranean during the Messinian: Constraints from integrated microfossil data of the Pissouri Basin (Cyprus)[END_REF]. Hydrological restriction of the Mediterranean led to an increase of salinity, a decrease in water column oxygenation and changes in biotic assemblages [START_REF] Bulian | Messinian West Alboran Sea record in the proximity of Gibraltar: Early signs of Atlantic-Mediterranean gateway restriction[END_REF]Cita et al., 1978;[START_REF] Kouwenhoven | Paleoenvironmental evolution of the eastern Mediterranean during the Messinian: Constraints from integrated microfossil data of the Pissouri Basin (Cyprus)[END_REF]. Thereafter, the sensitivity of the Mediterranean to precessional cycles became significantly more pronounced, and in turn greatly impacted the biogeochemistry of the basins [START_REF] Natalicchio | Paleoenvironmental change in a precession-paced succession across the onset of the Messinian salinity crisis: Insight from element geochemistry and molecular fossils[END_REF][START_REF] Sabino | The response of water column and sedimentary environments to the advent of the Messinian salinity crisis: insights from an onshore deep-water section (Govone, NW Italy)[END_REF]. This has been recorded as cyclical deposition of sapropel/diatomite-rich marls (Krijgsman et al., 1999) and marl/carbonate alternations [START_REF] Mancini | Calcareous Nannofossil Size and Abundance Response to the Messinian Salinity Crisis Onset and Paleoenvironmental Dynamics[END_REF] corresponding to periods of summer insolation maxima and minima, respectively. Lithological cyclicity is accompanied by variations in the size and abundance of fossils of marine fauna [START_REF] Mancini | Calcareous Nannofossil Size and Abundance Response to the Messinian Salinity Crisis Onset and Paleoenvironmental Dynamics[END_REF], as well as oscillations in the composition and isotopic signature of lipid biomarkers [START_REF] Natalicchio | Paleoenvironmental change in a precession-paced succession across the onset of the Messinian salinity crisis: Insight from element geochemistry and molecular fossils[END_REF][START_REF] Sabino | Climatic and hydrologic variability in the northern Mediterranean across the onset of the Messinian salinity crisis[END_REF]. Together, this evidence suggests the existence of a precessional-scale alternation between a highly productive/stratified water column and a vertically mixed water column.

The progressive intensification of hydrological restrictions between the Atlantic and the Mediterranean culminated at 5.97 Ma with the onset of evaporite precipitation that marked the beginning of the Messinian Salinity Crisis (MSC). During the first phase of the MSC, up to 220 m of cyclical alternations of gypsum and organic-rich deposits (characteristic of arid and humid periods, respectively) forming the Primary Lower Gypsum (PLG) unit (Fig. 1) accumulated in the shallow part of marginal to intermediate depth basins of the Mediterranean (Lugli et al., 2010;Raad et al., 2021;[START_REF] Roveri | The Messinian Salinity Crisis: past and future of a great challenge for marine sciences[END_REF]. During this time, marginal and intermediate depth basins still received significant marine water inputs, as indicated by the 87Sr/86Sr and S of gypsum and by the occurrence of marine biota (Lugli et al., 2007, Lugli et al., 2010 andreferences therein;[START_REF] Pellegrino | Tiny, glassy, and rapidly trapped: The nano-sized planktic diatoms in Messinian (late Miocene) gypsum[END_REF][START_REF] Reghizzi | Orbitally Forced Hydrological Balance During the Messinian Salinity Crisis: Insights From Strontium Isotopes (87Sr/86Sr) in the Vena del Gesso Basin (Northern Apennines, Italy)[END_REF].

During humid periods of precession cycles, non-evaporitic interbeds deposited in stratified basins where N2-fixing bacteria significantly contributed to primary productivity [START_REF] Isaji | Diazotrophy Drives Primary Production in the Organic-Rich Shales Deposited Under a Stratified Environment During the Messinian Salinity Crisis (Vena del Gesso, Italy)[END_REF] and anoxygenic phototrophic bacteria populated the upper water column (Sinninghe Damsté et al., 1995). Meanwhile, microbial mats dominated by colorless sulfide oxidizing bacteria (CSB) thrived at the seafloor (Dela [START_REF] Pierre | Flocculent layers and bacterial mats in the mudstone interbeds of the Primary Lower Gypsum unit (Tertiary Piedmont basin, NW Italy): Archives of palaeoenvironmental changes during the Messinian salinity crisis[END_REF]. Anoxygenic phototrophic bacteria and CSB mats develop in modern euphotic zone and benthic boundary respectively. They reoxidize H2S fluxes produced by microbial sulfate reduction (MSR) back to elemental sulfur or sulfate [START_REF] Brüchert | Biogeochemical and physical control on shelf anoxia and water column hydrogen sulphide in the Benguela coastal upwelling system off Namibia[END_REF][START_REF] Callbeck | Sulfur cycling in oceanic oxygen minimum zones[END_REF][START_REF] Jørgensen | SULFUR BIOGEOCHEMICAL CYCLE OF MARINE SEDIMENTS[END_REF]. The presence of sulfide oxidizing metabolisms in both the water column and at the sediment-water interface strongly suggest the prevalence of euxinic conditions during marls/organic-rich shale deposition. A. Sulfur δ34SSO4 and δ34SIS isotopic curve of the PLG succession of the VdG basin, calibrated with the 100ka orbital eccentricity and precession curves [START_REF] Laskar | La2010: a new orbital solution for the long-term motion of the Earth[END_REF]. Note the nearly invariant δ34SSO4 isotopic variations compared to δ34SIS. B. Panoramic view of the 220m thick VdG succession constituted by an alternation of thick gypsum beds and marl intervals (PLG unit). C. Distribution of Messinian evaporites in the Mediterranean and the location of the VdG basin (pink star) (see supplementary section 2.1 for references).

Precession maxima (arid periods) were characterized by rapid vertical accretion of bottom grown gypsum crystals [START_REF] Reghizzi | Orbitally Forced Hydrological Balance During the Messinian Salinity Crisis: Insights From Strontium Isotopes (87Sr/86Sr) in the Vena del Gesso Basin (Northern Apennines, Italy)[END_REF] -which allowed for the relatively good preservation of organic remains (Dela Pierre et al., 2015;Schopf et al., 2012) (Fig. 2). These bottom-grown gypsum crystals make up the bulk of the PLG unit, and contain abundant planktic marine diatoms, implying enhanced primary productivity in the upper water column during gypsum deposition (Dela Pierre et al., 2015;[START_REF] Pellegrino | Tiny, glassy, and rapidly trapped: The nano-sized planktic diatoms in Messinian (late Miocene) gypsum[END_REF]Schopf et al., 2012). The molecular fossil assemblages entrapped within Messinian gypsum crystals and their compound specific carbon isotopes are indicative of a diverse microbial community, including marine Thaumarchaeota communities living in the upper water column, and sulfate-reducing bacteria and CSB-forming microbial mats at the sea floor (Natalicchio et al., 2022). The latter benthic microbial mats are fossilized within the bottom-grown gypsum in the form of denselypacked networks of filamentous structures ("spaghetti-like" structures) together with pyrite and dolomite, implying that MSR activity occurred during gypsum formation (Dela Pierre et al., 2015;Natalicchio et al., 2022;Panieri et al., 2008;Schopf et al., 2012) (Fig. 2D-I).

Together, these petrographic and geochemical observations point towards a highly dynamic sulfur cycle during PLG formation, with sulfur cycled amongst microbial sulfate reducers and sulfide-oxidizing bacteria, as well as precipitated in its oxidized form as gypsum (CaSO4.2H2O).

Marine carbon and sulfur cycles are tightly linked, with MSR driving the remineralization of between 12 and 65% of the organic carbon flux to the seafloor [START_REF] Bowles | Global rates of marine sulfate reduction and implications for sub-sea-floor metabolic activities[END_REF][START_REF] Bradley | Widespread energy limitation to life in global subseafloor sediments[END_REF][START_REF] Canfield | Aquatic geomicrobiology[END_REF]. Constraining sulfur cycling during the formation of giant sulfate deposits is therefore critical to understand their significance for the carbon cycle. The accumulation of large gypsum deposits -comparable to the Mediterranean PLG-are known to greatly affect the long term carbon cycle and atmospheric oxygen levels by the removal of substantial amount of Ca2+ and SO42-ions from the ocean, (e.g. [START_REF] Halevy | Sulfate Burial Constraints on the Phanerozoic Sulfur Cycle[END_REF][START_REF] Shields | Evaporite weathering and deposition as a long-term climate forcing mechanism[END_REF]. However, the role of associated MSR and sulfide oxidation on carbon and sulfur cycles of basins where gypsum formed has been overlooked. This complex sulfur cycling during PLG formation is likely to have been widespread across Mediterranean marginal basins (Natalicchio et al., 2022;[START_REF] Pellegrino | Tiny, glassy, and rapidly trapped: The nano-sized planktic diatoms in Messinian (late Miocene) gypsum[END_REF]) and raises fundamental questions about the biogeochemical dynamic of these basins, including for instance: (i) What was the relative importance of microbially driven reduction of sulfate and reoxidation of sulfide, compared to the burial of sulfate through gypsum precipitation? (ii) How did sulfur cycling impact carbon cycling in Messinian Mediterranean marginal basins? (iii) How does this ancient system compare with modern marine environments?

In this study, we investigated the most complete PLG succession of the Mediterranean area in the Vena del Gesso basin (VdG) (Northern Apennines, Italy) (Fig. 1b). MSR and re-oxidative sulfur cycling result in specific multiple sulfur isotopes signatures, providing a tool to disentangle the relative influence of concurrent sulfur redox processes [START_REF] Johnston | Sulphur isotopes and the search for life: strategies for identifying sulphur metabolisms in the rock record and beyond[END_REF]Ono et al., 2006).

Therefore, to constrain the sulfur biogeochemical cycle during PLG deposition, we carried out measurements of multiple sulfur isotope ratios (34S/32S and 33S/32S) on gypsum sulfate and reduced sulfur minerals which co-precipitated with gypsum or that we extracted from marls, coupled to the measurement of oxygen isotopes of gypsum sulfate (18O/16O). This approach enabled us to estimate the relative fluxes of sulfide and sulfate in the basin and in turn to estimate the gross flux of organic carbon involved in MSR activity. In modern settings, this is typically achieved with incubation experiments and/or analyses of pore water geochemical profiles -tools that are not available when studying the rock record.

Geological setting

The PLG was deposited during the first stage of the MSC (5.97-5.60 Ma) in shallow to intermediate marginal basins of the Mediterranean, including southern Spain, Italy, Hellenic arc and Cyprus arc (Krijgsman et al., 1999;Lugli et al., 2010 and references therein) (Fig. 1c). This unit was deposited in approximately 380 kyrs [START_REF] Roveri | The Messinian Salinity Crisis: past and future of a great challenge for marine sciences[END_REF] and is formed by lithological cycles consisting of organic-rich (marls, shale and carbonate)/gypsum couplets (Fig. 1a). The deposition of these interbeds was controlled by precession-driven climate oscillations (Lugli et al., 2010;[START_REF] Roveri | The Messinian Salinity Crisis: past and future of a great challenge for marine sciences[END_REF]. Tuning of the lithological cycles to the astronomical reference curve (Fig. 1a) allows the precise dating of each cycle and bed by bed correlation of the PLG unit at the Mediterranean scale (see the pink dots in Fig. 1c) (Krijgsman et al., 1999;[START_REF] Roveri | The Messinian Salinity Crisis: past and future of a great challenge for marine sciences[END_REF]. Here we studied the VdG in the Northern Apennines (Italy) that exposes the most complete succession of PLG deposits in the whole Mediterranean area (Fig. 1b). Each PLG cycle of the VdG basin starts with a thin (dm to m thick) organic matter-rich (total organic carbon from 0.2 to 3.1%; [START_REF] Lugli | The Messinian 'Vena del Gesso' evaporites revisited: characterization of isotopic composition and organic matter[END_REF] marly bed (Fig. 2a).

The thin marly beds are followed in each PLG cycle by a 4 to 30 m thick bed of primary gypsum (Fig. 2a). The gypsum beds of the first three cycles (PLG1 to PLG3) are mostly composed of bottom-grown, vertically oriented, twinned selenite crystals (up to 2.5 m high) of the massive selenite lithofacies (Lugli et al., 2010).

Starting from the gypsum beds of the 4th PLG cycle, the massive selenite lithofacies is followed upwards by the banded selenite lithofacies composed of less than 10 cm-thick crusts of bottomgrown twinned gypsum crystals separated by thin carbonate laminae possibly reflecting cyclical fluctuation of the gypsum saturation interface (Lugli et al., 2010) (Fig. 2b).

From the 6th PLG cycle to the top of the section, massive and banded selenite are followed by the branching selenite lithofacies consisting of cm-sized clear bottom-grown selenite crystals growing inclined or horizontally and grouped together in nodules or lenses separated by thin fine grained carbonate or gypsum laminae (Lugli et al., 2010). The absence of karst features in any of these lithofacies suggests that PLG gypsum are fully subaqueous deposits formed in a basin permanently covered by waters supersaturated with respect to gypsum.

Isotopes to decipher sulfur cycling

δ34S and δ18O of sulfate

The sulfur and oxygen stable isotope compositions of pore water dissolved sulfate are frequently used to investigate the biogeochemical sulfur cycling in modern marine sediments [START_REF] Aller | Unsteady diagenetic processes and sulfur biogeochemistry in tropical deltaic muds: Implications for oceanic isotope cycles and the sedimentary record[END_REF][START_REF] Blonder | Impact of Aeolian Dry Deposition of Reactive Iron Minerals on Sulfur Cycling in Sediments of the Gulf of Aqaba[END_REF][START_REF] Gilhooly | A comprehensive sulfur and oxygen isotope study of sulfur cycling in a shallow, hyper-euxinic meromictic lake[END_REF]. These stable isotopic compositions are expressed in delta (δ) notations, representing the deviation of 34S/32S or 18O/16O ratios from international VCDT and SMOW standards respectively:

𝛿(‰) = ' ( )*+,-.*/ ( ,0+1/+./ -14 × 1000 (1)
With Rmeasured corresponding to the measured 34S/32S or 18O/16O ratios and Rstandard corresponding to the VCDT 34S/32S ratio or SMOW 18O/16O ratio.

Sulfate-reducing microbes preferentially metabolize 32S-and 16O-bearing dissolved sulfate.

Therefore, as dissolved sulfate concentration decreases and hydrogen sulfide concentration increases, residual pore-water sulfate is progressively enriched in 34S and 18O (Fig. 3) (Brunner et al., 2005).

However, the effect of MSR on sulfate and hydrogen sulfide concentration profiles and isotopic compositions can be obscured by re-oxidative processes; this happens because 34S-depleted sulfur metabolized by MSR is re-injected back into the 34S-enriched sulfate pool with little or no fractionation (Fig. 3). This results in a smaller net enrichment of dissolved sulfate 34S due to microbial sulfate reduction, and a lower net production of sulfide/consumption of sulfate.

Importantly, if all the hydrogen sulfide produced by MSR is entirely re-oxidized back to sulfate (Fig. 3) [START_REF] Aller | Unsteady diagenetic processes and sulfur biogeochemistry in tropical deltaic muds: Implications for oceanic isotope cycles and the sedimentary record[END_REF][START_REF] Blonder | Impact of Aeolian Dry Deposition of Reactive Iron Minerals on Sulfur Cycling in Sediments of the Gulf of Aqaba[END_REF][START_REF] Brüchert | Biogeochemical and physical control on shelf anoxia and water column hydrogen sulphide in the Benguela coastal upwelling system off Namibia[END_REF][START_REF] Callbeck | Sulfur cycling in oceanic oxygen minimum zones[END_REF], the sulfur cycle is "cryptic" -since there is no net isotopic effect on the 34S isotopic composition.

Furthermore, the shift in 18O sulfate isotopic values is often so small (<4‰) [START_REF] Aller | Unsteady diagenetic processes and sulfur biogeochemistry in tropical deltaic muds: Implications for oceanic isotope cycles and the sedimentary record[END_REF][START_REF] Blonder | Impact of Aeolian Dry Deposition of Reactive Iron Minerals on Sulfur Cycling in Sediments of the Gulf of Aqaba[END_REF] that it can be difficult to recognize the effect of cryptic sulfur cycling from natural variability of oceanic sulfate isotopic composition using this isotope system (Fig. 3) (Markovic et al., 2016). δ34SSO4 and δ18OSO4 isotopic composition of the dissolved sulfate ion in equilibrium with the PLG unit from VdG basin (dark pink) and the PLG unit of other marginal basins (light pink) (see Supplementary section 2.1 for references and 2.4 for calibration of δ34SSO4 and δ18OSO4 isotopic data from this study and from the literature). The range of isotope compositions of sulfate-bearing minerals and of sedimentary sulfides are represented as grey and light brown boxes, respectively (see Supplementary section 2.1 for references). The possible range of isotopic composition of riverine inputs is bracketed by the isotopic composition of sulfate-bearing minerals and sedimentary sulfides (see Supplementary section 2.1 for further details on riverine fluxes). The possible evolution of the isotopic composition of residual sulfate during MSR is indicated by the orange shaded area. The isotope composition of sulfate dissolved in top few centimeters of modern marine environments characterized by near complete reoxidation of H2S is represented by the blue box [START_REF] Aller | Unsteady diagenetic processes and sulfur biogeochemistry in tropical deltaic muds: Implications for oceanic isotope cycles and the sedimentary record[END_REF].

Multiple sulfur isotopes

Multiple sulfur isotopes (including 32S, 33S, 34S) have the potential to distinguish biogeochemically distinct hydrogen sulfide re-oxidation pathways occurring concomitantly with MSR -the signal of which might otherwise be obscured using traditional isotope systems [START_REF] Johnston | Sulphur isotopes and the search for life: strategies for identifying sulphur metabolisms in the rock record and beyond[END_REF]Ono et al., 2006;[START_REF] Zhang | Redox chemistry changes in the Panthalassic Ocean linked to the end-Permian mass extinction and delayed Early Triassic biotic recovery[END_REF].

The partitioning of 34S (34S/32S) relative to 33S (33S/32S) is predictable under thermodynamic equilibrium. However, intracellular biochemical branching processes will shift the minor 33S sulfur isotope partitioning (33S/32S) thus resulting in a deviation from the expected thermodynamic equilibrium relationship (Ono et al., 2006). This deviation is expressed as Δ33S

where:

∆ 88 𝑆 : = ;' == ( )*+,-.*/ == ( >?@A -14 × 1000B -C' =D ( )*+,-.*/ =D ( >?@A -14

E *F × 1000G (2)
The measurement of δ34S-Δ33S thus provides additional information as deviation from thermodynamic equilibrium can be linked to biological fractionation processes and mixing proportions of end-member sulfur species involved in sulfur cycling (Fig. 4) (Ono et al., 2006). Data points represent the Δ33S vs δ34S of sulfate and sulfide of gypsum (massive, banded, and branching selenite) and marl samples from the PLG unit of the VdG Basin. The dashed black arrow points to the isotopic composition of the first sulfide formed from oceanic sulfate by MSR considering the most depleted sulfide of marls measured in this study as a representative endmember. The curved arrow represents the calculated trajectory between this first sulfide and a pooled sulfide formed after the complete consumption of the initial sulfate reservoir. The blue line represents mixing of Messinian oceanic sulfate (SW) with the average isotopic composition of modern rivers (Freshwater), with friv =0 representing the marine endmember and friv=1 representing the riverine endmember. Note that the contribution of riverine SO4 to gypsum SO4 isotopic composition seems negligible.

The deviation from thermodynamic equilibrium is also expressed by the 33λ notation that defines the relationship between isotopic fractionation for 33S and 34S isotopes among two sulfur species:

88 𝜆 I = JK( == L M ) JK( =D L M ) (3)
with λ for thermodynamic equilibrium 33λeq = 0.515, and xαj = x/32Rproduct/x/32Rreactant, and x = 33 or 34.

Combining the 33λ and the classical 34Ɛ [= (34 xαj -1) × 1000] notations is particularly useful as λ and fractionation factors are dependent solely on the biogeochemical processes in question (e.g. MSR, biotic or abiotic sulfide oxidation). In addition, the combination of 33λ and 34Ɛ is not sensitive to the absolute isotopic value of the reacting sulfur species (Fig. 5) or the almost invariant temperatures relevant for sedimentary environments [START_REF] Eldridge | Theoretical estimates of equilibrium sulfur isotope effects in aqueous sulfur systems: Highlighting the role of isomers in the sulfite and sulfoxylate systems[END_REF]. Net fractionation factors 34ɛnet and 33λnet calculated for the gypsum (massive, banded, and branching selenite) and marl samples using the isotope compositions of the sulfate ion in equilibrium with gypsum and the sulfur isotope composition of iron sulfides (dots and crosses) compared with ranges of 34ɛ and 33λ values obtained in pure cultures of sulfur metabolizing microorganisms (shaded areas). For marls samples, where no gypsum precipitates occur, an isotope composition of the sulfate ion equal to the Messinian ocean value (Masterson et al., 2016) was set, based on the assumption that sulfides from marls were produced through the reduction of oceanic sulfate. See Supplementary section 4.1 for references on data from abiotic sulfide oxidation experiments (dark grey area with black dotted contour), data from pure cultures with sulfate reducers (light grey area), data from pure culture experiments with sulfur disproportionators and 34ɛnet and 33λnet fractionation factors measured in natural environments populated with biofilms dominated by chemolithotrophic sulfide oxidizers (white squares).

Material and methods

A total of 88 samples that are representative of each lithofacies of all cycles from the VdG succession were analyzed (Fig. 1A). Gypsum and marls from PLG cycle 1-2 were sampled in the Monticino quarry (44°13′29″N, 11°45′43″E) whereas cycles 3 to 16 were sampled in the Monte Tondo quarry (44°15′04″N, 11°40′13″E). Care was taken to avoid sampling the weathered surface and visible fractures in order to limit contamination.

Petrographic analyses

11 oriented samples were taken from the various gypsum lithofacies for petrographic observations. Fifteen petrographic thin sections and 10 thin slices (∼1 mm thick) were produced from the most representative samples by cleaving the selenite crystals along the main (010) cleavage planes with a razor blade. These were studied using transmitted, reflected, and UV-visible light microscopy, using a Leika DM 2700 P microscope with excitation filters 340-380

Nm at the Department of Earth Sciences of the University of Torino.

Oxygen and multiple sulfur isotope analyses

Sulfate from gypsum and iron sulfide from gypsum and marls were first reduced and transformed into Ag2S by wet chemistry. The multiple sulfur isotope composition was then measured using a MAT-253 mass spectrometer. Oxygen isotopic composition of sulfate from gypsum were measured on samples precipitated as BaSO4 (the full methods with references are detailed in supplementary material section 1).

Geochemical modeling

Description of the model

We interpret the multiple S isotope data using a steady-state, 3-box model based on a chemical and isotopic mass balance of sulfur in the SO42-, S° and H2S forms. In the model (see Fig. 6) sulfate is provided by oceanic (Fin) and riverine inputs (Friv) whereas it is removed from the basin through outflow into the open ocean (Fout) or gypsum precipitation (Fgyp). In addition, the sulfate can be reduced to hydrogen sulfide by microbial sulfate reduction (Fmsr), and hydrogen sulfide can be either buried as solid sulfide (Fis) or reoxidized (Fox). Sulfide reoxidation (biotic or abiotic) produces sulfate or an intermediate sulfur compound (elemental sulfur or sulfite) (Fox).

Finally, intermediate sulfur compounds can be reoxidized back to sulfate (Fox) or further disproportionated into hydrogen sulfide (Fdh) and sulfate (Fdo) (Fig. 6).

Fractionation factors

Sulfur redox or phase changes are often associated with a differential partitioning of sulfur isotopes. Therefore, the expressed fractionation Ɛnet and λnet between two co-occurring sulfur species of interest (typically sulfate and sulfide) results from the net effect of multiple processes involved in the cycling of these sulfur species (Fig. 5). To model the biogeochemical processes at play during the genesis of our samples, we assigned fixed fractionation factors established in pure culture experiments or environmental studies to each of the fractionating process considered in our model (see supplementary section 4.1), with the exception of sulfate reduction that potentially expresses a large range of isotopic fractionation [START_REF] Gomes | Sulfur isotope fractionation in modern euxinic systems: Implications for paleoenvironmental reconstructions of paired sulfate-sulfide isotope records[END_REF]Sim et al., 2011b;[START_REF] Wing | Intracellular metabolite levels shape sulfur isotope fractionation during microbial sulfate respiration[END_REF]. In the following paragraph, we detail the approach we chose to simulate fractionation of multiple sulfur isotopes during MSR.

Laboratory experiments and studies of natural environments suggest that there is a first-order control of the specific cellular rate of MSR (csSRR) on the extent of sulfur isotopic fractionation during MSR, with fractionation approaching equilibrium values (i.e. Ɛeq =-70‰ and λeq=0.515) when csSRR reaches minimum values (Sim et al., 2011a;Sim et al., 2011b;[START_REF] Stam | 34S/32S fractionation by sulfate-reducing microbial communities in estuarine sediments[END_REF]) (Fig. 5). The laboratory-observed csSRR at which these relationships have been established are often 1-3 orders of magnitude higher than the csSRR usually reported from natural marine environments [START_REF] Colangelo-Lillis | Diversity decoupled from sulfur isotope fractionation in a sulfate-reducing microbial community[END_REF]Jørgensen, 2021 and references therein). Some authors therefore argue that csSRR in marine environments are so low that cells always express fractionation factors close to equilibrium values (ƐMSR»-70‰ and λMSR»0.515) (e.g. [START_REF] Masterson | Interpreting multiple sulfur isotope signals in modern anoxic sediments using a full diagenetic model (California-Mexico margin: Alfonso Basin)[END_REF]. These authors propose that when apparent fractionation factors between sulfate and sulfide (Ɛnet and λnet) deviate from equilibrium values, the deviation is mainly a consequence of sulfide pooling produced by sedimentary or early diagenetic processes (see section 5.2.2.1 for further explanations on this mechanism) [START_REF] Masterson | Interpreting multiple sulfur isotope signals in modern anoxic sediments using a full diagenetic model (California-Mexico margin: Alfonso Basin)[END_REF].

Ɛnet>-70‰ and λnet<0.515 are often observed in systems where depositional controls are absent or negligible, implying that natural populations of microbial sulfate reducers can produce a large range of isotopic fractionations [START_REF] Colangelo-Lillis | Diversity decoupled from sulfur isotope fractionation in a sulfate-reducing microbial community[END_REF][START_REF] Gilhooly | A comprehensive sulfur and oxygen isotope study of sulfur cycling in a shallow, hyper-euxinic meromictic lake[END_REF]Gomes and Hurtgen, 2015 and references therein;litterature review in Sim et al., 2011a). In addition, modeling studies also suggest a first order control of csSRR on fractionation during MSR [START_REF] Bradley | Patterns of sulfur isotope fractionation during microbial sulfate reduction[END_REF][START_REF] Wing | Intracellular metabolite levels shape sulfur isotope fractionation during microbial sulfate respiration[END_REF]. Howbeit, secondary factors independent from csSRR might also influence the extent of isotopic fractionation expressed during MSR (e.g. [START_REF] Bradley | Patterns of sulfur isotope fractionation during microbial sulfate reduction[END_REF][START_REF] Brüchert | Controls on stable sulfur isotope fractionation during bacterial sulfate reduction in Arctic sediments[END_REF][START_REF] Colangelo-Lillis | Diversity decoupled from sulfur isotope fractionation in a sulfate-reducing microbial community[END_REF][START_REF] Stam | 34S/32S fractionation by sulfate-reducing microbial communities in estuarine sediments[END_REF][START_REF] Wing | Intracellular metabolite levels shape sulfur isotope fractionation during microbial sulfate respiration[END_REF].

Although it is not within the scope of this work to reconcile laboratory experiments with field observations and modelling studies, we point out that the correct estimation of csSRR in natural environments (in mol.cell-1.day-1) relies on: (i) a precise estimate of the number of physiologically active sulfate reducing microbial cells, and (ii) the absence of significative reoxidative sulfur cycling. Any overestimation of active MSR cell abundance [START_REF] Colangelo-Lillis | Diversity decoupled from sulfur isotope fractionation in a sulfate-reducing microbial community[END_REF] and/or underestimation of the rate of sulfide oxidation [START_REF] Jørgensen | SULFUR BIOGEOCHEMICAL CYCLE OF MARINE SEDIMENTS[END_REF][START_REF] Treude | Rapid sulfur cycling in sediments from the Peruvian oxygen minimum zone featuring simultaneous sulfate reduction and sulfide oxidation[END_REF] might lead to a large underestimation of csSRR. It is possible that the current observed range of csSRR in natural environments is biased towards low rates [START_REF] Jørgensen | SULFUR BIOGEOCHEMICAL CYCLE OF MARINE SEDIMENTS[END_REF], whereas csSRR is directly measured in pure culture experiments that allow the entire range of possible MSR-induced fractionations to be explored (Sim et al., 2011a;Sim et al., 2011b) (Fig. 5).

We consider that pure culture experiments linking fractionation factors ƐMSR and λMSR to csSRR provide a valid first order approximation of the behavior of sulfur isotope fractionation to MSR.

We used a nonlinear regression model in combination with pure culture experiment ƐMSR , λMSR and csSRR data to model the full range of MSR fractionation potentially expressed in natural environments (see supplementary 4.2 for detailed explanation). This approach does not consider the other environmental controls on MSR induced isotopic fractionation and should not be used to directly deduce a csSRR from our data but provide a simplistic way to model the ƐMSR and λMSR co-variation produced by MSR.

Results and discussion

Iron sulfides are found in all investigated gypsum and marl samples from VdG (Fig. 1A, supplementary section 3.1, and previously by Panieri et al. (2008) and Sinninghe Damsté et al. (1995)). Iron sulfides form from the reaction of hydrogen sulfide produced by MSR with reactive iron, thus providing evidence for an active biogeochemical sulfur cycle during the first phase of the MSC.

Marl interbeds

In marls, iron sulfides display the smallest variability, the lowest δ34SIS (average -32 ±5 ‰ vs CDT (1σ)) and the highest Δ33SIS (average 0.12±0.04‰ (1σ)) of all analyzed samples (Figure 1A, Figure 4, supplementary section 3.1). These isotopic values are typical of environments where the amount of available dissolved sulfate for MSR is effectively infinite (e.g. [START_REF] Gomes | Sulfur isotope fractionation in modern euxinic systems: Implications for paleoenvironmental reconstructions of paired sulfate-sulfide isotope records[END_REF][START_REF] Johnston | Sulphur isotopes and the search for life: strategies for identifying sulphur metabolisms in the rock record and beyond[END_REF]. We calculated net isotope fractionation factors 34ɛnet and 33λnet for marl samples considering a sulfate isotopic composition equal to that of the Messinian Ocean, based on the assumption that sulfide in marls was produced via the reduction of oceanic sulfate.

Marls fractionation factors fall within the MSR field (Figure 5.) and in fact represent the set of samples in this study that are the closest to the maximal fractionation expressed during MSR, with similar values to what is observed in modern organic-rich hypoxic to euxinic basins (e.g. [START_REF] Gomes | Sulfur isotope fractionation in modern euxinic systems: Implications for paleoenvironmental reconstructions of paired sulfate-sulfide isotope records[END_REF][START_REF] Johnston | Sulphur isotopes and the search for life: strategies for identifying sulphur metabolisms in the rock record and beyond[END_REF]. This is consistent with previous observations suggesting that these organic rich marls formed under stratified and frequently euxinic water columns (Dela [START_REF] Pierre | Flocculent layers and bacterial mats in the mudstone interbeds of the Primary Lower Gypsum unit (Tertiary Piedmont basin, NW Italy): Archives of palaeoenvironmental changes during the Messinian salinity crisis[END_REF][START_REF] Isaji | Diazotrophy Drives Primary Production in the Organic-Rich Shales Deposited Under a Stratified Environment During the Messinian Salinity Crisis (Vena del Gesso, Italy)[END_REF]Lugli et al., 2010;[START_REF] Sinninghe Damsté | Molecular indicators for palaeoenvironmental change in a Messinian evaporitic sequence (Vena del Gesso, Italy). I: Variations in extractable organic matter of ten cyclically deposited marl beds[END_REF] (Fig. 7a). (Natalicchio et al., 2022;Panieri et al., 2008;Schopf et al., 2012) and PLG deposits of other Mediterranean marginal basins (Dela Pierre et al., 2015;Natalicchio et al., 2022). The origin of the filaments has been disputed. They have been interpreted as cyanobacteria (Panieri et al., 2008), even if they are generally larger (~70µm diameter) than most modern cyanobacterial filaments (Schopf et al., 2012 and references therein) (Figure 2. E-I).

Gypsum beds

Moreover, the fact that filaments are mostly located both in the reentrant angle and along the vertical growth bands ((Figure 2.D-F) suggests that they formed a benthic mat located at the gypsum-bottom water interface and that filamentous microorganisms were entrapped within gypsum in their living position, rather than after having settled from the photic zone to the basin floor (Dela Pierre et al., 2015;Schopf et al., 2012). Finally, the presence of intracellular iron sulfide globules ((Figure 2.G-I) -interpreted as resulting from diagenetic transformation of sulfur globules (Dela Pierre et al., 2015;Schopf et al., 2012) -is a diagnostic feature of benthic giant CSB.

Such observations confirm the affiliation of the observed filament networks to chemotrophic microbial mats where both MSR and sulfide oxidation were active (Dela Pierre et al., 2015;Natalicchio et al., 2022;Schopf et al., 2012).

Isotopic composition of PLG gypsum deposits

5.2.2.1.

Influence of syn-depositionnal or early diagenetic processes on multiple sulfur isotopes

The sulfur isotopic composition of the dissolved sulfate ion in equilibrium with the VdG gypsum ranged between 20.73 and 24.11 ‰ for δ34SSO4, and between -0.001 and 0.049 ‰ for Δ33SSO4 (Fig. 1A,3,4 Supplementary Table S2). These values fall within a narrow range that is centered around the isotopic composition of the Messinian Ocean dissolved sulfate ion (δ34SMESS=22.2±0.2 ‰ and Δ33SMESS=0.035 ±0.013 ‰) (Masterson et al., 2016) 3,4).

(Figure 1.A,
In contrast to the almost invariant isotopic composition of the sulfate, the iron sulfides entrapped in the same gypsum samples display a large range of isotopic compositions (-36 < δ 34 SIS < +9‰ and -0.017 < Δ 33 SIS < 0.149), from high Δ 33 SIS -low δ 34 SIS, which is comparable to that of iron sulfides contained within marl deposits intercalated to the gypsum beds, to close to zero Δ 33 SIS and δ 34 SIS (Figure 1.A , 4 and Supplementary section 3.1 and 3.2).

The theoretical maximal sulfur isotopic fractionation induced by MSR (ƐMSR≈-70‰) generate sulfide isotope values of ~-50‰ in δ 34 SIS when the sulfate used during MSR has a marine isotopic composition (~+20‰). In some cases, the isotopic composition of iron sulfide differs substantially from this maximal estimate due to sulfide pooling controlled by local diagenetic or sedimentological processes (e.g. [START_REF] Masterson | Interpreting multiple sulfur isotope signals in modern anoxic sediments using a full diagenetic model (California-Mexico margin: Alfonso Basin)[END_REF]Pasquier et al., 2021b;[START_REF] Pasquier | Pyrite sulfur isotopes reveal glacial-interglacial environmental changes[END_REF]. Sulfide pooling occurs when the rate at which sulfate is consumed (and sulfide is produced) in a given environment overcomes the rate at which the sulfate pool is replenished by diffusion from the overlying sulfate-rich water column [START_REF] Fike | Rethinking the ancient sulfur cycle[END_REF]. This leads to the progressive enrichment in 34 S of the residual sulfate and produces sulfide pools following a Rayleigh-type distillation. The observed isotopic difference (Ɛnet) between oceanic sulfate values and sedimentary sulfide is therefore smaller than the actual MSR induced fractionation (ƐMSR) because the pool of sulfate used by MSR is not oceanic sulfate (and therefore replenishable) but an aliquot of it (and therefore finite). For example, sulfide pooling has been proven to control the isotopic composition of sulfides in cases where sedimentation rates are so high that porewater sulfate is disconnected from the overlying water column [START_REF] Pasquier | Pyrite sulfur isotopes reveal glacial-interglacial environmental changes[END_REF], or where high concentrations or reactivities of organic matter in sediments results in efficient sulfate depletion due to MSR (Pasquier et al., 2021b). In other cases, iron sulfides enriched in 34 S are observed due to sediment mixing by oceanic currents or bioturbation bringing together pooled 34 S-enriched sulfide with newly formed 34 S-depleted sulfide [START_REF] Fike | Rethinking the ancient sulfur cycle[END_REF]Pasquier et al., 2021a).

The petrographic evidence and isotopic signature of gypsum and co-occurring iron sulfides from the VdG deposits infers particular depositional settings that exclude the role of sulfide pooling in shaping the sulfur isotopic signatures measured in this study. We detail the reasoning in the following paragraphs.

In the δ 34 S vs Δ 33 S space, part of the gypsum iron sulfide samples fall along a curve (see the curved arrow in Figure 4.) consistent with the mixing of two sulfur reservoirs (Ono et al., 2006): (i) sulfide formed by MSR from an infinite, undistilled sulfate reservoir where the isotopic fractionation between sulfate and sulfide is close to equilibrium (corresponding here to the most depleted sulfide measured in this study with δ 34 S ~-40‰ and positive Δ 33 S ~+0.17‰) and (ii) what would seem to be pooled sulfide produced by the consumption of SO4 by MSR (and therefore characterized by an isotope composition that tends to approach that of initial, unreacted sulfate, as sulfate is consumed) (Ono et al., 2006;[START_REF] Zhang | Redox chemistry changes in the Panthalassic Ocean linked to the end-Permian mass extinction and delayed Early Triassic biotic recovery[END_REF].

To test the possible influence of sulfide pooling on the isotopic composition of our samples, we applied a Rayleigh type distillation to estimate the relative amount of sulfate necessary to be consumed to reproduce our data. We calculated that 55 to 95% of the sulfate reservoir should be consumed by MSR to reproduce the 34 S-enrichment in our samples in which δ 34 SIS falls between -20‰ and +10‰ (capturing 50% of the iron sulfides measured) (Figure 4.). We tested this with fractionation factors of -70‰ (i.e., the thermodynamic equilibrium fractionation between sulfate and hydrogen sulfide) and -58‰ (i.e., the net fractionation between the lowest measured δ 34 SIS and its corresponding δ 34 SSO4).

Such significant consumption of sulfate should either: (i) have led to the undersaturation of bottom waters with respect to gypsum, thus resulting in partial dissolution of the underlying gypsum, or (ii) have produced a substantial 34 S-enrichment of the isotopic signature of the residual sulfate pool preserved in gypsum. However, detailed petrographic studies of VdG gypsum here and in both [START_REF] Reghizzi | Orbitally Forced Hydrological Balance During the Messinian Salinity Crisis: Insights From Strontium Isotopes (87Sr/86Sr) in the Vena del Gesso Basin (Northern Apennines, Italy)[END_REF] and (Panieri et al., 2008) indicate that gypsum laminae have grown continuously without being interrupted by dissolution episodes (Figure 2.) and no substantial 34 S-enrichment of the gypsum phase (Figure 1., 3., and 4) have been observed in our samples, consistently with [START_REF] Lugli | The Messinian 'Vena del Gesso' evaporites revisited: characterization of isotopic composition and organic matter[END_REF].

Finally, a drastic drop in the sulfate concentration within CSB mats seems highly unlikely because (i) mats hosting CSB grew on gypsum lamina (Figure 2.) that would provide an effectively infinite source of sulfate upon gypsum undersaturation; however, no inter-laminae dissolution features are observed as stated above, and (ii) CSB can efficiently oxidize reduced sulfur to SO4 [START_REF] Brüchert | Biogeochemical and physical control on shelf anoxia and water column hydrogen sulphide in the Benguela coastal upwelling system off Namibia[END_REF][START_REF] Callbeck | Sulfur cycling in oceanic oxygen minimum zones[END_REF][START_REF] Treude | Rapid sulfur cycling in sediments from the Peruvian oxygen minimum zone featuring simultaneous sulfate reduction and sulfide oxidation[END_REF].

The combined action of biogeochemical production of SO4 by CSB, together with an infinite SO4 reservoir (represented by the gypsum substrate below the mat and the water column above the mat), exclude the possibility of a large reservoir effect due to sulfate consumption. We therefore conclude that pooling effects did not modify the isotopic composition of our samples.

Sulfide and sulfate mineral formation from the same sulfur pool

An alternative scenario is therefore needed to explain the large range of variation in δ 34 S isotopic composition of iron sulfides entrapped in gypsum (-58‰<Ɛnet<-14‰). Previous nano-SIMS studies have shown that the isotopic composition of sulfides measured within microbial mats varied with depth due to variable sulfate concentration and/or csSRR within the mat [START_REF] Fike | Micron-scale mapping of sulfur cycling across the oxycline of a cyanobacterial mat: a paired nanoSIMS and CARD-FISH approach[END_REF]. The isotopic compositions reported here, however, represent bulk measurements for each PLG gypsum cycle, and therefore are reflective of the most prominent processes governing the partitioning of sulfur species in the mat within several "generations" of CSB hosting mat (rather than local intra-mat conditions).

In the PLG deposits studied here, gypsum is a bottom-grown precipitate and is composed of an alternation of clear and turbid laminae ( Accordingly, iron sulfides co-precipitated with gypsum have multiple sulfur isotopic values (Δ 33 SIS -δ 34 SIS) that neither reflect the overall water column conditions of the basin nor the bottom water sulfur cycling, but instead provide information on the dynamic of sulfur cycling within the mats permineralized in gypsum.

We interpret the Δ 33 SIS -δ 34 SIS signature of iron sulfides as reflecting isotope fractionation processes associated with the biogeochemical sulfur cycle during gypsum precipitation (Figure 7). The same pool of sulfate was used to precipitate gypsum and to produce hydrogen sulfide

.2.3. Paired sulfate-sulfide isotopes as a fingerprint of reoxidative sulfur cycling

The combination of iron sulfide and gypsum co-precipitation and the absence of sulfide pooling influence here allows for the calculation of net fractionation factors ( 34 ɛnet -33 λnet) between sulfate from gypsum and their coeval iron sulfides in our samples.

These net sulfate-sulfide isotope fractionation factors fall between the range of isotopic fractionations observed in pure cultures of microbial sulfate reducers and those observed in pure cultures of sulfur disproportionators Figure 5). Therefore, we conclude that a process characterized by 33 λ higher than maximal 33 λMSR (>0.515) contributed, together with MSR, to define the isotopic composition of at least a part of the gypsum samples.

The candidate processes include biotic H2S oxidation ( 33 λH2S-SO4ox= 0.513-0.529) and intermediate oxidation state sulfur disproportionation ( 33 λH2S-SO4disp = 0.514-0.528) (see Table 2.1). In support of this idea, the ten most 34 S-enriched samples (δ 34 SIS between -20 and +8 ‰ and Δ 33 SIS between 0.01 and 0.05‰) deviate from the mixing curve ((Figure 4) that assumes sulfide originating from a reservoir affected by sulfide pooling. Indeed, the contribution from a pooled sulfide reservoir is highly unlikely (see section 5.2.2.1), implying that this deviation is due to an additional sulfide pool. This alternative sulfide pool has a higher Δ 33 SIS than seawater (0.035‰), and was possibly produced by disproportionation, as proposed in studies where a similar deviation was observed [START_REF] Zhang | Redox chemistry changes in the Panthalassic Ocean linked to the end-Permian mass extinction and delayed Early Triassic biotic recovery[END_REF].

This scenario is fully supported by the occurrence of two distinct sulfide pools in gypsum ((Figure 2., 4): the first one, associated with marine snow floccules ((Figure 2.J-L), was produced by MSR in an infinite sulfate reservoir; the second, associated with remains of sulfide oxidizing bacteria ((Figure 2G-I), results from the reoxidative part of the sulfur cycle.

In support of the significant occurrence of reoxidative sulfur cycling, the sulfur and oxygen isotopic composition of the dissolved sulfate ion in equilibrium with the VdG gypsum (20.73 < δ 34 SSO4 ‰ vs CDT < 24.11; -0.001 < Δ 33 SSO4 ‰ vs CDT < 0.049‰; 8.01 < δ 18 OSO4 ‰ vs SMOW < 10.84) (Figure 1., 3, 4 and supplementary Table S2 ) fall within a restricted range centered around the isotopic compositions of the dissolved sulfate ion of the Messinian Ocean (δ 34 SMESS=22.2±0.2

‰ and Δ 33 SMESS=0.035 ±0.013 ‰) (Masterson et al., 2016) and δ 18 OMESS=9±2‰ (Markovic et al., 2016)), implying negligible contributions from freshwater sulfate (Figure 1,3,4). Together with the very low iron sulfide concentration in gypsum (supplementary Figure S1), the lack of 34 S and 18 O enrichment of sulfate characteristic of MSR [START_REF] Kah | Geochemistry of a 1.2 Ga carbonate-evaporite succession, northern Baffin and Bylot Islands: implications for Mesoproterozoic marine evolution[END_REF]Melegy and Ismael, 2013) implies either that MSR activity during gypsum deposition was limited, or that the re-oxidation of sulfide into sulfate was nearly complete, i.e. the sulfur cycle was cryptic [START_REF] Aller | Unsteady diagenetic processes and sulfur biogeochemistry in tropical deltaic muds: Implications for oceanic isotope cycles and the sedimentary record[END_REF][START_REF] Blonder | Impact of Aeolian Dry Deposition of Reactive Iron Minerals on Sulfur Cycling in Sediments of the Gulf of Aqaba[END_REF]) (Figure 3.).

5.3. Modelling multiple sulfur isotopes in the Vena del Gesso basin.

Modeling strategy

To evaluate the importance of microbially-induced sulfide oxidation and to test the hypothesis of cryptic sulfur cycling, we quantified the relative contribution of biogeochemical processes in defining the multiple S isotope compositions of the VdG samples applying the numerical model described in section 4.3.

In order to find a best fit to our data we tested several scenarios of sulfur redox cycling (see supplementary section 9), each one characterized by different relative sulfur fluxes (controlled by dimensionless fx parameters defined in supplementary section 5) and by a combination of distinct isotope fractionations (see Table S3) according to the following scenarios:

(i)

MSR only (without reoxidation of sulfide).

(ii) MSR coupled to direct abiotic oxidation of sulfide to sulfate.

(iii) MSR coupled to phototrophic oxidation of sulfide to sulfate or elemental sulfur.

(iv) MSR coupled to partial oxidation of sulfide to elemental sulfur by chemolithotrophs and subsequent disproportionation of elemental sulfur.

(v) MSR coupled to partial oxidation of sulfide to sulfite and subsequent disproportionation of sulfite.

Supplementary section 9 provides full details on the modelling and results of the various biogeochemical scenarios. In all scenarios (i) to (v), the entire range (from 0 to 1) of fox was explored (fox refers to the fraction of H2S produced by MSR being reoxidized; see supplementary section 5 for a full description of the model equations and parameters). The isotope fractionation factors are fixed for all processes except for MSR, for which the full ƐMSR and λMSR ranges have been explored (see supplementary section 4.2 for further explanation).

Cryptic sulfur cycling during the formation of VdG gypsum layers

The best fit of the model results to the analytical dataset (Figure 8) is obtained with fmsr = 0.8 (representing the fraction of sulfur entering the system as SO4 and subsequently reduced by MSR, see Supplementary section 7), and scenario (iv), i.e. the partial oxidation of sulfide to elemental sulfur by chemolithotrophs and subsequent disproportionation of elemental sulfur (Figure 7.C and Figure 8). See Supplementary section 8 for further discussion on the reoxidative processes involved in sulfur cycling. In this best fit model run, 90% of the isotopic data in both diagrams (δ 34 SSO4 -Δ 33 SSO4, Figure 8.A; 34 ɛnet -33 λnet, Figure 8.B) plot within the area defined by 0.85 < fox < 1.

The combined sulfate and sulfide isotope dataset can therefore only be reproduced by considering that:

(1) at least 80% of the net incoming sulfate is reduced microbially to H2S rather than directly precipitated as gypsum (fmsr>0.8) (Fig. 8), and (2) nearly all of the MSR-produced H2S is reoxidized

Implications for carbon cycling and gypsum formation

To our knowledge, this is the first time that gypsum deposits are proven to be formed from sulfate sourced almost entirely from cryptic sulfur cycling. The average rate of vertical gypsum accumulation at VdG was 1.4 mm.year -1 [START_REF] Reghizzi | Orbitally Forced Hydrological Balance During the Messinian Salinity Crisis: Insights From Strontium Isotopes (87Sr/86Sr) in the Vena del Gesso Basin (Northern Apennines, Italy)[END_REF]. Given that our model results imply that more than 80% of sulfate constituting this gypsum was first reduced by MSR, we calculated that approximately 1.9 mmol.cm -2 .year -1 of sulfate was consumed by MSR. Then, by considering a rough C:S stoichiometry of 2:1 for the MSR reaction [START_REF] Jørgensen | SULFUR BIOGEOCHEMICAL CYCLE OF MARINE SEDIMENTS[END_REF], we estimated a gross benthic organic carbon flux of 3.8 mmol.cm -2 .year -1 remineralized by MSR during PLG gypsum formation. This is consistent with the large amounts of organic matter contained in the VdG interbeds. Marl deposits have a total organic carbon (TOC) content comparable to that of euxinic environments (Jovovic et al., 2020;[START_REF] Lugli | The Messinian 'Vena del Gesso' evaporites revisited: characterization of isotopic composition and organic matter[END_REF][START_REF] Sinninghe Damsté | Molecular indicators for palaeoenvironmental change in a Messinian evaporitic sequence (Vena del Gesso, Italy). I: Variations in extractable organic matter of ten cyclically deposited marl beds[END_REF], whilst TOC concentrations in gypsum range between 0.1 to 0.96% [START_REF] Lugli | The Messinian 'Vena del Gesso' evaporites revisited: characterization of isotopic composition and organic matter[END_REF], which is substantial considering the enormous dilution effect arising from gypsum depositional rates that are 25 to 250 times greater than marls [START_REF] Reghizzi | Orbitally Forced Hydrological Balance During the Messinian Salinity Crisis: Insights From Strontium Isotopes (87Sr/86Sr) in the Vena del Gesso Basin (Northern Apennines, Italy)[END_REF]). The benthic carbon flux involved in MSR we estimate for the VdG basin is comparable to that occurring in the world's most productive ocean regions, such as the Benguela and Peru-Chile margins, where the organic carbon remineralized through MSR ranges between 0.2 and 20mmol.cm -2 .yrs -1 [START_REF] Bowles | Global rates of marine sulfate reduction and implications for sub-sea-floor metabolic activities[END_REF][START_REF] Treude | Rapid sulfur cycling in sediments from the Peruvian oxygen minimum zone featuring simultaneous sulfate reduction and sulfide oxidation[END_REF].

In these areas, dense mats populated by microbial sulfate reducers and CSB, similar to the fossil microbial assemblages found in the VdG gypsum, fuel a cryptic sulfur cycle [START_REF] Brüchert | Biogeochemical and physical control on shelf anoxia and water column hydrogen sulphide in the Benguela coastal upwelling system off Namibia[END_REF][START_REF] Callbeck | Sulfur cycling in oceanic oxygen minimum zones[END_REF][START_REF] Treude | Rapid sulfur cycling in sediments from the Peruvian oxygen minimum zone featuring simultaneous sulfate reduction and sulfide oxidation[END_REF].

The above considerations do not necessarily imply that the VdG basin was characterized by similar primary productivity to modern oxygen minimum zones but, rather, by comparable benthic carbon fluxes and carbon and sulfur cycling rates. Indeed, aerobic consumption of organic carbon in the VdG basin may have been limited to the upper layer of the water column due to density stratification (Natalicchio et al., 2022), leading to high organic carbon export fluxes [START_REF] Pellegrino | Tiny, glassy, and rapidly trapped: The nano-sized planktic diatoms in Messinian (late Miocene) gypsum[END_REF] and large quantities of benthic organic carbon available for MSR (Figure 7.). Conversely, very high primary productivity controls the export flux of organic carbon in modern productive margins.

Our findings support a recent model of PLG deposition [START_REF] Aloisi | The geochemical riddle of "low-salinity gypsum" deposits[END_REF] which suggests that the activity of benthic CSB mats drove increases in dissolved sulfate concentration and triggering gypsum precipitation. The inhibition of gypsum precipitation in deeper, more distal lateral equivalents of PLG gypsum [START_REF] Natalicchio | Paleoenvironmental change in a precession-paced succession across the onset of the Messinian salinity crisis: Insight from element geochemistry and molecular fossils[END_REF][START_REF] Sabino | The response of water column and sedimentary environments to the advent of the Messinian salinity crisis: insights from an onshore deep-water section (Govone, NW Italy)[END_REF], together with our results

showing that more than 80% of the sulfate constituting the more marginal VdG gypsum has been cycled through a dynamic sulfur cycle, suggest that gypsum precipitation could have been controlled by the efficiency of the re-oxidative portion of the sulfur cycle, i.e. the cryptic cycling of sulfur.

Concluding remarks

Our work suggests a tight coupling of microbially driven sulfate reduction and sulfide oxidation in the Vena del Gesso basin.

Primary Lower Gypsum deposits accumulated during the first phase of the Messinian Salinity Crisis and characterized by identical lithofacies, stacking patterns, petrographic features and microbial fossils are widespread across the Mediterranean (see Figure 1.C) (Lugli et al., 2010;Natalicchio et al., 2022;[START_REF] Pellegrino | Tiny, glassy, and rapidly trapped: The nano-sized planktic diatoms in Messinian (late Miocene) gypsum[END_REF]. Whereas these deposits were previously thought to be restricted to semi-enclosed basins with a water column depth of no more than 200m (Lugli et al., 2010), recent studies have demonstrated that some were deposited at depths of up to 800m and in basins not necessarily restricted by a sill (Ochoa et al., 2015;Raad et al., 2021).

If the cryptic sulfur cycle we describe was also active during the deposition of gypsum elsewhere in the Mediterranean, this would imply that the restriction of the Atlantic connection at the end of the Miocene led to not only the widespread deposition of evaporites, but also the significantly altered the short-term carbon and sulfur cycles of the Mediterranean [START_REF] Waldeck | Calibrating the triple oxygen isotope composition of evaporite minerals as a proxy for marine sulfate[END_REF].

Further, the extensive occurrence of microbial mats populated by CSB on modern highproductivity continental margins or restricted basins, together with the results of the present study, suggest that cryptic sulfur cycling might have been prevalent in analogue evaporitic environments of the geological past.

standard using IAEA-S1 international standard that have a value of -0.3‰ VS-CDT. Long term reproducibility of IAEA-S1 (n=7) yielded δ 34 S=-0.30±0.2‰ and Δ 33 S=0.081±0.004. Error bars (1σ) on the whole procedure for sulfate and sulfide extraction and measurements were respectively calculated from replicate extraction and analyses of internal barite (n=7) and marls (n=3) standards. We obtained 0.2 and 0.004‰ of error on δ 34 S and Δ 33 S of sulfate and 0.4 and 0.008‰ for δ 34 S and Δ 33 S respectively for sulfides.

Oxygen isotopes of sulfate

Oxygen isotope analyses of BaSO4 (produced by re-precipitating the SO4 2-ion of gypsum in a BaCl2 solution) were conducted at University of Burgundy on a Elementar vario PYRO cube elemental analyzer in-line with an Isoprime 100 mass spectrometer in continuous flow mode;

analytical errors are ±0.4‰ (2σ) based on replicate analyses of the international barite standard NBS-127. The δ 18 O values were calibrated using the internal laboratory standard Edm-Tyr (δ 18 O = +8.6 ‰ V-SMOW) that was calibrated against the IAEA international reference NBS 127 (δ 18 O = +8.5 ± 0.2‰ V-SMOW).

Data from the literature used in main text figures

Figure 1A and 2A

In Figure 1.a ,3 and 4 of the main text we use δ 34 S, Δ 33 S and δ 18 O from other studies. This includes:

• Estimate of the δ 34 S and δ 18 O isotopic composition of Messinian oceanic dissolved SO4 2- (Markovic et al., 2016;Masterson et al., 2016;[START_REF] Turchyn | Oxygen Isotope Constraints on the Sulfur Cycle over the Past 10 Million Years[END_REF]) (Fig. 1a and 3). The δ 18 OSO4 of the Messinian ocean is taken from [START_REF] Turchyn | Oxygen Isotope Constraints on the Sulfur Cycle over the Past 10 Million Years[END_REF] and has a large error bar to account for the uncertainties pointed by (Markovic et al., 2016) linked to the measurement of δ 18 O in barite samples from which the curve of δ 18 OSO4 through time originates (Fig. 3).

• δ 34 S and δ 18 O of PLG deposits from the Spanish marginal basins of Sorbas (Evans et al., 2015;García-Veigas et al., 2018), Bajo Segura (García-Veigas et al., 2018), Murcia (Pierre, 1982) and Mallorca (García-Veigas et al., 2018), the Italian marginal Piedmont (El Kilany, 2018) and Caltanissetta basins (El Kilany, 2018;García-Veigas et al., 2018) and the Polemi and Kalavasos/Psematismenos basins of Cyprus (Pierre, 1982) (Fig. 3).

• δ 34 S and δ 18 O SO4 2-isotopic composition range produced by pyrite oxidation (Balci et al., 2007;Burke et al., 2018) (Fig. 3).

• Estimate of Mesozoic and Cenozoic evaporites SO4 2-δ 34 S and δ 18 O isotopic composition (Boschetti et al., 2011 and references therein;Claypool et al., 1980;Crockford et al., 2019;Ortí et al., 2010;Sakai, 1971;Utrilla et al., 1992) (Fig. 3).

• Modern sea water δ 34 S, Δ 33 S and δ 18 O SO4 2-isotopic composition [START_REF] Johnston | Placing an upper limit on cryptic marine sulphur cycling[END_REF](Fig. 1a,3 and 4).

Figure 1.c

In figure 1.c, modified after [START_REF] Lofi | Seismic Atlas of the Messinian Salinity Crisis markers in the Mediterranean Sea -Volume 2[END_REF], the pink dots represent the PLG units that were correlated at the Mediterranean scale [START_REF] Karakitsios | A record of the Messinian salinity crisis in the eastern Ionian tectonically active domain (Greece, eastern Mediterranean)[END_REF]Krijgsman et al., 1999;Lugli et al., 2010;[START_REF] Orszag-Sperber | The onset of the Messinian salinity crisis: Insights from Cyprus sections[END_REF]. In addition, numerous gypsum deposits have been (Hsü et al., 1978a;Hsü et al., 1978b, c;Hsü et al., 1978d;Hsü et al., 1973;[START_REF] Kastens | Site 654: Upper sardinian margin[END_REF][START_REF] Lugli | The deep record of the Messinian salinity crisis: Evidence of a non-desiccated Mediterranean Sea[END_REF] and gypsum from Albania [START_REF] Pashko | THE MESSINIAN EVAPORITES OF THE PREADRIATIC FORELAND BASIN (ALBANIA)[END_REF], Algeria [START_REF] Rouchy | The end of the Messinian salinity crisis: Evidences from the Chelif Basin (Algeria)[END_REF], Adriatic foredeep [START_REF] Manzi | The Messinian salinity crisis in the Adriatic foredeep: Evolution of the largest evaporitic marginal basin in the Mediterranean[END_REF], Balearic Promontory (Ochoa et al., 2015;Raad et al., 2021), Calabria [START_REF] Costanzo | Messinian twinned selenite from the Catanzaro Trough, Calabria, Southern Italy: field, petrographic and fluid inclusion perspectives[END_REF] in Italy, Cyprus [START_REF] Robertson | Depositional processes and basin analysis of Messinian evaporites in Cyprus[END_REF][START_REF] Rouchy | Stromatolites and Cryptalgal Laminites Associated with Messinian Gypsum of Cyprus[END_REF], Egypt [START_REF] Aref | Lithofacies characteristics, depositional environment and karstification of the Late Miocene (Messinian) gypsum deposits in the Northern Western Desert, Egypt[END_REF]Melegy and Ismael, 2013), Israel [START_REF] Cohen | Halite-clay interplay in the Israeli Messinian[END_REF][START_REF] Raab | The geological history of Messinian (upper miocene) evaporites in the Central Jordan Valley (Israel) and how strontium and sulfur isotopes relate to their origin[END_REF], Tunisia [START_REF] El | Messinian deposits and erosion in northern Tunisia: Inferences on Strait of Sicily during the Messinian Salinity Crisis[END_REF][START_REF] Harrab | Study of a Neogene basin dynamics: The "Bizerte basin", northeastern Tunisia: relevance to the global Messinian Salinity Crisis[END_REF][START_REF] Sghari | Un dépôt évaporitique du Miocène supérieur dans le Sud tunisien : une conséquence de la crise messinienne ? Géomorphologie : relief, processus[END_REF], Turkey [START_REF] Boulton | The Miocene-Pliocene boundary and the Messinian Salinity Crisis in the easternmost Mediterranean: Insights from the Hatay Graben (Southern Turkey)[END_REF][START_REF] Tekin | Sedimentology and paleoenvironmental evolution of Messinian evaporites in the Iskenderun-Hatay basin complex, Southern Turkey[END_REF] and Syria [START_REF] Hardenberg | Sedimentology of the NW margin of the Arabian plate and the SW-NE trending Nahr El-Kabir half-graben in northern Syria during the latest Cretaceous and Cenozoic[END_REF][START_REF] Rouchy | The end of the Messinian salinity crisis: Evidences from the Chelif Basin (Algeria)[END_REF].

Figure 5

In figure 5 of the main text we use data from laboratory experiments in order to compare them with the fractionation factors we calculated from our samples. Data from abiotic sulfide oxidation experiments comes from [START_REF] Eldridge | Rates and multiple sulfur isotope fractionations associated with the oxidation of sulfide by oxygen in aqueous solution[END_REF] (dark grey area with black dotted contour), data from pure cultures with sulfate reducers comes from [START_REF] Bertran | Deconstructing the Dissimilatory Sulfate Reduction Pathway: Isotope Fractionation of a Mutant Unable of Growth on Sulfate[END_REF][START_REF] Johnston | Multiple sulfur isotope fractionations in biological systems: A case study with sulfate reducers and sulfur disproportionators[END_REF][START_REF] Leavitt | Influence of sulfate reduction rates on the Phanerozoic sulfur isotope record[END_REF]Sim et al., 2011a;and Sim et al., 2011b) (light grey area) and data from pure culture experiments with disproportionators comes from [START_REF] Johnston | Multiple sulfur isotope fractionations in biological systems: A case study with sulfate reducers and sulfur disproportionators[END_REF] (white area with grey dashed contour). White squares represent the 34 ɛnet and 33 λnet measured in natural environments populated with biofilms dominated by chemolithotroph sulfide oxidizer [START_REF] Magnuson | Thiomicrorhabdus streamers and sulfur cycling in perennial hypersaline cold springs in the Canadian high Arctic[END_REF][START_REF] Zerkle | Sulfur isotope values in the sulfidic Frasassi cave system, central Italy: A case study of a chemolithotrophic S-based ecosystem[END_REF].

2.4

Calibration of δ 34 S, Δ 33 S and δ 18 O sulfate isotopic data from this study and from the literature δ 34 S and δ 18 O from this study have been set to the VSCDT and VSMOW scale respectively using recent calibration of NBS-127 reference material [START_REF] Brand | Assessment of international reference materials for isotope-ratio analysis (IUPAC Technical Report)[END_REF].

To compare these published data with ours, we corrected the reported data that were calibrated with obsolete values of NBS-127 by using the difference between the reported and newly calibrated δ 34 S and δ 18 O values of the standard following the method reported in [START_REF] Johnston | Placing an upper limit on cryptic marine sulphur cycling[END_REF]. In addition, because there is a small and almost temperature-independent isotopic fractionation between the dissolved SO4 2-ion and gypsum, we calculated the isotopic composition of the sulfate ion in equilibrium with gypsum by applying known fractionation factors (at earth surface temperature α 34 Gyp-SO4=1.00165 [START_REF] Raab | Sulfur isotopic variations during seawater evaporation with fractional crystallization[END_REF][START_REF] Thode | Sulfur-Isotope Geochemistry of Petroleum, Evaporites, and Ancient Seas[END_REF][START_REF] Van Driessche | Unraveling the Sulfate Sources of (Giant) Gypsum Crystals Using Gypsum Isotope Fractionation Factors[END_REF] ; 33 αGyp-SO4=1.00084 and α 18 Gyp-SO4=1.0035 [START_REF] Lloyd | Oxygen isotope behavior in the Sulfate-Water System[END_REF][START_REF] Van Driessche | Unraveling the Sulfate Sources of (Giant) Gypsum Crystals Using Gypsum Isotope Fractionation Factors[END_REF]). S1 and Fig. S1). The higher concentration of iron sulfide in marly interbeds are comparable to sulfide concentration found in typical ancient and modern euxinic environments (Jovovic et al., 2020). Instead, sulfide concentrations in gypsum samples is much lower and variable.

Isotopic data

Iron sulfide (IS) content and isotopic composition in gypsum and marls samples

Figure S1: Concentration of sulfides within gypsum and marls interbeds of the Vena del

Gesso succession as a function of their δ 34 SCRS isotopic composition 3.2 δ 34 S, δ 18 O and Δ 33 S isotopic compositions of sulfate from gypsum

Biogeochemical processes and their fractionation factors

Sulfur isotope composition of gypsum and pyrite accumulating in a marginal basin can vary as a function of (i) mixing between the ocean-and river-derived sulfate sources, (ii) the magnitude of gypsum and pyrite burial fluxes and (iii) the extent and mode of intra-basinal biogeochemical sulfur cycling and the associated isotope fractionation.

Within the marginal basin, sulfur is abiotically and microbiologically partitioned between the three sulfur pools of SO4 2-, S° and H2S via the following processes:

I. gypsum accumulation with rate Fgyp that extracts sulfate from the marginal basin:

SO4 2-+ Ca 2+ + 2H2O → Ca.SO4.2H2O (1)

II.

Dissimilatory microbial sulfate reduction (MSR) into sulfide, performed by sulfate reducing bacteria, with rate Fmsr:

SO4 2-+ 2CH2O → H2S + 2HCO3 - (2)
III. oxidation of sulfide produced by MSR to sulfate, with rate Fox, which can take place:

(iii.1) chemically, using molecular oxygen or metal oxides as electron acceptors [START_REF] Avrahami | The oxidation of the sulphide ion at very low concentrations in aqueous solutions[END_REF][START_REF] Zhang | The products from the oxidation of H2S in seawater[END_REF] (3.2)

H2S + 2O2 → SO4 2-+ 2H + (3) Or H2S + 3/2O2 → SO3 2-+ 2H + (3.1) and SO3 2-+ 1/2O2 → SO4 2-
Or 4FeOOH + 1/2HS -+ 8H + → 4Fe 2+ + 1/2SO4 2-+6H2O (Seitaj et al., 2015) (4)
Or 4MnO2 + HS -+ 7H + → 4Mn 2+ + SO4 2-+ 4H2O [START_REF] Aller | Complete oxidation of solid phase sulfides by manganese and bacteria in anoxic marine sediments[END_REF] (5) (iii.2) biologically, via phototrophic sulfide oxidation (green or purple sulfur bacteria): [START_REF] Findlay | A model of phototrophic sulfide oxidation in a stratified estuary[END_REF] (6) Or 2HS -+ CO2 + H + → 2S 0 + CH2O + H2O [START_REF] Findlay | A model of phototrophic sulfide oxidation in a stratified estuary[END_REF] (6.1) and

4CO2 + 2H2S + 4H2O → 4CH2O + 2SO4 2-+ 4H + (
3CO2 + 2S 0 + 5H2O → 3CH2O + 2SO4 2-+ 4H + (Findlay et al., 2017) (6.2) (iii.
3) biologically, via chemolithotrophic sulfide oxidation (colorless sulfur bacteria) using O2 or NO3 -as electron acceptors:

HS -+ NO3 -+ H + + H2O → SO4 2-+ NH4 + (7)
Or 4HS -+ NO3 -+ 6H + → 4S 0 + NH4 + + 3H2O [START_REF] Seitaj | Cable bacteria generate a firewall against euxinia in seasonally hypoxic basins[END_REF] (7.1) and

2S 0 + 3O2 + 2H2O → 2SO4 2-+ 4H + (Seitaj et al., 2015) (7.2)
In our model, the oxidation of H2S takes place either via processes 3,4,5,6 or 7 (with SO4 2-as an end product) or via processes 3.1, 6.1 or 7.1 (with S 0 or SO3 2-as an end product). While the rate of all of these four processes is denoted Fox, each of the four processes implies a distinct sulfur isotope fractionation, as defined in Table S3.

If not reoxidated into sulfate, the produced intermediate sulfur compound (S° or SO3 2-) can be involved in:

IV. disproportionation of S 0 or SO3 2-producing sulfate and sulfide, with rates Fdo and Fdh, respectively: (Bak and Cypionka, 1987;Thamdrup et al., 1993) (8) Or 4SO3 2-+ 2H + → H2S + 3SO4 2- (Bak and Cypionka, 1987) (

4S 0 + 4H2O → 3H2S + SO4 2-+ 2H +
Finally, aqueous sulfide produced by microbial reactions can be extracted from the sulfur cycling of marginal basins by V. solid iron sulfide formation and burial (Fcrs):

HS -+ Fe 2+ → FeS + H + (Eldridge et al., 2021) (10.1) HS -+ (x-1)S 0 → Sx 2-+ H + [START_REF] Eldridge | Theoretical estimates of equilibrium sulfur isotope effects among aqueous polysulfur and associated compounds with applications to authigenic pyrite formation and hydrothermal disproportionation reactions[END_REF] (10.2) FeS + H2S → FeS2 +H2 [START_REF] Eldridge | Theoretical estimates of equilibrium sulfur isotope effects among aqueous polysulfur and associated compounds with applications to authigenic pyrite formation and hydrothermal disproportionation reactions[END_REF] (11.1) FeS + Sx 2-→ FeS2 + Sx-1 2- [START_REF] Eldridge | Theoretical estimates of equilibrium sulfur isotope effects among aqueous polysulfur and associated compounds with applications to authigenic pyrite formation and hydrothermal disproportionation reactions[END_REF] (11.2)

Fixed fractionation factors

Sulfur isotope fractionation takes place during the biogeochemical processes introduced above as well as gypsum and pyrite burial fluxes. Fractionations factors α and λ are assigned fixed values (see Table S3), except for MSR for which isotope fractionation is a function of the cellular rate of microbial sulfate reduction (see section 4.2). In table S3, Rx are the isotope ratios ( x S/ 32 S, where x = 34 or 33) of the sulfur fluxes involved in the fractionation processes: Roc is the sulfur isotope ratio of the sulfate input from the ocean, Rriv is the sulfur isotope ratio of sulfate input from the rivers, Rmb is the sulfur isotope ratio of sulfate in the marginal basin, Rgyp is the sulfur isotope ratio of gypsum precipitated in the marginal basin, Rmsr is the sulfur isotope ratio of H2S produced by MSR, Rox is the sulfur isotope ratio of S° produced by oxidation of H2S, Rdo is the sulfur isotope ratio of SO4 2-produced by disproportionation of S° (or SO3 2-) and

Rdh is the sulfur isotope ratio of H2S produced by disproportionation of S°. 33 S and 34 S isotopes and for the different species of solid sulfide possibly generated through aqueous sulfide reaction with iron, it is, in our case, difficult to assess the exact nature of the iron sulfide preserved in our samples that we extracted chemically. Moreover, it is most likely that we are in presence of different types of solid sulfide compounds possibly formed by different processes. Therefore, and given that the fractionation factors remain very small and are most likely to have a negligible effect on the isotopic signature of the iron sulfide in our samples [START_REF] Eldridge | Theoretical estimates of equilibrium sulfur isotope effects among aqueous polysulfur and associated compounds with applications to authigenic pyrite formation and hydrothermal disproportionation reactions[END_REF][START_REF] Fry | Discrimination between 34S and 32S during bacterial metabolism of inorganic sulfur compounds[END_REF], we consider that the fractionation during solid sulfide burial as negligible (see Table S3).

Abiotic oxidation : Because there are no known fractionation factors for abiotic sulfide oxidation by iron and manganese oxides, we considered, like (Pellerin et al., 2015) , that these processes fractionate sulfur isotopes similar to abiotic oxidation of sulfide by molecular oxygen [START_REF] Eldridge | Theoretical estimates of equilibrium sulfur isotope effects among aqueous polysulfur and associated compounds with applications to authigenic pyrite formation and hydrothermal disproportionation reactions[END_REF].

Variable fractionation factor for microbial sulfate reduction

Multiple pure culture and environmental studies have found an inverse correlation between csSRR and the magnitude of isotopic fractionation (see section 4.3.2 of the main text for further discussion on the use of csSRR to model the range of variations expressed during MSR). [START_REF] Leavitt | Influence of sulfate reduction rates on the Phanerozoic sulfur isotope record[END_REF], applied a nonlinear regression model (equations ( 12) and ( 13 (

) 𝑐𝑠𝑆𝑅𝑅 = JK( == f `*+,-.*/ a =D f `b1 =D f `+c a =D f `b1 ) d == e ),. 12
where 34 ƐMeasured and 33 λMeasured represent the extent of fractionation measured for each individual SO4 2--H2S couple during the pure culture experiment, whereas 34 ƐMin , 33 λMin , 34 ƐMax and 33 λMax are the minimum and maximum empirical fractionations limits and 34 kmsr and 33 kmsr are pseudo first order rate constants obtained by fitting equations ( 12) and ( 13) to the experimental data. Leavitt et al., 2013, however, explored experiments with a variety of electron donors and sulfate concentrations, observing a wider range of isotope fractionations that approach the thermodynamic equilibrium between SO4 2- and H2S equal to 70‰ (Brunner and Bernasconi, 2005;[START_REF] Farquhar | Multiple sulphur isotopic interpretations of biosynthetic pathways: implications for biological signatures in the sulphur isotope record[END_REF][START_REF] Tudge | THERMODYNAMIC PROPERTIES OF ISOTOPIC COMPOUNDS OF SULPHUR[END_REF].

In our model, in order to cover the full range of isotopic fractionations that are possibly expressed by MSR activity, we expressed the cellular MSR rate of the Sim et al., 2011aSim et al., , 2011b experiments on a normalized scale, where csSRR= 0 corresponds to no sulfate reduction and csSRR=100 corresponds to the rate of sulfate reduction where the minimum isotope fractionation is expressed. We then refitted the data of Sim et al., 2011aSim et al., , 2011b to equations ( 14) and ( 15) (Fig. S2). This allowed us to obtain our own 34 ƐMin (-69) , 33 λMin (0.511) , 34 ƐMax (-5) , 33 λMax (0.515), 34 kmsr (0.0514) and 33 kmsr (0.1409) values.

In order to obtain the relationship that link 34 Ɛ and 33 λ to the normalized csSRR we rearranged equation ( 13).

88 𝜆 ghi = (𝑒 dkhl((× == e ),. (𝜆 gmn -𝜆 g:o ) + 𝜆 g:o (14

)
qq 𝛼 ghi = ( qs 𝛼 ghi )

qq E ),.

(15)

Model equations

As described in section 4.3.1 of the main text, the model is based on the chemical and isotopic mass balance for 34 S and 33 S sulfur isotope for SO4 2-,S 0 and H2S boxes (Main text Figure 6).

For the SO4 Rx are the isotope ratios ( x S/ 32 S, where x = 34 or 33) of the same fluxes.

To obtain output fluxes of gypsum-and pyrite-S expressed as a percentage of the sulfur fluxes added to the marginal basin from external sources, we define:

𝐹 :o + 𝐹 i:w = 100 (18)
Furthermore, we define parameters (fx, varying from 0 to 1) that express the relative magnitude of the different sulfur fluxes. Thus, fex represents the fraction of the incoming sulfur flux that is not involved in microbial reactions or gypsum precipitation:

𝑓 •n = € •-0 (€ b1 '€ .bƒ ) (19)
friv designates the fraction of the external sulfur inputs that is brought to the marginal basin by rivers:

𝑓 i:w = € .bƒ (€ b1 '€ .bƒ ) ( 20
)
and fmsr is the fraction of the net incoming sulfur flux being removed from the marginal basin through microbial sulfate reduction rather than by gypsum accumulation:

𝑓 ghi = € ),. (€ b1 '€ .bƒ d€ •-0 ) (21)
For the H2S box, the two mass balances are:

𝐹 ghi + 𝐹 u" = 𝐹 zy + 𝐹 vn (22)
and 𝐹 ghi × 𝑅 ghi + 𝐹 u" × 𝑅 u" = 𝐹 zy × 𝑅 zy + 𝐹 vn × 𝑅 vn (23)

Fpy is the flux of sulfur being buried as pyrite, Fox is the flux of sulfur being reoxidised to sulfate or intermediate sulfur compound, Fdh is the flux of sulfur being disproportionated to sulfide.

For the H2S box, we define fox as the fraction of the incoming sulfur fluxes that leaves the H2S box through H2S -oxidation (reactions 3, 4, 5 or 6) instead of being precipitated as pyrite:

𝑓𝑜𝑥 = € •c (€ ),. '€ / ‡ ) (24)
For the S° box, the two mass balances are:

𝐹 vn = 𝐹 u" + 𝐹 uv (25)
And

𝐹 vn × 𝑅 vn = 𝐹 u" × 𝑅 u" + 𝐹 uv × 𝑅 uv (26)
For the S° box, we define fdh as the fraction of sulfur added to the S° box by oxidation of H2S that is re-reduced to H2S through disproportionation:

𝑓𝑑ℎ = € / ‡ € •c (27)
Fluxes associated to biogeochemical transformations or mineral formation have specific fractionation factors 34 α and 33 λ defined in Table S3.

We define isotope fractionation factors for the considered processes that express a S-isotope fractionation:

𝛼 xyz = ( Š‹OE ( )• (28)
is the fractionation factor associated to the precipitation of gypsum,

𝛼 ghi = ( ),. ( )• (29)
is the fractionation factor associated to microbial sulfate reduction

𝛼 zy = ( OE‹ ( Ž•• (30)
is the fractionation factor associated to pyrite formation

𝛼 vn = ( •c ( Ž•• (31)
is the fractionation factor associated to H2S oxidation

𝛼 u" = ( / ‡ ( •° (32)
is the fractionation factor associated to reduction of S° to H2S through S° disproportionation

𝛼 uv = ( /• ( •° (33)
is the fractionation factor associated to oxidation of S° to SO42-through S° disproportionation

The set of 18 equations (22 to 33) and 18 unknowns (Fin,Fox,Friv,Fout,Fgyp,Fpy,Fmsr,Fdh,Fdo,RSO4,Rgyp,Rpy,RH2S,Rox,Rmsr,Rdo,Rdh,Rs°) was solved analytically with the Wolfram Mathematica 11 software, providing mathematical expressions for all sulfur fluxes and their isotope ratios (equations 37 to 54, below).

To facilitate reading of equations ( 39) to (56), we define:

𝐴 = 𝛼 xyz (𝑓 •n -1)(1 + 𝑓 ghi (𝑓 vn -1) -𝑓 u" 𝑓 vn ) (34) 𝐵 = -𝑓 •n + 𝛼 ghi 𝑓 ghi (𝑓 •n -1)(𝑓 u" 𝑓 vn -1) (35)
𝐶 = 𝛼 zy (𝑓 vn -1) (36) 𝐷 = (𝑓 u" 𝑓 vn -1)(𝑅 vk (𝑓 i:w -1) -𝛼 i:w 𝑅 i:w ) Following these new notations, solutions of equations ( 22) to (35) are given below

𝐹 :o = -100(𝑓 i:w -1) (37)
𝐹 i:w = 100𝑓 i:w (38)

𝐹 v{| = 100 × 𝑓 •n (39) 𝐹 ghi = -100𝑓 ghi (𝑓 •n -1) (40)
𝐹 vn = - € ),. -•c -/ ‡ -•c d- (41)
𝐹 u" = 𝐹 vn 𝑓 u"

𝐹 uv = -𝐹 vn (𝑓 u" -1) (43)

𝐹 xyz = - € •c (-),. '-/ ‡ -•c d-),. -•c d-) -),. -•c (44) 𝐹 zy = € ),. (-•c d-) (-/ ‡ -•c d-) (45)
𝑅 l˜s = ™[dL / ‡ -/ ‡ ›'L /• (-/ ‡ d-)(›dL •c -•c )] [L / ‡ -/ ‡ ›(•dž)] ' L /• (-/ ‡ d-)[(ž›'L •c -•c -*c (-/ ‡ -•c d-)d•(›dL •c -•c )] (46) 𝑅 ghi = 𝑅 l˜Ÿ 𝛼 ghi (47) 𝑅 xyz = 𝑅 l˜Ÿ 𝛼 xyz (48) 𝑅 ¡l = ™L ),. (-/ ‡ -•c d-)[L /• (-/ ‡ d-)dL / ‡ -/ ‡ ] [L / ‡ -/ ‡ ›(•dž)] ' L /• (-/ ‡ d-)[(ž›'L •c -•c -*c (-/ ‡ -•c d-)d•(›dL •c -•c )] (
𝑅 l°= -

( Ž•• L •c [L /• (-/ ‡ d-)dL / ‡ -/ ‡ ] (54)
To this set of equations were added equations ( 12), ( 14) and ( 15) that allowed to explore the full range of MSR fractionation.

Fixed Parameters

The isotopic composition of riverine input was set equal to the average values of modern rivers (δ 34 SSO4=4.8‰ (Burke et al., 2018) and Δ 33 SSO4=-0.0187‰ (Tostevin et al., 2014)). Based on the observation by Topper et al. (2011) that riverine inputs contributed for 10 to 80% of the total water inputs during PLG formation and on the fact that freshwater inputs are approximately 100 times less concentrated than oceanic waters in sulfate ion (Burke et al., 2018), we calculated that the fraction of incoming sulfate originating from riverine inputs on the total sulfate input, friv, ranged between 0 and 0.075, we fixed it here at 0.035. This is consistent with the blue mixing curve displayed on figure 4 that describe the mixing between oceanic sulfate and riverine inputs and from which can be inferred that the riverine inputs had a negligible influence on the sulfate isotopic composition of gypsum.

We then fixed the external basin sulfate isotopic composition equal to that of the sulfate from the first gypsum bed of the Vena del Gesso basin (δ 34 SSO4=21.78‰ and Δ 33 SSO4=0.033‰). This isotope composition is very close to that of sulfate dissolved in the global Messinian ocean (δ 34 SMESS=22.16±0.3 ‰ and Δ 33 SMESS= = 0.035 ±0.013 ‰ (Masterson et al., 2016)). Hydrochemical models of the MSC suggest that during deposition of PLG a two-way (inflow-outflow)

water exchange always existed between the Mediterranean Sea and the Atlantic Ocean, and between the marginal basins and the main Mediterranean basin (Topper et al., 2014) The main differences from our best fit results (where fex =0.9) concern the sulfate isotopic values, with the area covered being smaller when the choosen fex is larger, possibly requiring a larger fmsr to explain the whole range of values covered by our sulfate data (see supplementary section 6 for further explanation on this point).

The amount of intermediate sulfur compound being disproportionated into sulfide fdh is the complement of fdo, the amount of intermediate sulfur compound being disproportionated into sulfate. fdh was chosen depending on the scenario tested, with fdh=0 when the model was tested for the absence of reoxidation or disproportionation, fdh=0.25 when we tested the model for disproportionation of sulfite and fdh=0.75 when the elemental sulfur disproportionation scenario was choosen (see Pellerin et al., 2015 and references therein for stoichiometry of disproportionation reactions). Finally, the entire range of values for the fraction of sulfur reduced through MSR (0 < fmsr < 1) was tested in order to find the fmsr that would best fit our data (see supplementary information).

Constrains on fmsr

In the 34 ɛ-33 λ space, model results depend only on the magnitude of the isotope fractionation associated to MSR and on the fluxes (fox and fdh) and fractionation factors (αox, αdh and αdo)

controlling the reoxidative part of the sulfur cycle. Therefore, the 34 ɛnet-33 λnet dataset (built considering samples that have a paired sulfate-iron sulfide isotope composition), is used to constrain fox and ƐMSR. Most of the 34 ɛnet-33 λnet data is reproduced by the model only if fox > 0.85 and ƐMSR > -50‰ (Main text Fig. 8b).

Instead, in the δ 34 S-Δ 33 S space, the model results are a function of all sulfur fluxes and fractionation factors, and thus also of fmsr (Fig. 8a). Thus, with fox and ƐMSR constrained as explained above, fmsr can be constrained considering the δ 34 S-Δ 33 S dataset (Fig. S3). The sulfate isotope composition of VdG samples spans a small range from 20.7 to 24.1 ‰ in δ 34 S and from -0.001 to 0.0049 in Δ 33 S. In the δ 34 S-Δ 33 S space, the size of the area corresponding to the full considered range of fox and ƐMSR depends on fmsr: small values of fmsr produce a limited isotope reservoir effect, and thus a small deviation of model results from the initial basin sulfate composition, while elevated fmsr values result in a strong isotope reservoir effect and a large area defined by the model results (Fig. S3). With the constraints on fox and ƐMSR introduced above (fox > 0.85, ƐMSR >-50‰), only values of fmsr > 0.8 result in a good fit of the model to the sulfate data in the δ 34 SSO4-Δ 33 SSO4 space (Fig. S3).

Supplementary discussion on the reoxidative portion of the sulfur cycle

In the best fit model run (Main text Fig. 8), the use of fractionation factors and sulfur fluxes typical of elemental sulfur redox processing is in good agreement with the occurrence of iron sulfides associated with the remains of CSB in our samples, interpreted to result from early diagenesis of elemental sulfur. However, the nature of the microbial process driving the reoxidative portion of this cryptic sulfur cycle is still uncertain. The best fit model run is based on the combination of sulfide oxidation to S° by chemolitotrophs biofilms and subsequent microbial S° disproportionation (main text Fig. 8), because multiple sulfur isotope fractionation for these processes have been determined (Supplementary Table S3). However, the available isotope data for pure culture of sulfide-oxidizer is scant and we cannot exclude that other re-oxidative processes involving the S° intermediate produce a comparable isotope fractionation. In particular, the abundance of microbial fossils indicating the presence of CSB mats in Vena del Gesso gypsum makes these microorganisms the best possible alternative. To our knowledge, there are no experimental studies of sulfur fractionation by pure cultures of CSB. However, CSB being chemolithotrophs, and the biofilms from which comes the fractionation factors we used for our best fit for the oxidation of H2S to S 0 being dominated by chemolithotrophs organisms, it is likely that CSB produce a similar isotopic fractionation than the one we used while oxidizing H2S. Note that the occurrence of a reoxidative sulfur cycling involving the fractionation factors for elemental sulfur disproportionation has been previously suggested in porewaters displaying similar higher 33 λnet than maximal λMSR and in the presence of Beggiatoa CSB (Pellerin et al., 2015). Thus, our data and modelling imply that the cryptic sulfur cycle was driven by MSR and a re-oxidative process isotopically consistent with H2S oxidation to S 0 and subsequent elemental S° disproportionation, possibly involving the metabolism of CSB.

Other biogeochemical scenarios

In order to find a best fit to our data, we tested several scenarios involving different pathways for sulfur redox cycling. In all the tested scenarios the full range of fox (the relative amount of sulfide being oxidated into sulfate or intermediate sulfur compound) and ƐMSR (the extent of isotopic fractionation during MSR) were explored.

The amount of intermediate sulfur compound being disproportionated into sulfide fdh is the complement of fdo, the amount of intermediate sulfur compound being disproportionated into sulfate. fdh was chosen depending on the scenario tested, as explained below.

The occurrence of MSR without subsequent reoxidation was tested with attributing a fox=0, fdh=0 and varying fmsr and ƐMSR (Fig. S5). In Δ 33 S-δ 34 S diagram whereas our iron sulfide data does not fit the model, SO4 data could be explained by low values of ƐMSR and fmsr. In the 34 ɛnet and 33 λnet diagram, the model produces a curve that does not fit our data. Note that in the absence of reoxidation the model in the 34 ɛnet and 33 λnet is not an area but a curve following the different values of csSRR. This shows that the value of fmsr used in the 34 ɛ-33 λ space has no influence on the isotopic fractionation expressed. For each of the subsequent scenarios different values of fmsr (from 0 to 1) were tested in order to find a best fit to our data. This allowed us to understand what fraction of sulfate from the basin has been reduced through microbial sulfate reduction activity during PLG formation.

The direct reoxidation of hydrogen sulfide to sulfate was tested by combining the fractionation factors αox for the reoxidative flux fox attributed to equations 3, 4 or 5 (inorganic oxidation) or 6 (phototrophic oxidation) to fdh=0 (implying fdo=1) and no fractionation for the disproportionation flux αdh (Fig. S6). In both Δ 33 S-δ 34 S and 34 ɛ and 33 λ the model does not cover the full range of values expressed in our samples, supporting our conclusions that the filamentous fossils found in gypsum does not corresponds to cyanobacteria remains (Panieri et al., 2008) but to CSB remains(Dela Pierre et al., 2015;Natalicchio et al., 2022;Schopf et al., 2012). Note that in the model that consider phototrophic oxidation of hydrogen sulfide, low

ƐMSR and fox values in the model approach the values obtained for iron sulfide in marls, which is in good agreement with periodic euxinia during marl deposition (Fig. S5). This support the observation by Sinninghe Damsté et al. (1995) that have found abundant lipid biomarkers indicative of phototrophic green sulfur bacteria within the marls. We weren't able to test the hypothesis of a complete oxidation of sulfide into sulfate by chemolithotrophic bacteria because no fractionation factors are so far available in the literature (see section 4.1 above).

The disproportionation of elemental sulfur into sulfide was tested by combining the fractionation factors αox for the reoxidative flux attributed to equations 6.1 (phototrophic oxidation) or 7.1 (chemolithotrophic oxidation) to fdh=0.75 (implying fdo=0.25) (see the best fit Fig. 8 in the main text) in agreement with the stoichiometry of the reaction (Thamdrup et al., 1993) whereas the disproportionation of sulfite was tested by using fdh=0.25 (fdo=0.75) (Bak and Cypionka, 1987) (Fig. S7). As the fractionation factors of paired 33 S and 34 S isotopes for the oxidation of aqueous sulfide into sulfite remains unknown, when we tested our model for disproportionation of sulfite, the production of sulfite was done by testing all the fractionation factors available for aqueous sulfide oxidation (equations 3,4,5,6.1 and 7.1).

Whereas the areas covered by our model are larger than for previous scenarios (Fig. S7), allowing to fit a larger set of data point with the model, this scenario still doesn't reproduce the whole range of isotopic data we measured for the Vena del Gesso section. Moreover, the occurrence of remains of sulfur globules in our samples point its greater implication compared to sulfite as an intermediate sulfur compound in the genesis of the observed isotopic signatures of our samples. Mathematical relationship between the fluxes of sulfur Fx and the relative magnitude of these fluxes fx are given in section 5.

Figure 1 :

 1 Figure 1 : PLG in the Mediterranean and in the VdG basin.

Figure 2 :

 2 Figure 2: Macro-and microscopic characteristics of VdG PLG gypsum deposits. A. The basal part of cycle PLG5, composed of a basal dark marly interval followed by a thin whitish carbonate layer, in turn overlain by a gypsum bed consisting of dm-sized vertically oriented selenite crystals (massive selenite lithofacies). B. Close-up of the banded selenite facies (cycle PLG7). C. Close-up of a vertically oriented selenite crystal (massive selenite lithofacies). Note the turbid portion in the re-entrant angle of the twin. D. Photomicrograph (transmitted light) of a twinned selenite crystal showing the lamination in the re-entrant angle of the twin. E. Photomicrograph (transmitted light) of a gypsum twin with curved filaments aligned to the vertical growth bands. Lamination in the reentrant angle is also shown. F. Detail of (E) with filaments aligned along a vertical growth band. G, H, I. Photomicrographs in transmitted light (G and H) and UV light (I) showing the filaments inside gypsum. Pink arrows show tiny iron sulfide globules within the filaments. Note in (I) the high autofluorescence of the filaments. J, K, L. Photomicrographs in transmitted (J, K) and reflected (L) light of clay-rich aggregates in the turbid laminae of gypsum crystals. The aggregates are interpreted as marine snow floccules and contain abundant sulfide grains (arrows).

Figure 3 :

 3 Figure 3: Stable sulfate isotopic composition of VdG gypsum.

Figure 4 .

 4 Figure 4. Multiple Sulfur isotopic composition of sulfate and sulfide minerals constituting PLG deposits.

Figure 5 :

 5 Figure 5: Comparison between fractionation factors determined in laboratory experiments and net fractionations deduced from VdG sulfide and sulfate minerals.

Figure 6 :

 6 Figure 6: Structure of the numerical model simulating the biogeochemical cycling of sulfur in an idealized marginal basin. The three boxes correspond to the inventories of sulfate (SO42-), sulfide (H2S) and sulfur with an intermediate oxidation state (sulfite: SO32 or elemental sulfur: S°). The arrows represent the sulfur fluxes described in section 4.4.1.

  5.2.1. Petrographic evidence: co-occurrence of MSR and sulfide oxidation Petrographic analyses of gypsum samples show that iron sulfides are associated with either marine snow floccules (Figure 2.J-L) or filamentous fossils (Figure 2.I-G), which are dominant components of the VdG gypsum

  Figure 2. D-F). The mm-thick turbid laminae entrap filaments of CSB in close association with marine snow flocules rich in organic matter and iron sulfide minerals, (Figure 2. D-L). This suggests that sulfate reduction and sulfide reoxidation by CSB bacteria occurred in the same space, most probably in the mat adhering to the previously formed gypsum laminae (Figure 2. D-F, Figure 7. B,C). Since the surface underlying each successive layer of microbial mat was gypsum, the bottom boundary of the MSR zone was likely restricted to the base of the microbial mat, and sulfides were produced in-situ with no contribution of hydrogen sulfide diffusing from below (Figure7C).

  through MSR. It is important to note that the particular depositional conditions under which the VdG iron sulfide and gypsum minerals co-precipitated are drastically different from what is observed in traditional marine sediments settings ((Figure 7.B, C).

Figure 7 .

 7 Figure 7. Schematic representation of the biogeochemical conditions that accompanied PLG deposition. (A) Deposition of organic rich marls during humid phases (precession minima). (B) Deposition of gypsum in the margins and marls in the depocenter during arid phases (precession maxima). The relative size of the arrows represents the approximate magnitude of fluxes. (C) Details of a gypsum twin.e

( 0 .

 0 85<fox<1) (Figure 8.).We suggest that the relatively high ƐMSR and fox values inferred for gypsum compared to lower values for marls indicates the occurrence of alternations between drastically different biogeochemical environments associated with periods of marl and gypsum deposition (Figure 7 and 8.).

Figure 8 .

 8 Figure 8. Modelled isotope composition of VdG deposits compared to measured values. Data points in (A) represent the Δ 33 S vs δ 34 S of sulfate and sulfide of gypsum (massive, banded, and branching selenite) and marl samples from the PLG unit of the VdG Basin. Data points in (B) represent the isotopic fractionation calculated for paired sulfate and sulfide in the various gypsum lithofacies (massive, banded, and branching selenite). Individual curves represent the modelled isotopic compositions in (A) and net fractionations factors in (B) for different values of fox (solid lines) and ƐMSR (dashed lines) the extent of isotopic fractionation during MSR. ƐMSR approaches -70‰ when thermodynamic equilibrium is reached and fox=1, indicating that all the sulfide produced by MSR is reoxidized. Note that in both diagrams most of the data fall in the area defined by fox>0.85, indicating the occurrence of nearly complete reoxidation of sulfide, and supports the occurrence of ongoing cryptic sulfur cycling during gypsum formation. Most of the data fall in the area defined by ƐMSR <-50‰, indicating isotopic fractionation that deviated from thermodynamic equilibrium. (Error bars are 1σ).

  interpreted to corresponds to Messinian gypsum in other basins of the Mediterranean, even though a precise tuning and dating have not been realized on these sections (see the orange dots and green dots in main text Fig 1.c). They consist of PLG and resedimented PLG i.e Resedimented Lower Gypsum or Upper Gypsum deposits (formed at the end of the Messinian Salinity Crisis). This includes gypsum drilled during DSDP/IODP oceanographic expeditions LEG 13, 42 and 107 (see the green dots in main text Fig 1.c)

  )) on the results of a pure culture experiment with microbial sulfate reducer Desulfovibrio vulgaris hildenborough (DvH) under constant sulfate concentration and variable concentration of the electron donor lactate:

  only part of the potential fractionations possibly yielded by MSR communities because he only used lactate as an electron donor and set the sulfate concentration constant. On the other hand, Sim et al. (2011) carried out pure culture

Figure S2 :

 S2 Figure S2: Nonlinear regression model fits using the equations from Leavitt et al., 2013 to major and minor sulfur isotope fractionations measured in experiments by Sim et al 2011 a and b. In both panels A et B paired sulfate-sulfide major 34 Ɛmsr or minor 33 λmsr observed fractionation factors are plotted against the recast csSRR rate of MSR activity.

  and Fout are the exchange fluxes with the ocean, Friv is the riverine flux, Fgyp is the gypsum precipitation flux, Fmsr is the flux of sulfur reduced by MSR, Fdo is the flux of sulfur being disproportionated sulfate.

Figure S3 :

 S3 Figure S3: Sensitivity test for varying values of fex of 0.96 (green area) and 0.85 (brown area).

Figure S4 :

 S4 Figure S4: Figure 7: Measured and modelled Δ 33 SSO4 vs δ 34 SSO4 best fit with different values of

Figure S5 :

 S5 Figure S5: Measured and modelled Δ 33 SSO4 and Δ 33 SIS vs δ 34 SSO4 and δ 34 SIS data in A and 34 ɛnet

Figure S6 :

 S6 Figure S6: Measured and modelled Δ 33 SSO4 and Δ 33 SIS vs δ 34 SSO4 and δ 34 SIS data in A and 34 ɛnet

Figure S7 :

 S7 Figure S7: Measured and modelled Δ 33 SSO4 and Δ 33 SIS vs δ 34 SSO4 and δ 34 SIS data in A and 34 ɛnet

Table S1 :

 S1 Concentrations and multiple sulfur isotopic composition of Iron sulfide from theVena del Gesso gypsum and marly interbeds. Error is 1σ.

		Massive selenite Branching selenite	Ma III-01 Br X-04	0.046 0.074	8.88 -18.19	111 148
	4 11	Marl Marl	Pe IV-01 Pe XI-01	0.132 0.143	-35.37 -32.00	7614 6362
	4 11	Massive selenite Massive selenite	Ma IV-02 Ma XI-02	0.025 0.021	0.34 -0.26	236 191
	4 11	Banded selenite Banded selenite	Ba IV-05 BA XI-03	0.120 0.149	-32.56 -33.83	283 69
	5 5 11 12	Marl Massive selenite Branching selenite Marl	Pe V-01 Ma V-03 Br XI-04 Pe XII-03	0.142 0.018 0.099 0.133	-33.29 -11.16 -22.51 -34.31	8791 270 116 5543
	5 12	Banded selenite Massive selenite	Ba V-06 Ma XII-04	0.071 0.026	-27.89 -1.11	225 240
	6 6 12 6 12 6 13 7 13	Marl Massive selenite Banded selenite Banded selenite Branching selenite Branching selenite Banded selenite Marl Branching selenite	Pe VI-01 Ma VI-02 Ba XII-05 Ma VI-03 Br XII-06 Br VI-04 Ba XIII-03 Pe VII-01 Br XIII-04	0.161 0.037 0.127 0.082 0.128 0.034 0.127 0.119 0.016	-32.54 -13.90 -32.09 -21.86 -34.80 -12.74 -37.12 -31.12 -12.84	8590 112 76 234 251 502 236 1815 166
	7 14	Laminated gypsum Massive selenite	La VII-02 Ma XIV-02	0.120 0.036	-36.35 6.06	355 117
	7 14	Massive selenite Banded selenite	Ma VII-03 Ma XIV-03	0.068 0.023	-33.08 -4.98	111 219
	7 14	Banded selenite Branching selenite	Ba VII-04 Br XIV-04	0.034 0.061	-9.26 -29.46	35 460
	7 15 7 15 8 15	Branching selenite Marl Branching selenite Massive selenite Marl Banded selenite	Br VII-05 Pe XV-01 Ba VII-06 Ma XV-02 Pe VIII-01 Ba XV-03	0.054 0.121 0.126 0.029 0.125 -0.008	-17.83 -36.86 -27.16 -13.96 -32.08 -14.78	2699 97 74 351 1155 80
	8 15	Massive selenite Branching selenite	Ma VIII-02 Br XV-04	0.040 0.100	-26.00 -36.01	240 454
	8 16 8 16 9 16	Banded selenite Marl Branching selenite Massive selenite Marl Massive selenite	Ba VIII-03 Pe XVI-01 Br VIII-04 Ma XVI-02 Pe IX-01 Ma XVI-02 bis	0.086 0.014 0.089 0.031 0.130 0.009	-27.46 -17.84 -25.71 -2.73 -30.26 0.34	10197 52 89 29 2306 98
	PLG Cycle 9 16	Facies Massive selenite Banded selenite	Sample name Ma IX-02 Ba XVI-05	Δ 33 SIS ±0.008‰ 0.061 0.027	δ 34 SIS ±0.4‰ -21.34 -20.13	[IS] ±100ppm 147 117
	1 9 Iron sulfide concentrations calculated from our chromium reducible compound chemical Massive selenite Ma I-01 0.045 -3.63 15 Banded selenite Ba IX-03 0.122 -30.47 103
	2 10 2 10 3 10	Massive selenite Marl Massive selenite Massive selenite Marl Banded selenite	Ma II-02 Pe X-01 Ma II-07 Ma X-02 Pe III-01 Ba X-03	-0.017 0.132 0.025 0.045 0.143 0.047	-18.06 -31.15 0.36 -20.69 -33.87 -21.33	6190 37 171 101 13317 145

extractions range from 15 to 502 ppm in selenite samples and from 1155 to 13117 ppm in marls (Table

Table S2 :

 S2 Oxygen and multiple sulfur isotope composition of dissolved sulfate ion from the Vena del Gesso gypsum. Error is 1σ.

			Banded selenite Banded selenite	MT IV-05 Ba XIV-03	0.033 0.026		22.15 22.37	8.01 10.20
	5 14		Massive selenite Branching selenite	Ma V-03 Br XIV-04	0.028 0.038		21.17 21.89	8.43 10.08
	5 15		Banded selenite Branching selenite	Ba V-06 Br XV-04	0.037 0.043		21.83 21.99	8.53 10.11
	6 16		20 Massive selenite Massive selenite	Ma VI-02 Ma XVI-02	0.036 0.033		massive 21.92 21.88	9.24 n.d
	7 7 7 7 6 6 16	-20 -10 0 10 selenite Branching δ 34 Laminar selenite Massive selenite Banded selenite S (‰) vs CDT selenite Branching Banded selenite Massive selenite	Br VII-05 La VII-02 Ma VII-03 Ba VII-04 Br VI-04 Ma VI-03 Ma XVI-02 bis	0.023 0.033 0.043 0.038 0.037 0.034 0.043		23.47 21.38 23.68 21.85 banded pelite 23.05 branching 22.13 21.98 10.11 8.66 n.d n.d 9.05 8.80 8.88
	7		-30 Branching selenite	Br VII-06	0.026		20.73	8.76
	8 8		-40 Massive 1 selenite Banded selenite	10 Ma VIII-02 100 [Iron sulfide] (ppm) 1 000 0.039 Ba VIII-03 0.029	10 000 22.00 22.81	8.33 100 000 9.36
	8		Branching selenite	Br VIII-04	0.043		22.15	9.01
	9		Massive selenite	Ma IX-02	0.039		22.17	9.01
	9		Banded selenite	Ba IX-03	0.040		22.19	9.66
	Cycle 10 1 10 2 11 2 11 3 12 4 13 14	Facies Massive selenite Banded selenite Massive selenite Banded selenite Massive selenite Branching selenite Massive selenite Massive selenite Massive selenite Branching selenite Massive selenite Massive selenite	Sample name Ma X-02 Ba X-03 Ma I-01 Ba XI-03 MC II-02 Br XI-04 Ma II-07 Ma XII-04 Ma III-01 Br XIII-04 Ma IV-02 Ma XIV-02	Δ 33 SSO4 ±0.004‰ -0.001 0.049 0.033 0.039 0.037 0.036 0.033 0.039 0.029 0.030 0.029 0.020	δ 34 SSO4 ±0.2‰ 24.12 24.05 21.78 22.44 22.65 21.64 22.10 22.12 22.06 22.87 22.39 22.18	9.99 δ 18 OSO4 ±0.2‰ 9.29 8.24 9.81 9.26 9.39 n.d 10.39 8.67 10.08 8.32 10.09

Table S3 :

 S3 Sulfur fractionation factors used in the numerical model

	Reaction	Process	34 α	33 α	34 ε (‰)	33 λ	Reference
	Gypsum
	precipitation 𝒙 𝑹 𝒈𝒚𝒑	Equation (1)	1,00165 1,00085	1,65	0,515	(Thode and Monster, 1965)
	𝒙 𝑹 𝒎𝒃
	Solid Sulfide Burial 𝒙 𝑹 𝒑𝒚 𝒙 𝑹 𝑯𝟐𝑺	Equations (10), (11.1) and (11.2)	1	1	0	0,515	(Eldridge et al., 2021; Fry et al., 1986)*
		Inorganic
		H2S to SO4 2-					(Eldridge and
		Equations (3),	0,99440 0,99715	-5,6	0,508	Farquhar,
		(4) or (5)					2018)*
	Sulfide Oxidation 𝒙 𝑹 𝒐𝒙	Phototrophic H2S to SO4 2-Equation (6)	1,00154 1,00081	1,5	0,529	(Zerkle et al., 2009)
	𝒙 𝑹 𝑯𝟐𝑺	Phototrophic H2S to S 0 Equation (6.1)	1,00189 1,00099	1,89	0,535	(Zerkle et al., 2009)
		Chemolithotrophic H2S to S 0 Equations (7.1)	0,99600 0,99795	4	0,513	(Zerkle et al., 2016)
	Intermediate	Phototrophic S° to SO4 2-Equation (6.2)	0,99810 0,99906	-1,9	0,495	(Zerkle et al., 2009)
	compound	Disproportion					(Johnston et al.,
	Oxidation	S° to SO4 2-	1,01853 1,009584	18,5	0,519	2007; Johnston
	𝒙 𝑹 𝒅𝒐	Equation (8)					et al., 2005)
	𝒙 𝑹 𝑺𝒊	Disproportion					(Johnston et al.,
		SO3 2-to SO4 2-	1,01000 1,00527	10	0,528	2005; Pellerin
		Equation (9)					et al., 2015)*
	Intermediate Compound Reduction 𝒙 𝑹 𝒅𝒉 𝒙 𝑹 𝑺𝒊	Disproportion S° to H2S Equation (8) Disproportion SO3 2-to H2S Equation (9)	0,99382 0,99680 0,95500 0,97672	-6,2 -45	0,516 0,511	(Johnston et al., 2007; Johnston et al., 2005) (Johnston et al., 2005)*
	*see next paragraph
	Solid Sulfide Burial : Although Eldridge et al. (2021) have calculated fractionation factors for
	both

  2-box, the two mass balances are: 𝐹 :o + 𝐹 uv + 𝐹 i:w = 𝐹 xyz + 𝐹 ghi + 𝐹 v{| (16) and 𝐹 :o × 𝑅 vk + 𝐹 uv × 𝑅 uv + 𝐹 i:w × 𝑅 i:w = 𝐹 xyz × 𝑅 xyz + 𝐹 ghi × 𝑅 ghi + 𝐹 v{| × 𝑅 g}

  . At the acme of gypsum precipitation in a marginal basin connected to the Mediterranean sea,Topper and Meijer (2015) calculated that calcium sulfate concentration in the marginal basin is approximately 1.2 times greater than in the main Mediterranean basin, whereas the water outflow represents 69 to 77% of the water inflow from the main Mediterranean basin. This allowed us to estimate that the fraction of sulfate output from the Vena del Gesso basin to the Mediterranean represented between 87 and 96% of the incoming sulfate fluxes; thus, we set fex equal to 0.9 to represent sulfate exchanges with the external sulfate source during gypsum precipitation.

Supplementary Material

1. Isotopic measurements 1.1 Multiple sulfur isotopes Sulfate from gypsum and chromium reducible sulfides (iron sulfide) from gypsum and marls were reduced into H2S using two different reducing solutions following Jovovic et al. (2020).

Released H2S was trapped in a NaOH solution as HS -, transformed into Ag2S by adding a AgNO3 solution, purified with 3-5 drops of 68%HNO3 to dissolve co-precipitating Ag2O impurities [START_REF] Geng | A simple and reliable method reducing sulfate to sulfide for multiple sulfur isotope analysis[END_REF], rinsed three times with Milli-Q H2O and finally dried in an oven at 60°C. Ag2S powder was then reacted in Ni-reaction bombs with F2 at 250°C overnight to produce gaseous SF6 that was purified cryogenically and by gas chromatography before to be analyzed by a dual inlet mass spectrometer MAT-253 at Institut de Physique du Globe de Paris [START_REF] Geng | A simple and reliable method reducing sulfate to sulfide for multiple sulfur isotope analysis[END_REF]Ono et al., 2006). The δ 34 S values were calibrated with an in-house SF6