The microbiome of the Black Sea water column analyzed by shotgun and genome centric metagenomics

Background The Black Sea is the largest brackish water body in the world, although it is connected to the Mediterranean Sea and presents an upper water layer similar to some regions of the former, albeit with lower salinity and temperature. Despite its well-known hydrology and physicochemical features, this enormous water mass remains poorly studied at the microbial genomics level. Results We have sampled its different water masses and analyzed the microbiome by shotgun and genome-resolved metagenomics, generating a large number of metagenome-assembled genomes (MAGs) from them. We found various similarities with previously described Black Sea metagenomic datasets, that show remarkable stability in its microbiome. Our datasets are also comparable to other marine anoxic water columns like the Cariaco Basin. The oxic zone resembles to standard marine (e.g. Mediterranean) photic zones, with Cyanobacteria (Synechococcus but a conspicuously absent Prochlorococcus), and photoheterotrophs domination (largely again with marine relatives). The chemocline presents very different characteristics from the oxic surface with many examples of chemolithotrophic metabolism (Thioglobus) and facultatively anaerobic microbes. The euxinic anaerobic zone presents, as expected, features in common with the bottom of meromictic lakes with a massive dominance of sulfate reduction as energy-generating metabolism, a few (but detectable) methanogenesis marker genes, and a large number of “dark matter” streamlined genomes of largely unpredictable ecology. Conclusions The Black Sea oxic zone presents many similarities to the global ocean while the redoxcline and euxinic water masses have similarities to other similar aquatic environments of marine (Cariaco Basin or other Black Sea regions) or freshwater (meromictic monimolimnion strata) origin. The MAG collection represents very well the different types of metabolisms expected in this kind of environment. We are adding critical information about this unique and important ecosystem and its microbiome. Supplementary Information The online version contains supplementary material available at 10.1186/s40793-021-00374-1.


Background
The Black Sea is the innermost arm of the Mediterranean basin. Nearly severed from the rest by the tectonic movement of the African plate, it is only connected to the Mediterranean Sea by the narrow but deep strait of the Bosporus. The Black Sea has a positive hydric balance i.e. receives more freshwater than lost by evaporation and hence contains less salt (from 0.73% in epipelagic to 2.2% in meso-bathypelagic waters) than the Mediterranean proper (3.8%). In addition, the large watershed and riverine inputs lead to a richer nutrient status (meso-eutrophic) and permanent stratification with a colder, more saline deepwater mass that remains anaerobic and largely euxinic below 150-200 m [1][2][3]. All these properties make the Black Sea a unique brackish-marine environment. Its depth (average 1253 m with a maximum of 2212 m) makes this system much more stable than other brackish inland water bodies like the Baltic Sea in which the anaerobic compartment is only a recent development due to anthropogenic impacts [4].
Although a few studies have been carried out by metagenomics and metagenome-assembled genomes (MAGs) reconstruction, the information available in databases about this unique environment is scarce. A recent study showed for the first time the microbial structure of the sulfidic waters of 1000 m depth [5], mainly dominated by sulfate reducers (Desulfobacterota and Chloroflexi classes such as Dehalococcoidia or Anaerolinea), associated dissolved organic matter (DOM) degraders (Marinimicrobia, Cloacimonadota) and streamlined uncultured taxa such as Omnitrophota, Parcubacteriota or Woesearchaeota [6]. At the genome level, Black Sea microbes remain largely unknown. Only 10 MAGs were studied in the abovementioned work [5]. However, 179 MAGs from 50 to 2000 m have been recently deposited into Genbank (PRJNA649215) alongside two recent works that have described various Synechococcus phylotypes [7,8]. Most recently, members of widespread clades such as SUP05 (Ca. Thioglobus spp), Sulfurimonas bacteria, and uncultivated SAR324 and Marinimicrobia have been studied from this and other dysoxic environments [9]. Another habitat that could be comparable to the Black Sea in terms of water column stratification is the Cariaco Basin (Venezuela). Several studies on the redoxcline and euxinic strata on this system provided some of the first insights into microbial structuring in anoxic marine waters and Oxygen Minimum Zones (OMZs) [10,11]. Here, we present a genome resolved metagenomic study of different depths in the Black Sea adding a total of 359 high-quality MAGs. The epipelagic and DCM strata show an overall marine-brackish community composition with a predominance of microbes similar to the Mediterranean and the Caspian Seas [12,13]. The anaerobic compartment, which accounts for up to 80% of the total Black Sea volume, had a much more poorly understood microbiota including various members of the microbial "dark matter".

Black Sea shotgun metagenome analysis
Samples were collected along the Bulgarian coast at two stations (Fig. 1a). Sampling depths were chosen on the grounds of the physicochemical measurements (Fig. 1b) to cover representative temperature, oxygen, and chlorophyll-a values (Additional File 1). Thus, for St. 307, with a maximum depth of 1100 m, samples were from the near-surface at 5 m depth, the deep chlorophyll maximum (DCM) at 30 m, a sample from the redoxcline/pycnocline at 150 m, and finally, a sample at 750 m, corresponding to the euxinic water layer. Additionally, we collected a single near-surface sample (5 m) closer to the shore at station 301 with a maximum depth of 22.5 m to examine the coastal influence, but no significant difference was found with the surface offshore sample, and hence this sample will not be described any further although it was used to generate 87 MAGs.
We first performed an unassembled read analysis for each depth to obtain an overview of the taxonomic profile based on 16S rRNA gene fragments retrieved from the metagenomes that were compared to the SILVA database [14] (Fig. 1c) and the main predicted metabolic functions assessed by the SEED subsystems [15] (Fig.  1d). We have also made a distance-based redundancy analysis (dbRDA) to statistically assess the correlation among samples, environmental parameters, and the abundance of metabolic pathways (Fig. 2).
The oxic strata (surface and DCM) presented, at this large-grained level, a very similar taxonomic composition. We have used throughout the taxonomic nomenclature derived from the Genome Taxonomy Database (GTDB) [16]. In the photic zone (5 m and DCM), Alphaproteobacteria (orders SAR11, SAR116, Rhodobacterales, and Rhodospirillales), Gammaproteobacteria (mostly orders SAR86 and Pseudomonadales), and picocyanobacteria (order Synechococcales) were the most abundant groups, representing > 70% of the microbial community (assessed by 16S rRNA gene fragment abundance). It must be highlighted the complete absence of the genus Prochlorococcus that contrastingly predominates in the Eastern Mediterranean Sea [12]. This absence has been assessed through genome recruitment of all Prochlorococcus genomes available in the NCBI thus far, having obtained no significant mapped reads at > 70% identity and > 50 bp of alignment lengths. The predominant metabolic subsystems of the oxic layer ( Fig. 1d) were, as expected, associated with phototrophic lifestyles such as those from Synechococcales (photosystems/ phycobilisomes) or photoheterotrophy with type-1 rhodopsin pumps (generally present in SAR11, SAR86, or Flavobacteriales). Besides, ammonia and nitrite were the preferred N source, whilst DO, T, and pH were the main environmental variables correlating to oxic samples (Fig. 2).
The taxonomic composition derived from 16S rRNA gene fragments changed dramatically as we reached oxygen extinction in the redoxcline (150 m), where various taxa and microbial lifestyles coexisted, with a prevalence of anaerobic N and S related subsystems Each environmental measurement is color-coded. c Black Sea phylum level 16S rRNA gene taxonomic classification assessed with unassembled reads individually mapped into SILVA138 database. d Metabolic profiles assessed with abundance values obtained from reads and SEED subsystems and pathways, pecific proteins/genes or microbial systems. A blue row Z-score scale was applied to statistically assess abundance values differing between samples ( Fig. 1d and Fig. 2). Marinimicrobia (ca. 30% of 16S rRNA gene assigned reads) and Gammaproteobacteria (ca. 20%) were the dominant taxa (Fig. 1c). Also, significant percentages of various associated streamlined microbes such as Patescibacteria and Nanoarchaeota (ca. 4% each) were found in this layer. The redoxcline of the Black Sea presented the most metabolically diverse set of pathways among all analyzed samples ( Fig. 1d and Fig. 2). The main environmental variables that were statistically grouped with the redoxcline were total nitrogen (TN) and nitrate, which were associated with the different N biogeochemical cycle pathways. The highest abundance pathways were associated with denitrification (nitrogen gas as the final product), nitrate/nitrite ammonification, and dissimilatory nitrate reduction (with ammonium as the final product), but the N cycle was also completed with ammonia oxidation and N fixation pathways detected both in total reads (Fig. 2) and MAGs (see below). Nonetheless, various other metabolisms coexisted in this thin layer where oxygen disappears. Such pathways, assigned from global read abundance included fermentation and hydrogen metabolism, chemolithotrophic, anoxygenic photosynthesis, or sulfate reduction.
The euxinic waters at 750 m had an expected increase in sulfate reducers (Desulfobacterota), Dehalococcoidia/ Anaerolineae Chloroflexota. There were representatives from Omnitrophota and Kiritimatiellae (both classified inside Verrucomicrobiota according to SILVA standards [14], although Omnitrophota is classified as a single phylum according to GTDB [16]), Phycisphaerae Planctomycetota, Marinimicrobia, Nanoarchaeota, Patescibacteria, and Cloacimonadota. Finally, Halobacterota (Syntrophoarchaeia) and Crenarchaeota (Bathyarchaeia) minor representation (< 2% of total microbial biomass assessed by 16S rRNA gene) showed that methanogenesis coexisted with sulfate reduction in these euxinic waters, if in a much more reduced fraction. As for the SEED subsystems, the euxinic waters (750 m) showed an increase in fermentation, hydrogen metabolism, anaerobic respiratory reductases, or methanogenesis pathways ( Fig. 1d and Fig. 2). Overall, we observed a huge diversity of microbiota providing hydrogen and fermentation by-products that conformed a syntrophic network fueling the sulfate reducers at the redox end.

Comparison of the Black Sea and Cariaco Basin metagenomic datasets
To assess the representativity of our samples we compared the community structure by 16S rRNA gene (following the same abovementioned methodology) and the read similarity (clustering) among various similar datasets, including other Black Sea metagenomes available (PRJNA649215) [9] (see Fig. S1) and with a former sampling campaign from another marine euxinic system, the Cariaco Basin [10,11] (see Fig. S2).
As depicted from Fig. S1, we have observed both a clear 16S rRNA gene taxonomic profile (Fig. S1A) and read clusterization (Fig. S1B) between Black Sea samples obtained from different sampling campaigns. It is remarkable the stability observed in this system, even though samplings were taken 6 years apart, and in different locations. All the abovementioned taxa showed up in the compared samples as well, observing a clear sample clusterization between our DCM sample at 30 m and a 50 m former sample (PRJNA649215). Our redoxcline (150 m) sample appears in between 130 and 250 m samples from PRJNA649215, whilst the mesopelagic euxinic 750 m sample that we took clusters with the group of former meso/bathypelagic strata (500, 1000, and 2000 m).
Regarding the Cariaco Basin (Venezuela), we compared exclusively redoxcline metagenomes from this marine environment with our two redoxcline and anoxic datasets, observing a markedly different taxonomic profile (Fig. S2A) and read clusterization (Fig. S2B). As a general trend, Gammaproteobacteria, SAR324, Actinobacteriota, and Thermoplasmatota showed a higher % in the Cariaco Basin. The most similar samples to our (150 m) in terms of taxonomy were the 247 and 314 m.

Recovery of MAGs
Next, we have reconstructed metagenome-assembled genomes (MAGs) from each dataset separately. Automated binning followed by manual curation of the bins (see methods) generated 359 MAGs with > 50% completeness and < 5% contamination. Statistics of these MAGs are described in Table 1 (MAGs from oxic samples) and Table 2 (redoxcline/anoxic MAGs) and in individual detail in Additional File 2. We also plotted the  . Among all of the MAGs retrieved from this work, we selected the 50 most abundant MAGs (> 10 RPKGs in any of the recruited samples) from the oxic, redoxcline, and anoxic waters and recruited them at > 95% of identity (species level) on all metagenomes and showed their main distinguishing metabolic features (Fig. 4). The rest of the MAGs abundance in all datasets is shown in detail in Additional File 3.

MAGs from the epipelagic and DCM oxic strata
As expected, emblematic key players of the oxic waters harboring a photoheterotrophic/phototrophic lifestyle such as Ca. Pelagibacter, Ca. Actinomarina, SAR86, Synechococcus, or Flavobacteriaceae were detected in high numbers in all oxic metagenomic datasets. As noted above, we must highlight the complete absence of Prochlorococcus spp. genomic footprints (with no significant genome recruitment above 70% of identity with any existing Prochlorococcus genome), contrasting with a high abundance of Synechococcus MAGs that affiliated with the marine clades I, III, IV, VI and WPC1 including Fig. 3 Black Sea Distance-based redundancy analysis (dbRDA) between different samples (depths), environmental parameters and MAGs abundance values obtained from the RPKGs.. We provided MAG names for class (_c), order (_o), family (_f), genus (_g) or species (_s) according to the last GTDB classification scheme followed in Table 1 and Table 2. Each Black Sea sample is red square-coded. Environmental parameters are blue triangle-coded and MAGs are green triangle-coded  [18]. The main Actinobacteria MAGs retrieved presented relatively small genome sizes (1.2-2.2 Mb), among which we must highlight the presence of 5 novel Actinomarinales spp. (BS301-5 m-G7, BS307-5 m-G2, BS30m-G2/G3/G4) and a group of Ilumatobacteraceae genomes related to Caspian MAGs (Casp-actino5) [13]. The major SAR11 Alphaproteobacterial MAGs were eight novel Pelagibacterales that affiliated with the recently described groups Ia.1, IIaB/1and IIIa [19]. Remarkably, we obtained three novel MAGs from the order Rickettsiales, that showed similarities with groundwater aquifer MAGs (< 80% of ANI) and freshwaters from Lake Baikal [17]. Another relevant Alphaproteobacteria clade from which we obtained MAGs was SAR116, with six MAGs affiliated to the Puniceispirillum genus and five more were only classified as representatives of the family Puniceispirillaceae. A family that provided abundant MAGs in the Black Sea oxic waters is Flavobacteriaceae (23 MAGs), a group that was commonly detected in the Mediterranean [12] and Baltic Seas [20]. Various MAGs were related at GTDB genus level with MED-G11, MED-G14 MAGs, and at the species level (ANI > 95%) with MED-G20 Mediterranean Sea MAGs. Two MAGs also showed their closest relatives at the GTDB family level with Baltic Sea MAGs BACL11 and at the species level with BACL21. We also found five methylotrophic representatives from the clade OM43 (family Methylophilaceae) affiliating at the genus level to BACL14 Baltic Sea MAGs [20]. Eleven MAGs belonged to the cosmopolitan Gammaproteobacteria SAR86, so far only classified at this order level. Other Gammaproteobacteria that co-occurred in these samples were MAGs with similarity to Luminiphilus (11 MAGs) [21] and Litoricola (3 MAGs) [22] genera. Another relevant taxon from marine systems was the former marine group-II Euryarchaeota (Thermoplasmatota according to GTDB taxonomy). We retrieved six genomes affiliating to the family Poseidoniaceae and the other six to the genus Poseidonia. Only one genome was obtained affiliating to Thalassoarchaeaceae. Finally, three ultra-small (1 Mb of estimated genome size) Marinimicrobia MAGs were obtained from the oxic metagenomes, which so far are classified by the GTDB as genus Marinisoma.

MAGS from the redoxcline
Next, we showed the main ecological drivers of the redoxcline layer, which included chemolithotrophic S oxidizers and C fixers such as Ca. Thioglobus and dissimilatory nitrate reducers such as Sulfurimonas, both of which were recently described as members of the Black Sea microbiome [9]. The redoxcline also showed some other ecologically relevant nitrate reducers such as Bacteroidales MAGs, ammonia oxidizers such as Nitrosopumilus spp., sulfate reducers Desulfatiglandales and Desulfococcales novel species, and several Marinimicrobia representatives which are specialized in the H metabolism, denitrification, and mixed-acid fermentation. We noticed the presence of anoxygenic photosynthesis, exemplified by MAG BS150m-G13 showing > 99% of ANI with Chlorobium phaeobacteroides, a green sulfur bacterium (GSB) originally isolated from the Black Sea [23] (GCA_000020545.1), that was undergoing a nearly monoclonal bloom (Fig. S4). Chemoautotrophy was detected in Thioglobus sp. BS150m-G29 and G33 MAGs, both of which are novel representatives of the SUP05 clade which perform a wide variety of metabolisms including S oxidation and C fixation, and with only 80% of ANI with its closest relative (Ca. Thioglobus autotrophicus EF1) [24]. It appears that this is a case of a single species (recruiting at > 95% of nucleotide identity) abundant (> 70 RPKG, Fig. S5) in the Black Sea redoxcline. Methane oxidation (Methylobacter sp. BS150m-G31) and ammonia oxidation were also metabolisms observed in MAGs from this layer (Nitrosopumilus spp. BS150m-G38/39/40). Nitrite oxidation was detected in Nitrospinaceae BS150m-G45.
We also assessed the presence of our redoxcline MAGs with previously available metagenomes from the redoxcline from the Cariaco Basin (Venezuela) [10]. Overall, it seems that Marinisomatota/Marinimicrobia and Gammaproteobacteria chemolithotrophic groups are the most abundant key players of these two marine redoxclines, accounting for more than 50% of total microbial biomass (Fig. S2A). However, it must be noted that only a few MAGs from the Black Sea were detected there (Fig. S2C). Among them, two chemolithotrophic Gammaproteobacterial representatives (Ca. Thioglobus and a novel species BS150m-G30 classified only at the order level as o__GCA-2,400,775 by GTDB), sulfate reducers (Desulfatiglandales), denitrifying and hydrogenproducing Marinimicrobia (three species) and one Actinobacteriota (a novel species from the marine MedAcidi-G1 group). Apart from their metabolic potential fitting with microbial lifestyles from redoxcline layers, these species could play key roles in other marine redoxclines and oxygen minimum zones (OMZs), as their detection in two largely separated biomes with different salinities (ca. 2% in the Black Sea and 3.5% in the Cariaco Basin) indicate a widespread distribution in oxygen-depleted marine niches.

Euxinic MAGs and the "microbial dark matter"
The main environmental variables grouping with the mesopelagic sample of 750 m were PO 4 and Si, both of which are solubilized in anoxic layers and diffuse from the sediment layer (Figs. 2 and 3). Salinity also increased up to 2.2% in these euxinic waters. There is the expected predominance of sulfate reduction pathways, as carried out by Desulfobacterota MAG representatives (Desulfatiglandales BS750m-G47-G51 and BS750m-G54/G56, Desulfobacterales BS750m-G52/G53/G55), which perform the dissimilatory sulfate reduction pathway (dsr genes). MAGs that were detected at all euxinic strata included various other sulfate-reducers, their associated microbiota performing mixed-acid fermentations and H metabolism in syntrophism (Cloacimonadetes, Woesearchaeales, Ca. Aminicenantes) and the only MAG able to perform methanogenesis and ANME, a Syntrophoarchaeum that showed a remarkable abundance in 1000 and 2000 m metagenomic datasets, suggesting that methane metabolism indeed coexists with sulfate reduction and associated microbial fermenters and H 2 scavengers in a complex syntrophic network.
The key methanogenesis genes (mcr and fwd) were detectable in the abovementioned single MAG, Ca. Syntrophoarchaeum BS750m-G82. However, other genomes classified as Bathyarchaeota BS750m-G27/G28 showed most of the methanogenesis genes except for the key enzyme, Methyl-coenzyme M reductase (mcr). These enigmatic Bathyarchaeota representatives showed various mixed-acid fermentation pathways including the formation of H 2 and CO 2 (via the formate hydrogenlyase), formate (pyruvate-formate lyase), alcohol (alcohol dehydrogenase), or lactate (lactate dehydrogenase). It must be noted the potential capability of performing reverse methanogenesis, or anaerobic methane oxidation (ANME) by the abovementioned archaeon (MAG Ca. Syntrophoarchaeum BS750-G82). Heterodisulfide reductase genes (HdrABC), which are involved in electron transfer and the last step of methanogenesis by reducing CoB-CoM heterodisulfide, were detected in Bathyarchaeota and Syntrophoarchaeum MAGs. However, these genes were also found in Cloacimonadota, candidate division KSB1, Ca. Aminicenantes, Omnitrophota, Desulfobacterota, Planctomycetes, and Chloroflexi MAGs as well as in the unassembled reads (being completely absent from oxic datasets), suggesting that these electron transfer complexes are not exclusive to methanogens. As seen by the dbRDA, we also noted a global predominance of mixed-acid fermentation pathways (with ethanol, lactate, acetate, formate, or CO 2 /H 2 as products) and hydrogen uptake hydrogenases that couple with sulfate, fumarate, CO 2, or nitrate reduction, thus conforming a complex syntrophic network of microbes. This networking of syntrophic microbes (considered here as interspecies H transfer) includes the abovementioned uncultured taxa plus accompanying streamlined members of the "microbial dark matter" such as Omnitrophota, Patescibacteria (Ca. Microgenomates, Portnoybacteria, Paceibacteria) or Nanoarchaeota (Ca. Aenigmarchaeota, Woesearchaeota, Pacearchaeota), groups from which we also obtained MAGs (see Table 2). Various types of hydrogenases and hydrogen metabolism pathways grouped with the 750 m mesopelagic sample in the dbRDA plot (Fig. 2) and were found in the vast majority of microbes inhabiting this sulfide enriched waters, including NADreducing bidirectional (hox genes) and uptake hydrogenases (hup genes), NiFe (hyp genes) and FeFe (hym genes) hydrogenases, Coenzyme F420-reducing hydrogenases or carbon monoxide induced hydrogenases (CooHL genes), all of which showed the highest gene/ recA ratios (from 0.2 in hym genes to 1-2 for hyp and hoxF) in euxinic waters.
It was remarkable the presence of two Actinobacteriota MAGs (BS750m-G1/G2) in these sulfide-rich waters. These yet unclassified members have their highest resemblance with MAGs retrieved from groundwater aquifers (Actinobacteriota bacterium CG08_land_8_20_ 14_0_20_35_9, classified as UBA1414 by GTDB) and have very low GC content (31-34%) and predicted genome sizes (ca. 1.4-1.6 Mb). Their genomes presented various mixed-acid fermentative pathways associated with the production of ethanol (alcohol dehydrogenase), lactate (lactate dehydrogenase), formate (pyruvate-formate lyase), and H/CO 2 (formate hydrogen lyase). They also showed an active hydrogen metabolism with various NiFe hydrogenases including Coenzyme F420-reducing hydrogenase, hyp genes, and HyaA COG1740 355 Ni-Fe-hydrogenase I. Another remarkable group of microbes was Omnitrophota, from which we obtained 15 MAGs with variable estimated genomes sizes (from 1 to 3 Mb). For instance, the most abundant Omnitrophota MAG retrieved from our samples (BS750m-G77) presented a small predicted genome size (ca. 1.2 Mb) and was an obligate fermenter (mainly producing ethanol, H 2 /CO 2, and lactate). Another group of streamlined members of the microbial dark matter was Aenigmarchaeota (BS750m-G24/36/81/83/) and Nanoarchaeota (BS750m-G11/13/70) MAGs, which had estimated genome sizes of 1-1.5 Mb. Among their metabolic potentials, they were also mixed-acid fermenters, including lactate or H 2 /CO 2 as fermentation by-products, which would fuel the sulfate reducers, conforming a syntrophic network with the rest of the mixed-acid fermenters. Finally, another set of microbes of small genome sizes (0.6-1.6 Mb) were Patescibacteria (former Candidate Phyla Radiation). We must highlight the presence of the protein VirB4, associated with type IV secretion systems that work as injectors into host cells [26], in Ca. Microgenomates BS750m-G73/74, Ca. Paceibacteria BS750m-G71/75 and Ca. Portnoybacteria bacterium BS750m-G76. These proteins were only found in these microbes in the entire euxinic waters, which suggests a parasitic lifestyle in which these Patescibacteria could translocate nutrients, proteins, and DNA from or to a putative host [26].

Discussion
The first ecogenomic insights of Black Sea water column microbiota showed, at first sight, that this brackish microbiome has various taxa previously observed in other brackish systems such as Caspian [13] or Baltic Seas [20]. However, its interconnection with the Mediterranean Sea has allowed the coexistence of some other microbes not found in the other two [12]. However, the uniqueness of this enormous brackish water mass lies in its redoxcline and very deep euxinic strata, features that are not present in the Caspian or Mediterranean Seas.
The present study, together with others from the Black Sea [5,9] (bioproject PRJNA649215) and the Cariaco Basin in Venezuela [10,11] provide a first glimpse into the microbiome of anoxic marine water columns and Oxygen Minimum Zones (OMZs). The redoxclines of these habitats show a convergence of various coexisting metabolisms, among which we encounter anoxygenic photosynthesis such as that one observed in Chlorobium phaeobacteroides [27], ammonia oxidation by Nitrosopumilus spp., chemolithotrophic metabolisms carried by Gammaproteobacteria such Ca. Thioglobus, one of the most abundant and versatile players transitioning between oxic-anoxic regimes [24,28]. This microbe has been detected both in Cariaco and Black Sea basins [9] and its adaptive metabolism, which includes physiological adaptations to oxic-anoxic growth [28] and its broad metabolic potential has led it to colonize these zones with a large contribution to S and N biogeochemical cycles as a denitrifier, sulfur-oxidizing, and C fixer chemolithotroph. The simultaneous activity of sulfatereducing and sulfide-oxidizing microbes in these habitats has led to a term known as the "cryptic sulfur cycle" [25]. One of the most abundant microbes detected from these zones belongs to the Marinimicrobia, which are specialists of both oxic and anoxic waters [29,30], able to perform denitrification and various fermentations and H 2 metabolism and were recently labeled as organoheterotrophs with specific molybdoenzymes to generate energy from sulfur cycle intermediates [9].
As we approach the euxinic meso and bathypelagic waters, we encountered domination of sulfate reduction coupled with a complex variety of syntrophic networks that feed the ecosystem. In this sense, Desulfobacterota is a good example of a syntrophic phylum that could be able to accept electrons from organic electron donors, as noted previously in marine sediments [31]. It appears that the extremely high abundance of sulfate-reducers in the Black Sea has displaced methanogens, which are present in the water column but low in numbers. We found MAGs of Bathyarchaeota [32,33] representatives and a single example of a Syntrophoarchaeum in this study. This microbe could be performing reverse methanogenesis or anaerobic methane oxidation (ANME) rather than the standard methanogenesis. Recently, members of this newly identified species performed the complete oxidation of butane during the anaerobic methane oxidation process [34]. However, it appears that the competition between methanogens and sulfatereducers for acetate is dominated by the latter, which also take a fundamental role in the syntrophic networks by utilizing the H 2 produced by fermentation. Among all of the associated microbiota, we must highlight the presence of the Cloacimonadota phylum (previously known as WWE1), which have been detected in meromictic lakes as important carbon and sulfur recyclers [35] and appear to degrade propionic acid in syntrophic networks in bioreactors [36]. Novel microbes from lineages such as Cloacimonadota (WWE1), Marinisomatota (SAR406), Omnitrophicaeota (OP3), Bacteroidales, Kiritimatiellae, Anaerolinea/Dehalococcoidia formed a very complex syntrophic network where mainly mixed-acid fermentations with lactate, ethanol, formate, succinate, hydrogen, and CO 2 are the final products. Other microbes obtained as MAGs included members of the uncultured microbial dark matter, such as Patescibacteria and Nanoarchaeota [6,[37][38][39]. In particular, we have stumbled upon various members of Aenigmarchaeota [40] and Woesearchaeales with their usual streamlined genomes and hard to predict metabolic properties.
The study of microbes from this and other brackish water bodies could be useful in understanding better the reasons for the great divide between fresh and marine water bodies and shed light upon the physiological and genomic strategies that are responsible for marinefreshwater microbial transitions, considered infrequent [41] and demanding in terms of physiology and protein/ gene variation [42], that have been detected in several groups of aquatic microbes.

Methods
Sampling, DNA extraction, physical and chemical profiles measurement  29.001283 E, respectively). Up to 6.9 L of seawater from each sampling depth were filtered through a series of 20 μm Nylon Net filters (Millipore), 5 μm polycarbonate membrane filters (Millipore), and 0.22 μm SterivexTM Filter Units (Merck). DNA was then extracted only from these Sterivex filters using standard phenol-chloroform protocol [43]. In short, Sterivex filters were treated with CTAB lysis buffer and then treated with 1 mg ml − 1 lysozyme and 0.2 mg ml − 1 proteinase K (final concentrations). Then nucleic acids were purified with phenol/chloroform/isoamyl alcohol and chloroform/isoamyl alcohol.

Sequencing, assembly, and read annotation
The five samples were sequenced in one lane of Illumina HiSeq X 10 PE 2X150 bp (Novogene company), which provided ca. 180-200 million clean reads and 24 Gb of output for each sample. Samples were individually assembled using IDBA-UD [44] with the parameters --pre_correction, −-mink 50, −-maxk 140, −-step 10.Sub-assemblies of 20 million reads were done in each sample to retrieve some of the most abundant microbes that assembled poorly i.e. we reduced the total number of reads to obtain more complete bins of these representatives (e.g. Pelagibacterales, Actinomarinales, Synechococcus, Marinimicrobia or Thioglobus spp.). Annotation of contigs was assessed using Prodigal [45] for ORF prediction and then BLAST (nr database) using Diamond for functional annotation [46]. Proteins were annotated with customized nr, COG [47], and TIGFRAM [48] databases (last releases from January 2020) that served as a taxonomic proxy to understand what we obtained from our datasets. Features like tRNAs and rRNAs were detected with tRNAscan [49] and ssu-align [50], respectively.

Binning, classification, and MAG statistics
The Binning procedure was performed as follows: a first manual inspection was done assigning a hit (based on BLAST against Nr) to each CDS (coding sequence), which allowed us to classify contigs taxonomically into different phyla. This step, although not necessary for current binning approaches, has provided us good results in past studies of emblematic aquatic systems such as the Mediterranean Sea or the Lake Baikal [12,51]. Then, an initial binning step was applied for each set of contigs assigned to each phylum with METABAT2 [52] using coverage in the different samples. When META-BAT2 provided a single bin containing various species, we applied a manual inspection of the annotated contigs and we used other computational approaches based on, GC content, RPKGs abundance values, and tetranucleotide frequencies to refine and separate the bins [53,54]. Finally we only used MAGs with < 5% contamination and > 50% of completeness based on CheckM estimations [55]. MAGs were taxonomically classified according to the latest version of GTDB-tk and the database release89 [16] and whenever we could we used class, order, family, genus, or species names for all of them.
16S rRNA gene read classification, hierarchical cluster read analysis, read functionality The 16S rRNA gene reads were detected in a subset of 20 unassembled million reads from each metagenome. We first obtained candidates using USEARCH [56] with RefSeq 16S rRNA gene as database and then these putative 16S rRNA genes were confirmed using ssu-align [50]. Then, a BLASTN was performed against the SILVA database [14] (SILVA_138_SSURef_Nr99_Tax_silva from December 2019) to provide a taxonomic classification. Hierarchical cluster analysis (dendrograms) of different metagenomic samples with k-mer = 21 bp was assessed with SIMKA [57] and Bray-Curtis indexes of presence/ absence were obtained. Subsets of 20 million reads of each metagenomes were analyzed with BLASTX against SEED [15] database using Diamond [46], with parameters more-sensitive, max-target-seqs 1, e-value 0.00001 > 50 bp of alignment length and > 50% identity. The top hits were analyzed in search of specific genes and pathways based on the SEED database. Hits were normalized by the total number of hit counts for each sample and a row Z-score was calculated to assess statistical differences between samples for each metabolic pathway.

Relative abundance of MAGs
To estimate the relative abundance of the recovered genomes in various datasets we performed read recruitment and mapping, which was assessed considering BLASTN hits of the metagenomic reads against each MAG at > 95% identity and 50 bp of alignment length thresholds, as indication of belonging to the same species. A microbe was considered present in a metagenome if it was detected at > 1 RPKG (Reads per Kb of Genome per Gb of Metagenome). All relative abundances of our MAGs on Black Sea datasets are shown in Additional File 3. Datasets used for recruitment included Black Sea (PRJNA649215) and Cariaco Basin (PRJNA326482).

Redundancy analysis (RDA) of environmental variables, MAGs, and metabolic processes
Distance-based redundancy analysis (dbRDA) analysis was performed to describe the ordinations of the main MAGs and metabolic processes in an environmentally constrained space [58] and conducted with the R package vegan [59]. Environmental matrixes were constructed with 12 environmental variables for the 5 Black Sea samples. Each matrix was square-root transformed and normalized and subsequently transformed to Euclidean resemblance matrix. On the other hand, we constructed the other two matrixes with the recovered MAGs and metabolic processes with 30 selected metabolic processes. Both were standardized and square root transformed before performing a Bray-Curtis dissimilarity resemblance matrix. Two dbRDA were obtained, the first one comprising the recovered genomes (MAGs) matrix using environmental matrix as predictor variable, and a second one based on the metabolic processes matrix using environmental matrix as predictor variable.

Availability of data and materials
All metagenomes and reconstructed genomes derived from this work are publicly available under the NCBI Bioproject PRJNA638805. Metagenomes were deposited to NCBI-SRA with the accession numbers SRR12042682-SRR12042686.