Metabolites reflect variability introduced by mesoscale eddies in the North Pacific Subtropical Gyre

Mesoscale eddies significantly alter open ocean environments such as those found in the subtropical gyres that cover a large fraction of the global ocean. Previous studies have explored eddy effects on biogeochemistry and microbial community composition but not on the molecular composition of particulate organic matter. This study reports the absolute concentration of 67 metabolites and relative abundances for 640 molecular features, measured using liquid chromatography-mass spectrometry (LC-MS) following both targeted and untargeted approaches. This approach allowed us to better understand how mesoscale eddies impact the metabolome of the North Pacific Subtropical Gyre during two cruises in 2017 and 2018. We find that many metabolites track biomass trends, but metabolites like isethionic acid, homarine, and trigonelline linked to eukaryotic phytoplankton were enriched at the deep chlorophyll maximum of the cyclonic features, while degradation products such as arsenobetaine were enriched in anticyclones. In every analysis, the metabolites with the strongest responses were identified using LC-MS through untargeted metabolomics approaches, highlighting that the molecules most sensitive to environmental perturbation were not among the previously characterized metabolome. By analyzing depth variability (accounting for 20-40% of metabolomic variability across ~150 meters) and the vertical displacement of isopycnal surfaces (explaining 10-20% of variability across a sea level anomaly range of 40 centimeters and a spatial distance of 300 kilometers), this analysis constrains the importance of mesoscale eddies in shaping the chemical composition of particulate matter in the largest biomes on the planet.

1 Introduction

High frequency observations at Station ALOHA in the North Pacific Subtropical Gyre (NPSG) over the past 25 years have revealed temporal and spatial variability in what had previously been considered a relatively homogenous environment (Karl, 1999; Karl and Church, 2017; Karl et al., 2021). A major source of variability comes in the form of mesoscale eddies in which water is entrained into circular surface currents tens to hundreds of kilometers in diameter (McGillicuddy, 2016; Karl and Church, 2017).

Eddies can be observed via satellite altimetry, which measures anomalies in sea surface height. Those that have a positive sea level anomaly (SLA) typically indicate regions where deep water layers and isopycnal surfaces are depressed and are associated with clockwise rotation in the Northern Hemisphere. In contrast, a negative SLA corresponds to counterclockwise rotation in the Northern Hemisphere and the uplift of deep water layers into the sunlit region of the ocean (McGillicuddy, 2016). Mode-water eddies are an exception to the conventions established above (Sweeney et al., 2003; McGillicuddy et al., 2007) but here we only focus on the cyclonic (negative SLA) and anticyclonic (positive SLA) mesoscale features that are commonly observed near station ALOHA (Barone et al., 2019).

The uplift of deep, nutrient-rich seawater into the euphotic zone alters microbial communities (Rii et al., 2008, Rii et al., 2022). Measurements of chlorophyll in cyclones reveal a shallower and more intense maximum (Cornec et al., 2021; Barone et al., 2022) as a result of eukaryotic phytoplankton thriving in the higher nutrient concentrations while cyanobacterial biomass is reduced (Hawco et al., 2021). This growth of large eukaryotes corresponds to a net increase in biomass and primary productivity (Benitez-Nelson et al., 2007), especially near the deep chlorophyll maximum (Barone et al., 2022). In contrast, anticyclones produce conditions favorable to nitrogen fixation by accumulating diazotrophs such as Crocosphaera and Trichodesmium (Fong et al., 2008; Olson et al., 2015; Dugenne et al., 2020, Dugenne et al., 2023).

This shift in community structure and productivity should result in a corresponding shift in the composition of the particulate organic matter produced. Eukaryotic organisms have distinct chemical fingerprints from those of cyanobacteria and heterotrophic bacteria (Heal et al., 2021; Durham et al., 2022; Kuhlisch et al., 2023) and the relative contribution of these different taxa to total biomass should reflect their chemical composition. However, recent work from Harke et al. (2021) showed that overall community function may be robust across eddy types at the surface while Gleich et al. (2024) detected significant differences in protistan community composition and metabolic potential down to 250 meters, with heterotrophy-associated protistan transcripts enriched in the cyclone. The chemical composition of organic matter is one metric of community function since community metabolism and interactions are in part mediated through organic molecules, making their measurement useful for determining shifts in community dynamics.

In this work, we directly tested for differences in the chemical composition of particulate matter due to changes in eddy state in the NPSG. We sampled from cyclonic and anticyclonic eddies spatially near each other and used liquid chromatography-mass spectrometry (LC-MS) for both targeted and untargeted metabolomics to explore changes in the composition of small, polar molecules. We expected to find that overall metabolite abundance would reflect shifts in biomass over multiple depths. We also expected to see enrichment in the deep chlorophyll maximum of cyclonic features for those metabolites especially abundant in eukaryotic organisms. Finally, we predicted that these patterns would be robust across years and sampling regimes for a reliable way to link these ocean features with the chemical composition of organic matter in the open ocean.

2 Materials and Methods

2.1 Cruise information

Samples were collected from two cruises in the North Pacific Subtropical Gyre near Station ALOHA that targeted strong mesoscale eddy features as described in Dugenne et al. (2023) and Gleich et al. (2024) (Figure 1). Briefly, the 2017 MESO-SCOPE cruise (Microbial Ecology of the Surface Ocean-Simons Collaboration on Ocean Processes and Ecology, KM1709 on the R/V Kilo Moana) consisted of a transect across adjacent cyclonic and anticyclonic eddies as well as two long-term Lagrangian stations at the center of each eddy. The cyclonic eddy had a maximum negative sea level anomaly (SLA) of -20 centimeters and the anticyclonic eddy reached +24 centimeters. The transect samples were taken at various times of day as the ship transected the adjacent eddies while the eddy center samples were all collected between 5 and 8 pm. In 2018, the Hawaiian Eddy Experiment (HEE, FK180310 on the R/V Falkor) targeted new cyclonic and anticyclonic eddies in approximately the same location with samples taken at the center of each eddy (maximum negative anomaly in the cyclone = -15 cm, maximum positive anomaly in the anticyclone = +26 cm) (Figure 1).

Figure 1

Figure 1. Sampling during MESO-SCOPE and the Hawaiian Eddy Experiment (HEE). Cruise bounds are shown in the large central map with Station ALOHA colored as a point in red near the Hawaiian Islands. Yellow stars denote sampling locations and station numbers relative to the sea level anomaly contours in the background during 28 June 2017 for the MESO-SCOPE cruise (left) and the 6 April 2018 for the HEE cruise (right). Lagrangian sampling near the eddy centers of the MESO-SCOPE cruise was performed by following drifters deployed in proximity of Stations 6 and 12.

2.2 Biogeochemical data

Biogeochemical measurements were collected as described in Barone et al. (2022) and Dugenne et al. (2023) mostly following protocols used by the HOT program (http://hahana.soest.hawaii.edu/hot/methods/results.html). Briefly, all environmental data was collected via CTD rosette except SLA which was measured via satellite. Particulate carbon and particulate nitrogen were measured using an elemental analyzer. Nitrate + nitrite (N+N) was measured on an autoanalyzer except where concentrations were below 100 nM in which case they were measured using chemiluminescence. Soluble reactive phosphorus (phosphate, ${\text{PO}}_4^{3-}$) was measured on an autoanalyzer or following the magnesium induced coprecipitation method. Values were interpolated linearly to the depths at which metabolite samples were taken using all data collected at the given station. Two particulate carbon values were determined to be spurious, with concentrations 2-3 times higher than typical Station ALOHA values, and were instead estimated using a best-fit linear regression against beam attenuation. CTD data is available here in Supplementary Tables 1 and Supplementary Table 2.

2.3 Sample collection for particulate metabolites

Samples were obtained using the onboard CTD rosette to collect water from 15 meters, the deep chlorophyll maximum (DCM), and 175 meters during the MESO-SCOPE eddy transect and from the DCM ±10 and ±20 meters during the eddy center sampling in the MESO-SCOPE cruise and from 25 meters and the DCM during the HEE cruise. The DCM was determined visually from fluorometer data during the CTD downcast and Niskin bottles were tripped during the return trip to the surface. Seawater from each depth was sampled in triplicate by firing one Niskin bottle for each sample. Samples were brought to the surface and decanted into prewashed (3x with DI, 3x with sampled seawater) polycarbonate bottles for filtration. 10L samples were filtered by peristaltic pump onto 142mm 0.2 µm Durapore filters held by polycarbonate filter holders on a Masterflex tubing line. Pressures were kept as low as possible while still producing a reasonable rate of flow through the filter, approximately 250-500 mL per minute. Samples were then removed from the filter holder using solvent-washed tweezers and placed into pre-combusted aluminum foil packets that were then flash-frozen in liquid nitrogen before being stored at -80°C until extraction. A methodological blank was also collected by running filtrate through a new filter and then treated identically to the samples.

2.4 Metabolite sample extraction

Extraction followed a modified Bligh & Dyer approach as detailed in Boysen et al. (2018). Briefly, filters were randomized and added to PTFE centrifuge tubes with a 1:1 mix of 100 µm and 400 µm silica beads, approximately 2 mL -20°C Optima-grade dichloromethane, and approximately 3 mL -20°C 1:1 methanol/water solution (both also Optima-grade). Isotope labeled extraction standards were also added during this step (Supplementary Table 4). The samples were then bead-beaten three times, followed by triplicate extraction of the aqueous layer into a separate vial after spontaneous separation and replacement with fresh methanol/water mixture using a glass Pasteur pipette with additional bead-beatings in between. The aqueous fraction was then dried down under ultrapure nitrogen gas and warmed using a Fisher-Scientific Reacti-Therm module. Dry samples were reconstituted in 380 µL Optima-grade water and amended with 20 µL isotope-labeled injection standards (Supplementary Table 4) to measure the variability introduced by chromatography and ionization. The reconstituted fraction was syringe-filtered into precombusted glass inserts in liquid chromatography (LC) vials to remove any potential clogging material. Samples were then additionally diluted 1:1 with Optima-grade water to prevent overloading on the column and to reduce salt effects.

A pooled sample was created by combining 20 µL of each sample into the same LC vial, and a 1:1 dilution with water created a half-strength pooled sample from that to assess matrix effects and obscuring variation (Boysen et al., 2018). Also run alongside the environmental samples were samples of authentic standards split into two mixes (Supplementary Table 5). Standards were dissolved in Optima-grade water to help with identification and into an aliquot of the pooled sample to quantify matrix response factors. LC vials containing the samples were frozen at -80°C until thawing shortly before injection.

2.5 LC-MS methods

Separate LC-MS runs were used for the MESO-SCOPE eddy transect, the MESO-SCOPE eddy center, and the HEE cruise data. Eddy center samples were run in February 2018, eddy transect samples were run later during August of that year, and the HEE samples were run in July 2019. Each run was treated as a single batch and injected in sequence while maintaining the randomization that had occurred during extraction to minimize chromatographic shifts from solvent or column switching.

For each run, a Waters Acquity I-Class UPLC with a SeQuant ZIC-pHILIC column (5 µm particle size, 2.1 mm x 150 mm, from Millipore) was used with 10 mM ammonium carbonate in 85:15 acetonitrile to water (Solvent A) and 10 mM ammonium carbonate in 85:15 water to acetonitrile (Solvent B) at a flow rate of 0.15 mL/min. The column was held at 100% A for 2 minutes, ramped to 64% B over 18 minutes, ramped to 100% B over 1 minute, held at 100% B for 5 minutes, and equilibrated at 100% A for 25 minutes (50 minutes total). The column was maintained at 30°C. The injection volume was 2 µL for samples and standard mixes. When starting a batch, the column was equilibrated at the starting conditions for at least 30 minutes. To improve the performance of the HILIC column, we maintained the same injection volume, kept the instrument running water blanks between samples as necessary, and injected standards in a representative matrix (the pooled sample) in addition to standards in water. After each batch, the column was flushed with 10 mM ammonium carbonate in 85:15 water to acetonitrile for 20 to 30 minutes.

The Waters Acquity UPLC was coupled to a Thermo Q Exactive HF hybrid Orbitrap high resolution mass spectrometer equipped with a heated electrospray ionization source (H-ESI). The H-ESI voltage was set to 3.3 kV and sheath gas, auxiliary gas, and sweep gas flow rates were set at 16, 3, and 1, respectively. The capillary and auxiliary gas heater temperatures were maintained at 320°C and 100°C, respectively. Full scan analyses were performed with polarity switching and a scan range of 60 to 900 m/z at a resolution of 60,000. The instrument was mass calibrated at 200 m/z before each run began and once again during the eddy transect sample set to ensure calibrations were always performed within 3-4 days. All data files were then converted to an open-source mzML format and vendor centroided via Proteowizard’s msConvert tool (version 3.0.19297, Chambers et al. (2012)). All mzML files have been uploaded to Metabolomics Workbench under Project ID PR001738.

2.6 Metabolomic data analysis

2.6.1 Full scan feature extraction for known molecules

Compounds for which we had an authentic standard (Supplementary Table 5) were manually integrated using Skyline (versions 4.1 and 21.2.0.470, Adams et al. (2020)) and MSDIAL (version 4.36, Tsugawa et al. (2015)) with the standard mixes used to ensure the correct peak was integrated by matching retention time. Compounds were removed from the analysis if they were visually assessed to have poor peak quality in the samples or peak areas in the methodological blanks similar to those in the samples. Skyline was used for both the eddy center samples and the HEE cruise data while MSDIAL was used for the eddy transect data due to the larger number of files in that dataset. Raw peak areas were then normalized to their best-matched internal standard (Supplementary Table 4) per Boysen et al. (2018) except that all compounds were normalized to their best match no matter how minimal the improvement. While significant analytical drift was detected during the longer eddy transect sample run, the use of internal standards was able to negate the instrument’s general increase in sensitivity. Normalized peaks were then calibrated via a single-point standard curve as determined uniquely for each compound from the authentic standard mixes and used to convert normalized peak area into moles per µL injected and scaled to the amount of seawater filtered to provide a final estimate of environmental concentration for each compound in moles per liter of seawater filtered. Scripts used for this analysis are available at https://github.com/wkumler/MesoscopeMetabolomicsManuscript/tree/master/targeted.

2.6.2 Automated processing for detection of unknown features

Samples were also processed via the untargeted workflow detailed in Kumler et al. (2023). In short, XCMS (version 3.10.0, Smith et al. (2023)) was used for feature detection with CentWave peakpicking, Obiwarp alignment, and peak density correspondence/grouping. This was followed by manual inspection of the extracted ion chromatogram for each mass feature by an expert for quality with only those peaks classified as “Good” used in the analysis here. Features were then matched with and normalized to an internal standard (Supplementary Table 4) as described above. Scripts used for this analysis are available at https://github.com/wkumler/MesoscopeMetabolomicsManuscript/tree/master/untargeted.

2.7 Statistics

All statistics were run in R version 4.4.0. We used multivariate statistics provided by the vegan package (version 2.6-6.1, Oksanen et al. (2022)) to assess the overall impact of various metrics across the metabolome followed by univariate statistics to investigate individually significant responses. Non-parametric tests were used when data were unlikely to obey parametric assumptions and permutational statistics were preferred wherever possible. Multiple comparison problems were controlled where necessary using a false-discovery rate correction (Benjamini and Hochberg, 1995). For the analyses involving the eddy transect samples, sea level anomaly was treated as a continuous variable with the exception of the PERMANOVA, for which SLA was categorized as cyclonic if less than -5 cm and anticyclonic if greater than 5 cm, with a third category for stations in between (Anderson, 2001). Both the MESO-SCOPE eddy center and HEE analyses treated SLA as a categorical variable. Where time of day was included as an explanatory factor in the PERMANOVA, it was binned into 6-hour increments starting from 3AM (e.g. Midday contains casts collected between 9AM and 3PM). Marginal effects are reported for the PERMANOVA analyses to avoid influencing the results with the order of the terms in the model formula. Total sample size (n) for the eddy transect was 33 per depth (3 biological triplicates at each of 11 stations), while the eddy centers had an n of 15 per eddy (3 biological replicates at 5 depths) and the HEE data had an n of 24 (3 biological replicates at 4 stations (2 in the cyclone and 2 in the anticyclone) at 2 depths each).

3 Results

We explored the impacts of isopycnal uplift and depression on the composition of particulate organic matter across pairs of mesoscale eddies of opposite polarity in the North Pacific Subtropical Gyre (NPSG). The first eddy pair was the focus of the July 2017 cruise, Microbial Ecology of the Surface Ocean-Simons Collaboration on Ocean Processes and Ecology (MESO-SCOPE, KM1709 aboard the R/V Kilo Moana). During the first phase of our two-phase expedition, a transect was taken from the north to the south consisting of 11 stations across a pair of eddies (the same sampling sites as in Barone et al. (2022) and Dugenne et al. (2023), Figure 1). Samples for particulate metabolomic analysis were collected from 15 meters, the deep chlorophyll maximum (DCM), and 175 meters at each station onto 0.2 µm filters. The next phase of the MESO-SCOPE cruise focused on Lagrangian sampling in the two eddy centers, starting with the cyclone (SLA = -14 cm) before progressing to the anticyclone (SLA = 24 cm). There, metabolomics samples were taken from the DCM as well as 10 and 20 meters above and below it. The second eddy pair was targeted during the March 2018 Hawaiian Eddy Experiment (HEE, FK180310 aboard the R/V Falkor). During this followup cruise, a strong anticyclone (SLA = 21 cm) and a nearby cyclone (SLA = -13 cm) were targeted for a similar set of biogeochemical measurements described in Dugenne et al. (2023) and Gleich et al. (2024) (Figure 1). Metabolomics samples were taken from the center of each eddy as described above at both 25 meters and the DCM. We analyzed these three different datasets separately and discuss the findings from each in sequence below.

3.1 Metabolome variability across adjacent mesoscale eddies of opposing polarity during MESO-SCOPE explained by sea level anomaly variations

The metabolome clearly differed between cyclonic and anticyclonic samples, though the magnitude of this difference varied by depth. The largest differences in particulate matter composition were detected in the DCM and 175 meter samples, as shown by the NMDS in Figure 2.

Figure 2

Figure 2. Distribution of metabolites in multivariate space across adjacent eddies of opposite polarity during the MESO-SCOPE transect, broken down by depth. Top panels depict non-metric multidimensional scaling (NMDS) plots with individual samples colored based on their corrected sea level anomaly. NMDS stress values (s) have been reported in the bottom left corner, while PERMANOVA R² and p-values are reported in the top left. SLA trends are visible in the DCM and 175 meter samples, with dark blue circles consistently discriminating from the dark red circles along the first multidimensional axis. Bottom panels depict the direction and magnitude of this effect by plotting the mean value of all z-scored metabolites across three biological triplicates in color with the raw values in black behind.

At the DCM, approximately 14% of the variation in the metabolome could be explained by sea level anomaly (SLA) alone (PERMANOVA, p = 0.002). Additionally, SLA was strongly correlated with the first principal component (PC) of the metabolome with a Pearson’s r of 0.615 and 27% of the variance explained by PC1 (Supplementary Figure 1A), indicating that this was one of the largest sources of variation in the dataset. Notably, the samples taken from the exact center of the cyclonic eddy at Station 12 were highly distinct and likely drove much of the explained variance. In the 175 meter samples, 13% of the variance was explained by SLA (PERMANOVA, p = 0.004) and SLA was highly correlated with the second PC of the dataset (r = 0.697, % variance explained by PC2 = 15.4, Supplementary Figure 1A), though the first PC did not seem to have any visible pattern with metadata and appeared to largely capture variation between biological triplicates (Supplementary Figure 1B). At 15 meters, SLA trends are slightly less evident with a larger p-value (PERMANOVA, p=0.029) and a lower ${\text{R}}^2$ (0.10) (Figure 2). Time of day was also a significant factor in the PERMANOVA at the surface and DCM (p-values of 0.018 and 0.033, respectively) while the diel effects in the 175 meter samples were unsurprisingly much lower (p = 0.140). The significance of SLA as an explanatory factor was similar even when excluding time of day from the model expression (${\text{R}}_{SLA}^2$ and p-values within permutation error).

To characterize the difference between eddies in a lower-dimensional space, we plotted an average metabolite peak area per station at each depth. We performed z-score normalization on each mass feature to give them all equal weight while preserving the variance between samples then calculated the median z-score for each sample, resulting in large positive values when metabolites are more abundant in a given sample than in the median sample (Figure 2). The average metabolite collected at the exact center of the cyclonic eddy (Station 12) had higher values than other samples at this same depth, indicating that most metabolites were more abundant at the cyclone’s DCM than elsewhere in the transect. The 175 meter samples, on the other hand, showed the largest metabolite abundances in the anticyclone while the cyclone had consistently lower median z-scores. This method highlighted the negative correlation between the average metabolite and SLA at the DCM ($r_{DCM}=-0.716$) while the 15 meter and 175 meter samples showed the opposite trend ($r_{15m}=0.211$, $r_{175m}=0.830$).

This shift in metabolite abundance was largely driven by a corresponding shift in biomass. Particulate carbon and summed metabolite concentration were tightly correlated across all samples and depths (Type I linear regression $\beta =24.9\pm 1.82$ (SE), Pearson’s $r$ = 0.813, p-value< 0.001, Supplementary Figure 2). This overall trend was largely driven by higher values at 15 meters and the DCM than at depth, but when analyzed at each depth individually the general correspondence held (Supplementary Figure 2). The DCM and 175 meter samples had the stronger trends (${\beta }_{DCM}=22.6\pm 5.53$, $r=0.593$, p-value< 0.001; ${\beta }_{175m}=23.9\pm 7.36$, $r=0.503$, p-value = 0.003) while 15 meter samples showed no significant relationship (${\beta }_{15m}=3.4\pm 12.9$, $r=0.049$, p-value = 0.792). This meant that metabolites contributed a relatively fixed fraction of carbon to the particulate pool with our 53 quantified metabolites representing between 1.5% and 5% of the carbon in the system (Figure 3). Particulate nitrogen was tightly correlated with particulate carbon ($r=0.933$) but the fraction of particulate nitrogen represented by the quantified metabolites was slightly higher across the transect with contributions typically between 2% and 6%.

Figure 3

Figure 3. Differences in known metabolite concentration across the pair of adjacent eddies in the MESO-SCOPE transect separated by depth. Top panels depict concentrations of known compounds measured, where bar height corresponds to median triplicate concentration for each metabolite with the top 9 shown and the 44 other identified metabolites summed in grey. TMAO, trimethylamine N-oxide; HO-Ile, hydroxyisoleucine; GBT, glycine betaine; DCM, deep chlorophyll maximum. Center panels show the corresponding measurements of particulate carbon (PC), while the lower panels depict the fraction of total particulate carbon in the known metabolites. The mean value of three biological triplicates is shown in color and the raw values in black behind. Colors correspond to corrected sea level anomaly (Corr. SLA), with dark red indicating anticyclonic (positive SLA) and dark blue indicating cyclonic (negative SLA) eddy state.

Although the untargeted pool contains more molecular features, we expect that a majority of the signal has been captured by our 53 targeted compounds. Of the total peak area that passed quality control, approximately 67% was captured by these 53 molecules. An additional 15% of the total peak area was putatively annotated as various inorganic ions for which no authentic standard was available but could be matched by mass and isotope pattern. Additionally, our targeted list accounted for 13 out of the top 20 largest molecular features by total peak area with an additional 4 features annotated as the inorganic ions for 17/20 of the largest features known. However, these calculations should be taken as estimates because peak area does not correspond directly to environmental concentration and it is entirely possible for abundant environmental compounds to have small peak areas if they ionize poorly on the mass spectrometer or are removed during the sampling and extraction process.

Samples taken from the DCM of the cyclone center (Station 12, Figure 2) showed especially high particulate carbon (67% increase, rising from 1.71 μM across all other DCM samples to 2.85 μM at Station 12) and the highest measured total metabolite concentrations (73% higher, increasing from 10.0 nM to 17.2 nM). This explains the large separation observed between these samples and the rest of the data in the NMDS plot of Figure 2.

The most abundant intracellular metabolite quantified in the eddy transect samples was trimethylamine N-oxide, with an average concentration of 1.38 nM and a distinct DCM maximum (Figure 3), though a few samples had concentrations an order of magnitude higher for unknown reasons (Supplementary Figure 3). A majority of the known molecules (37/53) had a similar pattern with a subsurface maximum at the DCM, including glycine betaine, glutamate, and guanine. Thirteen molecules obeyed a different pattern and decreased monotonically with depth, including the molecules hydroxyisoleucine, gonyol, and dimethylsulfoniopropionate (DMSP). Two known molecules (arsenobetaine and O-acetylcarnitine) increased with depth, potentially representing particle degradation during export.

We also measured several phosphorus-containing compounds and their non-phosphorus equivalents (Supplementary Figure 4A). Phosphocholine is the polar headgroup for phosphatidylcholine, a class of intact polar diacylglycerols that has been previously described at this study site as showing phosphorus stress relief among eukaryotic organisms in response to the upwelling of phosphate-rich deep sea water (Bent et al., 2024). We measured this metabolite as well as choline (the non-phosphate containing equivalent) and glycerophosphocholine (GPC, which contains an additional glycerol group) with the expectation that phosphorus stress would result in higher choline:phosphocholine and choline:GPC ratios, as shown in Seelen et al. (in review). However, we saw no differences in this ratio across the eddy transect after accounting for depth effects (Supplementary Figure 4C). Additionally, choline:GPC ratios were highest at 175m where ${\text{PO}}_4^{3-}$ concentrations are highest, indicating that this metric may not accurately measure P-stress. Other phosphorus-containing lipid headgroups glycerophosphoglycerol and glycerophosphoethanolamine also showed no trend across the eddy at a given depth and instead showed trends functionally identical to the average metabolite in Figure 2 (Supplementary Figure 4B).

3.2 Particulate organic matter compositional shifts across eddy transect

In addition to absolute shifts in metabolite concentration due to biomass differences, we also explored compositional shifts in the metabolomes by assessing the fraction of the total metabolite pool contributed by each compound across the eddy transect. This allowed us to isolate the signal due to shifts in community composition and organismal response from those introduced by changes in biomass of the community as a whole.

In the 15 meter surface samples, only seven compounds had significant changes in relative peak area across the eddy transect, two of which were known: 5-oxoproline, which decreased from 0.2-0.3% of the total peak area at the center of the anticyclone to 0.05-0.1% in the cyclone; and a combined peak of sarcosine and beta-alanine which increased from 0.01% to 0.02-0.03%. The largest shift significant at an α = 0.05 level was an unknown mass feature with m/z = 189.12338 and a retention time around 8.6 minutes, the relative contribution of which increased from approximately 0.2-0.3% in the cyclone to ~1.5% in the anticyclone.

At the DCM, twenty-two compounds changed significantly in relative peak area across the transect. Seven of these were known metabolites, of which trigonelline and homarine were most enriched in the cyclone while arsenobetaine was the only metabolite with a larger fraction of the total peak area in the anticyclone. Homarine showed a very large and highly significant shift, representing about 2.5% of the total peak area in the anticyclone but 7-8% in the cyclone center. Trigonelline (N-methyl niacin), an isomer structurally very similar to homarine but biologically distinct, had a surprisingly strong correlation with homarine (Pearson’s $r=0.855$). Its peak areas were consistently around one third of homarine’s but still large enough to have the second-largest shift in relative peak area of the significantly different compounds. Similarly, arsenobetaine represented about 0.1% of the peak area in the cyclone and 0.5-0.6% of the total peak area in the anticyclone. The most significantly different compounds at the DCM in each direction, however, were both unknowns. The lowest p-value ($6.4\times {10}^{-6}$ after FDR correction) in the DCM data was a mass feature with an m/z of 173.09211, a retention time also around 8.6 minutes, and enrichment in the anticyclone with a putative chemical formula of [M+H] = C₇H₁₃N₂O₃ possibly glycylproline or prolylglycine. The next-lowest (p = $2.5\times {10}^{-4}$ after FDR correction) was enriched in the cyclone and had an m/z of 275.0712 and a retention time of 11.4 minutes with a putative formula of C₁₈H₁₁O₃.

Finally, in the 175 meter samples we detected 43 mass features with a significant association with SLA. Notably, all detected nucleobases (guanine, adenine, and cytosine) were positively associated with SLA (higher values in the anticyclone than the cyclone), though the strongest associations were found in O-acetylcarnitine, betonicine, and tyrosine. Both O-acetylcarnitine and betonicine effectively quintupled their contribution to total peak area, shifting from approximately 0.1% to ~0.5% over the transect. Acetylcholine was the only known compound more relatively abundant in the cyclone than in the anticyclone along with one unknown of m/z 131.0340. The unknown at m/z 173.0921 with the strongest trend at the DCM again had the strongest trend in the 175 meter samples, here with an FDR-corrected p-value of $2.0\times {10}^{-7}$.

3.3 High vertical resolution sampling near eddy center DCM reveals distinct metabolome between cyclone and anticyclone

We further investigated the metabolomic response to SLA at the DCM by collecting samples taken at high-resolution depth intervals around the DCM at the center of both eddy poles. Here, we found that this ~40 meter depth range and the transition between the eddies explained similar amounts of variation in the data (${\text{R}}_{depth}^2=0.253$, ${\text{R}}_{SLA}^2=0.238$) and were both highly significant factors (permutational p-values << 0.001, Figure 4). Additionally, both sample depth and SLA were correlated with the first principal component of the metabolite matrix ($r_{depth}=0.766$, $r_{SLA}=0.684$, fraction of variance explained by PC1 = 30.8%). The cyclonic samples in particular show a much larger spread than the anticyclonic ones, indicating larger sample variability in the cyclone DCM relative to the anticyclone. A depth gradient is visible along NMDS 2, with the deepest samples generally at the top right of the plot and the shallowest ones closer to the bottom (Figure 4).

Figure 4

Figure 4. NMDS plot of high-resolution depth sampling around the deep chlorophyll maximum (DCM, ~115 meters) at the two eddy centers during MESO-SCOPE. Red upward-pointed triangles are from the anticyclone and blue downward-pointed ones are from the cyclone. Shading intensity reflects the depth above or below the DCM. PERMANOVA estimates of the variance explained by depth and sea level anomaly (SLA) are noted in the upper left corner and the NMDS stress value is reported in the bottom left.

To characterize the observed differences in multivariate space we used k-means clustering as an unsupervised way to identify dominant trends within the data (Figure 5). We found that a majority of the metabolites (Clusters 1 and 2, 63% of the total) fell into clusters with larger peak areas in the cyclone while also detecting 32 metabolites that clustered such that the mean metabolite was enriched in the anticyclone (Cluster 4, Figure 5). Cluster 2 had a distinct decrease in relative peak area with depth, while the other clusters had much less clear depth trends. Cluster 1 also showed a DCM maximum for the samples from the anticyclone, correlating well with flow cytometry counts of picoeukaryotes (Pearson’s $r=0.851$) that was not present in the cyclone ($r=0.013$).

Figure 5

Figure 5. Distribution of metabolites in the high-resolution depth samples from the centers of each MESO-SCOPE eddy. The upper row of plots shows k-means clusters where points denote the average z-scored peak area for both known and unknown metabolites across the samples and are colored by the eddy from which they were taken. Clusters have been ordered by number of metabolites in each group and the total is denoted in the panel titles. Both depth trends (mostly a net decrease in metabolites with depth) and eddy effects (cyclonic enrichment in clusters 1 and 2, anticyclone enrichment in cluster 4) are observable. The lower plot shows the individual known and unknown metabolites where points correspond to the FDR-corrected p-value estimated by the nonparametric Mann-Whitney U test and the log₂ fold-change calculated with the average peak area in the cyclone divided by the average peak area in the anticyclone. Colors have been assigned using the k-means clusters and shapes have been assigned based on the status of the mass feature as either a known metabolite that was matched to an authentic standard or an unidentified metabolite. The dashed line across the figure represents the 0.05 level of significance as a visual cue for metabolites above which the differences between the eddies are unlikely to be due to chance. DMSP, dimethylsulfoniopropionate; DMS-Ac, dimethylsulfonioacetate.

The majority of the individual metabolites in these high-resolution depth profiles differed significantly between the two eddies (121/228, α = 0.05). Of those, 15 were enriched in the anticyclone and 106 were more abundant in the cyclone. As expected, all compounds enriched in the anticyclone were part of Cluster 4 and all those enriched in the cyclone belonged to either Cluster 1 or 2 (Figure 5).

Many of the known metabolites enriched in the cyclone function as osmolytes in the cell and might be enriched due to the increased eukaryotic phytoplankton biomass. However, some metabolites such as isethionate, the reduced sulfur osmolytes dimethylsulfoniopropionate (DMSP) and dimethylsulfonioacetate (DMS-Ac), and the isomers homarine and trigonelline, were enriched in excess of biomass (Figure 5). A few metabolites more abundant in the anticyclone were also given putative identifications based on RT and m/z matching with internal standards run at a later time on the mass spectrometer. Of these, the putative arsenobetaine was the most significantly different among these with a peak area at the anticyclone DCM nearly quadruple that of the cyclone (Figure 5). Of note, the 173.0921 m/z mass feature noted above as strongly enriched in the anticyclone was among the most significantly different between the two eddies, along with two isomers at an m/z of 170.1176 and retention times around 8-9 minutes that increased by a factor of 1.5 in the anticyclone (putative formula C₇H₁₄N₄O.

3.4 Hawaiian Eddy Experiment data reveals a different community and response to mesoscale eddies

Given the strong signals detected in the samples from the MESO-SCOPE cruise both across the eddy transect as well as in the high-resolution DCM sampling, we expected to find similar results in an analogous dataset. Samples were collected during a 2018 cruise on the R/V Falkor that again targeted both a cyclonic and an anticyclonic eddy in the North Pacific Subtropical Gyre near Station ALOHA as part of the Hawaiian Eddy Experiment (HEE).

The HEE data was characterized by high inter-replicate variability relative to the samples from the MESO-SCOPE cruise, with 25 meter samples in particular highly variable in multivariate space (Figure 6). Despite this, we still saw the importance of sea level anomaly as a significant explanatory factor in the dataset (PERMANOVA ${\text{R}}_{SLA}^2=0.090$, p-value = 0.017) as well as the larger depth differences between the surface (25 meters) and DCM (110 - 120 meters) (PERMANOVA ${\text{R}}_{depth}^2=0.251$, p-value<< 0.001). We saw no difference between the 6 AM and 6 PM samples (PERMANOVA ${\text{R}}_{time}^2=0.031$, p-value = 0.367) despite a few known compounds (trehalose and sucrose) showing large differences, likely due to a majority of the mass features demonstrating no diel effect (Supplementary Figure 5).

Figure 6

Figure 6. Non-metric multidimensional scaling (NMDS) plot from the Hawaiian Eddy Experiment cruise data, in which points correspond to individual samples and have been colored and shaped by their source eddy status and shaded by the depth from which they were collected. The NMDS stress value has been reported in the bottom left corner, while PERMANOVA ${\text{R}}^2$ and p-values are reported in the top left. Samples from 25 meters deep are visibly distinct from the deep chlorophyll maximum (DCM) samples and an SLA signal is visible in the 25 meter samples only.

When analyzed separately as distinct depths, the 25 meter samples had higher variances explained by eddy and a lower p-value (PERMANOVA, ${\text{R}}^2$=0.26, p-value=0.013), while the DCM samples were no longer likely to be distinct between the two eddies (PERMANOVA p-value = 0.171) (Figure 6, Supplementary Figure 6).

We were also able to use the HEE dataset to test whether the compounds detected as significantly different in the MESO-SCOPE dataset were also different in this eddy pair. We expected to find abundant TMAO and more hydroxyisoleucine at the surface than the DCM in addition to most compounds enriched in the cyclone DCM with biomass, especially the six compounds mentioned above in the results of Figure 5 (trigonelline, homarine, DMS-Ac, DMSP, taurine, and isethionic acid).

As expected, TMAO was again one of the most abundant metabolites detected in the particulate matter with concentrations around 1.2 nM at 25 meters and 0.4 nM at the DCM, second only to the high levels of the amino acid glycine (2 nM at 25 meters, 0.9 nM DCM), though the differences between 25 meters and DCM were not significant for either compound. We also detected a significant difference between depths for hydroxyisoleucine (25 meter mean = 0.32 nM, DCM mean = 0.21 nM, t-test p-value=0.022 with n=12 samples at each depth). These results imply that the overall depth structure of the metabolome of the gyre is relatively fixed for the most abundant compounds.

The eddy effects, on the other hand, were much less strong during this cruise. Of the six known metabolites with the strong enrichment in the cyclone at the DCM in the MESO-SCOPE cruise, only isethionic acid was significantly different in the HEE dataset (isethionic acid t-test p-value _FDR = 0.043, other five p-values > 0.25). This difference in isethionic acid concentration was only found after normalizing each sample to the sum of all metabolites in the sample to control for biomass. Arsenobetaine was again found to be strongly enriched in the anticyclone at the DCM both when normalizing to biomass and when not, just as in MESO-SCOPE. Surprisingly, the 173.0921 m/z mass feature was also slightly enriched but this time in the cyclone DCM with peak areas approximately 1.4 times the values in the anticyclone (t-test p-value = 0.067, Mann-Whitney p-value = 0.045). The 170.1176 m/z mass feature was not detected in the HEE samples.

4 Discussion

4.1 Multivariate approaches reveal metabolome-wide shifts across pairs of eddies of opposite polarity

We measured the metabolome of samples collected across two sets of adjacent eddies of opposite polarity to explore the effect of sea level anomaly (SLA) on metabolite composition and found that the altered biogeochemistry and microbial community composition between cyclonic and anticyclonic eddies explain a significant portion of the observed variations in their particulate metabolomes.

In all of the datasets analyzed here, we detected a significant difference in the composition of the metabolome between the adjacent eddies. This effect was strongest in the samples taken during the Lagrangian stations at the center of each eddy in the 2017 MESO-SCOPE cruise, with nearly a quarter of the total variance explained by the eddy from which the samples were taken. In the samples taken along the eddy transect, the largest effect was detected in the deepest samples taken from 175 meters and from the DCM, with much less of a response at 25 meters. This lack of SLA effect at the surface is unsurprising given the similarity in the biogeochemistry of the well-lit upper ocean previously reported (Barone et al., 2022). In the HEE samples, however, the opposite was true with a surprisingly larger SLA response at 25 meters than that at the DCM. In each case, the differences introduced by SLA were smaller than those of depth when compared directly, with even the high-resolution sampling around the DCM finding slightly more variance explained by the 40 meter difference in sampling depth than the 40 centimeter difference in SLA. This result agrees well with prior research showing the large effect of sampling depth on the metabolome (Heal et al., 2021; Kumler et al., 2023; Bent et al., 2024). Time of day was included as a significant factor in the MESO-SCOPE transect, with trehalose and sucrose in particular showing large differences between 6 AM and 6 PM as has been previously noted (Muratore et al., 2022), though this effect was diminished at 175 meters.

Although possible that the observed differences in metabolomes across adjacent eddies differing in polarity are due to latitudinal shifts or simply background variation in the gyre environment, the transect data in particular implies that this is unlikely to be the case. Samples taken when the absolute value of the SLA was less than 5 cm were more similar to each other than to those samples taken within the eddies despite larger differences in latitude. The transect sampled both outside of each eddy and between the two across a large spatial gradient (~4° latitude and 2° longitude), yet the non-eddy (absolute value of SLA< 5 cm) samples from all three locations grouped together, showed similar median metabolite concentrations, and were visually similar in the most abundant metabolite composition. These results imply that strong cyclonic and anticyclonic eddies represent endpoints in NPSG composition, in agreement with results from Barone et al. (2019) which found the most extreme values in the Hawaii Ocean Time series data typically detected when eddies passed over Station ALOHA. We also detected the largest differences in the DCM metabolome at the exact center of the eddy, with Station 12 highly distinct from even the nearest stations, while no such stark difference was found in the anticyclone (Figures 1 and 2). This indicates that the exact center of the cyclone has a unique metabolome at the DCM relative to the rest of the transect, rather than existing as a smooth gradient with SLA.

Shifts in total biomass were a large driver of metabolomic differences, with earlier work in the same location by Barone et al. (2022) showing that particulate carbon, chlorophyll, and beam attenuation were all 20-80% higher in the DCM of the cyclone relative to the anticyclone. We see this reflected particularly well in the targeted metabolites measured here, with clear depth trends and a strong correlation between total metabolite concentration and particulate carbon values. However, even after controlling for biomass effects by normalizing each sample to the sum of signal measured within it we still found a significant SLA effect, likely due to shifts in the community composition and in particular the increased eukaryote presence in the cyclone DCM.

The NMDS plots of the high-resolution depth samples were also illustrative of sample dissimilarity with clear SLA and depth trends. The anticyclonic samples grouped tightly and showed little variance with depth or triplicate when compared to the cyclonic samples with the exception of DCM triplicate B, which was highly distinct for unknown reasons. In contrast, the cyclonic samples were much more variable as is perhaps expected when the biomass is in the form of larger phytoplankton that are less homogenous in the environment. One notable aspect of their clustering was the way in which the deepest samples (DCM plus 20m) grouped most closely to the anticyclonic samples, perhaps indicating that below the DCM the metabolome rapidly approaches a uniform deep-water signal. This result was noted previously in the 175 meter samples of the eddy transect in Kumler et al. (2023), where samples from the deep euphotic zone showed greater intra-depth similarity than samples from the DCM or 25 meters.

Given the confluence of both depth and SLA signals and the way that the median metabolite plots of Figure 2 would confound the SLA signal with the depth variance, we used k-means clustering to group similar compounds and provide a reduced dimensionality space for visualization. This revealed four major patterns of metabolite response, with a majority of compounds responding to eddy state (Clusters 1, 2, and 4 in Figure 5). It is interesting to note that the sole cluster in which abundances in the anticyclone are greater than those of the cyclone is also the only cluster to generally increase in abundance with depth (Cluster 4), likely because compounds more concentrated in the anticyclone are the same kind of degradation products and recalcitrant carbon typically found at depth.

The Hawaiian Eddy Experiment (HEE) data upset many of the expectations we developed during the MESO-SCOPE cruise in the previous year. Most surprising was the large difference between the eddies detected at the surface, while the DCM samples were functionally indistinguishable. This contrasts directly with the results from MESO-SCOPE and the results in Gleich et al. (2024), who found eddy-driven shifts in protistan community composition to be larger at depth than at the surface. However, Dugenne et al. (2023) found differences in nitrogen fixation rate and nitrogen fixer composition varying widely throughout the upper water column. The large surface differences were especially surprising given the large inter-replicate differences between the HEE surface samples in this study, while the DCM samples tended to be much more consistent.

4.2 Univariate approaches highlight individual metabolites responding to lifted and depressed isopycnals

Several molecule- or pathway-specific narratives emerged from the metabolomic data. The untargeted approach used here allowed us to describe and characterize signals from small molecules whose identity was unknown, a particularly promising approach in open ocean gyres where the largest fraction of unknowns is found (Heal et al., 2021). We found that our list of authentic standards covered fewer metabolites in anticyclonic eddies at all depths, with the molecular diversity of the cyclone much better characterized. Many of the trends detected persisted even after normalizing within each sample, indicating that the shifts discussed have implications beyond simple scaling with biomass.

Of those molecules whose identity was known, the clearest response to SLA was that of the enigmatic osmolyte homarine at the DCM. This abundant compound increased approximately threefold in concentration from < 50 pM in the anticyclone to around 150 pM in the center of the cyclone in both the eddy transect and eddy center datasets. The pattern was weaker in the noisier HEE data but the largest concentrations were still detected in the cyclone and the median cyclone measurement was greater than the maximum value for the anticyclone. This molecule has a well-established role as an osmolyte in eukaryotic phytoplankton and cyanobacteria (Gebser and Pohnert, 2013; Dawson et al., 2020; Heal et al., 2021; Durham et al., 2022) as well as a documented decrease in abundance with depth (Heal et al., 2021), indicating that its response to SLA is potentially due to differences in physical and chemical attributes that in turn control biomass, community composition, and recycling rates. Curiously, the distribution of homarine was also tightly correlated with that of trigonelline. Although they have structural similarity, they are not known to have a relationship beyond their shared function as osmolytes.

Isethionate is a known osmolyte thus far found exclusively in eukaryotes (Durham et al., 2022) and strongly associated with a few diatoms in particular (Heal et al., 2021). This compound was enriched in the cyclone of the DCM in both the MESO-SCOPE and HEE cruises along with its precursor taurine, validating earlier results from Barone et al. (2022) that showed strong enrichment of eukaryotic phytoplankton at the DCM and Pseudo-nitzschia in particular. Additionally, the abundance other sulfur-containing compounds (e.g. dimethylsulfonioacetate, DMS-Ac and dimethylsulfoniopropionate, DMSP) in the cyclone may indicate links between the formation of eukaryote-rich cyclonic features and the production of dimethylsulfide, a climate-active volatile gas (Moran and Durham, 2019).

The strong SLA signal detected among the 175 meter samples was partially surprising to us given our expectation about the strongest effect at the DCM where eddy effects lift large concentrations of nutrients above the 1% light level. Instead, we found more compounds overall to be significantly different across the eddy transect in the 175 meter samples than we did at the DCM (α=0.05, 43 compounds at 175 meters vs 22 at the DCM). This was likely driven by the enhancement of organic matter and heterotrophic picoplankton below the DCM as seen in Barone et al. (2019), with the enrichment of nucleobases in the 175 meter samples additionally hinting at increases in bacterial biomass given that nucleobases tend to be highly abundant in bacterial cultures (Heal et al., 2021). The abundance of a mass feature putatively identified as acetylcholine in the 175 meter cyclone samples followed the same general trend as arsenobetaine, both enriched at depth and in the anticyclone, and may be another compound that results from heterotrophic degradation of phytoplankton metabolites (Durham et al., 2022). These results raise important questions about the depth at which the community is isolated from SLA effects, with additional data from Barone et al. (2022) and Gleich et al. (2024) indicating that even at 250 meters eddy effects are still discernable.

4.3 Metabolites link changes in biogeochemistry to community biology

The composition of particulate matter controls the turnover of organic carbon in the ocean, but it is still unclear how this composition is altered in the environment and the plankton community. Measuring metabolites can quantify the variability of the central molecules responsible for the flow of energy and matter through the marine ecosystem at many different scales. In these results, we focused on mesoscale eddies, which perturb the ocean ecosystem on relatively small horizontal scales on the order of tens to hundreds of kilometers.

Here, we reported how changes in the distribution of water layers of different density driven by mesoscale physical dynamics cascading effects in the biological and chemical characteristics of the eddies. The enrichment of eukaryotic organisms at the DCM of the cyclone mirrors other ocean processes where upwelling of deep, nutrient-rich water leads to increases in sulfur-rich metabolites representative of a shift towards organisms such as diatoms and dinoflagellates that undergo bloom dynamics (Durham et al., 2022; Kuhlisch et al., 2023). Barone et al. (2022) proposed that a “chemical wake” is formed after cyclone intensification when a transient phytoplankton community with net photosynthetic production induces export and leaves deficits and excesses of inorganic nutrients and oxygen, respectively. This conceptual model that can be extended to include many of the metabolites named here, as some metabolites are more bioavailable than others and the “flavor” of the chemical wake will change throughout the eddy life cycle.

The prevalence of unknown metabolites relative to known metabolites in the anticyclone may also in part explain why the increased biomass in the cyclone did not result in large increases in export (Barone et al., 2022; Bent et al., 2024). Our list of known compounds covers a large fraction of the small molecules common to most organisms, indicating that unknown molecules are more likely to be rare and therefore more difficult to degrade or a result of prokaryote dominance. If the increased biomass in the cyclone is dominantly in the form of highly bioavailable compounds, it is not surprising that much of it can be easily degraded. This hypothesis is supported by Gleich et al. (2024), who found that heterotrophy-associated transcripts were more abundant in the 2018 cruise cyclone than the anticyclone, indicating that there was fresh material for the heterotrophic bacteria to readily consume.

Almost all large phytoplankton classes in the 2017 eddy pair were enriched in the cyclone DCM relative to the anticyclone (particularly prymnesiophytes, dinoflagellates, diatoms, and green algae), while Prochlorococcus was significantly less abundant and Synechococcus was approximately equal (Barone et al., 2022). Thus, the metabolite data presented here does not help distinguish between large phytoplankton classes but instead provides additional cues about the shift from a metabolome dominated by cyanobacteria producers into one dominated by eukaryotic organisms. This shift occurs naturally over the life cycle of an eddy, with distinctions typically drawn between the eddy intensification, its mature phase, and eventual decay (Rii et al., 2008; McGillicuddy, 2016; Barone et al., 2022; Kuhlisch et al., 2023). Increased diapycnal mixing of nutrients in mature cyclones can alter the DCM community without large effects on integrated primary productivity or export even months after eddy intensification (Barone et al., 2022). We here demonstrated how this dynamic is linked with large changes in the chemical characteristics of the particulate material.

5 Conclusion

This study reports how changes in the biogeochemistry and community composition due to mesoscale eddies affect the metabolome and thereby the composition of particulate organic matter. A transect across adjacent mesoscale eddies of opposite polarity showed that many metabolites track closely with metrics of overall biomass, with eddy effects stronger at 175 meters than at 15 meters or the deep chlorophyll maximum (DCM). High-resolution depth sampling of the DCM at the center of each eddy elucidated several known and unknown biomarkers likely corresponding to the previously documented increase in eukaryotic phytoplankton in the cyclone. Metabolite clusters aligned well with expected trends in abundance with depth and between the two eddies, with convergence below the DCM towards a general deep-sea metabolite signal for many metabolites. By contrasting these effects with a follow-up analysis in the centers of a separate pair of mesoscale eddies of opposite polarity the following year, we learned that the impact of eddies may be influenced by their origin, age, and period of the year. Finally, the metabolites with the strongest responses to eddy effects were unknown, and a relatively smaller proportion of known metabolites in the anticyclone demonstrates the utility of untargeted metabolomics in exploring environmental variation and the large amount of information that is missed by only investigating known compounds. In this work, we have shown the sensitivity of the metabolome to various environmental and community factors and constrained the importance of eddy effects in models of elemental cycling.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: http://dx.doi.org/10.21228/M82719, Metabolomics Workbench Project ID PR001738.

Author contributions

WK: Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. WQ: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing – review & editing. RL: Conceptualization, Data curation, Investigation, Methodology, Writing – review & editing. BB: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing. LC: Data curation, Methodology, Project administration, Writing – review & editing. AI: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by grants from the Simons Foundation (SCOPE Award ID 329108 to AI, SF Award ID 385428 to AI, Postdoctoral Fellowship in Marine Microbial Ecology ID 548565 to WQ, and SCOPE Award ID 721264 to David Karl which supported BB).

Acknowledgments

The authors would like to thank the other members of the Ingalls Lab who provided assistance in sample processing, integration, and analysis. We are also grateful to the SCOPE science team (Tara Clemente and Tim Burrell) for sample collection during the Hawaiian Eddy Experiment. We would like to acknowledge the captain and crew of the R/V Kilo Moana and R/V Falkor during the KM1709 and FK180310 cruises for making this science possible. We would also like to acknowledge the analysts involved in the biogeochemical measurements for their willingness to provide data and answer questions (Eric Grabowski, Karin Bjorkman, and Rhea Foreman).

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2024.1481409/full#supplementary-material

Supplementary Figure 1 | Results of a principal component analysis performed separately on each depth of the eddy transect metabolome. (A) Regressions of SLA against various principal components, broken down by depth. % variance explained by each PC is noted in the upper right corner of the plot. Colors have been added corresponding to the SLA value for clarity but are redundant with the x-axis values. Lines of best fit and the associated standard error have been added behind the points in black and translucent grey, respectively. Strong associations exist between PC1 and SLA in the DCM samples and PC2 and SLA in the 175 meter samples. (B) Principal component 1 values for the 175 meter samples at each station showing large variance between triplicates rather than any external factor.

Supplementary Figure 2 | Type I linear regressions of particulate carbon and nitrogen against total nM carbon and nitrogen in known metabolites both as a whole and divided by depth. Points have been colored according to depth (top row of plots) or SLA (bottom row of plots) and have been placed on independent axes. Two outlier particulate carbon points have been interpolated from beam attenuation for 15 meter values at stations 10 and 11.

Supplementary Figure 3 | Stacked barplot of top metabolites by triplicate. Stacked barplots of known metabolite concentration as shown in main text (Figure 3) but with individual triplicates shown. The same data is shown in A) as B) but B) is rendered in relative contribution space instead of absolute concentration. TMAO = trimethylamine N-oxide, HO-Ile = hydroxyisoleucine, GBT = glycine betaine, DCM = deep chlorophyll maximum.

Supplementary Figure 4 | Evidence against phosphorus limitation across the MESO-SCOPE eddy transect. Top row of plots (A) correspond to extracted ion chromatograms for 5 metabolites useful for estimating the phosphorus stress of an environment, colored by the depth from which the sample was taken. Middle row of plots shows the integrated peak area for each metabolite across the 11 stations of the transect with sub-panels for each depth independently, with each point colored by the corrected sea level anomaly (SLA). Bottom row of plots shows the ratio of choline to phosphocholine (left) and choline to glycerophosphocholine (right) as a function of corrected SLA with no statistically significant trend visible at any of the three depths.

Supplementary Figure 5 | Boxplots showing the concentrations of six selected compounds between cruises and depths as well as time and SLA classes. The abscissa is separated by the time category into which sampling was binned, referring to 6-hour increments starting from 3AM (e.g. Midday contains casts collected between 9AM and 3PM). Colors denote the sea level anomaly bins with a ±5cm cutoff. For the eddy transect, six stations are included in the anticyclone, three in the cyclone, and two were neither. There were two morning stations, 4 midday stations, 2 evening stations, and 3 night stations. Diel effects are clearly visible for sucrose and trehalose while other compounds show very little diel differences. Trehalose was not quantified during the MESO-SCOPE transect for lack of an authentic standard to quantify the response factor.

Supplementary Figure 6 | Plots of the metabolome during the Hawaiian Eddy Experiment, broken down by depth. The top row of plots are non-metric multidimensional scaling (NMDS) plots, in which points correspond to individual samples and have been colored by their corrected sea level height anomaly. SLA trends are visible in the 25 meter samples, with dark blue circles consistently discriminating from the dark red circles. The bottom row of plots show the direction and magnitude of this effect by plotting the grand mean of the normalized metabolite peak areas in the stations taken between the adjacent eddies of opposite polarity.

References

Adams K. J., Pratt B., Bose N., Dubois L. G., St. John-Williams L., Perrott K. M., et al. (2020). Skyline for small molecules: A unifying software package for quantitative metabolomics. J. Proteome Res. 19, 1447–1458. doi: 10.1021/acs.jproteome.9b00640

PubMed Abstract | Crossref Full Text | Google Scholar

Anderson M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46. doi: 10.1111/j.1442-9993.2001.01070.pp.x

Crossref Full Text | Google Scholar

Barone B., Church M. J., Dugenne M., Hawco N. J., Jahn O., White A. E., et al. (2022). Biogeochemical dynamics in adjacent mesoscale eddies of opposite polarity. Global Biogeochemical Cycles 36, e2021GB007115. doi: 10.1029/2021GB007115

Crossref Full Text | Google Scholar

Barone B., Coenen A. R., Beckett S. J., McGillicuddy D. J., Weitz J. S., Karl D. M. (2019). The ecological and biogeochemical state of the north pacific subtropical gyre is linked to sea surface height. J. Mar. Res. 77, 215–245. doi: 10.1357/002224019828474241

Crossref Full Text | Google Scholar

Benitez-Nelson C. R., Bidigare R. R., Dickey T. D., Landry M. R., Leonard C. L., Brown S. L., et al. (2007). Mesoscale eddies drive increased silica export in the subtropical pacific ocean. Science 316, 1017–1021. doi: 10.1126/science.1136221

PubMed Abstract | Crossref Full Text | Google Scholar

Benjamini Y., Hochberg Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Society: Ser. B (Methodological) 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

Crossref Full Text | Google Scholar

Bent S. M., Muratore D., Becker K. W., Barone B., Clemente T., Fredricks H. F., et al. (2024). Lipid biochemical diversity and dynamics reveal phytoplankton nutrient-stress responses and carbon export mechanisms in mesoscale eddies in the North Pacific Subtropical Gyre. Front. Mar. Sci. 11. doi: 10.3389/fmars.2024.1427524

Crossref Full Text | Google Scholar

Boysen A. K., Heal K. R., Carlson L. T., Ingalls A. E. (2018). Best-matched internal standard normalization in liquid chromatography–mass spectrometry metabolomics applied to environmental samples. Analytical Chem. 90, 1363–1369. doi: 10.1021/acs.analchem.7b04400

PubMed Abstract | Crossref Full Text | Google Scholar

Chambers M. C., Maclean B., Burke R., Amodei D., Ruderman D. L., Neumann S., et al. (2012). A cross-platform toolkit for mass spectrometry and proteomics. Nat. Biotechnol. 30, 918–920. doi: 10.1038/nbt.2377

PubMed Abstract | Crossref Full Text | Google Scholar

Cornec M., Laxenaire R., Speich S., Claustre H. (2021). Impact of mesoscale eddies on deep chlorophyll maxima. Geophysical Res. Lett. 48, e2021GL093470. doi: 10.1029/2021GL093470

PubMed Abstract | Crossref Full Text | Google Scholar

Dawson H. M., Heal K. R., Torstensson A., Carlson L. T., Ingalls A. E., Young J. N. (2020). Large diversity in nitrogen- and sulfur-containing compatible solute profiles in polar and temperate diatoms. Integr. Comp. Biol. 60, 1401–1413. doi: 10.1093/icb/icaa133

PubMed Abstract | Crossref Full Text | Google Scholar

Dugenne M., Gradoville M. R., Church M. J., Wilson S. T., Sheyn U., Harke M. J., et al. (2023). Nitrogen fixation in mesoscale eddies of the north pacific subtropical gyre: Patterns and mechanisms. Global Biogeochemical Cycles 37, e2022GB007386. doi: 10.1029/2022GB007386

Crossref Full Text | Google Scholar

Dugenne M., Henderikx Freitas F., Wilson S. T., Karl D. M., White A. E. (2020). Life and death of Crocosphaera sp. In the Pacific Ocean: Fine scale predator–prey dynamics. Limnology Oceanography 65, 2603–2617. doi: 10.1002/lno.11473

Crossref Full Text | Google Scholar

Durham B. P., Boysen A. K., Heal K. R., Carlson L. T., Boccamazzo R., Deodato C. R., et al. (2022). Chemotaxonomic patterns in intracellular metabolites of marine microbial plankton. Front. Mar. Sci. 9. doi: 10.3389/fmars.2022.864796

Crossref Full Text | Google Scholar

Fong A. A., Karl D. M., Lukas R., Letelier R. M., Zehr J. P., Church M. J. (2008). Nitrogen fixation in an anticyclonic eddy in the oligotrophic North Pacific Ocean. ISME J. 2, 663–676. doi: 10.1038/ismej.2008.22

PubMed Abstract | Crossref Full Text | Google Scholar

Gebser B., Pohnert G. (2013). Synchronized regulation of different zwitterionic metabolites in the osmoadaption of phytoplankton. Mar. Drugs 11, 2168–2182. doi: 10.3390/md11062168

PubMed Abstract | Crossref Full Text | Google Scholar

Gleich S. J., Hu S. K., Krinos A. I., Caron D. A. (2024). Protistan community composition and metabolism in the North Pacific Subtropical Gyre: Influences of mesoscale eddies and depth. Environ. Microbiol. 26, e16556. doi: 10.1111/1462-2920.16556

PubMed Abstract | Crossref Full Text | Google Scholar

Harke M. J., Frischkorn K. R., Hennon G. M. M., Haley S. T., Barone B., Karl D. M., et al. (2021). Microbial community transcriptional patterns vary in response to mesoscale forcing in the north pacific subtropical gyre. Environ. Microbiol. 23, 4807–4822. doi: 10.1111/1462-2920.15677

PubMed Abstract | Crossref Full Text | Google Scholar

Hawco N. J., Barone B., Church M. J., Babcock-Adams L., Repeta D. J., Wear E. K., et al. (2021). Iron depletion in the deep chlorophyll maximum: Mesoscale eddies as natural iron fertilization experiments. Global Biogeochemical Cycles 35, e2021GB007112. doi: 10.1029/2021GB007112

Crossref Full Text | Google Scholar

Heal K. R., Durham B. P., Boysen A. K., Carlson L. T., Qin W., Ribalet F., et al. (2021). Marine community metabolomes carry fingerprints of phytoplankton community composition. mSystems 6, e01334–e01320. doi: 10.1128/mSystems.01334-20

PubMed Abstract | Crossref Full Text | Google Scholar

Karl D. M. (1999). Minireviews: A sea of change: Biogeochemical variability in the north pacific subtropical gyre. Ecosystems 2, 181–214. doi: 10.1007/s100219900068

Crossref Full Text | Google Scholar

Karl D. M., Church M. J. (2017). Ecosystem structure and dynamics in the north pacific subtropical gyre: New views of an old ocean. Ecosystems 20, 433–457. doi: 10.1007/s10021-017-0117-0

Crossref Full Text | Google Scholar

Karl D. M., Letelier R. M., Bidigare R. R., Björkman K. M., Church M. J., Dore J. E., et al. (2021). Seasonal-to-decadal scale variability in primary production and particulate matter export at station ALOHA. Prog. Oceanography 195, 102563. doi: 10.1016/j.pocean.2021.102563

Crossref Full Text | Google Scholar

Kuhlisch C., Shemi A., Barak-Gavish N., Schatz D., Vardi A. (2023). Algal blooms in the ocean: Hot spots for chemically mediated microbial interactions. Nat. Rev. Microbiol. 22, 138–154. doi: 10.1038/s41579-023-00975-2

PubMed Abstract | Crossref Full Text | Google Scholar

Kumler W., Hazelton B. J., Ingalls A. E. (2023). Picky with peakpicking: Assessing chromatographic peak quality with simple metrics in metabolomics. BMC Bioinf. 24, 404. doi: 10.1186/s12859-023-05533-4

PubMed Abstract | Crossref Full Text | Google Scholar

McGillicuddy D. J. (2016). Mechanisms of physical-biological-biogeochemical interaction at the oceanic mesoscale. Annu. Rev. Mar. Sci. 8, 125–159. doi: 10.1146/annurev-marine-010814-015606

PubMed Abstract | Crossref Full Text | Google Scholar

McGillicuddy D. J., Anderson L. A., Bates N. R., Bibby T., Buesseler K. O., Carlson C. A., et al. (2007). Eddy/wind interactions stimulate extraordinary mid-ocean plankton blooms. Science 316, 1021–1026. doi: 10.1126/science.1136256

PubMed Abstract | Crossref Full Text | Google Scholar

Moran M. A., Durham B. P. (2019). Sulfur metabolites in the pelagic ocean. Nat. Rev. Microbiol. 17, 665–678. doi: 10.1038/s41579-019-0250-1

PubMed Abstract | Crossref Full Text | Google Scholar

Muratore D., Boysen A. K., Harke M. J., Becker K. W., Casey J. R., Coesel S. N., et al. (2022). Complex marine microbial communities partition metabolism of scarce resources over the diel cycle. Nat. Ecol. Evol. 6, 218–229. doi: 10.1038/s41559-021-01606-w

PubMed Abstract | Crossref Full Text | Google Scholar

Oksanen J., Simpson G. L., Blanchet F. G., Kindt R., Legendre P., Minchin P. R., et al. (2022). Vegan: Community ecology package. Available online at: https://github.com/vegandevs/vegan (Accessed September 20, 2023).

Google Scholar

Olson E. M., McGillicuddy D. J., Flierl G. R., Davis C. S., Dyhrman S. T., Waterbury J. B. (2015). Mesoscale eddies and Trichodesmium spp. Distributions in the southwestern North Atlantic. J. Geophysical Research: Oceans 120, 4129–4150. doi: 10.1002/2015JC010728

PubMed Abstract | Crossref Full Text | Google Scholar

Rii Y. M., Brown S. L., Nencioli F., Kuwahara V., Dickey T., Karl D. M., et al. (2008). The transient oasis: Nutrient-phytoplankton dynamics and particle export in Hawaiian lee cyclones. Deep Sea Res. Part II: Topical Stud. Oceanography 55, 1275–1290. doi: 10.1016/j.dsr2.2008.01.013

Crossref Full Text | Google Scholar

Rii Y. M., Peoples L. M., Karl D. M., Church M. J. (2022). Seasonality and episodic variation in picoeukaryote diversity and structure reveal community resilience to disturbances in the north pacific subtropical gyre. Limnology Oceanography 67, S331–S351. doi: 10.1002/lno.11916

Crossref Full Text | Google Scholar

Smith C. A., Tautenhahn R., Neumann S., Benton P., Conley C., Rainer J. (2023). Xcms: LC-MS and GC-MS data analysis. Available online at: https://github.com/sneumann/xcms (Accessed September 20, 2023).

Google Scholar

Sweeney E. N., McGillicuddy D. J., Buesseler K. O. (2003). Biogeochemical impacts due to mesoscale eddy activity in the sargasso sea as measured at the Bermuda atlantic time-series study (BATS). Deep Sea Res. Part II: Topical Stud. Oceanography 50, 3017–3039. doi: 10.1016/j.dsr2.2003.07.008

Crossref Full Text | Google Scholar

Tsugawa H., Cajka T., Kind T., Ma Y., Higgins B., Ikeda K., et al. (2015). MS-DIAL: Data-independent MS/MS deconvolution for comprehensive metabolome analysis. Nat. Methods 12, 523–526. doi: 10.1038/nmeth.3393

PubMed Abstract | Crossref Full Text | Google Scholar