Proteomics of hyposaline stress in blue mussel congeners (genus Mytilus): implications for biogeographic range limits in response to climate change

Biogeographic distribution ranges of marine organisms are shifting due to climate change, specifically rising atmospheric and oceanic temperatures, increasing acidity of the ocean and more frequent and extreme precipitation events leading to greater hyposaline stress in estuaries and coastal waters (Harley et al., 2006; IPCC, 2007; Min et al., 2011; Pall et al., 2011). To assess which environmental stressor, either in isolation or combination, will affect the physiology of marine organisms the most and thus be the driving force for range shifts, we have to assess the physiological impacts of thermal, pH and hyposalinity stressors. The realization that extreme precipitation events may be a potential driving force for range shifts gives this research topic renewed urgency to improve our predictions of the ecological impacts of climate change.

It is now evident that rising temperatures affect rates of physiological processes and the integrity of the cell's macromolecular structure, and thereby contribute to shifting range limits (Hochachka and Somero, 2002; Pörtner, 2010; Tomanek, 2008; Tomanek, 2010). Although more extreme precipitation events due to higher atmospheric humidity levels associated with climate change have been documented (Groisman et al., 2005; Min et al., 2011), biologists are only now starting to evaluate the potential impacts of these events, e.g. greater levels of hyposaline stress, on species distribution ranges (Levinton et al., 2011). Extreme precipitation events will occur in a warmer world even if total precipitation levels do not increase (Karl and Trenberth, 2003). Analyses of regional past trends and projected future scenarios of precipitation and stream flow under different climate scenarios and their potential biological impacts are available for Chesapeake Bay. These suggest that winter flow will increase but summer flow will decrease, with an overall increase of acute hyposaline stress conditions (Najjar et al., 2010). An analysis of precipitation trends for the USA predicts an increase in extreme precipitation events for some coastal regions in California (Groisman et al., 2005), but does not state whether that will lead to heavier river flow rates.

To assess the effect of extreme precipitation events and their potential impacts on shifting distribution ranges, we decided to investigate the physiological responses to hyposaline stress of a pair of blue mussel species whose recent biogeographic changes have been documented and linked to changes in both temperature and salinity. One of the two blue mussel species is Mytilus galloprovincialis, which invaded southern California during the middle of the last century and has replaced the native M. trossulus from the southern part of its distribution range, from Baja California to central California (Braby and Somero, 2006a; Geller, 1999; McDonald and Koehn, 1988; Rawson et al., 1999). Although the range limits of these congeners are still in flux due to shorter climatic variations, e.g. the Pacific Decadal Oscillation, the main hybrid zone ranges roughly from Monterey Bay to San Francisco Bay, with small numbers of M. galloprovincialis hybrids found further north to Humboldt Bay (Braby and Somero, 2006a; Hilbish et al., 2010). Field surveys indicate that the distribution within the hybrid zone is determined by both temperature and salinity (Braby and Somero, 2006a; Schneider and Helmuth, 2007). Salinity seems to play a crucial role because M. trossulus occurs at sites with higher freshwater input that are warm enough to normally favor occurrence of the more warm-adapted M. galloprovincialis (Braby and Somero, 2006a). Based on their natural distribution, the Eastern Pacific M. trossulus seem to prefer colder temperatures and tolerate lower salinity levels, whereas the Mediterranean M. galloprovincialis is a warm-water species that prefers high salinity levels (Seed, 1992). Measurements of growth, heart rate and survival generally confirm these interspecific differences (Braby and Somero, 2006b; Schneider, 2008). One hypothesis for the underlying mechanistic differences is that M. trossulus may achieve tolerance to lower salinities by closing their shells, as indicated by a drop in heart rate (Braby and Somero, 2006b).

In this study, we have chosen to focus on the proteome to characterize the molecular mechanisms that set environmental tolerance limits, as changes in protein abundance represent modifications of the molecular phenotype of the cell and therefore functional changes (Feder and Walser, 2005). Mass-spectrometry-enabled proteomic analyses were first made possible with the completion of genome sequencing projects for model organisms (Aebersold and Mann, 2003; Mann et al., 2001). Through advances in mass spectrometry and the generation of expressed sequence tag (EST) libraries, proteomic studies on non-model organisms have constantly improved, leading to the generation of a number of new hypotheses about the stress responses of organisms to environmental change (Serafini et al., 2011; Tomanek et al., 2011; Tomanek, 2011; Tomanek, 2012).

By comparing proteomic responses to acute and chronic temperature stress in two closely related species of Mytilus that vary in distribution and invasiveness, we have generated several new hypotheses about how differently adapted congeners vary in their cellular responses to thermal stress and which cellular processes are involved in setting tolerance limits (Fields et al., 2012; Tomanek and Zuzow, 2010); simultaneously, our collaborators have focused on the transcriptomic responses of these congeners to acute heat and hyposaline stress (Lockwood et al., 2010; Lockwood and Somero, 2011). Here we exposed both blue mussel congeners to short exposures (4 h) of hyposaline stress (24.5 and 29.8 psu and a control of 35 psu), followed by a 0 and 24 h recovery at 35 psu, to mimic conditions typical for bays and coastal areas experiencing heavy freshwater input with a quick return to full salinity due to incoming tides and mixing with full-strength seawater. Our results in the current study indicate that the native M. trossulus is able to respond to a greater range of salinity variations than the invasive Mediterranean M. galloprovincialis. This increased plasticity with respect to salinity tolerance may better equip the native M. trossulus to compete with the invader in regions with warmer water and more frequent hyposaline stress despite the invaders increased heat tolerance. Our proteomic analysis implicates protein homeostasis, vesicle transport and cytoskeletal rearrangements as well as modifications in energy metabolism and oxidative stress response as cellular processes setting interspecific differences in salinity tolerance.

Mytilus trossulus Gould 1850 and M. galloprovincialis Lamarck 1819 were collected subtidally from Newport, OR, USA (44°38′25″N, 124°03′10″W), and Santa Barbara, CA, USA (34°24′15″N, 119°41′30″W), respectively. In a separate study, PCR was used to confirm that each site was occupied by only a single species (i.e. there were no hybrids present) (Lockwood et al., 2010). The experimental conditions were chosen to simulate temporary hyposaline stress conditions as they occur in estuaries and bays during heavy winter rains in California near the hybrid zone. However, these conditions are often quickly reversed due to incoming tides and dilution of freshwater.

Animals were kept for 4 weeks under constant immersion at 13°C in recirculating seawater tanks with a salinity of 35 psu and fed a phytoplankton diet (Phytofeast, Reed Mariculture, Campbell, CA, USA) every day. We employed two experimental treatments, 24.5 and 29.8 psu, and one control treatment of 35.0 psu. All treatments were kept at 13°C for the duration of the experiments. Animals were exposed for 4 h (or 0 h recovery), at which point we collected the first set of gill tissues (N=4–6 for all treatments). Another set was collected after a 24 h recovery period at 35.0 psu (N=6 for each treatment). The actual osmolalities measured with an osmometer (Advanced Instruments, Norwood, MA, USA) were 750, 858 and 979 mOsm kg⁻¹ for the 24.5, 29.8 and 35.0 psu treatments, respectively. The first time point was chosen because it coincides with the time of collection of the samples used for the transcriptomic analysis (Lockwood and Somero, 2011), the second one because it allowed the organism to respond to the stress by translating proteins in high enough abundances and assessed the proteomic response to a hyperosmotic stress (relative to 24.5 and 29.8 psu) upon return to control conditions (35.0 psu). One possible behavioral response of Mytilus to hyposaline stress is shell closure to avoid direct contact with the medium (Braby and Somero, 2006b), which would be difficult to control. To avoid this confounding variable, we placed a small cork (5 mm diameter) between the shells to characterize the cellular response of gill tissue to the three salinity treatments. Mussels were immediately dissected on chilled aluminum foil and tissues were kept frozen at −80°C until processing.

Sample preparation followed procedures outlined previously (Tomanek and Zuzow, 2010). Briefly, gill tissue was lysed in homogenization buffer [7 mol l⁻¹ urea, 2 mol l⁻¹ thiourea, 1% amidosulfobetaine-14, 40 mmol l⁻¹ Tris-base, 0.5% immobilized pH 4–7 gradient (IPG) buffer (GE Healthcare, Piscataway, NJ, USA) and 40 mmol l⁻¹ dithiothreitol] at a ratio of 1:4. After centrifugation at 20°C for 30 min at 16,100 g, the proteins were precipitated by adding four volumes of ice-cold 10% trichloroacetic acid in acetone and incubating the solution at −20°C overnight. The precipitate was centrifuged at 4°C for 15 min at 18,000 g, the supernatant was discarded, and the protein pellet was washed with ice-cold acetone and centrifuged again at 4°C. After air-drying, the pellet was re-suspended in rehydration buffer (7 mol l⁻¹ urea, 2 mol l⁻¹ thiourea, 2% cholamidopropyl-dimethylammonio-propanesulfonic acid, 2% nonyl phenoxylpolyethoxylethanol-40, 0.002% Bromophenol Blue, 0.5% IPG buffer and 100 mmol l⁻¹ dithioerythritol). The protein concentration was determined with the 2D Quant kit (GE Healthcare), according to the manufacturer's instructions.

Prior to isoelectric focusing, IPG strips (pH 4–7, 11 cm; BioRad Laboratories, Hercules, CA, USA) were passively rehydrated with 200 μl of 2.5 μg μl⁻¹ protein in rehydration buffer in wells for 13 h. Isoelectric focusing was conducted using the following protocol: 250 V for 15 min, gradient voltage increase to 8000 V for 1 h, 8000 V for 3 h 45 min, and reduced to 500 V (Ettan IPGphor3, GE Healthcare).

To prepare for second-dimension SDS-PAGE electrophoresis, strips were incubated in equilibration buffer (375 mmol l⁻¹ Tris-base, 6 mol l⁻¹ urea, 30% glycerol, 2% SDS and 0.002% Bromophenol Blue) for two 15 min intervals, first with 65 mmol l⁻¹ dithiothreitol and then with 135 mmol l⁻¹ iodoacetamide. IPG strips then were placed on top of 11.8% polyacrylamide gels, which were run (Criterion Dodeca, BioRad Laboratories) at 200 V for 55 min at 10°C. Gels were subsequently stained with colloidal Coomassie Blue (G-250) and destained with Milli-Q water for 48 h. The resulting gels were scanned with an Epson 1280 transparency scanner (Epson, Long Beach, CA, USA).

Digitized images of two-dimensional (2-D) gels were analyzed using Delta2D (version 3.6, Decodon, Greifswald, Germany) (Berth et al., 2007). Spot boundaries were detected on a fused composite 2-D gel image and transferred back to the original gel images. After background subtraction, the relative amount of protein in each spot (i.e. spot volume) was quantified by normalizing against total spot volume of all proteins in the image.

To determine which spots changed significantly in response to salinity (24.5, 29.8 and 35.0 psu) and recovery time (0 and 24 h), we used a two-way ANOVA (P<0.02) with salinity and recovery time as the main effects and the effect of time on the response to salinity as the interaction effect. We generated a null distribution for the two-way ANOVA (1000 permutations) to account for unequal variance and non-normal distributions of the response variables. In the Results and Discussion, we include all identified proteins of certain functional categories and indicate whether they showed significance for one or both of the main effects and the interaction. The complete data set, separated by main and interaction effects, is available in supplementary material Figs S1–S4 and Tables S1, S2. Because there is only limited overlap between the proteome maps of the two congeners, as well as uncertainty whether overlapping proteins were orthologous or paralogous homologs, a two-way ANOVA comparing species was not possible. Following the two-way ANOVA, post hoc testing to compare treatments was conducted using Tukey's analysis (P<0.05) in Minitab (version 15, Minitab, State College, PA, USA), to support conclusions about differences in protein abundances (single-protein graphs are not shown).

Proteins that changed abundance in response to temperature acclimation were excised from gels and prepared for analysis by mass spectrometry (MS) following previously published protocols (Fields et al., 2012; Tomanek and Zuzow, 2010).

We obtained peptide mass fingerprints (PMFs) using a matrix-assisted laser desorption ionization tandem time-of-flight mass spectrometer (Ultraflex II, Bruker Daltonics, Billerica, MA, USA). We selected a minimum of six and a maximum of 20 peptides for tandem MS in order to obtain information about their b- and y-ions.

Analysis of peptide spectra followed previously published procedures (Fields et al., 2012; Tomanek and Zuzow, 2010). We used flexAnalysis (version 3.0, Bruker Daltonics) to detect peptide peaks (with a signal-to-noise ratio of 6 for MS and 1.5 for MS/MS). Porcine trypsin (Promega, Madison, WI, USA) was used for internal mass calibration.

To identify proteins we used Mascot (version 2.2, Matrix Science, Boston, MA, USA) and combined PMFs and tandem mass spectra in a search against two databases. One database is an EST library that represents 12,961 and 1688 different gene sequences for M. californianus and M. galloprovincialis, respectively (Lockwood et al. 2010). The other was NCBI, with 77,410 total nucleotide sequences with Mytilus as the taxonomic restriction, and 5,266,919 sequences under Metazoa. Oxidation of methionine and carbamidomethylation of cysteine were our only variable modifications. Our search allowed one missed cleavage during trypsin digestion. For tandem MS we set the precursor-ion mass tolerance to 0.6 Da, the default value in Mascot. The molecular weight search (MOWSE) score that indicated a significant hit was dependent on the database: scores higher than 40 and 51 were significant (P<0.05) for a search in the Mytilus EST and NCBI database, respectively. However, we only accepted positive identifications that included two matched peptides regardless of the MOWSE score.

To associate proteins with similar changes in abundance across samples, we employed hierarchical clustering with average linking (Delta2D), using a Pearson correlation metric. To further assess the importance of specific proteins in differentiating the proteomes of mussels exposed to different salinities, we employed principal component analysis (PCA; Delta2D) based on proteins whose abundances changed significantly during and after exposure to hyposaline stress (two-way ANOVA, P<0.02). Component loadings, which quantify the contribution of each protein in the separation of samples along a given component, are reported in supplementary material Figs S2–S4 only if greater than ±1.0.

Proteins from gill tissue of mussels exposed to three salinity levels (24.5, 29.8 and 35 psu) for 4 h and mussels that were given a chance to recover from these salinities at control conditions (35.0 psu) for 24 h were separated by 2-D gel electrophoresis and yielded 336 and 310 distinct protein spots in M. trossulus and M. galloprovincialis, respectively (Fig. 1). Of the total protein spots, 39% in M. trossulus and 29% in M. galloprovincialis changed in response to hyposaline conditions.

In M. trossulus, principal component 1 (PC1) explains 27.4% of the variation of protein abundance data and separates out the mussels exposed to 24.5 psu for 4 h followed by exposure to 35.0 psu for 24 h during recovery (24.5 psu +24 h) (Fig. 2A). Along the y-axis, PC2 explains 12.3% of the variation, and it separates the control from the 29.8 psu treatment (Fig. 2A). These two PCs show that the greatest variation in the data is found in the response of M. trossulus 24 h after an acute exposure to 24.5 psu (PC1), followed by the variation between the control and 29.8 psu treatments (PC2; Fig. 2A). These patterns suggest that the broadest proteomic adjustments of gill tissue occur during recovery from 24.5 psu.

In M. galloprovincialis, the contributions of PC1 and PC2 (26.6 and 14.0%, respectively) to explaining the variation in protein abundance in response to hyposaline treatment and recovery conditions are similar to those of M. trossulus. But in contrast to M. trossulus, it is the 29.8 psu +24 h treatment that is separated the furthest from the other treatments along PC1 (Fig. 2B). PC2 separates 29.8 psu 0 h from the remaining treatments. Thus, PC1 and PC2 indicate that the 29.8 psu hyposaline treatment causes the largest variation in protein abundance in M. galloprovincialis, more so after +24 h than 0 h recovery.

Despite explaining similar levels of variation in protein abundance in both species, the PCAs reveal differences in how the two species vary in their response to hyposaline stress. In summary, M. trossulus responds strongest +24 h into the recovery from a 24.5 psu exposure, whereas M. galloprovincialis responds strongest +24 h into the recover from a 29.8 psu exposure while showing limited changes in protein abundance to 24.5 psu. PC2 clusters three hyposaline conditions close to each other, with the exception of 24.5 psu +24 h, and placed them in opposition to the control treatments in M. trossulus, suggesting that these conditions require a similar proteomic response and thus do not differ among each other enough to represent greatly differing stress levels. In M. galloprovincialis, it is only the 29.8 psu 0 h exposure that is placed in this position along PC2.

To assess the acute response as well as the recovery from hyposaline stress, we collected samples at two time points, 0 h and +24 h into recovery (at 35.0 psu). This scenario mimics the effect of a heavy rain event diluting full-strength into more brackish seawater, just to return to full-strength seawater after the cessation of the rain event, or more accurately, the down-flow of a freshwater surface layer through an estuary or near coastal waters. Of the 336 proteins detected in M. trossulus, 27% (91 spots) changed during recovery; in M. galloprovincialis, 29% (89 of the 310 spots) changed.

In M. trossulus, PC1 explains 17.1% of the variation and separates the 24.5 psu + 24 h from the 0 h time point, with the 29.8 and 35.0 psu +24 h treatments in between (Fig. 3A). PC2 explains 9.9% of the variation in M. trossulus and mainly separates the 29.8 and 35.0 psu +24 h treatments (negative on y-axis) from all other treatments.

In M. galloprovincialis, the first PC explained 28% of the variation in protein abundance (Fig. 3B), approximately 11% more than in M. trossulus. Overall, PC1 separates all +24 h from all 0 h treatments. The group separated the furthest along PC1 is 29.8 psu +24 h. The other +24 h recovery treatments, 24.5 and 35.0 psu, are separated from 29.8 psu and overlap. PC2 in M. galloprovincialis explains 12.9% of the variation in protein abundance and separates the 24.5 and 35.0 psu 0 h and 29.8 psu +24 h treatments (positive range of PC2) from all others (Fig. 3B).

The PCAs show a clear separation with decreasing salinities during recovery in M. trossulus along PC1 (Fig. 3A). There is little separation among the 0 h groups. This suggests that most of the proteomic changes occur during recovery and are greater with decreasing salinity.

Mytilus galloprovincialis shows a similar pattern of separation along PC1, but with 29.8 instead of 24.5 psu +24 h being the treatment with the greatest separation and 24.5 and 35.0 psu +24 h overlapping (Fig. 3B).

Both species show species-specific changes in the abundance of chaperones that are localized to the mitochondria, prohibitin (Liu et al., 2009) and the endoplasmic reticulum (ER), e.g. 78 and 94 kDa glucose regulated protein (GRP78 or BiP and GRP94), protein disulfide isomerase (PDI) and translocon-associated protein β (part of the Sec61 channel to translocate proteins during translation into the lumen of the ER) (Araki and Nagata, 2012). GRP94 is a heat shock protein (HSP) 90 homolog that facilitates folding of secreted and membrane proteins and holds misfolded proteins until they can be transported out of the ER for further degradation (Araki and Nagata, 2012; Eletto et al., 2010). It also is a major calcium binding protein in the lumen of the ER and its upregulation is considered an indicator of ER stress, mainly because of its activation of insulin-like growth factors, which facilitate recovery from ER stress while blocking apoptosis (Eletto et al., 2010). GRP78 or BiP may precede GRP94 as a folding catalyst (Melnick et al., 1994).

Although M. trossulus showed the highest abundances of two GRP94 isoforms and cystatin-B at 24.5 psu 0 h (Fig. 4A; cluster T_CHA), M. galloprovincialis increased abundances of GRP94, two GRP78 isoforms, heat shock cognate (HSC) 70 and PDI at 29.8 psu 0 h (Fig. 4B; cluster G_CHB). These interspecific differences parallel our results from the PCAs (Fig. 2) and suggest that M. trossulus is able to tolerate greater acute hyposaline stress than M. galloprovincialis before disruption of proteostasis in the ER.

In addition, abundances of T-complex protein 1 (TCP-1), a tubulin- and actin-folding chaperone, decreased during hyposaline treatments (0 h) in M. galloprovincialis only, suggesting that proper folding of cytoskeletal elements, such as building blocks for cilia, was disrupted (Fig. 4B) (Sternlicht et al., 1993). One small HSP whose main function is to stabilize cystoskeletal elements (Haslbeck et al., 2005) showed overall higher levels at all salinities at 0 h than after 24 h of recovery in M. galloprovincialis (spot 41 was also identified as a small HSP but has a much higher than expected molecular mass and thus may not be a small HSP). Together, these data suggest that proteostasis, especially in the ER, of the cytoskeleton and possibly cilia, is important in setting species-specific limits to hyposaline stress in Mytilus gill tissue. Protein folding in the ER is important for secreted proteins, especially as part of the mucus that is transported across the ventral grove of the gill to capture food particles that will be transported towards the mouth through ciliary movements.

The ER maintains an oxidizing environment that facilitates the formation of disulfide bonds (Araki and Nagata, 2012; Csala et al., 2010). As a consequence, protein folding in the ER is closely linked and sensitive to changes in the redox environment. For example, abundance changes in GRP94 and PDI, members of a subfamily of the thioredoxin-like proteins (Funato and Miki, 2007), represent key indicators for the disruption of proteostasis in the ER (Eletto et al., 2010; Feige and Hendershot, 2011). Importantly, reactive oxygen species (ROS) cannot only interrupt disulfide bonds but are actually generated by the oxidation of sulfhydryl groups in the ER, specifically hydrogen peroxide, and may make up as much as a quarter of all ROS produced in the cell (Araki and Nagata, 2012; Csala et al., 2010; Malhotra and Kaufman, 2007). Furthermore, a cluster of three proteins in M. galloprovincialis may play important roles in protein folding or ROS scavenging in the ER: thioredoxin-like protein [a protein disulfide reductase (Holmgren and Lu, 2010)], nucleoredoxin [a putative thioredoxin (Funato and Miki, 2007)] and superoxide dismutase (Fig. 5B; cluster G_EC). Their abundances decreased during recovery from 24.5 psu +24 h, possibly indicating a downregulation of protein folding activity and protein synthesis in the ER in response to hyposaline stress, which would explain why the proteomic response at 24.5 psu in M. galloprovincialis was closer to 35.0 psu than 29.8 psu (Fig. 2B).

Two proteins that are part of cluster G_EC (Fig. 5B), NADH dehydrogenase [complex I of the electron transport chain (ETC)] and superoxide dismutase (SOD), are shared between the congeners. Although abundances of NADH dehydrogenase overall were lower during recovery, they were comparatively higher at 29.8 psu +24 h in comparison to 24.5 psu and the control +24 h treatments in both congeners (Fig. 5). However, SOD showed decreasing abundances during recovery from 24.5 psu +24 h only in M. galloprovincialis, in contrast to M. trossulus, which decreased SOD at 24.5 and 29.8 psu +24 h. Isoforms of NADP-dependent isocitrate dehydrogenase (NADP-ICDH) are part of this cluster (Fig. 5A, cluster T_EA) and showed reduced abundances at 24.5 and 29.8 psu +24 h in M. trossulus. We have hypothesized that all three proteins may play a role in regulating oxidative stress, through ROS production (NADH dehydrogenase), ROS scavenging (SOD) or maintenance of high levels of reduced glutathione for ROS scavenging in Mytilus in the mitochondria during acute heat stress and acclimation to cold (NADP-ICDH) (Fields et al., 2012; Tomanek and Zuzow, 2010). Of these three, at least NADP-ICDH has been shown to reside in the ER (Margittai and Bánhegyi, 2008) and could contribute to ROS scavenging in the ER. SOD could scavenge the hydrogen peroxide normally produced during protein folding in the ER.

The picture that emerges is one of protein unfolding in the ER during the acute phase of hyposaline stress, as indicated by the upregulation of the molecular chaperones GRP78 and GRP94, with species-specific abundance patterns (e.g. 29.8 and 24.5 psu in M. galloprovincialis and M. trossulus, respectively) and proteins (e.g. PDI in M. galloprovincialis), followed by a reduction in protein synthesis and folding during recovery, as indicated by reduced abundances of a subset of the same proteins [e.g. GRP78 and GRP94 at 29.8 and 24.5 psu in M. galloprovincialis and M. trossulus (GRP94 spot 36 only), respectively]. The proposed reduction in protein synthesis and protein folding in the ER would cause a reduction in the production of ROS, specifically H₂O₂, which may be indicated by the lower abundances of proteins involved in ER redox regulation in M. galloprovincialis (thioredoxin-like and nucleoredoxin at 24.5 psu +24 h). The lower abundances of additional oxidative stress proteins (SOD and NADP-ICDH), possibly located in the ER or the nearby cytosol, during recovery also supports an inference of lower levels of oxidative stress. Further support for the notion of reduced protein synthesis may come from two proteins, HSC71 and translocon-associated protein, from cluster T_CHB in M. trossulus (Fig. 4A), both of which showed increasingly lower abundances with lower salinities during recovery (between 35 and 24.5 psu for HSC71 and between 35 or 29.8 and 24.5 psu for translocon-associated protein) and are indicators of chaperone activity of newly synthesized proteins that are processed through the ER (Araki and Nagata, 2012). Although the comparison between the congeners is suggestive, a more comprehensive characterization is necessary before we can discern that differences in regulating the link between ER-localized protein maturation and ROS production contribute to setting tolerance limits to hyposaline stress.

Finally, proteases break down irreversibly denatured proteins and thereby remove them from a pool of possibly toxic aggregates that could interact with other functioning proteins (Wong and Cuervo, 2012). In contrast to acute heat stress, where we identified a number of proteasome isoforms (Tomanek and Zuzow, 2010), in the present study we identified only one proteasome α-type subunit in M. trossulus that showed higher abundance at 24.5 psu +24 h (Fig. 4A). Cystatin-B is a protease inhibitor, especially of cysteine proteases, which binds irreversibly to proteases and thereby protects cells from their activity (Chapman et al., 1997). We identified three isoforms of cystatin-B, with higher abundances during acute stress (spot 16), control conditions (spot 14) and recovery from extreme hyposaline stress (spot 15), with only minor shifts in molecular mass, thus possibly suggesting a role for PTMs in regulating their activity. Interestingly, cystatin-B together with fatty acid binding protein (FABP) (see below; Fig. 6) have both been suggested to be urinary biomarkers for acute kidney injury (Vaidya et al., 2008).

Because the production of ROS is closely linked to the ETC and therefore to energy metabolism, we cover both functional categories together (Murphy, 2009). Proteins involved in energy metabolism and those indicating oxidative stress showed more pronounced changes during recovery than during acute hyposaline stress in M. trossulus (Fig. 5A). The hierarchical clustering showed two main patterns: one cluster with abundances decreasing at 24.5 and/or 29.8 psu during recovery (T_EA), and another with increasing abundances mainly at 24.5 psu +24 h (T_EB). Proteins of cluster T_EA [with the exception of ATP synthase (spot 6) and NADH dehydrogenase] showed decreasing abundances in response to hyposaline stress during recovery. Proteins of this cluster represent the pyruvate dehydrogenase (PDH) reaction [dihydrolipoyl dehydrogenase (DLDH) is part of the PDH complex] as well as the Krebs cycle [mitochondrial malate dehydrogenase and NADP-ICDH] and ATP production (ATP synthase). With the exception of the latter enzyme, they were all hypothesized to respond to increased ROS production by decreasing ROS-generating NADH-producing pathways while increasing ROS-scavenging NADPH-producing pathways, in the case of NADH-ICDH, during acute heat stress in M. trossulus (Tomanek and Zuzow, 2010). A similar response may be seen here during recovery from hyposaline stress, with the exception that abundances of NADP-ICDH did not increase. A possible reason for this may be that we were not able to distinguish between the cytosolic (and ER) and the mitochondrial isoforms of NADP-ICDH (Margittai and Bánhegyi, 2008). However, given that three typical oxidative stress proteins, DyP-type peroxidase (a catalase) (Sugano, 2009), SOD and the mitochondrial isoform of aldehyde dehydrogenase (ALDH) (Ellis, 2007) reduced abundances in parallel to the decreasing abundances of NADH-producing enzymes, this suggests that the changes in proteins involved in energy metabolism may be linked to reduced ROS production.

The complementary cluster (T_EB) mainly showed increasing abundances at 24.5 psu +24 h (Fig. 5A). The two ALDH isoforms are involved in the detoxification of different species of aldehydes, which are produced in part through other ROS interacting with the double bonds of polyunsaturated fatty acids and thus lipid peroxidation (Ellis, 2007). The electron transfer flavoprotein-α transfers electrons that are made available through the β-oxidation of fatty acids via FADH₂ to the ETC (Salway, 2004). Finally, propionyl CoA carboxylase plays a role in the metabolic pathways of valine, methionine and threonine oxidation to succinyl CoA (Salway, 2004). These changes also indicate increasing levels of one specific type of oxidative stress, possibly limited to a specific group of macromolecules, e.g. lipids, as well as possible alternative strategies to regulate energy metabolism to reduce ROS production.

Mytilus galloprovincialis gill tissue showed three clusters: one with decreasing abundances at 24.5 psu +24 h (Fig. 5B, G_EC), similar to the one discussed for M. trossulus (Fig. 5A, T_EA), one with lower abundances at 29.8 psu +24 h (G_EA) and one with higher abundances (G_EB) at 29.8 psu +24 h. Proteins involved in producing (PDH) and oxidizing NADH (NADH dehydrogenase), as well as SOD as a scavenger of hydrogen peroxide, and nucleoredoxin, a thioredoxin and therefore a disulfide reductase (Funato and Miki, 2007), all showed lower abundances at 24.5 psu +24 h. In a direct comparison of the same proteins (SOD and PDH or DLDH), M. trossulus showed lower abundances at 24.5 and 29.8 psu +24 h. These results are suggestive of an important role for a reduction in energy metabolism, e.g. metabolic depression, in setting limits to hyposaline conditions, possibly through the reduced production of ROS.

Clusters G_EA and G_EB are complementary and indicate that while ATP synthase abundance is up, the abundances of oxidative stress proteins, such as DyP-type peroxidase and peroxiredoxin 5, are down at 29.8 psu +24 h (Fig. 5B). The cytosolic paralog of malate dehydrogenase showed two isoforms in both clusters, suggesting a possible PTM, e.g. acetylation, regulating its activity (Zhao et al., 2010). This pattern suggests that there may be a transitory increase in energy demand during recovery from 29.8 psu in M. galloprovincialis.

In summary, during recovery from hyposaline stress, metabolic pathways involving NADH production and oxidation are downregulated to a greater extent in M. trossulus, including exposure to both 24.5 and 29.8 psu +24 h, than in M. galloprovincialis, which showed decreasing abundances only at 24.5 psu +24 h. These changes are paralleled by decreasing abundances of oxidative stress proteins, with some proteins likely localized to the ER, where we hypothesize that they showed decreasing abundances due to a decrease in protein synthesis and folding of proteins with disulfide bridges, which in turn may lower the production of ROS. This link between reduced protein synthesis and folding and lower levels of ROS production could be an underappreciated reason for the translational arrest during stress (Holcik and Sonenberg, 2005). Two additional themes distinguished the proteomic response of the congeners. First, M. trossulus showed changes in proteins indicating an upregulation of metabolic pathways (β-oxidation and metabolism of branched amino acids) at 24.5 psu +24 h that were not seen in M. galloprovincialis and could indicate alternative metabolic pathways used by M. trossulus during hyposaline stress. Second, M. galloprovincialis showed increasing abundances of ATP synthase but lower abundances of oxidative stress proteins at 29.8 psu +24 h, possibly indicating a transient increase in energy demand that M. trossulus did not show.

Proteins constituting the cytoskeleton or elements of cilia, actin binding and regulatory proteins as well as small GTPases involved in vesicle formation and transport showed three major clusters in M. trossulus: one in which five actins, one α-tubulin and gelsolin, an actin severing protein (Silacci et al., 2004), showed higher abundances at mild (29.8 psu +24 h) but, in the case of some proteins, lower abundances at extreme (24.5 psu +24 h) hyposaline stress during recovery (Fig. 6A; cluster T_CC). A complementary cluster showed higher abundances at extreme hyposaline stress and included three actins, a β-tubulin, F-actin capping protein β, G-protein β and Rab1-GDP dissociation inhibitor (Rab1-GDI; cluster T_CA). Both clusters contain an isoform of the Na⁺/H⁺ exchange regulatory factor (NHE-RF). A third cluster is characterized by lower abundances at one or both hyposaline stress conditions during recovery (+24 h) and contains an actin, α-tubulin, actophorin (a cofilin or actin depolymerization factor) and Ras-like GTPase Sar1 (cluster T_CB). Clusters T_CC and T_CB both contain an isoform of FABP.

The distinct changes in clusters that mainly contain actin isoforms during recovery with different levels of hyposaline stress (T_CA and T_CC) may be explained in part by actin-binding and regulatory proteins that are also part of these clusters. For example, at 24.5 psu +24 h, abundances of actophorin and gelsolin are lower, while abundance of the F-actin capping protein is higher (Fig. 6A). Lower abundances of the former proteins indicate that ‘treadmilling’ of actin or the growth of actin filaments, a process that can expand the cell membrane and therefore cell volume, is inhibited upon return to control conditions following extreme hyposaline stress (Le Clainche and Carlier, 2008). This hypothesis is further supported by the simultaneously higher abundances of F-actin capping protein, which would prevent actin filaments from growing.

We also identified two small GTPases – Ras-like GTPase Sar1, which recruits membrane coat proteins that facilitate vesicle formation, and Rab1-GDI, a protein that inhibits Rab1 – which regulate vesicle transport from the ER to the Golgi apparatus (Di Ciano-Oliveira et al., 2006; Marks et al., 2009). Thus, the simultaneously higher abundance of Rab1-GDI and lower abundance of Ras-like GTPase Sar1 during recovery from extreme hyposaline stress may be hypothesized to indicate a downregulation of vesicle formation and transport from the ER, possibly reversing the activation of these processes during acute hyposaline stress (0 h).

Two isoforms of NHE-RF also changed in opposite clusters (T_CA and T_CC). NHE-RF can be phosphorylated by protein kinase A and affects the signaling of G-protein coupled receptors in addition to transporters (e.g. Na⁺/H⁺ exchanger), ion exchangers and signaling proteins (Ardura and Friedman, 2011). Some NHE-RFs have a C-terminal binding domain that connects them to the cytoskeleton, suggesting a role in sensing cytoskeletal modifications and, by extension, cell volume (Thelin et al., 2005). Given the difference in molecular mass between the isoforms (13 kDa), they may present different orthologs rather than PTMs (Ardura and Friedman, 2011).

Finally, the role of the two FABP isoforms is unclear. Their abundance changes are complementary, possibly because of PTMs (Fig. 6A). They may be involved in the synthesis of lipids, including phospholipids, in the ER to modify membranes that may be transported to the outer cell membrane (Storch and Thumser, 2000).

To understand the changes associated with the cytoskeleton and vesicle transport, it is important to recall that the PCAs for M. galloprovincialis (Figs 2, 3) showed limited proteomic changes for the extreme hyposaline stress conditions. During acute stress (0 h), proteins represented in cluster G_CB showed higher abundances at 29.8 psu only (Fig. 6B). The majority of those are five isoforms of α-tubulin, one isoform of β-tubulin, three actins and Rab1-GDI, which would indicate that vesicle formation and transport are inhibited during the early response to mild hyposaline stress. At least during 0 h, cluster G_CD included proteins with lower abundances at 29.8 psu, such as: radial spoke head 9 (RSH9), a cilia protein; Rho-GDI, an inhibitor of the small GTPase Rho; a β-tubulin; profilin, an actin-binding protein that increases the rate and affects the direction of actin treadmilling as well as prevents G-actin aggregation, depending on its PTMs (Le Clainche and Carlier, 2008); and two isoforms of myosin light chain 1, which may be connecting actin and myosin near the periphery of the cell (Estévez-Calvar et al., 2011). The cluster is in some way complementary to G_CB, at least during the acute phase of the stress.

During recovery (+24 h), proteins of cluster G_CA showed higher abundances at mild hyposaline stress (Fig. 6B). They include α- and β-tubulins, intermediate filament and two isoforms of RSH. This cluster is similar to G_CB (higher abundances at 0 h) in that it contains several tubulin isoforms. Cluster G_CC showed the opposite patterns during recovery (+24 h) and contains three actins, β-tubulin and F-actin capping protein.

Although species-specific patterns of protein abundance exist, namely the greater number of tubulin isoforms changing abundance in M. galloprovincialis but not M. trossulus (Fig. 6) and specific proteins that were only identified for one of the congeners, e.g. Ras-like GTPase Sar1 and FABP in M. trossulus, these differences and the proteins the congeners have in common point to a related cellular response to hyposaline stress. This response includes vesicle formation and transport in response to osmotic cell swelling (van der Wijk et al., 2003), represented in part by the small GTPases known to affect this process (Di Ciano-Oliveira et al., 2006; Marks et al., 2009). In addition, vesicle transport, with a close connection to modifications to cilia architecture, occurs with the help of tubulin, and depends on radial spokes (Silverman and Leroux, 2009). The other set of proteins associated with cell-volume regulation includes the actin-based cytoskeleton, specifically those proteins that are involved in actin ‘treadmilling’ (Le Clainche and Carlier, 2008), which seems to be regulated during recovery (Fig. 6). The species-specific patterns point to a role for tubulin, and possibly its PTMs (specifically acetylation), as an important process in affecting vesicle transport and cytoskeletal rearrangements (Perdiz et al., 2011), and thereby reduced tolerance towards hyposaline conditions in M. galloprovincialis. This hypothesis is further supported by the observation of decreasing abundances of Rho-GDI, an inhibitor of the small GTPase Rho, which has been shown to control this process (Destaing et al., 2005), in M. galloprovincialis during mild hyposaline stress. Rho also affects several downstream protein kinases, which in turn either indirectly, through additional kinases, or directly affect myosin light chains and thereby cell volume, the cellular stress response, several actin-binding proteins and the formation of actin stress fibers (Di Ciano-Oliveira et al., 2006; Marks et al., 2009). These changes, in addition to those directly linked to vesicle formation and transport, suggest that small GTPase-mediated processes contribute to setting species-specific limits to hyposaline conditions.

The proteomic response of both Mytilus congeners to hyposaline stress showed common themes: ER molecular chaperones indicate protein unfolding during the acute phase; vesicle transport and cytoskeletal modifications suggest adjustments in cell volume, especially during recovery; and proteins involved in energy metabolism and ROS scavenging indicate that a reduction in energy demand may be accompanied by reduced ROS production, also during recovery. However, the differences in protein abundances suggest that M. trossulus can respond to a greater hyposaline challenge (24.5 psu) than M. galloprovincialis (29.8 psu), specifically during recovery. It is possible that a reduction of protein folding in the ER during recovery may be linked to decreased oxidative stress in the ER, thereby lowering ROS production and, as a possible consequence, protein denaturation (Dalle-Donne et al., 2003), more so in M. galloprovincialis than in M. trossulus. Both vesicle transport and cytoskeletal modifications play a role in the response to hyposaline stress. While in M. trossulus, the abundances of a number of actin-binding regulatory proteins changed, a number of tubulin isoforms changed in M. galloprovincialis. Although the former may be linked to adjustments in cell volume, the latter may be linked to the transport of membrane vesicles, possibly to first increase cell volume during acute hyposaline stress and then to retrieve membranes during recovery. Changes in proteins involved in energy metabolism indicate an overall reduction in energy metabolism upon return to control conditions in both congeners, with an indication of a transient increase in energy metabolism at mild hyposaline stress (29.8 psu) during recovery in M. galloprovincialis, suggesting species-specific differences in time course and scope of adjustment in energy metabolism. In general, abundances in oxidative stress proteins parallel changes of proteins involved in energy metabolism.

Abundance changes of ER chaperones in response to osmotic stress have also been observed in proteomic analyses of mouse embryonic stem cells and kidney cells (Dihazi et al., 2005; Mao et al., 2008). Proteins involved in small GTPase and cytoskeletal pathways were enriched in osmoregulatory tissues of sharks (Lee et al., 2006). Several of the proteins representing energy metabolism in Mytilus were also found in the rectal glands of sharks in response to a feeding-associated salt load (Dowd et al., 2008), but shark gill tissue showed a number of proteasome isoforms in response to salinity change (Dowd et al., 2010), a response that was almost absent in Mytilus. Our results indicate that these cellular processes play an important role in setting tolerance limits towards hyposaline stress. Furthermore, the number of actin-binding regulatory proteins and tubulin isoforms potentially associated with vesicle transport provide novel insights into the cellular processes contributing to salinity tolerance limits, especially in gill tissue, which excretes proteins as part of the mucus needed to trap food. A comparison of the proteomic responses of Mytilus gill tissue to acute heat stress and temperature acclimation with the current data set shows some stressor-specific cellular processes, e.g. protein degradation during acute heat stress, as well as responses that are common to all of the stressors, e.g. a trade-off between energy metabolism and oxidative stress (Tomanek, 2012). Together, these studies emphasize the importance of oxidative stress, and the comparisons between Mytilus congeners suggest that ROS-induced physiological tolerance limits play an important role in setting biogeographic distribution limits.

Finally, unlike our proteomic analysis, the transcriptomic analysis of gill tissue of Mytilus specimens from the same experiment (but limited to the 35 and 29.8 psu +0 h treatments) showed very few changes between the congeners (Lockwood and Somero, 2011). In addition, there is almost no overlap between the transcript and our protein abundance changes, suggesting that interspecific differences at the level of the proteome are crucial to setting tolerance limits to hyposaline stress. Some of the proteomic changes observed here are likely based on PTMs, e.g. FABP in M. trossulus (Fig. 6A), a conclusion that is supported by changes in protein kinase activities during hyposaline stress in Mytilus (Evans and Somero, 2010).

Thus, the comparison of the proteomic responses of gill tissue of both congeners to hyposaline stress conditions shows that, at the level of the molecular phenotype, the warm-adapted M. galloprovincialis may be limited in its expansion north by an increase in precipitation events and increased freshwater input near coastal waters. Moreover, it is significant to note that this study illustrates possible molecular level mechanisms to predict the results of closely related species competition in response to climate change.

We thank Daniel D. Magee, Jeremy K. LaBarge and Brent L. Lockwood for their assistance in conducting the original experiment. The experimental design was created in collaboration with Brent L. Lockwood and George N. Somero of Stanford University (Lockwood and Somero, 2011). Peter Field, Jennifer Oquendo and Shelley Blackwell provided helpful editorial suggestions.

 
• ALDH

  aldehyde dehydrogenase

 
• BiP

  binding immunoglobulin protein

 
• CCT

  chaperonin containing TCP-1

 
• DLDH

  dihydrolipoyl dehydrogenase

 
• ER

  endoplasmic reticulum

 
• EST

  expressed sequence tag

 
• ETC

  electron transport chain

 
• FABP

  fatty-acid binding protein

 
• F(G)-actin

  filamentous (globular)-actin

 
• FADH₂

  flavin adenine dinucleotide dihydrogen

 
• GRP

  glucose-regulated protein

 
• G(D)TP

  guanosine 5’-(di-)triphosphate

 
• G_CA

  M. galloprovincialis; cytoskeleton-associated proteins; cluster A

 
• G_CHA

  M. galloprovincialis; protein chaperoning/degradation; cluster A

 
• G_EA

  M. galloprovincialis; energy metabolism; cluster A

 
• HSC

  heat shock cognate

 
• HSP

  heat shock protein

 
• ICDH

  isocitrate dehydrogenase

 
• IPG

  immobilized pH gradient

 
• MOWSE

  molecular weight search

 
• NAD(H)

  nicotinamide adenine dinucleotide (reduced form)

 
• NADP(H)

  nicotinamide adenine dinucleotide phosphate (reduced form)

 
• NHE-RF

  Na⁺/H⁺ exchange regulatory factor

 
• PC

  principal component

 
• PCA

  principal component analysis

 
• PDH

  pyruvate dehydrogenase

 
• PDI

  protein disulfide isomerase

 
• PMF

  peptide mass fingerprint

 
• PTM

  post-translational modification

 
• Rab-GDI

  Rat Brain (small GTPase)-GDP dissociation inhibitor

 
• Ras

  Rat-sarcoma (small GTPase)

 
• Rho

  Ras-homology (small GTPase)

 
• ROS

  reactive oxygen species

 
• RSH

  radial spoke head

 
• SAR

  Secretion-associated Ras-like (small GTPase)

 
• Sec61

  ER protein transport protein

 
• SOD

  superoxide dismutase

 
• TCP-1

  T-complex protein 1

 
• T_CA

  M. trossulus; cytoskeleton-associated proteins; cluster A

 
• T_CHA

  M. trossulus; protein chaperoning/degradation; cluster A

 
• T_EA

  M. trossulus; energy metabolism; cluster A