Prevalence of Heterotrophic Methylmercury Detoxifying Bacteria across Oceanic Regions

Microbial reduction of inorganic divalent mercury (Hg2+) and methylmercury (MeHg) demethylation is performed by the mer operon, specifically by merA and merB genes, respectively, but little is known about the mercury tolerance capacity of marine microorganisms and its prevalence in the ocean. Here, combining culture-dependent analyses with metagenomic and metatranscriptomic data, we show that marine bacteria that encode mer genes are widespread and active in the global ocean. We explored the distribution of these genes in 290 marine heterotrophic bacteria (Alteromonas and Marinobacter spp.) isolated from different oceanographic regions and depths, and assessed their tolerance to diverse concentrations of Hg2+ and MeHg. In particular, the Alteromonas sp. ISS312 strain presented the highest tolerance capacity and a degradation efficiency for MeHg of 98.2% in 24 h. Fragment recruitment analyses of Alteromonas sp. genomes (ISS312 strain and its associated reconstructed metagenome assembled genome MAG-0289) against microbial bathypelagic metagenomes confirm their prevalence in the deep ocean. Moreover, we retrieved 54 merA and 6 merB genes variants related to the Alteromonas sp. ISS312 strain from global metagenomes and metatranscriptomes from Tara Oceans. Our findings highlight the biological reductive MeHg degradation as a relevant pathway of the ocean Hg biogeochemical cycle.


S2
within the same genome, the designed primers only covered the gene/s copies with the highest similarity between them and between the different Alteromonas species, since evidence of horizontal gene transfer has been detected in mer genes, usually in plasmids, transposons or genomic islands 5 (Supplementary Table S1). On the other hand, as Marinobacter merA sequences differ greatly from species to species and most of the Marinobacter isolates present in our culture collection taxonomically classified as Marinobacter hydrocarbonoclasticus, Marinobacter aquaeolei or Marinobacter salarius (Supplementary Table S2), primers were designed to cover the merA sequence variants of these species (Supplementary Table S1).
Otherwise, as merB genes were only found in Alteromonas mediterranea DE strain (CP003917) and Marinobacter aquaeolei VT8 strains (NC_008740) (also named Marinobacter hydrocarbonoclasticus VT8), merA + merB sets of primers were designed specifically for these mentioned strains. In this case we used the online tool Primer-BLAST of the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) 6 . The input sequence for the generation of these primers was a concatenate nucleotide sequence of the merA and merB genes that were co-localized one next to the other in the mer operon (Supplementary Figure S1A). All the designed primers had to meet the following requirements: optimum polymerase chain reaction (PCR) product size of 1200 bp for merA or 2100 bp for merA + merB (referred hereafter as merAB), annealing temperature around 57 ºC, primers length 20 bp, and 50 % GC content. Designing optimal PCR primer sequences following the above mention parameters allows a successful PCR with sensitive, specific, and reproducible results avoiding crosshybridization with primers and other gene sequences present in the reaction mixture 7

Phylogenetic trees with aminoacid merA and merAB sequences
On the other hand, phylogenetic trees were also constructed with the amino acid sequences of the amplified merA and merAB genes. Using Kyoto Encyclopedia of Genes and Genomes (KEGG) identifiers K00520 (for merA) and K00221(for merB) we retrieved the merA and merB genes present in UniProtKB. Those sequences were then submitted to BLAST against our detected merA and merB genes, and best hits were included in the tree as reference sequences. Sequences were aligned with ClustalW of the Geneious software v.11.0.5 8 with the Gonnet substitution matrix and default gap extension and opening penalties as described previously 9 . For merA, the dihydrolipoamide dehydrogenase protein sequences from Magnetospirillum magneticum AMB-1 (WP_011386317.1) and Pseudomonas fluorescens Pf0-1 (WP_011336663.1) served as outgroups. Likewise, we trimmed the N-terminal region of the aligned sequences as in Barkay et al. (2010) 10 . In the case of merAB, outgroup sequences were not included. Phylogenetic trees were constructed using maximum-likelihood inference with RAXML-NG 0.9.0 11 and the LG evolutionary model with optimization in the among-site rate heterogeneity model and the proportion of invariant sites (LG+G+I), and 100 bootstrap replicates.

Measurement of the biotic and abiotic degradation of MeHg
In order to characterize the MeHg degradation rates of our bacteria caused by the action of the merA and merB genes, 2 ml samples were taken from the 1 µM and 5 µM growth curves of the most tolerant strain, Alteromonas sp. ISS312, at times 0, 6, 12, 24 and 48 h. Besides, in order to check the possibility that MeHg was being abiotically removed, we measured the MeHg concentrations from samples taken from multiwell plates experiments. These included one well of a liquid culture from the ISS312 strain in Zobell broth (initial O.D. at 600 nm of 0.05) amended with 5 µM CH 3 HgCl and incubated at RT during 72 h in the dark, as well as two different controls to detect the possible abiotic degradation of MeHg: (i) medium control (CH 3 HgCl at concentrations of 1 µM and 5 µM with no strain added), and (ii) killed control (CH 3 HgCl at concentrations of 1 µM and 5 µM added to autoclaved liquid cultures of the strain). Three replicates of 2 ml samples from each well, containing the CH 3 HgCl treatment and the different controls, were taken at times 0 and 72 h. Collected samples were immediately frozen at -80 ºC.
Concentration of Hg species was measured by direct derivatization of the culture samples with sodium tetraethyl borate and injection into a hyphenated system consisting of a gas chromatograph coupled to an atomic fluorescence detector via pyrolysis (GC-pyro-AFS) as previously described elsewhere 12 . Briefly, 2 ml samples were used for derivatization. The pH of the extracts was adjusted to 3.9 by adding 5 ml of 0.1 M acetic acid-sodium acetate buffer and ammonia (20 %) if necessary. Then, 2 ml of hexane and 250 µl of sodium tetraethyl borate (6 %, S4 w/v) were added and the mixture was manually shaken for 5 min. The sample was centrifuged for 5 min at 600 g. The organic layer was transferred to a chromatographic glass vial and stored at -18 ºC until analysis. When Hg species were not detectable, the organic layer was preconcentrated under a gentle stream of nitrogen to a low volume (50-100 µL) just before the measurement. The procedural detection limits, after pre-concentration, were 0.19 and 0.23 nM for MeHg and inorganic Hg, respectively.
Strain ISS312 genome sequencing DNA from strain ISS312 was extracted using the DNeasy Blood & Tissue kit (Qiagen), following the manufacturer's recommendations. A genome library was prepared with a Celero TM DNA-seq library system and sequencing was performed with paired-ended 300 bp long reads by IGA-Tech with a MiSeq Illumina machine. The sequence data was filtered to remove the adapters and the unpaired reads with cutadapt v.1.16 13 and the quality was assessed before and after with fastaqc v.0.11.7. The clean data were used to do the assembly with Spades v.3.12.0 and then optimized with the same program. QUAST v.5.0.1 14 and ALE were used to assess the quality of the assemblies and the best scores were selected. K-mer 121, 125 and 127 were selected for the optimized-combined assembly and the quality was assessed again. The annotation was done with Prokka v1.13 15 and the completeness of the genome was checked with CheckM v1.0.18 16 . In order to search for plasmids within our contigs, we used the database PLSDB 17 .

Fragment recruitment analysis of the genome of strain ISS312 and MAG-0289 in bathypelagic metagenomes
The abundance of ISS312 across the global bathypelagic ocean was assessed thanks to its complete genome. Fragment Recruitment Analysis (FRA) was performed by mapping the metagenomic reads of 58 bathypelagic microbial metagenomes from 32 stations 18 from Malaspina Expedition, including free-living (0.2-0.8 µm) and particle-attached (0.8-20 µm) bacterial communities. Analyses were done with BLASTn v2.7.1+ 3 using the following alignment parameters: -perc_identity 70, -evalue 0.0001. Only those reads with more than 90 % coverage and mapping identity equal to or higher than 95 % were kept for analysis. In order to remove possible false mapping hits to the conserved regions of rRNA genes, reads aligning to the regions annotated as ribosomal genes were not considered for the analysis. Read counts from mapped reads from each metagenome were corrected by their sequencing depth to make them comparable through samples.
The abundance of the MAG-0289 across the Malaspina bathypelagic dataset was done as described above, but with a previous subsampling step to the shallower sequencing depth (4,175,346 read pairs) with bbtools reformat.sh v.38.08 (https://sourceforge.net/projects/bbmap/).
Here, when reads could map with the same probability to any of the genomes (same e-value, same alignment length and identity), they were assigned at random. The average nucleotide identity (ANI) of ISS312 against MAG-0289 was calculated with fastANI v1.2 as both were related to the Alteromonas mediterranea species.

Preparation and observation of samples for transmission electron microscopy
In order to study potential effects of MeHg on strain ISS312 cells morphology, the isolate was grown separately during 24 h in Zobell broth not amended with CH 3 HgCl and amended at a final concentration of 5 µM with shaking in the dark. The overnight culture was centrifuged at 1000 g during 15 min and the supernatant was discarded. The pellet was fixed with paraformaldehyde 2 % final concentration during 30 min at RT. After fixation, the pellet was processed as previously described 19 to finally obtain thin sections of the samples that were examined by transmission electron microscopy (TEM, JEM-1400 plus, JEOL). Visualizations were done by the microscopy service of the Universitat Autònoma de Barcelona (http://sct.uab.cat/microscopia/en/content/inici).

Exploring the MARINHET culture collection
We detected a total of 20 different bacterial taxa in the IMG/JGI database that matched at the genera level with the taxonomic assignation of the MARINHET culture collection isolates, and therefore, with potential for carrying merA and merB genes (Supplementary Table S9). The comparative analyses between the 16S rRNA gene sequences of these 20 genera containing the targeted merA and merB genes and the partial 16S rRNA sequences of our isolates revealed a total of 352 strains that were, at least, 99 % identical to one of the putative candidates genera. These  Table S10). From these, we selected a total of 244 strains affiliating to Alteromonas sp. and 46 strains to Marinobacter sp. (Supplementary   Table S2) for performing the functional screening of merA and merB genes through PCR.

Ribosomal phylogenetic reconstruction of strains codifying merA and/or merB genes
The phylogenetic analyses of the 16S rRNA gene from all of the detected strains harboring merA and merAB genes, together with the rest of the screened strains, revealed some patterns (Supplementary Figure S7). First, most of the isolated Alteromonas strains with merA were related to Alteromonas australica and Alteromonas mediterranea, and only strains affiliating to the last one presented both genes merA and merB (merAB). We also detected one strain with merAB genes affiliated to Alteromonas macleodii. Secondly, Marinobacter strains displaying merAB were related to Marinobacter hydrocarbonoclasticus, Marinobacter salarius and some uncultured Marinobacter strains. We are aware that our primers do not match the whole diversity within Alteromonas since there is a large uncultured Alteromonas cluster related to NW Mediterranean, Indian Ocean and North and South Atlantic samples and from photic and aphotic layers that was not covered by our primer-set (Supplementary Figure S7). Moreover, this phylogeny showed that positive strains for merA and merB genes clustered together with strains which do not harbor these genes (Supplementary Figure S7). It has been described that some bacterial species codify different sequence variants of the merA gene 9,20 including the Alteromonas genus 5 . Therefore, it is possible that some of the Alteromonas and Marinobacter strains within the MARINHET collection tested presented different sequence variants, although we were only able to detect the ones targeted by the primers designed. On the other hand, the operon mer can be either codified in the chromosome 21 or in plasmids 22,23 , and usually, mer genes are components of transposons 24 , and integrons 25,26 . Thus, it is not surprising to find some strains within the same species without the mer operon.

Results from merA and merAB phylogenies
The amino acid merA phylogeny, which theoretically could include the different gene sequence variants covered by the primers designed, unveiled that all the Alteromonas merA sequences Marinobacter spp. Alteromonas spp.

5' TTG CTT SAC RTC CTT GGT 5' 5' CCT TCG TGT GCT AAG ATT TG 3'
Marinobacter spp. Alteromonas spp.    Tables   Table S1. Downloaded merA and merB sequences from the JGI/IMG database (2016) to design the PCR primers. In green are selected those merA and merB copies that should be recognized by the primers designed. Table S2. Information of the Alteromonas sp. and Marinobacter sp. strains used for PCR screening of the merA and merB genes. Table S3. Characteristics of the marine seawater samples from five different cruises used for the isolation of bacteria screened for Hg resistance genes and results from PCR screening. ATP09: Arctic Tipping Points cruise in 2009; BBMO: Blanes Bay Microbial Observatory. Table S4. Summary of the number of positive strains per oceanographic region. Table S5. Summary of the results for the MIC determination assays to HgCl 2 . Table S6. Concentrations of MeHg and inorganic mercury (IHg) at different time points during the growth curves at 1 µM and 5 µM. SD: standard deviation; LOD: limit of detection. Table S7. Biotic (culture) and abiotic controls (killed and medium alone) of the 5 µM MeHg degradation performed by the ISS312 Alteromonas strain. SD: standard deviation; LOD: limit of detection. Table S8. Taxonomic classification of the 54 merA and 6 merB gene homologous to merA and merB genes from ISS312 strain found in global metagenomes and metatranscriptomes from Tara Oceans expedition. Table S9. Candidate genera for mercury bioremediation obtained from the search of merA and merB genes in the KEGG database and from the IMG/JGI database. Table S10. Summary of the putative isolates within the MARINHET culture collection that could harbor merA and merB genes based on the BLASTn search between the partial 16S rRNA genes sequences of the isolates and the 16S rRNA gene sequences of the candidate genera extracted from the JGI/IMG database.