Introduction

The continued increase in anthropogenic carbon dioxide emissions is leading to ocean acidification, with the ocean absorbing an average of 25% of human-caused emissions annually1. If atmospheric carbon dioxide concentration continues to rise at the current rate, the pH of oceans is predicted to fall 0.3–0.4 units by the end of the century2,3. This pH drop could exacerbate conditions in the U.S. Pacific Northwest, where the ocean pH is lower than that of the global ocean due to natural oceanographic processes including regional upwelling, and could pose a greater challenge to the marine inhabitants already coping with this lower pH. An additional and compounding factor to ocean acidification is ocean deoxygenation, which co-varies with pH and temperature. Given that global ocean temperatures will rise with global warming from continued greenhouse gas emissions, hypoxic zones are expected to increase in duration, intensity, and frequency4. It is not certain how future ocean acidification and deoxygenation environmental stress might affect important Pacific Northwest fisheries like the Dungeness crab fishery, which is the most lucrative and valued at more than $200 million annually5.

Ocean acidification is predicted to have negative indirect effects on Dungeness crab through loss of prey directly affected by ocean acidification6, but knowledge of how ocean acidification might directly impact Dungeness crab is limited. It has been postulated that Dungeness crab may exhibit limited tolerance for acid-base disturbances given past observations of Dungeness crab as weak osmoregulators7,8,9,10 and that acid-base balance and osmoregulation are tightly coupled in decapod crustaceans11. However, it was shown that following a brief two-week exposure to a future-predicted seawater pH of 7.4, adult Dungeness crab are able to acclimate by increasing hemolymph ion levels (bicarbonate, calcium, chloride, sulfate, and sodium), and by decreasing both oxygen consumption and nitrogen excretion12. On the other hand, Dungeness crab larvae have shown reduced survival and development rate in response to a 45-day exposure to the projected future seawater pH of 7.5 and 7.113. Moreover, adult Dungeness crab have shown behavioral responses to declining oxygen conditions (21–1.5 kPa pO2 over a 5-hour period), including reduced feeding and a preference for the area with the highest pO2 level when placed in a seawater oxygen gradient (2.5–10.5 kPa pO2 for 1 hour)14. Dungeness crab also have shown physiological responses to declining oxygen conditions (18–3 kPa pO2 over a 6-hour period), including redistributing hemolymph to high-energy-demand tissues15. Despite these clear biological responses to low pH and oxygen, the biochemical mechanisms underlying Dungeness crab response to combined pH and oxygen stress have not yet been defined.

One way to gain broad insight into the biochemical processes underlying the physiological status of organisms is through surveying metabolomes, which is now possible for nearly any species with the advancement of high-throughput metabolite profiling, also known as metabolomics. Exploratory untargeted metabolomics approaches can offer unbiased analyses of the composition of all detectable metabolites for the rapid and quantitative detection of stress responses, which often leads to the development of targeted approaches and identification of stress-indicating biomarkers16,17. Functional analyses using pathway inference can subsequently be performed by integrating metabolomics data with databases and other “-omics” datasets using bioinformatics tools to ultimately establish causal networks between different experimental conditions and outcomes. Untargeted metabolomics in the context of understanding ocean acidification has recently been applied to reef-building coral, and revealed metabolite profiles that were predictive of primary production activity and molecules that could be used as potential biomarkers of ocean acidification18.

To better understand which biochemical pathways might be altered in the response of Dungeness crab to ocean acidification and given that early life stages of marine organisms will likely be more sensitive to ocean acidification than adults19, we examined the metabolomes and lipidomes of individual Puget Sound juveniles exposed to current pH (7.85) and future pH (7.45) conditions for an average of 32 days starting as megalopae. To account for the dissolved oxygen (DO) levels that naturally co-vary with pH, we included both ambient oxygen (8.9 mg/L or saturated O2) and low oxygen (3.0 mg/L or 33% O2 saturation) treatments with our pH conditions in a factorial design to understand how pH, oxygen, and/or their interaction might influence metabolite abundances.

Results

Juvenile Dungeness crab general metabolome and lipidome composition

Untargeted lipidomics and metabolomics were carried out on 60 individual juvenile crabs exposed to ocean acidification treatments (Table 1 and Supplementary Table 1) from the end of megalopal stage through the entire duration of their first juvenile instar until 2 days after molting to their second juvenile instar (Supplementary Table 2). There are three main technologies used for untargeted metabolomics, liquid chromatography-mass spectrometry (LC-MS), gas chromatography-mass spectrometry (GC-MS), and nuclear magnetic resonance (NMR) spectroscopy. Because the goal of our study was to identify primary metabolites and complex lipids in the most comprehensive way, we employed GC-MS and LC-MS to general metabolites and lipids respectively. Both GC-MS and LC-MS have better sensitivity and resolution over NMR, with GC-MS having the advantage of identifying both soluble and volatile compounds and LC-MS having the advantage of identifying insoluble compounds20. From the lipidomics and metabolomics profiles generated by the West Coast Metabolomics Core using in-house open-access bioinformatics pipelines21,22 (see Methods for additional details), a total of 3113 lipids and 651 general metabolites were detected of which 88% (3320/3764 total detected compounds) had spectra that did not match LipidBlast23 or BinBase24 records and were classified as unknown compounds (Supplementary Table 3 and 4). Eighteen lipids were detected in fewer than 50% of individual profiles and were therefore excluded from all downstream analyses. The MS identification techniques found 29/281 metabolite classes among the 246 compounds detected and listed in Human Metabolome Database (HMDB)25 and 13/77 lipid classes among 195 compounds detected and listed in LIPID MAPS Structural Database26 (Supplementary Fig. 1). Among the largest classes represented in all identified compounds were LIPID MAPS classes glycerophosphocholines, triradylglycerols, and fatty acids, and HMDB classes organic acids, carboxylic acids, and organic oxygen compounds, consistent with previously observed Cancer magister body composition proportions of these compounds27,28.

Table 1 Summary of water chemistry for experimental treatments showing the mean ± standard deviation for each parameter.
Full size table

A combination of statistical methods reveal condition-specific responsive compounds

In general, individual metabolite and lipid profiles showed high variation regardless of treatment group, with no obvious trends revealed by clustering the relative abundance of all compounds by treatment groups (Supplementary Fig. 2). To identify specific compounds affected by treatments, we applied both univariate and multivariate statistical approaches, combining the strengths of different methods29. We did not detect a significant difference in molting times of crabs in different treatment groups (Supplementary Fig. 3) and therefore did not include molting time or stage duration as effects in our statistical models. To assess individual effects as well as interaction effects from pH and DO treatments on individual compounds, we used a two-way analysis of variance (ANOVA) on the processed abundance data (Supplementary Table 3 and 4). Prior to performing ANOVA, we verified compounds showed homogeneity of variances across treatment groups (>96% compounds showed a Levene test P value > 0.05), and that metabolite and lipid abundances were mostly (on average 69% compounds showed a Shapiro-Wilks test P value > 0.05) normally distributed within treatment groups (Supplementary Table 5). Although about one-third of compounds violated the ANOVA normality assumption, ANOVA can be considered robust to violations of this assumption when datasets have more than 10 samples per treatment group30,31. Of all compounds analyzed, 56/651 metabolites (including 24/160 known metabolites) and 98/3095 lipids (including 7/284 known lipids) showed overall model significance at P < 0.1 and at least one model term (pH, DO, and/or pH × DO interaction) significant at P < 0.05 (Supplementary Figs 4 and 5, and Supplementary Table 6). Because no metabolite or lipid overall model P values passed their respective 10% false discovery rate Benjamini-Hochberg P value threshold of 1 × 10−4 and less than 3 × 10−5, we chose to use their uncorrected P values.

To identify discriminatory compounds and model the relationship between metabolite profiles and exposure treatments using a multivariate linear regression approach, we used partial least square discriminatory analysis (PLS-DA). For facilitating comparisons of metabolite composition among treatment groups, compound abundances were centered around the mean and scaled by the reference group (ambient pH and DO) compound standard deviations prior to multivariate analysis32. The PLS-DA model for the metabolite data shows partial separation of low pH-treated samples from ambient pH-treated samples by the first and second components accounting for a total of 30.9% (Fig. 1a). Moreover, the metabolite data shows partial separation of low DO-treated samples in the second and third components, explaining an additional 6.8% of variation (Fig. 1b). For the lipid data, the PLS-DA model shows less clear treatment group separation on the first component (Fig. 1c) but shows partial separation of low pH-treated samples from ambient pH-treated samples on the second and third components explaining 14.7% of variation (Fig. 1d). These trends persisted even after removing stringently defined outlier samples (Supplementary Fig. 6a,b,f,g). Because the overall PLS-DA model had weak predictive power (Supplementary Fig. 6c–e,h–j), we used importance thresholds defined by the point of diminishing returns in the PLS-DA component loadings plots (Supplementary Fig. 7 and Methods), and identified a total of 45 metabolites (including 14/160 known metabolites) and 18 unknown lipids as compounds important in discriminating between treatment groups.

Figure 1
figure 1

Scores plots of PLS-DA analysis. Metabolite data show partial separation of low pH treatment groups (blue and purple) from low DO treatment groups (pink and green) by (a) Components 1 and 2 and (b) components 2 and 3. Lipid data show partial separation of treatment groups by (c) components 1 and 2 and (d) components 2 and 3. Shapes correspond to treatment groups (ambient pH and DO, Δ; ambient pH and low DO, +; low pH and ambient DO, x; low pH and low DO, ◊. Numbers refer to the sample number. Treatment group ellipses are the 95% pairwise confidence regions of component 1 and 2 variance76.

Full size image

Finally, we used a multivariate nonlinear-based supervised random forest classification to predict treatment classes from the metabolomic and lipidomic data (Supplementary Fig. 8a,b). A total of 9 metabolites (including 4/160 known metabolites) and 15 lipids (including 1/284 known lipids) were considered compounds important to predicting treatments based on how much their removal from the model decreased prediction accuracy (Supplementary Fig. 8c,d, Supplementary Table 8, and Methods). While no compounds were commonly identified by all 3 statistical methods as was expected given the fundamental differences of the methods, 7 metabolites were commonly identified by ANOVA and PLS-DA, 8 metabolites were in common between ANOVA and random forest, 1 metabolite overlapped between random forest and PLS-DA, and 4 lipids overlapped between ANOVA and random forest (Fig. 2a,b). This combination of analyses yielded a comprehensive set of compounds affected by treatments that no one method was capable of capturing on its own.

Figure 2
figure 2

Venn diagrams of all (a) metabolites and (b) lipids identified by ANOVA, PLS-DA, and/or random forest as significantly changed by exposure treatments.

Full size image

Low oxygen more strongly affects compound abundance than pH

Heatmaps of the statistically identified 94 (including 35 known) metabolites and 127 (including 7 known) lipids, show several compounds respond to low pH and oxygen treatments by commonly increasing or decreasing abundance relative to the ambient treatment (Fig. 3a,b). This is also summarized by violin plots of known compound abundances in Fig. 4a–c. For example, 5-methoxytryptamine, butyrolactam, cysteine, cysteine, homoserine, pipecolinic acid, piperidone, ribose, and xanthine commonly show a decrease in abundance in response to low DO treatments relative to ambient DO treatments. Glutamic acid, maltose, and maltotriose commonly show an increase in abundance in response to low DO treatments relative to ambient DO treatments (Fig. 4a). Overall, metabolite abundance tends to decrease rather than increase in response to low DO treatments relative to ambient DO treatments, exemplified by the mostly blue color of compound abundance averages for the low DO treatment group (green) in Fig. 3a. Specific to the low pH treatment, more metabolites show average abundances similar to ambient treatment signifying pH has a less dramatic effect on metabolites compared to low DO treatment (Figs 3a,b and 4a,b). This is also apparent in the combined low pH and low DO treatment, where metabolites and lipids show average abundances similar to the low DO treatment signifying low DO has a more dominant effect on compounds than low pH (Fig. 3a,b).

Figure 3
figure 3

Heatmap plots showing the average compound abundance (average peak intensity) for each treatment group for the (a) 94 general metabolites and (b) 127 lipids selected by multivariate and univariate statistical methods used to evaluate treatment effects. Individual known and unknown (“unk”) compound names are listed over the columns and treatment groups are listed over the rows (orange, ambient pH, ambient O2; green, ambient pH, low O2; purple, low pH, ambient O2; and blue, low pH, low O2). Compound average abundances are shown as auto-scaled within each compound (“relative abundance”).

Full size image
Figure 4
figure 4

Compound abundance levels (peak intensities) of known metabolites and lipids selected by univariate and multivariate statistical methods used to evaluate treatment effects. Violin plots with boxplot insets show the distribution of compound abundances within each treatment group for compounds statistically showing (a) DO treatment effect, (b) pH treatment effect, and (c) DO:pH interaction effect. Binbase names for each known compound are listed across the top of each plot. Abundance levels (peak intensities, noted as “quant”) are listed on the y-axis while treatments are listed along the x-axis. Statistical methods that each known compound was identified by are noted by shapes in the upper corners of each plot (ANOVA, triangle; PLS-DA, square; and random forest, diamond).

Full size image

Pathway analysis reveals potential mechanisms of pH and oxygen stress tolerance

To explore the physiological relevance of treatment-responsive compounds, we performed biochemical pathway analysis on all known compounds. While compounds classified as sugars and fatty acids were affected by the low pH and DO treatments, most of the significantly affected known compounds were part of the amino acids class, indicating amino acid metabolism was most significantly altered by the treatments. This is consistent with previous findings in NMR-derived metabolomes in adult green shore crabs exposed to varying degrees of hypercapnia33. Abbreviated biochemical pathway networks focusing on affected amino acids are shown in Fig. 5, and the complete biochemical pathway networks affected by low pH, low DO, and the combined low pH and DO treatments can be found in Supplementary Figs 9–11. The trends in amino acid abundances resulting from either low pH, low DO, osr the combined low pH and DO treatment (Fig. 5a–c), suggest different amino acid metabolic pathways are affected in response to each factor or combination of factors.

Figure 5
figure 5

Differential amino acid network effects from (a) low DO, (b) low pH, and (c) combined low pH and low DO. Compounds are clustered by chemical similarity. Node fill color is colored by Cohen D effect size comparing the treatment group to the ambient pH, ambient DO group, and node borders are colored by the treatment group mean fold change relative to the ambient pH, ambient DO group. Node shapes indicate the statistical method(s) from which the compound was classified as important. Node shape and label size are enlarged if the compound was identified as important by a statistical method. Gray edges indicate nodes sharing a KEGG biochemical pathway(s).

Full size image

Specific to low DO treatment response compared to the normoxic treatment, energy conservation pathways appear to be most affected (Fig. 5a). The increased lysine with decreased pipecolinic acid and the piperidine derivative piperidone abundance is suggestive of downregulation of the lysine degradation pathway34,35. The decreased abundance of both cysteine and the cysteine homodimer, cystine, suggests that cysteine catabolism is upregulated. This coupled with the increased abundance of glutamic acid suggests the glutathione synthesis pathway could be downregulated, as both cysteine and glutamic acid are precursor molecules in glutathione synthesis36,37. Glutathione itself was not detected, but is typically challenging to detect by MS techniques in marine animals due to its reactivity38. The glycogen intermediates maltose and maltotriose show an increase in abundance consistent with previously observed glycogen synthesis pathway upregulation during hypoxic stress39 and potentially linked to gluconeogenesis upregulation in crabs exposed to oxidative stress40. Purine and pyrimidine metabolic intermediates orotic acid, ribose, and xanthine show decreased abundance suggesting purine and pyrimidine metabolic pathways could be downregulated41.

Specific to low pH treatment response compared to ambient pH treatment, ATP generation pathways appeared most affected (Fig. 5b). The decreased abundance of aspartic acid and the subtle decreased abundance of maleic acid may suggest the citric acid cycle intermediates that they form (oxaloacetate, malic acid, and fumaric acid)42,43 are favored. Although, oxaloacetate was not detected (likely due to the instability of the alpha keto acid compound)44, and malic acid and fumaric acid did not show a significant difference across treatments. However, citric acid did show an increase in abundance. Taken all together, the altered abundance of these compounds suggest that citric acid cycle activity could be upregulated in response to low pH45. An increase in glutaric acid in response to low pH suggests the catabolism of its parent molecule, glutaryl CoA, which produces ATP in addition to glutaric acid46. Glutaryl CoA catabolism also supports citric acid cycle activity because glutaryl CoA can uncompetitively inhibit the alpha-ketoglutarate dehydrogenase complex that facilitates a rate limiting step in the citric acid cycle47.

Among the compounds affected by the combined low pH and DO treatment, phosphoethanolamine and phophatidylethanolamine (PE(p-34:2) or PE(o-34:3)) show increased abundance compared to ambient treatment (Supplementary Fig. 11), suggesting that low pH and DO might alter the glycerophospholipid synthesis pathway of which phosphatidulethanolamine is a product and phosphoethanolamine is a substrate intermediate in a rate limiting step48. Alanine shows a decrease in abundance, suggesting that alanine synthesis is decreased and that its substrate for synthesis, pyruvate, may be limited49 and/or that alanine catabolism, which has previously been observed in crabs under oxidative stress40 may be occurring. The increased abundance of cystathionine and subtle increase of methionine suggest cysteine and homoserine synthesis are likely stalled50 in response to combined low pH and low DO exposure. This overlaps with the decrease in homoserine abundance and decrease in cysteine synthesis products (cysteine and cystine) observed in low DO treatment (Fig. 5b,c).

Discussion

We used metabolomics to explore how the Dungeness crab might respond to the simultaneous change in oxygen and carbon dioxide in predicted future climate change scenarios in the Pacific Northwest. Although the untargeted metabolomics and lipidomics approaches rapidly identified hundreds of different molecular species, the majority of compounds identified in our study, including most compounds selected by statistical methods, lack annotations in existing metabolite databases. Thus, our understanding of their role in response to low pH and low DO treatment is limited without extensive further validation. Still, through shedding light on the simultaneous activity of hundreds of known compounds, we were able to observe a dynamic range of metabolic responses among individuals within treatment groups that indicates that Dungeness crab have flexibility in how their biochemistry compensates for environmental change. While we attempted to control for variation between individuals by collecting animals from the same location within a 2-month period, we were unable to control for prior environmental exposure or genetic background of the wild-caught animals we used. Where prior studies have found high genetic variation among individuals within one sampling site51,52, we suspect that both genetics and prior environmental exposure likely contributed to the high variation in metabolite abundances among individuals within treatment groups, which may have obscured the different treatment effects. Within-sample variation could be reduced by sampling from a specific tissue rather than using whole animals, however this may be challenging in early stage juvenile Dungeness crabs given their small size.

We applied three different statistical methods (ANOVA, PLS-DA, and random forest) for selecting important metabolites in order to combine the strengths of powerful univariate and multivariate analyses29. Although using univariate statistics in the strictest sense with an FDR correction to control for false discoveries showed no statistically significant variation for compounds across treatments, not all compounds were identified with the same confidence level and were subject to the sensitivity of the mass spectrometry method used. Applying an FDR correction to all detected compounds assumes that all compounds had equal chance for discovery, which can mask important biology in highly variable datasets53. For this reason and due to our use of wild-caught whole-animal samples, we chose a more liberal focus on all compounds showing an overall ANOVA model P value < 0.1 and at least one model term P value of <0.05, or identified as important in PLS-DA or random forest models in our pathway analysis.

While we acknowledge a level of uncertainty in our pathway analysis due to liberal selection of important treatment-responsive compounds, the resulting proposed affected pathways are consistent with previous observations of pH and hypoxia effects on different organisms. In general, amino acid metabolism is a well-documented mechanism for stress tolerance54,55. This class of compounds contains versatile chemical structures that serve as buffering molecules, antioxidants, signaling molecules, and chemical building blocks for the synthesis of proteins important in stress response (i.e., heat shock proteins, unfolded protein response proteins, ion channels). Under low oxygen conditions, it is in the best interest of the animal to limit non-essential energy consuming pathways56. One way Dungeness crab may do this is through reducing the activity of evolutionarily conserved gamma-glutamyl cycle that synthesizes glutathione and consumes ATP57,58 by limiting cysteine availability through catabolism. Interestingly, glutathione reduction in response to hypoxia has been observed in multiple mammalian cell lines59,60,61,62. Also under low oxygen conditions in nature, food can be scarce, and it is theorized that glycogen storage during low oxygen can help prepare cells for low nutrient conditions. In multiple mouse and human cancer cell lines39,63,64, increased glycogen storage has been observed in response to hypoxia, and is induced by highly evolutionarily conserved hypoxia-inducible factor transcriptional signalling39. It seems that Dungeness crab may also adopt this strategy in combating low oxygen conditions. Similar to the metabolomes observed in the Hammer et al. (2012) study of green shore crab exposed to hypercapnia33, we also observed slight increase in the branched chain amino acids (BCAAs) leucine and isoleucine in the low pH group, which may reflect increased osmo- and/or pH regulation given BCAAs are implicated in osmotic stress signaling65 and that osmo- and pH regulation are coupled in decapods66,67. However, unlike Hammer et al.33, we did not observe differences in lactate, taurine, sarcosine, or glycine between our control or low pH samples, possibly due to the inherent differences between our studies (i.e., species, life stage, pH, exposure length, metabolomics methods, etc.). Under low pH, adult Dungeness crab initially develop hypercapnia which then abates over time via elevated hemolymph bicarbonate likely generated through gill restructuring and the upregulation of energy consuming ion-exchange proteins12. Our results indicate that this process may occur in early juveniles in low pH conditions given that we found metabolite profiles that support energy generation via increased citric acid cycle activity.

Having parsed out effects from low oxygen and low pH, we found that in combined low oxygen and pH conditions low oxygen has a more dramatic effect on metabolite abundance. However, the animals in this treatment group still showed pathway alterations similar to the low pH treatment group including increased abundance of citric-acid-cycle-related metabolites to potentially increase ATP production. It is not yet clear what the longer-term consequences are of these metabolic adjustments or how long these responses could be sustained. Future avenues of research to expand on these findings should include targeted metabolomics to confirm the compounds identified in this study as well as capture more antioxidants and citric acid cycle intermediates that were not detected in this study. Targeted expression profiling would also be helpful in confirming the biochemical pathway activity predictions from this study. Ultimately, longer-term exposure experiments on crabs reared over generations might best reveal the maximal duration that these metabolic responses can be sustained and any long-term consequences of low pH and oxygen exposure. This exploratory metabolomics and lipidomics analysis uncovered potential biochemical pathways affected by experiments simulating ocean acidification and hypoxia, and can now serve as preliminary hypotheses for deeper investigations of how the Dungeness crab (and even other crustaceans) may tolerate global ocean change.

Methods

Animals

Cancer magister megaolopae were collected from a single site in Puget Sound (47.950232, −122.301784) on several days between the June and September 2016 using light traps that were set overnight. The contents of each trap were immediately transferred into a 5-gallon bucket that was, within 5 minutes, pooled into a cooler with ice packs and an air bubbler. After transferring contents from 7 traps in approximately 1 hour, the cooler containing megalopae was brought within 5 minutes into the lab and megalopae were individually transferred onto the water flow-through system described below. The total time for transferring megalopae from the cooler onto the seawater flow-through system was about 2 hours.

CO2 exposure experiments

Crabs were held in individual 250 mL customized jars on Mobile Ocean Acidification Treatment Systems (MOATS). These systems flowed one-micron-filtered, UV-sterilized, Puget Sound seawater maintained at 12 °C. Prior to flowing through jars, seawater was degassed and oxygen, nitrogen and carbon dioxide gases were resupplied to finely control dissolved gas levels. Temperature, pH, and dissolved oxygen were continuously monitored throughout the duration of the experiment by Omega thermistors, Honeywell Durafet III probes, and Vernier optical dissolved oxygen probes, respectively. The pH was additionally validated by periodic sampling of water for dissolved inorganic carbon and total alkalinity, and by bi-weekly spectrophotometric pH measurements using an Ocean Optics USB 230 2000+ Fiber Optic Spectometer with SpectraSuite software and a 5 mM solution of Sigma Aldritch m-cresol purple indicator dye. Both Durafet probes and spectrophotometric pH method have previously been shown to have the level of precision required for the pH treatments targeted68,69. MOATS seawater pH and dissolved oxygen were automatically adjusted through a data-driven feedback system. Megalopae were fed Artermia salina (San Francisco Bay brand) at a target concentration of 1 nauplius per mL every 3 days. Once megalopae transitioned to the first juvenile instar, they were fed small pieces of squid. Upon transitioning to the second juvenile instar, crabs were held on the MOATS for an additional 48 hours post-molt in attempt to reduce variation due to potential stochastic physiological processes associated with molting70,71. Juveniles were then immediately frozen and stored at −80 °C after lightly blotting with a paper towel to remove excess seawater. To reduce variation from length of exposure, which in total ranged from 20–65 days, only 60 juvenile crabs with the average exposure time (30–33 days) were chosen for metabolomics analysis with 15 individuals from each treatment group.

Water chemistry analysis

Mean and standard deviation calculations for pH, DO, and temperature (Table 1) were based on logger data (measured every 10 minutes) with the exception of pH and DO data from MOATS 5 and 12, where Spectrophotometric pH and Presens DO measurements were instead used due to incorrectly calibrated pH and DO probes in these MOATS. Temperature logger readings above 25 °C and below 4 °C were excluded from the mean and standard deviation analysis because these were out of the achievable range for the equipment used and were indicative of a thermistor malfunction. Total alkalinity (TA) and salinity were discrete measurements, with TA ranging from 2003.4 to 2036.1 μmol/kg and salinity ranging from 29.42 to 29.96 ppt throughout the course of the experiment. The complete data table is included as Supplementary Table 1.

Sample preparation

Frozen juvenile crabs were sent to the West Coast Metabolomics Center, Davis, CA for sample preparation, metabolomics and lipidomics profiling. Whole animals were thawed, weighed, and derivatized as previously described72,73 (Supplementary Table 2). Briefly, samples were extracted at −20 °C with 2 mL of degassed acetonitrile/isopropanol/water (3:3:2) solution and solvents were evaporated to complete dryness with a Labconco Centrivap cold trap concentrator. Membrane lipids and triglycerides were subsequently removed from dried samples with 50% acetonitrile, and samples were again concentrated to complete dryness. 15 mg of each sample preparation was used for metabolomic profiling. For lipidomic profiling, 15 mg of the same sample preparation was also used to which internal standards, C8-C30 fatty acid methyl esters were added. Aliquoted samples were derivatized with methoxyamine hydrochloride (Sigma-Aldrich) in pyridine (Acros Organics) and then by N-methyl-N-(trimethylsilyl) trifluoroacetamide (Sigma-Aldrich) for trimethylsilylation of acidic protons.

Metabolite and lipid data acquisition

General metabolite and lipid abundances were quantified from derivatized samples by gas-chromatography, time-of-flight mass spectrometry (GC-TOF/MS) and charged-surface, hybrid-column, electrospray-quadrupole, time-of-flight mass spectrometry (CSH-ESI QTOF MS/MS), respectively. For metabolites, an Agilent 6890 gas chromatograph (Santa Clara, CA) was used with a Leco Pagasus IV time-of-flight mass spectrometer running Leco ChromaTOF software 2.32 (St. Joseph, MI). The following temperature profile was used: 50 °C to 275 °C final temperature at a rate of 12 °C/s and hold for 3 minutes. Injection volume was 0.5 μl with 10 μl/s injection speed on a splitless injector with a purge time of 25 seconds. Liner (Gerstel #011711-010-00) was changed after every 10 samples, (using the Maestro1 Gerstel software vs. 1.1.4.18). Before and after each injection, the 10 μL injection syringe was washed 3 times with 10 μL ethyl acetate. For gas chromatography, a 30 m long, 0.25 mm i.d. Rtx-5Sil MS column (0.25 μm 95% dimethyl 5% diphenyl polysiloxane film) with additional 10 m integrated guard column was used (Restek, Bellefonte PA). 99.9999% pure Helium with a built-in purifier (Airgas, Radnor PA) was set at a constant flow of 1 mL/minute. The oven temperature was held constant at 50 °C for 1 minute and then ramped at 20 °C/minute to 330 °C at which it is held constant for 5 minutes. The transfer line temperature between gas chromatograph and mass spectrometer was set to 280 °C. Electron-impact ionization at 70 V was employed with an ion source temperature of 250 °C. Acquisition rate was 17 spectra/second, with a scan mass range of 85–500 Da.

For positively charged lipids, an Agilent 6530 QTOF mass spectrometer with resolution 10,000 was used and for negatively charged lipids, an Agilent 6550 QTOF mass spectrometer with resolution 20,000 was used. Electrospray ionization was used to ionize column elutants in both positive and negative modes. Compounds were separated using a Waters Acquity ultra-high-pressure, liquid-chromatography charged surface hybrid column C18 (100 mm length × 2.1 mm internal diameter; 1.7 um particles) using the following conditions: mobile phase A (60:40 acetonitrile:water + 10 mM ammonium formiate + 0.1% formic acid, mobile phase B (90:10 isopropanol:acetonitrile + 10 mM ammonium formiate + 0.1% formic acid), 65 °C column temperature, a flow rate of 0.6 mL/minute, an injection volume of 3 uL, an injection temperature of 4 C, and a gradient of 0 minutes 15%, 0–2 minutes 30%, 2–2.5 minutes 48%, 2.5–11 minutes 82%, 11–11.5 minutes 99%, 11.5–12 minutes 99%, 12–12.1 minutes 15%, and 12.1–15 minutes 15%. The capillary voltage was set to + 3.5 and −3.5 kV, and the collision energy to 25 or 40 eV for positive and negative modes. Mass-to-charge ratios (m/z) were scanned from 60 to 1700 Da and spectra acquired every 2 seconds.

Spectral data processing

Acquired metabolite data were processed using UC Davis’s BinBase workflow, which performs data processing including peak detection at signal-to-noise levels of 5:1 throughout the chromatogram. Resulting apex masses are reported for use in the BinBase algorithm to facilitate metabolite identification and quantification74. Metabolites were identified through comparison to the BinBase database24 and peak heights were normalized to total metabolite content75.

For acquired lipid data, MassHunter (Qual v. B05.00) was used to find peaks from the raw data in up to 300 chromatograms. These peaks were imported into Mass Profiler Professional for alignment to determine which peaks occur in at least 30% of the chromatograms. These peaks were then quantified with MassHunter. Resulting accurate mass data and tandem MS/MS spectra were compared to LipidBlast libraries for compound identification23. All spectra mapping to previously identified compounds are accessible in the public database MassBank of North America (http://mona.fiehnlab.ucdavis.edu/). Spectra that did not map to previously identified compounds are accessible by searching their BinBase identifiers listed in Supplementary Tables 3 and 4 in the BinVestigate database (http://binvestigate.fiehnlab.ucdavis.edu).

Statistical analyses

All statistical analyses excluded 18 compounds that were detected in less than half of the individual crabs surveyed, and were carried out in R, except for PLS-DA and random forest analyses which were carried out using the MetaboAnalyst web interface.

Univariate statistics

Prior to applying a univariate statistical test, normality and heteroscedasticity of compound abundances within treatment groups was first assessed using the Shapiro-Wilks test and the Levene test, respectively (Supplementary Table 5). After confirming normality and heteroscedasticity assumptions were correct, a two-way ANOVA was then applied to each compound. A Benjamini-Hochberg FDR correction was applied to ANOVA P values to correct for multiple testing, but corrected P values were not used in selecting important compounds. Metabolites and lipids were selected for pathway analyses if they had an overall model ANOVA P value less than 0.1 and an effect P value with less than 0.05 without an FDR correction applied. Results from ANOVA tests are reported in Supplementary Table 6.

Multivariate statistics

Prior to applying multivariate testing, data were normalized by mean-centering and scaling by the reference-group standard-deviations of each compound in order to make compounds more comparable to one another32. Normalized data was then uploaded to MetaboAnalyst76 with no further normalizations applied. The PLS-DA function was run with default settings, and score plots were generated that showed sample distribution across pairwise components along with the 95% component-pairwise confidence region for samples in each treatment (Fig. 1). Next, the effect of outlier samples (samples falling outside the 95% component-pairwise confidence region) was checked. Score plots (Fig. 5) were first checked for common outlier samples across all plots. Since no common outlier samples were detected across all plots, outlier sample effects were secondarily checked by redefining outlier samples more stringently (samples falling outside or near the 95% component-pairwise confidence region in any plot). After removing stringently redefined outliers (metabolite samples 17, 20, and 56 and lipid samples 6, 11,16, 20, 37, 38, and 43), rerunning the PLS-DA function, and observing the same treatment group separation trends (Supplementary Fig. 6a,b, f,g), the effect from outlier samples was determined to minimal so all samples were retained. A table of component loadings for each score plot from the initial PLS-DA analysis (Fig. 1) was exported from MetaboAnalyst (Supplementary Table 7). Compound loadings were ordered by from greatest variance explained to lowest, and plotted for PLS-DA components. Importance thresholds were placed at the point of diminishing returns in the plot curves (Supplementary Fig. 7). Specifically for metabolites, loadings were plotted PLS-DA components 1, 2, and 3 (Supplementary Fig. 7a–c) since these PLS-DA components gave the greatest separation between treatment groups (Fig. 1a,b). The points of diminishing returns that importance thresholds were placed in the metabolite loadings plots were as follows: component 1, loadings threshold > 0.05; component 2, loadings threshold > 0.1; and component 3, loadings threshold > 0.05. Specifically for lipids, loadings were plotted for PLS-DA components 2 and 3 (Supplementary Fig. 7d,e) since these gave the best separation between treatment groups, while component 1 could only distinguish the ambient from the altered treatment groups (Fig. 1c,d). The points of diminishing returns that importance thresholds were placed in the lipids loadings plots were as follows: component 2, loadings threshold > 0.05; and component 3, loadings threshold > 0.05. For the random forest analysis, the Random Forest function was run with 2000 decision trees and either 25 predictors for the metabolite dataset or 56 predictors for the lipid dataset, conforming to the default classification value being the square root of the number of variables77. The complete table of important variables and their mean decrease in random forest accuracy prediction exported from MetaboAnalyst (Supplementary Table 8). The mean decrease in prediction accuracy for each compound was ordered from largest to smallest and plotted (Supplementary Fig. 8c,d). Importance thresholds were drawn from the plots at the points of diminishing returns, which for metabolites was a mean decrease in accuracy > 0.001 and for lipids was a mean decrease in accuracy > 0.00047.

Heatmaps in Fig. 3 and Supplementary Fig. 2 were made using the heatmap function in the MetaboAnalyst web interface with Euclidean distance, ward clustering, plotting auto-scaled compound abundance averages for each treatment group.

Pathway analysis

Pubchem and KEGG identifiers were obtained for metabolites and lipid international chemical identifier keys provided by the West Coast Metabolomics Core using Chemical Translation Service78 and Pubchem Identifier exchange79. Metamapp80 was then used to map metabolites by chemical and biochemical relationships. The generated SIF file (Supplementary File 1) was imported into Cytoscape81. For visualizing the pH, DO, and pH:DO interaction effects, both mean abundance fold change and Cohen D effect size values were calculated for each compound within a treatment group relative to the ambient treatment group. Cohen D effect size \((\frac{{\mu }_{treatment}-{\mu }_{ambient}}{{S}_{pooled}})\) was calculated for each compound processed abundance values within a treatment group relative to the ambient treatment group using the R package effsize82. The treatment group −log2 fold change values and the Cohen D effsize output from R were added to a node attribute file (Supplementary Table 9) that was also imported into Cytoscape. Low pH, low DO, and low pH:low DO networks were styled using the following settings: node shapes were set according to statistical method used to select the compound (ANOVA, triangle; PLS-DA, square; random forest, diamond; both PLS-DA and ANOVA, hexagon; not selected as important by any method, circle); node border was color by the treatment group mean compound abundance fold change relative to ambient (1.25X–2X fold change, yellow gradient; −1.25X–1.25X fold change, grey; −1.25X–10X fold change, blue gradient); node fill was colored by the Cohen D effect size of treatment on compound abundance relative to ambient (effect size of 0.5–1, yellow gradient; −0.5–0.5, white; and −0.5–1, blue gradient); node shape size and label size were set to fixed values (150 height and width with size 50 font for statistically selected compounds, and 30 height and width with size 12 font for compounds not identified by statistics (Supplementary Fig. 9–11).