Competition-cooperation in the chemoautotrophic ecosystem of Movile Cave: first metagenomic approach on sediments

Background Movile Cave (SE Romania) is a chemoautotrophically-based ecosystem fed by hydrogen sulfide-rich groundwater serving as a primary energy source analogous to the deep-sea hydrothermal ecosystems. Our current understanding of Movile Cave microbiology has been confined to the sulfidic water and its proximity, as most studies focused on the water-floating microbial mat and planktonic accumulations likely acting as the primary production powerhouse of this unique subterranean ecosystem. By employing comprehensive genomic-resolved metagenomics, we questioned the spatial variation, chemoautotrophic abilities, ecological interactions and trophic roles of Movile Cave’s microbiome thriving beyond the sulfidic-rich water. Results A customized bioinformatics pipeline led to the recovery of 106 high-quality metagenome-assembled genomes from 7 cave sediment metagenomes. Assemblies’ taxonomy spanned 19 bacterial and three archaeal phyla with Acidobacteriota, Chloroflexota, Proteobacteria, Planctomycetota, Ca. Patescibacteria, Thermoproteota, Methylomirabilota, and Ca. Zixibacteria as prevalent phyla. Functional gene analyses predicted the presence of CO2 fixation, methanotrophy, sulfur and ammonia oxidation in the explored sediments. Species Metabolic Coupling Analysis of metagenome-scale metabolic models revealed the highest competition-cooperation interactions in the sediments collected away from the water. Simulated metabolic interactions indicated autotrophs and methanotrophs as major donors of metabolites in the sediment communities. Cross-feeding dependencies were assumed only towards 'currency' molecules and inorganic compounds (O2, PO43−, H+, Fe2+, Cu2+) in the water proximity sediment, whereas hydrogen sulfide and methanol were assumedly traded exclusively among distant gallery communities. Conclusions These findings suggest that the primary production potential of Movile Cave expands way beyond its hydrothermal waters, enhancing our understanding of the functioning and ecological interactions within chemolithoautotrophically-based subterranean ecosystems. Supplementary Information The online version contains supplementary material available at 10.1186/s40793-022-00438-w.

The Movile Cave (SE Romania) was the first to be mentioned as defying the conventional view of subterranean ecosystems as supported by aboveground photosynthesis. Movile Cave is a small-surfaced, closed chemoautotrophic system [8] driven by in-house sulfur and methane oxidation and CO 2 fixation as primary production processes. The cave, located in south-eastern Romania (Dobrogea region), developed in oolitic and fossil-rich limestone of the Sarmatian age (Late Miocene) and was sealed off during the Quaternary by a thick and impermeable layer of clays and loess [17]. Movile Cave has a complex geological evolution with an ongoing speleogenesis driven by two main processes: the sulfuric acid corrosion in the partially submerged lower cave level; and the condensation-corrosion processes active in the upper level of the cave [18]. The upper gallery (approx. 200 m long) is dry, whereas the lower gallery (approx. 40 m long) is partially flooded by sulfidic hydrothermal waters (T ~ 20.9 °C). The two cave levels converge in the so-called Lake Room (Fig. 1). Air pockets (Air Bells) are present in the submerged gallery (Fig. 1).
Here, an active redox interface is created on both the water's surface and cave walls, colonized by floating microbial mats and biofilms. Most of our knowledge of this ecosystem comes from studies of trophic chain components such as microbial mats associated with the surface of the hydrothermal water and its close proximity [19][20][21][22]. In addition, the rich endemic invertebrates of Movile Cave live in the proximity of the sulfidic water. The fauna is abundant in the Lake Room [15,19], particularly at the water surface in the Air-Bells [23], where the redox potential is high.
Except for a snapshot investigation of the rock's surface collected at about 2 m away from the water [24]-which revealed the presence of sulfur oxidizers and methylotrophic bacteria-Movile Cave has not yet been explored for the microbial life associated with sediments distant from the sulfidic waters. Also, no metagenome-assembled genomes (MAGs) reconstruction studies have been carried out, and the genomic information available in databases about this unique environment is limited to a Ca. Methylomonas sp. LWB [19] and a Thiovulum sp. [22] genome assemblies from the microbial mat.
Considering the unique features of Movile Cave as a model chemoautotrophic subterranean ecosystem and its overlooked microbial communities away from the sulfidic water, we hypothesized that the scarcity of high-energy organic substrates beyond the primary energy source of sulfidic water would strengthen the ecological interactions among community members to optimize the flow of nutrients and ultimately, the energy gain. To address this hypothesis, genome-resolved metagenomics was employed to infer the diversity and abundance of spatially distinct sediment-associated microbiomes and to predict their metabolic abilities, ecological interactions and roles within the cave's ecosystem. Our genome resolved metagenomic study adds 106 higher-quality MAGs to the Movile Cave microbial genomes. The results broaden our understanding of sediment microbial communities' role within the complex food web of Movile Cave, which supports 52 endemic invertebrate species [15].

Study site and sample collection
Sediment samples were collected from seven sites located in the cave's upper and lower levels (galleries) (Fig. 1) in December 2019. Sampling sites were selected based on their physical characteristics (detailed in Table 1 and Additional file 1: Fig. S1) and positions related to the Lake Room. Three sampling sites were near the water, in the Lake Room (PMV1, PMV3, PMV4), and the other four (PMV2, PMV6, PMV7, PMV8) in the upper, dry galleries. The PMV5 sample consisted of sulfidic water and was not the subject of this study. Stable environmental conditions typical to a subterranean habitat combined with hydrothermal activity drive a constant temperature of the water and air at ≈ 20.9 °C within Movile Cave. Aseptically collected top sediments (up to 50 g) were kept at 4 °C after sampling and while transporting, then at − 20 °C until processing.

Mineralogy and geochemistry measurements
Powdered X-ray diffraction analyses were performed on sediments in order to establish their mineralogy. Samples were analyzed with a Rigaku Ultima IV diffractometer in parallel beam geometry equipped with CuKα radiation (wavelength 1.5406 Å). The XRD patterns were collected in 2Θ range between 5 to 80 with a speed of 2°/min and a step size of 0.02°. PDXL software from Rigaku, connected to the ICDD database was used for phase identification. The quantitative determination was made using the RIR (Reference Intensity Ratio) methodology.
The pH and electric conductivity were measured in 1:5 solid to water extract using a Seven Excellence (Mettler Toledo) multimeter. The N, C and H were measured on freeze-dried samples using a Flash 2000 (Thermo Scientific, Thermo Fisher Scientific, Waltham, MA, United States) analyzer. For measurement of metals, samples were digested using a 1:3 mixture of 65% HNO 3 and 37% HCl at reflux conditions. Major elements (Na, K, Ca, Mg, Al, Fe, S, P) in the digested samples were determined by inductively coupled plasma optical emission spectrometry (ICP-OES) using a 5300 Optima DV spectrometer (Perkin Elmer, Shelton, CT, United States), while trace metals and metalloids by inductively coupled plasma mass spectrometry (ICP-MS) using an ELAN DRC II spectrometer (Perkin Elmer, Waltham, MA, United States). Additional measurements were done for the aqueous fraction of the PMV4 sample. Total nitrogen (TN) was measured in unfiltered water, while dissolved nitrogen (DN) was measured in water samples filtered through 0.45 µm pore size PTFE syringe filters by catalytic combustion and chemiluminescence detection using a Multi N/C 2100 (Analytik Jena, Jena, Germany) analyzer. Total carbon (TC) and total inorganic carbon (TIC) were measured from unfiltered samples, while dissolved ) were measured by ion-chromatography, using a 761 IC compact ion chromatograph (Metrohm, Herisau, Switzerland), equipped with a Metrosep 5-100/4 column and a Metrosep A Supp 4/5 mm guard column after filtration of samples through 0.45 μm cellulose acetate membrane filters. Na, Mg, K and Ca were measured in the filtered samples after acidulation with 65% HNO 3 to a pH < 2 by ICP-OES.

DNA extraction and sequencing
Community DNA from sediments was extracted using the DNeasyPowerMax Soil Kit (Qiagen, Gaithersburg, MD, United States); precipitated with 1:10 (v:v) of 3 M Na-Acetate (pH 5.2) and 1:1 (v:v) of 100% ethanol; air-dried and finally resuspended in 100 µL of sterile TE buffer. DNA extracts were quantified with Qubit TM dsDNA BR Assay Kit and a Qubit 4.0 fluorometer (Invitrogen, Thermo Fisher Scientific, Waltham, MA, United States). Novogene Europe Co., Ltd (Cambridge, UK) performed the library preparation and shotgun metagenomics sequencing. For library preparation, Illumina protocol NEBNext ® Ultra ™ DNA Library Prep Kit (New England BioLabs, Ipswich, MA, United States) was used starting from 200 ng DNA. Sonication was used to randomly fragment the DNA to a size of around 350 bp. Sequencing indexes were added to each sample. The PCR amplification was conducted after the DNA fragments were end-polished, A-tailed, and ligated with the adapter. AMPure XP system was used for PCR product purification. Agilent 2100 Bioanalyzer was used to analyze the size distribution of the DNA library, and qPCR was performed to quantify DNA according to conventional Illumina library quantification protocol. The samples were sequenced and paired-end reads were generated (PE150) on the Illumina NovaSeq 6000 platform.

Reads assembly and annotation
Metagenomic raw reads (fastq) were adapter-trimmed using Trimmomatic v0.39 [26] and, their quality was evaluated using FastQC v0.11.9 [27]. The resulting sequences were assembled into contigs using MEGA-HIT v1.2.9 [28], with the following parameters: no-mercy, k-min = 31, k-max = 101, k-step = 10. The assembly of reads was done separately, for each metagenome (sample). Assembled contigs (> 1 kb) were binned using two binners: MaxBin v2.0 [29] and MetaBAT2 v2.12.1 [30]. DAS Tool (v1.1.2) was applied to integrate the bins generated from the two methods [31]. Only bins with over 70% completion and less than 10% contamination based on the results of CheckM v1.1.3 [32] were analyzed. Protein-coding sequences (CDS) of selected metagenome-assembled genomes (MAGs) were predicted and annotated with Prodigal v2.6.1 [33] and DIAMOND v2.0.5 [34] against Uniprot databases. Furthermore, the functional annotation of CDS and the prediction of metabolic pathways related to sulfur, nitrogen, carbon, and methane metabolism was carried out using the KEGG's Ghost Koala tool [35]. For confirmation, CDS of interest were identified by blastp against NCBI nr database. In addition to the results obtained by functional annotation, the shotgun preprocessed reads were mapped to a set of reference sequences of nitrogen and carbon fixation hallmark genes with Minimap2 v2.21 with the short reads settings (-sr option).

Taxonomic assignment and abundance of MAGs
Taxonomic classification of MAGs was assessed using Genome Taxonomy Database Toolkit (GTDB-Tk v1.5.0, reference data for GTDB R06-RS202) [36] based on 120 bacterial and 122 archaeal marker genes. For phylogenetic analysis, closely related MAGs, reference genomes, or type strains genomes were selected from GTDB (July 2021). The marker genes alignment was concatenated and then imported into R using the readAA Multiple Alignment function from the package Biostrings. Phylogenetic distances were calculated using the function dist. alignment from the package seqinr, and a phylogenetic tree was built using the neighbor-joining tree estimation as implemented in the function nj in the R package ape.
The tree was visualized using the online iTOL tool [37]. Finally, for MAGs coverage, the filtered shotgun reads were mapped back to contigs belonging to each MAG with Minimap2 v2.21 [38], and the coverage was summarized with the coverage command from SAMtools v1.11 [39]. The relative abundance of each MAG in each metagenome was given by the number of reads mapped per Kb of MAG divided by Gb of corresponding metagenome (RPKG) [40].

Genome-scale metabolic reconstruction and interaction models
Genome-scale metabolic models derived from MAGs (metaGEMs) were reconstructed using CarveMe v1.5.1 [41], with the default CPLEX solver (v12.8.0.0) (IBM, 2019). Two sets of metaGEMs were reconstructed: one set without gap-filling for any particular medium, based only on genetic evidence; and another set gap-filled for a custom minimal medium that would guarantee the models' growth in a nutrient-poor environment as in Movile Cave. MEMOTE tests suite [42] was applied to the gap-fill models initialized for the minimal medium to ensure that the metaGEMs could generate biomass and reproduce growth in a scarce environment. The custom minimal medium composition is detailed in Additional file 2: Table S1. Community simulations were performed with SMETANA v1.1.0 [43] using the reconstructed metaGEMs. For the simulations, the samples were divided into two groups based on their distance from the sulfidic water: samples from the lower gallery in the Lake Room (PMV1, PMV3, PMV4) and samples from the upper gallery (PMV2, PMV6, PMV7, PMV8). Metabolic interactions were assessed for each sample, and metabolic dependencies were calculated for each condition (lower vs. upper gallery). Global metabolic interactions algorithms MRO (metabolic resource overlap) and MIR (metabolic interaction potential) were applied for both sets of reconstructed metaGEMs assuming complete or minimal media. MRO/MIR scores were used as a measure of how much the species compete for the same metabolites (MRO) and how many metabolites they can share to decrease their dependency on external resources (MIP). Only the gap-filled models were used for metabolic dependency, calculated assuming a complete medium. The complete medium composition relies on the BIGG Database [44] (http:// bigg. ucsd. edu/). By doing so, we imposed the growth of community members on the minimal medium as a constraint. However, we relaxed the community simulation to all possible metabolic exchanges of a complete environment (to identify the metabolic exchanges essential for the community's survival). The obtained SMETANA scores estimate the strength of metabolic coupling (metabolites exchange) in a community. Only a SMETANA score of 1 was considered for analysis. The interactions characterized by a score of 1 are viewed as essential metabolic exchange (essentialities) for the community, and this metabolic coupling is considered a "dependency". SMETANA score was used to highlight the cross-feeding interactions and the keystone MAGs (phyla). Significantly different (p < 0.001) compound exchanges across conditions were identified with Wilcoxon rank-sum test using the compare_means function of ggpubrR package. Model's reconstruction and community simulations parameters and the simulation codes are summarized and detailed in Additional file 2: Tables S2, S3.

Statistical analysis
We used multidimensional scaling (MDS) to represent the sediment samples in a bi-dimensional space. For the MDS we used a proximity matrix of similarity between the objects (the chemical characteristics of the sediment samples) to the coordinates of these same objects in a 2-dimensional space so that the objects may be visualized easily. Stress is the goodness-of-fit statistic that MDS tries to minimize. Stress varies between 0 and 1, with values near 0 indicating better fit. The analysis was done in XLSTAT 2021.4.1 (Addinsoft, New York).

Physicochemical and mineralogical characterization
The sediments had a slightly alkaline pH, in the range of 7.9-9.2, and low (PMV1, PMV7 and PMV8) or high (PMV2, PMV3, PMV4 and PMV6) electrical conductivity ranging from 60.6 to 1738.8 µS cm −1 ( Table 2). The physicochemical characteristics of the sediment samples are shown in Table 2. Total carbon (TOC) and total nitrogen (TN) content, measured for the aqueous fraction of the sample collected from the water edge (PMV4) were low (TOC = 14.9 mg L −1 , TN = 3.7 mg L −1 ) with a TOC/TN of ~ 4 mg L −1 (Additional file 3: Table S4).
In the upper cave gallery, the identified minerals are calcite, dolomite, aragonite, quartz, muscovite, gypsum, and in small amounts, goethite. Near the Lake, in the lower gallery, the association is dominated by dolomite and calcite, followed by quartz, gypsum and muscovite ( Table 2). Carbonates are represented mainly by calcite derived from the limestone bedrock, dolomite and aragonite. Dolomite is present in substantial amounts in both levels of the cave. The bedrock mineralogy reveals the presence of magnesium-rich carbonate, which could be the explanation for the occurrence of aragonite, a signature of PMV8 mineralogy. Samples mineralogy is shown in Table 2 and Additional file 3: Table S5.
A multidimensional scaling (MDS; Fig. 2) analysis shows that PMV6 and PMV8 are more dissimilar from the other samples, PMV6 having the higher electrical conductivity and higher S, Na, Mg, Ca combined with low K, Al and Fe content. The presence of gypsum (sulfate mineral) in that sample (Table 2) explains the high content of S and Ca. In contrast, PMV8 had the lowest electrical conductivity and content in S, Ca (no gypsum in this sample), with a higher concentration of Al, Li, and Na. PMV3 and PMV4 were also separated, while PMV1, PMV2 and PMV7 were more similar along the MDS 1st Dimension. PMV3 had low Na and Mg, and high Fe (the only sample where goethite-iron(III) oxide-hydroxide-is well represented), V, As, Cr, etc. PMV4 had the highest content of Ca, Ni, and Ba, which probably gives the black color to the sample. The other samples (PMV1, PMV2, and PMV7) were less dissimilar, with muscovite (hydrated phyllosilicate mineral of aluminum and potassium) and quartz (SiO 2 ) in their mineralogical composition.

MAGs recovery and phylogeny
A total of 106 medium-to high-quality MAGs (> 70% complete and < 10% contamination) [48] were recovered from the sediment metagenomes of Movile Cave. The taxonomic assignment performed using GTDB-tk indicated that the recovered MAGs span 19 bacterial and three archaeal phyla. Out of those, 4 were candidate phyla Ca. Patescibacteria, KSB1, Krumholzibacteriota and Zixibacteria. The percentages of reads mapped to MAGs varied among samples from ca. 80% to 13% for MAG recovery efficiency (summarized in Additional file 4: Table S6). The low binning efficiency in some of the samples points toward a very diverse community (feature specific to sediment samples) where the low abundance species were unable to be assembled into MAGs as a consequence of low or partial coverage. Based on the relative abundance of MAGs in the sample (expressed as RPKG), the most abundant phyla were Acidobacteriota, Chloroflexota, Proteobacteria and Planctomycetota (> 20 RPKGs in any of the samples) (Fig. 3A.). In the lower gallery, Chloroflexota was abundant in PMV1 and PMV3, while Proteobacteria (class Gammaproteobacteria) dominated PMV4. In the upper, dry gallery, Acidobacteriota was the most abundant phylum for PMV8 and PMV7. For the phyla with a lower relative abundance (< 20 RPKGs), Ca. Patescibacteria was the signature phylum for PMV2, Thermoproteota for PMV3, Methylomirabilota for PMV7 and Ca. Zixibacteria for PMV8 (Fig. 3B.). Out of the 106 MAGs, 25 were classified as high-quality (> 90% complete and < 5% contamination) (Fig. 3C.), but only two were assigned to species level (PMV3.39 and PMV4.117 affiliated to Chromohalobacter marismortui) with available genomes in the NCBI database (gANI > 95%). One-third of the MAGs (32 of 106) was classified by Relative Evolutionary Divergence (RED) as new species of established GTDB taxa [49]. Given the taxonomic novelty of recovered MAGs, RefSeq type material genomes (GCF) and representatives MAGs from GeneBank (GCA) were incorporated into the phylogenetic tree construction for a better representation (Fig. 3C.) MAGs attributes, extended taxonomy and abundance are summarized in Additional file 4: Tables S7-S9.

Potential for biogeochemical cycling of S and N, CH 4 oxidation and CO 2 fixation Sulfur metabolism
Microbial sulfur cycling has been proposed as a driving force for bacterial proliferation in microbial mats of Movile Cave [24] but not investigated in the cave's sediments. Here, we examined the presence of the marker genes for sulfur cycle in the metagenomes of Movile Cave's sediments. For sulfur oxidation (Fig. 4), the complete canon- ) reduction (sat + aprAB) was present in the PMV4 dataset, in a MAG affiliated to order Thermodesulfovibrionales (phylum Nitrospirota) (PMV4_ maxbin.018). An NCBI blast (blastx against RefSeq Select database) showed the AprAB sequences similarities to  /S 0 ) disproportionation and tetrathionate respiration were partially encoded (phsA, phsC; ttrA/B) in different phylogenetic groups, but a Planctomycetota-affiliated MAG (PMV1.61) that might encode both as ttrB and phsA were present in the bin. The production of hydrogen sulfide from polysulfides by sulfhydrogenase complex (hyd(G,B),(A,D)) was partially encoded, especially by MAGs affiliated to Acidobacteriota, Chloroflexota and Planctomycetota.

Nitrogen metabolism
Analogous to sulfur, we addressed the potential for biogeochemical N cycling across sediments of Movile Cave. Primary producers capable of obtaining energy for autotrophy by nitrogen oxidation (nitrification) were questioned. The first step of nitrification, the ammonia NH 3 + to nitrite NO 2 − oxidation, was suggested in Movile Cave sediments as ammonia monooxygenase AMO (amoABC) was encoded in PMV1 to PMV3 datasets. All three subunits of AMO's were encoded in an order Methylococcales-affiliated MAG (PMV2.70_sub). However, the amoABC sequences were homologous (> 65% similarity) to the sequence encoding particulate methane monooxygenase pMMO in the genus Methyloterricola. The methanotrophs can oxidize both substrates (NH 3 and CH 4 ) but grow only on their characteristic substrate [50]. Interestingly, AmoA-like subunit was annotated in an archaeon MAG affiliated to Nitrososphaera genus (phylum Thermoproteota) (PMV3_maxbin.65_sub). For the second step of nitrification, the oxidation of NO 2 − to NO 3 − , none of the bins encoded the nitrite oxidoreductase (nxrAB) found in known nitrite-oxidizing genera. Also, no marker genes (hzsA and hzo) for the anaerobic ammonium oxidation (anammox) (NH 3 to N 2 H 4 and then to N 2 ) were annotated in the Movile datasets.
In the chemoautotrophic cave environment, the nitrogen demand can be met by converting inorganic nitrogen to a biologically useful form by nitrogen reduction. The conversion involves microbial dinitrogen fixation or nitrate assimilation. Surprisingly, the N 2 fixation pathway was not detected in any sediment datasets, as no nitrogenase genes (nif) for the reduction of atmospheric molecular nitrogen (N 2 ) to ammonia (NH 3 + ) was annotated in MAGs. This result was reinforced by mapping the metagenomics reads against a set of nif reference sequences (listed in Additional file 5: Table S10) as no nifH/D were found in the PMV datasets except PMV4 (Additional file 5: Fig. S2). The assimilatory process (NO 3 − , NO 2 − conversion to NH 3 + ), catalyzed by nitrate and nitrite reductases (nasA/B, narB; nirA), was fully encoded in PMV3, PMV6, and PMV8 datasets. The bins that encoded the nitrate-assimilatory enzymes were taxonomically diverse as the ability is however widely distributed among bacteria.
Methanogenesis potential was absent in the Movile sediment datasets as suggested by the absence of methyl coenzyme-M reductase (mcrA) marker gene across investigated MAGs.

CO 2 fixation
The CO 2 fixation is presumably critical in a chemoautotrophic ecosystem such as Movile Cave. To verify the autotrophic potential of the sediment communities, the presence of genes for the key enzymes of CO 2 fixation pathways, namely the Calvin-Benson-Bassham (CBB) cycle (key enzyme: ribulose 1,5-bisphosphate carboxylase/oxygenase (RubisCO), genes cbbL and cbbS), the reverse tricarboxylic acid (rTCA) cycle (key enzyme: ATP citrate lyase, genes aclA and aclB) and the reductive acetyl-CoA, or Wood-Ljungdahl (WL) pathway (key enzyme: CO dehydrogenase/acetyl-CoA synthase (CODH/ACS complex), genes acsA and acsB) was investigated.
In addition to MAG analysis, mapping the metagenomics reads against a set of reference CO 2 fixation gene sequences (listed in Additional file 5: Table S10) highlighted the CBB cycle key genes as present in all datasets, including PMV6 and PMV7, and the WL pathway marker genes as hallmark genes for only PMV4 dataset (Additional file 5: Fig. S2).
Interestingly, none of the genes of interest for biogeochemical cycling of S and N, CH 4 oxidation and CO 2 fixation were annotated in MAGs assigned to Ca. Patescibacteria or Myxococcota phyla.
All the functional genes implicated in the sulfur and nitrogen cycling, methane oxidation and CO 2 fixation annotated in MAGs from Movile Cave sediment are listed in Additional file 6: Tables S11-S14 and an overview of the genes, encoding MAGs and pathways are shown in Fig. 4.

Potential metabolic interactions and dependencies within the communities
To gain a deeper insight into the microbial metabolic webs in Movile Cave sediments, we used the MAGs to construct metagenome-scale metabolic models (metaGEMs) for simulation and prediction of potential metabolic interactions and dependencies within the communities. Noteworthy, no metabolic models could be constructed for three MAGs belonging to Ca. Patescibacteria (PMV2.56, 2.61 and 2.72) thus were excluded from the analysis. The inability to construct metabolic models for MAGs belonging to Ca. Patescibacteria can be a consequence of the size of their reduced genome (usually < 1 Mb). The rest of the metaGEMs could generate biomass and reproduce growth in the simulation conditions, as evidenced by MEMOTE tests results included in the Additional file 7: Table S15. The competitioncooperation potential predicted by SMETANA (Fig. 5A.; Additional file 7: Table S16) was highest in the communities associated with the upper gallery (e.g. PMV7). This was evident, especially in the case of the necessary resources overlap (competition). The pattern persisted regardless of metaGEMs reconstruction or community simulation parameters. When community simulation assumed minimum nutrient availability, as presumed for the Movile Cave environment, the community members of the upper gallery (samples PMV6 to PMV8) exhibited the highest similarities of metabolic requirements (competition) and highest potential for community self-sufficiency (cooperation) (Fig. 5A.b.).
For an insightful analysis of the microbial metabolic webs in the Movile Cave sediments, we grouped the samples based on distance, namely samples from the lower, Fig. 5 Competition-cooperation landscape of each sample and cross-feeding interactions across wet and dry galleries. A. The competition (MRO) and cooperation (MIP) scores (divided by the numbers of MAGs in the community) are shown for different reconstruction and simulation parameters: a. metaGEMs reconstructed only on genetic evidence and community simulation on complete medium (unconstrained environment); b. metaGEMs reconstructed by gap-filling for minimal media and community simulation on minimal media (constrained environment). B. and C. Alluvial diagrams showing compounds exchanged (SMETANA score = 1) in each condition (lower/wet and upper/dry gallery) between the donor (left) and receiver (right) phyla. The colors are used only to distinguish distinct components of the alluvial diagrams 'wet' gallery (PMV1, PMV3, PMV4) and samples from the upper, 'dry' gallery (PMV2, PMV6, PMV7, PMV8). We examined the cross-feeding interactions and keystone phyla across those two conditions ('wet' and 'dry'). Our simulation results (SMETANA score = 1) suggested that only 'currency' molecules and inorganic ions (O 2 , PO 4 3− , H + , Fe 2+ , Cu 2+ ) were exchanged in the lower gallery. Here, the donor MAGs belonged to archaea Ca. Thermoplasmatota (PMV4.75), bacteria candidate division NC10 (Methylomirabilota) (PMV4.88) and the bacteria genus Methylocella (PMV1.12) (Fig. 5B.). In the upper gallery, besides inorganic ions (Fe 2+ , NO 3− ), hydrogen sulfide (H 2 S) and methanol (CH 4 O) were also exchanged. Here, the essential donors were affiliated to the Alphaproteobacteria within order Dongiales (PMV8.34_sub) and Planctomycetota within order Phycisphaerales (PMV2.11) (Fig. 5C.). For both conditions, the taxonomic affiliation of receiver MAGs was diverse, represented by MAGs belonging to 8 and 17 phyla. The MAGs that interact in each simulation condition and their role as donors or recipients of exchanged metabolites are listed in Additional file 7: Tables S17 and S18.
Out of the predicted essential compounds (SMETANA score = 1) that are readily exchanged among MAGs, H 2 S, NO 3 − , methanol, and H + exchanges were significantly different (Wilcoxon rank-sum test, BH adjusted p-values < 0.001) across upper and lower galleries (Additional file 7: Table S19).

Discussions
Our current understanding of microbial life within the Movile Cave ecosystem was limited to the hydrothermal waters, as only microbial mats, water samples and lake sediments were previously investigated [18,19]. Our study gives a first insight into the cave's chemoautotrophic and primary production potential beyond the hydrothermal waters, namely the cave's sediments.
Based on physicochemical and mineralogical characterization, sampled sediments showed a very different chemical composition, even if located just meters away from one another. PMV6 and PMV8 differed the most from the other samples, followed by PMV3 and PMV4. Mineralogy partially supports such chemical discrepancies, but some of the best-represented elements are of different, unknown origins, probably linked to the history of the region where the cave is located and the input of the underground sediments from the Miocene to the more recent climatic events.
The taxonomic novelty of recovered MAGs was high, with over 60% unaffiliated to a genus and 30% classified by RED as new species of established GTDB taxa. The community was diverse with 22 microbial phyla, mainly Acidobacteriota, Chloroflexota, Proteobacteria, Planctomycetota, Ca. Patescibacteria, Thermoproteota, Methylomirabilota, and Ca. Zixibacteria based on the relative abundance of MAGs in samples. Despite the high novelty at lower taxonomic levels, the overall phyla diversity of sediment MAGs was typical of sulfidic and nonsulfidic cave microbiology [51][52][53]. Something worth noticing is the lack of Betaproteobacteria among recovered MAGs, since it was previously demonstrated for the Movile Cave ecosystem [16,21].
Functional analyses suggested the presence of chemolithoautotrophic primary producers in the cave sediments from both galleries. As expected, sulfur oxidation, proposed as the ecosystem's driving force, was mainly present in the lakeside dataset PMV4. The complete oxidation pathways (Sox and Dsr pathways), RuBisCO and partial nitrate/nitrite respiration and denitrification were encoded by the Thiohalomonadales-affiliated MAG (class Gammaproteobacteria). As Thiohalomonas members are obligate chemolithoautotrophic facultative anaerobic sulfur-oxidizing bacteria [54], this MAG might play a key role as a primary producer within the ecosystem. It uses reduced sulfur compounds as energy sources, nitrate as an electron acceptor, and assimilating CO 2 via the Calvin-Benson-Bassham cycle. Another primary producer inferred from the lake edge sample was an Arcobacteraceae-related MAG (phylum Campylobactero ta/Epsilonproteobacteria) that encoded the capacity for sulfur oxidation (Sox) potentially coupled with nitratereduction (NapB) and CO 2 fixation via the reverse TCA cycle (aclAB). Those are typical metabolic features of Epsilonproteobacteria chemolithotrophic primary producers from the deep-sea hydrothermal vent [55]. Similar metabolic traits were found for the recently characterized Thiovulum sp. that dominates Movile Cave's hypoxic Air Bells microbial mat [22]. Moreover, Arcobacteraceae family was found dominant in the white filaments from the thermal sulfidic spring of Fetida Cave, Italy [56]. Sulfur oxidation was not, however, limited to the lakeside sediments. Genes coding for Sox pathway (soxCD) were found in PMV2, 6 and 7 datasets of the upper, dry gallery across MAGs affiliated to Alpha-, Gammaproteobacteria and Gemmatimonadota phylum. The Gemmatimonadota role in the sulfur cycle is uncertain, but MAGs with similar attributes were identified in the Siberian soda lakes [6]. In the sediments from around the lake, sulfur respiration was also evident, which is less critical in self-sustaining ecosystems and a characteristic of heterotrophs. In addition to previous reports [16,21] describing Deltaproteobacteria sulfate reducers in Movile Cave water and microbial mat samples, we assembled a hypothetical thermophilic sulfate reducer MAG affiliated to Nitrospirota (order Thermodesulfovibrionales) possessing DSR pathway genes (aprAB). These results might indicate that sulfate reduction is important in and around the sulfidic water but not away from it.
Considering the relatively high (0.2-0.3 mM) ammonium concentrations in the cave waters [57], it was formerly implied that nitrogen oxidation might also be a driving force for the ecosystem [21]. A few Nitrospirotaaffiliated MAGs were identified in the upper and lower galleries, but none of the nitrite or ammonia-oxidation genes were annotated. No nitrifying bacteria were identified based on gene annotation, but an ammonia-oxidizing Nitrososphaera archaeon encoding for an AmoA-like subunit was assembled from the Lake Room PMV3 sediments. Nitrososphaera archaea are facultative chemolithoautotrophs with potential metabolic flexibility [58,59]. Our metabolic annotation results suggest that the nitrogen cycle throughout the Movile Cave sediments is mostly driven by nitrate/nitrite respiration and denitrification. Those nitrogen reduction processes are linked to taxonomically diverse bacteria and spread among all investigated locations. Assimilatory and dissimilatory nitrate reduction to ammonium were predicted in all datasets, emphasizing the ability of the microbial communities to provide bioavailable N for the other trophic links in the Movile Cave ecosystem. Methanotrophs were also postulated as primary producers in the water and microbial mats [19][20][21]24, 60] as 1-2% methane concentration was highlighted, especially in the Air-Bells [57]. In this study we were able to assemble methanotrophic MAGs from both lower and upper galleries. These were affiliated to the uncultured UBA1147 genus of Methylococcales (in PMV1, PMV2, PMV6) and Methylocella genus of Rhizobiales (in PMV1, PMV3, PMV6) and encoded for subunits of pMMO monooxygenase. Although no methane-oxidation gene could be annotated, Methylomirabilales-affiliated (candidate division NC10) MAGs were found in PMV2 (order Rokubacteriales), PMV4 and PMV7 (order Methylomirabilales) datasets. Interestingly Methylomirabilales-affiliated MAG was one of the most abundant MAGs in PMV7. The presence of Methylococcales (Methylococcus, Methylomonas, Methylocaldum) and Rhizobiales (Methylocystis/Methylosinus, Methylocella) methanotrophs was previously postulated in the microbial matsbased on pMMOA-sequence clones from a CH 4 -enriched culture [19,20]. The genus Methylocella comprises facultative methanotrophs that utilize multicarbon compounds (acetate, pyruvate, succinate, malate, and ethanol) [61]. While the order Rokubacteriales might comprise non-methanotrophic bacteria [62], Methylomirabilales are known autotrophic (CO 2 fixing via CBB cycle), denitrifying methanotrophs that can oxidize methane anaerobically. The molecular oxygen needed for methane oxidation is generated by reducing nitrate to dinitrogen gas and O 2 and bypassing the nitrous oxide formation [63][64][65]. This is considered the fourth biological pathway known to produce oxygen besides photosynthesis, chlorate respiration, and the detoxification of reactive oxygen species [63].
The key enzymes for CO 2 fixation were found in 5 out of 7 sediment samples (except PMV6, PMV7), including PMV8 collected farthest from the lake. Noteworthy, in PMV8 MAGs both the CBB and rTCA enzymes were annotated. WL pathway of anaerobic CO 2 fixation was postulated only in the water edge dataset PMV4. A sulfur-respiring (aprAB) MAG affiliated to the order Thermodesulfovibrionales (phylum Nitrospirota) encoded for the WL pathway key enzymes (acsAB). However, the Thermodesulfovibrionales order includes members (genus Thermodesulfovibrio) that are chemoorganotroph, fermentative and dissimilatory sulfate-reducing bacteria [66]. Therefore, the MAG encodes the WL pathway in reverse to break down acetate to CO 2 and H 2 . Uncultured Nitrospirota MAG from PMV8 dataset encoded key genes for rTCA cycle, confirming early evidence revealed by stable-isotope probing (SIP) of water and microbial mat samples that Nitrospirota might be able of CO 2 fixation [21]. CBB cycle key enzymes were detected in most sediment samples, except for PMV6 and PMV7. The archaeal MAG encoding CBB capability found in the water's edge sample (PMV4) pertains to Ca. Thermoplasmatota's class EX4484-6 that includes MAGs retrieved from marine hydrothermal vent sediments (BioProject PRJNA362212). Thermoplasmata members were also identified in the snottites (thick snot-like biofilm) in the Frasassi caves and in the gypsum moonmilk of Fetida Cave, both in Italy [51,67,68]. The presence of RuBisCO (cbbL) and sulfhydrogenase (hydBG) subunits in the genome may point toward an autotrophic organism that couples sulfur (S 0 , S n ) reduction to hydrogen oxidation. Among the bacteria encoding RuBisCO, a MAG affiliated to the order Dongiales from the upper gallery PMV8 stands out. GTDB's Dongiales order comprises the members of the formerly known family Rhodospirillaceae. The RuBisCO enzyme presence and the lack of genes for the photosynthetic reaction center (pufLM) are intriguing for Rhodospirillaceae-like bacteria. Rhodospirillaceae are known as photoautotrophic, photoheterotrophic, and chemoheterotrophic bacteria [69], but not as chemoautotrophs. Since the sulfite dehydrogenase (soeA) subunit was also annotated, the Dongiales-affiliated MAG may be a chemolithoautotrophic bacterium that fix CO 2 and uses sulfur compounds as electron sources. Another interesting phylum encoding autotrophic features was Ca. KSB1 phylum, which currently consists of a single class termed UBA2214. The two Ca. KSB1-related MAGs recovered in this study encode oxic-(CBB) and anoxic (WL) CO 2 fixation potential. Both MAGs also seem to carry nitrate reduction by DNRA and denitrification, which can serve as a nitrogen retainer (DNRA) and a nitrogen remover (denitrification) in the environment. N 2 and CO 2 fixation prediction obtained by MAGs functional analysis was mostly reinforced by metagenomic read mapping on the marker genes (Additional file 5: Figs. S2, S3). This reinforcement could also be considered a validation of the obtained MAGs-based metabolic prediction even if a low binning efficiency (percentages of reads mapped to MAGs) was the case for some samples.
The metabolic interactions were investigated to extend our understanding of the metabolic potential in sediment-associated microbial communities. Microbial communities in the lower, wet and upper dry galleries were analyzed using MAG-based metabolic models (meta-GEM). The community metabolic modeling approach using metaGEM reconstruction and in silico simulation is only recently applied in different fields [70][71][72][73][74]. The community MRO and MIP distribution patterns support the expectation of lower nutrient availability in the dry gallery vs. the wet one. In scant nutrient environments, the microbes compete over the available nutrients (high MRO), and its members might need to have complementary biosynthetic capabilities to decrease their dependency on the scarce external resources (high MIP). In caves, it is known that selfish competition for resources can be replaced by cooperative and mutualistic associations, such as the ones seen in biofilms [75], maintaining bacterial communities with diverse metabolic pathways, interdependent and cooperative [76].
On the other hand, the MRO/MIP distribution pattern supports the supposition of higher nutrient availability in the areas where invertebrates were present. Those communities (PMV1, PMV2, PMV3, PMV4) had lower MRO, and MIP than the communities where no invertebrates were identified (PMV6, PMV7, PMV8). The simulated cross-feeding interactions highlighted the key donor MAGs for each condition as Methylomirabilales and Methylocella methanotrophs and Thermoplasmatota autotrophic archaeon for the lower, wet gallery; and the autotrophic Dongiales and a Phycisphaerales-related MAG for the dry gallery. Little could be deduced for the Phycisphaerales-affiliated donor MAG. None of the marker genes for sulfur, nitrogen, methane metabolism, or CO 2 fixation were assembled or annotated. This is not surprising since the phylum Planctomycetes has the highest values (35-65%) of protein sequences with unknown functions among bacterial phyla [77].
Moreover, Phycisphaerales MAG is the second of its genus (SLJD01). The first was assembled from the surface sediments of a hypersaline soda lake in Siberia (GCA_007135295). As an observation, MAGs belonging to Ca. Patescibacteria could not be modeled or were on the receivers' side of cross-feeding interactions. This is typical for Ca. Patescibacteria members with a reduced genomes size (usually < 1.5 Mb) that lack essential biosynthetic capacities have metabolic dependence and may have a parasitic or symbiotic lifestyle [78][79][80].
In the absence of light, the chemolithoautotrophic ecosystems are fueled by the oxidation of reduced compounds such as H 2 S (HS − ), CH 4 , NH 3 (NH 4 + ), Fe 2+ and H + . The metabolic modeling community simulations pointed to distinct metabolic dependencies of simple compounds in lower and upper galleries. We postulate that the compound accessibility influences the established microbial dependencies in the environment, as H 2 S, NO 3 − (a result of ammonia oxidation), and CH 4 O (a result of methanotrophy) are more available in the lower gallery than in the upper one. Hence the organism's dependencies on those compounds (i.e., a likelihood of species A growth depending on metabolite X from species B) are established in the upper gallery, where those compounds are lacking in the environment. Similarly, in the case of O 2 , the dependencies appear in the lower gallery where O 2 concentration is reduced. Ferrous iron (Fe 2+ ) dependencies are characteristics of both conditions, probably because under oxidizing conditions, iron is found mainly in Fe(III) (oxyhydr) oxide minerals (i.e. goethite).

Conclusions
The present work addressed the diversity, biogeochemical potential and ecological interaction of sedimentassociated microbiome located near and distant from the sulfidic hydrothermal waters feeding the chemoautotrophic-based Movile Cave. Metagenomic-based approaches have indicated that the diversity of the microbiomes detected in Movile Cave sediments spans a wide taxonomic range and is likely to have a high degree of novelty. This study pinpoints chemolithoautotrophy as an essential metabolic asset in an organic carbon-poor Movile Cave environment that is not confined to the high-redox potential hydrothermal water and nearby sediments but to as far as the most distant locations in the dry gallery. This assumption is supported by the recovery of autotrophic MAGs encoding CO 2 fixation ability via at least three different pathways (CBB, rTCA, WL). Sulfur oxidation was predicted for microbiomes detected nearby the sulfidic water, whereas ammonia-oxidation might be active in cave sediments in contrast with apparently absent nitrification.
Additionally, methanotrophy has been inferred across all sampled sediments. Therefore, it seems to play a key role in the primary production along the entire Movile Cave, not only in the water proximity. Despite simulation simplifications, our modeling approach postulates that nutrient scarcity is the driving force of competitioncooperation patterns across Movile Cave. The metabolic annotation and simulations point towards the autotrophic and methanotrophic MAGs as key donors in the sediments. Cross-feeding interactions can reveal the limiting compounds in the environment and the notable differences between the lower and the upper galleries.
Our findings point to the potential ecological roles and interactions of the sediment-associated microbiome and add to the previous microbiological investigations focused on the sulfidic waters in Movile Cave, thus comprehensively expanding our understanding of the peculiar chemoautotrophically-based subterranean ecosystems. Nevertheless, prospective direct metabolic quantitative assessments adjoined by multi-omics and isolation and culturing efforts are needed to further unveil the full microbiological picture of this intriguing cave ecosystem.