Changes in protein expression in the salt marsh mussel Geukensia demissa: evidence for a shift from anaerobic to aerobic metabolism during prolonged aerial exposure

Residents of the intertidal zone experience cyclical immersion and aerial exposure, and the stark contrast in physical conditions between these two environmental states leads to unique physiological and biochemical stresses. Nearly all intertidal species have arisen from marine rather than terrestrial ancestors, and so stresses associated with the tidal cycle usually occur during aerial exposure. These stresses may include hypoxia, desiccation, heat or cold, UV exposure, and hyposaline stress due to freshwater runoff (Sagarin et al., 1999; Helmuth et al., 2006; Somero, 2012). Indeed, such physiological stresses often determine the upper vertical limit of a species' intertidal distribution (Tomanek and Somero, 1999; Stillman, 2002) and many intertidal species experience environmental conditions that approach their physiological limits (Stillman, 2003; Jost and Helmuth, 2007). Because of the magnitude and rapidity of changes in the physical environment of the intertidal zone, intertidal organisms have been used for many years as model organisms by researchers interested in physiological and biochemical adaptation to stressful environments.

Geukensia demissa (Dillwyn) (Mytilidae), the ribbed mussel widespread in salt marshes along the east coast of North America from southern Florida to the Gulf of St Lawrence (Blackwell et al., 1977; Gosner, 1978), is one such species that must withstand extreme variation in environmental conditions across a tidal cycle, throughout a day, and over the course of a year. It is the bivalve found highest in the intertidal zone throughout most of its range (Kuenzler, 1961; Lent, 1969), likely due to predation pressure from the blue crab, Callinectes sapidus (Lin, 1989), and is often emersed for longer periods than it is immersed (Kuenzler, 1961), spending up to 80% of the tidal cycle in air (Hilbish, 1987). It thus is repeatedly exposed to a wide range of physical stressors, and is remarkably tolerant of them. Geukensia demissa survives temperatures of 45°C in the summer in the warmer parts of its range (Jost and Helmuth, 2007), as well as temperatures as low as −22°C in the winter in the north (Kanwisher, 1955), perhaps through the use of bacteria as ice-nucleators in its pallial fluid (Loomis and Zinser, 2001). Long periods of aerial exposure also lead to significant hypoxia, and G. demissa is known to air-gape, likely to allow aerial respiration (Kuenzler, 1961; Lent, 1968; Hilbish, 1987); however, gaping leads to desiccation and G. demissa has been reported to survive the loss of 38% of its soft-tissue water content (Lent, 1968).

Geukensia demissa also is an ecologically important and abundant member of the salt marsh community along the east coast of North America. The mussel is found embedded among the rhizomes of the ubiquitous Spartina alterniflora cordgrass that dominates the low marsh, and the byssus of the mussel helps stabilize the organic network that allows the mudflat to accrete; in the absence of the mussel the salt marsh erodes (Bertness, 1984; Smith and Frey, 1985). Jordan and Valiela calculated that G. demissa filter a volume of water greater than the tidal volume of the marsh each tidal cycle (Jordan and Valiela, 1982), and Bertness found that the nitrogenous wastes deposited by the mussel increase cordgrass productivity by approximately 50% (Bertness, 1984). Given the importance of salt marsh habitat as a nursery for estuarine and marine species, as well as its function as a buffer against shoreline erosion, the substrate stabilization and enhancement of productivity offered by G. demissa are ecologically and economically significant.

In the present study we have taken a proteomic approach to more fully describe the physiological mechanisms by which G. demissa withstand long periods of aerial exposure. Although temperature and UV stresses are correlated with emersion, we held these factors constant to focus on aerial exposure alone, comparing responses to increasing periods of emersion and decreasing periods of post-emersion recovery. We chose experimental conditions to mimic short-duration exposure during neap tides (6 h exposure, 18 h immersion recovery; hereafter referred to as 6E/18R), longer term exposure (12 h exposure, 12 h recovery; 12E/12R), and exposure conditions that might be experienced by mussels highest in the marsh during extreme spring tides (18 h exposure, 6 h recovery; 18E/6R).

We chose to vary exposure and recovery periods concomitantly for two reasons. First, the study design follows the patterns of exposure mussels experience in the field, where emersion and immersion periods are inversely related, but the duration of the entire tidal cycle is constant. Second, it is likely that substantial stress-related protein expression changes occur during the aerial exposure period as well as during recovery, and if exposure or recovery time alone were varied, the total treatment time (i.e. exposure + recovery) would differ among the groups. This would give the groups different amounts of time to initiate stress-related transcription and translation, likely affecting protein abundance and protein expression patterns (PEPs) and thereby obscuring the effects of emersion stress.

Within the gill tissue of G. demissa, our results show a clear signal of increased anaerobic capacity after the shortest period of aerial exposure, but as emersion duration lengthens, evidence of upregulation of fermentative metabolism disappears, and the proteomic response shifts toward oxidative stress protection. After the longest exposure, changes in cytoskeletal components predominate, which may be associated with repair or modification of gill structures after air gaping and desiccation.

We prepared two-dimensional polyacrylamide gels for 24 mussels, six replicates from each of the three treatment groups plus six controls that did not experience aerial exposure. After fusing the gel images of the 24 mussel gill protein samples and manually removing spots associated with electrode smears and the molecular weight ladder, Delta2D software detected 800 protein spots (see supplementary material Fig. S1 for representative gels from each group). Once we had filtered spots using a low-volume cut-off (below 0.05% average spot volume in all groups), 363 spots remained, which were used for the analyses described below.

To identify those protein spots that most significantly contribute to the differences in protein expression among groups after the aerial exposure treatments, we performed principal components analysis (PCA) on the 24 samples using data from all 363 protein spots (Fig. 1). The four treatment groups cluster separately, indicating that the treatments significantly and replicably alter G. demissa gill PEPs. The 6E/18R group is most separated from the others on the first principal component (PC1; horizontal axis in Fig. 1A,B). By definition, PC1 accounts for the greatest amount of variation in protein expression among the samples: 15.2% of variation in the PEPs of all 363 spots is contained within PC1. PC2, the vertical axis in Fig. 1A, accounts for 11.8% of variation in PEPs, and this component most clearly separates 12E/12R from 18E/6R. Along PC3 (vertical axis in Fig. 1B; 8.3% of expression variation), the control group is separated from the three treatment groups, especially 12E/12R and 18E/6R, which cluster loosely together. Along higher PCs, samples do not cluster by treatment group (data not shown); these higher components are associated with inter-individual differences in protein expression. The first three PCs explain 35.2% of the variation in gill protein expression among the mussels, and although this percentage may appear low, we note that all 363 detected proteins are included in the analysis, not only those that were shown to change significantly in abundance in response to treatment (cf. Tomanek and Zuzow, 2010).

We used the component loadings derived from PCA to rank individual protein spots and thus determine their relative contributions to the differences apparent among the treatment groups in Fig. 1. Ordering by positive component loading along PC1 reveals which proteins contribute most to separating 6E/18R from the other three groups (Fig. 1A). Similarly, along PC2, proteins most affected by 12E/12R treatment have high positive loading values, while proteins whose abundance changes most after 18E/6R treatment have the most negative loadings. Finally, those protein spots with greater abundance in control mussels relative to exposure treatments have higher negative loading values along PC3 (Fig. 1B), where the control group is most clearly differentiated from the others.

To better understand how protein expression changes led to the separation of treatment groups in Fig. 1, we subjected the 20 spots with the highest component loadings from each group to hierarchical cluster analysis (HCL); Fig. 2 shows the expression profiles of these 80 proteins. Mussels within each treatment group cluster together (along the top of the heat map), and there are discrete protein clusters associated with each mussel treatment group (left); the function of proteins in each cluster is described below. Of the 80 spots analyzed, 70 were found to differ significantly in abundance (ANOVA; P>0.02; significance determined by permutation) (Fig. 2).

After trypsin digestion and MS/MS analysis of the 20 spots with highest loadings from each group, we were able to identify 66 proteins. [Six of the 14 unidentified spots produced significant matches to our G. demissa EST library, but produced no match (‘no significant similarity’) or a match to a hypothetical protein in the NCBI nr database.] Proteins are listed in Tables 1, 2, 3 and 4 according to component loading; more detailed information about the identification of each protein can be found in supplementary material Table S1. The 70 spots that differed significantly in abundance among the groups also are indicated in Tables 1, 2, 3 and 4.

Positive component loadings on PC1 are most closely associated with proteins increasing in abundance in group 6E/18R (Fig. 1), and proteins associated with catabolic processes are well represented in this group (Table 1, Fig. 3). Seven of the 18 identified proteins with high positive loadings on PC1 fall into this category, including succinate dehydrogenase/fumarate reductase (SucDH/FRD) subunit B, propionyl-CoA carboxylase (PCC), fumarate hydratase (fumarase), phosphoenolpyruvate carboxykinase (PEPCK) and cytosolic malate dehydrogenase (cMDH). Notably, these five enzymes are central to succinate/propionate fermentation (Fig. 3A) (see Discussion). The final two enzymes associated with catabolic processes in the 6E/18R group are a putative methylmalonate-semialdehyde dehydrogenase and mitochondrial glutamate dehydrogenase (Fig. 3B). After 6 h of aerial exposure, only two oxidative stress proteins were found among the 20 with the highest positive component loadings on PC1, mitochondrial peroxiredoxin (Prx)-5 and a transketolase-like protein (Fig. 3B). Three proteins associated with the cytoskeleton also are present – adenylate cyclase-associated protein 1, actin-2 and tubulin α1C (Fig. 3B).

Proteins significantly increasing in abundance in the 12E/12R group can be identified through high positive loadings on PC2 (Fig. 1A). Eight of the 18 spots identified among this group are oxidative stress proteins (Table 2, Fig. 4), including glutathione S-transferase (GST) 1, GST2, GSTσ3 and Prx-6, as well as 6-phosphogluconate dehydrogenase, which, as an enzyme of the pentose phosphate pathway (PPP), produces reducing equivalents in the form of NADPH to be used in repairing oxidatively damaged biomolecules. It is likely as well that aldehyde dehydrogenase (AldDH; two spots) and type 3 alcohol dehydrogenase, both with relatively high positive loadings on PC2, are involved in oxidative damage repair (see Discussion). In addition, three cytoskeletal proteins (adenylate cyclase-associated protein 1, tubulin α1C and tektin-4) are among the 18 identified proteins with the highest positive loading values on PC2 (Fig. 4). Finally, although three of 18 proteins with high positive loadings on PC2 are functionally associated with energy metabolism (Table 2) – triose phosphate isomerase (TPI), glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and arginine kinase (AK) – only GAPDH differs significantly in abundance among treatment groups (Fig. 4).

Proteins that increase most significantly in abundance in response to the 18E/6R treatment have high negative loadings on PC2 (Fig. 1A). Here, cytoskeletal proteins dominate, with seven of 17 identified proteins belonging to this functional category (Table 3, Fig. 5). These include two spots each identified as actin-interacting protein 1 and tektin-3, as well as tektin-4, β-tubulin and a tubulin-polymerization promoting protein. Energy metabolism and oxidative stress proteins are relatively less abundant, with two spots representing GAPDH identified, and one μ-class GST.

Negative loadings on PC3 (Fig. 1B) separate the control samples from the three treatment groups, and thus indicate those proteins that are relatively high in abundance in the control mussels compared with the exposed mussels. Here, 18 of 20 proteins were identified (Table 4), with BLAST searches of three spots returning the same hypothetical protein (homologous to a protease). Three of the identified proteins are associated with energy metabolism; however, only one, AK, is significantly increased in abundance (Fig. 6). Another, catalase, is an oxidative stress protein, and three are cytoskeletal or cell-adhesion proteins.

In summary, after G. demissa experience aerial exposure followed by immersion, changes in abundance of gill proteins associated with catabolism/fermentation, the cytoskeleton and oxidative stress predominate: 16, 11 and 16 of 66 identified proteins, respectively, fall into these categories across the three treatment groups. However, the duration of exposure and recovery clearly affects expression patterns: relatively short exposure followed by extended recovery (6E/18R) leads to increased abundance of catabolic and especially fermentation proteins; the 12E/12R treatment induces increases in the abundance of oxidative stress proteins more than the other treatments; and 18E/6R appears to impact cytoskeletal structure most significantly.

One of the hallmarks of bivalve anaerobic metabolism is the use of fermentative pathways that produce succinate and propionate as end products, rather than lactate as is more commonly found in vertebrates (Hochachka and Mustafa, 1972; Kluytmans et al., 1977; Ho and Zubkoff, 1982). Succinate/propionate fermentation appears to have evolved in eukaryotes such as helminth parasites (Barrett, 1991; Müller et al., 2012) and bivalves that undergo extended periods of hypoxia because it has the potential to produce more ATP – up to 6 ATP molecules per glucose molecule consumed, compared with 2 ATP molecules for lactate fermentation (Hochachka and Somero, 2002; Storey and Storey, 2005). Perhaps more importantly, in this type of fermentation, acidosis is relatively delayed, because lactate fermentation produces 1 mole ATP per mole H⁺ released, whereas propionate fermentation produces 3 moles ATP per mole H⁺ (Hochachka and Somero, 2002).

Work performed in the 1970s described the pathways involved in succinate/propionate fermentation in bivalves, using both blue mussels (genus Mytilus) and G. demissa to examine metabolic processes in the gill (Kluytmans et al., 1977; Malanga and Aiello, 1972), adductor mussel (Hochachka and Mustafa, 1972) or the entire organism (De Zwaan and Zandee, 1972; De Zwaan and van Marrewijk, 1973; Kluytmans et al., 1977; Ho and Zubkoff, 1982). The alternative aerobic and anaerobic pathways used by Mytilid mussels are shown in Fig. 3A. Phosphoenolpyruvate (PEP) in the cytosol represents an important branch point between aerobic and anaerobic metabolism, controlled by reversible inactivation of pyruvate kinase (Grieshaber et al., 1994; Hochachka and Somero, 2002; Müller et al., 2012). In G. demissa after 6 h of aerial exposure and 18 h of immersed recovery, we found that PEPCK was increased in abundance (Fig. 3A), suggesting the shunting of glycolytic metabolites into the succinate/propionate fermentative pathway. Furthermore, the next three enzymes in the pathway leading to succinate, namely cMDH, fumarase (in the mitochondrion) and sucDH/FRD (subunit B), were all higher in abundance in the 6E/18R group relative to control, indicating a concerted upregulation of this pathway. SucDH/FRD is the source of two further ATP molecules produced per glucose molecule, and the succinate generated can be decarboxylated to proprionyl CoA, gaining two further ATP molecules (Hochachka and Somero, 2002; Storey and Storey, 2005). We observed that PCC, the central enzyme in this second stage of the fermentative process, is also increased in abundance in G. demissa gill after 6E/18R treatment. In summary, we have found a clear signal – five proteins out of 18 identified – that fermentative capacity is increased in these mussels after a moderate (6 h) period of hypoxia induced by aerial exposure.

Remarkably, however, these five enzymes associated with anaerobic metabolism are not increased in abundance compared with control in mussels experiencing the longer duration aerial exposures, 12E/12R and 18E/6R (Fig. 3A). This result seems counter-intuitive, because G. demissa are well known for their ability to withstand long periods of aerial exposure (Kuenzler, 1961; Lent, 1969; Hilbish, 1987), and anaerobic metabolism combined with metabolic depression are considered the main tools with which bivalves survive hypoxia associated with aerial exposure (Gäde, 1983; Hochachka, 1986; Hochachka and Somero, 2002). However, as discussed below, this conundrum may be resolved by linking proteomic changes we found at later time points with behavioral responses of the mussel, specifically air gaping, after long periods of exposure.

A number of other energy metabolism proteins were identified among the 80 proteins having the highest component loadings across all four treatment groups. Two proteins associated with protein catabolism were also increased in abundance in the 6E/18R group (Fig. 3B): methylmalonate semialdehyde dehydrogenase is a component of the branched-chain amino acid degradation pathway, and glutamate dehydrogenase is central to the extraction of NH₃ from amino acids to provide carbon skeletons for catabolic processes. Both increase in abundance in the 6E/18R group in a pattern similar to that of the succinate/proprionate fermentation enzymes described above.

Notably, GAPDH is the only energy metabolism enzyme that showed significant abundance changes in the 12E/12R and 18E/6R groups (Figs 4, 5); this enzyme appears to have multiple roles in the cell beyond glycolysis, acting as a transcription factor (Zheng et al., 2003), a pro-apoptotic signal (Hara et al., 2005) and a sensor of oxidative stress (Chuang et al., 2005), and perhaps modifying cytoskeletal structure (Tisdale, 2002). Examination of representative gel images from each treatment group (supplementary material Fig. S1) reveals that the GAPDH identified in the 12E/12R and 18E/6R mussels occurs as a series of spots that have the same M_r but different isoelectric point (pI) values. Such a pattern is consistent with post-translational modifications (PTMs) to the enzyme that lead to different charge states, and thus horizontal positions on these gels, and indeed GAPDH is known to undergo a wide variety of PTMs (e.g. Tisdale, 2002; Sen et al., 2008; Seo et al., 2008). However, further research is necessary to describe the nature and function of these apparent PTMs to GAPDH.

In the control group, AK is the only identified energy metabolism protein to differ significantly in abundance relative to the exposure treatments (Fig. 6). AK plays an important role in recovery from anaerobiosis by generating the phosphagen phosphoarginine, which is able to rephosphorylate ADP during hypoxia without the generation of H⁺ (Ellington, 2001; Grieshaber et al., 1994). High loading values for AK in control as well as changes in the horizontal position of AK on treatment gels (supplementary material Fig. S1), which indicates a change of pI of AK after exposure relative to control, are consistent with PTM. Dawson and Storey have shown that crayfish AK is phosphorylated during hypoxia (Dawson and Storey, 2011), perhaps favoring phosphoarginine utilization over arginine phosphorylation. Analogously, modification of phosphoarginine levels via changes in AK activity may be a strategy G. demissa employs when preparing for extended bouts of aerial exposure.

A nearly inescapable correlate of hypoxia exposure is oxidative damage, due in part to reactive oxygen species (ROS) production as electrons leak from the slowed electron transport chain, as well as further ROS production upon re-oxygenation after severe hypoxia (Hermes-Lima and Zenteno-Savín, 2002; Welker et al., 2013). Here, we found a clear response to (or preparation for) oxidative damage in G. demissa as a consequence of prolonged aerial exposure. The strongest response, however, is not after the longest exposure period, 18E/6R, but after intermediate exposure, 12E/12R. Proteins with the highest loadings in this treatment group include five that are clearly associated with oxidative stress – three isoforms of GST (GST-1, GST-2 and GST-3), a peroxiredoxin (Prx-6) and the PPP enzyme 6-phosphogluconate dehydrogenase (Fig. 4). GSTs comprise a large and diverse superfamily of enzymes whose functions include the transfer of glutathione (GSH) groups to initiate repair of oxidatively damaged biomolecules (Hayes et al., 2005). Notably, in a similar study by Almeida et al. examining the effects of aerial exposure followed by re-immersion in the brown mussel Perna perna (Almeida et al., 2005), the authors found that of the six antioxidant enzymes assayed – catalase, superoxide dismutase, glutathione peroxidase, glutathione reductase, glucose 6-phosphate dehydrogenase and GST – only GST increased in activity after aerial exposure. Prx-6 is expressed in the cytosol and uses GSH, rather than thioredoxin like other Prx isoforms, to reduce and detoxify hydrogen peroxide and lipid peroxides (Schremer et al., 2007; Fisher, 2011). Interestingly, this GSH is transferred to the Prx-6 active site cysteinyl residue by GSTs (Zhou et al., 2013). Finally, 6-phosphogluconate dehydrogenase is also linked to the glutathione system, because as a component of the PPP it is central to the production of NADPH in the cytosol, which is used to regenerate GSH from oxidized glutathione.

Notably, the enzymes AldDH and type-3 alcohol dehydrogenase (ADH-3) also are both increased in abundance at 12E/12R relative to other treatments (Fig. 4). [Aldehyde dehydrogenase was identified twice, although examination of 12E/12R gels (supplementary material Fig. S1) indicates that these identifications may represent a single protein spot.] Aldehyde dehydrogenase acts to reduce aldehyde groups formed from sources such as lipid peroxides (Kotchoni et al., 2006) during oxidative stress, and is an enzyme of the ‘minimal stress proteome’ (Kültz, 2005). In addition, it has been shown to change in abundance in response to acute heat and xenobiotic stress in Mytilid mussels (Tomanek and Zuzow, 2010; Apraiz et al., 2006). Comparably, the ADH-3 identified in the 12E/12R group uses the reducing power of glutathione to eliminate formaldehyde (Danielsson et al., 1994) and S-nitrosothiols (Liu et al., 2001). It is worth noting that ADH-3 is the only ADH form found in cephalopod mollusks (Kaiser et al., 1993). These results suggest a role for AldDH and ADH-3 in protection from oxidative damage, and, assuming this evaluation is correct, eight of 18 proteins with the highest component loadings in 12E/12R are associated with protection from or repair of oxidative damage.

In contrast, the 6E/18R mussels only had one protein clearly associated with oxidative stress among the 18 identified, mitochondrial Prx-5, as well as a second protein, transketolase, that may be associated with oxidative stress because of its role in the non-oxidative phase of the PPP (Fig. 3B). The 18E/6R group only included one oxidative stress protein, GST μ, a cytoplasmic isoform, despite experiencing the longest period of exposure (Fig. 5). These results support the interpretation that the shortest duration of exposure, 6E/18R, is not sufficient to lead to oxidative damage requiring a significant increase in the activity of oxidative stress proteins, whereas the 12E/12R treatment is. In our model, the absence of a robust oxidative stress response in the 18E/6R group is associated with their shift of strategy from anaerobiosis to air gaping.

Interestingly, one important oxidative stress protein, catalase, appears to decrease significantly in abundance in all treatment groups relative to control (Fig. 6). Catalase is a ubiquitous enzyme central to the detoxification of hydrogen peroxide, so an apparent reduction in abundance during oxidative stress is counter-intuitive, although there is some evidence for such a response in human mesangial cells (Venkatesan et al., 2007).

As mentioned above, an unexpected finding of this study is that the mussels experiencing the longest period of aerial exposure, 18E/6R, evinced no signal of increased fermentative capacity and only minimal upregulation of oxidative stress protection and repair. We hypothesize that these mussels have a weak proteomic response to hypoxia because they have altered their strategy of coping with aerial exposure from one of anaerobiosis, which appears robustly at 6E/18R, to one of aerobic metabolism through exposure of gills to air.

Air-gaping behavior in G. demissa has been described repeatedly (Lent, 1968; Lent, 1969; Malanga and Aiello, 1972; Widdows et al., 1979), and it has been shown that air-gaping G. demissa are able to undergo aerobic metabolism (Widdows et al., 1979; Nicchitta and Ellington, 1983). Indeed, air gaping is not limited to G. demissa, and has been described in a number of other high-intertidal mussels (Nicastro et al., 2010; Dowd and Somero, 2013) as a response to both hypoxia and high temperature (via evaporative cooling) (Lent, 1968). The present experiment was not designed to examine the impact of gaping on protein expression, and so we measured neither the time of onset nor the duration of gaping in mussels of different treatment groups. However, we were able to determine straightforwardly whether gaping had occurred, because during aerial exposure the mussels were positioned partially buried in artificial seawater (ASW)-saturated sand. When they had opened their valves, the water trapped within was clearly apparent, pooled above the sand in the bottom of the air-filled jar. All of the mussels in the 18R/6E group were gaping at the end of the experiment, with trapped water pooled in the jar; in the 6E/18R group there was no water in the jars, indicating the mussels never gaped during the exposure. The 12E/12R group was intermediate – water present in the jars indicated that mussels had gaped, but not all mussels were gaping at the end of the exposure; thus, we cannot quantify the extent of the behavior in this group.

Although gaping allows G. demissa to resume aerobic metabolism, there are trade-offs to the air-gaping strategy. When the mussel's valves are separated, the water held within is lost, which clearly leads to the risk of desiccation, and early work showed that G. demissa is remarkably tolerant to tissue water loss (Lent, 1968). In the present experiment, desiccation stress was minimized by ensuring air in the jars was water saturated. A second stress that may occur during gaping is a loss of physical support for the delicate gill tissue. When the valves are closed and filled with water, the gills are effectively weightless; however, when the water is lost, the gills collapse and adhere. We hypothesize that the significant changes in expression of gill cytoskeletal proteins in the 18E/6R group are a response to physical damage that occurs as the gills lose structural support when exposed to air, or a preparation for future air exposure during the next tidal cycle. In this group at least seven of 17 identified proteins with high component loadings are associated with the cytoskeleton (Fig. 5). Notably, three are tektin isoforms (two tektin-3 and one tektin-4); tektins are integral microtubule proteins of cilia and flagella (Pirner and Linck, 1994; Amos, 2008). The presence of tektins is noteworthy because the numerous cilia of mussel gills are responsible for creating the water currents that allow the mussel to feed and respire. Any damage to these cilia, which are directly exposed to the external medium, would require rapid repair to maintain feeding and metabolic processes. In addition to the tektins, β-tubulin increased in abundance, as did tubulin-polymerization protein family member 3, suggesting repair or reorganization of microtubule assemblies, perhaps including those related to ciliary function. The actin cytoskeleton also appears to be modified in 18E/6R mussels, as actin-interacting protein 1 also increased in abundance; this protein works together with actin and cofilin to help control actin filament disassembly (Ono, 2003), thus affecting cell morphology and migration, cell division and endocytosis (Lin et al., 2010).

In contrast, the 6E/18R and 12E/12R groups both showed weaker cytoskeletal responses, with only three of the identified proteins in each group associating with this function (Fig. 3B and Fig. 4, respectively). Again, changes to the tubulin cytoskeleton – tubulin α-1c in 6E/18R and tubulin α-1c and tektin-4 in 12E/12R – may indicate restructuring of cilia. There were three cytoskeleton-associated proteins with high component loadings in the control group; that is, occurring with relatively low abundance in all treatment groups relative to control. These included a tektin isoform, an actin, and talin-2, which is a major component of focal adhesion plates linking the internal actin cytoskeleton to integrin and altering its affinity for extracellular matrix proteins (Calderwood, 2004; Critchley and Gingras, 2008).

To summarize, we argue here that the changes in cytoskeletal protein abundance in the 18E/6R group are a result of structural damage due to air gaping. However, we recognize that further work examining changes to gill structure during gaping in G. demissa and other intertidal mussels is needed to provide support for our hypothesis that air gaping results in a trade-off between the benefits of aerobic metabolism and the costs of damage to gill tissue.

Molecular chaperones help stabilize or refold other proteins during denaturing stresses, and are central to the minimal cellular stress response (Kültz, 2005). Thus, it is noteworthy that chaperones do not represent a larger proportion of the proteins whose abundances change significantly in response to aerial exposure. Chaperones with high loading values according to PCA include calreticulin (12E/12R; Table 2), an endoplasmic reticulum chaperone, heat shock protein (HSP) 70 (18E/6R; Table 3), a ubiquitous cytosolic chaperone, and HSP24.1 (control, Table 4), one of the family of small HSPs that interacts with and stabilizes the actin cytoskeleton during stress (Arrigo, 1998; Dalle-Donne et al., 2001). However, these last two chaperones were not found to differ significantly in abundance by ANOVA. The small number of chaperones with high loading values, the lack of statistical significance in their abundance changes and the lack of a clear pattern among the treatment groups suggest that protein denaturation is not a significant consequence of the aerial exposures experienced by these mussels, although greater desiccation stress during air gaping might lead to greater chaperone expression.

In conclusion, our results demonstrate that G. demissa, and perhaps other intertidal bivalves, do not rely on a single strategy through time to cope with the stress of aerial exposure. Instead, the evidence we present here indicates that exposures of increasing duration lead to the initiation of a variety of responses. First, fermentative capacity is upregulated, but prolonged hypoxia followed by re-oxygenation may lead to oxidative stress, so longer term exposure is marked by increases in oxidative stress proteins. Finally, after the longest duration exposures, fermentation is abandoned as a survival strategy, and air gaping is employed to allow aerobic metabolism to resume, with concomitant physical damage to gill tissue. Further research is necessary to describe how correlated stresses, including desiccation, UV exposure and high or low temperature may affect the biochemical and behavioral survival strategies hypothesized here.

Adult G. demissa (shell length 7–10 cm) were collected from the salt marsh south of Stone Harbor, NJ, USA (39.04N, 74.77W) in May 2012. Mussels were chilled and transported to Franklin and Marshall College, Lancaster, PA, USA, where they were kept immersed in ASW (33–35 ppt; Instant Ocean, Marineland Labs, Mentor, OH, USA) in recirculating aquaria. Experimental animals were allowed to acclimate for 4 weeks constantly immersed at 15±0.2°C, and were fed, per mussel, 0.1 ml Instant Algae shellfish diet 1800 (Reed Mariculture, Campbell, CA, USA) three times per week.

Experimental mussels were divided into four groups of six individuals. Control mussels experienced no emersion before they were killed. Treatment mussels were exposed to water-saturated air in 800 ml jars in groups of three; within the jars they were buried approximately ¼ of their shell length in an ASW–sand slurry, in the same vertical position they maintain in the salt marsh. Jars were loosely covered and held at 15°C in a temperature-controlled chamber (±0.2°C) for the duration of exposure, and mussels were returned to the 15°C recirculating aquaria for the indicated recovery times before they were killed. No mortality occurred during treatment and recovery. The three treatment groups experienced 6 h of aerial exposure followed by 18 h of immersed recovery (6E/18R), 12 h exposure and 12 h recovery (12E/12R) or 18 h exposure and 6 h recovery (18E/6R).

Protein extraction, purification and separation methods followed those described elsewhere (Tomanek and Zuzow, 2010). At the end of the recovery period, gill tissue was excised and immediately homogenized (1:4 w/v) in denaturing homogenization buffer [7 mol l⁻¹ urea (Bio-Rad, Hercules, CA, USA), 2 mol l⁻¹ thiourea, 1% ASB-14 (Calbiochem, San Diego, CA, USA), 40 mmol l⁻¹ Tris base, 0.5% ampholyte solution (pH 3–10, GE Healthcare, Piscataway, NJ, USA), 1.2% Destreak Solution (GE Healthcare) and 0.001% Bromophenol Blue] using ground glass homogenizers. (Unless otherwise specified, chemicals were purchased from Sigma Chemical Co., St Louis, MO, USA.) Homogenates sat at room temperature (RT) for 1 h to solubilize proteins and were centrifuged (16,000 g, 30 min, RT) to pellet insoluble material.

Proteins were precipitated from homogenate supernatants using 10% trichloroacetic acid in acetone (1:4 v/v supernatant: TCA/acetone). After overnight incubation (−20°C), precipitated proteins were pelleted by centrifugation (16,000 g, 30 min, 4°C), supernatant was removed, and pellets were washed three times by vortex-mixing with ice-cold acetone. Protein was resuspended using a volume of rehydration buffer (7 mol l⁻¹ urea, 2 mol l⁻¹ thiourea, 2% w/v CHAPS, 2% NP40 substitute, 0.6% pH 3–10 IPG buffer, 1.2% Destreak Solution and 0.01% Bromophenol Blue) equal to the original homogenization buffer volume. Protein was vortex-mixed and resolubilized overnight at 4°C. Samples were centrifuged (16,000 g, 30 min, RT) and the protein concentration of each supernatant was determined using the 2D Quant Kit (GE Healthcare). Supernatants were used immediately for electrophoresis or stored at −80°C.

Proteins were separated on polyacrylamide gels by pI in the first dimension and by molecular mass in the second (Görg, 2004). Protein samples (200 μl of 1 mg ml⁻¹) were loaded into a single well of an 11 cm isoelectric focusing (IEF) tray (Bio-Rad IEF Cell). An immobilized pH gradient (IPG) IEF strip (11 cm, pH 3–10; Bio-Rad) was laid over each sample, and covered with mineral oil to prevent evaporation. Strips rehydrated passively for 6 h and actively (50 V) for 12 h; water-saturated wicks were placed between the electrodes and the strips before focusing [12,000 volt hours (V h); maximum current 50 μA]; wicks were replaced and focusing continued for a total of 24,000 V h. Oil was then blotted from the strips and they were either frozen (−80°C) or used immediately for second-dimension electrophoresis.

In the second dimension, proteins were separated by mass using precast SDS-polyacrylamide gradient gels (8–16%; Bio-Rad Criterion). Proteins in the IPG strips were equilibrated with SDS by agitation in equilibration buffer [6 mol l⁻¹ urea, 375 mmol l⁻¹ Tris-HCl pH 8.8, 30% glycerol, 2% SDS (Bio-Rad) and 0.01% Bromophenol Blue]. Cysteinyl residues were covalently capped with carbamidomethyl groups to prohibit formation of disulfide bridges, by adding 10 mg ml⁻¹ dithiothreitol (DTT) (Bio-Rad) to the buffer used for the first 15 min equilibration, followed by 15 min in equilibration buffer without DTT but containing 25 mg ml⁻¹ iodoacetamide (Bio-Rad). Strips were placed atop SDS-PAGE gels and held in place with 0.8% agarose. Gels were run simultaneously in groups of six or 12 in Tris-glycine running buffer (25 mmol l⁻¹ Tris base, 190 mmol l⁻¹ glycine, 0.1% SDS) using the Bio-Rad Dodeca Cell (200 V, ~1 h) cooled to 10°C. Gels were stained overnight with colloidal Coomassie Blue, and destained with multiple washes of 18 MΩ water.

Gel images were digitized at 600 dpi resolution using an Epson model 700 V transilluminating scanner. Images were preprocessed using the histogram tool within the Epson software to maximize contrast, setting γ=1 to maintain linearity in pixel density. Images were cropped to the gel boundaries using Photoshop software (Adobe Systems). Twenty-four 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 (Berth et al., 2007). The fusion image was used to detect spots employing default software settings, and spot boundaries were transferred back to the original 24 gel images. The 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 in response to each treatment.

Before analysis, we applied a low-volume cut-off filter to remove spots that had a mean relative spot volume less than 0.05% in every group, because protein identification via tandem mass spectrometry (MS/MS; described below) often is unsuccessful when protein volume is below this threshold. We employed PCA within the MeV statistics suite of Delta 2D to determine which protein spots contributed most to the separation of treatment groups, by ranking according to component loading along the first three component axes. We chose the 20 spots from each treatment group (i.e. 80 spots in total) with the highest loadings from the PCA analysis for subsequent identification via MS/MS. We also performed hierarchical cluster analysis (HCL) (Pearson correlation; average linking) on these 80 spots to better understand relationships in PEPs among spots.

We employed one-way ANOVA to determine which of the 80 proteins changed significantly in abundance in response to exposure treatments, utilizing a null distribution (5000 permutations) to account for potential unequal variances and non-normal distributions of protein abundance (Tomanek and Zuzow, 2010). In order to limit Type 1 errors we used a P<0.02 rather than a multiple comparison correction (Tomanek and Zuzow, 2010; Fields et al., 2012). We applied a Tukey's HSD post hoc test (IBM SPSS Statistics, version 20) to identify treatments with differing abundance for each protein shown to be significant via one-way ANOVA.

The 80 spots identified as contributing most significantly to differences in PEPs among the treatment groups were excised from gels and destained with 25 mmol l⁻¹ ammonium bicarbonate in 1:1 (v/v) acetonitrile:water. Gel pieces were dehydrated by soaking in 100% acetonitrile, and were re-swelled in digestion buffer (25 mmol l⁻¹ ammonium bicarbonate) containing 10 ng μl⁻¹ modified porcine trypsin (Promega). After 4 h incubation (37°C) with shaking, the digest solution was collected for MS analysis.

Tryptic peptides were separated and analyzed using a high-pressure liquid chromatography–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 non-polar C8 column (Agilent) before injection into the trap. Agilent software detected peptides via MS, which were automatically analyzed via MS/MS to allow determination of amino acid sequences. 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 constructed using RNA isolated from gill tissue of four G. demissa individuals (MPG mRNA purification kit, PureBiotech, Middlesex, NJ, USA) collected from southern (Darien, GA, USA) and northern (Nummy Island, NJ, USA) populations. cDNA was synthesized at the Pennsylvania State University Genomics Core Facility and sequenced using a Roche 454 FLX system. ESTs were assembled de novo using iAssembler (Zheng et al., 2011). After assembly, the library contained 47,736 unique entries that we used 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.

We thank Dr I. Albert of the Bioinformatics Consulting Center at Pennsylvania State University for valuable help in assembling the G. demissa gill EST library. We also express our gratitude to Dr K. Hess for his help in the mass spectrometry portion of this study. Our discussions with Dr L. Tomanek on analysis of proteomic data were invaluable.