Microscale C-O-S isotope signatures record microbial sulfate reduction in ancient deep subsurface pyrite–calcite–barite assemblages

The deep biosphere harbors diverse microbial communities amounting to a substantial fraction of Earth's biomass for much of its history. Analyses of fracture-hosted minerals (e.g., pyrite, calcite) and their C-O-S isotope signatures are beginning to establish a fossil record of subsurface microbial communities and their relationship with major geological events. Fracture-coating pyrite from deep within the Fennoscandian Shield records the largest ³⁴S-enrichment observed on Earth to date (δ³⁴S = +147 ‰) – likely the result of late-stage Rayleigh distillation during closed-system microbial sulfate reduction (MSR). This implies even heavier S isotope values for the complementary sulfate reservoir during pyrite formation, possibly recorded in coeval sulfate minerals such as barite [BaSO₄]. However, barite has been poorly explored as an archive of ancient deep subsurface biosignatures, particularly in samples with extreme ³⁴S-enrichment in pyrite. Here, we compiled published microscale δ¹³C_calcite, δ¹⁸O_calcite, and δ³⁴S_pyrite data with new secondary ion mass spectrometry (SIMS) δ¹⁸O_barite and δ³⁴S_barite analyses from two localities on the Fennoscandian Shield (Forsmark and Laxemar, Southeastern Sweden), aiming to constrain the record of microbial processes in barite. We identified two groups of barites distinguished by their distinct ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends. Comparison of these trends with δ¹³C values of cogenetic calcite reveals variable associations with sulfate-dependent anaerobic oxidation of methane (AOM) and organotrophic sulfate reduction (OSR). Strongly ¹³C-depleted calcite (δ¹³C_calcite < −40 ‰), indicating a relationship with AOM, almost exclusively occurs in the steeper-trending Group 1. Evidence for barite affected by OSR (δ¹³C_calcite > −30 ‰) occurs in both groups and dominates the shallower-trending Group 2. The different trends likely resulted from an interplay of MSR pathways (AOM vs. OSR) and variations in sulfate reduction rates (SRR), which we attribute to the relative abundances of sulfate and electron donors in the parent fluid of barite. Lower sulfate and electron donor abundances likely permitted methanogenesis and AOM at lower SRR, resulting in a steeper ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trend in Group 1 barite. Conversely, higher sulfate and electron donor abundances inhibited methanogenesis and AOM, favoring OSR at a higher SRR, which is reflected by the shallower ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trend in Group 2 barite. Comparison between the δ³⁴S_pyrite range (−53.9 ‰ to +131.7 ‰) and the δ³⁴S_barite range (+7.6 ‰ to +52.0 ‰) across both localities demonstrates that the sulfate reservoir corresponding to extremely ³⁴S-enriched pyrite is not recorded in barite. The discrepancy between barite and pyrite could reflect their formation in separate microenvironments or barite undersaturation at high degrees of Rayleigh distillation. Our study highlights the need for coupled microscale C-O-S analyses in fracture-hosted mineral assemblages to understand metabolic processes in the ancient deep biosphere, while stressing that these records are affected by local hydrogeochemical conditions.

1 Introduction

The average depth of the 122 °C isotherm in the continental crust is 4.8 km, defining a vast rock volume that is cooler than the upper known temperature limit of life (Takai et al., 2008). For this reason, fractures and rock pores in the continental crust represent one of the largest habitable environments on Earth (Magnabosco et al., 2018). This “deep biosphere” hosts diverse communities of chemolithoautotrophic primary producers and heterotrophic microorganisms, including methanogens, methanotrophs, and sulfate reducers (Chivian et al., 2008; Kieft et al., 2018; Ruff et al., 2024; Stevens and McKinley, 1995; Wu et al., 2016). Recent estimates of the total biomass of the deep biosphere vary (Bar-On et al., 2018; Magnabosco et al., 2018; McMahon and Parnell, 2014), but its significant impact on global biogeochemical cycles is generally accepted. Before the rise of land plants, the deep subsurface may have even held the majority of Earth's total biomass (McMahon and Parnell, 2018). Indeed, portions of the continental crust have been continuously habitable since 1.1 billion years ago (Ga) and may have experienced habitable periods since 1.9 Ga (Drake and Reiners, 2021). Understanding the diversity and biogeochemical impact of the deep biosphere throughout Earth's history requires robust and geologically stable biosignatures. However, the deep biosphere fossil record remains patchy and poorly understood (Gustafsson and Drake, 2025; Ivarsson et al., 2020).

Minerals and their stable isotope ratios can represent valuable biosignatures in deep time (Runge et al., 2023). For instance, the largest magnitudes of S isotope fractionation between sulfate and sulfide (³⁴ϵ = δ³⁴S_sulfate − δ³⁴S_sulfide) are associated with kinetic isotope effects (KIE) during microbial sulfate reduction (MSR). Depending on the microbial strain and environmental parameters, such as the temperature, or the availability of sulfate and electron donors (organic matter or H₂), ³⁴ϵ induced by MSR can vary between 0 ‰ and ∼ 70 ‰ (Bradley et al., 2016; Chambers et al., 1975; Farquhar et al., 2003; Habicht et al., 2002; Ries et al., 2009; Sim et al., 2011; Wing and Halevy, 2014; Wortmann et al., 2001). These signatures can be recorded in ³⁴S-depleted authigenic sulfide minerals (e.g., pyrite) (Halevy et al., 2023; Pasquier et al., 2025; Strauss, 1997). Microbial sulfate reduction can also be coupled with anaerobic oxidation of methane (AOM), which is well supported in the deep subsurface (Bell et al., 2022; Drake et al., 2015, 2017; Kotelnikova, 2002). Microbially produced methane is ¹³C-depleted relative to CO₂ by up to 95 ‰ (Whiticar, 1999), hence methane oxidation via AOM can result in the precipitation of strongly ¹³C-depleted calcite (Aloisi et al., 2002; Campbell et al., 2002; Feng et al., 2010; Lin et al., 2011; Natalicchio et al., 2012; Reitner et al., 2005). Calcite precipitation resulting from organotrophic sulfate reduction (OSR) typically reflects the less ¹³C-depleted organic carbon substrate. Notably, dating of fracture-coating minerals (adularia $\mathrm{Rb}/\mathrm{Sr}$, calcite $\mathrm{U}/\mathrm{Pb}$, and hematite $\mathrm{U}\text{-}\mathrm{Th}/\mathrm{He}$) co-genetic to biosignatures in pyrite and calcite records ancient MSR deep within the Fennoscandian Shield between the Ordovician (460 ± 30 Ma) and the Neogene (13.7 ± 0.6 Ma) (Drake et al., 2017, 2018, 2021). Thus, combining S and C isotope ratios is a powerful approach for reconstructing ancient MSR pathways.

During the early stages of MSR, δ¹⁸O_sulfate increases via a minor KIE and an equilibrium isotope effect (EIE) between H₂O and intermediate S species (e.g., sulfite) before their re-oxidation to sulfate (Antler et al., 2013; Böttcher et al., 1999; Brunner et al., 2005; Fritz et al., 1989; Mizutani and Rafter, 1973; Wortmann et al., 2007) (Fig. 1). The slopes of ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends are strongly controlled by the sulfate reduction rate (SRR). At a lower SRR, the equilibration of intermediate S species with water is more pronounced, and δ¹⁸O_sulfate rises more rapidly relative to δ³⁴S_sulfate (i.e., a steeper slope). At a higher SRR, the rise in δ¹⁸O_sulfate relative to δ³⁴S_sulfateis less pronounced (i.e., a shallower slope). In addition, AOM and OSR can result in different slopes at similar SRR (Antler et al., 2014), and the slopes resulting from AOM differ between methane diffusion-limited settings (e.g., marine sediments: steeper trend) and methane-in-excess environments (e.g., cold seeps, stratified estuaries: shallower trend) (Antler et al., 2015) (Fig. 1). Owing to the preferential consumption of ${}^{\mathrm{32}}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ by microorganisms, MSR in a closed system (i.e., when the SRR exceeds the resupply of sulfate) can produce ³⁴S-enriched sulfate in the residual porewater (Anderson et al., 1987; Jørgensen et al., 2004; Rudnicki et al., 2001; Zaback et al., 1993). This process has resulted in δ³⁴S_pyrite values of up to +147 ‰ in deep fractures of the Fennoscandian Shield, the most extreme ³⁴S-enrichment yet reported (Drake et al., 2018). This implies an even more ³⁴S-enriched complementary sulfate reservoir during pyrite formation. Barite [BaSO₄] is a robust archive of δ¹⁸O_sulfate and δ³⁴S_sulfate values in its parent fluid because its precipitation is associated with minimal stable isotope fractionation and it is thermodynamically stable and poorly soluble (Crockford et al., 2019; Griffith and Paytan, 2012; Han et al., 2025; Paytan et al., 1998, 2002; Turchyn and Schrag, 2004, 2006; Wood et al., 2021). Thus, the co-variation of δ¹⁸O_barite and δ³⁴S_barite values may record past and extant environmental conditions and microbial S cycling, complementing the C-O-S isotope records in pyrite and calcite. However, pyrite–calcite–barite assemblages are poorly explored archives of ancient deep subsurface microbial processes, particularly in samples containing extremely ³⁴S-enriched pyrite.

Figure 1 Controls on ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends during the early stages of microbial sulfate reduction (MSR). The intracellular reduction of sulfate to sulfide (HS⁻) proceeds via multiple reversible steps, including the activation of adenosine 5^′-phosphosulfate (APS) and reduction to sulfite (${\mathrm{SO}}_{\mathrm{3}}^{\mathrm{2}}$). The kinetic isotope effect (KIE) associated with MSR produces ³⁴S-depleted sulfide and ³⁴S-enriched residual sulfate. Concerning O, the KIE associated with MSR is small, and the equilibrium isotope exchange (EIE) between sulfate and water is kinetically inhibited. However, intermediate S species such as sulfite readily exchange O isotopes with water and become relatively ¹⁸O-enriched. Intracellular re-oxidation of this sulfite produces ¹⁸O-enriched sulfate. The ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends in residual sulfate are controlled by the expression of the EIE associated with sulfite-water exchange relative to the sulfate reduction rate (SRR). This expression becomes more pronounced when the SRR is lower, or when anaerobic oxidation of methane (AOM) occurs under methane (CH₄) limited conditions, resulting in a steeper slope. After Wortmann et al. (2007), Antler et al. (2013, 2014, 2015), Han et al. (2025).

Here, we compile previously published microscale δ¹³C_calcite,δ¹⁸O_calcite, and δ³⁴S_pyrite data from multiple deep drill cores taken at two localities on the Fennoscandian Shield (Forsmark and Laxemar/Äspö, Southeastern Sweden) with new secondary ion mass spectrometry (SIMS) δ¹⁸O_barite and δ³⁴S_barite analyses from the same drill cores. In most cases, barite occurred in fractures that contained associated pyrite and/or calcite. Although we did not detect extreme ³⁴S-enrichment in barite, our results demonstrate that barite records an ancient archive of OSR and AOM, while highlighting that these archives are affected by local hydrogeochemical conditions.

2 Methods

2.1 Materials and sampling

The sampled barite-bearing fractures are from Proterozoic crystalline bedrock at Forsmark, Laxemar, and Äspö (Southeastern Sweden) (see Drake et al. (2018) for details on the geological setting). Laxemar and Äspö are situated on the Baltic Sea coastline, and only ca. three kilometers apart; hence we consider them as one site (i.e., Laxemar) in this study. Drill cores were mainly obtained by triple-tube drilling conducted by the Swedish Nuclear Fuel and Waste Management Co. This drilling technique preserves fragile minerals on fracture surfaces, facilitating the investigation of ancient biological activity in deep fracture systems. Barite was sampled from fractures in a total of 40 drill cores. The samples represented depths ranging from 60 to 1660 m below sea level. The exact locations of the boreholes are provided in Runge et al. (2026). Sample characterization and mineral identification were performed directly on the fracture surfaces using a FEI QUANTA FEG 650 Scanning Electron Microscope (SEM) equipped with an integrated energy-dispersive spectroscopy (EDS, Oxford T-Max 80 detector) under low-vacuum conditions. The acceleration voltage was 20 kV, and the instrument was calibrated with a cobalt standard. The barite grains were hand-picked from the fracture surface under a microscope, embedded in several epoxy mounts (2.5 cm diameter), and polished to expose cross-sections of the grains. The polished grains were examined using SEM prior to SIMS analysis to detect growth zonation and avoid impurities and cracks.

Table 1 Stable isotope data (δ³⁴S, δ¹⁸O, δ¹³C) of fracture-hosted minerals (barite, pyrite, calcite) from Forsmark and Laxemar (Fennoscandian Shield, Sweden). Barite data were analyzed in this study. Calcite and pyrite data were obtained from Drake et al. (2017, 2018).

2.2 Secondary-ion mass spectrometry (SIMS) analysis of barite (δ¹⁸O and δ³⁴S)

Secondary-ion mass spectrometry (SIMS) was conducted on a CAMECA IMS1280 ion probe at the NordSIMS lab, Swedish Museum of Natural History, Stockholm. The analytical routines for S and O isotope analyses of barite follow those described by Drake et al. (2024), which are based on the approaches by Whitehouse (2013) for S and Heinonen et al. (2015) for O. The samples were sputtered using a ¹³³Cs⁺ primary beam with 20 kV incident energy and a beam current of ∼ 3 nA, which was rastered over a 5 µm × 5 µm area during analysis. Each analysis comprised 40 s of pre-sputtering to remove the gold coating over a 15 µm × 15 µm area, centering of the secondary beam in the field aperture, and data acquisition in 12 four-second integration cycles. Secondary ion signals (³²S and ³⁴S or ¹⁶O and ¹⁸O) were detected simultaneously using two Faraday detectors with a common mass resolution of 4860 ($M/\mathrm{\Delta }M$) for S and 2430 ($M/\mathrm{\Delta }M$) for O. Data were normalized for instrumental mass fractionation using the S0327 reference material, with a conventionally determined δ³⁴S value of +22.0 ± 0.3 ‰ and δ¹⁸O value of +11.0 ± 0.1 ‰ (Liseroudi et al., 2021). All results are reported as δ³⁴S_V−CDT (Ding et al., 2001) for S and δ¹⁸O_V−SMOW (Coplen, 1995) for O. The typical uncertainties after propagating the within-run and external uncertainties from the standard measurements were ± 0.2 ‰ (1σ) for δ³⁴S and between ± 0.2 ‰ and ± 0.6 ‰ (1σ) for δ¹⁸O values. The previously published calcite δ¹⁸O_V−PDB values were converted to δ¹⁸O_V−SMOW for comparability with barite using the formula δ¹⁸O_V−SMOW = 1.03092 × δ¹⁸O_V−PDB + 30.92. All SIMS data, including previously published data, are compiled in Runge et al. (2026).

2.3 Statistical analysis

Statistical analysis of 366 paired ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{barite}}/{\mathit{\delta }}^{\mathrm{34}}{\mathrm{S}}_{\mathrm{barite}}$ data was conducted using a Gaussian Mixture Model (GMM), which is a useful tool for classifying multivariate stable isotope data (Göhring et al., 2016). Notably, GMM is highly sensitive to outliers; hence, they need to be excluded before GMM to prevent adverse effects on the clustering performance. We identified outliers in the paired δ³⁴S–δ¹⁸O data using density-based spatial clustering of applications with noise (DBSCAN; Ester et al., 1996). DBSCAN identifies data points that do not belong to any dense region as noise without assuming a predetermined number or shape of clusters, making it useful for outlier detection in bivariate isotope datasets (Valentine et al., 2015). The parameters for DBSCAN are the minimum number of points (minPts) required to form a dense region and the neighborhood radius (ϵ). We tested the sensitivity of our outlier detection to these parameters over a range of minPts values (2 to 15). For each minPts value, ϵ was set at the point of maximum curvature of the sorted k-distance curve (k = minPts) (Ester et al., 1996) (Fig. S1 in the Supplement). The ${\mathit{\delta }}^{\mathrm{34}}\mathrm{S}/{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ data were converted to z-scores before clustering to ensure neither variable dominated ϵ, and DBSCAN was run for each parameter pair (Fig. S2 in the Supplement). The stability of the identified parameter combinations was quantified using the Jaccard index (the number of outliers identified in two neighboring minPts settings divided by the number identified in each), and the number of dense clusters was tracked across settings (Fig. S3 in the Supplement). The number of clusters was stable at 2 for minPts ≥ 6, and neighboring minPts settings from 13 to 14 identified identical points as outliers (Jaccard = 1.00), indicating that DBSCAN using minPts values of 13 and 14 yields stable outlier detection results. Therefore, we used minPts = 13 and the associated ϵ = 0.586. This parameter combination identified 21 of 366 data points as outliers that were excluded before GMM.

Figure 2 Depth distribution and violin plots of stable isotope data for fracture-hosted minerals from Forsmark. (A): δ³⁴S_barite, δ³⁴S_pyrite. (B): δ¹⁸O_barite, δ¹⁸O_calcite. The δ¹⁸O_calcite data were converted from V-PDB to V-SMOW for comparison with δ¹⁸O_barite. (C): δ¹³C_calcite. Each data point represents one SIMS analysis spot ($\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ = meters above sea level). Barite data were analyzed in this study. Calcite and pyrite were obtained from Drake et al. (2017, 2018).

We performed GMM on the outlier-cleaned data set using the Python module Scikit-learn (Pedregosa et al., 2011). To select the appropriate number of clusters for GMM in our ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{barite}}/{\mathit{\delta }}^{\mathrm{34}}{\mathrm{S}}_{\mathrm{barite}}$ dataset, we identified the best-fit model by calculating the Bayesian Information Criterion (BIC) as a function of the number of clusters. The minimum BIC value indicated our dataset contained two clusters (Fig. S4A in the Supplement). Random state optimization for GMM was performed over a range of zero to ten, yielding a random state of zero at the minimum BIC. York regression was used to calculate the slopes and r² values of the resulting clusters (Fig. S4B) (York et al., 2004). Notably, our dataset contains data points from two localities (i.e., Laxemar and Forsmark), whereas Laxemar provides the dominant fraction (323 vs. 21 for S). Deriving robust and meaningful clusters requires consistent behavior of these clusters across localities. Hence, we tested the sensitivity of our GMM approach to differences between localities by excluding Forsmark and repeating the same GMM procedure described for the entire dataset. After excluding the Forsmark data, the BIC-guided decision for a two-cluster model remained robust (Fig. S4C). Additionally, the resulting slopes of the individual clusters changed minimally after excluding Forsmark data (Fig. S4B and D). This suggests our clusters were independent of the localities. Thus, we used the cluster assignment and slopes resulting from GMM of the whole dataset to interpret our dataset. In concert, these observations strongly support a two-cluster GMM approach for classifying our ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{barite}}/{\mathit{\delta }}^{\mathrm{34}}{\mathrm{S}}_{\mathrm{barite}}$ dataset. Nevertheless, the resulting clusters should not be regarded as globally representative and may differ for barite from other areas.

Figure 3 Depth distribution and violin plots of stable isotope data for fracture-hosted minerals from Laxemar. (A): δ³⁴S_barite, δ³⁴S_pyrite. (B): δ¹⁸O_barite, δ¹⁸O_calcite. The δ¹⁸O_calcite data were converted from V-PDB to V-SMOW for comparison with δ¹⁸O_barite. (C): δ¹³C_calcite. Each data point represents one SIMS analysis spot ($\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ = meters above sea level). Barite data were analyzed in this study. Calcite and pyrite were obtained from Drake et al. (2017, 2018).

3 Results and discussion

3.1 Pyrite–calcite–barite assemblages record ancient microbial sulfate reduction in the deep subsurface

Coupled δ³⁴S_pyrite and δ¹³C_calcite signatures provide evidence of ancient organotrophic sulfate reduction (OSR) and sulfate-dependent anaerobic oxidation of methane (AOM) in deep subsurface fractures at Forsmark and Laxemar (Drake et al., 2013, 2015, 2018, 2021). The compiled data from Forsmark and Laxemar showed a δ³⁴S_pyrite range of −53.9 ‰ to +131.7 ‰, a δ¹⁸O_calcite range of +10.5 ‰ to +28.1 ‰, and a δ¹³C_calcite range of −93.1 ‰ to +36.5 ‰ (Drake et al., 2017, 2018) (Table 1, Figs. 2 and 3). A prominent role of thermochemical sulfate reduction (TSR) in these localities has been ruled out because of the large apparent ³⁴ϵ between plausible δ³⁴S_sulfate values of fracture fluids in these systems and pyrite (55 ‰ to 66 ‰: Drake et al., 2021). Moreover, TSR generally occurs at temperatures > 100 °C (Machel, 2001), which is inconsistent with fluid inclusion homogenization temperatures and clumped isotope thermometry indicating formation temperatures < 100 °C of coeval calcite (Drake et al., 2018; Herlambang et al., 2023). Thus, the wide δ³⁴S_pyrite range (−53.9 ‰ to +131.7 ‰) was interpreted to reflect Rayleigh distillation during microbial sulfate reduction (MSR) in semi-closed fracture systems, producing ³⁴S-depleted sulfide and extremely ³⁴S-enriched sulfate (Drake et al., 2018). The overall ranges of δ¹⁸O_barite (−0.9 ‰ to +19.7 ‰) and δ³⁴S_barite (+7.6 ‰ to +52.0 ‰) across both localities are consistent with MSR, which can produce ¹⁸O- and ³⁴S-enrichment in the residual sulfate pool (Antler et al., 2013, 2014, 2015; Brunner et al., 2005; Drake et al., 2018; Han et al., 2025; Jørgensen et al., 2004; Rudnicki et al., 2001; Wortmann et al., 2007; Zaback et al., 1993). These observations strongly suggest MSR could be recorded in co-genetic barite from Forsmark and Laxemar.

The extremely high δ³⁴S_pyrite values at Forsmark and Laxemar imply the presence of an even more ³⁴S-enriched sulfate reservoir potentially archived in co-genetic barite. To detect this archive and investigate its relationship with OSR and AOM, we analyzed δ³⁴S_barite and δ¹⁸O_barite in the same drill cores that were previously sampled for δ³⁴S_pyrite and δ¹⁸O_calcite, and δ¹³C_calcite analysis (Drake et al., 2017, 2018) (Figs. 2 and 3; Table 1). The δ³⁴S_sulfate values of the paleo-fracture fluids are unknown and likely changed over time. However, the δ³⁴S_sulfate values of the current fracture fluids (+14 ‰ to +27 ‰: Drake et al., 2013) are similar to the median δ³⁴S_barite values at Forsmark (+26.6 ‰) and Laxemar (+18.3 ‰), suggesting dissolved sulfate in such fluids is a plausible source for barite. The total range of δ³⁴S_barite (+7.6 ‰ to +52.0 ‰) is consistent with processing of such fracture fluids by MSR under closed-system conditions (Paytan et al., 2002). However, this range is much narrower than that of δ³⁴S_pyrite (−53.9 ‰ to +131.7 ‰) across Forsmark and Laxemar (Figs. 2 and 3; Table 1). The ancient barite from Forsmark and Laxemar does not record the extreme ³⁴S-enrichment expected from the δ³⁴S_pyrite data. This agrees with barite that recently precipitated on drill core equipment within a deep borehole at Äspö (Fichtner et al., 2024), and carbonate-associated sulfate from the Forsmark and Laxemar drill cores (Fichtner et al., 2026). Therefore, a record of extremely ³⁴S-enriched sulfate in the deep subsurface remains undetected.

Figure 4 Gaussian Mixture Modeling (GMM) of barite δ¹⁸O-δ³⁴S data and correlation with published calcite δ¹⁸O-δ¹³C data. Each data point represents one spot with paired SIMS δ¹⁸O and δ³⁴S analyses. (A): δ¹⁸O and δ³⁴S analyses of barite and GMM-derived cluster assignments of data points. The two-component GMM of barite δ¹⁸O and δ³⁴S data after exclusion of outliers (black) yielded “Group 1” and “Group 2” barite. York regression yielded slopes of 0.95 for Group 1 and 0.55 for Group 2. The ranges of the regression lines equal the data range used for GMM and cannot be extended to the outliers. (B): δ¹⁸O and δ³⁴S analyses of barite color-coded by δ¹⁸O of calcite in the same sample. There was no apparent correlation between the Group 1 and Group 2 trends and their corresponding δ¹⁸O_calcite values, demonstrating that the trends arise from fractionation processes affecting the sulfate reservoir, rather than fluid mixing. (C): δ¹⁸O and δ³⁴S analyses of barite color-coded by δ¹³C_calcite in the same sample. The δ¹³C_calcite data suggest that AOM barite occurs almost exclusively in Group 1. In contrast, OSR-related barite occurred in Groups 1 and 2. The color coding in (B) and (C) follows recommendations for scientifically derived color maps (Crameri et al., 2020).

Table 2 Statistical parameters of the groups resulting from the Gaussian Mixture Model.

3.2 Variable sulfate reduction rates and pathways explain microscale δ¹⁸O-δ³⁴S-δ¹³C data

The co-variation of the δ¹⁸O_barite and δ³⁴S_barite values exhibited substantial scatter (Fig. 4A), implying a mixture of different sulfate sources and/or an overlap of different fractionation processes. Since the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ systematics of dissolved sulfate in the low-temperature groundwater systems at Forsmark and Laxemar are affected by OSR and AOM (Drake et al., 2017, 2018; Fichtner et al., 2024), we aimed to identify whether our dataset contains groups of barite associated with either OSR or AOM. To identify barite associated with these groups, we classified our dataset using Gaussian Mixture Modeling (GMM) and compared the linear regression slopes of the resulting clusters with known sulfate ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends in MSR-affected settings. We identified a steeper trending “Group 1” cluster (slope: 0.95 ± 0.01; r² = 0.39) and a shallower trending “Group 2” cluster (slope: 0.55 ± 0.01; r² = 0.38) (Fig. 4A; Table 2). The clusters resulting from GMM are statistical entities and cannot be interpreted as continuous fractionation trends in our study. This is because our dataset contains multiple barite growth generations from different locations and records different fractionation trends, non-coeval sulfate reservoirs between growth zones, or local fluid mixing. Therefore, our statistical analysis can only detect average ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends of overlapping fractionation or mixing processes involving fracture water sulfate that was processed by OSR and AOM at different points in space and time. These overlapping trends likely explain the low r² values. There was no correlation between Group 1 or Group 2 barite and the δ¹⁸O values of co-genetic calcite (Fig. 4B). This suggests that the difference between Group 1 and Group 2 does not simply reflect isotopically distinct water sources (Pedersen et al., 1997). Instead, they most likely record different processes affecting the respective sulfate reservoirs from which barite precipitated.

Previous studies established that sulfate reduction rates (SRR) control the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ slope in residual sulfate during the early stages of MSR (Brunner et al., 2005; Han et al., 2025; Wortmann et al., 2007). At a lower SRR, δ¹⁸O rises rapidly relative to δ³⁴S leading to a steeper slope, whereas at a higher SRR, δ¹⁸O rises slowly relative to δ³⁴S, leading to a shallower slope (Fig. 1). Thus, the difference in ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ slopes could reflect a lower SRR in Group 1 than in Group 2 barite (Fig. 4A). Moreover, differences in metabolic pathways (OSR vs. AOM) can lead to variable ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends in sulfate at similar SRR (Antler et al., 2013, 2014, 2015). To clarify the effects of SRR and MSR pathways on ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ systematics recorded in barite at Forsmark and Laxemar, we correlated our GMM-classified δ¹⁸O_barite and δ³⁴S_barite data (i.e., excluding outliers) with the δ¹³C values of calcite in the same samples (Figs. 2C, 3C, and 4C). When petrographic correlation of barite with individual calcite generations was not possible, we used average δ¹³C_calcite values from the barite-bearing samples. Out of all classified barite data, 71 spots in Group 1 and 81 spots in Group 2 had associated δ¹³C_calcite data. Strongly ¹³C-depleted calcite (δ¹³C < −40 ‰) was almost exclusively associated with Group 1 barite (28 out of 71) (Fig. 4C). The δ¹³C values of this calcite (−70 ‰ to −40 ‰) cannot be explained by oxidation of organic matter via OSR alone and suggest variable input of inorganic C from the oxidation of primarily microbial methane (Drake et al., 2015, 2017). This implies that some Group 1 barite associated with such calcite was affected by AOM. In contrast, only 2 out of 81 spots in Group 2 barite were associated with strongly ¹³C-depleted calcite. Instead, Group 2 was almost exclusively associated with calcite that has a δ¹³C ≈ −30 ‰ or higher (Fig. 4C), suggesting methane was not the main C source (Whiticar, 1999). This calcite likely reflects a more prominent contribution from organic matter rather than methane oxidation, implying that the associated barite was affected by OSR (Drake et al., 2017). Notably, such OSR-related barite not only dominates Group 2, but also widely occurs in the AOM-affected Group 1 (Fig. 4C). Thus, it seems unlikely that the different MSR pathways are the only control on the different slopes of both groups. Instead, we interpret the different ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends to result from an interplay of MSR pathways and SRR: Group 1 reflects AOM and OSR at a relatively low SRR, whereas Group 2 dominantly reflects OSR at a relatively high SRR.

The outliers identified via DBSCAN (n = 21) are strongly heterogeneous and do not follow a common ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trend (Fig. 4A, Table 2). Thus, the outliers do not form a coherent group and likely reflect highly variable fluid compositions and fractionation processes. Moreover, a relationship of the outliers with fluid sources and MSR pathways is unclear because only four outliers have associated δ¹⁸O_calcite and δ¹³C_calcite values (Fig. 4B and C).

3.3 Hydrogeochemical controls on sulfate reduction rates and pathways

The major control on SRR is the availability and quality of electron donors (i.e., H₂, CH₄, or organic matter) (Leavitt et al., 2013; Lyons and Gaudette, 1979; Westrich and Berner, 1984). SRR is high when electron donors are more bioavailable and abundant, and vice versa. The abundances of electron donors in the parent fluids from which our ancient barite precipitated are unknown. However, low SRR due to nutrient limitation was inferred based on the large apparent ³⁴ϵ associated with pyrite formation at these sites (Drake et al., 2018), consistent with the common free energy limitation in deep subsurface habitats in general (Hoehler and Jørgensen, 2013). Thus, a relatively lower availability of electron donors in Group 1 could have caused the lower inferred SRR than in Group 2 (Fig. 5A). Moreover, the typical range of sulfate concentrations in modern anoxic fracture fluids at Forsmark and Laxemar (< 1 to 10 mM: Hallbeck and Pedersen, 2008; Laaksoharju et al., 2008; Drake et al., 2013) is much lower than in marine environments (28 mM in seawater: Jørgensen et al., 2019). Such sulfate concentrations may additionally limit SRR in the deep subsurface (Boudreau and Westrich, 1984; Pallud and Van Cappellen, 2006). Sulfate concentrations in these subsurface habitats also vary around the threshold at which sulfate reducers outcompete methanogens for substrate utilization (ca. 1–10 mM: Iversen & Jorgensen, 1985; Glossner et al., 2016; Sela-Adler et al., 2017). If sulfate concentrations exceed this threshold, OSR can be favored, inhibiting methanogenesis and AOM. Thus, δ¹³C evidence for AOM indirectly supports sulfate limitation in Group 1 (Fig. 4c). Thus, a combined effect of electron donor limitation and low sulfate concentrations could have favored methanogenesis and AOM at relatively low SRR (Fig. 5A). Conversely, a higher abundance of electron donors could explain a higher inferred SRR in Group 2 than in Group 1 (Fig. 5B). At higher sulfate availability, methanogenesis is inhibited, possibly resulting in a lack of evidence for AOM in Group 2 (Fig. 5B). Instead, OSR is favored and SRR is high, resulting in a shallower trend in Group 2 (Fig. 5B). Thus, we interpret the relative abundances of sulfate and electron donors to control the observed ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends in ancient subsurface barite by modulating SRR and MSR pathways.

Figure 5 Conceptual model of the controls on δ¹⁸O_barite and δ³⁴S_barite systematics in the continental deep subsurface. (A): At low sulfate and electron donor abundances, microbial methanogenesis occurs and microbial sulfate reduction proceeds via both sulfate-dependent anaerobic oxidation of methane (AOM) and organotrophic sulfate reduction (OSR) at a relatively low sulfate reduction rate (SRR). This results in the steep ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{barite}}/{\mathit{\delta }}^{\mathrm{34}}{\mathrm{S}}_{\mathrm{barite}}$ trend recorded in Group 1 barite. (B): At higher sulfate and electron donor abundance, microbial methanogenesis is inhibited, and MSR proceeds predominantly via OSR at relatively high SRR. This results in the shallow ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{barite}}/{\mathit{\delta }}^{\mathrm{34}}{\mathrm{S}}_{\mathrm{barite}}$ trend recorded in Group 2.

The stronger ³⁴S-enrichment in the OSR-dominated Group 2 than in the AOM-affected Group 1 contrasts with observations from marine sediments, where ³⁴S-enriched barite is a commonly used proxy for AOM at the SMTZ (Torres et al., 1996; Wood et al., 2021; Yang et al., 2023; Zan et al., 2022; Zhou et al., 2015). This discrepancy could be due to the higher sulfate concentration in seawater (28 mM: Jørgensen et al., 2019) compared to the subsurface at Forsmark and Laxemar (< 1 to 10 mM), allowing sulfate to become strongly ³⁴S-enriched before barite becomes undersaturated (Fichtner et al., 2024). The relatively steep slope of the AOM-affected Group 1 is comparable to some previous observations from marine methane-diffusion-limited AOM settings, such as sedimentary SMTZs, or OSR-dominated methane-devoid settings (Antler et al., 2013, 2015; Blake et al., 2006; Böttcher et al., 2006). At the same time, Group 2 shows little clear evidence of AOM but a shallow slope similar to marine methane-in-excess environments (i.e., the methane flux exceeds methane consumption via oxidation), such as cold seeps (Aharon and Fu, 2000, 2003; Antler et al., 2013, 2014; Rubin-Blum et al., 2014). These observations highlight a further discrepancy between marine environments and the deep subsurface. In marine environments, steeper ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ slopes tend to be found in settings affected by OSR or methane-limited AOM, whereas shallower slopes result from AOM at high methane concentrations. In our subsurface barite, steeper slopes tend to be more associated with evidence for AOM, whereas shallower slopes tend to occur in OSR-dominated settings. The relatively shallow slope of OSR-related barite in the subsurface implies a lesser extent of equilibrium isotope exchange of S intermediates with water compared to marine settings. This could be achieved via a higher overall SRR in the subsurface than in marine settings; however, this seems unlikely given the oligotrophic nature of subsurface habitats (Küsel and Dopson, 2026). Additional factors influencing the different behavior of ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ trends in the deep subsurface may be the increasing temperature and pressure with depth or factors within microbial communities, such as low cell densities or adaptations to extreme nutrient limitation. However, the impact of these factors on ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ relationships in residual sulfate during MSR remains speculative. Irrespective of this unresolved problem, the differences in ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ systematics between marine and subsurface barite highlight that interpretations cannot be extended from one system to the other without additional constraints.

Figure 6 Backscatter electron micrographs of calcite-pyrite–barite assemblages in drill core samples investigated in this study. (A, B): Sample KAS02 (802 m depth) from Äspö (i.e., Laxemar). (A): Barite (Brt) intergrown with calcite (Cal), demonstrating the coeval formation of barite and calcite. Calcite in this sample has an average δ¹³C of −72.1 ‰, demonstrating calcite precipitation after the oxidation of biogenic methane via AOM (Drake et al., 2017). (B): Euhedral pyrite (Py) intergrown with calcite, demonstrating that pyrite is co-genetic with the outermost calcite generation. (C, D): Sample KLX01 (220 m depth) from Laxemar. (C): Barite intergrown with calcite, demonstrating their co-genetic relationship. (D): Pyrite intergrown with calcite, demonstrating the coeval formation of pyrite with calcite and, thus, barite. The calcite generation coeval with barite has a $\mathrm{Rb}/\mathrm{Sr}$ age of 394 ± 14 Ma and δ¹³C_calcite = −11.6 ‰ (Drake et al., 2017).

3.4 Barite and pyrite record distinct hydrogeochemical regimes

The absence of extreme ³⁴S-enrichment in barite from Forsmark and Laxemar is unexpected considering the corresponding δ³⁴S_pyrite values of up to +131.7 ‰, predicting an even more ³⁴S-enriched sulfate reservoir during pyrite formation. The δ³⁴S_barite values of up to +52.0 ‰ observed among the outliers suggest that the SRR at times exceeded the sulfate resupply, decreasing the steady-state sulfate concentration (i.e., closed-system conditions) (Fig. 4A). Hence, Rayleigh distillation is locally recorded in subsurface barite, although not to the extreme degree observed in pyrite (Table 1). This discrepancy between the barite and pyrite records has been suggested to result from (i) the precipitation of barite and pyrite under separate microenvironmental conditions, involving the isolation of evolved ³⁴S-enriched fluids, or (ii) late-stage barite undersaturation which may inhibit sulfate precipitation or favor its post-formation dissolution before the formation of a ³⁴S-enriched sulfate reservoir (Fichtner et al., 2024, 2026). Barite undersaturation may arise from sulfate depletion at high degrees of Rayleigh distillation (Fichtner et al., 2024; Wood et al., 2021). Alternatively, barite precipitation during earlier stages may deplete Ba²⁺ (Aloisi et al., 2004), which in crystalline bedrock fractures is mainly supplied by feldspar dissolution. However, none of the δ³⁴S analyses of carbonate-associated sulfate in the samples analyzed here detected this reservoir (Fichtner et al., 2026). This suggests that the mechanism limiting the precipitation of extremely ³⁴S-enriched sulfate is more likely associated with the S cycle, rather than cation availability. Sampling biases cannot be excluded but have been mitigated by the large number of δ³⁴S analyses (n = 505) of barite from multiple fractures at two localities. Notably, the strongest ³⁴S-enrichment observed here occurs in the outlier group and is not associated with cogenetic calcite data (Fig. 4A–C; Table 2). This makes an interpretation of the associated MSR pathway and SRR difficult for these cases. However, it seems likely that these barites represent relatively late-stage precipitates that possibly formed in an environment where barite saturation was less limiting.

Figure 7 Microscale stable isotope variations in individual barite and pyrite grains. (A): SEM image of a polished barite grain with three distinct growth zones in sample KAS02 (802 m depth) from Äspö (i.e., Laxemar). The white dotted lines highlight the boundaries between growth zones, whereas the white dashed lines highlight analytical transects. (B, C): SIMS analytical transects of barite and associated pyrite in the same sample (A). Transects across three pyrite crystals from the same sample show a δ³⁴S_pyrite variation of ca. 65 ‰ across two of the crystals. (D): SEM image of a polished barite grain with two distinct growth zones and SIMS analytical transects of barite and associated pyrite in sample KLX01 (220 m depth) from Laxemar. The overlain heatmap was interpolated via cubic convolution and does not reflect measured δ³⁴S values beyond the analyzed spots. (E, F): SIMS analytical transects of barite and associated pyrite in the same sample (D()). The transect across a zoned barite crystal (E) shows co-variation of δ³⁴S_bariteand δ¹⁸O_barite with higher ³⁴S- and ¹⁸O-enriched values in the later generation (rim). The transect across five pyrite grains shows a δ³⁴S_pyrite variation of ca. 20 ‰ across one grain (F). The color coding of the analysis spots in (A) and (D) reflects their assignment to Group 1 (blue) and Group 2 (orange) according to Fig. 4A. White dots could not be assigned to a group because corresponding δ¹⁸O_barite values were lacking. The error bars associated with pyrite analyses are smaller than the symbols.

Microscale observations of pyrite–calcite–barite assemblages could help investigate this interpretation further (Figs. 6 and 7). Sample KAS02 (802 m) can serve as an example of AOM-associated barite due to its association with calcite that has an average δ¹³C of −72.1 ‰ and pyrite with δ³⁴S values spanning −23.5 ‰ to +42.5 ‰ (Fig. 6A and B). Additionally, all SIMS spots on barite in this sample belong to Group 1 (Fig. 7A). Together, these observations are consistent with a sulfate pool affected by AOM at relatively low SRR. Notably, this barite exhibited distinct crystallographic and isotopic zonation with only minor δ³⁴S and δ¹⁸O variabilities between growth zones of ∼ 4 ‰ and ∼ 0.8 ‰, respectively (Fig. 7B). This low variability suggests relatively open system conditions during barite growth (Fichtner et al., 2024; Wood et al., 2021). In contrast, the much more pronounced core-to-rim δ³⁴S increase in the co-genetic pyrite from −23.5 ‰ to +42.5 ‰ indicates MSR in a more closed system, meaning that SRR exceeded the resupply of sulfate during pyrite growth (Fig. 7C). Apparently, barite and pyrite in this sample record precipitation in distinct hydrogeochemical regimes, possibly at different stages of fracture fluid evolution or in separate microenvironments.

Sample KLX01 (220 m) serves as an example of OSR-related barite due to its association with calcite lacking a substantial δ¹³C depletion (−11.6 ‰) (Fig. 6C and D). Notably, barite in this sample exhibited a massive core and a porous rim (Fig. 7D). These texturally distinct growth zones were associated with a sharp core-to-rim δ³⁴S and δ¹⁸O increase of ∼ 12 ‰ and ∼ 6 ‰, respectively (Fig. 7D and E). This sharp increase suggests that barite in the rim precipitated from sulfate that had already undergone ³⁴S-enrichment via sulfate reduction elsewhere in the fracture system (e.g., Drake et al., 2018). Additionally, SIMS spots related to Group 1 predominated in the core, whereas spots related to Group 2 exclusively occurred in the rim (Fig. 7D). This implies that the processing of the respective sulfate reservoirs reflected by core and rim occurred under different SRR and/or MSR pathways. Thus, the distinct textural and isotopic differences demonstrate that the generations reflect two barite precipitation events, rather than continuous barite growth. The porosity in the rim suggests barite dissolution (Fig. 7D), possibly due to undersaturation via rapid sulfate depletion at the high SRR inferred for Group 2 (Fichtner et al., 2024). Thus, a relatively high SRR in a sulfate-poor environment may result in a preservation bias against extremely ³⁴S-enriched barite in deep subsurface environments. The relatively minor core-to-rim δ³⁴S increase of ca. 0 ‰ to 20 ‰ in co-genetic pyrite implies that sulfate did not become strongly depleted during pyrite formation in this microenvironment (Fig. 7F). While there was no evidence for AOM from δ¹³C_calcite in this sample, it cannot be resolved petrographically whether the calcite and pyrite are co-genetic to the barite core or rim. Therefore, it remains unclear whether the shift from Group 1 to Group 2 during barite growth is due to different MSR pathways or SRR, or both. Nevertheless, the differences between the samples KAS02 (802 m) and KLX01 (220 m), as well as their evolution during barite growth, demonstrate a strong control of hydrogeochemical conditions on their isotopic composition in different microenvironments and over time. This highlights that the associated barite and pyrite record different environmental conditions and should be correlated with caution.

4 Conclusions

We compiled previously published microscale δ¹³C_calcite, δ¹⁸O_calcite, and δ³⁴S_pyrite data with new SIMS δ¹⁸O_barite and δ³⁴S_barite analyses to investigate biosignature records in ancient deep subsurface mineral assemblages from Forsmark and Laxemar/Äspö (Southeastern Sweden). We identified two groups of barites recording microbial sulfate reduction (MSR) associated with distinct ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/{\mathit{\delta }}^{\mathrm{34}}\mathrm{S}$ slopes. Comparison with coeval δ¹³C_calcite data revealed that evidence for sulfate-dependent anaerobic oxidation of methane (AOM) is mostly restricted to the steeper-trending Group 1, whereas Group 2 was dominated by organotrophic sulfate reduction (OSR). We interpret the different slopes to result from an interplay of variable sulfate reduction rates (SRR) and MSR pathways (AOM vs. OSR): Group 1 reflects AOM and OSR at relatively low SRR, whereas Group 2 dominantly reflects OSR at relatively high SRR. Variations in SRR and MSR pathways likely reflect the relative sulfate and electron donor abundances in the barite's parent fluids. Barite records δ³⁴S values up to +52.0 ‰, which is well below the extreme ³⁴S-enrichment observed in some pyrite from the same locality (δ³⁴S_pyrite = +131.7 ‰). The discrepancy between barite and pyrite records could reflect their formation in separate microenvironments or may be due to undersaturation caused by sulfate and/or Ba²⁺ depletion at high degrees of Rayleigh distillation. In concert, microscale C-O-S isotope signatures in fracture-hosted pyrite–calcite–barite assemblages at Forsmark and Laxemar record biosignatures of AOM and OSR deep within the Fennoscandian Shield since at least 394 Ma, but these records are affected by local hydrogeochemical conditions. Signatures of metabolic processes in barite varied over time within the same microenvironment and between different fractures. This highlights that understanding deep subsurface geobiology requires integrated microscale analyses of mineral assemblages distributed over large crustal volumes.