Metagenomic Analysis Shows the Presence of Bacteria Related to Free-Living Forms of Sulfur-Oxidizing Chemolithoautotrophic Symbionts in the Rhizosphere of the Seagrass Zostera marina

Seagrasses play an important role as ecosystem engineers; they provide shelter to many animals and improve water quality by filtering out nutrients and by controlling pathogens. Moreover, their rhizosphere promotes a myriad of microbial interactions and processes, which are dominated by microorganisms involved in the sulfur cycle. This study provides a detailed insight into the metabolic sulfur pathways in the rhizobiome of the seagrass Zostera marina, a dominant seagrass species across the temperate northern hemisphere. Shotgun metagenomic sequencing revealed the relative dominance of Gamma- and Deltaproteobacteria, and comparative analysis of sulfur genes identified a higher abundance of genes related to sulfur oxidation than sulfate reduction. We retrieved four high-quality draft genomes that are closely related to the gill symbiont of the clam Solemya velum, which suggests the presence of putative free-living forms of symbiotic bacteria. These are potentially highly versatile chemolithoautotrophic bacteria, able to alternate their metabolism between parallel pathways of sulfide oxidation (via sqr and fcc), nitrate reduction (denitrification or DNRA) and carbon fixation (via CBB or TCA cycle), depending on the environmental availability of sulfide. Our results support the hypothesis that seagrass meadows might function as a source of symbionts for invertebrates that inhabit within or around seagrass meadows. While providing ideal conditions for the proliferation of these free-living forms of symbionts, seagrasses would benefit from their genetic versatility, which contributes to sulfide detoxification and ammonium production, the seagrasses' preferred nitrogen source.


INTRODUCTION
Global carbon and sulfur biogeochemical cycles are tightly coupled in marine sediments, mostly through sulfate reduction, which is responsible for approximately 50% of the organic carbon remineralization in anoxic coastal shelf sediments (Jørgensen, 1982). Although marine sediments in general harbor active biogeochemical cycles, these seem to be further boosted in seagrass-vegetated areas (Devereux, 2005). The rhizosphere of seagrasses is a thin layer of sediment surrounding the roots, which is highly enriched with dissolved organic matter (OM), mostly originated from photosynthetic products that are released through the roots (Holmer and Nielsen, 1997;Pérez et al., 2007). These root exudates create microniches that are distributed along the rhizosphere according to the availability of the most energetically favorable terminal electron acceptors for microbial respiration (Capone and Kiene, 1988;Devereux, 2005), creating hotspots of microbial activity (Holmer and Nielsen, 1997;Blaabjerg and Finster, 1998;Donnelly and Herbert, 1999). The high abundance of sulfate in seawater (Capone and Kiene, 1988) coupled with OM enrichment in the rhizosphere, enhance the activity of sulfate reducing bacteria (SRB, Holmer and Nielsen, 1997;Blaabjerg and Finster, 1998;Donnelly and Herbert, 1999). This results in the production of high levels of sulfide, which pose a strong threat to seagrass health and survival (Borum et al., 2005;Pérez et al., 2007). Detoxification of sulfide in the rhizosphere occurs by different chemical and biological processes. Chemical oxidation can occur by oxygen released from the seagrass roots during photosynthesis (Jørgensen and Nelson, 2004;Borum et al., 2005;Frederiksen and Glud, 2006;Holmer et al., 2006), or by binding to iron resulting in iron sulfide and pyrite (Jørgensen and Nelson, 2004). Although the majority of studies focus on the chemical oxidation of sulfide, its biological oxidation by sulfideoxidizing bacteria (SOB) also plays an important role in coastal marine sediments (Jørgensen and Nelson, 2004;Frigaard and Dahl, 2008).
Zostera marina, also known as eelgrass, is the most widespread seagrass species in the world (Green and Short, 2003). It is found in subtidal areas up to 15 meters depth (Borum and Greve, 2004) along temperate coasts (Green and Short, 2003) throughout the northern hemisphere. Eelgrasses, as seagrasses in general, are ecosystem engineers. They provide valuable ecosystem services, such as a habitat, refuge and nursery ground for many animals, improve water quality through sediment and organic matter retention, and by filtering out nutrients and contaminants (Gacia et al., 1999;Short et al., 2000). Recently it was shown that seagrasses might even control the abundance of potential pathogens of marine animals and humans (Lamb et al., 2017).
In a previous study in which we used 16S rRNA gene sequencing, we found that bacteria involved in the sulfur cycle are abundant in the rhizosphere of different European seagrass species (Cúcio et al., 2016). Although the importance of biogeochemical sulfur cycling in coastal marine sediments has been widely studied, to our knowledge, no information is available about the molecular pathways involved in sulfate reduction and sulfide oxidation in the seagrass rhizosphere. In the present study, we investigated the diversity, structure and gene content of bacterial communities involved in dissimilatory sulfur cycling in the rhizosphere of Zostera marina, and hypothesized that the rhizosphere harbors a diverse community of bacteria involved in the sulfur cycle, which is enriched in sulfide oxidizers over sulfate reducers. Furthermore, we reconstructed and characterized four gammaproteobacterial draft genomes. Shotgun metagenomic sequencing allowed us to investigate for the first time in detail the molecular pathways involved in sulfate reduction and sulfide oxidation in the rhizosphere.

Sample Collection and Preparation
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.

Metagenomic Sequencing
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.
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 precorrection 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).

Phylogenetic Reconstruction of Sulfur Genes
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 maximumlikelihood 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.

Reconstruction and Analysis of Draft Genomes
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 2 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). Nitrogenaseencoding 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.

RESULTS
Metagenomic sequencing resulted in approximately 98 and 91 million high-quality paired-end reads from the rhizospheres sampled in Portugal (ZmPt) and France (ZmFr), respectively. The metagenome ZmPt yielded better quality results than the metagenome ZmFr, i.e., more than twice as many contigs larger than 1 kb and a N50 of 1,150 bp, compared to a N50 of 1,097 bp in ZmFr ( Table 2). The co-assembly resulted in a total of 76,983 contigs larger than 1 kb and a N50 of 1,768 bp, higher values than those observed for ZmPt and ZmFr ( Table 2). The average G+C content of the contigs from ZmPt, ZmFr and the co-assembly was 55.61, 50.79, and 54.20%, respectively ( Table 2).

Bacterial Community Structure
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  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 coassembly 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.
Diversity of Sulfur Genes in the Rhizobiome of Z. marina FIGURE 1 | Community composition of the microbiomes from the rhizosphere of Zostera marina from Portugal (ZmPt) and France (ZmFr). Taxonomic classification of metagenomic reads of bacterial classes (bar charts), and orders of the two most abundant classes, the Gamma-and the Deltaproteobacteria (pie charts). All unassigned reads, as well as those unclassified at higher taxonomic ranks were excluded from the analysis. Taxonomy assignment was performed on raw reads in the software program Kaiju (Menzel et al., 2016).
FIGURE 2 | Mapping of reads to the co-assembled metagenomes. Visualization of the contigs from Portugal (blue, ZmPt) and France (green, ZmFr) mapped against the co-assembled metagenomes, represented in a log scale with a maximum of 2.2 for the former and 3.2 for the latter, respectively. The GC content is also represented, with a maximum of 73%. The dendogram was calculated based on sequence composition and differential coverage of raw reads.
sulfur pathways dominated ZmFr in terms of coverage (average coverage of 20x vs. 24x for oxidative and reductive sulfur genes, respectively, Figure 3).
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).

Reconstructed Draft Genomes
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).

Bacterial Community Structure
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 Frontiers in Marine Science | www.frontiersin.org 6 May 2018 | Volume 5 | Article 171 FIGURE 4 | Diversity of dissimilatory sulfite reductase (dsrAB) genes. Phylogenetic tree of dsrAB marker genes present in the rhizobiomes of Zostera marina from Portugal (red, prefix PT) and France (blue, prefix FR). Contigs labeled with suffix "-a" correspond to sequences annotated to dsrA, "-b" to dsrB, and "-ab" to concatenated dsrAB. Sequences were concatenated whenever dsrA and dsrB were assembled in the same contig. Contigs were aligned at the protein level to the dsrAB reference alignment published by Müller et al. (2015) using the ARB alignment tool. The phylogenetic tree was inferred using the approximate maximum-likelihood method in FastTree2 (Price et al., 2010). The green background indicates sequences involved in oxidative processes, the blue background sequences involved in reductive processes. Bootstrap values are indicated. Scale bar indicates percentage sequence difference.  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.

Diversity of Sulfur Genes in the Rhizobiome of Z. marina
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 reverseoperating homologous gene of the aprA found in sulfatereducing 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).
FIGURE 5 | Phylogenetic affiliation of reconstructed draft genomes. Draft genomes were aligned to genomes in the CheckM reference tree (Parks et al., 2015). Subsequently, sequences were selected and the final phylogeny was calculated using the approximate maximum-likelihood method in FastTree2 (Price et al., 2010). Bootstrap values are indicated. Scale bar indicates percentage sequence difference.

Symbiont-Related Bacteria in the Rhizosphere of Z. marina
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 2 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 sulfuroxidizing 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. (1996Gros et al. ( , 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 FIGURE 6 | Conceptual model of sulfur, nitrogen and carbon fixation pathways in one of the reconstructed genomes, CO_1, obtained from the rhizosphere of Zostera marina. At aerobic conditions and low sulfide concentrations fcc is used to oxidize sulfide with oxygen, while at high sulfide concentrations sqr is used to oxidize sulfide. In both cases, the Calvin cycle (CBB) is used to fix CO 2 . However, at microaerobic and anaerobic conditions and low sulfide concentrations fcc is oxidizing sulfide via denitrification (purple), while at high sulfide concentrations sqr is oxidizing sulfide via DNRA (dissimilatory nitrate reduction to ammonium represented in blue). At these conditions carbon is fixed using the TCA-cycle. The ammonium produced in the DNRA is used as a nitrogen source by the seagrasses. The dashed arrow represents the proposed interactions between the sulfur and nitrogen cycle in which sulfide reduces nitrite to nitric acid.
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-living state.
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.

Reconstructed Symbiont-Related Draft Genomes
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 2 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 freeliving 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 NH + 4 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 2 S and NO − 2 .
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.

DATA AVAILABILITY
The raw sequence reads have been deposited as dataset SRP126211 in the NCBI Sequence Read Archive (SRA).

AUTHOR CONTRIBUTIONS
CC, AE, and GM designed the study and collected the samples; CC performed all lab experiments, and data analysis together with LO; CC wrote the manuscript and all authors contributed to the discussion of the results and to the final version of the manuscript.

ACKNOWLEDGMENTS
We thank Francisco Rodríguez-Valera, Mário Lopez-Lopez, and Rohit Ghai for an introduction into metagenomic analysis, Emily D. Melton for helpful discussions about microbial sulfur metabolism, Tom Berben for very helpful contribution to scripting, and Nicole Dubilier for insightful comments on the final version of the manuscript.