Corinne Biderre-Petit

Damien Courtine

Claire Hennequin

Pierre E Galand

Stefan Bertilsson

Didier Debroas

Arthur Monjot

Cécile Lepère

Anna-Maria Divne

Corentin Hochart

A pan-genomic approach reveals novel Sulfurimonas clade in the ferruginous meromictic Lake Pavin

Keywords: Meromictic lake, freshwater, sulfur oxidation, Sulfurimonas, pan-genome

The permanently anoxic waters in meromictic lakes create suitable niches for the growth of bacteria using sulfur metabolisms like sulfur oxidation. In Lake Pavin, the anoxic water mass hosts an active cryptic sulfur cycle that interacts narrowly with iron cycling, however the metabolisms of the microorganisms involved are poorly known. Here we combined metagenomics, single-cell genomics, and pan-genomics to further expand our understanding of the bacteria and the corresponding metabolisms involved in sulfur oxidation in this ferruginous sulfide-and sulfate-poor meromictic lake. We highlighted two new species within the genus Sulfurimonas that belong to a novel clade of chemotrophic sulfur oxidizers exclusive to freshwaters. We moreover conclude that this genus holds a key-role not only in limiting sulfide accumulation in the upper part of the anoxic layer but also constraining carbon, phosphate, and iron cycling.

Introduction

Sulfur cycling is an important driver of microbial growth and influences element conversions in most environments [START_REF] Berben | Diversity and Distribution of Sulfur Oxidation-Related Genes in Thioalkalivibrio, a Genus of Chemolithoautotrophic and Haloalkaliphilic Sulfur-Oxidizing Bacteria[END_REF]. In the absence of oxygen, the sulfide produced as part of this cycle [START_REF] Booker | Sulfide Generation by Dominant Halanaerobium Microorganisms in Hydraulically Fractured Shales[END_REF][START_REF] Grein | Unifying concepts in anaerobic respiration: Insights from dissimilatory sulfur metabolism[END_REF][START_REF] Ravot | Haloanaerobium congolense sp. Nov., an anaerobic, moderately halophilic, thiosulfate-and sulfur-reducing bacterium from an African oil field[END_REF] is further oxidized mostly in environments with low oxygen levels and/or limited free transition metals. In some aquatic environments, biological rates of sulfide oxidation can be orders of magnitude higher than rates of abiotic sulfide oxidation. Microbial oxidation can thereby be of vital importance, as sulfide is potently toxic to aerobic organisms and consumption in anoxic waters can effectively prevent its efflux to the aerobic waters. Bacteria involved in sulfur oxidation (SOBs) are phylogenetically and metabolically diverse and spread across various environments.

Moreover, although sulfur-based metabolism appeared almost 3.5 billion years ago, many enzymes and enzymatic pathways involved in this process remain to be discovered [START_REF] Meier | Niche partitioning of diverse sulfur-oxidizing bacteria at hydrothermal vents[END_REF][START_REF] Nosalova | Cold Sulfur Springs-Neglected Niche for Autotrophic Sulfur-Oxidizing Bacteria[END_REF][START_REF] Wacey | Microfossils of sulphurmetabolizing cells in 3.4-billion-year-old rocks of Western Australia[END_REF].

In most anaerobic sub-surface systems, the dominant chemolithotrophic SOBs belong to the co-occurring Sulfurovum, Sulfuricurvum and Sulfurimonas-related species (Campylobacteria formely Epsilonproteobacteria). Sulfurimonas representatives thrive in a broad spectrum of environments [START_REF] Han | The globally widespread genus Sulfurimonas: Versatile energy metabolisms and adaptations to redox clines[END_REF] including pelagic redox zones [START_REF] Grote | Genome and physiology of a model Epsilonproteobacterium responsible for sulfide detoxification in marine oxygen depletion zones[END_REF][START_REF] Hamilton | Coupled reductive and oxidative sulfur cycling in the phototrophic plate of a meromictic lake[END_REF][START_REF] Henkel | Candidatus Sulfurimonas marisnigri sp. Nov. And Candidatus Sulfurimonas baltica sp. Nov., thiotrophic manganese oxide reducing chemolithoautotrophs of the class Campylobacteria isolated from the pelagic redoxclines of the Black Sea and the Baltic Sea[END_REF][START_REF] Salmaso | Effects of Habitat Partitioning on the Distribution of Bacterioplankton in Deep Lakes[END_REF], hydrothermal vents [START_REF] Hu | Sulfurimonas indica sp. Nov., a hydrogen-and sulfur-oxidizing chemolithoautotroph isolated from a hydrothermal sulfide chimney in the Northwest Indian Ocean[END_REF][START_REF] Sikorski | Sulfurimonas autotrophica type strain (OK10)[END_REF][START_REF] Takai | Sulfurimonas paralvinellae sp. Nov., a novel mesophilic, hydrogen-and sulfur-oxidizing chemolithoautotroph within the Epsilonproteobacteria isolated from a deep-sea hydrothermal vent polychaete nest, reclassification of Thiomicrospira denitrificans as Sulfurimonas denitrificans comb. Nov. And emended description of the genus Sulfurimonas[END_REF]S. Wang, Jiang, et al., 2021;S. Wang, Shao, et al., 2021), petroleum reservoirs [START_REF] Lahme | Metabolites of an Oil Field Sulfide-Oxidizing, Nitrate-Reducing Sulfurimonas sp[END_REF][START_REF] Tian | Compositions and Abundances of Sulfate-Reducing and Sulfur-Oxidizing Microorganisms in Water-Flooded Petroleum Reservoirs with Different Temperatures in China[END_REF], groundwaters [START_REF] Anantharaman | Thousands of microbial genomes shed light on interconnected biogeochemical processes in an aquifer system[END_REF][START_REF] Probst | Differential depth distribution of microbial function and putative symbionts through sediment-hosted aquifers in the deep terrestrial subsurface[END_REF] and mud volcanoes [START_REF] Ratnikova | Sulfurimonas crateris sp. Nov., a facultative anaerobic sulfuroxidizing chemolithoautotrophic bacterium isolated from a terrestrial mud volcano[END_REF]. They are particularly abundant in oxygen minimum zones (15 to 30% of all microorganisms [START_REF] Callbeck | Sulfur cycling in oceanic oxygen minimum zones[END_REF][START_REF] Van Vliet | The bacterial sulfur cycle in expanding dysoxic and euxinic marine waters[END_REF]) and deep-sea hydrothermal vents (up to 70% of all microorganisms [START_REF] Campbell | The versatile ε-proteobacteria: Key players in sulphidic habitats[END_REF][START_REF] Han | The globally widespread genus Sulfurimonas: Versatile energy metabolisms and adaptations to redox clines[END_REF][START_REF] Wang | Diversity of Total Bacterial Communities and Chemoautotrophic Populations in Sulfur-Rich Sediments of Shallow-Water Hydrothermal Vents off Kueishan Island, Taiwan[END_REF]). Their widespread occurrence is largely explained by their versatile metabolisms linked to the use of a wide variety of electron donors (e.g. sulfide, thiosulfate, S 0 , H2) and acceptors (e.g. NO3 -, O2, MnO2) and by their adaptative abilities (e.g. aerobic to anoxic growth conditions, osmotic stress, etc) [START_REF] Götz | Targeted metabolomics reveals proline as a major osmolyte in the chemolithoautotroph Sulfurimonas denitrificans[END_REF][START_REF] Molari | A hydrogenotrophic Sulfurimonas is globally abundant in deep-sea oxygen-saturated hydrothermal plumes[END_REF][START_REF] Rogge | Success of chemolithoautotrophic SUP05 and Sulfurimonas GD17 cells in pelagic Baltic Sea redox zones is facilitated by their lifestyles as K-and r -strategists: SUP05 and Sulfurimonas in sulfidic redox zones[END_REF][START_REF] Taylor | Temporal shifts in dominant sulfur-oxidizing chemoautotrophic populations across the Cariaco Basin's redoxcline[END_REF][START_REF] Zeng | Microorganisms from deep-sea hydrothermal vents[END_REF].

In addition to chemolithoautotrophic oxidation of sulfur compounds, they also perform chemolithoautotrophic disproportionation of thiosulfate and elemental sulfur that may represent a previously unrecognized primary production process (S. [START_REF] Wang | Disproportionation of Inorganic Sulfur Compounds by Mesophilic Chemolithoautotrophic Campylobacterota[END_REF]. Sulfurimonas spp. have mostly been isolated from hydrothermal vents and marine sediments.

However, high-scale DNA sequencing showed that they are widespread in fresh waters such as aquifers [START_REF] Anantharaman | Thousands of microbial genomes shed light on interconnected biogeochemical processes in an aquifer system[END_REF][START_REF] Probst | Differential depth distribution of microbial function and putative symbionts through sediment-hosted aquifers in the deep terrestrial subsurface[END_REF][START_REF] Tian | Compositions and Abundances of Sulfate-Reducing and Sulfur-Oxidizing Microorganisms in Water-Flooded Petroleum Reservoirs with Different Temperatures in China[END_REF], springs [START_REF] Hahn | Microbial Diversity and Sulfur Cycling in an Early Earth Analogue: From Ancient Novelty to Modern Commonality[END_REF], and lakes [START_REF] Diao | Spatio-temporal dynamics of sulfur bacteria during oxic-Anoxic regime shifts in a seasonally stratified lake[END_REF][START_REF] Fang | Compositional and Metabolic Responses of Autotrophic Microbial Community to Salinity in Lacustrine Environments[END_REF], especially in the anoxic zone of meromictic (permanently stratified) systems [START_REF] Berg | Microbial diversity involved in iron and cryptic sulfur cycling in the ferruginous, low-sulfate waters of Lake Pavin[END_REF][START_REF] Block | Vertical structure of the bacterial diversity in meromictic Fayetteville Green Lake[END_REF][START_REF] Hamilton | Coupled reductive and oxidative sulfur cycling in the phototrophic plate of a meromictic lake[END_REF][START_REF] Salmaso | Effects of Habitat Partitioning on the Distribution of Bacterioplankton in Deep Lakes[END_REF]. In the meromictic Lake Pavin, the co-occurring Sulfuricurvum and Sulfurimonas-related species comprise up to 12% of the total bacterial community at the oxicanoxic boundary [START_REF] Biderre-Petit | Identification of sulfur-cycle prokaryotes in a low-sulfate lake (Lake Pavin) using aprA and 16S rRNA gene markers[END_REF]. However, compared to most other euxinic meromictic lakes (both anoxic and sulfidic), Lake Pavin contains very low free sulfide concentrations (<15µM) in the upper part of its monimolimnion, thus representing a good model for studying the ecology and metabolic traits of Sulfurimonas.

The small surface (44 ha) and great depth (92m) of Lake Pavin impede light penetration below the oxycline, resulting in the fixation of carbon predominantly by chemolithoautotrophic microorganisms in the anoxic bottom layer (i.e. monimolimnion). In this sulfide and sulfatepoor (<20 µm) water mass, about 80% of total sulfide detected is in the electroactive FeS colloidal form [START_REF] Bura-Nakić | Electrochemical and Colorimetric Measurements Show the Dominant Role of FeS in a Permanently Anoxic[END_REF]. [START_REF] Berg | Microbial diversity involved in iron and cryptic sulfur cycling in the ferruginous, low-sulfate waters of Lake Pavin[END_REF] showed that the cryptic sulfur cycle in these waters is tightly interacting with the biogeochemical iron cycle, suggesting that the microorganisms involved in these two cycles contribute to mineral formation through biomineralization processes [START_REF] Miot | Mineralogical Diversity in Lake Pavin : Connections with Water Column Chemistry and Biomineralization Processes[END_REF]. The bacteria and metabolisms involved remain however poorly known.

The main goal of our study was to test the hypothesis that Sulfurimonas strains that thrive in Lake Pavin are specific to freshwater systems, and act as key-players in tightly linked sulfur, carbon, phosphate, and iron cycles. We therefore investigated the distribution of Sulfurimonas spp. and their adaptation processes along the anoxic water column. Ribotypes and metabolisms of this genus were studied via metabarcoding, genome-resolved metagenomics, and single cell genomics, and compared with publicly available data from all environments shown to host Sulfurimonas spp by using a pan-genomic approach.

Materials and methods

Sample collection

Water samples for Single Amplified Genomes (SAGs) were collected from Lake Pavin on December 2015 at 65 and 80m depth, while samples for metagenomes were collected on October 2015 at 50, 65, 70 and 80m. All samples were recovered from a platform positioned above the deepest point of the lake (i.e. 92m deep) using a Van Dorn bottle. Samples for SAGs were cryopreserved in 1X TE, 5% glycerol (final concentration) for later cell-sorting, and biomass for metagenome sequencings (called MetaG2, N=4) were concentrated by filtration onto 0.22μm pore size polycarbonate (47-mm diameter, TSTP, Millipore, Billerica, Massachusetts, USA) and stored at -80°C for later DNA extraction.

Single amplified genomes (SAGs)

Single-cell sorting, whole genome amplification and SAG screening were mainly performed as previously described [START_REF] Alneberg | Genomes from uncultivated prokaryotes: A comparison of metagenome-assembled and single-amplified genomes[END_REF], for details see Supplemental material). Thirtyone SAGs that carried 16S rRNA sequences, among which one affiliated to the genus Sulfurimonas, were selected for genome sequencing. Their DNA content was quantified with Qubit dsDNA HS assay kit (ThermoFisher Scientific, Waltham MA) and diluted to 0.2ng/µl in nuclease-free water. MiSeq sequencing libraries were prepared with Nextera XT DNA Library Preparation Kit and Nextera XT Index Kit (Illumina, San Diego CA), and run according to the manufacturer's instructions. After Qubit quantification, the libraries were pooled in equimolar concentration and sequenced on a MiSeq with the V3 2x300bp strategy. Subsequently, we used the bbduck2 and bbnorm programs from the BBTools package (release 12.10.2015) (http://jgi.doe.gov/data-and-tools/bbtools/) to remove sequencing adapters and discard the low quality reads (option: trimq = 28, maq = 20, maxns = 0, and minlength = 25) as well as to normalize the quality-controlled reads (parameters target = 40, and mindepth = 0). The remaining reads were then assembled with SPAdes v3.5.0 with --careful and --sc options [START_REF] Bankevich | SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing[END_REF]. Completeness and contamination were estimated by CheckM v1.0.5 [START_REF] Parks | CheckM : Assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes[END_REF]. Gene prediction for the SAGs was performed using Prodigal V2.6.3 [START_REF] Hyatt | Prodigal: Prokaryotic gene recognition and translation initiation site identification[END_REF].

Metagenome sequencing, assembly, and binning

For MetaG2 metagenomes (N=4), genomic DNA was extracted using a standard phenolchloroform method as previously described (Biderre-Petit, Jézéquel, et al., 2011). DNA libraries were prepared and sequenced by GenoScreen (Lille, France), with an Illumina HiSeq2500 sequencer to generate 2x150bp paired-end reads. After sequencing, reads were preprocessed to remove Nextera adapters using the bbduck program as described above. Prior to assembly, raw paired-end reads were filtered using Trimmomatic v.0.33 (option: LEADING=28, TRAILING=28, SLIDINGWINDOW=4:15 and MINLEN=30) [START_REF] Bolger | Trimmomatic: A flexible trimmer for Illumina sequence data[END_REF].

Then, IDBA-UD v1.1.1 was used to assemble each sample [START_REF] Peng | IDBA-UD: A de novo assembler for singlecell and metagenomic sequencing data with highly uneven depth[END_REF] with the default parameters and both pair-end reads (-r entry) plus single-end reads (--long entry).

To further assemble additional (low coverage or less abundant) contigs, we performed an additional co-assembly of all Lake Pavin metagenomes samples, i.e. MetaG2 obtained in the present study and MetaG1 obtained from a previous study [START_REF] Biderre-Petit | New insights into the pelagic microorganisms involved in the methane cycle in the meromictic Lake Pavin through metagenomics[END_REF]. The MetaG1 (N=4) and MetaG2 (N=4) datasets were co-assembled using MEGAHIT v.1.1.1 with the 'meta-sensitive' preset option (D. [START_REF] Li | MEGAHIT: An ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph[END_REF] followed by the filtering of the resulting contigs by length (>2.5kb, i.e. the minimum contig length generally used for the initial binning [START_REF] Kang | MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities[END_REF]).

To bin the contigs, the process was as follows: (i) metagenomic reads were individually mapped against contigs using the Burrows-Wheeler-Algorithm (BWA) v.0.7.17 (Li and Durbin 2009) and stored as BAM files using Samtools, (ii) each BAM file was processed by anvi-profile as implemented in anvi'o to get the coverage and abundance information for each contig, (iii) anvi-merge was used to combine single profiles into a merged profile database, (iv) contigs were clustered via the automatic binning algorithm CONCOCT v.1.0.0 [START_REF] Alneberg | Binning metagenomic contigs by coverage and composition[END_REF] with a maximal number of clusters set to 100 to minimize 'fragmentation errors ', finally (v) the resulting bins were manually curated through the Anvi'o interactive interface using anvirefine v6.1 [START_REF] Eren | Anvi'o: An advanced analysis and visualization platform for 'omics data[END_REF] and as detailed by the Anvi'o tutorial (https://merenlab.org/2015/05/11/anvi-refine/). The manually curated bins will hereafter be referred to as Metagenome-Assembled Genomes (MAGs). MAG abundance across MetaG1 and MetaG2 datasets was assessed with the anvi-summarize program from Anvi'o.

Quality evaluation and MAG taxonomic analysis

The SAGs and MAGs genome quality was evaluated using both Anvi'o and the 'lineage workflow' of CheckM v1.0.12 [START_REF] Parks | CheckM : Assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes[END_REF]. Completeness/completion and contamination/redundancy metrics were then used to compute the mean completeness (mcpl) and the mean contamination (mctn) values to classify genome quality according to established community standards [START_REF] Bowers | Minimum information about a single amplified genome (MISAG) and a metagenomeassembled genome (MIMAG) of bacteria and archaea[END_REF]. Finally, quality scores were calculated as follows: Q-score = mcpl - 5 × mctn; and Q′-score = mcpl - 5 × mctn + mctn × (Checkm strain heterogeneity)/100 + 0.5 × log [N50].

The Genome Taxonomy Database (GTDB, release r202) [START_REF] Parks | GTDB: An ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy[END_REF] and associated taxonomic classification toolkit (GTDB-Tk, v1.7.0, default parameters) [START_REF] Chaumeil | GTDB-Tk: A toolkit to classify genomes with the Genome Taxonomy Database[END_REF] were adopted to perform MAGs taxonomic classification on the basis of a set of 120 bacterial and 122 archaeal gene markers.

To produce Sulfurimonas MAG/SAG hybrid assemblies of superior contiguity and accuracy, we tested two different pipelines allowing genome assembly, reconciliation, and merging, namely CISA v1.3 [START_REF] Lin | CISA: Contig Integrator for Sequence Assembly of Bacterial Genomes[END_REF][START_REF] Lin | CISA: Contig Integrator for Sequence Assembly of Bacterial Genomes[END_REF]GARM v0.7.5 (Soto-Jimenez et al., 2014). Hybrid genome quality was evaluated as described before for SAGs and MAGs.

Sulfurimonas MAGs and SAG abundance

Prediction of 16S rDNA genes was performed on Sulfurimonas SAG and MetaG2 contigs, using HMM search [START_REF] Huang | Identification of ribosomal RNA genes in metagenomic fragments[END_REF] with the SILVA database [START_REF] Yarza | Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences[END_REF] in rna_hmm tools (https://github.com/meb-team/Public_Tools). Predicted sequences were checked with BLASTN+ v.2.2.31+ [START_REF] Altschul | Basic local alignment search tool[END_REF] searches against the SILVA database (option : perc_identity 97% and evalue 1e -10 ) [START_REF] Yarza | Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences[END_REF]. 16S rDNA sequences were clustered using cd-hit-est (v4.6) (-c 0,95) (W. Li & Godzik, 2006, p.). Quality-controlled metagenomic reads from all samples were individually mapped against the clustered sequences. Mapping was performed using BWA (v.0.7.17) with default parameters and relative abundances were calculated with BAM-Tk (https://github.com/meb-team/BAM-Tk).

Comparative genomic analysis of Sulfurimonas MAGs, SAGs and reference genomes

For performing a comparative genomic analysis between the Sulfurimonas MAGs and SAG (identified here) and the reference genomes, we first downloaded the 17 Sulfurimonas genomes available from NCBI Genbank (December 2022). These genomes include 15 species among which one is represented by two strains (S. xiamenensis 1-1N and S. xiamenensis CVO), as well as one strain (Sulfurimonas sp. SWIR-19) (Table S1). We onwards refer to these 17 datasets as Sulfurimonas reference genomes for convenience. Sulfurimonas sp. HSL-3221 (CP087998.1) was excluded from the study because of its high GC% (>57%) and low 16S rRNA gene sequence similarity with Sulfurimonas representatives (GC% ranging from 33.3 to 38.8%; up to 93% 16S rRNA identity with S. lithotrophica), suggesting that it belongs to another genus.

We also downloaded 55 MAGs and one additional SAG affiliated with Sulfurimonas (completness >70%, redundancy <4%, and GC% <40%) from NCBI genome database and IMG/M (https://img.jgi.doe.gov) web portals (Table S2).

The average nucleotide identity (ANI) and amino acid identity (AAI) enabled delineation of operational units at the genome level. ANIb was calculated using pyANI v0.2.12 [START_REF] Pritchard | Genomics and taxonomy in diagnostics for food security: Soft-rotting enterobacterial plant pathogens[END_REF]. It was estimated by aligning fragments of 1020nt [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF] with BLASTN+ [START_REF] Altschul | Basic local alignment search tool[END_REF][START_REF] Camacho | BLAST+: Architecture and applications[END_REF] and averaging the sequence identity between pairs of genomes. AAI was calculated using the program CompareM v0.1.2 (https://github.com/dparks1134/comparem) and the parameters aai_wf --proteins -sensitive. The matrices were then plotted as a heatmap and clustered using clustermap from the Seaborn statistical data visualization package, v0.11.2, using default parameters [START_REF] Waskom | seaborn: Statistical data visualization[END_REF].

Analyses of acquired resistance genes were performed by ResFinder v4.1 [START_REF] Zankari | Identification of acquired antimicrobial resistance genes[END_REF].

Insertion sequence (IS) elements were manually annotated according to the ISfinder database (http://www-is.biotoul.fr/). Icefinder v2 (M. Liu et al., 2019) was used to identify putative integrative and conjugative elements (ICE) and integrative mobilizable elements (IME). Finally, VFDB (B. Liu et al., 2019) was used, with the genus Helicobacter as reference, to highlight putative virulence factors and TAfinder v2 [START_REF] Xie | TADB 2.0: An updated database of bacterial type II toxin-antitoxin loci[END_REF] with the option "4.2 unannotated complete genome" to identify the presence of toxin-antitoxin (TA) loci. Contigs were also manually analyzed to complete automatic annotations performed by the different tools, using both BLAST searches and predictions performed by Prodigal [START_REF] Hyatt | Prodigal: Prokaryotic gene recognition and translation initiation site identification[END_REF].

Pan-genomics analysis

The pan-genome was computed to identify the clusters of shared genes (GC) according to amino acid sequence homology. Firstly, anvi'o v7.1 [START_REF] Eren | Anvi'o: An advanced analysis and visualization platform for 'omics data[END_REF] was used to determine the pan-genome for the 17 complete reference genomes using Diamond (--sensitive) for sequence alignment (Buchfink et al., 2021) and the Markov cluster algorithm (MCL) for clustering (Van Dongen, 2008). The other parameters were kept as default except the inflation, set to 8. Gene collections were classified as follows: core genome included GCs present in all genomes; accessory genome included two collections -dispensable (present in at least two genomes, hereafter referred to as shell GCs) and unique (singletons, present in just one genome, hereafter referred to as cloud GCs). Next, these results were used to assign a category to the proteins predicted from the Lake Pavin SAG and MAGs assemblies using a BLASTP v2.9.0+ sequence alignment. Proteins without BLASTP hit were compared in an allagainst-all BLASTP, e-value <10 -5 and query coverage >80%, to extract GCs present only in a unique Lake Pavin strain. GC functional annotation was performed using eggNOG-mapper (version emapper-2.1.9) based on eggNOG orthology data [START_REF] Huerta-Cepas | eggNOG 5.0: A hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses[END_REF].

Phylogeny

A phylogenomic tree was generated with the program GToTree v1.8.2 (Lee, 2019) using a set of 74 markers spanning all bacterial genomes. The representative Sulfurimonas genomes were downloaded from NCBI RefSeq (Tables S1 andS2). GToTree identified target genes from genomes and MAGs using HMMER3, aligned each gene set using MUSCLE program v5.1 [START_REF] Edgar | Muscle5: High-accuracy alignment ensembles enable unbiased assessments of sequence homology and phylogeny[END_REF] and performed automated trimming with trimAl v1.4.1 ("-automated1") [START_REF] Capella-Gutiérrez | trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses[END_REF]. Individual gene alignments were concatenated to construct a species tree using FastTree 2.1.11 the default parameters [START_REF] Price | FastTree 2 -Approximately Maximum-Likelihood Trees for Large Alignments[END_REF].

Protein sequences to the NapA, PiT, SulP, YeeE, ferritin, sulfide: quinone reductases, and DsrE proteins of Sulfurimonas genomes were collected using accession numbers found in literature and for those not present in the literature, by searching the NCBI non-redundant database by protein BLAST search. A phylogenetic tree was generated for each protein, from aligned protein sequences (MUSCLE program [START_REF] Edgar | MUSCLE: a multiple sequence alignment method with reduced time and space complexity[END_REF])) using MEGA software package version 7.0 [START_REF] Kumar | MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets[END_REF]. Maximum likelihood (ML) method was used based on the Whelan and Goldman model (JTT model, discrete Gamma distribution (5 categories)) [START_REF] Whelan | A General Empirical Model of Protein Evolution Derived from Multiple Protein Families Using a Maximum-Likelihood Approach[END_REF], with a bootstrap iteration of 1000. All positions with less than 95% site coverage were eliminated.

Results

Reconstruction of Sulfurimonas genomes from Lake Pavin and species delineation

The metagenome assembly from all depths represented ≃1000 bins (Table S3) among which two binned MAGs were affiliated to the genus Sulfurimonas, one (Bin58) being >94% complete and high-quality and another (Bin76), >57% complete and medium quality (Table S4). One SAG (733-F2, 27.5% completeness, 0.2% single-marker-gene redundancy) recovered from the water column was also affiliated to this genus (Table S5). Across the genome regions of Bin58 and SAG 733-F2 showing homology, the within-cluster sequence identity was close to or 100% (sequence inversions were observed at some contig extremities), suggesting that both represent the same genomic population (hereafter referred to as LP strain 1). However, merging SAG 733-F2 and Bin58 genomic data led to lower estimates of completeness, higher contamination, and lower Q-score (Table S6), these two genomes were therefore kept separate in further analyses. For Bin76, the within-cluster sequence identity dropped below 85%, suggesting that this MAG corresponds to a second strain (referred to as LP strain 2).

Bin58 and SAG 733-F2 had ANIb and AAI of 99.8% and 98.7% respectively, and >99% with MAG LacPavin_0419_WC90 (Jaffe et al., 2022). The corresponding numbers for Bin76 fall below 84%, while it shares >97% with MAG LacPavin_0818_WC45 (Jaffe et al., 2022) (Fig. 1). These results confirm that Lake Pavin hosts two distinct Sulfurimonas species (threshold used for species demarcation >95% [START_REF] Kim | Towards a taxonomic coherence between average nucleotide identity and 16S rRNA gene sequence similarity for species demarcation of prokaryotes[END_REF][START_REF] Thompson | Microbial genomic taxonomy[END_REF]). Additionally, these values did not exceed 83% when compared to publicly available data for the genus Sulfurimonas (Fig. S1 andS2), suggesting that LP strains 1 and 2 represent novel species. This conclusion was supported by analyses of the 16S rRNA sequence from SAG 733-F2 with the closest database hit being S. gotlandica at 96.3% sequence identity. Moreover, the phylogenomic analysis assigned the LP strains to a clade which includes mostly strains from natural freshwater systems (Fig. 2A).

Unlike SAG 733-F2, no 16S rRNA gene was found in the recovered MAGs. By contrast, we detected three 16S rRNA sequences in the metagenome assemblies that exhibited 100% homology with the one from SAG 733-F2 (Fig. 2B). Regarding 5S rRNA genes, a copy was present in Bin58. We also isolated contigs with 5S rRNA genes from metagenome assemblies that were 100% identical between each other and with Bin58. These rRNA genes were flanked by various genomic regions (Fig. 2B).

Abundance profiles based on mapping reads from the Lake Pavin metagenomes revealed a similar vertical pattern for both SAG 733-F2 and Bin58, with a maximum relative abundance between 65m and 70m (Fig. 3A andB). A similar pattern was observed for a 16S rRNA OTU with 100% sequence similarity to the SAG 733-F2 (i.e. Sulfu. 1, Fig. 3C) based on amplicons from both DNA and RNA samples previously obtained for the Lake Pavin water column [START_REF] Biderre-Petit | Distribution of Dehalococcoidia in the Anaerobic Deep Water of a Remote Meromictic Crater Lake and Detection of Dehalococcoidia-Derived Reductive Dehalogenase Homologous Genes[END_REF]. This OTU also showed a significant abundance in the upper part of the anoxic zone (60-63m), suggesting that LP strain 1 (i.e. SAG 733-F2/Bin58) expand across a wider zone along the depth profile. In comparison, Bin76 had a marked peak in relative abundance at 65m (Fig. 3A). Its distribution pattern was similar to the one of the second most abundant Sulfurimonas OTU (Sulfu. 2, Fig. 3C), suggesting that they represent the same phylotype.

Among the 17 Sulfurimonas reference genomes, S. sediminis S2-6 and S. sp. strain SWIR-19 have ANIb and AAI ≥ 96% (Fig. 1), and 16S rRNA gene sequences that are 99.3-99.6% identical (considering only the most similar copies of SWIR-19, Table S7). This suggests that they belong to the same species. The question of heterotypic synonymy for S. sediminis S2-6 and S. hydrogeniphila might also arise. Indeed, they share similar genome wide identity (94% AAI, 92% ANI) and 16S rRNA genes that are 99.3% identical (Table S7) along with many shared biochemical and physiological features (S. Wang, Jiang, et al., 2021;S. Wang, Shao, et al., 2021).

Metabolic potential of the LP strains

Lake Pavin SAG and MAGs were assessed for the presence of specific genome-encoded functional traits involved in numerous processes, including mediation of carbon, nitrogen, sulfur, and dihydrogen (H2) cycles, oxidative stress, motility, vitamin biosynthesis and transport (Table S8, Fig. 4).

Sulfur metabolism.

Bin58 had the potential to use bisulfide (HS -) as electron donor via either sulfide:quinone reductase (Sqr) or flavocytochrome c sulfide dehydrogenase (Fcc) pathways. Like in most Sulfurimonas representatives, Type II (SqrB), IV (SqrD) and VI (SqrF) (Fig. S3, Table S9) were identified, with two copies of Type II as in S. marina, S. hongkongensis and S. gotlandica.

Moreover, the presence of a FccB homolog suggests that this pathway could also be used to catalyze initial sulfide oxidation (Fig. 5A). The less complete Bin76 genome had only one Type II and one Type IV Sqr. Only the soxCDY2Z2H gene cluster for oxidation of reduced sulfur compounds was found in both LP strains. In S. xiamenensis CVO and S. denitrificans, that display this phenotype, S 0 oxidation correlates to high expression of an OprD-like porin (CVO_00715 and Suden_1917, respectively, both sharing 57% amino acid sequence identity), for which a homolog (55% amino acid similarity) was found in Bin58. The elemental sulfur reduction in LP strains could also be mediated by the polysulfide reductase-like molybdoenzymes (psrBCA) (Fig. 5A, Table S8).

Among the three major protein families identified for the transport of oxidized sulfur compounds, two (SulP and PiT) clustered with clades involved in the bicarbonate:Na + symport and inorganic phosphate transport, respectively (Fig. S4). The only similar sequence to our PiT protein was identified in a MAG (90% protein identity, accession number OHE15791.1). By contrast, the third family of transporters, encoding the YeeE/YedE proteins, clustered with thiosulfate transporters [START_REF] Tanaka | Crystal structure of a YeeE/YedE family protein engaged in thiosulfate uptake[END_REF] and is predicted to be involved in thiosulfate ion transfer into the cytoplasm. In Bin76, the yeeE genes were tandemly organized in an operon with other genes involved in sulfur metabolism, i.e. rhodanese (rhd)-like and thioredoxin genes (Fig. S5). Only Ca. S. marisnigri contained yeeE homologs with a similar tandem organization (HUE87_05335 and HUE87_05340, > 72% amino acid sequence identity). Finally, TauE-like exporters might transport sulfite from cytoplasm to the periplasm [START_REF] Gwak | Sulfur and methane oxidation by a single microorganism[END_REF].

Both LP strains possessed a single homolog of the DsrE-like protein (Fig. S6A), involved in cytoplasmic sulfur trafficking during the sulfur oxidation. By comparison, zero to five DsrE-like proteins were detected in the full range of Sulfurimonas genomes and MAGs, with the lowest number typically found in freshwater strains (Table S10). Among all proteins identified in Sulfurimonas genomes, only one grouped in a cluster previously determined by L.-J. [START_REF] Liu | Thiosulfate Transfer Mediated by DsrE/TusA Homologs from Acidothermophilic Sulfur-oxidizing Archaeon Metallosphaera cuprina[END_REF], i.e. DsrE2A in a MAG (SAMN18877703), while the others formed eight additional clusters, herein referred to as DsrEa to DsrEh (Fig. S6B). Their distribution may not be random, as some of them are absent in for instance the freshwater strains (Table S10). The LP strain DsrE proteins belonged to the cluster DsrEa and contained three strictly conserved cysteine residues (Fig. S6A).

Regarding other cytoplasmic processes, a sulfur metabolizing pathway for assimilatory sulfate reduction was also present in the LP strains with homologs to Cys proteins (CysD, E, K, I) and a sulfate adenylytransferase (Sat, Fig. 5A) [START_REF] Han | The globally widespread genus Sulfurimonas: Versatile energy metabolisms and adaptations to redox clines[END_REF].

Nitrogen metabolism.

All components required for the denitrification pathway were detected in LP strains (Fig. 5B, Table S11), with a periplasmic nitrate reductase catalytic subunit (NapA) organized in a well conserved operon. Only S. pluma is currently known to lack this operon [START_REF] Molari | A hydrogenotrophic Sulfurimonas is globally abundant in deep-sea oxygen-saturated hydrothermal plumes[END_REF].

Moreover, as previously shown, NapA from hydrothermal vent strains represent a monophyletic group distinct from genomes obtained from other ecosystems (Fig. S7). The LP strains also contain genes encoding ammonium (atmB1 and atmB2 genes) and nitrate (narK homolog) uptake. These compounds can subsequently be funneled into amino acid metabolism via the ferredoxin-nitrite reductase (nirA) and the GS-GOGAT/GDH pathways (Fig. 5B). Unlike Bin58 and most other Sulfurimonas genomes, Bin76 contains a high-affinity nitrate ABC transporter (NrtABC) for nitrate/nitrite uptake for which homologs have only been found in Ca. S. marisnigri, S. litotrophica, S. marina, S. pluma and S. gotlandica. Bin76 genes also encode key proteins for N2 fixation (nifHD, nifA, nifT) with homologs only detected in Ca. S. marisnigri, Ca. S. baltica, S. litotrophica and S. marina. However, while both Candidatus strains can fix N2 [START_REF] Henkel | Candidatus Sulfurimonas marisnigri sp. Nov. And Candidatus Sulfurimonas baltica sp. Nov., thiotrophic manganese oxide reducing chemolithoautotrophs of the class Campylobacteria isolated from the pelagic redoxclines of the Black Sea and the Baltic Sea[END_REF], S. marina was seemingly unable to grow with N2 as the sole nitrogen source (Z. [START_REF] Wang | Sulfurimonas marina sp. Nov., an obligately chemolithoautotrophic, sulphur-oxidizing bacterium isolated from a deep-sea sediment sample from the South China Sea[END_REF].

H2, O2, oxidative stress and central carbon metabolisms.

Most Sulfurimonas can grow with H2 as an energy source. The LP strains contained [NiFe]hydrogenases of groups 1a and 1b that mediate respiratory H2 uptake, and group 4e, representing a multisubunit membrane-bound Energy-Conserving Hydrogenase (Ech) that couples ferredoxin oxidation to H2 evolution (Table S8). Unlike group 1 [NiFe]-hydrogenases, which is ubiquitous in Sulfurimonas spp., group 4e is present only in marine pelagic species (S. gotlandica, Candidatus S. marisnigri and Ca. S. baltica) and in both LP strains.

Like other Sulfurimonas genomes [START_REF] Molari | A hydrogenotrophic Sulfurimonas is globally abundant in deep-sea oxygen-saturated hydrothermal plumes[END_REF], LP strains carry various systems for uptake and intracellular availability control of iron and manganese. Among them, a gene encodes a DNA-binding protein belonging to the ferritin-like diiron-carboxylate protein superfamily (known as Dps). This gene is found in almost half of the Sulfurimonas genomes (Fig. S8). In contrast, a second member within this gene superfamily (i.e. a gene encoding for a 24-subunit non-heme ferritin) was only found in four Sulfurimonas reference genomes (Fig. S8). In addition to the ferritin, LP genomes have genes encoding a superoxide dismutase (SOD), producing H2O2 and O2 as byproducts, and a superoxide reductase (SOR), producing only H2O2 (Table S8, Fig. 4A), both designed to specifically minimize stress and cell damage from reactive oxygen species (ROS) [START_REF] Johnson | Distribution of reactive oxygen species defense mechanisms across domain bacteria[END_REF][START_REF] Lushchak | Adaptive response to oxidative stress: Bacteria, fungi, plants and animals[END_REF].

In assessing the potential for aerobic respiration, both LP strains have a gene coding for the cbb3-like high-O2 affinity cytochrome, as all other Sulfurimonas spp. This cytochrome is active in anaerobic processes such as dinitrogen (N2) fixation and denitrification [START_REF] Bell | Active sulfur cycling in the terrestrial deep subsurface[END_REF]. In addition, two hemerythrins, non-heme intracellular O2-binding proteins, were present in LP strains, one with a single-domain and one with a multi-domain that belongs to the hemerythrin-GGDEF family (Table S8). We also detected an Aer-PAS domain primary energy sensor for motility as well genes coding for a flagellum and chemotaxis (Fig. 4A, Table S8).

Finally, both LP strains contain all essential genes for carbon fixation via the reverse tricarboxylic acid (rTCA) cycle, a common autotrophic carbon fixation pathway of all Sulfurimonas spp. (Fig. 4A).

Mobilome and defense mechanism analysis

MGE in LP strains

MGE is a generic term for a multitude of genomic sequences (plasmid, prophage, ICE, IME, insertion sequences, etc.) transmitted vertically with cell division or through horizontal transfer (HGT). Most of such genomic features were identified in the LP strain genomes via either automatic or manual annotations (Table S12). An IME was found in both LP strains while a prophage was present only in Bin58. The IME encodes a tyrosine integrase and harbors a compact mobilization module with two genes encoding relaxases (one of which is a canonical MobC) and one gene encoding a non-canonical coupling protein (Type IV secretory system conjugative DNA transfer family protein) (Fig. 6). Only one MAG affiliated to Arcobacter sp., reconstructed from cold oxic subseafloor [START_REF] Tully | A dynamic microbial community with high functional redundancy inhabits the cold, oxic subseafloor aquifer[END_REF], have been shown to host a similar element with a conserved gene order (Fig. 6). With regards to the prophage, genomic elements with a similar organization and content have only been reported for Sulfurimonas MAGs (Fig. S9). Both the IME and prophage carried many flexible genes, mostly encoding proteins of unknown function. Therefore, their role in the acquisition of novel adaptive phenotypes remains unknown.

Defense systems

Among all defense mechanisms identified, Toxin-Antitoxin (TA) systems were prominent, representing small genetic modules with implications in a multitude of bacterial lifestyle adaptation [START_REF] Fraikin | Type II Toxin-Antitoxin Systems: Evolution and Revolutions[END_REF]. Among the six TA types (types I-VI) currently classified, type II, that function via protein-protein interactions, is the most well characterized and the most prevalent across bacterial genomes [START_REF] Kamruzzaman | Biological Functions of Type II Toxin-Antitoxin Systems in Bacteria[END_REF][START_REF] Yeo | GNAT toxins of bacterial toxin-antitoxin systems: Acetylation of charged tRNAs to inhibit translation: GNAT toxins of bacterial toxin-antitoxin systems[END_REF]. Here, we identified a repertoire of 51 different type II TA systems in the Sulfurimonas genomes involving 16 toxin families that are mostly involved in translational arrest (Table 1). Each of the 17 reference genomes had 1 to 28 TA pairs through its chromosome. For those genomes featuring plasmids, one (S. paralvinellae) and five additional loci (S. aquatica) were identified. For LP strains, we identified 12 TA pairs in Bin76, nine in Bin58 and eight in SAG 733-F2, as well as the Clp proteases involved in their degradation (Table S12). Overall, the TA pattern was similar for all species of a specific environment even if the TA pair number per species may varied (Table 1).

In contrast, both TA number and pattern are kept across strains of a same species, regardless if they are from the same (S. sediminis and SWIR-19) or different habitats (S. xiamenensis strains). Interestingly, two species (S. paralvinellae and S. marina) have only one TA pair.

However, compared to the species that live in the same environment, these species seem to be restricted to more specific niches. Indeed, S. paralvinellae, in hydrothermal vents, lives associated with a host, unlike any of the other species and S. marina, in marine sediments, thrives at 2000m depth while the other species were retrieved from less than 250m deep. In addition to type II, a putative type IV TA system (AbiE), widespread in bacteria in which it ensures phage resistance (Dy et al., 2014), was present in seven genomes with an occurrence per genome ranging from one to five (Table 1).

Resistome

The resistome analysis identified resistance to several antibiotic classes within the LP strain 1 genomes including beta-lactams, quinolones, macrolides and fusaric acid (FA). We found a gene encoding an AmpC beta-lactamase (Table S12), with the only Sulfurimonas homolog being a gene found in a MAG from groundwater (Finland, SAMN12221366) that exhibits the same genomic context, i.e. associated with two amino acid ABC transporter substrate-binding proteins and a HAMP domain-containing histidine kinase. Interestingly, LP strain 1 is the only known Sulfurimonas representative harboring genes involved in FA resistance. The fuaABC operon consists of a transcription factor of the family MarR, an inner membrane efflux pump, a periplasmic efflux transporter and an outer membrane protein. Our strain also has other efflux pump systems belonging to MATE, ABC (MAcAB-TolC) and RND (AcrAB-TolC) families that could contribute to the generation of resistant phenotypes (Table S8).

Virulome

Virulome prediction identified genes involved in adherence, lipid A synthesis and biofilm formation (Table S12). LP strain 1 features all traits required for biofilm formation capacity (chemotaxis, flagellum-assisted motility, EPS production and export). Moreover, we identified a pgaBCD locus that is required for a biofilm adhesin synthesis (poly-β-1,6-GlcNAc or PGA) [START_REF] Wang | The pgaABCD Locus of Escherichia coli Promotes the Synthesis of a Polysaccharide Adhesin Required for Biofilm Formation[END_REF], but where the pgaA subunit encoding the outer membrane protein is lacking. Homologous genes were only detected in S. gotlandica and S. denitrificans, also devoid of the pgaA subunit.

Pan-genome analysis of the genus Sulfurimonas

Comparative analysis of the 17 Sulfurimonas reference genomes revealed a pan-genome composed of 9021 gene clusters (GCs) among which 915 core GCs, 3508 shell GCs and 4598 cloud GCs (Fig. 7A), which is consistent with previous studies [START_REF] Molari | A hydrogenotrophic Sulfurimonas is globally abundant in deep-sea oxygen-saturated hydrothermal plumes[END_REF]S. Wang, Jiang, et al., 2021). With regards to the LP genomes (SAG 733-F2, Bin58 and Bin76), they contained 272 to 842 core GCs, 460 to 1152 shell GCs, and 57 to 192 cloud GCs (Fig. 7B).

Among the shell GCs identified in LP strains, 340 to 946 belonged to the shell GCs when analyzing the 17 Sulfurimonas reference genomes while 96 to 206 belonged to the cloud GCs (Fig. 7B). This is what is generally observed with the addition of new genomes in pan-genomic analyses: certain GCs firstly identified as cloud are subsequently shared by more than one genome after adding additional genomes resulting in their updated classification as shell GCs.

In this study, 316 distinct GCs classified as cloud GCs when analyzing the 17 Sulfurimonas reference genomes were reclassified into shell GCs when adding LP genomes (Table S13). The prediction of their function via assignment to the Clusters of Orthologous Groups (COGs) classification shows that 39% (123 GCs) could not be categorized due to the absence of common features with other proteins. The 61% remaining clusters fell into 18 COGs with the categories « function unknown » (S, 18%), « signal transduction mechanisms » (T, 27.3%) and « cell wall, membrane, envelope biogenesis » (M, 10.8%) being the biggest three clusters (Table S13).

Comparative analysis based on the cloud GCs identified in LP strains revealed 104 GCs found only in SAG 733-F2, with 127 and 38 exclusive to Bin58 and Bin76, respectively (Fig. 7C). The vast majority of these were not categorized (76%, 66.2% and 71.1%, respectively). For the classified ones, nearly a quarter fall into the categories « S » and « R (only general predicted function) » (Table S14). The other most common COGs were related to « mobilome (X) » and « defense mechanisms (V) » (Fig. 8). For LP strain 1, proteins related to « amino-acid transport and metabolism (E) », « replication, recombination, and repair (L) » and « T » were also identified in lower proportions. In addition, a few cloud GCs (22 in total, Fig. 7C) were shared by both strains, most of which (95.4%) had a COG classification. The categories « inorganic ion transport and metabolism » (P), « cell cycle control, cell division, chromosome partitioning (D) », « intracellular trafficking, secretion, and vesicular transport (U) » and M were the main ones (Fig. 8, Table S9). Among these proteins, some could be relevant for the LP strain fitness, as for instance, a chromate-resistance protein, involved in the reduction of chromate accumulation (element with toxic and mutagenic effects), porins of the OprB family, central in carbohydrate transport, and the Min system proteins (MinE, MinC, MinD) (Fig. 4). LP strains differ from all Sulfurimonas genomes by the presence of the MinE protein (Fig. S5) that jointly with MinD enables oscillatory behavior of the Min system. In its absence, the system remains static [START_REF] Taviti | Bacterial Min proteins beyond the cell division[END_REF]. In addition to its involvement in cell division, this system is also crucial for cellular processes such as bacterial motility, RNA decay, colonization, and virulence.

Discussion

In this comparative study of Lake Pavin microbiota, we identified two uncharacterized Sulfurimonas species inhabiting the bottom anoxic zone using both MAGs and SAGs data. The clear difference we observed in level of completeness between SAG and MAGs affiliated to this genus is congruent with that was previously observed by [START_REF] Alneberg | Genomes from uncultivated prokaryotes: A comparison of metagenome-assembled and single-amplified genomes[END_REF] and attributed this to constraints related to the WGA during single-genome amplification. The MAG quality was ensured by low level of redundancy, rejecting the possibility of contaminant metagenome contigs included in the binning and assembly. Moreover, our abundance profiles suggest that the two LP strains may differ in their ecological preferences related to oxygen (O2) and electron acceptor/donor availability. Such coexistence of Sulfurimonas species within a single ecosystem with overlapping niches has never been described before.

The presence of various genes in the vicinity of the ribosomal genes for the LP Lake Pavin strain 1 implies multiple copies per genome, as observed for most Sulfurimonas members (Table S1). Furthermore, phylogenomic analysis as well as the diversity and distribution pattern of specific genes (e.g. dsrE genes) suggest that LP strains form, with other freshwater representatives, a distinct freshwater clade within this genus. Regarding all the other Sulfurimonas strains analysed in this study, S. sediminis S2-6 and S. sp. strain SWIR-19 exhibit a high degree of genome relatedness. This finding allows us to propose renaming the latter to S. sediminis strain SWIR-19, similarly to the recent renaming of S. sp. strain CVO (renamed S. xiamenensis strain CVO [START_REF] Henkel | Candidatus Sulfurimonas marisnigri sp. Nov. And Candidatus Sulfurimonas baltica sp. Nov., thiotrophic manganese oxide reducing chemolithoautotrophs of the class Campylobacteria isolated from the pelagic redoxclines of the Black Sea and the Baltic Sea[END_REF]). The high genetic similarity, also observed with the S. hydrogeniphila strain, could be explained by the fact that they all originate from the Indian Ocean (Table S1) with the most closely related S. sediminis S2-6 and SWIR-19 isolated from the same area (GPS coordinates 37.78S 49.64E and 37.28S 49.65E, respectively (S. Wang, Shao, et al., 2021)), while the third, S. hydrogeniphila, was recovered from a more distant region (6.21N 60.31E (S. Wang, Jiang, et al., 2021)). This is consistent with the identification of ecotypes specific to an environment (as defined by [START_REF] Molari | A hydrogenotrophic Sulfurimonas is globally abundant in deep-sea oxygen-saturated hydrothermal plumes[END_REF]; e.g. brackish water samples, hydrothermal plume samples, etc…) defined on the basis of oligotypes (amplicon sequence variants distinguished on the basis of information entropy determined for each nucleotide position [START_REF] Eren | Anvi'o: An advanced analysis and visualization platform for 'omics data[END_REF]).

From a functional point of view, like other Sulfurimonas-spp. from terrestrial subsurface environments, the LP strains display only the soxCDY2Z2H genes and not the soxXY1Z1AB genes [START_REF] Lahme | Comparison of sulfide-oxidizing Sulfurimonas strains reveals a new mode of thiosulfate formation in subsurface environments[END_REF]. The products of this gene cluster are involved in S 0 oxidation with thiosulfate as metabolic end product. Cell attachment to S 0 particles are required for their oxidation and different mechanisms were suggested, depending on the Sulfurimonas spp. (e.g. flagellar proteins, Type-IV pilin-like proteins). Recently, an outer membrane protein (OprD-like porin) was also implicated in the S 0 oxidation in strain CVO for which a homolog was found in LP strains. Here, the question whether this protein is directly involved in degradation of S 0 or facilitates its transport into the cell remains open. Data available for other proteins involved in the sulfur cycle are also poor. This is notably the case for the DsrE proteins for which investigations are needed to identify not only their function but also the proteins they may interact with. For instance, the DsrE2A protein, organized in a rhd-tusA-dsrE2A gene cluster, was shown to mediate sulfur transfer and sulfite oxidation in cytoplasm [START_REF] Dahl | Cytoplasmic sulfur trafficking in sulfur-oxidizing prokaryotes[END_REF][START_REF] Liu | Thiosulfate Transfer Mediated by DsrE/TusA Homologs from Acidothermophilic Sulfur-oxidizing Archaeon Metallosphaera cuprina[END_REF]. Since the gene encoding DsrEa in LP strains is in a slightly different cluster (tusA gene replaced by the protein-disulfide reductase (dsbD gene)), its true function in the cell could be cryptic. We also detected in the LP strain genomes both a PiT transporter homolog, closer to inorganic phosphate transporters than to sulfur compounds [START_REF] Marietou | Sulfate Transporters in Dissimilatory Sulfate Reducing Microorganisms: A Comparative Genomics Analysis[END_REF], and proteins involved in polyphosphate formation and degradation (Ppk). This suggests the LP strains have the capacity to accumulate polyphosphate under favorable conditions as a cellular storage of energy, analogous to the Sulfurimonas subgroup GD17/GD1T isolated from Baltic Sea [START_REF] Möller | Sulfurimonas subgroup GD17 cells accumulate polyphosphate under fluctuating redox conditions in the Baltic Sea: Possible implications for their ecology[END_REF].

The presence of various systems for intracellular iron availability control in LP strain genome suggests its potential role as electron donors. Also, though essential, iron can be a harmful element for the cell by catalyzing the formation of ROS by reacting with H2O2 in a Fenton reaction. Almost half of the Sulfurimonas genomes, including LP strains, encode a Dps-like protein known to minimize damage from such oxidative stress. Indeed, this protein quenches both Fe 2+ and H2O2 while physically protecting genomic DNA and hence protects the cell against Fenton-reaction [START_REF] Bradley | Bacterial iron detoxification at the molecular level[END_REF].

Numerous genes implicated in biofilm formation and motility are present in LP strain genomes. Some of these genes were found only in a very limited number of species (e.g. pga, hemerythrin genes). With regard to one of the two hemerythrins found, the GGDEF domain located at its C terminus is believed to be involved in the synthesis of a second messenger that contributes to the control of motility, virulence and biofilm formation [START_REF] French | Diversity and distribution of hemerythrin-like proteins in prokaryotes[END_REF][START_REF] Römling | Cyclic di-GMP: The First 25 Years of a Universal Bacterial Second Messenger[END_REF]. This suggests a link between O2 levels and the cyclic-di-GMP signaling network in LP strains. However, how the interaction of the hemerythrin-like domain with O2 alters the activity of the associated signaling domain is still not known and should be investigated. Furthermore, the importance of O2 in Sulfurimonas strains motility is reinforced by the presence, in the few strains where this hemerythrin where detected (including LP strains), of an Aer-PAS domain primary energy sensor for motility, known to respond to changes in O2 concentration. Hence, in cooperation with chemotaxis, such traits could confer an ability for the cell to position themselves in a microenvironment appropriate for their metabolic needs [START_REF] Maschmann | Redox properties and PAS domain structure of the Escherichia coli energy sensor Aer indicate a multistate sensing mechanism[END_REF].

For the PGA adhesin, since it has been hypothesized that the PgaA subunit mediates its translocation and/or docking to the cell surface, it would be interesting to investigate how PGA is released at the surface in these strains in its absence. Moreover, both LP strains and S. gotlandica all thrive in the pelagic redox zone of the water column, whether in freshwater lake or the Baltic Sea. Therefore, although this has not yet been demonstrated experimentally, capacity to live in biofilm at such redox transition zones might confer ecological advantages, such as efficient nutrient cycling, stress tolerance and resistance to grazing. As S. gotlandica is a cultivated species, it would be possible to elucidate some of these issues and determine the importance of pga and hemerythrin genes in biofilm formation.

We suggest that MGEs might play important roles in the ability of Sulfurimonas spp. to disperse. Among them, IMEs have never been reported in Sulfurimonas before and, interestingly, the gene associations within IME in LP strains differs from the mandatory ones described by [START_REF] Coluzzi | A Glimpse into the World of Integrative and Mobilizable Elements in Streptococci Reveals an Unexpected Diversity and Novel Families of Mobilization Proteins[END_REF] in Streptococci. For TAs, the wide variety of associated antitoxins suggests that each TA type II pair may function independently. Moreover, since variations in their patterns are observed both between environments and between species, further work is needed to investigate how such genomic features contribute to Sulfurimonas strain acclimation to variable conditions and particular environment. Indeed, these characteristics may, in part, explain the low number of TA pairs observed in certain strains.

The LP strains also harbor a broad arsenal of virulence and resistance determinants that may demonstrate their capacity for genetic recombination and enable rapid and efficient adaptation to the local environment. Among the resistance determinants, the FA operon shows an organization that is consistent with the finding that Gram-negative bacteria resistance to high FA would be conferred by an efflux pump [START_REF] Crutcher | Microbial Resistance Mechanisms to the Antibiotic and Phytotoxin Fusaric Acid[END_REF].

Nevertheless, the ability of our strain to grow when exposed to FA remains to be experimentally confirmed. Moreover, the similarity of the MarR sequences with those of Sulfuricurvum kujiense (Table S14) suggests an acquisition via HGT. Therefore, further studies are needed to shed light on the evolutionary and ecological forces underpinning the early evolution and differentiation of these freshwater strains, as well as their ability to disperse globally.

The upper boundary of the Lake Pavin monimolimnion, where Sulfurimonas live, is dynamic and uncertain as it responds to weather-induced mixing. Our results suggest that LP strain 1, like S. gotlandica, has a strategy based on polyphosphate storage capacity that can be used as an energy source to sustain cell integrity and reach favorable niches when substrate availability changes. This capacity to accumulate polyphosphate has also been demonstrated for magnetotactic bacteria in the monimolimnion of Lake Pavin [START_REF] Rivas-Lamelo | Magnetotactic bacteria as a new model for P sequestration in the ferruginous Lake Pavin[END_REF], suggesting that this process might be more widely used than expected in this ecosystem.

Furthermore, the LP strain 1 was detected in a zone of the water column where conditions varied from anoxic, ferruginous with elemental sulfur (≈60m, S: 3.3µM) to slightly sulfidic (≈70m, H2S: 15µM). Together, this suggests Sulfurimonas could be key-players in the interconnected cycling of sulfur, carbon, phosphate, and iron cycles in Lake Pavin. Hence this lake may be a promising environment to gain new knowledge about freshwater Sulfurimonas strains and to elucidate their interactions with the environment.

Funding

The work was supported by the Agence Nationale de la Recherche (ANR) through the project EUREKA (ANR-14-CE02-0004-01).

OTU with seed-sequence 100% identical to the SAG 733-F2 sequence; Sulfu2: the second most abundant OTU affiliated to the genus Sulfurimonas and with less than 97% identity with the SAG 733-F2 sequence. S9, and for nitrogen pathways, Table S10).  Table S1: Characteristics of the 17 Sulfurimonas reference genomes and strains from the water column of Lake Pavin.

Table S2: List of Sulfurimonas MAGs and SAGs downloaded from the public databases and used in this study.

Table S3: List of Bins produced by de novo assembly of pooled reads from all metagenomes (MetaG1 and MetaG2) collected along the Lake Pavin water column.

Table S4: Statistics on the assembly of SAG and MAGs of Sulfurimonas genus obtained from Lake Pavin, completeness and contamination estimates.

Table S5: Taxonomic affiliation of SAGs obtained at 65 and 80 m depth in Lake Pavin water column based on their 16S rRNA gene sequence (Sanger sequencing).

Table S6: Quality evaluation of hybrid assembly genomes and of Bin58 MAG and SAG_733-F2 used to produce them.

Table S7:

Multiple sequence alignement of the different 16S rRNA gene copies from Sulfurimonas sp. SWIR19, S. sediminis S2-6 and S. hydrogeniphila using Clustal Omega tool (https://www.ebi.ac.uk/Tools/msa/clustalo/). Table S8: Specific functions involved in numerous processes, including: carbon, nitrogen, sulfur, hydrogen (H2) cycling, oxidative stress, motility, vitamin biosynthesis, transport, etc, identified in Lake Pavin strains. Table S9: List of genes involved in sulfur cycling identified in Lake Pavin strains.

Table S10: DsrE-like protein diversity and distribution in Sulfurimonas strains. Table S11: List of genes involved in nitrogen cycling identified in Lake Pavin strains.

Table S12: Predicted mobilome, resistome and virulome in MAGs and SAG sequences (* genes found only in Lake Pavin strains).

Table S13: List of cloud gene clusters found both in Sulfurimonas reference genomes and Lake Pavin strains to be reclassified into shell GC and their classification into COG categories.

N A D H H 2 O 2 H 2 O c h e m o ta x is HCO 3 -+ H + CO 2 + H 2 O 2O 2 .-+ 2H + H 2 O 2 + O 2 SOD H 2 O 2 H 2 O tpx,trx Dpsl iron storage 2Fe 2+ + H 2 O 2 + 2H 2 O 2O 2 .-+ 2H + + e - H 2 O 2 SOR Fep FeoB Exb D Exb D TA M CA m la A B C A B C M a c A B p q iC

p a r a l v i n e l l a e S . a u t o t r o p h i c a S . i n d i c a S . h y d r o g e n i p h i l a S . s e d i m i n i s ‡

S . s p _ S W

I R 1 9 ‡

S . a q u a t i c a S . g o t l a n d i c a S . C a n d i d a t u s b a l t i c a S . C a n d i d a t u s m a r i s n i g r i S . m a r i n a S . d e n i t r i f i c a n s S . h o n g k o n g e n s i s S . l i t o t r o p h i c a S . x i a m e n e n s i s

Hypothetical proteins

Known proteins

Bin 58 K119_515439 (36.9 kb) 1 2 3 4 5 6 7 8 9

Figure 4 :

 4 Figure 4: Overview of the metabolism of the Lake Pavin Sulfurimonas strains inferred from their genome. A) General functions. Transporter families (ABC, MFS, P-type ATPase, symporter, etc) are represented by different colors and forms. B) Functions only carried by Lake Pavin (LP) strains (strain 1, strain 2 or both -black boxes) or shared by one or both LP strains with only one other Sulfurimonas genome (blue box).

Figure 5 :

 5 Figure 5: Overview of the putative pathways of oxidative and assimilation metabolisms for sulfur (A) and of reductive and assimilation metabolisms for nitrogen (B) within the Lake Pavin strains (for details about genes involved in sulfur pathways see TableS9, and for

Figure 6 :

 6 Figure 6: Modular structure of integrative and mobilizable element (IME) found in Sulfurimonas strains from Lake Pavin. The gene blocks conserved between SAG 733-F2 and Bin58, SAG 733-F2 and Bin76 and SAG 733-F2 and Arcobacter sp. are shown. Proteins not found in other Sulfurimonas genomes appear in bold; those found only in a few genomes, in bold italics.

Figure 7 :

 7 Figure 7: Pan-genome analysis of the 17 Sulfurimonas reference genomes used in the present study and the Lake Pavin strains. A) Flower plot displaying the numbers of core gene clusters (GCs) (in the center), flexible (shell) GCs (in the annulus) and unique (cloud) GCs (in the petals) for the 17 Sulfurimonas genomes. B) Number of core, flexible (shell), and unique (cloud) GCs for SAG 733-F2, Bin58 and Bin76 genomes. C) Venn diagram displaying the number of cloud GCs shared by lake Pavin strains or specific to each. The number of genes corresponding to the GC number is indicated in parentheses in order: SAG 733-F2, Bin58, Bin76.

Figure 8 :

 8 Figure 8: Clusters of Orthologous Groups (COGs) classification of the cloud GCs identified in Lake Pavin strains.Table 1: Classification and distribution of Type II and IV toxin-antitoxin systems in Sulfurimonas species.

Figure S1 :

 S1 Figure S1: Pairwise average nucleotide identity (ANIb) heatmap of Sulfurimonas genomes, MAGs and Bin58.

Figure S2 :

 S2 Figure S2: Pairwise average amino acid identity (AAI) heatmap of Sulfurimonas genomes, MAGs and Bin58.

Figure S3 :

 S3 Figure S3: Maximum likehood phylogenetic tree of deduced protein sequences of the sulfide:quinone reductase (Spr). Number at nodes represent bootstrap values (1 000 replicates of the original dataset), only number above 50% are shown. The scale bar represents the average number of substitutions per site. The accession number for the sequences are given at the right side of species name. The protein sequences from Lake Pavin strains are indicated by a dot for Bin58 and a circle for Bin76, with colors depending on the sqr type: type II in red (SqrB), type IV in green (SqrD) and type VI in blue (SqrF).

Figure S4 :

 S4 Figure S4: Maximum likehood phylogenetic trees of deduced protein sequences of two potential sulfate transporters (PiT and SulP). On the left: PiT protein sequence tree. The two protein clades are surrounded with different colors and given a letter label (A, B). Reference sequences: PITA, inorganic phosphate transporter (UniProt P0AFJ7), PITB: inorganic phosphate transporter (UniProt P0AFJ7), Cys Pit: sulfate transporter (BSUA_01689). On the right: SulP protein sequence tree. The three protein clades are surrounded with different colors and given a letter label (A, B, C) while separate protein groupings are denoted by subscript number (i, ii). Reference sequences: BICA (Aii), bicarbonate transporter (UniProt Q14SY0), DAU_Salm (B): dicarboxylic transporter (ACY88617), Riv1739c and SulP_Pseudo (PA1647) (C): sulfate transporter. Number at nodes represent bootstrap values (1 000 replications of the original dataset), only number above 50% are shown. The scale bar represents the average number of substitutions per site. The accession number are given for all sequences. Bin58 and SAG 733-F2 sequences are identified by a red dot.

Figure S5 :

 S5 Figure S5: Schematic representation of yee gene genomic context in the Bin76 genome. All genes potentially involved in sulfur cycling are in ochre-brown color. Each box represents genes in the contig for which a homolog was present only in a restricted number of Sulfurimonas genomes and MAGs. Their genomic context in the Lake Pavin strain 1 is also shown (Node for SAG 733-F2 and locus tag, i.e. K119_, for Bin58). Ochre arrows: genes encoding YeeE proteins; green arrows: Min operon; black blue arrow: gene encoding a high potential iron-sulfur protein.

Figure S6 :

 S6 Figure S6: Sequence alignment of DsrE homologs (A) and phylogenetic tree (B). In A, the sequences are identified by their locus tags. The active site cysteine conserved in all sequences is labeled with an asterisk and highlighted in cyan. For clarity, only those parts of the alignment are shown to contain conserved cysteine residues (N and C terminus, cysteines highlighted in yellow). In blue: sequence identifiers of Sulfurimonas genomes. Bins are highlighted in red and their motifs in bold italics. In B, the evolutionary history was inferred using the maximum likelihood method. Number at nodes represent bootstrap values (1 000 replicates of the original dataset), only number above 50% are shown. The scale bar represents the average number of substitutions per site. The accession number are given for all sequences. Bin sequences are identified by a red dot, Sulfurimonas genome sequences by a blue diamond. Subclusters DsrE2A, DsrE3A, DsrE3B, DsrE4 and DsrE5 were those previously described by Liu et al (2014). DsrE corresponds to the classical DsrE/TusD subunit from the DsrEFH operon that acts as cytoplasmic acceptor for sulfur. C  corresponds to the active site cysteine conserved in all sequences.

Figure S7 :

 S7 Figure S7: Maximum likehood phylogenetic trees of deduced protein sequences of periplasmic nitrate reductase catalytic subunit (NapA). Number at nodes represent bootstrap values (1 000 replicates of the original dataset), only number above 50% are shown. The scale bar represents the average number of substitutions per site. The accession number are given for all sequences. Each dot color corresponds to an environmental category for Sulfurimonas genomes: blue, freshwater (Bin58 sequence); yellow: marine sediment; orange: marine water; light yellow: brackish water, green: oil field water, gray: mud volcano and light brown: hydrothermal vent.

Figure S8 :

 S8 Figure S8: Maximum likehood phylogenetic tree of deduced protein sequences of the ferritin-like diiron-carboxylate protein superfamily. Number at nodes represent bootstrap values (1 000 replicates of the original dataset), only number above 50% are shown. The scale bar represents the average number of substitutions per site. The accession numbers for the sequences are given at the right side of species name. For the Sulfurimonas genomes, the ferritin belonging to the DNA-binding protein from nutrient starved cells (Dps) family is identified by a green dot while that belonging to the 24-subunit non-heme ferritin family, by a red dot.

Figure S9 :

 S9 Figure S9: Modular structure of integrated prophage found in Lake Pavin strain 1. The gene blocks conserved in the two MAGs affiliated to the Sulfurimonas genus (UBA and RIFOXYD2) are shown. Proteins found in no other Sulfurimonas genomes appear in red while those found only in a few genomes, in blue.

  Figure3

①⑧

  Figure4

  Figure 7

  Figure S3

  Sulfurimonas MAG sequences in NR database (no Sulfurimonas species or other taxa sequences)

Table S14 :

 S14 Cloud gene clusters identified in Lake Pavin strains for which a COG classification was found.

		Anvi'o MetaG1	Anvi'o MetaG2	MetaG2 16S	MetaB 16S gDNA	MetaB 16S cDNA
	4m	0	0
	50m			0.03	0.16	0.54
							0	0
	53m	0.03	0.35
							0.77	0	0	0
							6.48	0.04	1.26	0.04
							7.74	0.07	6.14	0.38
	65m	2.25	0.87	12.7	4.77	5.63
							13.08	3.64	7.2	2.22
	70m			12.66	2.57	4.55
	80m	0.21	0.37	5.5	1.13	1.47	0.78	0.71	0.94	0.39
	Figure 1 Figure 2 Bin58 Bin76	Bin58	Bin76 SAG 733-F2	Sulfu. 1	Sulfu. 2	Sulfu. 1	Sulfu. 2

Acknowledgements

The authors would like to acknowledge support from the Genomics infrastructure services at Science for Life Laboratory (https://www.scilifelab.se). Single Cell Sorting, whole genome amplification, screening and library preparation were performed at the Microbial Single Cell Genomics Facility (MSCG). Sequencing was performed by the SNP&SEQ Technology Platform.

The facilities are part of the Uppsala node at the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory. The SNP&SEQ Platform is also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. We are also grateful to the Mésocentre Clermont Auvergne University (https://mesocentre.uca.fr/) for providing help, computing, and storage resources. We thank Maria Hammond (MSCG, Upsala), Hélène Gardon (PhD, Université Clermont-Auvergne) and Vincent Dardot (M2 student, Université Clermont-Auvergne) for efficient technical assistance and expertise.

Data accessibility

The raw sequences of the metagenomic datasets MetaG1, previously described in [START_REF] Biderre-Petit | New insights into the pelagic microorganisms involved in the methane cycle in the meromictic Lake Pavin through metagenomics[END_REF], are archived in MG-RAST under accession numbers mgm4580845.3 to mgm4580848.3 for the metagenomic reads and the accession numbers mgm4557752. 3.100, mgm4557753.3.100, mgm4581104.3.100 and mgm4581106.3.100 for the assembled contigs.

Metabarcoding data are available at the NCBI SRA database under the BioProject accession number PRJNA975855 and SRA accession numbers SRX20521566 to SRX20521576. Raw sequences of the metagenomic datasets MetaG2 as well as the Sulfurimonas SAG and MAGs produced as part of the present study are available at NCBI SRA database under the BioProject accession PRJNA983208.

Conflict of interests

The authors declare that they have no known competing interests.

Author contributions AM conceived figures, involvement in writing-review and editing; CBP was responsible for conceptualization, supervision, data curation, data analysis and interpretation and writingoriginal draft; Che was involved in formal analysis and writing -review and editing; CHo was responsible for metagenome assembly and analysis, data curation, methodology and writing -review and editing; CL was involved in writing -original draft , DC participate to the pangenome analysis, data curation, methodology, writing -review and editing; PEG was responsible for sampling and funding acquisition (ANR EUREKA), writing -review and editing; DD was responsible for sampling and funding acquisition (ANR EUREKA); SB, AMD and MH were responsible for SAG preparation and sequencing, writing -review and editing. All authors contributed to the manuscript and approved the final version. and cDNA samples (MetaB 16S cDNA) collected along the Lake Pavin water column. Sulfu1: HicA HicB Chromosomal location: l: 1 copy, : 2 copies, t: 3 copies, n: 4 copies, : 5 copies, v: 7 copies, : orphan toxin (1 copy).

Figure titles and legends