All chemicals used were of analytical grade and sourced from Merck KGaA (Darmstadt, Germany). De-ionized water (Millipore, Darmstadt, Germany) was utilized for preparing the mineral medium, the analytical standards, and all other solutions.

A bench-scale anaerobic bioreactor (total volume: 250 mL) was filled with 180 mL of mineral medium and inoculated with 20 mL of groundwater collected from a petrochemical-contaminated aquifer in Italy. The composition of the mineral medium is reported in Supplementary Table 1. Upon setup, the bioreactor was purged with pure nitrogen gas (N₂) to establish anaerobic conditions. Toluene was the sole carbon and energy source and periodically spiked to reach a final concentration of 0.15 mM at the beginning of each feeding cycle. Each feeding cycle ended when toluene was almost completely degraded. Sulfate was also periodically added to a final concentration of 2 mM per cycle as the terminal electron acceptor. The culture was maintained under at room temperature (22 ± 3 °C) under continuous mixing, with pH kept constant at 7.2 ± 0.3.

In order to confirm that the removal of toluene was coupled with the reduction of sulfate, microcosm experiments were set up in serum bottles (vol: 60 mL), containing 15mL of mineral medium and inoculated with 10mL of the microbial culture described in the previous paragraph after 150 days of enrichment. Oxygen was eliminated by N₂ purging, and toluene and sulfate were added at a final concentration of 0.15mM and 2mmol/L respectively. Two kinds of treatments were implemented: in the first type, sodium molibdate was added (final concentration of 20Mm), while in the second (i.e., control) no molibdate was added. Molybdate is a well-known selective inhibitor reduction (Isa, 2004). The microcosms were kept under constant mixing by means of a magnetic stirrer and incubated at room temperature (22 ± 3 °C). Throughout the whole study, the pH was kept constant at 7.2 ± 0.3. Sulfate and toluene concentrations were monitored throughout the experiment.

The presence of O₂, H₂, and CH₄ was monitored by analyzing gaseous samples with a gas-chromatograph (Agilent 8860, GC system USA) equipped with a thermal conductivity detector (TCD). Toluene concentrations were measured by analyzing gaseous samples with a gas-chromatograph (Agilent 8860, GC system USA) equipped with a flame ionization detector (FID). Acetate concentrations were determined by injecting 1 μL of filtered (0.2 μm) and acidified (10 % vol/vol with 0.3 M oxalic acid) liquid samples, into a gas-chromatograph (Agilent 8860, GC system USA) equipped with a flame ionization detector (FID). The methods used for GC analysis, are reported in the Supporting Information (Tab. S2). Sulfate concentrations were determined by injecting a liquid sample (filtered, 0.22 μm porosity) into an ion chromatograph (IonPac AS14 analytical column, Dionex DX-100 system, Dionex Corp., Sunnyvale, CA, USA).

Liquid sample (50 mL) of the culture were collected at steady state with a sterile syringe, filtered through hydrophilic polycarbonate membranes (0.2 μm pore size, 25 mm diameter, Millipore, Italy) and immediately processed for DNA extraction using the DNeasy PowerLyzer PowerSoil Kit (QIAGEN), following the manufacturer's protocol. DNA concentration (>100 ng/μL) and purity were assessed with the Qubit dsDNA HS Assay Kit and on the NanoDrop One device (both from Thermo Fisher Scientific, USA). DNA size distribution was analyzed using Genomic DNA ScreenTapes on the Agilent Tapestation 4200 (Agilent, USA). DNA was then used for long amplicon sequencing through 16S rRNA targeted and untargeted metagenomics. Amplicon libraries targeting the bacterial 16S rRNA gene were prepared using the Oxford Nanopore 16S Barcoding Kit (SQK-16S024, Oxford Nanopore Technology, ONT), according to the manufacturer's instructions. The genomic library was loaded onto a FLO-MIN106D flow cell (R9.4.1 chemistry) and sequenced using the MinION Mk1B device (Oxford Nanopore Technology, ONT). For shotgun metagenomics, barcoded SQK-NBD114.96 DNA libraries were prepared with minor modifications to the manufacturer's protocol (Oxford Nanopore Technologies, United Kingdom). The libraries were loaded onto primed FLO-PRO114M (R10.4.1) flow cells and sequenced on a PromethION P24 device running MinKNOW Release 24.02.19. Signal data was basecalled and demultiplexed using Dorado basecall server v. 7.3.11 (Oxford Nanopore Technologies) with the super-accurate (SUP) algorithm (dna_r10.4.1_e8.2_400bps_5khz_sup.cfg). Residual adapters and low-quality reads were trimmed using nanoq v. 0.10.0 with parameters -q 15, -l 1000, -S 150 and -E 150 to remove the first and last 150 bp of each read, along with any reads below q15 and length 1000 bp (Steinig, 2022).

The raw sequencing reads and metagenome-assembled genomes (MAGs) generated in this study have been deposited in NCBI under BioProject accession PRJNA1250843.

Reads from 16S rRNA long amplicon sequencing were basecalled and demultiplexed using MinKNOW Guppy Dorado 7.3.9, employing the super-accurate basecalling algorithm (configurations: dna_r10.4.1_400bps_sup.cfg and dna_r9.4.1_450bps_sup.cfg). Bacterial 16S rRNA reads were filtered by length and quality using Chopper v0.8.0 (https://github.com/wdecoster/chopper), applying a minimum quality score of 9 and selecting reads between 1400 and 1700 bp. Taxonomic assignment and abundance calculations were performed using EMU v3.4.5 (Curry et al., 2022) with the Silva NR99 database v138.2.

Metagenome assemblies were generated using flye v. 2.9.4 with the parameters –meta and –extra-params min_read_cov_cutoff=5 (Kolmogorov et al., 2020). The assemblies were polished once using medaka v. 1.11.3 (https://github.com/nanoporetech/medaka) and contigs shorter than 1,000 bp were removed using seqkit v. 2.8.2 (Shen et al., 2016).

Metagenome-assembled genomes (MAGs) were binned with SemiBin v. 2.1.0 (–self-supervised, –sequencing-type long_read, –minfasta-kbs 500) using a self-supervised contrastive learning algorithm (Pan and Zhao, 2023). Inter-assembly MAG dereplication was performed using the drep v. 3.5.0 dereplicate workflow (-comp 50 -con 9.99 -p -l 500000) alongside checkm2 v. 1.0.2 (Olm et al., 2017; Chklovski et al., 2023). Dereplicated MAGs with a minimum completeness of 50% and contamination below 10% were functionally annotated with DRAM v.0.1.2 via KBase (Shaffer et al., 2020).

The tool gtdbtk v. 2.4.0 was used to taxonomically classify MAGs against the Genome Taxonomy Database (GTDB) release 220 [https://doi.org/10.1093/bioinformatics/btz848]. Maximum likelihood phylogenomic tree was generated with GTDB-tk de-novo workflow through the identification, alignment, and concatenation of 120 phylogenetically informative marker genes (bac120). The resulting alignment file was finally analyzed in IQ-TREE (10.1093/molbev/msu300) to compute a maximum likelihood phylogenetic tree under the Whelan and Goldman (WAG) model with 1000 pseudo-bootstrap replicates (Aroney et al., 2025; Chaumeil et al., 2020; Nguyen et al., 2015).

The temporal dynamics of toluene and sulfate concentrations over 170 days of operation are shown in Figure 1. After an initial startup phase (feeding cycle 1), toluene degradation rates increased significantly and remained stable in subsequent feeding cycles, indicating successful acclimation of the microbial consortium. Throughout the experiment, toluene degradation was consistently coupled to sulfate reduction, confirming the role of SO4²⁻ as the terminal electron acceptor for anaerobic toluene oxidation. During the 8th feeding cycle, the toluene and sulfate loads were doubled to test the robustness of the system. As a result, this perturbation caused a temporary decline in degradation rates, most likely due to toluene toxicity effects. However, from cycle 9 onward, when the initial substrate concentrations were restored, both toluene biodegradation and sulfate reduction fully recovered, returning to pre-inhibition levels. Excluding the startup phase (cycle 1) and the inhibited phase (cycle 8), the toluene removal rate remained stable, averaging 19.1 ± 1.2 μmol L⁻¹ d⁻¹ (Figure 1C). Throughout the entire experimental period, O₂ was always below the detection limit, confirming strictly anaerobic conditions, and H₂ and CH4 were not detected or remained at negligible levels. Additionally, acetate concentrations remained consistently below 5 mg L⁻¹, indicating that it did not accumulate to a significant extent (data not shown).

From stoichiometric analysis of toluene consumption and sulfate reduction, we found that the observed cumulative sulfate consumption was approximately 1.9 ± 0.3 times higher than the theoretical value expected for complete oxidation of toluene to CO₂ coupled to sulfate reduction. This imbalance can be attributed to “futile” redox cycles that dissipate the excess of reducing potential by consuming electron acceptors without contributing directly to substrate oxidation (Sharma and Khandelwal, 2024). It was recently reported that futile sulfur metabolism can be employed as strategy to achieve rapid sulfide detoxification while maintaining intracellular redox homeostasis under substrate-excess conditions (Jia et al., 2026).

To further confirm that the removal of toluene was coupled with the reduction of sulfate, we conducted a targeted experiment where molibdate was employed as selective inhibitor of sulfate reduction in a microcosm inoculated using the original enriched culture. As a result, negligible degradation was observed in a microcosm containing molibdate, while clear toluene removal occurred in the molibdate-free control (Supplementary Figure 1).

Long-read 16S rRNA sequencing revealed a selected and low-complexity microbial consortium, co-dominated by Desulfoprunum (~34%) and Sulfurovum-affiliated sulfur oxidizers (Sulfurovaceae, ~34%), with minor contributions from Stenotrophomonas and Achromobacter spp. (< 10% each) and Pseudomonas, Rhizobium, Bacillus, Variovorax, and Brevundimonas (< 2%) (Figure 2). Such a co-dominance pattern is unusual, as sulfate-reducing enrichments are typically dominated by a single SRB lineage, such as Desulfobacula or Desulfosarcina in marine sediments, or Desulfobulbaceae in hydrocarbon-contaminated sites (Beller et al., 1996; Widdel, 2001; Kleinsteuber and Schleinitz, 2012; Rabus et al., 2016). Interestingly, the co-occurrence of SRB (Desulfoprunum), responsible for toluene degradation via dissimilatory sulfate reduction, and Sulfurovum-affiliated taxa, known to oxidize sulfide to sulfur, suggests a potential stabilizing mechanisms that may help mitigate the toxic effects of H₂S, a byproduct of sulfate reduction, thereby supporting long-term toluene degradation.

Genome-centric metagenomic analysis using long-read sequencing supported this hypothesis, enabling the reconstruction of seven medium to high-quality MAGs (Supplementary Table 2): BN30871 sp. TD_0, Achromobacter pulmonis TD_1, Desulfoprunum sp000769715 TD_2, UBA9337 sp. TD_4, Stenotrophomonas acidaminiphila_A TD_6, Stenotrophomonas maltophila_O TD_3, and Stutzerimonas stutzeri TD_5. Taxonomically, these MAGs are affiliated with the Patescibacteriota order Moranbacterales (UBA9337), Gammaproteobacteria (S. acidaminiphila_A, S. maltophila_O, A. pulmonis, S. stutzeri), Desulfocapsaceae (Desulfoprunum sp000769715 TD_2), and Sulfurovaceae (BN30871) (Figure 3).

MAG-resolved analyses revealed that Desulfoprunum sp000769715 TD_2 was the sole community member carrying the genes of the fumarate addition pathway (bssABC, bbsEFGH), which are required for the anaerobic activation of toluene and its conversion to benzoyl-CoA, together with the dissimilatory sulfate reduction genes (aprAB, dsrAB) (Winderl and Schaefer, 2007; Rabus et al., 2016). The genetic profile of Desulfoprunum sp000769715 TD_2 is consistent with that of a canonical SRB capable of mineralizing aromatic hydrocarbons under anaerobic conditions. Indeed, previous studies have reported the involvement of Desulfoprunum species in anaerobic degradation of aromatic compounds under sulfate-reducing conditions, including toluene, further supporting its role as the key toluene degrader in this microbial consortium (Junghare, 2015; Irianni-Renno et al., 2024; Hudari, 2025).

Interestingly, genes associated with other sulfur transformations, including assimilatory sulfate reduction (ASR), were more evenly distributed across the remaining MAGs (Figure 4). In particular, the Gammaproteobacteria MAGs (A. pulmonis TD_1, S. stutzeri TD_5, S. acidaminiphila_A TD_6) likely contributed to the sulfate demand, as they encoded nearly complete sets of cys genes required for assimilatory sulfate reduction to H₂S (Kushkevych et al., 2020). In particular, S. stutzeri TD_5 harbors the complete set of cysNCDHJI genes required for assimilatory sulfate reduction (ASR) through the cysteine biosynthesis pathway (Figure 3). ASR is essential for the production of sulfur-containing amino acids such as cysteine and methionine. A key step shared with dissimilatory sulfate reduction (DSR) is the ATP-dependent activation of sulfate to adenosine-5′-phosphosulfate (APS), followed by its conversion to 3′-phosphoadenosine-5′-phosphosulfate (PAPS). PAPS is then reduced to sulfite and subsequently to sulfide by sulfite reductase, after which sulfide is incorporated into L-cysteine via O-acetyl-serine(thiol)-lyase. The synthesis of L-cysteine from inorganic sulfate thus represents the main route for sulfur incorporation into biomass under these conditions (Kushkevych et al., 2020). Additionally, A. pulmonis TD_1 and S. stutzeri TD_5, together with the Sulfurovaceae MAG (BN30871 sp. TD_0), encode sulfide:quinone oxidoreductase (sqr), suggesting the potential for sulfide oxidation to elemental sulfur (S°), with electrons transferred to the quinone pool. While the presence of sqr indicates a pathway for sulfide detoxification, the ultimate electron acceptor in this anaerobic system remains unknown. Further studies are needed to clarify the electron transfer pathways involved in this process. This process links sulfide detoxification to energy conservation via reduction of the quinone pool (Landry and Ballou, 2021). Finally, the Patescibacteriota MAG UBA9337 sp. TD_4 (Moranbacterales) encodes a putative sulfhydrogenase (HydAB), which may catalyze the reduction of elemental sulfur (S°) back to hydrogen sulfide (H₂S) (Pedroni et al., 1995), thereby completing the cryptic sulfur cycle within the consortium.

Together, these complementary sulfur transformations outline a potentially distributed redox network that links toluene oxidation to internal sulfur turnover (Figure 4). In this network, DSR by Desulfoprunum couples' toluene oxidation with H₂S release, while co-occurring taxa oxidize H₂S to elemental sulfur (S°) via Sqr, thereby mitigating sulfide toxicity and regenerating electron acceptors. In parallel, other community members perform ASR to sustain biosynthetic demands, which likely accounts for the sulfate consumption (Figure 4).

Although similar sulfur cycling dynamics have been described in sedimentary systems, here we propose a conceptual framework in which tightly coupled sulfur redox interactions contribute to anaerobic hydrocarbon degradation in a microbial consortium. This framework is supported by genome-resolved evidence and community composition, rather than by direct measurements of intermediate sulfur species. In the enrichment studied here, toluene degradation does not appear to rely exclusively on canonical sulfate reduction but may be stabilized by a distributed sulfur-cycling network involving sulfate reducers, sulfide oxidizers, and sulfur assimilators (Figure 4).

Such coordinated interactions could enable internal sulfide detoxification and redox buffering, thereby sustaining long-term hydrocarbon degradation under sulfidic conditions.

From an ecological perspective, these findings reinforce the view of anaerobic hydrocarbon degradation as a community-driven process enabled by sulfur-cycling interactions, rather than the activity of a single taxonomic group. By integrating contaminant turnover with internal sulfur transformations, microbial consortia may enhance metabolic stability and resilience in anoxic environments.

Genome-resolved analyses reveal complementary metabolic roles across community members, providing a mechanistic basis for proposing cryptic sulfur cycling as a stabilizing feature of hydrocarbon-degrading communities. This model will benefit from future studies that directly quantify sulfur intermediates and functional activity (e.g., dissolved sulfide, elemental sulfur, proteomics), which will be important to confirm and refine our understanding of sulfur-coupled hydrocarbon degradation dynamics.

The enrichment culture described here demonstrates how anaerobic toluene degradation can be sustained through cooperative sulfur cycling involving multiple microbial guilds. Within the consortium, Desulfoprunum performs toluene oxidation via the fumarate addition pathway (Bss/Bbs) coupled to dissimilatory sulfate reduction (AprAB/AsrAB), releasing H₂S. This sulfide is subsequently detoxified and partially oxidized to elemental sulfur (S°) by co-occurring taxa, including Sulfurovaceae, Achromobacter, Stenotrophomonas, and Stutzerimonas carrying the sqr gene, while UBA9337 (Moranbacterales), encoding HydAB, may reduce S° back to H₂S (Figure 4), potentially closing a cryptic sulfur cycling. This distributed sulfur transformation network likely acts as a redox-buffering and detoxification mechanism, preventing sulfide accumulation. By partitioning key sulfur processes across specialized taxa, the community maintains redox balance and minimizes sulfide inhibition. Although our enrichment culture lacks the large-scale spatial stratification observed in natural environments (e.g., marine sediments or oxygen minimum zones), the co-occurrence of sulfate-reducing and putative sulfur-oxidizing taxa, together with their respective functional genes (e.g., dsrAB, aprAB, sqr, hydAB), suggests potential metabolic complementarity. This raises the hypothesis that microscale redox heterogeneity, possibly arising within microbial aggregates or biofilms, may enable the partial coexistence of these energetically opposing sulfur transformations even in a well-mixed bioreactor system (Jia et al., 2026). Nonetheless, we emphasize that this model is based on genome-resolved inference and future work will be needed to validate these interactions, including direct measurements of sulfur species and spatially resolved analyses to confirm the presence of redox gradients and the activity of the implicated pathways.

From an applied perspective, this cooperative sulfur metabolism offers a mechanistic explanation for the functional resilience and robustness of anaerobic hydrocarbon-degrading consortium. Such internal sulfur cycling could enhance the performance of in situ bioremediation systems, including engineered bioreactors, by mitigating sulfide toxicity and sustaining electron flow, provided that sulfate remains available.

Ecologically, these findings extend the relevance of cryptic sulfur cycling, previously described in marine sediments, to hydrocarbon-fed consortia in contaminated aquifers. In this context, hydrocarbons appear to fuel a distributed sulfur redox network involving reducers, oxidizers, and recyclers. Such multi-guild cooperation may stabilize degradation processes, enhances community resilience under stress conditions (e.g., contamination), and embed contaminant turnover within broader carbon–sulfur biogeochemical coupling. If widespread, hydrocarbon-driven cryptic sulfur cycling could represent a general ecological mechanism that sustains microbial functioning and contaminant attenuation in anthropogenically impacted anoxic environments.