Metagenomic survey of methanesulfonic acid (MSA) catabolic genes in an Atlantic Ocean surface water sample and in a partial enrichment

Seawater sample collection and metagenomic DNA isolation

Atlantic Ocean surface water was collected on Dec 3rd, 2014 along the coast of Leça da Palmeira, Portugal (approximate coordinates 41.226956, −8.720528). Approximately 15 L of seawater were collected off the rocky shore at rising tide into clean bottles, which were immediately transported to the lab in an isothermal bag with ice packs. Before starting the procedures, the bottles were shaken in order to homogenize the samples. Three fractions of 0.5 L were filtered through 1.2 µm glass fiber filters, which were immediately wrapped in aluminum foil and stored at −80°C for later quantification of chlorophyll. The determination of chlorophyll a concentration was performed through a spectrophotometric method as previously described (Inag, 2009). Chlorophyll a concentration was calculated through the Jeffrey and Humphrey equation (Jeffrey & Humphrey, 1975). Measures for pH and conductivity were also performed. Two fractions of 5 L of the seawater sample were filtered in parallel through 1.2 µm, 0.45 µm and 0.2 µm filters in succession. The filters from one 5 L fraction were immediately used for DNA extraction using the PowerWater DNA Isolation Kit (MO BIO Laboratories, Inc.) according to the manufacturer’s instructions.

Marine enrichment with MSA and metagenomic DNA isolation

The biomass retained on the filters from the second 5 L fraction was resuspended in the last 200 mL of the ocean water sample. The filters were removed and 1 mL of alkaline (pH 8) sodium methanesulfonate (MSA) 1 M was added to the suspension (5 mM final concentration). Similar amounts of alkaline MSA were used to spike the enrichment at days 7, 9 and 12. The suspension was incubated aerobically at room temperature (ca. 20 °C) in the dark, continuously mixed by a magnetic stirrer. This schedule was deliberately maintained during 16 days, with no subculturing, in order to produce just a partially enriched culture. The biomass at the end of the enrichment process was collected by centrifugation and DNA was extracted with PowerWater DNA Isolation Kit (MO BIO Laboratories, Inc., with adaptations to accommodate biomass in a pellet rather than on a filter). For convenience, the seawater sample (time 0) and the corresponding enrichment (E) culture will be designated SCD0 and SCDE, respectively. This study was limited to the observation of the evolution of the microbial community of a single sample and as such the following results are to be considered exploratory.

Whole metagenome sequencing

The metagenomes (from SCD0 and SCDE samples) were sequenced at Molecular Research LP (Shallowater, TX, USA). Paired-end sequencing libraries were prepared (2 × 101 bp) and sequencing was performed using the Illumina HiSeq2500 platform. The libraries were prepared using Nextera DNA Sample preparation kit (Illumina) following the manufacturer’s user guide. The initial concentration of DNA was evaluated using the Qubit^® dsDNA HS Assay Kit (Life Technologies). The samples were then diluted accordingly to achieve the recommended DNA input of 50 ng at a concentration of 2.5 ng/µL. Subsequently, the samples underwent fragmentation, addition of adapter sequences and PCR amplification (5 cycles) during which a unique index was added to each sample. The average library size was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies). The libraries were then pooled in equimolar ratios at 2 nM and 5 µL of the library pool was clustered using the cBot (Illumina) and sequenced paired-end for 200 cycles using the HiSeq 2500 system (Illumina).

The quality of the sequencing reads from both libraries was assessed using FastQC on Galaxy web-based platform (https://usegalaxy.org/). The libraries were also checked for human contamination with Kraken Metagenomics (Zaharia et al., 2011; Wood & Salzberg, 2014) at Illumina BaseSpace (http://basespace.illumina.com/home/index). As the resulting values for the presence of human sequences were very low (0.15% for sample SCD0 and 0.04% for sample SCDE), no filtering was performed before the subsequent steps of analysis (general statistics are described in Table S1).

Sequencing data for the two samples, SCD0 and SCDE, were submitted to the European Nucleotide Archive (http://www.ebi.ac.uk/ena) under project number PRJEB9018 and sample accession numbers ERS700852 and ERS700853, respectively. The analyses of the metadata were performed through the EBI Metagenomics service pipeline (https://www.ebi.ac.uk/metagenomics/pipelines/2.0 (Mitchell et al., 2015a)) that includes a quality control step, a taxonomic analysis step based on 16S rDNA sequences and a functional analysis of predicted protein coding sequences using the InterPro resource (Mitchell et al., 2015b).

BIOM (see http://biom-format.org/) files containing phylogenetic classification information provided by EBI metagenomics were used to construct rarefaction curves using MEGAN (version 5.10.6; Huson, Mitra & Ruscheweyh, 2011). Phylogenetic data were also used to estimate the alpha diversity of the two samples (Shannon index (Shannon, 1948), evenness (Mulder et al., 2004) and Chao species estimator (Chao, 1984)). For beta diversity analysis, Jaccard (Jaccard, 1912), Kulczynski (Faith, Minchin & Belbin, 1987), and Chao (Chao, Chazdon & Shen, 2005) indices were calculated through the Vegan package (Oksanen et al., 2015). Bray-Curtis dissimilarity index (Bray & Curtis, 1957) was calculated starting from relative abundances. Only significant differences in phylogenetic or functional composition were retained (Fisher’s exact / χ² test with a significance level of 0.05 with a Bonferroni correction on the number of comparisons, in order to minimize false positives; code was adapted from Metastats (White, Nagarajan & Pop, 2009)).

Assembled metagenomic data from both samples were also submitted to the DOE Joint Genome Institute’s Integrated Microbial Genome Metagenomic Expert Review (IMG/MER) annotation pipeline (http://img.jgi.doe.gov/) for functional and taxonomic annotation (Markowitz et al., 2014) (SCD0: GOLD Analysis Project Id Ga0069134/biosample ID Gb0111627; SCDE: GOLD Analysis Project Id Ga0069135/biosample ID Gb0111630). Before submission, reads were quality checked (FastQC for quality control at Galaxy for trimming and filtering): only sequences with quality scores equal or higher than 20 over 95% or more of the nucleotides were kept. For each sample, forward and reverse sequence files were merged and assembled using Megahit (Li et al., 2015) (general statistics are reported in Table S2). A further analysis was performed at MG-RAST (Meyer et al., 2008) using the subsystems approach (Overbeek et al., 2005): metagenome SCD0 was submitted under sample no 4698364.3 and SCDE under sample no 4698363.3. Binning was performed using MetaBat v0.25.4 (Kang et al., 2015) and the bins obtained were analyzed by MG-RAST.

A flowchart with the major steps of the analysis is available in Fig. S1.

msmA and msmE amplicon surveys

In order to obtain DNA amounts sufficient for the subsequent analyses, the REPLI-g^® MiniKit (QIAGEN) was used to amplify the DNA from each metagenomic sample, according to the manufacturer’s instructions. Negative controls starting from ultrapure nuclease-free water yielded no product.

The amplified metagenomes were tested with primers directed to genes msmA and msmE. As amplicons with less than 400 bp were needed for Ion Torrent sequencing, primer pairs SarA139fwd/SarA488rev and SarE322fwd/SarE704rev (Table S3) were employed in order to obtain amplicons with 349 bp and 382 bp, respectively. However, when amplification was performed directly with these primers, the signals obtained were too weak or absent. As such, a nested PCR approach was carried out. The first round was performed using primer sets SarA124fwd/SarA1053rev or SarE133fwd/SarE1125rev for msmA and msmE amplicons, respectively (Table S3). Some parameters had to be adjusted in order optimize the reactions, and the successful conditions are reported in Table S4. Negative controls received PCR water instead of DNA. Positive controls contained DNA from Sargasso Sea Metagenome clone EF103447 (Leitão, Moradas-Ferreira & De Marco, 2009). Only primers directed to the low-GC Sargasso Sea Metagenome msm sequences had previously been successful at amplifying these genes from a seawater metagenomic sample from the same location (Henriques & De Marco, 2015a). Accordingly, the primers used in this work were based on the known Sargasso Sea Metagenome msm sequences. Amplification products were submitted for sequencing at Stabvida Lda. (Caparica, Portugal). After determining the exact PCR product concentrations with Qubit 2.0 Fluorometer (Invitrogen) and Qubit dsDNA BR Kit (Invitrogen), platform-specific barcoded adapters were added to each sample in the preparation of libraries with KAPA Library Preparation Kit (Kapa Biosystems) and NEXTflex™ DNA Barcodes for Ion PGM (Bioo Scientific). Samples were sequenced on an Ion Torrent Personal Genome Machine (PGM, Life Technologies, Thermo Fisher Scientific) using Ion PGM™ Hi-Q Sequencing Kit reagents. Reads were provided free from barcodes and adapters sequences.

Sequencing data from msmA amplicons, SCD0-A and SCDE-A sets, and msmE amplicons, SCD0-E set, were submitted to the European Nucleotide Archive (http://www.ebi.ac.uk/ena) under sample accession numbers ERS954926, ERS954925 and ERS954927, respectively, within the metagenomic study PRJEB9018.

The overall quality scores of the output data from the Ion Torrent sequencing was assessed using FastQC on Galaxy web-based platform. Sequencing reads were then trimmed and filtered based on length and quality. For the msmA amplicon, reads ranging from 300 to 365 bp were selected; for the msmE amplicon, sequences with length between 335 to 420 bp were retained. In the quality filtering step, sequences with quality lower than 20 over more than 15% of the nucleotides were discarded. Chimeric sequences were removed with Chimera Check (Edgar et al., 2011) available in the FunGenePipeline (Fish et al., 2013). The resulting files were submitted to FrameBot (Wang et al., 2013) for translation and frameshift correction (FunGenePipeline). Sequences containing stop and/or undefined codons were discarded (general statistics are reported in Table S5). For each dataset, protein sequences were aligned with hmmalign (HMMER3 (Eddy, 2011)) and clustered by complete linkage clustering (mcClust (Loewenstein et al., 2008), FunGenePipeline). Rarefaction (FunGenePipeline (Fish et al., 2013)) was based on data from the clust file at 0.03 distance. Conservation analysis of the MsmA and MsmE predicted protein sequences was performed on alignments generated by hmmalign (HMMER3) and analyzed by Jalview (Waterhouse et al., 2009). Alpha diversity analysis (Shannon index, evenness and Chao species estimator) was performed on clustering values obtained at 0.03 distance for each set of results. Beta diversity (Bray-Curtis, Jaccard, Kulczynski, and Chao indices) between SCD0-A and SCDE-A sequence sets was computed on the clustered (0.03 distance) aligned MsmA predicted sequences. For Bray-Curtis, relative cluster abundances were employed. The representative sequence for each cluster was obtained using Representative Sequences (FunGenePipeline). In this case, distance cutoff values were chosen in order to obtain in the region of 20–30 clusters from each dataset (0.15 for SCD0-A/SCDE-A and 0.10 for SCD0-E). The nucleotide sequences corresponding to these cluster-representative sequences were retrieved and used to infer phylogenetic trees by Maximum Likelihood with 100 bootstrap iterations (as implemented in PhyML (Guindon et al., 2010) at Phylogeny.fr (Dereeper et al., 2008)). A flowchart with the major steps of this analysis is available in Fig. S2.

Characterization of the samples

The seawater sample collected for this study was at a temperature of 15.5°C and 3.34% salinity (conductivity 47.9 mS/cm). The pH value was 8.07 and the concentration of chlorophyll a was 1.06 mg/m³.

Optical density (at 600 nm) of the filtered sample was 0.606 at time 0 (SCD0), dropped to 0.415 at day 2 and rose towards the last days of incubation to 0.758 (SCDE). During the enrichment process the suspension progressively lost its original greenish hue.

Analysis of phylogenetic composition

Despite the fact that sample SCDE underwent just a partial enrichment with no dilution or subculturing, this was enough for the G + C content of the community’s metagenome to shift from 43.30% (SCD0) to 53.29% (SCDE).

Phylogenetic composition analysis was performed by EBI Metagenomics on the metagenomic data, based on prokaryotic rRNA gene sequences. Rarefaction curves were generated from these results for both metagenomes (Fig. S3). These curves show that a reasonable coverage of the communities’ composition was achieved: however, neither reaches a definite plateau meaning that deeper sequencing levels would have been desirable.

Phylogenetically, the community before the enrichment was composed primarily by taxa typical of surface seawater, namely Alphaproteobacteria of the Pelagibacteraceae and the Rhodobacteraceae, Gammaproteobacteria of the Halomonadaceae and Piscirickettsiaceae, Archaea of the Crenarchaeota (with prominence for Nitrosopumilus) and the Euryarchaeota (Thermoplasmata Marine group II), Bacteroidetes of the Flavobacteriaceae and the uncultured marine phylum SAR406 (Figs. 1 and 2).

Clearly, some of these groups declined or disappeared during our enrichment program (Figs. 1, 2 and 3): all the Archaea and SAR406, the Pelagibacteraceae, Halomonadaceae, Flavobacteriaceae (a family known to harbor methylotrophic species (De Marco et al., 2004; Moosvi et al., 2005a; Boden et al., 2008; Madhaiyan et al., 2010)) and Variovorax paradoxus (Comamonadaceae) a known methylotroph (Anesti et al., 2005). Other groups thrived: the Rhodobacteraceae and the Hyphomonadaceae within the Alphaproteobacteria, the Oceanospirillaceae, the Piscirickettsiaceae, two families of Alteromonadales and genus Alcanivorax within the Gammaproteobacteria, the Saprospiraceae within the Bacteroidetes. It is worth noticing that: fam. Rhodobacteraceae includes several methylotrophic species and namely marine MSA-degrader Marinosulfonomonas (Thompson, Owens & Murrell, 1995); fam. Piscirickettsiaceae includes the genus Methylophaga, a known methylotrophic taxon relevant in surface sea waters (Neufeld et al., 2008). Indeed, the signal for Methylophaga itself at the genus level shot up from below detection to 82 hits; genus Alcanivorax, which contains marine representatives that have previously been found in a dimethylsulfide-based enrichment (Schäfer, 2007) and one known potentially methylotrophic species (A. borkumensis strain SK2—http://www.ebi.ac.uk/biomodels-main/BMID000000083743). Despite the fact that to our knowledge no methylotrophic species is known among the Hyphomonadaceae, the considerable increase in its abundance after the enrichment (225x) seem to suggest a possible involvement of members of this group in the turnover of MSA. Other rarer taxa that showed sizeable abundance increases after the enrichment were the Cryomorphaceae (Bacteroidetes), among the Gammaproteobacteria several branches within the Alteromonadales, several groups within the Oceanospirillales like the Oceanospirillaceae (and species Neptuniibacter caesariensis), the Oleiphilaceae, genus Halomonas, and rather surprisingly the Enterobacteriaceae. Notably, known methylotrophic-containing taxa like the Hyphomicrobiaceae (Alphaproteobacteria) and Methylotenera mobilis (Betaproteobacteria) also increased significantly.

Metagenomic data binning by MetaBat failed to generate discrete bins with sample SCD0: this is not surprising with a diverse natural community under no definite selective pressure where no single species dominate. However, three bins were obtained with sample SCDE which by and large matched the results obtained by direct metagenomics analysis. Bin 1, assigned to an Alcanivorax sp. strain, contains genes coding for diagnostic enzymes for the assimilation of C₁ carbon through the serine cycle (serine hydroxymethyltransferase, serine-glyoxylate aminotransferase and isocitrate lyase). Bin 2, tentatively classified as the genome of a Hyphomonas sp. strain, also contains genes for the serine cycle (serine hydroxymethyltransferase) and for the anaplerotic ethylmalonyl-CoA and glyoxylate pathways (propionyl-CoA carboxylase, ethylmalonyl-CoA mutase and isocitrate lyase). Bin 3 showed phylogenetic markers mostly associated with genus Methylophaga: it contains genes diagnostic for the ribulose monophosphate (RuMP) cycle (3-hexulose-6-phosphate synthase and 6-phospho-3-hexuloisomerase), an alternative metabolic route for the assimilation of C₁ carbon, which is the typical metabolic route found in Methylophaga species.

Alpha and beta diversity estimates were based on the phylogenetic data obtained by the EBI Metagenomics analysis. The values are listed in Tables 1 and 2. The effects of the enrichment were an increase in the number of taxa observed, internal diversity (H’), and evenness. This was due to the sharp decline or disappearance of some of the major taxa composing the original metagenome (all Archaeal groups, algal chloroplast sequences, the Pelagibacteraceae, and the SAR406 lineages) accompanied by the appearance of a large number of less representative taxa, which shows that the enrichment process was halted at a very early stage, as intended.

Overall analysis of the functional annotation

Significant shifts between the functional annotation of samples SCD0 and SCDE obtained by EBI Metagenomics were observed.

The analysis revealed somewhat surprising large increases in categories GO:0000103 “sulfate assimilation” and GO:0008272 “sulfate transport.” Among the other categories, significantly rising in hit numbers were GO:0016846 “carbon-sulfur lyase activity,” GO:0016705 “oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen” and GO:0004497 “monooxygenase activity,” GO:0051537 “2 iron, 2 sulfur cluster binding,” GO:0006730 “one-carbon metabolic process,” GO:0004488 “methylenetetrahydrofolate dehydrogenase (NADP⁺)” and GO:0008864 “formyltetrahydrofolate deformylase,” all of which can be associated to the catabolism of MSA. The latter two may suggest that a selection for organisms that metabolize formaldehyde through the condensation with tetrahydrofolate (rather than with tetrahydromethanopterin) may have occurred. Oddly, though, GO:0046653 “tetrahydrofolate metabolic process” and GO:0046654 “tetrahydrofolate biosynthetic process” were found among the significantly decreased categories. Fittingly, several category associated with oxygenic photosynthesis and GO:0004329 “formate-tetrahydrofolate ligase activity” (incorporation of formate into biomass in anaerobic conditions) also decreased in abundance.

Among the monooxygenase genes significantly increased in number during the enrichment, 5 hits for gene msuD/ssuD were found in the assembled metagenomic data of the enriched sample (SCDE), while no homolog was present in sample SCD0.

Since no hits for msm genes were observed in the metagenomic results by either EBI Metagenomics or IMG/MER, the sequence data were also screened locally by hmmsearch (HMMER3) and tblastn (McGinnis & Madden, 2004). The only hits obtained with protein MsmA (expectation value cutoff 10⁻⁵) contained short Rieske-associated motifs (spacers varying from 16 to 18 in SCD0, and from 16 to 23 in SCDE) and were low scoring. As for protein MsmE, no hits were found at E-value ≤10⁻⁵.

In order to complement the metagenomic data, msmA and msmE genes were investigated by an amplicon survey experiment. Direct amplification was weak or absent, so a nested PCR approach was adopted for both genes, which witnesses to the overall low concentration of these sequences in the metagenomes. Amplification of gene A was achieved from both samples. On the contrary, despite repeated and insistent efforts, we were unable to amplify gene E from the enrichment sample (SCDE). For convenience, the three amplicon sets will be designated as follows: SCD0-A (msmA amplicons from sample SCD0), SCDE-A (msmA amplicons from sample SCDE) and SCD0-E (msmE amplicons from sample SCD0). The three amplicon pools were sequenced by ion Torrent technology. The G + C content of the processed reads was close to 40% in all 3 cases (Table S5). Namely, while community metagenome grew from 43.30% to 53.29% G + C content, there was almost no difference in base composition between of the msmA reads before and after the enrichment. The reads were then translated and checked for frameshift-generating sequencing errors and all sequences containing an undetermined codon or a stop codon were eliminated. Statistics relative to each case are shown in Table S5.

Rarefaction curves constructed from clustering data at a 0.03 distance are represented in Fig. 4. Although a definite plateau is not achieved, gene A looks sampled at satisfactory depth in both SCD0-A and SCDE-A sets. On the contrary, the level of amplification and sequencing of gene E in set SCD0-E was clearly insufficient. These results suggest that a deeper sequencing would be needed. However, it is also true that a huge part of the data was deleted at the quality control step (Table S5) so better sequencing quality rather than deeper sequencing might be the solution.

The data from clustering of MsmA and MsmE sequences at 0.03 cutoff were analyzed for alpha and beta diversity. Clearly there was an evolution in MsmA sequences during the enrichment process which lead to the loss of some sequence diversity (Table 3). The shift in MsmA allele composition is very significant as illustrated by the high values of the beta diversity indexes (Bray Curtis = 0.95; Jaccard = 0.98; Kulczynski = 0.95; Chao index = 0.76).

Conservation analysis was performed with the predicted MsmA and MsmE sequences (Figs. 5 and 6). The edges of the sequences, corresponding to the PCR primers, should be 100% conserved. The low conservation levels seen in these regions are artifacts due to the occurrence of a few shorter sequences in the datasets and lower quality at the end of reaction. In the case of MsmA (Fig. 5), the levels of amino acid conservation are generally higher and more constant along the central region of the fragment. If the primers regions are not taken into account, all the conservation values are ≥92.35%, with the exception of a one-amino acid trough (a Serine at position 30) with 58% conservation. Particularly, all the amino acids in the Rieske associated motif and spacer demonstrate conservation levels higher than 99%. Regarding MsmE (Fig. 6), it is possible to observe regions with high levels of conservation, which in some cases are greater than 99%. However, the overall conservation is much less constant than for MsmA, with 4 positions below 85%. The lowest conservation value (61%) was found for a Serine at position 10.

Phylogenetic trees of the predicted MsmA and MsmE sequences

Phylogenetic trees were constructed based on the nucleotide data corresponding to the representative protein sequences of clusters obtained with distance cutoff values 0.15 (for MsmA) or 0.10 (for MsmE). In the process, sequences from this study were aligned with metagenomic seawater homologs obtained in a previous work (Henriques & De Marco, 2015a), with sequences EF103447, EF103448 and EF103448 from SSM clones (Leitão, Moradas-Ferreira & De Marco, 2009), with sequences from the GOS project (Rusch et al., 2007) and with sequences from cultured strains. Phylograms of msmA and msmE sequences are shown in Figs. 7 and 8. Although some degree of caution is required in the analysis of these sequences due to the fact that they were produced by three rounds of amplification (whole-genome amplification followed by nested PCR), the results obtained seem to fit logical expectations. In the case of gene msmA, it is possible to observe a clear separation (95% bootstrap value) between two principal branches: one consists only of metagenomic sequences and includes all the msmA reads generated in this study; the other group (omitted as root in Fig. 7) contains the sequences from cultured strains (Alpha, Beta or Gammaproteobacteria) of both marine and soil origin. In a similar way, the phylogram for gene msmE shows an unambiguous split between two major groups (94% bootstrap value): one comprising mainly metagenomic sequences, including all the msmE reads from this study, sequences from the GOS project, and two from cultured marine strains C. Puniceispirillum marinum IMCC1322 and C. Filomicrobium marinum str. Y; plus a second group containing all other sequences from cultured strains (Alpha and Betaproteobacteria).