Transcriptomic responses to heat stress in invasive and native blue mussels (genus Mytilus): molecular correlates of invasive success

Invasive species pose numerous threats to ecosystems, including displacement of native species (Mooney and Cleland, 2001; Sax et al., 2005). Much work has focused on the anthropogenic causes of species introductions (Carlton and Geller, 1993) and the species-level interactions that contribute to the spread of invasive species (Holway, 1999; Huxel, 1999; Jensen et al., 2002). However, little is currently known about the role played by physiological factors in facilitating invasions. For example, an invader may displace native species if it is better adapted to the abiotic conditions in the invaded environment than native competitors (Blackburn and Duncan, 2001). With increasing temperatures due to climate change, it is conceivable that more warm-adapted species may be able to enter new habitats and replace similar species that have lower tolerance of high temperature (Somero, 2010). In fact, it has been shown that global warming generally will favor invasions in marine habitats (Stachowicz et al., 2002).

To further analyze the role that history of thermal adaptation might play in conferring a competitive advantage to a warm-adapted invasive species, we initiated a broad comparative study of two congeners of blue mussels: Mytilus trossulus, a species native to the West Coast of North America, and Mytilus galloprovincialis, a species that originated in the Mediterranean Sea and which has successfully invaded coastal habitats at several sites in Africa, Asia and North and South America. Along the coastline of California, M. galloprovincialis has replaced the native congener at sites south of Monterey Bay. A mosaic hybrid zone containing both congeners and their hybrids is found in the San Francisco and Monterey Bays, and the native M. trossulus remains dominant north of Bodega Bay (Suchanek et al., 1997; Rawson et al., 1999; Braby and Somero, 2006a). These two congeners are closely related, having diverged approximately 3.5 million years ago when M. trossulus populations in the Northern Pacific and Northern Atlantic were separated due to an expansion of the polar ice cap (Vermeij, 1991; Seed, 1992). This separation gave rise to a second species of blue mussel, Mytilus edulis, which spread to European waters from North America. Subsequently, populations of M. edulis in the Mediterranean Sea became isolated when the Strait of Gibraltar closed approximately 2 million years ago, giving rise to a third species of blue mussel, M. galloprovincialis (Seed, 1992). This species has evolved in an environment with higher and more stable temperatures than the other two blue mussels, and a number of physiological, biochemical and molecular studies have shown it to be more warm-adapted than its north temperate congener Mytilus trossulus. For example, heart function in M. galloprovincialis can be sustained at higher temperatures than in the native species whereas M. trossulus can maintain heart function down to lower temperatures (Braby and Somero, 2006b). The kinetic properties of the orthologs of cytosolic malate dehydrogenase in the two congeners are consistent with M. galloprovincialis being a more warm-adapted species (Fields et al., 2006). Lastly, induction of heat shock proteins (Hsps) has been shown to be more pronounced in M. galloprovincialis (Hofmann and Somero, 1996), suggesting a more robust ability of this species to initiate transcription of genes encoding heat-shock proteins.

Because the ability to respond physiologically to acute heat stress is based in large measure on capacities for broad-scale changes in gene expression (i.e. in the transcriptome) (Podrabsky and Somero, 2004; Gracey et al., 2004; Buckley et al., 2006), we sought to examine how the transcriptomes of M. trossulus and M. galloprovincialis were affected by acute increases in temperature. To this end, we developed a novel oligonucleotide microarray that contained 11,506 genes based on sequences available for two congeners of Mytilus, M. galloprovincialis and the ribbed mussel, Mytilus californianus. Mytilus californianus is an out-group species that is approximately 7.6 million years diverged from M. galloprovincialis and M. trossulus (Ort and Pogson, 2007). Hybridization to the array thus involved both homologous and heterologous hybridizations, leaving 8874 of the probes, representing 4488 different genes, with sufficient binding to be studied. Mytilus galloprovincialis and M. trossulus exhibited strikingly similar transcriptional responses to acute heat stress, particularly among molecular chaperone genes; however, they showed significantly different responses among a subset of genes involved in oxidative stress, proteolysis, energy metabolism, ion transport, cell signaling and cytoskeletal reorganization. These differences may be partially responsible for the physiological differences in heat tolerance noted in earlier studies and may play a role in determining the future success of the invasive M. galloprovincialis in warming marine habitats.

We collected mussels from floating boat docks: Mytilus trossulus (Gould 1850) from Newport Harbor, Newport, OR, USA (44°38′25″N, 124°03′10″W) and Mytilus galloprovincialis (Lamarck 1819) from Santa Barbara Harbor, Santa Barbara, CA, USA (34°24′15″N, 119°41′30″W). All specimens thus were acclimatized to conditions of continuous immersion. Care was taken to collect animals of similar size (5–7 cm) to avoid potential artifacts from size-related influences on gene expression. Specimens were transported to the laboratory of Dr Lars Tomanek at California Polytechnic State University, San Luis Obispo, and placed in a recirculating seawater tank in which they were common-garden acclimated at 13°C for 4 weeks, in order to minimize the effects of past thermal history. We chose both collection sites to be well outside of the hybrid zone (McDonald and Koehn, 1988; Braby and Somero, 2006a) to ensure no hybrid individuals were included in our study. We confirmed the species identification by genotyping each mussel using published methods (Heath et al., 1995; Rawson et al., 1996). The 13°C acclimation temperature mimicked an ecologically realistic thermal environment for both species (Braby and Somero, 2006a).

All mussels remained submerged in seawater throughout the heat-stress experiment. At the start of the experiment, six mussels of each species from the 13°C common garden were killed and their gill tissues dissected and frozen on liquid nitrogen. Treatment animals were then transferred from the common-garden aquarium to a tank in which they were exposed to a heat ramp of +6°C h^–1, from 13°C to 32°C, to simulate thermal heating on a hot day during falling tide (Tomanek and Somero, 1999). Control animals were simultaneously transferred to a separate tank that remained at 13°C. At designated temperature points on the heat ramp (24°C, 28°C and 32°C), six treatment mussels of each species were transferred to a separate tank and maintained at the respective temperature for an additional hour. Treatment and control mussels corresponding to each temperature point were then killed and their gill tissues dissected and frozen on liquid nitrogen. We used gill tissues for this study because (1) stress-induced protein expression has been shown to be responsive to environmental stress in Mytilus gill (Hofmann and Somero, 1995; Hofmann and Somero, 1996), (2) Mytilus californianus has been shown to have highly similar stress-induced transcriptomes between gill, mantle and adductor muscle (Gracey et al., 2008), indicating that gill can be used as a representative tissue, and (3) the microarray was designed from expressed sequence tag (EST) sequences generated from the gill tissue of the out-group species, M. californianus (see Microarray development section below). Tissues were stored at –70°C prior to RNA extraction. For the companion proteomic study (Tomanek and Zuzow, 2010), mussels from the experiment described above were returned to the 13°C common-garden aquarium for a recovery period of 24 h before being killed.

Microarrays were designed using approximately 26,000 publicly available ESTs from M. californianus [ES7325872–ES738966; ES408175–ES387463 (Gracey et al., 2008)] and M. galloprovincialis [AJ623313–AJ626468; AJ516092–AJ516921 (Venier et al., 2003)]. Both sets of ESTs were sequenced from libraries constructed using Suppression Subtractive Hybridization (SSH) (Diatchenko et al., 1996) to enrich for transcripts expressed in response to various stressors.

Sequences for each species were separately screened for vector contamination, clustered, assembled and loaded into a PostgresQL database using the PartiGene bioinformatics pipeline (Parkinson et al., 2004), yielding a total of 12,961 and 1688 putative gene transcripts (clusters) for M. californianus and M. galloprovincialis sequences, respectively. Clusters were annotated for identity against the UniRef 90 protein database using NCBI BLAST (basic local alignment search tool), with an e-value cut-off of 10^–8. Clusters were also annotated for Gene Ontology (GO) and enzyme commission number (EC) using the program annot8r (Schmid and Blaxter, 2008), again with an e-value cut-off of 10^–8.

Long oligo (60-mer) probes were designed against the M. californianus clusters using the probe design program YODA (Nordberg, 2005) under default parameters, with a maximum of 5 probes per cluster. This yielded 43,969 total unique probes to 11,506 clusters, with a mean of 2.6 probes per cluster. To measure the effects of interspecific sequence variation on probe performance, the above probes were compared with M. galloprovincialis clusters using BLAST. Probes designed to the M. californianus library that successfully aligned with the smaller M. galloprovincialis library had counterparts designed using the M. galloprovincialis sequence. This resulted in 556 pairs of probes that had matched M. californianus and M. galloprovincialis sequences, with a mean number of 4.6 divergent nucleotide bases per probe. Altogether, there were 44,524 unique probes that were duplicated or triplicated randomly to fill a 105,000-feature microarray that was in situ synthesized by Agilent Technologies, Inc. (Santa Clara, CA, USA).

Total RNA was extracted from gill tissues with a TissueLyser homogenizer (Qiagen, Valencia, CA, USA) and RNeasy kit (Qiagen). Following RNA extraction, samples were treated with DNase using the Ambion TURBO DNA-free kit (Applied Biosystems, Foster City, CA, USA). The DNase was removed by incubation with resin beads (included in the kit), which bind DNase, followed by application to Ambion QuikSpin columns (Applied Biosystems) to filter out the DNase-bound beads. Total RNA quality was assessed on 1% agarose gels. RNA quality was also assessed and concentration determined for this and all subsequent steps with a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). 1 μg of total RNA was reverse transcribed and amplified with the Ambion AminoAllyl MessageAmp II kit (Applied Biosystems). To reduce problems in transcript binding due to sequence differences between species, we employed a reference design that used separate pools of reference RNA for each species. Reference amplified RNA (aRNA) was created for each species by pooling RNA before and after amplification. The reference pool was made up of RNA from seven different samples: one control sample from the beginning of the experiment, one treatment sample from each temperature point, and one control sample from each temperature point. 5 μg of aRNAs from reference and experimental samples were separately labeled with Alexa Fluor fluorescent dyes (Invitrogen, Carlsbad, CA, USA), 555 nm for reference and 647 nm for control or treatment, and dye incorporation was assessed with a NanoDrop spectrophotometer (NanoDrop Technologies). 750 ng of labeled aRNAs from reference and experimental samples were competitively hybridized to the microarray. We followed the recommended protocols from Agilent Technologies, Inc. for all hybridization and washing steps. Microarrays were scanned with an AXON GenePix 4000B scanner (Axon Instruments, Molecular Devices, Sunnyvale, CA, USA).

Fluorescent signal intensities for each spot on the microarray were extracted and Lowess normalized within each microarray using Feature Extraction Software (version 9.5.3.1; Agilent Technologies, Inc.). Only spots that Feature Extraction Software measured to have good hybridization among all samples were included in further analysis. Because each gene was represented by a set of one to five unique probe sequences, we summarized the expression level of each gene in each sample by computing the geometric mean of the normalized signal intensities for each probe set on a single microarray. We then computed the log₂-ratios of the experimental channel divided by the reference channel and normalized these expression values to the median log₂-ratio of the control samples of the same species from the same temperature point on the heat ramp (i.e. 13°C, 24°C, 28°C and 32°C). Log₂-ratios of the samples from the 13°C temperature point were normalized to the median log₂-ratio within that group, separately for each species.

We were interested in genes that showed a significant change in expression in response to heat stress and genes that showed different responses to heat stress between the two species. Therefore, we conducted a two-way analysis of variance (ANOVA), in which temperature and species were modeled as fixed effects, and focused on genes that were significant for the temperature effect or the species × temperature interaction. We ignored the species term from the ANOVA as this effect could have highlighted genes that differentially bound to probes on the microarray due to differences in sequence homology, thus not reflecting true differences in gene expression. In accordance with statistical convention, all genes with a significant species × temperature interaction were deemed not to have significant temperature effects, even if the temperature term from the ANOVA had a low P-value. We also conducted a separate one-way ANOVA to compare the transcriptomes of control samples from different time points in the heat ramp. Hereby we were able to determine if there was significant temporal variation in gene expression that was not due to temperature treatment. All P-values were adjusted with a false discovery rate (FDR) correction for multiple testing by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). All genes with FDR-corrected P-values less than 0.01 were considered significant. Analyses were conducted in the R statistical programming environment (R Development Core Team, 2009). Microarray data were submitted to the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/; GEO Accession: GSE19031).

We used principal component analysis (PCA), to illustrate the relative contribution of experimental variables (i.e. temperature and species) to the variation in gene expression. PCA identifies axes of maximal variance in the data, defined as principal components (PCs), revealing explanatory variables in large data sets. We evaluated the results of the PCA using loading plots; the spatial clustering of experimental samples on loading plots reflects the degree of similarity between transcriptomes. We used the normalized log₂-ratio data of all genes with good hybridization among all samples (see above). The PCA was computed using the correlation matrix, and analyses were conducted in the R statistical programming environment (R Development Core Team, 2009).

To describe the effect of temperature on biological function, we conducted a functional gene group correlation analysis (GCA) in the ErmineJ software package (Lee et al., 2005). We used the log₂-ratio data of all genes found to be significant (ANOVA, P<0.01) for the temperature effect and defined functional groups based on GO biological process (Ashburner et al., 2000). GCA tests whether genes in a particular functional category, or gene set, share a similar expression pattern. We used the following parameters in the analysis: (1) a minimal gene set size of 5 genes, (2) a maximal gene set size of 100 genes, and (3) the model was run for 10,000 iterations. ErmineJ uses the Benjamini–Hochberg method (Benjamini and Hochberg, 1995) for FDR correction, and we considered gene sets with FDR-corrected P-values less than 0.1 to be significant, following the authors' recommendations (Lee et al., 2005).

To compare the expression patterns of genes with species-specific responses to temperature, we used the log₂-ratio data from genes found to have significant species × temperature interactions and conducted a hierarchical clustering analysis. Hierarchical clustering groups genes based on expression pattern across experimental conditions, thereby identifying clusters of genes with common expression profiles. We performed the hierarchical clustering algorithm using the Spearman correlation dissimilarity matrix and the Ward agglomerative linkage method. Analyses were conducted using the gplots package (Warnes et al., 2009) in the R statistical programming environment (R Development Core Team, 2009).

Clusters defined in the hierarchical clustering analysis were functionally characterized by a functional over-representation analysis (ORA) in ErmineJ (Lee et al., 2005). ORA tests whether a particular functional category, or gene set, is significantly enriched among genes in a particular cluster versus the gene sets represented among genes in all other clusters. We used a minimal gene set size of 2 genes and a maximal gene set size of 100 genes. As in the GCA method described above, ORA uses the Benjamini–Hochberg method (Benjamini and Hochberg, 1995) for FDR correction, and we considered gene sets with FDR-corrected P-values less than 0.1 to be significant. In contrast to the GCA, for the ORA we defined gene sets based on hypothesized gene functions that we formulated using information from the GO database (Ashburner et al., 2000), the UniProt Knowledge Database (The UniProt Consortium, 2010) and the relevant published literature. These multiple sources of functional information provided a more detailed annotation than the GO database alone. Based on these criteria, we defined the following 22 functional categories: alternative splicing, amino acid metabolism, anti-apoptosis, cell cycle arrest, cell signaling, cytoskeletal reorganization, energy homeostasis, energy metabolism, immune response, ion transport, lipid metabolism, membrane stabilization, molecular chaperone, oxidative stress, post-translational modification, protein targeting, proteolysis, RNA splicing, transcription, transcriptional regulation, vesicle formation and unknown. This process of annotation was only feasible for the relatively small number of genes in this subset but it provided a more complete functional characterization of the genes that showed species-specific responses to heat-stress (see supplementary material Table S1 for gene annotations and references).

Heterologous hybridization to the microarray rendered 8874 probes that represented 4488 genes whose spots reliably hybridized in all samples, and thus performed well for both species. Only data from these 4488 genes were used in the transcriptomic analyses. 556 pairs of probes had matched sequences between M. californianus and M. galloprovincialis ESTs. Among these matched pairs of probes, we found no significant effect of sequence divergence on log₂-ratio of M. galloprovincialis or M. trossulus samples (data not shown).

Mytilus galloprovincialis and M. trossulus showed large shifts in gene expression in response to the heat-ramp treatment. Strikingly, the two species showed overall similar transcriptional responses (Fig. 1A), and 1531 genes had significant temperature effects (34% of all gene assayed; ANOVA, FDR-corrected P<0.01 for temperature effect). We found no evidence of temporal variation in gene expression, as there were no genes significantly differentially expressed between control groups from different time points (ANOVA, FDR-corrected P>0.99).

The PCA illustrated that most of the variance in gene expression was due to heat-shock treatment and not due to species differences. Principal component 1 (PC1), which explained 20% of the variance in expression, clustered the 13°C control samples separately from the 24°C, 28°C and 32°C heat-shocked samples but showed no separation between the species (Fig. 1B). Principal components 2, 3 and 4 explained 10%, 5% and 4% of the variance, respectively, and showed no clear clustering of any particular experimental group (data not shown). Interestingly, principal component 5 (PC5), which explained 3% of the variance in expression, clustered the two species separately for the heat-shocked samples (24°C, 28°C and 32°C; Fig. 1B).

The GCA revealed that among the 1531 genes with significant temperature effects, there were two GO biological processes whose genes showed consistent responses to heat stress: protein folding (GO:0006457; Fig. 2A; GCA, FDR-corrected P=0.035) and proteolysis (GO:0006508; Fig. 2B; GCA, FDR-corrected P=0.024). The protein folding biological process consisted of 22 genes that were mostly up-regulated at all heat-shock temperatures (24°C, 28°C and 32°C; Fig. 2A). Notably, in this group there were two genes of DnaJ (an ortholog of Hsp40), two genes of Hsp90, Hsp60, and eight genes encoding subunits of chaperonin TCP1. The proteolysis biological process consisted of 31 genes, most of which showed no change in expression or were down-regulated at 24°C and 28°C and then were up-regulated at 32°C (Fig. 2B). Many of the genes in the proteolysis category were specific to ubiquitin-mediated proteolysis, such as ubiquitin carboxyl-terminal hydrolase 5, five genes encoding subunits of the proteasome, ubiquitin fusion degradation 1-like protein, ubiquitin-conjugating enzyme E2 and ubiquitin specific proteases 14 and 52. Another gene of note in this category was leucine aminopeptidase (Lap), which was also down-regulated at 24°C and 28°C and up-regulated at 32°C (Fig. 2B). Interestingly, despite the overall similar patterns of expression between the two species, among genes in the proteolysis category M. trossulus showed an average of 1.3-fold higher up-regulation than M. galloprovincialis at 32°C (Fig. 2B).

The largest changes in expression in response to heat-stress were among genes that encoded the molecular chaperone heat shock protein 70 (Hsp70). In total, there were seven Hsp70 genes that were significantly differentially expressed, all of which were strongly up-regulated at all heat-shock temperatures (24°C, 28°C and 32°C) in M. galloprovincialis and M. trossulus. One Hsp70 gene showed an up-regulation of more than 256-fold. As can be seen in Fig. 3, which shows the mean log₂-ratio among all seven Hsp70 genes, there was a strikingly similar response between the species and the expression levels were constant across the three heat-shock temperatures. The mean ± s.e.m. log₂-ratio among the seven Hsp70 genes, across the three heat-shock temperatures, was 4.97±0.64 in M. galloprovincialis and 5.10±0.52 in M. trossulus (Fig. 3). This translates to a 31-fold average induction of Hsp70 mRNAs in M. galloprovincialis and 34-fold average induction in M. trossulus.

96 genes showed significantly different responses to heat stress between M. galloprovincialis and M. trossulus, 2% of all the genes assayed on the microarray (ANOVA, FDR-corrected P<0.01 for species × temperature interaction). Based on their expression profiles, the hierarchical clustering analysis identified four clusters among these genes (Fig. 4). Cluster 1 contained 28 genes, most of which were down-regulated in M. galloprovincialis and up-regulated in M. trossulus at 24°C. Cluster 2 contained 7 genes, all of which showed strong up-regulation in M. galloprovincialis at 24°C, 28°C and 32°C but showed no clear pattern in M. trossulus. Cluster 3 consisted of 21 genes, which showed no clear pattern in M. galloprovincialis but were consistently down-regulated in M. trossulus at 32°C. Cluster 4 contained 40 genes, which showed no clear pattern in M. galloprovincialis but were consistently up-regulated in M. trossulus at 32°C.

The ORA revealed key functional categories that were significantly enriched in all four clusters (Table 1). However, Cluster 1 was only significantly enriched for the ‘unknown’ functional category, indicating that, based on information from the relevant literature, Cluster 1 was functionally unclassifiable. Cluster 2 was enriched for the oxidative stress functional category, and Cluster 3 was enriched for the ion transport category (Table 1). Cluster 4 had four enriched functional categories, more than any other cluster, one of which was proteolysis. Of the four genes in the proteolysis category, three of them fell into Cluster 4. The proteolysis category was also identified in the GCA of genes with significant temperature effects (see above) and showed a similar trend in which genes were more highly up-regulated in M. trossulus than in M. galloprovincialis at 32°C (Fig. 2B and Fig. 4). The ORA of Cluster 4 demonstrated an even clearer pattern than that of the GCA: genes in the proteolysis category were strongly induced in M. trossulus but showed no change in M. galloprovincialis at 32°C.

Eight out of 34 (23.5%) of the genes associated with proteolysis, which either had significant temperature effects or species × temperature interactions, encoded subunits of the proteasome (Fig. 2B and Fig. 4). Fig. 5 shows the average expression among these eight proteasome genes, and indicates that M. trossulus showed a dramatic rise in the level of induction of proteasomal mRNAs with increased temperature whereas M. galloprovincialis showed a muted response.

The gene that showed the biggest difference between the species in response to heat stress was the molecular chaperone gene small heat shock protein 24 (Hsp24). This gene fell into Cluster 2 and showed an average of a 64-fold induction at all temperatures in M. galloprovincialis but only a small change in M. trossulus (Fig. 4). Interestingly, when compared with the levels of expression of other molecular chaperone genes (i.e. Hsp70, Hsp90, Hsp60, DnaJ and chaperonin TCP1), the differential induction of Hsp24 stands out as an important characteristic of the species-specific responses to heat stress (Fig. 6).

The ability for an organism encountering stress to regulate cellular processes by transcriptional control can allow it to cope with stress-induced damage of macromolecular components and to avoid extreme cellular damage that can lead to cell death. Here, we compared two species of mussels, the more warm-adapted invasive species M. galloprovincialis and the more cold-adapted non-invasive M. trossulus in their transcriptomic responses to acute thermal stress. Our results suggest that minor differences in transcriptomes seem to underlie the well-established differences in thermotolerance of these species, as relatively few genes exhibited species-specific transcriptional changes in response to acute heat-stress.

Our data demonstrate that, in response to acute heat stress, M. galloprovincialis and M. trossulus are highly similar at the level of the transcriptome. In particular, these species are virtually indistinguishable in their gene expression of major components of the cellular stress response. The increased expression of molecular chaperones is a primary component of the cellular stress response (Kültz, 2005) and a key indicator of environmental stress (Dahlhoff, 2004). The two species were similar in both the timing and intensity of expression of molecular chaperone genes, such as Hsp70 and genes belonging to the protein folding GO biological process (i.e. Hsp90, Hsp60, DnaJ and chaperonin TCP1). These results were unexpected because previous studies of marine mollusks that have compared differentially thermally adapted species in response to acute thermal stress have found differences in the timing (Tomanek and Somero, 1999) and/or intensity (Dong et al., 2008) of expression of Hsp70 protein. In particular, Hofmann and Somero showed that M. trossulus had a high level of constitutive expression of Hsp70 but had a lesser ability to induce this protein (Hofmann and Somero, 1996). One potential explanation for the discrepancy between our data and those of Hofmann and Somero (Hofmann and Somero, 1996) is that we did not assay protein expression but rather extrapolated from mRNA levels. However, while transcript levels are only a proxy for protein expression, other studies have shown a tight correlation between the transcription and expression of Hsps (Buckley et al., 2006). The high levels of Hsp gene hybridization observed here, as expected, suggests that our results are not an artifact of the microarray design. We note that a related proteomic study of Tomanek and Zuzow (Tomanek and Zuzow, 2010) also showed that Hsp70 induction, as measured by protein concentration, in response to acute heat stress was similar between M. galloprovincialis and M. trossulus. Therefore, further study is needed to reconcile the current findings with previous work.

Because the two species showed similar expression of molecular chaperone genes, and molecular chaperone genes showed the largest changes in expression, the more subtle changes in gene expression that differentiated the two species were relatively obscure. This is why the first four principal components (PCs 1, 2, 3 and 4) in the PCA did not reveal species-level differences. The fifth principal component (PC5) showed a clear distinction between the species, among the heat-shocked samples.

Despite the similarities in the transcriptomes of M. galloprovincialis and M. trossulus, our results also show that there are important differences. The subset of genes that showed significant species-specific responses to heat stress exhibited consistent patterns of expression, allowing them to be clustered into distinct blocks of genes that shared similar expression profiles, as well as definable functions (Fig. 4 and Table 1).

Mytlius galloprovincialis' response to acute heat stress can be primarily distinguished from that of M. trossulus by the higher induction of three oxidative stress genes and one molecular chaperone gene, Hsp24. These genes have been shown to be important elements of stress response and stress resistance in other organisms and therefore could be key factors that increase thermotolerance in M. galloprovincialis. We consider first the potential consequences of oxidative stress genes and discuss Hsp24 in the next section. Environmental stress causes cellular membranes to become compromised, which can lead to the release of reactive oxygen species (ROS) in the form of oxygen radicals from the mitochondria (Kültz, 2005; Murphy, 2009). ROS are toxic to the cell because, like other forms of stress, they cause damage to macromolecules (Kültz, 2005). Therefore, antioxidant systems that scavenge ROS are an important coping mechanism. Three out of the four oxidative stress genes that were differentially expressed between the species showed high induction in M. galloprovincialis with little or no response in M. trossulus: pyridoxal phosphate phosphatase PHOSPHO2 and β-1,3-glucuronyltransferase 2, both of which were up-regulated at all heat-shock temperatures, and glutamate decarboxylase 1, which was up-regulated in M. galloprovincialis and down-regulated in M. trossulus at 32°C. Pyridoxal phosphate phosphatase PHOSPHO2 metabolizes antioxidant compounds (Bilski et al., 2000), β-1,3-glucuronyltransferase 2 detoxifies oxidatively-damaged molecules for excretion from the cell (Stakhiv et al., 2006), and glutamate decarboxylase 1 has been shown to confer oxidative stress resistance in the yeast Saccharomyces cerevisiae (Coleman et al., 2001). The one oxidative-stress-related gene that was more highly expressed in M. trossulus (up-regulated at 32°C) than in M. galloprovincialis (slightly up-regulated at 24°C and down-regulated at 32°C) was cadmium metallothionein precursor, which scavenges toxic heavy metal ions (Nordberg, 1989) and has been shown to be up-regulated in the protist Tetrahymena thermophila in response to hydrogen peroxide-induced oxidative stress (Dondero et al., 2004). The increased expression of oxidative stress genes is an indication that both species probably encountered oxidative damage in response to acute heat stress. However, the more robust response of M. galloprovincialis in three of these genes might signify its greater ability to cope with heat-stress-induced macromolecular damage.

Indeed, the M. trossulus-specific transcriptomic response to heat stress is, in part, characterized by the up-regulation at 32°C of genes involved in proteolysis, an indicator of severe cellular stress (Dahlhoff, 2004; Kültz, 2005). Proteolysis is the directed degradation of proteins, and in the context of environmental stress serves to degrade and remove permanently denatured proteins (Travers et al., 2000; Kültz, 2003; Kültz, 2005). Mytilus galloprovincialis and M. trossulus showed similar timing of induction of many genes involved in proteolysis. However, the level of induction was much higher in M. trossulus, and M. trossulus showed strong induction of three subunits of the proteasome for which M. galloprovincialis showed no induction at all. The involvement of proteolysis deduced from our data is consistent with previous studies that have shown higher levels of protein ubiquitination, a signal of proteolysis, in M. trossulus than in M. galloprovincialis under common-garden conditions in the laboratory (Hofmann and Somero, 1996) and during summer months in the field (Dutton and Hofmann, 2008). The companion proteomic study of Tomanek and Zuzow (Tomanek and Zuzow, 2010) also showed evidence for increased proteolysis in M. trossulus but not in M. galloprovincialis.

Increased levels of stress at 32°C in M. trossulus are also indicated by the up-regulation of energy metabolism genes. Many processes involved in the cellular stress response require ATP (Kültz, 2005), so energy metabolism is crucial to mounting a robust stress response. Other species have shown induction of energy metabolism genes in response to heat stress, including the yeast S. cerevisiae (Gasch et al., 2000), the goby Gillichthys mirabilis (Buckley et al., 2006), the Antarctic fish Trematomus bernacchii (Buckley and Somero, 2009) and the mussel M. californianus (Gracey et al., 2008). Of the three genes in the energy metabolism category that showed species-specific expression, all were strongly induced in M. trossulus at 32°C. Two of these genes encode subunits of mitochondrial enzymes that are critical for ATP production: NADH dehydrogenase (ubiquinone) 1 β subcomplex subunit 8 and pyruvate dehydrogenase E1 component subunit α– somatic form. NADH dehydrogenase is the first enzyme in the electron transport chain, suggesting that heat stress induced an increase in aerobic respiration in M. trossulus. However, the companion proteomic study of Tomanek and Zuzow (Tomanek and Zuzow, 2010), showed the opposite trend (i.e. down-regulation of NADH dehydrogenase at 28°C and 32°C), suggesting downstream regulation of this enzyme. Interestingly, the NADH dehydrogenase complex has also been shown to produce ROS (Murphy, 2009), which might explain the need for the suppression of this protein despite its increased mRNA expression. The third gene in this category, patatin-like phospholipase domain-containing protein, encodes an enzyme that catalyzes the first step in triglyceride hydrolysis in lipid droplets (Zimmerman et al., 2004). This process makes stored fat available as an energy source, and thus is important in maintaining energy homeostasis when ATP is in demand (Rajala and Scherer, 2003).

Another characteristic of the M. trossulus response to heat stress is the down-regulation at 32°C of genes involved in ion transport. The differential regulation of ion transport genes suggests that gill cells in M. trossulus might have undergone changes in osmotic balance and/or cell volume. Heat stress can cause damage to membrane lipids (Swan and Watson, 1999), which can increase the permeability of the cell membrane (Kültz, 2005). Therefore, the maintenance of solute balance can be challenging for cells under conditions of environmental stress. Indeed, temperature stress has been shown to induce osmotic changes in fish gill (Stanley and Colby, 1971; Nordlie, 2009). In this context, it is noteworthy that the gene leucine aminopeptidase (Lap) was induced at 32°C in both species but more strongly in M. trossulus (Fig. 2B). Lap encodes a protein in the proteolysis category that in mussels also functions to provide free amino acids, which can be used as organic osmolytes to regulate cell volume (Koehn and Hilbish, 1987).

It is important to note that at 32°C M. trossulus up-regulated genes involved in cell signaling. The induction of cell signaling genes is difficult to interpret because these genes encode proteins that are involved in triggering a cascade of events in multiple cellular pathways. Two of these genes encode 14-3-3 proteins, which interact with other proteins by binding to phosphorylated serine and threonine groups, affecting a diverse array of cellular processes, including metabolism, protein trafficking, signal transduction, apoptosis and cell cycle regulation (Morrison, 2009). The third gene in the cell signaling functional category encodes ras-related c3 botulinum toxin substrate 1 (Rac1), which promotes actin polymerization (Wennerberg and Der, 2004) and has been shown to be involved in cell adhesion (Fukata and Kaibuchi, 2001) and phagocytosis (Nauseef, 2004). Further research is needed to fully characterize the functional relevance of these trends in heat-stress-induced gene expression of cell signaling genes.

Hsp24 showed by far the largest difference between the species in stress-induced expression. This disparity, as well as the degree to which Hsp24 was induced in M. galloprovincialis, suggests that it probably played an important role in the stress response of that species. Hofmann and Somero similarly discovered a protein of the same size class as Hsp24, which they did not identify, that was induced in M. galloprovincialis in response to heat stress and not in M. trossulus (Hofmann and Somero, 1996). In addition, the companion study of Tomanek and Zuzow (Tomanek and Zuzow, 2010) found that M. galloprovincialis and M. trossulus showed different responses to heat stress in the expression of four small Hsp isoforms. Homologs of Hsp24 belong to the alpha-crystallin family of small Hsps (de Jong et al., 1993), and have been shown to be important chaperones in the stress responses of a diverse array of animals, including humans (Horwitz, 1992), brine shrimp Artemia franciscana (Clegg et al., 1999), the annual killifish Austrofundulus limnaeus (Podrabsky and Somero, 2004) and the mussel M. californianus (Gracey et al., 2008). Small Hsps have been shown to inhibit apoptosis (Beere, 2004) and to increase resistance to oxidative stress in mammalian cells (Lee et al., 2004). Mechanistically, Hsp24 is likely to act as a chaperone for maintaining cytoskeletal structural elements during stress (Derham and Harding, 1999), as well as preventing the aggregation of denatured proteins (Horwitz, 1992; Horwitz et al., 1992; Merck et al., 1993). Protection of the cytoskeleton is likely to be a critical element of thermotolerance, particularly in gill tissue. In fact, our data show that one characteristic of the heat-stress response in M. trossulus is the induction at 32°C of genes involved in cytoskeletal reorganization. Previous studies on fish gills have shown that environmental stress has major effects on cytoskeletal structure and organization (Buckley et al., 2006; Evans and Somero, 2008), probably due to the effects on the cellular osmotic balance and volume (Di Ciano et al., 2002). Therefore, the fact that M. trossulus showed a muted response of Hsp24 might explain the need for this species to up-regulate genes in the cytoskeletal reorganization category, as well as ion transport genes. This is a further indication that heat stress had heightened effects on the molecular and cellular integrity of M. trossulus gill tissues.

The comparative approach we have employed here has given us a unique opportunity to gain insight into the evolution of transcriptional regulation, and more specifically the evolution of thermotolerance. A large fraction of the heat-stress-induced transcriptome is conserved between M. galloprovincialis and M. trossulus, which indicates that the majority of the cellular stress response is constrained, i.e. evolving under stabilizing selection (Whitehead and Crawford, 2006). This is consistent with previous work that has showed conserved transcriptional patterns across diverse taxa among genes associated with other fundamental biological processes, such as aging and development (McCarroll et al., 2004). Conversely, a small subset of genes exhibit species-specific patterns of heat-stress-induced expression. Given our experimental design, these trends in gene expression could be the result of either evolutionary divergence between the species or developmental differences arising from the distant locations of our two field collection sites. We chose our collection sites to minimize the chance of encountering hybrid individuals. Common-garden acclimation probably removed most of these site-specific environmental effects but the effect of rearing environment cannot be ruled out. Nevertheless, the fact that there are relatively few differences in heat-stress-induced gene expression between congeners that have well-documented differences in thermotolerance suggests that the differential regulation of a small subset of genes has played an adaptive role in the evolution of these species.

Mytilus galloprovincialis and M. trossulus show remarkably similar transcriptional responses to acute heat stress. However, species-specific transcriptional responses begin to elucidate differences in thermotolerance, particularly among genes involved in oxidative stress, proteolysis and cytoskeletal reorganization and the molecular chaperone gene Hsp24. The reasons for M. galloprovincialis' greater thermotolerance probably go beyond the realm of transcriptional regulation. Much of the regulation and activation of proteins occurs post-translationally via phosphorylation and acetylation events, and this will be an important next step in elucidating key differences in the physiologies of these two species.

To our knowledge this is the only study to date that has compared the transcriptomes of an invasive and a native marine species in response to thermal stress. Thus, not only does our study provide a window into the evolution of thermotolerance, it also demonstrates that molecular physiology is an important component of invasion biology. While it was surprising to find that relatively few differences in gene expression are associated with well-established differences in the thermal physiologies of these two congeners, it may help to explain how such closely-related organisms can have such dramatically different invasive capacities: while M. galloprovincialis has established populations around the world, (Braby and Somero, 2006a; Schneider and Helmuth, 2007; Schneider, 2008), there are currently no known invasive populations of M. trossulus (Seed, 1992). It remains to be seen whether invasive populations of M. galloprovincialis will continue to spread, as warmer sea surface temperatures due to global climate change will cause the expansion of ideal thermal habitat for the invasive species (Stachowicz et al., 2002; Schneider, 2008). The complexity of the natural environments in which these species live make it crucial to consider other forms of environmental stress, such as extremes of salinity, that might be limiting their biogeographical ranges.

We thank Dr Lars Tomanek of California Polytechnic State University, San Luis Obispo, CA, USA, for collecting M. galloprovincialis and maintaining mussels in his laboratory. We also thank Dr Tomanek, along with Daniel Magee, Jeremy LaBarge and Jacob Valenzuela, for assistance with the thermal stress experiment. We thank Dr Gavin Sherlock for helpful advice on the microarray hybridization design. We thank Laurie Kost for genotyping the mussels used in the study. We thank two anonymous reviewers for helpful comments on the manuscript. This study was conducted in collaboration with Dr Lars Tomanek and Marcus J. Zuzow (Tomanek and Zuzow, 2010).

 
• ANOVA

  analysis of variance

 
• aRNA

  amplified ribonucleic acid

 
• BLAST

  basic local alignment search tool

 
• EC

  enzyme commission

 
• ESTs

  expressed sequence tags

 
• FDR

  false discovery rate

 
• GCA

  gene group correlation analysis

 
• GEO

  Gene Expression Omnibus

 
• GO

  Gene Ontology

 
• Hsp

  heat shock protein

 
• Hsp24

  small heat shock protein 24 gene

 
• Hsp24

  small heat shock protein 24 protein

 
• Hsp60

  heat shock protein 60 gene

 
• Hsp70

  heat shock protein 70 gene

 
• Hsp70

  heat shock protein 70 protein

 
• Hsp90

  heat shock protein 90 gene

 
• Lap

  leucine aminopeptidase gene

 
• ORA

  over-representation analysis

 
• PC

  principal component

 
• PCA

  principal component analysis

 
• Rac1

  ras-related c3 botulinum toxin substrate 1

 
• ROS

  reactive oxygen species

 
• SSH

  suppression subtractive hybridization