Skip to main content

Investigating the molecular basis of local adaptation to thermal stress: population differences in gene expression across the transcriptome of the copepod Tigriopus californicus



Geographic variation in the thermal environment impacts a broad range of biochemical and physiological processes and can be a major selective force leading to local population adaptation. In the intertidal copepod Tigriopus californicus, populations along the coast of California show differences in thermal tolerance that are consistent with adaptation, i.e., southern populations withstand thermal stresses that are lethal to northern populations. To understand the genetic basis of these physiological differences, we use an RNA-seq approach to compare genome-wide patterns of gene expression in two populations known to differ in thermal tolerance.


Observed differences in gene expression between the southern (San Diego) and the northern (Santa Cruz) populations included both the number of affected loci as well as the identity of these loci. However, the most pronounced differences concerned the amplitude of up-regulation of genes producing heat shock proteins (Hsps) and genes involved in ubiquitination and proteolysis. Among the hsp genes, orthologous pairs show markedly different thermal responses as the amplitude of hsp response was greatly elevated in the San Diego population, most notably in members of the hsp70 gene family. There was no evidence of accelerated evolution at the sequence level for hsp genes. Among other sets of genes, cuticle genes were up-regulated in SD but down-regulated in SC, and mitochondrial genes were down-regulated in both populations.


Marked changes in gene expression were observed in response to acute sub-lethal thermal stress in the copepod T. californicus. Although some qualitative differences were observed between populations, the most pronounced differences involved the magnitude of induction of numerous hsp and ubiquitin genes. These differences in gene expression suggest that evolutionary divergence in the regulatory pathway(s) involved in acute temperature stress may offer at least a partial explanation of population differences in thermal tolerance observed in Tigriopus.


Thermal environments often vary spatially and temporally across the geographic range of a species and such variation can be an important selective force leading to genetic divergence of local populations [1]. In particular, selection resulting from the extremes in local temperature variation can have a broad impact on an organism’s genome, due to the important role of temperature on a range of biochemical and physiological processes [2]. Several studies have shown that tolerance of intertidal invertebrates to thermal stress is correlated with local temperature regimes; much of this work has focused on the tolerance limits of different species, both among unrelated taxa [3] and among congeners [46]. Results of these studies typically indicate that the geographic and ecological ranges of species are determined by thermal constraints, but the genetic basis for thermal tolerance is often not known and in some cases the observed phenotypic variation is driven by plasticity in traits rather than genetic differences among populations [7, 8]. In order for local adaptation to evolve across environmental gradients, there must be underlying genetic variation in the phenotypic trait and moderate to strong environmental selection [9]. Additionally, rates of interpopulation gene flow must remain low, or gene swamping will override the effect of diversifying selection [10].

Several empirical studies have found evidence of the genetic fixation of alternative alleles related to local adaptation of populations along latitudinal clines [11, 12] and altitudinal gradients [13, 14]. These studies suggest that local adaptation sometimes involves structural changes in single protein-coding loci with large effects [15, 16]. Point mutations in protein coding sequences can result in allelic gene products with differing thermal optima [17] and directly impact protein function (e.g., transporter affinity or enzyme activity) or alter the thermal stability of proteins, i.e., their propensity to lose native conformation at warmer temperatures. However, not all environmental adaptations involve changes in protein structure, as genes can alternatively evolve new patterns of expression. Natural selection on regulatory variation might play a major role in evolution of local adaptation [18, 19], particularly since a single change in a master regulator gene can result in widespread changes in gene expression across the transcriptome [20]. Of course, there is also the possibility that genetic changes underlying local adaptation involve a combination of both structural changes and changes in gene regulation. Efforts to differentiate these alternatives are only now becoming feasible with next generation sequencing (permitting the analysis of regulatory elements on a genome-wide scale) and studies of transcriptome responses to the environment.

Until relatively recently, assessing structural variation of genes and gene expression during thermal stress in non-model organisms has often been limited to few genes. Heat shock proteins (Hsps) have been the most important targets for investigation because they play an important role in maintaining cellular functions under environmental stress [4, 2126]. With the recent accessibility of genomic tools for use in non-model organisms, several studies have taken advantage of microarray technology for profiling expression of hundreds to thousands of genes simultaneously. These investigations have shown that the gene-regulatory response of intertidal animals to thermal stress is much more complex than previously predicted. For instance, Lockwood et al. [5] showed that different species of Mytilus mussels rely, at least partially, on different sets of genes (including Hsps) during heat shock experiments. Moreover, microarray studies have revealed classes of genes that were previously unmonitored during thermal stress, but that appear to play a significant role in the thermal response [5, 2729]. Next-generation RNA sequencing technology (RNA-seq) [30] now allows us to simultaneously assemble transcriptomes and quantify gene expression across tens of thousands of genes without any a priori genomic information. Another useful aspect of the RNA-seq approach is that the data consist of actual transcript sequences, permitting accurate determination of expression levels of closely related paralogs in gene-rich families (such as the Hsps) that might otherwise be pooled in microarray studies [31], and identification of any structural variation among copies of ortholgous genes found in different populations.

Here we initiate investigations of the genetic mechanisms that underlie population differences in thermal tolerance observed in Tigriopus californicus[32]. Besides experiencing extreme daily and seasonal temperature fluctuations (annual range of 6-33°C at one site) in their local shallow intertidal pools [33], this species is distributed over 5000 km of coastline, from Baja California, Mexico to Alaska, spanning ~35 degrees of temperate latitude. Under moderate temperature conditions, populations of T. californicus show no consistent differences in survivorship [34]. Willett [32] and Kelly et al. [35] independently demonstrated that T. californicus exhibits pronounced interpopulation differences in tolerance to high temperature stress. With large local populations, restricted gene flow, and short generational time, we hypothesize that differences in high temperature tolerance are driven by local adaptation. As a first step in understanding the molecular mechanisms that might underlie adaptation to local thermal regimes among T. californicus populations, we employed RNA-seq to characterize divergence in transcriptome-wide gene expression profiles between two populations with different thermal tolerances. We additionally examine the structural variation underlying members of the hsp70 gene family to examine whether sequence variability might account for differences between these populations.


Preparation of animals

At least 2000 Tigriopus californicus were collected from each of two natural populations, San Diego (SD: N 32.746842, W −117.254650) and Santa Cruz (SC: N 36.949589, W −122.047036), and maintained in the laboratory at 20°C with a 12-hour light:dark cycle for 1–3 generations before initiation of experiments. Two replicate samples of over 300 individual copepods of all life stages were pooled from laboratory stocks of SD and SC and flushed repeatedly in filtered (0.22 micron) seawater over two days. One sample from each population served as a control treatment, while the other sample served as the acute heat stress treatment. Pooling of large numbers of individuals into single RNA samples allows for biological averaging, and for small experimental designs, pooling has been shown to approximate non-pooled, replicated designs [36]. To initiate an acute heat stress, samples were placed in 50-ml Falcon tubes and immersed in a 35°C water bath for one hour, allowed to recover for 1 hour at 20°C, and then immediately transferred into Tri Reagent (Sigma) for RNA extraction. Samples in the control treatment were maintained at 20°C and otherwise treated identically. Samples in Tri-Reagent were disrupted with a tissue homogenizer and total RNA was then extracted following the manufacturer’s protocol.

Illumina sequencing

The Illumina mRNA sequencing kit (RS-100-0801) and multiplexing kit (PE-400-1001) were used to isolate mRNA and prepare sequencing libraries. Matching barcodes were assigned to the SD control and heat stress sample (TGACCA), as well as the SC control and heat stress samples (GCCAAT), to allow for multiplexing on two lanes of an Illumina GA-II genome analyzer flow cell. This was done to remove possible PCR bias in comparisons of control versus heat stress samples. For sequencing, the two barcoded control samples were pooled in one lane, while the two barcoded heat stress samples pooled in an adjacent lane of the flow cell. Previous experiments have suggested that there is relatively little between-lane bias [37]. Each sample was sequenced for single-end 100 basepair read lengths.

Preparation of reference assemblies

Reference transcriptomes for SD and SC were previously developed from 454 pyrosequencing of each population [38]. To improve coverage and representation of rare gene transcripts in these reference assemblies, we performed de novo assembly using the previously published 454 data and the new Illumina data. Separate assemblies were made for SD and SC sequencing reads. CLC Genomics Workbench v4.5 (CLC Bio) was used to remove adaptors and trim sequences of low quality reads, and the reads were assembled into contigs with the following assembly parameters: minimum fraction length of read overlap = 0.5, minimum sequence similarity = 0.9, minimum contig length 60 basepairs.

Assembled contigs from each population were identified as putative gene products using BLASTX searches in NCBI’s non-redundant (nr) protein database and annotated for gene function using the Blast2GO software [39, 40]. The highest scoring BLAST hit, at a threshold E-value of 10-3, was used to assign a gene name to each assembled contig. For contigs with a positive BLAST hit, Gene Ontology (GO) [41] terms were retrieved at a threshold of 10-6.

Orthology between SD and SC contigs was assessed by a reciprocal best-hit approach. We used NCBI’s executable scripts to create BLAST-searchable databases with each set of contigs and then performed BLASTN searches between the two data sets. We parsed the output keeping only pairs of sequences that were each other’s best hit, with E ≤10-20, as putatively homologous.

Read mapping and expression analysis

Illumina reads from SD and SC samples were mapped to their respective reference assembly using CLC Genomics Workbench (read mapping parameters: minimum fraction length of read overlap = 0.5, minimum sequence similarity = 0.99). The expression of each gene was determined by calculating the sum of the reads mapping to each contig, and then normalized by dividing the number of observed reads by the length of the contig and by the total number of reads sequenced across the transcriptome (reads per kilobase per million reads, or RPKM). To compare levels of expression between control and heat stressed samples in each population, we used Kal’s Z-test of proportions [42]. The Z-test calculates the difference in the proportion of reads observed between two treatments and compares whether each sample was drawn from the same distribution, using the normal distribution as an approximation of the binomial distribution. Due to the large number of tests, statistical significance was assessed using p-values adjusted by Bonferroni correction. In some cases Bonferroni can be overly conservative, so we consider tests as significant with a Bonferroni p-value ≤ 0.10. This was more conservative than applying a false discovery rate at p < 0.05, which detected many more (>10x) contigs. We then filtered the dataset of significant loci and retained only genes where at least one of the orthologous contigs had >10 mapped reads. For graphical purposes, RPKM values were log10 transformed due to the large range in expression values across genes. Using the GO terms for each transcriptome, we examined whether the group of genes showing a significant change in expression in each population were enriched for particular gene ontology terms relative to their reference transcriptome. Statistical significance was assessed using Fisher’s exact tests implemented in the Blast2GO software.

Evolutionary analysis of the Hsp70 family

As discussed below, RNA-seq results identified 18 loci annotated as hsp70. Because these genes were prominent in the genes responsive to thermal stress, we carried out a more detailed analysis of this gene family. We performed additional BLASTX against the Drosophila melanogaster and Caenorhabditis elegans genomes in order to identify T. californicus contigs that were paralogous hsp70 loci. For each identified hsp70 copy, we aligned the orthologs between SD and SC and performed comparative analyses of sequence divergence. To look for evidence of molecular adaptation, coding sequences were compared to calculate the ratio of non-synonymous (dN) to synonymous (dS) nucleotide substitution rates [43]. Typically, dN/dS ratios much greater than one indicate strong diversifying selection, whereas values close to one are compatible with neutral evolution and values significantly less than one are consistent with balancing or purifying selection. To examine relationships among the Hsp70 proteins, orthologs of Drosophila melanogaster [GenBank: NM_169441.1, NM_079632.5, NM_079615.3, NM_080059.2, NM_169469.1, NM_080188.2, NM_141952.2, NM_167306.1, NM_079017.2], Caenorhabditis elegans [EMBL: F26D10, C15H9, F43E2, C37H5, F44E5.4, F44E5.5, C12C8, F11F1], and Escherichia coli [UniProt: P77319] were obtained as out-group sequences and amino acids were aligned using the ClustalW2 algorithm [44]. A neighbor-joining tree was constructed in Mega5 [45], using the p-distances model with pairwise deletion of gaps and missing data. Support for topological relationships in the phylogram was evaluated by bootstrapping with 1,000 replicates.

Validation of RNA-seq expression levels

In order to evaluate the accuracy of the RNA-seq expression profiles, we measured the expression of a set of heat shock genes using quantitative PCR (qPCR). We first quantified expression using the same RNA samples prepared for RNA-seq, as a means of technical replication. cDNA from the transcriptome samples was diluted to ~100 ng/ul and 1ul was amplified over 25 cycles with 200 nM of each gene specific primer (Additional file 1) using Brilliant III Ultra-Fast SYBR Green QPCR Mastermix (Agilent Technologies). Several housekeeping genes were chosen to normalize the expression of target genes, including beta-actin, glyceraldehyde 3-phosphate dehydrogenase, tubulin, and myosin. Ct values for each sample were normalized by geometric averaging of the housekeeping genes, implemented in the program GeNorm [46].

We assessed biological accuracy further by replicating the thermal stress experiment with additional samples (biological replicates). For each replicate, thirty individuals were placed in 15-mL tubes containing filtered seawater, and the experimental thermal stress was repeated as above. RNA was extracted using Tri-Reagent, with tissue homogenization performed using a Mini-BeadBeater 8 (BioSpec) with silica beads (0.175 mg) at full speed for 30s. RNA was quantified by spectrophotometry, and 500 ng of RNA were used for cDNA synthesis with the High Capacity RNA-to-cDNA kit (Applied Biosystems). The experiment was performed with three replicates per temperature for each population. Quantitative PCR was performed in 20-μL reactions with the iTaq Fast SYBR Green Supermix (Bio-Rad), and contained 2 μL of a 1:10 dilution of cDNA and 250 nM of each primer. Fluorescence was measured in a Stratagene Mx 3000P System (Stratagene), with the following thermal profile: 95°C for 3 min, followed by 40 cycles of 95°C for 10 s and 59°C for 30 s. At the end of thermal cycling, a melting curve analysis was performed to confirm the presence of single amplification products. A total of 14 target genes were surveyed, and their expression levels were normalized to the expression of housekeeping genes.


Illumina sequencing, de novo assembly, and RNA-seq mapping

The sequencing data from this study have been submitted to the NCBI Gene Expression Omnibus under accession no. GSE38546. Single-end 100 basepair Illumina sequencing of the SD samples resulted in 9,292,573 reads, with 6,794,534 reads in the heat stress sample and 2,498,039 reads in the control sample. Sequencing of the SC samples yielded 10,794,996 reads, with 6,205,028 reads in the heat stress sample and 3,256,758 in the control sample. We combined samples from each population with previous 454 sequencing data [38] and assembled the reads into 66,943 unique SD contigs and 60,648 unique SC contigs. The average length of the SD contigs was 498 basepairs (range: 49–9,830 bp) and the average length of the SC contigs was 520 bp (range: 36–9,003 bp). These contigs formed the reference transcriptomes used in our RNA-seq gene expression analysis. Approximately 29% and 33% of the SD and SC assemblies, respectively, were annotated using BLAST, and a total of 19,622 contigs were identified as orthologs based on reciprocal BLAST to each transcriptome. Approximately 78% of the total reads were uniquely mapped in the SD control sample, while ~82% reads were uniquely mapped in the SD heat stress sample. In the SC control sample, ~71% reads were uniquely mapped, while ~80% reads were uniquely mapped in the SC heat stress sample. A distribution of the RPKM values shows that most of the contigs in each treatment had values around 10 RPKM, with slightly lower values in the heat stress samples (see Additional file 2).

Differential gene expression

In the comparison of SD control and heat shock samples, a total of 590 contigs (0.88% of 66,943 contigs) showed evidence of a statistically significant difference in expression (Bonferroni corrected p <0.10). Of these 590 contigs, 234 were removed due to a low number of mapped reads, leaving 356 (0.53%) as significantly up or down regulated. In the comparison of SC control and heat shock samples, a total of 903 contigs (1.49% of 60,648 contigs) showed evidence of a statistically significant difference in expression (Bonferroni corrected p <0.10). Of these 903 contigs, 36 were removed due to a low number of mapped reads, leaving 867 (1.43%) as significantly up or down regulated. Most of the contigs identified as significant had RPKM values between 10 and 100, with better coverage than the overall transcriptome (see Additional file 2).

We examined the functional classes of genes that showed the significant changes in expression and tested for enrichment of gene ontology terms relative to their frequency in each transcriptome. Among the gene ontology terms that were significantly enriched (Table 1), genes responding to stress and genes involved in cuticle formation were the most overrepresented in SD and ribosomal genes and lipid transport genes were the most overrepresented in the SC population. We also focused on the pattern of expression in four classes of genes that showed significant fold changes in both populations: heat shock proteins, genes associated with proteolysis, mitochondrial genes, and genes associated with cuticle formation (see Additional file 3 for frequency of other gene ontology terms). We plotted the log-transformed expression of genes that were significantly up or down-regulated in the heat-shock treatment relative to control samples (Figure 1). Plots of the SD heat stress versus control sample (Figure 1A) show that variation in up-regulated gene expression was greater than in down-regulated genes, with a group of 125 genes exhibiting a >10X-fold up-regulation relative to the control sample. Thirteen genes with >10X-fold up-regulation were identified as heat shock proteins (Hsps). Nineteen genes involved in the degradation of damaged proteins (proteolysis associated genes) were up-regulated, as were 20 genes involved in cuticle formation. Down-regulated genes include eight mtDNA-encoded genes involved in cellular respiration, as well as some other classes of genes (e.g. lipid transport and hydrolases, not shown). In comparison, plots of the SC heat stress versus control sample (Figure 1B) show relatively even distribution of up-regulated and down-regulated genes, with much less variation in expression than the SD treatment. Only 15 genes show a >10X-fold up-regulation, and two of these are hsps. Seventeen proteolytic genes were up-regulated and eighteen were down-regulated in the SC comparison. Thirteen mitochondrial genes in SC were down-regulated and show roughly the same response as in SD. In striking contrast to the SD comparison, thirteen SC genes involved in cuticle formation were largely down-regulated.

Table 1 Enrichment of gene ontology (GO) terms in genes that were differentially expressed during heat stress
Figure 1

Statistically significant up/down-regulated genes in A) San Diego and B) Santa Cruz heat stress treatments. Scatterplots of A) San Diego Control versus San Diego Heat Stress and B) Santa Cruz Control versus Santa Cruz Heat Stress. Expression values are normalized for contig length and per million reads and log-transformed. Open circles are unidentified genes, closed black circles are BLAST-classified genes, and closed colored circles represent specific groups of classified genes.

Although the direction of change in expression under heat stress is the same in both populations across most genes (with cuticle-associated genes being an exception), the magnitude of up- or down-regulation is markedly different. In the SD sample, several genes showed >100-fold change in expression, while the greatest change in the SC sample was 25-fold (Table 2; Figure 1). An examination of the identity of genes with statistically significant changes in expression revealed that 63 were shared in both populations (Table 2). Twenty-five of these were up-regulated in both populations, including eight heat shock proteins. In all cases, the heat shock proteins showed higher fold changes in the SD population. Fourteen genes were down-regulated in both populations, including 7 mitochondrial genes with very similar levels of change in both populations. Finally, 24 genes showed population-specific direction of responses, 12 of which belonged to two functional classes showing opposite expression profiles. Seven cuticle proteins were up-regulated in SD (mean = 16.1-fold) and down-regulated in SC (−3.1-fold), while 5 hydrolases were down-regulated in SD (−2.9-fold) and up-regulated in SC (2.4-fold) (Table 2).

Table 2 Genes found in both populations with statistically significant changes in expression after thermal stress

When comparing transcriptional profiles of hsps, some population-specific responses to heat shock were also observed (Table 3). A total of 29 putatively unique hsp transcripts were identified as orthologs, only nine of which showed significant change in the same direction (up-regulation) in both samples. For each of 10 SC hsp70 paralogs that can be matched to an SD ortholog, the patterns of induction during thermal stress are similar in the two populations. However, the amplitude of the stress response is markedly different, with the SD orthologs consistently showing a higher level of stress induction in the more thermally tolerant SD population. This pattern holds across other hsp genes, with hsp induction invariably higher in SD. For 13 putative hsps, an orthologous copy was not identified in one population through reciprocal BLAST. These genes, missing from one transcriptome, most likely reflect an incomplete transcriptome assembly; alternatively a gene may have been lost in one of the populations. These hypotheses not yet been tested via PCR amplification from genomic DNA.

Table 3 RNA-seq quantification of gene expression in heat shock proteins

Sequence divergence in the hsp70 gene family

Average sequence divergence of the 29 orthologous hsp70 pairs is 2.9% and ranges from 0.0- 8.7% (Table 3). The ratio of rates of non-synonymous to synonymous substitutions (dN/dS) is less than one in all pairs, providing no evidence for divergent selection. Based on full-length transcript sequences, we unambiguously identified eight hsp70 paralogs in the T. californicus transcriptome. We examined the relationship between these paralogs and Hsp70 proteins from Drosophila melanogaster and Caenorhabditis elegans. A neighbor-joining tree based on protein alignment (Figure 2) shows that five pairs have an ancient origin among Metazoans, while 3 pairs have diversified from an ancestral copy of Hsp70 #2. In all cases, orthologous copies in SD and SC are more closely related to one another than to other Hsp70 paralogs. Notably, the relatively closely related group of paralogs (#1, #2, #5, and #16) show markedly different response to thermal stress.

Figure 2

Neighbor-joining phylogram of Hsp70 proteins from Tigriopus californicus and out-groups. Numbers next to internal nodes represent bootstrap values from 1,000 replicates. For T. californicus, numbers next to population labels show the fold-change during thermal stress (‘n.s.’ = not significantly up- or down-regulated), as quantified by RNA-seq, while numbers associated with each pair of orthologs show the percent divergence in amino acid sequence between populations. NCBI accession numbers for T. californicus Hsp sequences can be found in Table 3.

Validation of RNA-seq gene expression levels

We performed qPCR assays on a subset of genes using the same RNA samples prepared for RNA-seq (a type of technical validation), as well as using RNA from additional samples (for biological validation). In nearly all comparisons (75%), fold-changes in expression measured by qPCR were in the same direction as those estimated by RNA-seq (Table 4). Moreover, the magnitude of these changes were also highly congruent between the two methods, as tested by linear regressions of log2 fold-differences (technical replicates: r = 0.85, F1,34 = 87.4, P < 0.001; biological replicates: r = 0.84, F1,26 = 68.2, P < 0.001).

Table 4 Comparison of RNA-seq and qPCR measurements of fold change in gene expression in selected genes


Population divergence in acute heat stress

Populations of Tigriopus californicus show strikingly different survivorship following acute thermal stress [32], with populations in southern California (i.e. San Diego) exhibiting higher survivorship. Here we tested the hypothesis that differences in patterns of gene expression play a role in adaptation to thermal stress. Our results, based on an RNA-seq approach, demonstrated that the two populations under study have genetically diverged in their transcriptome response to acute thermal stress. While the number of reads obtained for each treatment in each population limited our power to detect rare transcripts and incomplete assemblies could result in uneven representations of genes in each population, our qPCR data confirmed that the RNAseq results are robust measurements of gene expression changes in the tested populations. Our experiment showed that statistically significant changes in expression involved a small number of genes across the transcriptome of T. californicus (356 in San Diego and 867 in Santa Cruz, accounting for less than 1.5% of their respective assemblies). Among the subset of genes with statistically significant changes in expression, only 63 were responsive in both populations, suggesting that the two populations have diverged in the identity of affected loci. Furthermore, population differences were evident in both the direction and magnitude of expression changes in these 63 genes. The functional groups represented by these genes also differed between the two populations, as different gene ontology terms were enriched relative to the transcriptome. The San Diego population from southern California is more tolerant of thermal stress and exhibited a much larger up-regulation of stress-related genes following heat shock. In particular, this included a large number of genes that are members of heat shock protein families. Although there is up-regulation of hsp70 paralogs in Santa Cruz, the fold change is much smaller than observed in the San Diego population. These results contrast with numerous other studies over latitudinal gradients that have not detected interpopulation expression differences in hsps[26]. For example, even at the interspecific level, Lockwood et al. [5] found that the magnitude of expression of seven hsp70 genes was strikingly similar between a northern and a southern species of Mytilus mussels. Only one hsp (hsp24) exhibited a strong difference in expression between the two geographically separated species.

We also examined whether structural variation in Hsp orthologs might explain differential survivorship following acute thermal stress. Sequence analyses of the Hsp gene families suggest population differentiation at the structural level is not elevated compared to the rest of the genome. We identified 42 members of Hsp families across the San Diego and Santa Cruz populations and 29 were assigned to a matching ortholog. While some differences are evident in these orthologous copies, the amount of sequence divergence (average 2.9%) is close to the median level of 2.7% sequence divergence found across all genes in a comparison of the San Diego and Santa Cruz transcriptomes [38]. The orthologous pair with the largest amount of sequence divergence in our analysis (8.7%, hsp67b2 #1) is not an outlier considering the range of sequence divergence evident across the entire transcriptome [38]. Finally, the lack of evidence that non-synonymous substitution rates are elevated in these genes, suggests that positive selection on the coding region may have a limited role in the diversification of heat shock protein function. We note that non-synonymous substitution rates were measured across entire genes; more rigorous tests of selection based on codon-substitution models will require sampling individuals across multiple populations of T. californicus[47].

In addition to the Hsp genes, several other functional categories of genes showed significant response to thermal stress. One similarity between the two populations was the pattern of down-regulation of genes encoded in the mtDNA. Since certain reactions in the electron transport chain are known to be major sources of reactive oxygen species [48], a decrease in expression of mtDNA genes may represent a means of reducing stress. In contrast to the mtDNA genes, most orthologous genes showing significant change in expression in response to thermal stress have diverged either in the magnitude of expression or direction of regulation. Perhaps the most striking difference observed was in genes involved in cuticle formation; these were enriched and up-regulated in San Diego, while being down-regulated in Santa Cruz. Cuticle genes are known to be up-regulated during diapause in many invertebrates [49], presumably to reduce cuticle permeability, but their specific role during acute thermal stress is not known. Genes functioning in proteolysis, which are known to play an important role in protein degradation following stress, also showed a differential response with significant up-regulation of some genes in southern but not in northern populations of Tigriopus. This pattern is opposite to the pattern found in Mytilus[5], where northern populations exhibit an up-regulation of proteolytic genes.

In sum, our results show a variety of significant differences between populations in gene regulation following acute thermal stress and we hypothesize that these contribute to the population differences in survivorship following stress. While mRNA transcript abundance is not equivalent to protein abundance [50], other studies have shown that gene expression changes in experimental treatments are relevant predictors of organismal physiology [51], and that there is often a high correlation between gene expression and protein abundance, including the Hsp gene family [52].

Evidence of local adaptation in Tigriopuscalifornicus

Numerous species exhibit morphological and life history variation along latitudinal clines [53, 54]. These traits typically involve changes in body size, phenology [11, 12, 55], or temperature adaptation [56]. Recently, Willett [32] and Kelly et al. [35] provided evidence that the intertidal copepod Tigriopus californicus exhibits variation in temperature tolerance along a latitudinal cline in California. Willett showed that populations in northern California exhibit reduced survivorship when exposed to temperatures above 36°C for one hour, while survivorship remained high in southern California populations. Additionally, southern California populations exposed to chronic temperature stress survived longer than populations from northern California. Based on observed tradeoffs in fitness across a range of temperatures, Willett argued that this provided evidence of local adaptation to temperature in T. californicus. Kelly et al. [35] used selection experiments to demonstrate that thermal tolerance of T. californicus responded to selection over 10 generations, providing direct evidence for additive genetic variation for tolerance to thermal stress. Furthermore, Kelly et al. found that total quantitative genetic variation was almost completely partitioned among populations, rather than within populations, suggesting strong local adaptation.

Our results support the conclusions of both of the above studies and provide a possible molecular basis for local adaptation, namely that the pathways that control the regulation of numerous genes involved in response to acute temperature stress have diverged. Other studies have shown that stabilizing selection often acts to prevent genetic differentiation in expression patterns among populations and among species [57, 58]. The clearest pattern in our study emerges from the dramatic differences in expression in the hsp70 gene. Heat shock genes are known to be important components of temperature adaptation and resistance to thermal stress [1, 59], as they play a critical role in preventing the degradation of intracellular proteins. Whole transcriptome studies have shown that Hsps play a dominant role in heat stress response in other species [60], as well as responding to other forms of environmental stress [61, 62]. Rapid adaptation to novel thermal environments by altering expression of Hsps, and that of other chaperone proteins, may be a common form of local population adaptation in diverse taxa [5].


In this study we contrasted the whole genome response in gene expression between two T. californicus populations that exhibit differences in tolerance to thermal stress. Broad changes in gene expression were observed, with the more thermal tolerant San Diego population showing more pronounced up-regulation of hsp and ubiquitin pathway genes than the more sensitive Santa Cruz population. Although this up-regulation involves a large number of heat shock proteins, the most dramatic response involved the hsp70 gene family. Orthologous copies of the Hsps do not show evidence of divergent selection for structural changes, suggesting, at least in part, that regulatory variation between populations may play an important role in their differential sensitivity to acute heat stress.


  1. 1.

    Angilletta MJ: Thermal adaptation: a theoretical and empirical synthesis. 2009, Oxford: Oxford University Press

    Google Scholar 

  2. 2.

    Clarke A: Costs and consequences of evolutionary temperature adaptation. Trends Ecol Evol. 2003, 18: 573-581. 10.1016/j.tree.2003.08.007.

    Article  Google Scholar 

  3. 3.

    Davenport J, Davenport JL: Effects of shore height, wave exposure and geographical distance on thermal niche width of intertidal fauna. Mar Ecol Prog Ser. 2005, 292: 41-50.

    Article  Google Scholar 

  4. 4.

    Sorte CJB, Hofmann GE: Changes in latitude, changes in aptitudes: Nucella canaliculata (Mollusca: Gastropoda) is more stressed at its range edge. Mar Ecol Prog Ser. 2004, 274: 263-268.

    Article  Google Scholar 

  5. 5.

    Lockwood BL, Sanders JG, Somero GN: Transcriptomic responses to heat stress in invasive and native blue mussels (genus Mytilus): molecular correlates of invasive success. J Exp Biol. 2010, 213: 3548-3558. 10.1242/jeb.046094.

    PubMed  CAS  Article  Google Scholar 

  6. 6.

    Kuo ESL, Sanford E: Geographic variation in the upper thermal limits of an intertidal snail: implications for climate envelope models. Mar Ecol Prog Ser. 2009, 388: 137-146.

    Article  Google Scholar 

  7. 7.

    James AC, Azevedo RBR, Partridge L: Genetic and environmental responses to temperature of Drosophila melanogaster from a latitudinal cline. Genetics. 1997, 146: 881-890.

    PubMed  CAS  PubMed Central  Google Scholar 

  8. 8.

    Ayrinhac A, Debat V, Gibert P, Kister A-G, Legout H, Moreteau B, Vergilino R, David JR: Cold adaptation in geographical populations of Drosophila melanogaster: phenotypic plasticity is more important than genetic variability. Funct Ecol. 2004, 18: 700-706. 10.1111/j.0269-8463.2004.00904.x.

    Article  Google Scholar 

  9. 9.

    Kawecki TJ, Ebert D: Conceptual issues in local adaptation. Ecol Lett. 2004, 7: 1225-1241. 10.1111/j.1461-0248.2004.00684.x.

    Article  Google Scholar 

  10. 10.

    Lenormand T: Gene flow and the limits to natural selection. Trends Ecol Evol. 2002, 17: 183-189. 10.1016/S0169-5347(02)02497-7.

    Article  Google Scholar 

  11. 11.

    Li Y, Huang Y, Bergelson J, Nordborg M, Borevitz JO: Association mapping of local climate-sensitive quantitative trait loci in Arabidopsis thaliana. Proc Natl Acad Sci USA. 2010, 107: 21199-21204. 10.1073/pnas.1007431107.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  12. 12.

    Ikeda H, Setoguchi H: Natural selection on PHYE by latitude in the Japanese archipelago: insight from locus specific phylogeographic structure in Arcterica nana (Ericaceae). Mol Ecol. 2010, 19: 2779-2791. 10.1111/j.1365-294X.2010.04700.x.

    PubMed  CAS  Article  Google Scholar 

  13. 13.

    Storz J, Sabatino S, Hoffmann F, Gering E, Moriyama H, Ferrand N, Monteiro B, Nachman M: The molecular basis of high-altitude adaptation in deer mice. PLoS Genet. 2007, 3: e45-10.1371/journal.pgen.0030045.

    PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Dahlhoff EP, Rank NE: Functional and physiological consequences of genetic variation at phosphoglucose isomerase: Heat shock protein expression is related to enzyme genotype in a montane beetle. Proc Natl Acad Sci USA. 2000, 97: 10056-10061. 10.1073/pnas.160277697.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  15. 15.

    Hoekstra HE, Drumm KE, Nachman MW: Ecological genetics of adaptive color polymorphism in pocket mice: geographic variation in selected and neutral genes. Evolution. 2004, 58: 1329-1341.

    PubMed  CAS  Article  Google Scholar 

  16. 16.

    Schulte PM: Environmental adaptations as windows on molecular evolution. Comp Bioch Phys B. 2001, 128: 597-611. 10.1016/S1096-4959(00)00357-2.

    CAS  Article  Google Scholar 

  17. 17.

    Arnold FH, Wintrode PL, Miyazaki K, Gershenson A: How enzymes adapt: lessons from directed evolution. Trends Biochem Sci. 2001, 26: 100-106. 10.1016/S0968-0004(00)01755-2.

    PubMed  CAS  Article  Google Scholar 

  18. 18.

    McDonald JF, Chambers GK, David J, Ayala FJ: Adaptive response due to changes in gene regulation: a study with Drosophila. Proc Natl Acad Sci USA. 1977, 74: 4562-4566. 10.1073/pnas.74.10.4562.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  19. 19.

    Barrier M, Robichaux RH, Purugganan MD: Accelerated regulatory gene evolution in an adaptive radiation. Proc Natl Acad Sci USA. 2001, 98: 10208-10213. 10.1073/pnas.181257698.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  20. 20.

    Carroll SB: Endless forms: the evolution of gene regulation and morphological diversity. Cell. 2000, 101: 577-580. 10.1016/S0092-8674(00)80868-5.

    PubMed  CAS  Article  Google Scholar 

  21. 21.

    Hofmann G, Somero G: Evidence for protein damage at environmental temperatures: seasonal changes in levels of ubiquitin conjugates and hsp70 in the intertidal mussel Mytilus trossulus. J Exp Biol. 1995, 198: 1509-1518.

    PubMed  CAS  Google Scholar 

  22. 22.

    Hofmann G, Somero G: Interspecific variation in thermal denaturation of proteins in the congeneric mussels Mytilus trossulus and M. galloprovincialis: evidence from the heat-shock response and protein ubiquitination. Mar Biol. 1996, 126: 65-75. 10.1007/BF00571378.

    CAS  Article  Google Scholar 

  23. 23.

    Sagarin RD, Somero GN: Complex patterns of expression of heat-shock protein 70 across the southern biogeographical ranges of the intertidal mussel Mytilus californianus and snail Nucella ostrina. J Biogeogr. 2006, 33: 622-630. 10.1111/j.1365-2699.2005.01403.x.

    Article  Google Scholar 

  24. 24.

    Dutton JM, Hofmann GE: Biogeographic variation in Mytilus galloprovincialis heat shock gene expression across the eastern Pacific range. J Exp Mar Biol Ecol. 2009, 376: 37-42. 10.1016/j.jembe.2009.06.001.

    CAS  Article  Google Scholar 

  25. 25.

    Hammond LM, Hofmann GE: Thermal tolerance of Strongylocentrotus purpuratus early life history stages: mortality, stress-induced gene expression and biogeographic patterns. Mar Biol. 2010, 157: 2677-2687. 10.1007/s00227-010-1528-z.

    PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Feder ME, Hofmann GE: Heat-shock proteins, molecular chaperones, and the stress response: evolutionary and ecological physiology. Annu Rev Physiol. 1999, 61: 243-282. 10.1146/annurev.physiol.61.1.243.

    PubMed  CAS  Article  Google Scholar 

  27. 27.

    Gracey AY, Chaney ML, Boomhower JP, Tyburczy WR, Connor K, Somero GN: Rhythms of gene expression in a fluctuating intertidal environment. Curr Biol. 2008, 18: 1501-1507. 10.1016/j.cub.2008.08.049.

    PubMed  CAS  Article  Google Scholar 

  28. 28.

    Place SP, O'Donnell MJ, Hofmann GE: Gene expression in the intertidal mussel, Mytilus californianus: physiological response to environmental factors on a biogeographic scale. Mar Ecol Prog Ser. 2008, 356: 1-14.

    CAS  Article  Google Scholar 

  29. 29.

    Stillman JH, Tagmount A: Seasonal and latitudinal acclimatization of cardiac transcriptome responses to thermal stress in porcelain crabs, Petrolisthes cinctipes. Mol Ecol. 2009, 18: 4206-4226. 10.1111/j.1365-294X.2009.04354.x.

    PubMed  CAS  Article  Google Scholar 

  30. 30.

    Wang Z, Gerstein M, Snyder M: RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009, 10: 57-63. 10.1038/nrg2484.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  31. 31.

    Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Meth. 2008, 5: 621-628. 10.1038/nmeth.1226.

    CAS  Article  Google Scholar 

  32. 32.

    Willett CS: Potential fitness trade-offs for thermal tolerance in the intertidal copepod Tigriopus californicus. Evolution. 2010, 64: 2521-2534. 10.1111/j.1558-5646.2010.01008.x.

    PubMed  Article  Google Scholar 

  33. 33.

    Powlick JJ: Habitat characters of Tigriopus californicus (Copepoda: Harpacticoida), with notes on the dispersal of supralittoral fauna. J Mar Biol Assoc UK. 1999, 79: 85-92. 10.1017/S0025315498000095.

    Article  Google Scholar 

  34. 34.

    Edmands S, Deimler JK: Local adaptation, intrinsic coadaptation and the effects of environmental stress on interpopulation hybrids in the copeod Tigriopus californicus. J Exp Mar Biol Ecol. 2004, 303: 183-196. 10.1016/j.jembe.2003.11.012.

    Article  Google Scholar 

  35. 35.

    Kelly MW, Sanford E, Grosberg RK: Limited potential for adaptation to climate change in a broadly distributed marine crustacean. Proc R Soc Lond B Biol Sci. 2012, 279: 349-356. 10.1098/rspb.2011.0542.

    Article  Google Scholar 

  36. 36.

    Kendziorski C, Irizarry RA, Chen KS, Haag JD, Gould MN: On the utility of pooling biological samples in microarray experiments. Proc Natl Acad Sci USA. 2005, 102: 4252-4257. 10.1073/pnas.0500607102.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  37. 37.

    Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y: RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008, 18: 1509-1517. 10.1101/gr.079558.108.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  38. 38.

    Barreto FS, Moy GW, Burton RS: Interpopulation patterns of divergence and selection across the transcriptome of the copepod Tigriopus californicus. Mol Ecol. 2011, 20: 560-572. 10.1111/j.1365-294X.2010.04963.x.

    PubMed  Article  Google Scholar 

  39. 39.

    Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005, 21: 3674-3676. 10.1093/bioinformatics/bti610.

    PubMed  CAS  Article  Google Scholar 

  40. 40.

    Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talón M, Dopazo J, Conesa A: High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008, 36: 3420-3435. 10.1093/nar/gkn176.

    PubMed  PubMed Central  Article  Google Scholar 

  41. 41.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, The Gene Ontology Consortium, et al: Gene Ontology: tool for the unification of biology. Nat Genet. 2000, 25: 25-29. 10.1038/75556.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  42. 42.

    Kal AJ, van Zonneveld AJ, Benes V, van den Berg M, Koerkamp MG, Albermann K, Strack N, Ruijter JM, Richter A, Dujon B, et al: Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sources. Mol Biol Cell. 1999, 10: 1859-1872.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  43. 43.

    Miyata T, Yasunaga T: Molecular evolution of mRNA: a method for estimating evolutionary rates of synonymous and amino acid substitutions from homologous nucleotide sequences and its applications. J Mol Evol. 1980, 16: 23-36. 10.1007/BF01732067.

    PubMed  CAS  Article  Google Scholar 

  44. 44.

    Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, et al: Clustal W and Clustal X version 2.0. Bioinformatics. 2007, 23: 2947-2948. 10.1093/bioinformatics/btm404.

    PubMed  CAS  Article  Google Scholar 

  45. 45.

    Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar SF: MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011, 28: 2731-2739. 10.1093/molbev/msr121.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  46. 46.

    Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3: research0034.0031-research0034.0011.

    Article  Google Scholar 

  47. 47.

    Schoville SD, Flowers JM, Burton RS: Diversifying selection underlies the origin of allozyme polymorphism at the phosphoglucose isomerase locus in Tigriopus californicus. PLoS One. 2012, 7: e40035-10.1371/journal.pone.0040035.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  48. 48.

    Murphy MP: How mitochondria produce reactive oxygen species. Biochem J. 2009, 417: 1-13. 10.1042/BJ20081386.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  49. 49.

    MacRae T: Gene expression, metabolic regulation and stress tolerance during diapause. Cell Mol Life Sci. 2010, 67: 2405-2424. 10.1007/s00018-010-0311-0.

    PubMed  CAS  Article  Google Scholar 

  50. 50.

    Feder ME, Walser JC: The biological limitations of transcriptomics in elucidating stress and stress responses. J Evol Biol. 2005, 18: 901-910. 10.1111/j.1420-9101.2005.00921.x.

    PubMed  CAS  Article  Google Scholar 

  51. 51.

    Maier T, Güell M, Serrano L: Correlation of mRNA and protein in complex biological samples. FEBS Lett. 2009, 583: 3966-3973. 10.1016/j.febslet.2009.10.036.

    PubMed  CAS  Article  Google Scholar 

  52. 52.

    Buckley BA, Gracey AY, Somero GN: The cellular response to heat stress in the goby Gillichthys mirabilis: a cDNA microarray and protein-level analysis. J Exp Biol. 2006, 209: 2660-2677. 10.1242/jeb.02292.

    PubMed  CAS  Article  Google Scholar 

  53. 53.

    Endler JA: Geographic Variation, Speciation, and Clines. 1977, Princeton, NJ: Princeton University Press

    Google Scholar 

  54. 54.

    Vermeij GJ: Biogeography and Adaptation: Patterns of Marine Life. 1978, Cambridge: Mass.: Harvard University Press

    Google Scholar 

  55. 55.

    Ragland GJ, Kingsolver JG: Influence of seasonal timing on thermal ecology and thermal reaction norm evolution in Wyeomyia smithii. J Evol Biol. 2007, 20: 2144-2153. 10.1111/j.1420-9101.2007.01428.x.

    PubMed  CAS  Article  Google Scholar 

  56. 56.

    Sarup P, Frydenberg J, Loeschcke V: Local adaptation of stress related traits in Drosophila buzzatii and Drosophila simulans in spite of high gene flow. J Evol Biol. 2009, 22: 1111-1122. 10.1111/j.1420-9101.2009.01725.x.

    PubMed  CAS  Article  Google Scholar 

  57. 57.

    Lemos B, Meiklejohn CD, Cáceres M, Hartl DL: Rates of divergence in gene expression profiles of primates, mice, and flies: stabilizing selection and variability among functional categories. Evolution. 2005, 59: 126-137.

    PubMed  CAS  Article  Google Scholar 

  58. 58.

    Whitehead A, Crawford DL: Variation within and among species in gene expression: raw material for evolution. Mol Ecol. 2006, 15: 1197-1211. 10.1111/j.1365-294X.2006.02868.x.

    PubMed  CAS  Article  Google Scholar 

  59. 59.

    Sørensen JG, Kristensen TN, Loeschcke V: The evolutionary and ecological role of heat shock proteins. Ecol Lett. 2003, 6: 1025-1037. 10.1046/j.1461-0248.2003.00528.x.

    Article  Google Scholar 

  60. 60.

    Sørensen JG, Nielsen MM, Kruhøffer M, Justesen J, Loeschcke V: Full genome gene expression analysis of the heat stress response in Drosophila melanogaster. Cell Stress Chaperon. 2005, 10: 312-328. 10.1379/CSC-128R1.1.

    Article  Google Scholar 

  61. 61.

    Chapman RW, Mancia A, Beal M, Veloso A, Rathburn C, Blair A, Holland AF, Warr GW, Didinato GUY, Sokolova IM, et al: The transcriptomic responses of the eastern oyster, Crassostrea virginica, to environmental conditions. Mol Ecol. 2011, 20: 1431-1449. 10.1111/j.1365-294X.2011.05018.x.

    PubMed  Article  Google Scholar 

  62. 62.

    Rhee J-S, Raisuddin S, Lee K-W, Seo JS, Ki J-S, Kim I-C, Park HG, Lee J-S: Heat shock protein (Hsp) gene responses of the intertidal copepod Tigriopus japonicus to environmental toxicants. Comp Biochem Phys C. 2009, 149: 104-112.

    Google Scholar 

Download references


We would like to acknowledge The Scripps Research Institute DNA Array Core Facility for conducting the Illumina Sequencing. We would also like to thank the anonymous reviewers for their helpful comments. This study was funded in part by grants from the National Science Foundation (DEB 0717178, DEB 1051057) to RSB.

Author information



Corresponding author

Correspondence to Sean D Schoville.

Additional information

Competing interests

The authors’ have no competing interests to declare.

Authors’ contributions

SDS carried conducted the heat shock treatments, carried out the transcriptome analysis, and drafted the manuscript. FSB assisted in the transcriptome assembly, qPCR analysis and helped to draft the manuscript. GWM prepared the transcriptome samples. AW sequenced the HSP genes and assisted with qPCR analysis. RSB initiated the design and coordination of the study, obtained funding, and helped to draft the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and Permissions

About this article

Cite this article

Schoville, S.D., Barreto, F.S., Moy, G.W. et al. Investigating the molecular basis of local adaptation to thermal stress: population differences in gene expression across the transcriptome of the copepod Tigriopus californicus. BMC Evol Biol 12, 170 (2012).

Download citation


  • Next-generation sequencing
  • RNA-seq
  • Heat shock proteins
  • Diversification
  • Temperature stress
  • Climate