Pacific geoduck aquaculture is a growing industry, however, little is known about how geoduck respond to varying environmental conditions, or how the industry will fare under projected climate conditions. To understand how geoduck production may be impacted by low pH associated with ocean acidification, multi-faceted environmental heterogeneity needs to be included to understand species and community responses. In this study, eelgrass habitats and environmental heterogeneity across four estuarine bays were leveraged to examine low pH effects on geoduck under different natural regimes, using targeted proteomics to assess physiology. Juvenile geoduck were deployed in eelgrass and adjacent unvegetated habitats for 30 days while pH, temperature, dissolved oxygen, and salinity were monitored. Across the four bays, pH was lower in unvegetated habitats compared to eelgrass habitats. However this did not impact geoduck growth, survival, or proteomic abundance patterns in gill tissue. Temperature and dissolved oxygen differences across all locations corresponded to differences in growth and targeted protein abundance patterns. Specifically, three protein abundance levels (trifunctional-enzyme β-subunit, puromycin-sensitive aminopeptidase, and heat shock protein 90-⍺) and shell growth positively correlated with dissolved oxygen variability and inversely correlated with mean temperature. These results demonstrate that geoduck may be resilient to low pH in a natural setting, but other abiotic factors (i.e. temperature, dissolved oxygen variability) may have a greater influence on geoduck physiology. In addition this study contributes to the understanding of how eelgrass patches influences water chemistry.

Protein Preparation: Relative protein abundance was ultimately assessed in a two-phase proteomics approach using Selected Reaction Monitoring (SRM), with targets identified using Data Independent Acquisition (DIA). Tissues were prepared separately for DIA and SRM, both following the protocol in Timmins-Schiffman et al. (2014) with a few exceptions. Tissue was homogenized with sterile plastic pestle in 100 µl lysis buffer (50 mM NH4HCO3, 6M urea solution) and sonicated with Sonic Dismembrator (Fisher Scientific, Model 120) at 50% amplitude for ten seconds, three times. Protein concentration was quantified via Pierce™ BCA Protein Assay Kit (ThermoFisher Scientific, Waltham, MA USA). Mini-Trypsin Digestion: Aliquots of protein (30 µg for DIA, 100 µg for SRM) were suspended in Lysis Buffer (50 mM NH4HCO3 + 6 M urea solution) to a total volume of 100 µl followed by: 1) a 1 hour incubation at 37°C with 200 mM Tris(2-carboxyethyl)phosphine (2.5µl) and 1.5 M Tris at pH 8.8 (6.6 µl); 2) 1 hour at room temperature in dark with 200 mM iodoacetamide (20 µl); 3) 1 hour at room temperature with 200 mM diothiothreitol (20 µl); 4) 1 hour at room temperature with 2 ug/µl Lysyl Endopeptidase (Lys-C, Wako Chemicals) (3.3 µg); 5) overnight at room temperature in 25 mM NH4HCO3 (800 µl) + high pressure liquid chromatography grade methanol (200 µl) + Pierce Trypsin Protease, MS Grade (1 ug/µl, Thermo Scientific) at 1:30 enzyme:protein ratio (3.3 µg). Samples were evaporated to near dryness at 4°C using a CentriVap Benchtop Vacuum Concentrator. Desalting: Samples were desalted to isolate peptides using MacroSpin Columns (Nest Group, 50-450 µl, Peptide Protein C18). Peptides were reconstituted in 5% acetonitrile + 0.1% trifluoroacetic acid (TFA) (100 µl), then 10% formic acid (70 µl) was added to achieve pH ≤2. Columns were washed with 60% acetonitrile + 0.1% TFA (Solvent A, 200 µl) four times, then equilibrated with 5% acetonitrile + 0.1% TFA (Solvent B, 200 µl) three times. Peptides were bound to the column by running the digest through the column twice, followed by peptide elution with two additions each of Solvent A (100 ul). Columns were spun for 3 minutes at 3000 rpm on VWR Galaxy 16DH digital microcentrifuge at each stage. Samples were evaporated to near dryness at 4°C, then reconstituted in the Final Solvent (3% acetonitrile + 0.1% formic acid) (60 µl for 0.5 µg/µl final concentration of protein, and 50 µl for 2 µg/µl final concentration for DIA & SRM, respectively). Peptide sample preparation and internal standard: Final mixtures for mass spectrometry included 3.13 fmol/µl Peptide Retention Time Calibration mixture (PRTC), 0.33 µg/µl and 0.5 µg/µl peptides for DIA and SRM, respectively, in Final Solvent for 15 µl total volume. To confirm that peptides were quantified correctly in SRM, 10 µg from 5 randomly selected geoduck peptide samples were pooled, and 8 dilutions were prepared by combining with oyster peptides at known percentages of total protein content (10%, 13.3%, 20%, 40%, 60%, 80%, 87.7%, 90%) and analyzed alongside other samples. DIA was performed to assess global protein abundance patterns and to identify consistently detectable peptides for SRM targets. Eight samples, one per deployment location, were analyzed in technical duplicates via liquid chromatography tandem mass spectrometry (LC-MS/MS) with the Thermo Scientific™ Orbitrap Fusion Lumos™ Tribrid™. Prior to sample analysis, the 30 cm analytical column and 3 cm trap were packed in-house with with C18 beads (Dr. Maisch HPLC, Germany, 0.3 μm). For each sample, 3 µl of geoduck peptides (1.0 µg) + PRTC peptides was injected and analyzed in MS1 over 400-900 m/z range, in 5 m/z overlapping isolation windows from 450-850 m/z with 15K resolution in MS2. Final Solvent blanks were run between each geoduck peptide injection to ensure against peptide carry-over. Lumos MS/MS method and sequence files are available in the project repository (GitHub). SRM samples were analyzed on a Vantage Triple-Stage Quadrupole Mass Spectrometer (Thermo Scientific, San Jose, CA, USA), and injected by a nanoACQUITY UPLC® system (Waters, Milford, MA, USA) at random in two technical replicates. For each sample, 2 µl of peptides + PRTC solution containing 1.0 µg of geoduck peptides was injected, trapped on a 3 cm pre-column and separated on a 30 cm analytical column using a chromatography gradient of 2-60% acetonitrile over 60 minutes. Columns were prepared as in DIA (above). Samples were injected in randomized groups of 5, followed by a Peptide Retention Time Calibration (PRTC) plus bovine serum albumin peptides (BSA) standard, then Final Solvent blank. Vantage MS sequence and method files are available in the project repository (GitHub). DIA: Proteins were inferred using an assembled, translated, and annotated P. generosa gonad transcriptome (combined male and female) (Timmins-Schiffman et al. 2017, BioProject Accession: PRJNA316216). Transcriptome peptides were queried against those detected by the Lumos MS/MS using PEptide-Centric Analysis (PECAN) (Ting 2017) to create a peptide spectral library (.blib type file). DIA raw files were first demultiplexed using MSConvert (ProteoWizard version 3.0.10738, 2017-04-12) (Chambers et al. 2012) with filters set to vendor centroiding for msLevels [2,3] ( --”peakPicking true 1-2”), and optimizing overlapping spectra (“demultiplex optimization-overlap only”). The transcriptome fasta file was tryptic digested in silico in Protein Digestion Simulator (version 2.2.6471.25262), set to Fully Tryptic from 400-6000 fragment mass range, 5 minimum residues allowed, 3 maximum missed cleavages and peak matching thresholds set to 5 ppm mass tolerance, and 0.05 ppm NET tolerance. Skyline version 3.7.0.11317 (MacLean et al. 2010) automatically selected transition peaks and quantified peptide abundances using peak area integration. All PRTC peptide peak selections were manually verified and corrected. Skyline peak selection error rate was calculated by manually checking chromatograms from 100 proteins across all DIA samples. Auto-selected peaks were assigned correct or incorrect selection based on transition retention time alignment across replicates, using PRTC peptides as a reference. Transition peak area, which is assumed to correlate to peptide transition abundance, was exported from Skyline for analysis in R version 2.4-3 (R Core Team 2016). Abundance was normalized by the total ion current (TIC) for each injection. Thirteen proteins were selected for SRM targets (Table 2). First, candidate proteins (~200) from DIA were selected based on biological function listed in the Universal Protein Knowledgebase (Apweiler et al. 2004) and evidence of stress response in bivalves from the scientific literature. Candidate proteins were screened for detectability using DIA results. Selected proteins were required to have ≥3 high quality peptides, each with ≥3 transitions, present in all DIA biological and technical replicate data. Quality peptides had uniform peak morphology and retention time in Skyline across replicates. A total of 49 peptides were selected for SRM: 39 to quantify 13 target proteins (116 transitions), and 10 for internal standard (30 transitions). A full list of transition targets are published on panoramaweb.org and available in the project repository (GitHub). SRM: Peptides were identified and quantified via Skyline-daily version 3.7.1.11357 (MacLean et al. 2010). Raw SRM files were imported into a Skyline-daily project along with the target protein transitions, and the spectral library (.blib file) created previously in the DIA Protein Identification step. SRM peptides were verified by regressing PRTC peptide retention time (RT) in SRM against retention time in DIA. A fitted model from PRTC peptides predicted RT of protein target peptides. Where necessary, peak selection and boundaries were manually adjusted for SRM peptide chromatograms, and actual RT were regressed against predict RT to confirm correct selection (F(1,38): 5768, p-value: < 2.2e-16, Adjusted R-squared: 0.9933) (Supplemental Figure 5). Transition peak area, defined henceforth as abundance, was exported from Skyline for further analysis in R (R Core Team 2016). Abundance results from the separate serial dilution samples were used to remove peptides that did not adhere to the dilution curve. Briefly, dilution abundances (exported from Skyline) for each transition were normalized by the most dilute sample abundance, then regressed against predicted ratios. Peptides with slope coefficient 0.2<x<1.5 and adjusted R2 >0.7 were included in analysis. Ten of the 39 peptides were discarded from the dataset based on dilution standards results (Supplemental Figure 6). To determine and remove disparate technical replicates, NMDS analysis was performed as described above. Technical replicates with ordination distance >0.2 were removed, and only samples with two technical replicates were preserved for analysis (Supplemental Figure 7). Thirteen technical replicates from different samples and all replicates from three sample were discarded, for 84% technical replicate and 94% biological replicate retention. Within samples, transitions with coefficients of variation (CV) > 40% between technical replicates were also discarded (2% of all transitions across 21 samples). In final dataset for differential analysis, 10 proteins, 26 peptides, and 77 transitions were retained. Mean transition abundance was calculated for replicates, with zero in the place of n/a values, which Skyline generates for replicates without peaks. Transition abundances within each peptide were summed for a total peptide abundance before analyzing for differential abundance. Data are available via ProteomeXchange with identifier PXD012266

Juvenile geoduxk were collected during low tide and transferred on wet ice to shore where mortality and size were recorded. Live animals were dissected, and ctenidia tissue was isolated and flash-frozen in an ethanol-dry ice bath. Ctenidia was selected for proteomic analysis due to its direct interaction with the environment, importance in gas and ion regulation, and its implication in environmental stress response (Timmins-Schiffman et al., 2014, Matozzo et al. 2013, Zhang et al. 2015, Thompson et al. 2015). During sampling all instruments were sterilized between samples with bleach then ethanol, and rinsed with nanopure water. Samples were held on dry ice while transported back to the lab and stored at -80ºC. For DIA, 8 ctenidia tissue samples were analyzed, one sample from each location and habitat: FB-eelgrass (G048), FB-unvegetated (G058), PGB-eelgrass (G077), PGB-unvegetated (G068), CI-eelgrass (G010), CI-unvegetated (G018), WB-eelgrass (G131), WB-unvegetated (G119). For SRM, new ctenidia samples were examined, 12 individuals per bay (Fidalgo Bay, Port Gamble Bay, Case Inlet, Willapa Bay), with 6 from each habitat (eelgrass, unvegetated) for 48 samples total.