Metagenomes of polyamine-transforming bacterioplankton along a nearshore–open ocean transect

Dissolved organic nitrogen (DON) is a well-recognized pool of labile N for marine microorganisms (Sipler and Bronk 2015). A large number of studies on bacteria-mediated DON transformation have been focused on amino acids and urea. However, increasing attention has been given to a group of understudied DON compounds [i.e., polyamines (PAs)] (Liu et al. 2015; Lu et al. 2014, 2020), which are produced by organisms of all three life domains (Tabor and Tabor 1984) and mediate a number of essential cellular functions (Igarashi and Kashiwagi 1999, 2000). PAs share many biogeochemical features with amino acids, such as C: N ratio, wide distribution, and low concentrations (nmol/L level) in seawater (Lee and Jørgensen 1995; Lu et al. 2014). Putrescine is thought to be one of the most abundant PAs in aquatic environments (Höfle 1984). Using putrescine, spermidine, and spermine as models, PAs were found to fulfill nearly 10% of the N demand of marine bacterioplankton (Liu et al. 2015). Studies have also found that PA concentrations and turnover can be highly dynamic, observed daily and varying seasonally and spatially in marine environments (Liu et al. 2015; Lu et al. 2014); these concentrations and turnover rates were sometimes comparable to dissolved free amino acids (Liu et al. 2015; Lu et al. 2014).

In addition to direct measurements, the importance of PAs to aquatic microbial communities is also evident by the wide distribution of genes involved in PA transport through polyamine ABC transporters (Igarashi and Kashiwagi 2010) and degradation via γ-glutamylation, transamination, or spermidine cleavage pathways in bacterial genomes (Mou et al. 2010), metagenomes (Mou et al. 2011), metatranscriptomes (Lu et al. 2020), and metaproteomes (Sowell et al. 2008). Most of the field studies on PAs have been performed in nutrient-rich coastal seawaters (Lu et al. 2015; Mou et al. 2011; Nishibori et al. 2001, 2003); many of these studies have consistently suggested that a diverse group of marine bacteria are involved in PA transformation (Mou et al. 2011, 2014), with roseobacter and SAR11 of Alphaproteobacteria as the most dominant taxa (Lu et al. 2015; Mou et al. 2011). These studies have also suggested that among the three known PA transformation pathways, transamination outweighs γ-glutamylation and spermidine cleavage during bacteria-mediated PA degradation in coastal seawater (Lu et al. 2015; Mou et al. 2011). However, studies that examined PA transformation in oligotrophic areas are scarce and even fewer have assessed the potential differences among transformations of individual PA compounds (Lu et al. 2015, 2020).

In the present study, metagenomes that were responsive to additions of PA compounds were compared among bacterioplankton collected along a nearshore-offshore-open ocean transect that covers a gradient of nutrient loads in the Gulf of Mexico. In addition to putrescine, another two recently recognized common marine polyamines, i.e., spermidine or spermine were used. Therefore, in addition to examining whether PA-transforming genes and taxa may vary along a nutrient gradient, our study also allows the comparison of PA transformations among different PA compounds.

Surface water samples were collected at three stations in the Gulf of Mexico along a transect, i.e., one nearshore (NS), one offshore (OS), and one open ocean (OO) station (Fig. 1). The measured in situ environmental variables varied among NS, OS, and OO (Lu et al. 2020; Supplementary Table S1). From the NS site that is distant from the shore (OS and OO), the salinity increased, while nitrate (NO₃^-), nitrite (NO₂^-), dissolved organic carbon (DOC), dissolved nitrogen (DN), and water temperature (Temp) decreased (one-way ANOVA, P < 0.05). Opposite from the other nutrients, concentrations of ammonium (NH₄⁺) increased and had the highest values at OS (one-way ANOVA, P < 0.05). HPLC was used to detect putrescine (PUT), spermidine (SPD), and spermine (SPM) (Lu et al. 2014; Fig. 2A). SPD (3.9–5.6 nmol/L) and SPM (4.4–5.5 nmol/L) were measured at all three sites with no significant spatial difference (one-way ANOVA, P > 0.05; Fig. 2A). PUT was not detected at any site in the present study. The individual and total concentrations of PAs here were similar to previously reported values from marine environments (Lu et al. 2014; Nishibori et al. 2001, 2003) but appeared to be less diverse.

Figure 1.  The sampling sites of NS, OS, and OO in the Gulf of Mexico in May, 2013. The depth of water column at each site is listed in the parentheses

Figure 2.  A The concentrations of individual PA compounds in situ environment, and B the degradation rate of individual PA compounds in microcosm experiment. The error bar represents standard error

Duplicate microcosms were set up on the deck of the sampling cruise ship using NS, OS, and OO samples and amended with 200 nmol/L (final concentration) of PUT, SPD, SPM, or an equal volume of diH₂O (serving as controls, CTL). After 48 h of incubation at ambient temperature and light, residual PAs in each of the microcosms were measured and used to calculate PA degradation rates (Fig. 2B). Significant PA degradation (50–100%) was observed in all microcosms. Comparable rates of PA degradation by marine bacteria have been observed in the Atlantic Ocean (Höfle 1984; Lu et al. 2015). No significant difference of PA degradation rate was observed among microcosms of different sites and PA amendments (Fig. 2B).

Metagenomes of samples from each of the incubated microcosms (48 h) were sequenced, which yielded a total of 6, 700, 391 Illumina Mi-Seq sequences with an average length of 363 bp (Table 1). Samples for PUT metagenomes in the OO sites were lost during transportation and were not included in the analysis. Ribosomal RNA (rRNA) sequences accounted for 0.7–1.5% of the metagenomic sequences, which is typical for metagenomes (Mou et al. 2008). A total of 6, 564, 670 sequences were identified as putative protein-coding sequences, 62.0% of them further received annotations to the gene level. The gene sequences were then functionally classified into 1742–2289 unique clusters of orthologous groups (COGs, 26.8–50.8%), and 150–184 unique Kyoto encyclopedia of genes and genomes (KEGG) pathways (19.2–33.6%).

NMDS and ANOSIM analyses were performed based on the relative abundances of major (≥97%) COGs among PUT, SPD, and SPM metagenomes. They consistently showed that metagenomes of NS, OS, and OO bacterioplankton were well separated (r_ANOSIM=0.65, P < 0.05; Fig. 3). The distinct COG contents between NS and OS/OO sites suggest that there may be a variety of metabolic strategies among sites with different environmental conditions (Lu et al. 2020).

Figure 3.  The non-metric multidimensional scaling (NMDS) plot based on the relative abundance of major COGs in metagenomes of nearshore (NS; triangle), offshore (OS; square), and open ocean (OO; circle) in the Gulf of Mexico

Xipe-TOTEC analysis identified 129–197 COGs that were significantly enriched in the PUT metagenomic libraries relative to the CTR libraries in the NS site. The number of COGs increased to 249–296 in OS and then to 289–357 in the OO site. A total of 12 of the COGs have been shown to be involved in polyamine transformation steps, such as PA synthesis (3 COGs), transport (6 COGs), and degradation (3 COGs) (Fig. 4; Supplementary Table S3). The distribution of these PA-related COGs was site- and PA compound-specific. Specifically, the OS site had the highest number (up to 10 in total) of PA-related COGs among the three sites. In line with overall COG data, this again indicates that PA compounds may be metabolized via different pathways (Chou et al. 2008; Dasu et al. 2006; Mou et al. 2010) by marine bacterioplankton in a given location. In addition, most of the PA-related COGs were found in SPM (7 COGs) metagenomes (Fig. 4), indicating the importance and universal distribution of SPM along the studied transect.

Figure 4.  Fold changes of previously identified PA-transforming COGs in PA vs. CTR metagenomes based on Xipe-TOTEC analysis. AA amino acid, BioA adenosylmethionine-8-amino-7-oxononanoate aminotransferase

The taxonomic affiliations of all enriched COGs, i.e., PA-responsive COGs, in metagenomic libraries, were significantly different among sampling sites and PA compounds (i.e., PUT, SPM, and SPD) at the family level, as revealed by NMDS (data not shown) and ANOSIM (r_ANOSIM=0.87, P < 0.05) analyses. In the NS metagenomes, Rhodobacteraceae (Alphaproteobacteria) was generally dominating in PUT (14.4% in average), SPD (15.0%), and SPM (21.6%) libraries (Fig. 5). Rhodobacteraceae-affiliated roseobacters are known for their high abundance and strong ability in processing plankton-derived DOC compounds (González et al. 2000; Hahnke et al. 2013). These results are consistent with a previous metatranscriptomic study, suggesting the importance of roseobacters in PA transformations in nearshore marine environments (Lu et al. 2020).

Figure 5.  Taxonomic binning of the protein-encoding sequences in significantly enriched COGs at bacterial family level in the PA amended metagenomes (PUT, SPD, and SPM). Rows were organized using the hclust method in R

In the OS metagenomes, PA-responsive COGs were most affiliated with three Gammaproteobacteria families, i.e., Alteromonadaceae (14.6%), Pseudomonadaceae (13.6%, ), and Alcanivoracaceae (13.2%) in the PUT library (Fig. 5). In the SPD library, unclassified Chroococcales (12.4%; Cyanobacteria) and Planctomycetaceae (7.6%; Planctomycetes) were the most abundant PA-responsive families (Fig. 5). In SPM, Prochlorococcaceae (13.0%; Cyanobacteria), Rhodobacteraceae (9.2%, Alphaproteobacteria), and Comamonadaceae (9.0%; Betaproteobacteria) were predominant (Fig. 5).

In the OO metagenomes, PA-responsive COGs were again mostly affiliated with Gammaproteobacteria but with different families from those in the OS metagenomes; these included Idiomarinaceae (13.7%, ) and Shewanellaceae (13.0%) in the SPD library, and Pseudoalteromonadaceae (55.4%, ) in the SPM library (Fig. 5). The importance of various Gammaproteobacteria in transforming PAs in non-coastal sites has also been suggested in a study performed at the South Atlantic Bight (Lu et al. 2015, 2020). No PA-responsive COGs were found to be affiliated with SAR11, although this taxon dominates open ocean environments and has been suggested to be involved in PA transformations in a metaproteomic study (Sowell et al. 2008).

Variations in the relative abundance of responsive bacterioplankton to different PA compounds (ANOVA, P < 0.05) were observed in this study; this agrees with the findings of a previous study, which identified different PA-responding bacterial families between PUT and SPD transformation (Lu et al. 2015). However, a metatranscriptomic study performed at the same Gulf of Mexico sites indicates that PUT, SPM, and SPD may share similar taxa for their transformation (Lu et al. 2020). The metatranscriptomic study also identified that PA-responsive taxa varied between the OS and OO sites, which is also identified in the present study. This discrepancy may partly be explained by method differences, such as length of incubation time (days for the metagenomic study vs. hours for the metatranscriptomic study). More importantly, the two studies sequenced different nucleic acids, i.e., community DNAs for the metagenomic study and community mRNA for the metatranscriptomic study (Aguiar-Pulido et al. 2016; Shakya et al. 2019). Response of mRNA synthesis to external stimuli can be observed in seconds to minutes, while changes of DNA require a much longer time (days) (Aguiar-Pulido et al. 2016). However, mRNA is far less stable than DNA and some mRNA can be degraded before sample collection (Shakya et al. 2019). Moreover, different bacteria may have different sensitivity to PA compounds, some can respond in seconds, and others may take longer times.

Genes that are known to participate in processing PAs, including transporter genes (potABCDEFGH), γ-glutamylation genes (puuABCDE), transamination genes (spuC, kauB, and GabT), and SPD cleavage genes (spdH and gltA) were identified, and their relative abundance was calculated for each metagenomes (Supplementary Table S3). Among PA uptake- and degradation-related genes, the ones that were significantly enriched in PA metagenomes relative to the corresponding CTRs were designated as PA-responsive.

Responsive PA transporter genes were only identified in SPM libraries of OS and OO (OR > 1, P < 0.02; Fig. 6A). The taxonomic affiliations of these transporter genes were primarily assigned to Rhodobacteraceae (26.9% of the total assigned putative PA genes) and Alteromonadaceae (66.7%), respectively (Fig. 6B).

Figure 6.  A The enriched PA degradation pathways; B the relative abundance of diagnostic PA uptake/degradation genes and their taxonomic affiliations in CTR, PUT, SPD, and SPM

Responsive PA degradation genes, however, were identified in all three PA metagenomes (PUT, SPD, and SPM) of all three sites (OR > 1, P < 0.02). At the NS site, putative γ-glutamylation genes were enriched in SPD metagenomes (OR > 1, P < 0.02; Fig. 6B), and were mostly affiliated with Rhodobacteraceae (6.0%) (Fig. 6B). In contrast, putative SPD cleavage genes were enriched in the PUT and SPM metagenomes of NS (OR > 1, P < 0.02; Fig. 6B), and the taxonomic binnings of these genes were primarily assigned to Methylophilaceae (2.8%; Betaproteobacteria) and SAR11 clade (3.3%; Alphaproteobacteria), respectively (Fig. 6B).

At the OS, putative γ-glutamylation genes were enriched in PUT and SPD metagenomes (OR > 1, P < 0.02; Fig. 6B), and were primarily binned to Alteromonadaceae (1.0%) and Plantomycetaceae (1.6%), respectively (Fig. 6B). Putative SPD cleavage genes also showed enrichment in SPD metagenomes (OR > 1, P < 0.02; Fig. 6B) and they were mainly assigned to Planctomycetaceae (1.8%) and Alteromonadaceae (1.4%) (Fig. 6B).

At the OO, putative γ-glutamylation genes were enriched in SPD metagenomes (OR > 1, P < 0.02; Fig. 6B) and the majority of the sequences were affiliated with OMG group (2.0%) (Fig. 6B). Putative SPD cleavage genes were enriched in SPM metagenomes of OO (OR > 1, P < 0.02; Fig. 6B) and were taxonomically binned to Rhizobiaceae (1.8%) and Shewanellaceae (1.8%) (Fig. 6B). Putative transamination genes were enriched also in SPD metagenomes (OR > 1, P < 0.02; Fig. 6B) and were mostly affiliated with Alteromonadaceae (3.2%) (Fig. 6B). No genes of the PA transamination pathway were identified as PA-responsive in either NS site or OS site.

Results of PA-responsive gene analysis largely agree with the COG gene results on the diversity of bacterial taxa that are potentially involved in PA transformation (Fig. 5), including families of Proteobacteria (alpha, beta, and gamma-lineages), Planctomycetes and Cyanobacteria. At the family level, Rhodobacteraceae of Alphaproteobacteria consistently appeared to be the most important in PA transformation at all sites and for all tested PA compounds. In addition, these two sets of results both indicate that PA-responsive bacterioplankton taxa of individual PA compounds can be site-specific. The PA-responsive gene analysis, recovered genes of SAR11 (Alphaproteobacteria) from most PA metagenomes and these genes were mostly affiliated with transporter genes. The indicated importance of SAR11 in PA transformation is consistent with the previous studies (Lu et al. 2020; Sowell et al. 2008).

The absence and relatively low abundance of PA transamination genes in the PA-responsive gene groups at the NS and OS sites suggests a minor importance of this pathway in PA transformation in the coastal regions of the Gulf of Mexico. This result contrasts with the finding of a previous metatranscriptomic study of coastal PA-transforming bacterioplankton in the South Atlantic Bight, in which transamination dominated the PUT and SPD degradation (Mou et al. 2011). This reveals large-scale spatial differences in PA transformation.

Similar to the COG-gene-based analysis, results of PA-responsive gene analysis of metagenomes had some different findings from the previous metatranscriptomic study (Lu et al. 2020). For example, in this study, PA transporter genes were responsive to the PA addition in the OS and OO sites; however, they were only responsive in the NS site in the metatranscriptomic study. In addition, in the NS site, this study found that γ-glutamylation was most important in transforming PAs. However, the metatranscriptomic study identified the transamination pathway as the most important at the same site. We again attributed these different results mostly to method differences described in the above section. Nonetheless, most results of this metagenomic and previous metatranscriptomic studies were consistent (Lu et al. 2020) and allow a better understanding of PA transformations of individual compounds at coastal and open ocean environments.

Using metagenomic approaches, we identified variations in PA-transforming bacterioplankton genes and taxa among different marine systems in the Gulf of Mexico. Transamination, γ-glutamylation, and SPD cleavage all participated in PA degradation but SPD cleavage might be more important in all tested sites. PA-transforming bacteria covered a diverse group of bacteria from phyla of Actinobacteria, Bacteroidetes, Cyanobacteria, Frmicutes, Planctomycetes, and Proteobacteria. At the nearshore site, Rhodobacteraceae of Alphaproteobacteria played a key role in driving PA transformation, while at offshore and open ocean sites, bacterial families of Gammaproteobacteria were the predominant PA-transforming bacterial taxa. Although PA-transforming bacterioplankton of individual PAs can be site-specific, they showed little variations among PA metagenomes of the same site, suggesting that PA-transforming taxa are capable of processing multiple PA compounds.

Surface water (0–2 m) was collected onboard the R/V Pelican using 12 L Niskin bottles mounted on a rosette sampling system (Sea-Bird Electronics, Bellevue, WA) in May 2013 in the Gulf of Mexico. Sampling was performed at three stations, i.e., one nearshore (NS), one offshore (OS), and one open ocean (OO) station (Fig. 1). In situ temperature (T) and salinity (S) of the sampling sites were measured by a conductivity–temperature–depth (CTD) water column profiler (Sea-Bird Electronics, Bellevue, WA, USA) at the same time as sampling.

After collection, part of water samples was immediately filtered through 3 μm pore-size membrane filters (EMD Millipore Corp., Billerica, MA, USA) to remove larger organisms and debris. Eight carboys, each containing 18.9 L water filtrates, were set up as microcosms; each received additions of either 200 nmol/L (final concentration) PUT, SPD, SPM, or an equal volume of diH₂O (serving as controls, CTL). Microcosms were incubated on board at ambient air temperature and light for 48 h. Samples (500 ml) were taken from the microcosm at the beginning and end of the experiment and filtered through 0.2 μm pore-size polycarbonate membrane filters (Pall life sciences, Ann Arbor, MI), and the resulting filtrates were immediately frozen at - 20 ℃ onboard and transported on dry ice back to the lab, where they were kept at - 80 ℃ before being used for nutrient analysis and PA concentration measurements.

After 48 h of incubation, bacterial cells in microcosms were collected onto 0.2 μm pore-size isopore membrane filters (EMD Millipore Corp., Billerica, MA, USA) by filtration. The total filtering time for each sample was kept within 30 min. Filters were immediately stored in a liquid nitrogen and transported back to the lab, where they were transferred into a-80 ℃ freezer until DNA extraction.

All plastic bottles were acid-washed and all glassware was combusted at 500 ℃ for at least 6 h and rinsed with diH₂O prior to use.

Concentrations of dissolved organic carbon (DOC), dissolved nitrogen (DN), nitrate (NO₃^-), nitrite (NO₂^-), soluble reactive phosphorus (SRP), ammonium (NH₄⁺), and individual PAs (PUT, SPD, and SPM) were measured following a procedure described previously (Lu et al. 2020). Briefly, DOC and DN were measured using a TOC/TN analyzer (TOC-V_CPN; Shimadzu Corp., Tokyo, Japan); NO₃^- (Jones 1984) and NH₄⁺ (Strickland and Parsons 1968) and SPR (Murphy and Riley 1962) were measured spectrophotometrically; NO₂^- were measured based on a microplate technique (Hernández-López and Vargas-Albores 2003). Polyamine concentrations were measured using a Shimadzu 20A high-performance liquid chromatography (HPLC; Shimadzu Corp., Tokyo, Japan) that was equipped with a 250×4.6 mm i.d. 5 µm particle size, Phenomenex Gemini-NX C₁₈ column (Phenomenex, Torrance, CA, USA). Samples were derivatized with o-phthaldialdehyde, ethanethiol, and 9-fluorenylmethyl chloroformate before running on the HPLC machine (Lu et al. 2014).

Genomic DNA was extracted from 0.2 μm pore-size membrane filters using the Qiagen DNeasy DNA extraction kits (Qiagen, Chatsworth, CA, USA). An additional step of bead beating with 0.1 mm size glass beads (0.2 g/filter; BioSpec, Bartlesville, OK, USA) for 10 min at 3000 r/min was added after enzymatic lysis with lysozyme and proteinase K during DNA extraction (Hunt et al. 2013). The quantity of DNA was determined with the Quant-iT PicoGreen ds DNA Assay Kits (Life technologies, Carlsbad, NY, USA). DNA extracts of treatment replicates were pooled before sequencing. DNAs of the PUT microcosms at the OO site were lost during processing. DNA library of each treatment sample was prepared with TruSeq Nano DNA Sample Prep Kits and sequenced using the Illumina MiSeq 2×250 pair-end technology (Illumina Inc., San Diego, CA, USA) at the University of Minnesota Genomics Center.

The raw DNA sequences were deposited in the Sequence Read Archive of NCBI under accession no. SRP049693.

Quality control and annotations of Illumina sequences (paired-end) were performed in the Metagenome Rapid Annotation using Subsystem Technology (MG-RAST) v3 (Meyer et al. 2008). Putative protein-coding sequences were identified and annotated using the sBLAST analysis against a non-redundant protein database (M5NR), which integrates a number of databases such as SEED, KEGG, and RefSeq. Taxonomic affiliation of all the putative protein-coding sequences was performed based on cut-off of E value≤10^–20, percent identity≥40%, and alignment length≥69 bp, which is approximately corresponding to bit score≥40 (Mou et al. 2008). Homologs to empirically identified genes that are involved in polyamine transport and degradations (Supplementary Table S2; Chattopadhyay et al. 2009; Chou et al. 2008; Dasu et al. 2006; Lu et al. 2002; Mou et al. 2010) were putatively identified by tBLASTn with a cut-off value of bit score≥40 (Mou et al. 2008).

A non-metric multidimensional scaling (NMDS) analysis was performed to ordinate CTR, PUT, SPD, and SPM metagenomic libraries using PRIMER-v5 (Plymouth Marine Laboratory, Plymouth, UK; Clarke and Warwick 2001). The similarity matrix was calculated based on normalized and square-root-transformed relative abundances of major COGs (collectively accounted for≥97% of total sequences) using the Bray–Curtis algorithm. The robustness of NMDS grouping patterns was statistically evaluated by ANOSIM (analysis of similarity), which is an analogue of the standard univariate ANOVA (analysis of variance). The ANOSIM index r_ANOSIM was calculated on a scale of 0–1. When P < 0.05, the sample groups were identified as well separated when r_ANOSIM > 0.75, clearly different but overlapping when 0.5 < r_ANOSIM≤0.75, or barely separable when r_ANOSIM < 0.25 (Clarke and Warwick 2001).

Pair-wise comparisons were performed to compare the gene contents between PA amended (PUT, SPD, or SPM) and CTR metagenomes using the Xipe-TOTEC program (Rodriguez-Brito et al. 2006). Resampling (20, 000 times) was undertaken for each comparison based on a sample size equal to the smaller number of sequences in the two metagenomic libraries. Significant differences were reported when P < 0.01.

Xipe-TOTEC analysis tests proved to be insensitive to variations in sequences of low abundance (Mou et al. 2008). To address this potential bias, odds ratios (OR) and bionomical distribution probabilities of variations in the relative abundances of the putative PA uptake and degradation genes were also calculated (Gill et al. 2006). The equation for OR calculation is [n_p/(N_p - n_p)]/[n_c/(N_c - n_c)], where n_p is the number of target gene sequences in the polyamine (PUT, SPD, or SPM) metagenomes, n_c is the number of target gene sequences in the controls (CTR) metagenomes; N_p and N_c are the total number of sequences in the polyamine and control metagenomes, respectively. The binomial distribution was presumed in each library and probability (P) was calculated with the [n_c/(N_c - n_c)] as the expected sequence frequency. False discovery rates (FDR) of P_binomial values were calculated using the Benjamini-Hochberg method (1995). PA gene was considered as significantly enriched when the corresponding OR > 1 and FDR < 5%.

Environmental variables and bacterial taxa of enriched COGs and PA diagnostic genes were compared between or among samples using t test or ANOVA in the R software package (R Core Development Team 2005). Significant differences were reported when P < 0.05.

The online version contains supplementary material available at https://doi.org/10.1007/s42995-021-00114-x.

We appreciate the help of stuff working on R/V Pelican for assisting sample collection and providing CTD data, and ZF Liu for measuring the in situ concentrations of nitrate, ammonium, and soluble reactive phosphorus. This study was supported by National Science Foundation Grants (OCE1029607 to X.M) and Kent State University.

XL performed experiments and bioinformatic analysis, and interpreted the data; KW performed bioinformatic analysis and interpreted the data; XM designed the experiment and interpreted the data. XL, KW, and XM wrote the manuscript.

The authors declare no conflict of interest.

No animal and human right are involved in this research.