Rapid proteomic responses to a near-lethal heat stress in the salt marsh mussel Geukensia demissa

Temperature, through its impact on the function and structure of biomolecules, has a pervasive influence on the distribution and survival of organisms (Pörtner, 2002; Hochachka and Somero, 2002; Angilletta, 2009). To survive in thermally variable habitats, or to invade new ones, eurythermal ectotherms must compensate for the changes in macromolecular stability, especially in proteins, caused by variation in temperature (Fields, 2001; Fields et al., 2015). The type of response to temperature stress an organism mounts will vary depending on the time scale involved: over evolutionary time periods, changes mainly occur at the genetic level (i.e. adaptation), leading to alterations in protein structure and function (Somero, 1995, 2004). However, when temperature changes occur relatively rapidly with respect to the lifespan of the individual, organisms will rely on acclimatory processes to alter physiological sensitivity to temperature (Peck, 2011; Tattersall et al., 2012), for example, by altering abundance of particular proteins (Hazel and Prosser, 1974; Shaklee et al., 1977; Somero, 2004), expressing alternate protein isoforms (Baldwin and Hochachka, 1970; Zakhartsev et al., 2007) or modifying protein function via post-translational modifications (PTMs) (Walsh, 2006). Despite the large body of research that has shown that acclimatory changes in protein expression in metazoans can occur quickly, studies of proteomic responses to stress often examine changes in protein expression hours after cessation of the stress (e.g. Gardeström et al., 2007; Tomanek and Zuzow, 2010; Fields et al., 2014). Thus, relatively little is known about the earliest proteomic responses to acute heat exposure in ectotherms.

To more fully describe the early time course of changes in protein expression in response to acute heat stress, we have examined proteomic responses of gill tissue from Geukensia demissa (Dillwyn 1817) (formerly Modiolus demissus; Mytilidae), a mussel found in the high intertidal zone of salt marsh habitats along the east coast of North America (Gosner, 1978). Many temperate intertidal invertebrates regularly must withstand large and sometimes unpredictable shifts in temperature as they are repeatedly exposed to submersed then aerial conditions, and thus these organisms are excellent study systems for examining biochemical and physiological responses to rapid temperature change. Geukensia demissa is especially suited to studies of environmentally induced stress, including heat, cold, hypoxia and salinity, because in its habitat it often is emersed for longer periods than it is immersed (Kuenzler, 1961). For example, G. demissa can withstand temperatures as high as 45°C in the summer (although with some mortality; Jost and Helmuth, 2007), and temperatures at least as low as −22°C in the winter (Kanwisher, 1955), indicating that the species has considerable capacity to respond to temperature variability.

To better understand how G. demissa, and eurythermal poikilotherms in general, survive rapid, large changes in temperature, we acutely exposed mussels in the laboratory to a temperature near the upper thermal limit determined in the field, and allowed them to recover for varying periods of time, from 0 to 24 h, before they were euthanized. We then performed an analysis of protein expression changes in gill tissue. We describe here the rapid changes in protein expression profiles (PEPs) that occur after a near-lethal heat stress, as well as the classes of proteins that are predominantly expressed at different times during recovery. Our results provide insight into the biochemical and cellular mechanisms that allow organisms in thermally variable environments to quickly protect themselves from heat-induced damage, to repair cellular and tissue damage caused by rapid temperature increases, and to prepare for future temperature fluctuations.

Geukensia demissa were collected from Stone Harbor, NJ, USA (latitude 39.039; longitude −74.776), in May 2012, under New Jersey Division of Fish and Wildlife Scientific Collecting Permit no. 1233. Animals were chilled and transported to Franklin & Marshall College, where they were kept immersed in artificial seawater (ASW; 33–35 ppt; Instant Ocean, Marineland Labs, Mentor, OH, USA) in recirculating aquaria at 15°C. Mussels were fed Instant Algae Shellfish Diet 1800 (Reed Mariculture, Campbell, CA, USA) at a volume of 0.1 ml per individual, three times per week. Mussels were allowed to acclimate for 3 weeks before acute heat exposure. There was no mortality during collection, transport and acclimation. No approval of the research protocol was required from the Franklin & Marshall College Institutional Animal Care and Use Committee.

We monitored body temperature of mussels by drilling a small hole in one valve, approximately 1 cm from the posterior edge, and inserting a K-type thermocouple wire into the mantle cavity; the wire was held in place by cyanoacrylate glue. Treatment mussels were divided into groups of six and acutely exposed to 45°C in 800 ml jars filled with approximately 4 cm ASW-saturated sand. To mimic heat stress as experienced by mussels in the field, animals were exposed to air and buried vertically to approximately one-quarter of their body length, anterior down. Jars were capped and immersed in pre-heated aquaria, with the water temperature controlled so that mussels reached 45±1°C in 25 min. Target body temperature was maintained for 1 h, and mussels were either immediately euthanized (0 h group) or re-immersed in the 15°C aquaria to recover for 3, 6, 12, 18 or 24 h before they were euthanized. A seventh group of mussels acted as a control, and was kept immersed at 15°C until euthanization.

Extraction, purification and separation of proteins were performed as described previously (Fields et al., 2014). Briefly, gill tissue was excised from each mussel at its specified time point, and each sample was processed and analyzed individually. Proteins were extracted and purified using trichloroacetic acid/acetone precipitation, and were resolubilized in a denaturing buffer. Proteins were separated using two-dimensional (2D) gel electrophoresis. In the first dimension, proteins were separated by isoelectric point (pI) using isoelectric focusing along an 11-cm pH 3–10 gradient. In the second dimension, after carbamidomethylation of cysteinyl residues, proteins were separated by relative molecular mass (M_r) using gradient sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE; 8–16%). Gels were stained with colloidal Coomassie Blue.

Gel images were recorded at 600 dpi resolution using an Epson model 700V transilluminating scanner. Images were preprocessed using the histogram tool within the Epson software to maximize contrast, setting γ=1 to maintain linearity in pixel density, and were cropped to the gel boundaries using Photoshop software (Adobe Systems, New York, NY, USA). Gel images were imported into Delta2D image analysis software (DECODON, Greifswald, Germany), where images were sorted into treatment groups, warped (i.e. corresponding spots from separate gels were associated) and fused. The fusion image was used to detect spots, and after manual editing to remove smears and false spots (e.g. dust), spot outlines were transferred back to the original gel images. (See Fig. S1 for images of representative gels from each group.) The normalized relative spot volumes of these proteins (i.e. the percentage of the total protein volume on each gel ascribed to each spot) were used for statistical analysis of changes in protein abundance among treatments.

Within the multi-experiment viewer module of Delta2D, an ANOVA test (P<0.05; significance determined by permutation; false positive identifications not exceeding 5%) was used to identify spots changing significantly in abundance across all treatment groups, and Dunnett's multiple comparisons test (GraphPad Prism, GraphPad Software, La Jolla, CA, USA) was used to determine which proteins differed significantly in relative spot volume compared with controls. Hierarchical cluster analysis (HCL), using spots identified as changing significantly via ANOVA and utilizing average linking and Pearson's correlation coefficient, allowed clustering of protein spots according to PEPs. Clustering of mussels by PEPs was further examined via principal components analysis (PCA).

Those protein spots that were determined to have changed significantly in abundance were digested with modified porcine trypsin (Promega, Madison, WI, USA) following standard protocols (Fields et al., 2014). Tryptic peptides were separated and analyzed using HPLC–electrospray ionization ion trap MS/MS (Agilent 1100 series SL ion trap LC/MS; Santa Clara CA, USA). The mobile phase consisted of an increasing linear gradient of acetonitrile in water, acidified with formic acid (0.1%); peptides were separated on a nonpolar C8 column (Agilent) before injection into the trap. Peptide mass lists were submitted to Mascot MS/MS Ions search software (version 3.1; Matrix Science Inc., Boston, MA, USA) for identification. Carbamidomethylation of cysteinyl residues (fixed) and oxidation of methionyl residues (variable) were the only peptide modifications employed during the searches; precursor mass tolerance was set to 1.2 Da and MS/MS tolerance was 0.6 Da, the defaults recommended by Mascot.

For bioinformatic analysis, we used a G. demissa gill expressed sequence tag (EST) library (Fields et al., 2014) containing 47,736 unique entries as the database against which Mascot searched for matches. A protein identification was accepted if the molecular weight search (Mowse) score exceed the significance threshold (P<0.05) and the match included at least two non-overlapping peptide fragments. ESTs that matched MS/MS peptide sequences were then identified via BLAST search (tblastn) of mollusk (taxon ID: 6447) nucleotide sequences in the National Center for Biotechnology Information (NCBI) database.

Acute exposure of G. demissa to 45°C, both in the field and the laboratory, has been shown to cause mortality (Jost and Helmuth, 2007), and during the recovery period of the current experiment, five heat-exposed mussels died (approximately 12%). This left 37 mussels for analysis, distributed across the experimental groups as follows: control, 6; 0 h, 6; 3 h, 5; 6 h, 5; 12 h, 4; 18 h, 5; and 24 h, 6.

After spot editing to remove spurious marks from dust, molecular weight markers and streaks associated with the isoelectric focusing electrodes, Delta2D software detected 1957 protein spots on the fusion image. To remove from analysis those spots that were too faint to be identified by mass spectrometry, we excluded any spot that had a mean spot volume less than 0.05 in each group, after which 933 protein spots remained.

To examine how protein expression in G. demissa gill changed over time after acute heat exposure, we conducted a one-way ANOVA on the expression profiles of the 933 protein spots. A total of 115 spots (12.3%) were found to differ significantly among groups, and a post hoc (Dunnett's) test indicated that many differed from the control at more than one time point. The numbers of spots changing significantly relative to the control in each group are given in Table 1, and it is clear that the magnitude of proteomic response increases to a maximum at 6 h before attenuating. This is corroborated by an HCL of the 115 spots using Pearson correlation and average linkage clustering (Fig. 1), where the 6 h group has the largest number of proteins with changes in abundance (both increasing and decreasing), followed by the 12 and 3 h groups.

To better understand the magnitude of differences in protein expression among groups, we performed PCA on the 115 spots that changed abundance significantly according to ANOVA, in order to simplify the multi-dimensional data set while preserving variance among samples. In Fig. 2A, PC1 (which explains 24.3% of the variation present) clearly separates the 3, 6 and 12 h groups from all other treatments, while the 24 h group is less fully separated in the opposite direction. In addition, PC2 (16.3% of variation) separates the 6 h group from the combined 3 and 12 h groups. The pattern shown in Fig. 2A confirms the results of the HCL (Fig. 1), i.e. that the 6 h recovery mussels showed the most dramatic changes in PEPs relative to the control, followed by the 3 and 12 h groups. In contrast, control, 0 h and 18 h groups cluster together, indicating smaller differences in PEPs among these mussels. Along PC3 (Fig. 2B; 10.8% of variation), the 3 and 12 h groups are differentiated, with the 12 h mussels being most divergent. PC4 (6.0%) separates the 0 h mussels from others, while higher PCs are associated with individual rather than group variation (data not shown). PCA using all 933 detected spots (i.e. including the 818 spots not selected by ANOVA as differing significantly in abundance in response to recovery time), although more strongly affected by inter-individual variation in protein expression than the PCA shown in Fig. 2, supports the separation of the 6 and 12 h groups from the other treatments (Fig. S2).

To explore the biochemical strategies used by G. demissa to avoid or recover from heat-induced damage, we attempted to identify each of the 115 protein spots that changed significantly in abundance after heat exposure. Geukensia demissa is a non-model organism with limited sequences available in public protein and nucleic acid databases. Using a G. demissa gill EST library combined with homology searching of mollusk genomic and EST databases at NCBI, we were able to identify 85 of the 115 spots (73.9%), which are listed in Table S1. Fourteen proteins were identified more than once, including identical GenBank accession numbers, and a total of 35 spots shared an identification with at least one other spot. In cases where spots with synonymous identifications were contiguous on the 2D gels, the spots were combined to form a single spot. In some cases, however, the same protein was identified in spots at separate locations on the gel. In most cases, the M_r of the spots with matching identifications was the same (within the resolution of the gel), but the pI values varied. In these cases, which we interpret as evidence of a PTM altering the charge of the protein, the protein abundances were analyzed separately for each spot. After coalescing contiguous spots and applying Dunnett's multiple comparisons test to find proteins with abundances significantly different from the control, the 85 identified proteins with significant changes in abundance were reduced to 53 unique spots (Table S1).

As is evident from HCL and PCA (Figs 1, 2), few proteins showed changes in abundance immediately after heat exposure, likely because of the short duration and extreme nature of the stress. Nevertheless, 11 spots changed significantly in abundance at 0 h (i.e. during the stress, rather than the recovery period), of which seven were identified, corresponding to three proteins: prohibitin, voltage-dependent anion channel 2-like protein (VDAC2) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) (Fig. 3). Note that the latter two are represented on the 2D gels by three spots with differing pI values, suggesting that changes in PEPs of these proteins likely are due at least in part to PTMs. As described in the Discussion, a significant change in abundance of each of these proteins at 0 h suggests that the earliest response to damaging heat stress is induction of pro- or anti-apoptotic responses.

The number of proteins that changed in abundance grew substantially at 3, 6 or 12 h of recovery. At 3 h, 28 spots showed significant changes in abundance relative to control, of which 18 spots corresponding to 14 proteins were identified, and five of these can be categorized as chaperones or foldases (Fig. 4A): 78 kDa glucose regulated protein (GRP78), small heat shock protein 24.1 (sHSP24.1), peptidylprolyl cis-trans isomerase (PPIase), protein deglycase DJ-1 and universal stress protein A-like (USP-A). In addition to these stress proteins, other proteins changing in abundance for the first time at 3 h (Fig. 4B) are actin 2, mitochondrial H⁺ ATPase subunit α and nacre protein.

The 6 h recovery mussels show the greatest change in protein expression, as evidenced by HCL and PCA (Figs 1, 2); 44 spots changed significantly in abundance at this time point relative to controls. Among the 23 proteins identified here are a further subset of stress proteins – peroxiredoxin 6 (Prx-6), dyp-type peroxidase, nucleoredoxin and NADP-dependent isocitrate dehydrogenase (Fig. 5A) – suggesting that during recovery from heat stress, oxidative damage is addressed after protein misfolding and aggregation are dealt with. Furthermore, at 6 h, a number of glycolytic and citric acid cycle (CAC) enzymes change in abundance (Fig. 5B), including fructose-bisphosphate aldolase (ALD), triosephosphate isomerase (TPI), mitochondrial malate dehydrogenase (mMDH) and citrate synthase (CS), as well as arginine kinase (AK), which in mollusks produces the energy storage phosphagen phosphoarginine. It should be noted as well that GAPDH, which remains increased in abundance at 6 h (Fig. 3), also is a glycolytic enzyme. Other proteins increased in abundance at 6 h are hemicentin-1 and profilin, N(G), N(G) dimethylarginine dimethylaminohydrolase 1 (DDAH1), vitelline membrane outer layer protein 1 (VMO-1), heterogeneous nuclear ribonucleoprotein 27C (hnRNP27C) and a C1q-motif containing protein (Fig. 5C).

Of the 34 significantly different spots detected at 12 h, 22 were identified, and grouped into 14 proteins. Four of these identified proteins show significant changes in abundance at 12 h for the first time (Fig. 6), including putative perlucin-4 and perlucin-like protein, each of which is a component of the protein matrix of the shell; universal stress protein-like isoform 2 (USP-2); and alcohol dehydrogenase 3 (ADH3). It should be emphasized, however, that a number of proteins described at earlier time points maintain significant changes in abundance at 12 h, for example, GAPDH (Fig. 3), PPIase, DJ-1 and USP-A (Fig. 4), and Dyp-type peroxidase, hnRNP27C and VMO-1 (Fig. 5).

At 18 h and beyond, the number of spots showing significant differences in abundance relative to controls decrease substantially, suggesting that much of the proteomic response to acute heat stress has been accomplished by this time, and PEPs are reverting to the pre-stress pattern. At 18 h, 13 spots differed in abundance from control, but only two new proteins show significant changes in abundance (Fig. 7), although again a number of proteins that were described earlier maintain significant differences in abundance at 18 h. These two proteins are cathepsin B, a cysteine protease with broad specificity that is normally associated with the lysosome (but see Discussion), and flotillin-2, a protein associated with microdomains of the plasma membrane.

The number of spots changing in abundance at 24 h is comparable to that at 18 h (14 versus 13). The five proteins identified with changes in abundance for the first time in this group (Fig. 8) were the mitochondrial oxidative stress protein Prx-5, nucleoside diphosphate kinase (NDPK), 15-hydroxyprostaglandin dehydrogenase (15-PGDH), 5′-AMP-activated protein kinase (AMPK) and eukaryotic translation initiation factor 5A-1-like (eIF5A). Interestingly, PCA suggests that the PEP of the 24 h group is somewhat different from that of the other groups: PC1 (Fig. 2A) separates 24 h mussels from an undifferentiated cluster of groups including control, 0 h and 18 h groups. One potential explanation for the separation of the 24 h group, discussed below, is that changes in protein expression at this time are signaling a return to an anabolic state as the stress response diminishes.

We believe that the present study helps fill a gap in our understanding of proteomic responses to acute environmental stress in animals. Whereas earlier studies examining recovery from stress often have examined changes in protein abundance at a single time point (e.g. Gardeström et al., 2007; Tomanek and Zuzow, 2010; Serafini et al., 2011; Fields et al., 2012b; Garland et al., 2015; but see Dowd et al., 2010), here we have examined the progression of changes in abundance over a series of time points after acute heat stress. Our results, we argue, may help to form a clearer and more complete description of how organisms cope with potentially damaging stress at a physiological and biochemical level. In G. demissa, changes in protein expression induced by an acute and near-lethal heat stress occur rapidly. Even at 0 h post-stress, a small number of potentially significant changes in protein abundance are apparent. At 3 h, the number of proteins whose abundance changes relative to control levels increases, and at 6 h the proteomic response to stress appears maximal (Table 1, Figs 1, 2). After 6 h, the magnitude of proteomic response decreases substantially. Importantly, not only is there a pattern of increasing changes in protein expression over the first 6 h of recovery, but there also appears to be a series of coordinated biochemical responses to the stress as time progresses.

The identified proteins that changed in abundance at the 0 h time point (Fig. 3), and thus which must have been synthesized or modified during the period of heat stress, provide a signal of changes to apoptotic signaling pathways. Prohibitin, which increased 2.1-fold at 0 h before returning to control levels at 3 h, is found in a number of cellular compartments and has been ascribed diverse roles, including as a mitochondrial chaperone (Nijtmans et al., 2000) and in regulation of transcription and cell cycle progression (Artal-Sanz and Tavernarakis, 2009). Most importantly, decreased prohibitin protein levels appear to induce apoptosis (Sánchez-Quiles et al., 2010), perhaps through greater sensitivity to oxidative damage (Theiss et al., 2006), and the 3′-untranslated region of prohibitin mRNA has been identified as a tumor suppressor with anti-proliferative activity (Song et al., 2014). Based on these findings, we suggest that the rapid increase in prohibitin abundance in G. demissa during heat stress may be a protective response to delay or avoid heat-induced apoptosis.

Like prohibitin, the abundance of VDAC2 also changes rapidly during heat exposure, with significant effects observable at 0 h. However, unlike prohibitin, VDAC2 appears on the 2D gels as three spots with similar M_r values but differing pI values, and each of the spots has a unique expression profile (Fig. 3). The spot with the highest pI increases 1.7-fold in abundance at 0 h relative to control, but by 3 h has decreased to control levels. In contrast, the spot with the lowest pI shows significantly increased abundance at 3 h only, and the spot with an intermediate pI value is significantly increased at 6 h only. This pattern of shifting pI with time is likely a result of PTMs that alter protein charge while minimally impacting M_r. VDAC2 is an important component of the mitochondrial outer membrane permeability pore, which allows movement of ions and metabolites between the cytosol and the intermembrane space (Mertins et al., 2014). [Note that despite the nomenclature, it appears that bivalves have only one VDAC isoform, whereas mammals have three (Shoshan-Barmatz et al., 2010).] During periods of cellular stress, however, VDAC impacts apoptotic processes (Shimizu et al., 1999), although there is controversy regarding whether the main role of VDAC is pro- or anti-apoptotic. Growing evidence in mammals suggests that VDAC overexpression leads to apoptosis, perhaps through induction of the mitochondrial permeability transition and release of pro-apoptotic cytochrome c into the cytosol (Shoshan-Barmatz et al., 2010). In contrast, recent work in the oyster Crassostrea gigas indicates that bivalve VDAC2 inhibits UV-induced apoptosis, likely through binding and inhibition of the pro-apoptotic Bak protein (Li et al., 2016).

The third protein increased in abundance at 0 h is GAPDH, most familiar as a glycolytic enzyme involved in ATP generation. However, GAPDH is known to have a number of ‘moonlighting’ roles beyond glycolysis, including as an oxidative stress sensor whose inactivation can shunt glucose into the pentose phosphate pathway, as a transcription factor in the nucleus, and as a pro- or anti-apoptotic factor depending on its location (Tristan et al., 2011). Like VDAC2, GAPDH in G. demissa appears as three discrete spots with differing pI values (Fig. 3), which may represent PTMs that activate different functions. The highest-pI isoform decreases significantly in abundance at 0 h and remains lower than control for the entire 24 h recovery period. In contrast, the middle- and lowest-pI isoforms increase in abundance, but not until 3 or 12 h, respectively.

Given the roles of prohibitin and VDAC in apoptosis, it is interesting to note that GAPDH has also been found to impact apoptotic processes in a variety of cellular locations. For example, cells undergoing apoptosis show increased levels of GAPDH in the nucleus (Dastoor and Dreyer, 2001), where GAPDH may act as part of a complex that causes proteolytic cleavage of nuclear proteins, among other pro-apoptotic functions (Chuang et al., 2005). However, GAPDH also protects telomeres from rapid degradation during periods of oxidative stress, thereby inhibiting apoptosis (Sundararaj et al., 2004; Demarse et al., 2009). In the mitochondria of mammals, GAPDH has been shown to accumulate during stress, interact with VDAC1 and induce mitochondrial permeabilization, thus leading to apoptosis (Tarze et al., 2007).

In summary, the PEPs of prohibitin, VDAC2 and GAPDH strongly indicate that apoptotic signaling pathways are activated immediately upon exposure to near-lethal heat stress. However, because of the complexity of apoptotic pathways, and the varying roles that have been ascribed to each of the proteins, it is not clear whether the changes in the abundances represent a pro- or anti-apoptotic signal, or perhaps both, depending on the amount of molecular damage experienced by different cells.

The universal response to cellular stress involves expression of chaperones to protect and refold proteins undergoing denaturation, or to send irretrievably damaged proteins into the ubiquitin-proteasome pathway for degradation. In addition, it has become clear that protection from oxidative stress is a component of the minimal stress proteome (Kültz, 2005), and oxidative stress proteins are often increased in abundance after heat stress, including in bivalves (Tomanek and Zuzow, 2010; Fields et al., 2012a). Stress proteins form a substantial portion of the proteomic response to acute heat exposure in G. demissa, but interestingly, changes in abundance do not appear immediately (i.e. at 0 h), instead becoming apparent at the 3 h time point.

Among the stress proteins increased in abundance at 3 h are GRP78, HSP24.1 and PPIase (Fig. 4A). GRP78 is an HSP70 analog found in the endoplasmic reticulum, where it has a central role in the unfolded protein response (Lai et al., 2010) and may help limit the effects of oxidative stress (Liu et al., 1997). Small heat shock proteins such as HSP24.1 are found in the cytosol, where they interact with and stabilize the actin cytoskeleton during stress (Arrigo, 1998; Dalle-Donne et al., 2001; Haslbeck et al., 2005), as well as inhibit apoptosis (Beere, 2004) and prevent the aggregation of denatured proteins (Horwitz, 1992; Merck et al., 1993). PPIase, which is represented by multiple spots on the 2D gels, is a foldase also found in the endoplasmic reticulum that helps refold partially denatured proteins (Wang and Heitman, 2005). One PPIase spot (low pI) increases significantly in abundance at 3, 6 and 12 h while a second decreases in abundance at 3 h, suggestive of PTM. Combined with the increased abundance of prohibitin at 0 h, which has been shown to have chaperoning activity in the mitochondrion (Nijtmans et al., 2000), the changes in expression of GRP78, PPIase and HSP24.1 indicate a rapid increase in chaperoning capacity in a variety of compartments in the cell. Notably, three of these chaperones return to control levels relatively quickly, prohibitin at 3 h, GRP78 at 6 h and HSP24.1 at 12 h. Perhaps surprisingly, one protein whose abundance was not found to change significantly at any time point was cytosolic HSP70, which often is detected in proteomic studies of heat stress in bivalves (e.g. Tomanek and Zuzow, 2010; Fields et al., 2012a), including G. demissa (Fields et al., 2012b). Further research is necessary to determine under what conditions HSP70 levels increase, for example, in response to acute but relatively mild heat stress, versus the near-lethal stress that mussels experienced in this study.

Oxidative stress also elicits a clear proteomic response in G. demissa after acute heat exposure, but the oxidative stress response appears to be delayed relative to the chaperone/protein folding response. Most oxidative stress proteins are first significantly increased in abundance at 6 h, and have returned to control levels by 12 h (Fig. 5A). These include cytosolic Prx-6, which detoxifies lipid peroxides (Schremer et al., 2007; Fisher, 2011); Dyp-type peroxidase, which removes H₂O₂ (Sugano, 2009); and nucleoredoxin, a thioredoxin-family protein that, despite its name, is predominantly cytosolic (Funato and Miki, 2007). In addition, one of the two NADP⁺-dependent isocitrate dehydrogenase spots first increases in abundance at 6 h. Unlike NAD⁺-dependent isocitrate dehydrogenase, which is a CAC enzyme, the NADP⁺-dependent isoform is associated with protection of matrix proteins from oxidative stress, as it shuttles electrons from the CAC intermediate isocitrate to NADPH, which then donates electrons to oxidatively damaged macromolecules (Jo et al., 2001). Taken together, these changes in abundance indicate that G. demissa gill strongly increases oxidative repair capacity in cytosol and matrix, but with the response reaching a maximum only after 6 h of recovery.

In contrast, three proteins associated with oxidative stress response have significant changes in abundance at times other than 6 h. One, protein deglycase DJ-1, is significantly reduced in abundance at 3 h, and its abundance remains significantly below that of the control at 6 and 12 h as well (Fig. 4A). Deglycase DJ-1 repairs methylglyoxal- or glyoxal-glycated (Richarme et al., 2015) and sulfenylated proteins (Gautier et al., 2012) and defends cells against oxidative damage (Wilson, 2011). The decrease in DJ-1 abundance is counterintuitive; however, it is possible that the change in expression represents a PTM, with a concomitant increase in abundance in another spot on the 2D gel that was not detected as significantly increasing in abundance. At 12 h, two ADH3 spots are upregulated (Fig. 6); ADH3 acts as an antioxidant enzyme, using electrons donated by glutathione to eliminate formaldehyde (Danielsson et al., 1994) and S-nitrosothiols (Liu et al., 2001) in proteins that have been oxidatively damaged. Finally, Prx-5, which is increased in abundance at 24 h, detoxifies H₂O₂ in the mitochondrial matrix (Fig. 8).

Two other, more enigmatic stress-associated proteins were identified amongst those changing in abundance in G. demissa after heat exposure. USP-A showed a significant decrease in abundance at 3 and 12 h (Fig. 4A), and USP-2 was decreased at 12 h (Fig. 6). USPs are known to respond to a variety of stresses in prokaryotes (Sousa and McKay, 2001), and more recently it has been shown that they possess chaperoning activity in plants, especially under conditions of oxidative or heat stress (Jung et al., 2015). However, there currently is little information on the function of USPs in animals.

In addition to clear proteomic responses to acute heat exposure in apoptotic pathways and stress proteins, proteins associated with energy homeostasis were strongly affected in G. demissa. Most of these proteins first showed significant changes in PEPs at 6 h (Fig. 5B); for example, among glycolytic enzymes, ALD and TPI are increased in abundance only at 6 h. In addition, the glycolytic enzyme GAPDH shows changes in expression, as noted above (Fig. 3). Although we argue that its initial alteration in expression profile suggests it plays a role in apoptotic pathways during heat stress, the complexity of the changes in GAPDH expression shown in Fig. 3 suggest that it indeed may be taking on a number of different roles, including ATP generation, and some of these changes continue to be significant at 6 h and later.

The CAC enzymes increasing in abundance at 6 h are CS and mMDH (Fig. 5B); both return to control levels by 12 h. One component of oxidative phosphorylation also was detected, mitochondrial H⁺-ATPase α, but unlike the other energy metabolism proteins described here, it alters abundance at 3 h (Fig. 4B), with two separate spots both following similar PEPs. The PEP of AK also changes significantly at 6 h, with a high-pI spot increasing in abundance while a low-pI spot decreases. The pattern reverts to control at 12 h, strongly suggesting a reversible PTM.

Two other proteins identified here that function in cellular energy homeostasis are noteworthy, although their changes in abundance do not occur until 24 h post-stress (Fig. 8). AMPK is known as a ‘master regulator’ of energy homeostasis (Hardie, 2007), and responds to low ATP/AMP ratios by activating ATP-generating pathways while inhibiting anabolic pathways. In heat-stressed G. demissa gill, AMPK levels rise slightly, but not significantly, at 6 and 12 h. However, AMPK levels decrease substantially at 24 h, suggesting that gill tissue is shifting from a catabolic to an anabolic state. Concomitantly, NDPK levels increase (Fig. 8). NDPK transfers the γ-phosphate from a nucleotide triphosphate (NTP) to a nucleotide diphosphate, maintaining balance among fully phosphorylated NTPs in the cell for both catabolic and anabolic purposes. Because ATP is the most prevalent NTP, however, the activity of NDPK can be deleterious during times of stress as it removes needed ATP. NDPK is phosphorylated and inactivated by AMPK in order to favor maintenance of ATP levels during stress (Onyenwoke et al., 2012), and the reciprocal changes in PEPs of AMPK and NDPK in G. demissa gill may indicate a return to ‘normality’ in energy homeostasis at 24 h.

The actin cytoskeleton is one of the earliest targets of heat-induced oxidative stress (Dalle-Donne et al., 2001), and in this experiment three spots on the 2D gels were identified as actin-2, one at a relatively high pI and two others at a lower pI, but with different M_r values (Fig. 4B). The actin spot with the highest expression level increased in abundance at 3 h, but returned to control levels by 6 h. The other two spots appear to have reciprocal expression patterns, with the high-pI form increasing significantly in abundance at 18 h while the low-pI form decreased significantly at 24 h. In addition to actin, other proteins associated with the cytoskeleton suggest that recovery from heat stress includes repair of the actin cytoskeleton and cell proliferation. Small HSP24.1, as described above, increases transiently at 3 h (Fig. 4A), and stabilizes the actin cytoskeleton by regulating actin filament dynamics and preventing aggregation (Dalle-Donne et al., 2001; Wettstein et al., 2012). At 6 h, profilin is increased 2.8-fold (Fig. 5C); profilin also plays a role in stabilization and synthesis of the actin skeleton by affecting rates of actin polymerization and depolymerization (Yarmola and Bubb, 2006). Interestingly, profilin is necessary for abscission during cytokinesis (Böttcher et al., 2009), and also helps control cell motility and migration (Witke, 2004; Pernier et al., 2016).

Confirming the importance of cell proliferation and motility during recovery from heat stress, four other proteins, flotillin-2, hemicentin-1, cathepsin B and DDAH1, link intracellular processes affecting cytoskeletal dynamics to the extracellular matrix (ECM) and cell migration. Flotillin-2, which increases significantly in abundance at 18 h (Fig. 7), has been associated with cholesterol-rich microdomains in the plasma membrane (Simons and Toomre, 2000), where it interacts with the actin cytoskeleton (Langhorst et al., 2007), is involved in cell adhesion and movement through the ECM (Neumann-Giesen et al., 2007), and is essential for appropriate targeting of motile cells, for example during development in Drosophila (Langhorst et al., 2005). A second protein, hemicentin-1, which is increased in abundance at 6 h (Fig. 5C), is an ECM protein that, like profilin, is involved in cell division (specifically, maturation of the cleavage furrow) (Xu et al., 2013). Cathepsin B (18 h) is a cysteine protease often associated with the lysosome, but also has been shown to be expressed extracellularly, where it degrades structural proteins including collagen, laminin and fibronectin, and is central to controlled digestion of the ECM (Buck et al., 1992), thus enabling cell migration via podosomes (Tu et al., 2008; Schachtner et al., 2013). Finally, the increase in abundance of DDAH1 at 6 h is a notable marker of cell proliferation and motility. DDAH1 metabolizes asymmetric dimethylarginine, which itself is an inhibitor of all isoforms of nitric oxide synthase (Tran et al., 2003). Thus, DDAH1 favors synthesis of the signaling molecule nitric oxide, and in mammalian systems has been shown to have pro-proliferative effects on endothelial cells (Fiedler et al., 2009) and enhance motility of trophoblasts (Ayling et al., 2006). Interestingly, DDAH1 has been observed to colocalize with F-actin and nitric oxide synthase in lamellipodia of migrating endothelial cells (Wojciak-Stothard et al., 2007), indicating its central role in enabling cell motility. The changes in abundance of these four proteins as well as profilin, then, provides strong evidence that cell proliferation and migration within gill tissue is occurring during later stages of recovery from heat stress, from 6 h to 18 h.

More speculatively, we suggest that a number of proteins changing in abundance with heat exposure in G. demissa gill, which do not have obvious connections to an abiotic stress response, in fact are associated with the ECM and also are involved in tissue reorganization and regeneration. These include nacre protein (3 h; Fig. 4B), perlucin-4 and perlucin-like protein (12 h; Fig. 6), all of which typically are associated with the organic matrix of the shell (Blank et al., 2003) but here are found in gill tissue; VMO-1 (6 h; Fig. 5C), a component of the egg ECM (Mann, 2007); and putative C1q domain containing protein (Fig. 5C), which usually is ascribed a role in the innate immune response (Gerdol et al., 2011). None of these proteins normally are associated with gill tissue or abiotic stress, but they all play roles (or have domains that play a role, in the case of C1q domain containing protein) in extracellular networks in tissues outside the gill. For example, perlucins have previously been shown to be expressed in gill (e.g. Moreira et al., 2015), and structurally have been found to contain the carbohydrate-recognition domain of a C-type lectin, thus suggesting that perlucins are involved in cell–cell recognition and adhesion (Mann et al., 2000). Similarly, C1q proteins have been associated with shell biomineralization in mussels (Yin et al., 2005), suggesting that these proteins also may play a structural role in the ECM. Thus, we suggest that despite their identification with processes or structures outside the gill, here we have found these proteins, or structurally homologous proteins, which as part of the ECM of the gill are changing in abundance along with other structural proteins as gill tissue is repaired and remodeled.

As noted above, the PEPs of two proteins, AMPK and NDPK, suggest that as far as energy homeostasis is concerned, at 24 h G. demissa gill may be returning to an anabolic state after a stress-induced switch to ATP-generating pathways. Two other proteins that increase in abundance at 24 h, eIF5A and 15-PGDH (Fig. 8), are more difficult to categorize. Although the function of eIF5A is not completely understood, the protein either is found in the cytosol, where its activity as an elongation factor is required for cell proliferation, or it localizes to the nucleus and helps regulate apoptosis, depending on the presence or absence of a unique hypusine PTM (Caraglia et al., 2001; Lee et al., 2009). Interestingly, decreased levels of active (hypusinated) eIF5A inhibit cell growth (Nishimura et al., 2005), and we propose that the changes in eIF5A we have found in G. demissa at 24 h may indicate a relief of that inhibition as gill recovers from heat stress. 15-PGDH is responsible for the first step in the degradation of prostaglandins (PGs) by oxidizing the hydroxyl group at C15 (Tai et al., 2002), for example, in PGE₂. In mollusks, PGs are most often associated with sexual maturation and gamete release (Deridovich and Reunova, 1993), and it is not clear why 15-PGDH abundance would change after heat stress if this was their only function. However, it has been shown in G. demissa that PGE₂ levels are significantly increased during osmotic stress (Freas and Grollman, 1980), suggesting that this PG also plays a role in stress response, and thus the increase in abundance of 15-PGDH found in the present study may signal the attenuation of this stress-mediated hormonal signal once the tissue has recovered.

The relative timing of changes in abundance of proteins associated with apoptosis, protein folding, oxidative stress and energy metabolism is intriguing given our current understanding of translational control during acute stress. It is clear that severe heat stress strongly inhibits protein translation (Lindquist, 1980), and it is notable that changes in protein abundance at 0 h, particularly of VDAC2 and GAPDH, are consistent with PTMs, rather than protein synthesis (Fig. 3). Although the translation of most proteins is suppressed during severe stress, some proteins, notably chaperones, continue to be synthesized at a relatively high rate (Hartl et al., 2011), which may account for the early appearance of chaperones in our data (Fig. 4A). A number of mechanisms are responsible for the preferential expression of stress-related proteins, including selective arrest of elongation of nascent polypeptides on the ribosome (Shalgi et al., 2013) and selective inhibition of splicing of mRNAs in the nucleus (Shalgi et al., 2014). Notably, it appears that proteins specifically related to stress response, and especially those involved in protein folding, oxidative stress and energy production, are efficiently spliced and exported from the nucleus during severe heat stress, while other mRNAs stall in the pre-splicing stage (Shalgi et al., 2014). These strategies of transcriptional control may help explain the time course of protein abundance changes we have found, and the relative pre-dominance of chaperones, oxidative stress proteins, and glycolytic and CAC enzymes, at least during the first 6 h of recovery.

Equally importantly, our study design may have impacted the timing of expression of certain classes of stress proteins. This is because the mussels were stressed in air, but were immersed for recovery. Although protein denaturation would be expected to occur during the period of heat stress, and cease upon return to the lower recovery temperature, thus necessitating rapid expression of chaperones, oxidative stress might be relatively delayed. This is because a substantial portion of oxidative damage could occur upon return to the recovery tank, when the sudden exposure to a pulse of oxygenated water after a period of hypoxia might lead to reperfusion-like oxidative stress (Welker et al., 2013), helping to explain why changes in abundance of oxidative stress proteins lag behind those of chaperoning proteins.

In conclusion, the data presented here suggest that G. demissa reacts to acute heat stress over time with a series of discrete and coordinated proteomic responses. We have found that the earliest response involves activation or inhibition of apoptotic pathways, followed by macromolecular repair. Later, oxidative stress proteins change in abundance along with ATP-generating pathways. Still later, tissue-level repair is evident through changes in proteins associated with the cytoskeleton, cell proliferation and migration, and the ECM. Finally, 24 h after stress gill tissue appears to shift from the stress-responsive, catabolic state to an anabolic state. In future studies it will be interesting to determine whether a similar series of proteomic changes occur in other ectotherms, and in response to stresses other that heat.

We would like to thank Dr Ken Hess for his help in the mass spectrometry portion of this study.