The rhizosphere of the seagrass Zostera marina was sampled at Culatra Island (ZmPt-Faro, Portugal, 36°59′56.0″ N 7°49′31.7″ W) and Pointe de Cléguer (ZmFr-Roscoff, France, 48°43′37.1″ N 3°58′35.9″ W), in 2013 as described and used in Cúcio et al. (2016). Briefly, sampling consisted of randomly collecting cores (15 cm diameter) of Z. marina (n = 5), from which 4–6 shoots were carefully separated and shaken to remove loose sediment. The rhizosphere was retrieved by washing the roots in 0.2 μm-filtered seawater, and was subsequently transported to the laboratory in cool boxes. Upon arrival in the laboratory, meiofauna and plant detritus were removed whenever present, and samples were treated in a homogenizer (Stomacher 80 Laboratory Blender, Seward Medical) to loosen the bacteria from the sediment particles. After 3 homogenizing cycles of 1 min at normal speed, the supernatant was centrifuged for 30 min at 10,000 g (Costa et al., 2006), and the resulting pellet was used for DNA extractions.

Extraction of genomic DNA from the rhizosphere bacteria was performed using the PowerSoil DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA), according to manufacturer's instructions. Part of the DNA resulting from this extraction was used in our previous publication as 5 individual replicates per location (Cúcio et al., 2016), and due to the highly comparable community composition among the samples within location (Cúcio et al., 2016), all replicates of each location were pooled to create two rhizosphere samples for shotgun metagenomic analysis, one from Portugal (Pt) and another from France (Fr). The quality of the genomic DNA was determined by agarose electrophoresis, and the quantification was performed with a dsDNA HS Assay Kit on a Qubit 2.0 Fluorometer.

The genomic DNA of the two samples was sequenced by the Beijing Genomics Institute (BGI, Hong Kong, China) in one lane of paired-end sequences with an insert size of 170 bp, using an Illumina HiSeq 2000 PE 100 sequencer.

The quality of the metagenomic reads was checked with FastQC (Babraham Bioinformatics, http://www.bioinformatics.babraham.ac.uk/projects/fastqc), and trimmed using CLC Genomics Workbench (v8.5.1, CLCbio) with a quality score limit of 0.05, and a maximum allowance of 2 ambiguous nucleotides.

Taxonomic assignment of metagenomic reads was carried out using Kaiju (http://kaiju.binf.ku.dk), using the non-redundant protein database NCBI BLAST nr as reference database (Menzel et al., 2016). “Greedy” run mode was applied with a minimum match score of 70, and an allowance of five mismatches (Menzel et al., 2016). High-quality metagenomic paired-end reads were assembled using IDBA-UD (v1.1.1, Peng et al., 2012), with pre-correction for improved efficiency of highly uneven sequencing depth, and k-mer sizes from 40 to 100. To obtain more insights into the eelgrass rhizobiome and to increase the amount of information for coverage and binning, both metagenomes were further co-assembled using MEGAHIT (v1.0.3-6, Li et al., 2015), using a minimum k-mer size of 25 incremented in steps of 10 on each iteration. Hereinafter, results obtained from the individually assembled metagenomes are referred to by the name of their sampling location (i.e., “Portugal” or “ZmPt”, and “France” or “ZmFr”), and results obtained from the co-assembled metagenomes are referred to as “co-assembly”.

Protein-coding genes were predicted with Prodigal (Hyatt et al., 2010), and tRNAs with tRNAscan-SE (Lowe and Eddy, 1997). Annotation of protein-coding genes was performed with BLAST against COGs (Tatusov et al., 2003), TIGRfams (Haft et al., 2001), and NCBI's non-redundant database. In parallel, annotation of functional genes was also performed using the KEGG (Kyoto Encyclopedia of Genes and Genomes) Orthologs (KO) database in the KEGG's online annotation tool, GhostKOALA (Kanehisa et al., 2016). For those genes that can participate in both oxidative and reductive pathways of the sulfur cycle (Table 1), the function was confirmed by phylogenetic placement in reference trees, as well as according to the percentage of similarity to the closest hit attributed by BLAST during annotation, using a minimum identity of 75%.

The presence of specific eukaryotic sequences was investigated based on raw reads, using the MG-RAST server (v4, Meyer et al., 2008).

Databases of the marker genes for dissimilatory sulfur processes (“sulfur” genes) adenylylsulfate reductase (aprA) and sulfate adenylyltransferase (sat) were built using reference alignments of the ortholog clusters ENOG4105CFK and ENOG4107RIQ, respectively, retrieved from the eggNOG phylogenomic database (v4.5, Huerta-Cepas et al., 2016). A database of sulfide:quinone oxidoreductase/flavocytochrome c (sqr/fcc) orthologs was generated with sequences used in previous work performed by Marcia et al. (2010) and Han and Perner (2015). To complement the databases with genetic information about symbiotic bacteria, amino acid sequences of the proteins aprA, sat, sqr and fccB from symbionts were obtained from the Nucleotide database of NCBI (https://www.ncbi.nlm.nih.gov) by searching the terms “<gene of interest> AND symbiont.” A reference alignment for these proteins was performed on the online multiple sequence alignment tool Clustal Omega (Sievers et al., 2011). To build the phylogenetic tree for the sulfur gene dissimilatory sulfite reductase dsrAB, we used the reference database published by Müller et al. (2015).

Contigs annotated as orthologs of each gene were manually aligned to the reference alignments of their respective database in the ARB alignment tool (Ludwig et al., 2004), and phylogenetic trees were calculated using an approximately maximum-likelihood method in FastTree v2.1.3 (Price et al., 2010). Contigs that did not align to the reference databases used were not included in this analysis.

Assembled contigs (>1 kb) from ZmPt, ZmFr and co-assembly, were binned according to their pentanucleotide frequency using VizBin (Laczny et al., 2015). Taxonomic affiliation, completeness, contamination, strain heterogeneity and coverage of draft genomes were determined with the software CheckM (v1.0.5, Parks et al., 2015). Furthermore, assembled genomes were manually inspected to verify that the operons/gene cassettes are consistent in their structure. We retrieved about 20 draft genomes, but only selected those that had a completeness of > 70% and a contamination of < 5%. To identify the potential for CO₂ fixation, we searched for the carbon fixation marker genes cbbL, cbbS and cbbM from the Calvin-Benson-Bassham cycle (CBB), and korAB, frdA, porAB, and aclAB from the reverse tricarboxylic acid cycle (rTCA). The presence of nitrate-reducing and anammox bacteria was investigated through the presence/absence of genes involved in denitrification (narG, napAB, nirS, nirK, norC, and nosZ), dissimilatory nitrate reduction to ammonium (DNRA; nrfA), and anaerobic ammonium oxidation (Anammox; hzo, and hao). Nitrogenase-encoding genes nifH, nifD and nifK were also searched for, in order to identify nitrogen-fixing bacteria. Furthermore, sulfur metabolism was assessed by investigating the presence of the sulfur genes mentioned in the previous section, and summarized in Table 1. The search for these genes was performed based on their KO identifiers and gene/protein names.

Based on the most abundant bacterial classes identified among classified reads, the community structure, as revealed by taxonomic assignment of raw sequence reads, was highly similar between Portugal and France. The percentage of unclassified reads was similar in both metagenomes (37.36% in ZmPt and 38.57% in ZmFr). Nearly 98% of the classified reads from both metagenomes could be assigned to Bacteria, whereas < 2% of the reads were assigned to Archaea, and 0.1–0.2% of the reads had a viral origin. Based on classified reads, the bacterial communities were dominated by Gammaproteobacteria (25.3% for ZmPt, and 21.6% for ZmFr), followed by Deltaproteobacteria (20.0% for ZmPt, and 16.3% for ZmFr), and Alphaproteobacteria (7.9% for ZmPt, and 12.2% for ZmFr) (Figure 1, Supplementary Table 1). The most abundant classes were comparable between sites. However, a larger difference was found in the number of reads from the less abundant classes, such as Flavobacteriia, Anaerolineae, Gemmatimonadetes, Epsilonproteobacteria, Caldinae, Acidimicrobiia, Dehalococcoidia, Ardenticatenia, Mollicutes, and Thermoflexia (Supplementary Table 1). The order Desulfobacterales equally dominated the Deltaproteobacteria in both sites, whereas Chromatiales, the most dominant order of the Gammaproteobacteria in the metagenomes, harbored nearly two-fold more reads in Portugal than in France (Figure 1, Supplementary Table 2).

To provide genome and ecosystem level insight in the general rhizosphere of Z. marina we co-assembled the metagenomes from Portugal and France. By mapping the reads from both metagenomes to the co-assembly, we identified a higher coverage in Portugal than in France (Figure 2). Furthermore, although most reads were not shared between sites, three clusters of the co-assembly formed a clear exception, as their coverage was high for both metagenomes (clusters 1, 2, and 3, highlighted in Figure 2). Assembly of the reads of each cluster followed by phylogenetic placement with CheckM allowed the identification of cluster 1 as a member of the family Desulfobacteraceae within the Deltaproteobacteria and cluster 2 and 3 as Gammaproteobacteria.

A total of 263 high-quality sulfur-related gene copies were retrieved from the rhizosphere of ZmPt, from which 90% were assigned to the oxidative sulfur cycle with an average total coverage of 125x, in contrast to reductive genes, which only had an average total coverage of 35x (Figure 3).

The most abundant oxidation pathway in ZmPt was sulfide oxidation via sqr with a total coverage of 187x, followed by sulfite oxidation via aprAB, and thiosulfate oxidation via the sox complex (Figure 3). On the other hand, aprAB dominated the pool of reductive genes, with a total coverage of 70x (Figure 3). Only 65 sulfur-related gene copies were identified in the rhizosphere of ZmFr (Figure 3), of which 77% were assigned to genes involved in oxidative pathways. Nevertheless, reductive sulfur pathways dominated ZmFr in terms of coverage (average coverage of 20x vs. 24x for oxidative and reductive sulfur genes, respectively, Figure 3).

Phylogenetic analysis allowed the distinction of dsrAB (Figure 4), aprAB (Supplementary Figure 1) and sat (Supplementary Figure 2) genes from oxidative and reductive microorganisms. Contig sequences aligning to sequences in the reverse operating dsrAB database were mostly affiliated to the Gammaproteobacteria, whereas those aligning to the direct operating dsrAB were affiliated to the Deltaproteobacteria. Several contigs were grouping with symbiotic bacteria, such as a cluster of five contigs closely related to the endosymbiont of the polychaete Oligobrachia haakonmosbiensis (contigs FR-20896a, PT-3377ab, PT-61199b, PT-1051a, and PT-10561ab). Eight other contigs were directly related to endosymbionts of the oligochaete worms Olavius algarvensis and O. ilvae (contigs PT-42780ab, PT-47632ab, PT-7484ab, PT-55619ab as Gammaproteobacteria, and contigs FR-942b, FR-17686a, FR-8850b, and PT-27360ab as Deltaproteobacteria). Moreover, three other contigs were related to Candidatus Thiobios zoothamnicoli (contig PT-954ab, PT-562ab, and PT-16897a), an ectosymbiont of the ciliate Zoothamnium niveum (Figure 4).

As the taxonomic resolution obtained from the dsrAB tree is not extremely high due to the large number of uncultured bacteria to which our sequences aligned, sister lineages containing species such as Thioalkalivibrio thiocyanoxidans, Thiothrix nivea, and Allochromatium vinosum, allowed the placement of our sequences in the orders Chromatiales and Thiotrichales (Figure 4). On the reductive side of the phylogeny, however, a more detailed identification was possible for members of the families Desulfovibrionaceae, Desulfobulbaceae and Desulfobacteraceae (Figure 4).

Metagenomic binning resulted in 4 draft genomes isolated from ZmPt, which ranged in genome size from ~1.9 to ~3.1 Mb, and GC content between 53.6 and 62.7% (Table 3). With a coverage of 16x, Pt_3f was the most abundant organism among the recovered draft genomes, whereas Co_1 was the least abundant with a 10x coverage (Table 3).

Genes for the oxidation of reduced sulfur compounds were identified in all draft genomes, and additionally Co_1 and Co_2 contained sgpA, one of the genes encoding proteins for sulfur globule envelopes. Co_1 and Pt_3f contained the genes necessary for the complete oxidation of sulfide to sulfate (Table 4). Furthermore, Co_1 was the draft genome with the highest number of genes involved in nitrogen cycling, including the following two pathways: i) denitrification to dinitrogen, and ii) dissimilatory nitrate reduction to ammonium (DNRA). Nitrogenase-encoding genes nifH, nifD and nifK were absent in all draft genomes. Genes for two carbon fixation pathways were identified in Co_1, Co_2, and Pt_2d (cbbL and cbbS from the CBB cycle, and korAB, frdA, sdhA, sdhC, and porAB from rTCA cycle), whereas in Pt_3f genes for the CBB cycle were absent (Table 4). Comparative sequence analysis attributed these draft genomes to the family Chromatiaceae, and placed them in close proximity to the Solemya velum gill symbiont (IMG_2502171186; Figure 5).

Marine plants anchored in anoxic sediments harbor a steep redox gradient around their roots that creates microniches for different types of bacteria (Holmer and Nielsen, 1997; Blaabjerg and Finster, 1998; Devereux, 2005). By using 16S rRNA amplicon sequencing, we have found that the seagrass rhizosphere microbiome (“rhizobiome”) is enriched by Gamma- and Deltaproteobacteria (Cúcio et al., 2016). Although Gammaproteobacteria contains an extremely large and metabolically versatile group of bacteria, their importance in the sulfur cycle is indicated by the high abundance of members of the orders Chromatiales and Thiotrichales in the metagenomes. Members of the order Chromatiales are mainly phototrophic sulfur oxidizers (Imhoff, 2005). Even though most of these purple sulfur bacteria (PSB) require light as an energy source, some are able to grow chemotrophically under micro-oxic conditions in the absence of light (Imhoff, 2005). Chromatiales have been previously shown to be dominant in the rhizosphere of seagrasses (Cúcio et al., 2016) and salt marsh plants (Thomas et al., 2014). Some of the most abundant bacteria among the Chromatiales (Thiohalocapsa sp. and Lamprocystis purpurea) identified in the rhizosphere of Zostera marina from Portugal, are bacteria that can grow in the dark as chemolithoautotrophs or chemoorganotrophs using sulfide/thiosulfate and pyruvate, respectively, as electron donors (Eichler and Pfennig, 1988; Caumette et al., 1991). While profiting from the oxygen that is released from the roots during daytime (Borum et al., 2005; Frederiksen and Glud, 2006), these bacteria can actively consume the sulfide from the rhizosphere, and so contribute to a detoxification of the root area.

Nearly half of the class Deltaproteobacteria in our metagenomes was composed by members of the order Desulfobacterales. Ten percent of this order in ZmPt (corresponding to 5% of the Deltaproteobacteria and 2% of the Proteobacteria) were assigned to the species Desulfosarcina cetonica. The genus Desulfosarcina was among the most abundant taxa of the Deltaproteobacteria in the rhizosphere of the seagrasses Zostera marina, Z. noltii and Cymodocea nodosa (Cúcio et al., 2016). Desulfosarcina cetonica is a complete oxidizer that is able to utilize sulfate, thiosulfate and/or elemental sulfur as electron acceptors for the oxidation of a wide range of organic compounds, including acetone, benzoate, ethanol, lactate, and acetate (Galushko and Rozanova, 1991). As previously discussed by Cúcio et al. (2016), the presence of SRB capable of ethanol oxidation in the rhizosphere might represent a trade-off between plant and bacteria, in which the latter, although it produces sulfide, contributes to another mean of detoxification of the rhizosphere by consuming the ethanol released by the roots. Another deltaproteobacterium present in ZmPt was the (uncultured) PSCGC 5451, which was sequenced in a project targeting hydrocarbon degradation in marine sediments (GOLD–JGI). This strain, like several other SRB found (namely those belonging to the family Desulfobacteraceae), have been described to degrade hydrocarbons (Kleindienst et al., 2014). Desulfobacteraceae were also dominant among other Deltaproteobacteria on belowground compartments of the tropical seagrass Halophila stipulaceae (Mejia et al., 2016).

The overall community composition at the class level was similar in both locations, contradicting previous findings in salt marshes and seagrass meadows (Thomas et al., 2014; Cúcio et al., 2016). We consider that this discrepancy is likely to be associated with three reasons, i) different sequencing depth between the samples, as indicated by large variation in coverage (visible in Figure 2), (ii) biases inherent to PCR amplification in the former studies, and iii) to the high number (ca. 38%) of unclassified reads in the present study, which masks the presence, c.q. dominance of unknown bacteria. Furthermore, shotgun metagenomic sequencing does not feature the biases inherent to PCR amplification used in previous studies. Furthermore, the difference in the GC content found between both metagenomes (55.6 and 50.8% in ZmPt and ZmFr, respectively, Table 2), suggests that they harbor different microbial communities. Notwithstanding this fact, some sequences were abundant in both metagenomes (see clusters 1, 2, and 3, and encircled regions in Figure 2), highlighting an overrepresentation of some bacteria in these metagenomes, in particular in ZmFr. Cluster 1, in particular, was identified as a member of the Desulfobacteracae. Previous determination of the core rhizobiome of seagrasses (Cúcio et al., 2016) identified the Desulfobacteraceae as the most abundant family within the core rhizobiome; therefore it is not surprising that these bacteria are highly represented in both locations.

Although the importance of biological oxidation of sulfide in the rhizosphere of seagrasses has not yet been fully understood (Sayama et al., 2005; Preisler et al., 2007; Lenk et al., 2011), Luther et al. (2011) have demonstrated that biological sulfide oxidation rates are three or more orders of magnitude higher than abiotic rates, indicating the importance of sulfide oxidizing bacteria in the rhizosphere of seagrasses.

Overall, we observed a higher diversity and abundance of sulfur genes involved in oxidation than in reduction, such as the gene sqr encoding the enzyme sulfide:quinone reductase (Figure 3, Supplementary Figure 3). The gene sqr and its homolog fccAB, are essential for sulfide oxidation. They can be present in one or multiple copies per genome (e.g., Reinartz et al., 1998; Chan et al., 2009), and are characterized and divided into six different types, which exhibit different affinities to sulfide (Marcia et al., 2010). According to phylogenetic placement of the contig sequences annotated as SQR, we identified a majority of sequences clustering with the SQR types VI, I, and III. (Supplementary Figure S3). SQR type VI, responsible for growth on sulfide as the only electron donor (Chan et al., 2009), was the most abundant sulfur gene in ZmPt. This type has been shown to be transcribed under sulfide concentrations up to 8 mM in Chlorobium tepidum (Chan et al., 2009), indicating that sulfide-oxidizing bacteria present in the rhizosphere of Z. marina might be adapted to perform sulfide oxidation under high sulfide concentrations (Han and Perner, 2015). Nevertheless, the presence of other types of SQR in the rhizosphere of Z. marina with affinity to sulfide at the micromolar level, as well as a high diversity of FCC (Brune, 1995; Marcia et al., 2010), suggest that sulfide oxidizing bacteria possess a diverse set of genes which allow them to thrive under a large range of sulfide concentrations. Such versatility enables these bacteria to successfully explore the complex redox gradient and settle in specific microniches of the rhizosphere.

The oxidation of sulfide via SQR and/or FCC, as well as the oxidation of thiosulfate via the Sox multi-enzyme complex, result in the formation of sulfur globules as an intermediate product (Dahl and Prange, 2006; Eichinger et al., 2014). Although sulfur globule formation is not yet fully understood, the pathway for thiosulfate oxidation via SoxB seems to result in the accumulation of sulfur intra- or extracellularly (Dahl and Prange, 2006; Hensen et al., 2006; Eichinger et al., 2014). A study on endosymbiotic bacteria suggested that the formation of sulfur globules actively contributes to the detoxification of sulfide from the cells of the host (Eichinger et al., 2014), and functions as a reservoir of energy for the bacterium itself. Furthermore, reduced sulfur compounds are preferentially utilized as electron donors rather than organic compounds (Grimm et al., 2011). Intra- and extracellular sulfur globules are oxidized via the reverse-operating dissimilatory sulfite reductase, encoded by the genes rdsrAB (Dahl et al., 2005; Dahl and Prange, 2006; Hensen et al., 2006; Müller et al., 2015), which were abundant in ZmPt (Figure 3). The well represented presence of all these genes allowed us to speculate that bacteria in the seagrass rhizosphere are well adapted to varying concentrations and different types of reduced sulfur compounds, and furthermore, the high abundance of rdsrAB genes might represent an advantageous side mechanism for SOB to obtain energy (Müller et al., 2015).

Another abundant gene involved in oxidation, as determined by phylogenetic analysis (Figure 3, Supplementary Figure S1), was the APS reductase-encoding gene aprA. This is the reverse-operating homologous gene of the aprA found in sulfate-reducing prokaryotes (Meyer and Kuever, 2007a,b), which is essential for the oxidation of sulfite to an APS intermediate and is further oxidized to sulfate by the ATP sulfurylase (Sat).

The functional genes dsrAB, aprAB, and soxB are commonly used to infer the phylogenetic affiliation of SRB and SOB (Meyer and Kuever, 2007b; Meyer et al., 2007; Müller et al., 2015). Whereas aprAB and soxB genes have suffered several events of lateral gene transfer (LGT) among SRB and SOB (Meyer and Kuever, 2007b; Meyer et al., 2007), the dsrAB gene is highly conserved in both groups, and its phylogeny is congruent with that of the 16S rRNA gene (Müller et al., 2015). Despite that a majority of sequences of dsrAB clustered with those of uncultured bacteria (Figure 4), indicating the presence of a high diversity of unknown bacteria in the rhizosphere of Z. marina, another number of sequences analyzed in our study clustered in close proximity to those of well-known photo- and chemolithoautotrophic bacteria, such as Allochromatium vinosum (Kämpf and Pfennig, 1980), and (endo) symbionts of several marine invertebrates (Figure 4). Likewise, sequences of other target genes involved in oxidation processes, such as dsrAB, aprA, sat, sqr, and fcc also clustered closely with sequences from chemolithoautotrophic and endosymbiotic bacteria (Figure 4, Supplementary Figures S1–S3). In particular, several assembled sequences closely clustered with those of the gamma- and deltaproteobacterial endosymbionts of the gutless worm Olavius algarvensis (Figure 4). The gammaproteobacterial endosymbionts are SOB that have genes for autotrophic CO₂ fixation and accumulation of sulfur globules, and provide their gutless host with their vital nutrition (Woyke et al., 2006).

Three out of the four draft genomes isolated from the rhizosphere of Z. marina clustered in proximity to a sulfur-oxidizing gill symbiont of the clam Solemya velum. Gros et al. (2003) also found multiple 16S rDNA sequences identical to Codakia orbicularis symbiont. Nevertheless, no eukaryotic genes from this clam or from O. algarvensis were detected in our metagenomes, suggesting that these bacteria might be free-living forms of these symbionts. Free-living forms of endosymbionts of the clam Codakia orbicularis (Lucinidae) were previously identified in seagrass meadows (Gros et al., 2003; König et al., 2016). Furthermore, Gros et al. (1996, 2012) demonstrated that this clam is able to acquire its chemoautotrophic gill endosymbionts from the environment, although the release of these bacteria from the host to the environment has, to our knowledge, never been observed for these clams (Brissac et al., 2009). Symbionts of members of the family Solemyidae were initially thought to be exclusively vertically transmitted (Krueger et al., 1996), nevertheless recent evidence suggests a very dynamic mode of transmission, which further includes horizontal/environmental acquisition (Dmytrenko et al., 2014). Though little information is available about the transmission of chemoautotrophic symbionts of gutless oligochaetes, such as O. algarvensis, some evidence supports both vertical (Giere and Langheld, 1987) and environmental transmission (Dubilier et al., 2006). Hence, there seem to be several indications suggesting that seagrass meadows could harbor a pool of fauna endosymbionts in a free-livingstate.

Some seagrasses can, at least partly, cope with high levels of sulfide through a three-stage symbiosis between a Lucinidae clam and their chemolithoautotrophic gill symbionts (Van der Heide et al., 2012). However, this does not include meadows devoid of clams. Due to the large proportion of reads annotated as symbiotic bacteria and the absence of invertebrate hosts in our metagenome, we support the hypothesis of Brissac et al. (2009), which states that seagrass meadows function as a source of endosymbionts. Horizontal (environmental) symbiont transmission requires an inoculum of free-living bacteria (reviewed in Bright and Bulgheresi, 2010), which might originate by detaching from their symbiotic association upon the host's death, as recently shown in Riftia pachyptila by Klose et al. (2015). We hypothesize that these bacteria could accumulate and proliferate in the rhizosphere of seagrasses, until they are brought to the water column by bioturbators, for instance, and inoculate putative hosts in the periphery of the seagrass meadow.

Like the gill symbiont of Solemya velum, all four draft genomes recovered in our study possessed genes necessary for carbon fixation, as well as a repertoire of genes for the complete oxidation of sulfide to its most oxidized form sulfate (Table 4; Stewart et al., 2011), thereby confirming that these organisms can function as chemolithoautotrophs. Furthermore, they are also genetically equipped to perform nitrate reduction (Table 4), indicating that they may be able to couple sulfide oxidation to nitrate reduction, as previously described for bacteria such as Thioploca sp. and Beggiatoa sp. (Fossing et al., 1995; Otte et al., 1999; Sayama, 2001; Sayama et al., 2005).

A vast range of genes involved in the sulfur and nitrogen cycles, as well as carbon fixation, were identified in Co_1 (Table 4). The presence of alternative pathways in sulfur reactions (sulfate, thiosulfate, sulfur oxidation), carbon fixation mechanisms (CO₂ fixation via the CBB and carbon fixation via the TCA cycle) and nitrate reduction pathways (denitrification, DNRA) in this draft genome, are likely an advantageous trait in a dynamic environment as the seagrass rhizosphere. For this reason, we propose a conceptual model for the presence and function of Co_1-like free-living forms of symbiotic bacteria, in the rhizosphere of seagrasses (Figure 6). Having such a complete set of genes in its genome potentially gives this bacterium the ability to exploit a large variety of substrates, which is particularly important if it needs to thrive as free-living as well as in symbiotic associations. For instance, this bacterium is able to synthesize organic matter under oxic via the CBB cycle, and under micro-oxic/anoxic conditions via the TCA cycle (Figure 6). To our knowledge, the presence of these two pathways simultaneously has only been described in the endosymbionts of the hydrothermal vent worms Riftia pachyptila and Tevnia jerichonana, and represents an adaptation to transient conditions marked by the amount of energy available (Markert et al., 2007; Gardebrecht et al., 2012). Similar to the gill symbiont of Solemya velum (Stewart et al., 2011), this microorganism is capable of oxidizing sulfide at low and high environmental concentrations using FCC and SQR, respectively (Brune, 1995; Marcia et al., 2010). Moreoever, it might be able to store sulfur intracellularly as suggested by the presence of the gene sgpA. Sulfide concentrations in the rhizosphere of seagrasses can be very dynamic and particularly respond to light/dark regimes, mainly due to the lower oxygen release in the dark, compared to light conditions (e.g., Pagès et al., 2012). Given its ability to thrive in such conditions, this bacterium also has the genetic setup necessary to switch their usage between nitrogen compounds, according to the availability of sulfide. Brunet and Garcia-Gil (1996) found evidence that under low sulfide concentrations denitrification is favored, whereas under high levels of sulfide, DNRA is more likely to proceed (as sulfide inhibits the final two steps of denitrification, Figure 6). In this regard, high sulfide concentrations may benefit the seagrass (below the threshold of sulfide detrimental for the fitness of the plant), by inducing the production of ammonium, which is the seagrasses' preferred nitrogen source (Alexandre et al., 2015). Like sulfide, ammonium concentrations in the rhizosphere increase during dark periods (Pagès et al., 2012). Although these authors attributed higher NH4+ concentrations to a decrease in seagrass uptake, according to our conceptual model, high sulfide concentrations sustain DNRA, which results in the release of ammonium during the night.

Notwithstanding the fact that nitrate reduction via DNRA is favored by environmental conditions (Sayama, 2001), Co_1 possesses genes that encode nitric oxide and nitrous oxide reductases. The presence of these genes indicates that this organism could have the potential to completely denitrify nitrate to dinitrogen, however we did not find the genes that encode the nitrite reductase responsible for the reduction of nitrite to nitric oxide (nir, Figure 6). Either nir genes were not captured due to the limited sequencing depth (this genome is missing approximately 20% of its genes), or they are naturally not present in Co_1. While we cannot rule out the first possibility, complete denitrification with the concomitant absence of nir genes in Co_1 could be bridged by the reduction of nitrite to nitric oxide through interaction with sulfide (Figure 6, Grossi, 2009; Cortese-Krott et al., 2015). This abiotic reaction could be of major importance for the seagrasses, because it would allow the double detoxification of two toxic components, H₂S and NO2-.

The present study unveiled the so-far unknown diversity of sulfur genes present in the rhizosphere of the seagrass Zostera marina, and supported the hypothesis that the rhizobiome is enriched with sulfur-oxidizing bacteria. We found indications supported by phylogenetic inference that the rhizosphere hosts highly versatile bacteria related to symbionts of marine invertebrates that are capable of exploiting a wide range of sulfur and nitrogen compounds, using alternate pathways that are favorable under different environmental conditions. We suggested a conceptual model for the sulfur, nitrogen and carbon fixation metabolism of these organisms in the rhizosphere of seagrasses. The draft genomes recovered further supported the presence of chemolithoautotrophic bacteria closely related to free-living forms of symbionts, that are able to couple the complete oxidation of sulfide to sulfate with nitrate reduction to ammonium. The dominance of these bacteria should be determined by other techniques, such as fluorescent in situ hybridisation (FISH) or qPCR, and their contribution to the well-being of seagrasses needs to be unraveled in future studies by characterization of their activity.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.