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

1Thermal 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. 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. 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 and among congeners. 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. 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. Additionally, rates of interpopulation gene flow must remain low, or gene swamping will override the effect of diversifying selection.

2Several empirical studies have found evidence of the genetic fixation of alternative alleles related to local adaptation of populations along latitudinal clines and altitudinal gradients. These studies suggest that local adaptation sometimes involves structural changes in single protein-coding loci with large effects. Point mutations in protein coding sequences can result in allelic gene products with differing thermal optima 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, particularly since a single change in a master regulator gene can result in widespread changes in gene expression across the transcriptome. 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.

3Until 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. 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. 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. Next-generation RNA sequencing technology (RNA-seq) 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, and identification of any structural variation among copies of ortholgous genes found in different populations.

4Here we initiate investigations of the genetic mechanisms that underlie population differences in thermal tolerance observed in Tigriopus californicus. Besides experiencing extreme daily and seasonal temperature fluctuations (annual range of 6-33°C at one site) in their local shallow intertidal pools, 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. Willett and Kelly et al. 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.

5At 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. 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.

6The 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. Each sample was sequenced for single-end 100 basepair read lengths.

7Reference transcriptomes for SD and SC were previously developed from 454 pyrosequencing of each population. 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.

8Assembled 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. 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) terms were retrieved at a threshold of 10-6.

9Orthology 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.

10Illumina 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. 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.

11As 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. 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. A neighbor-joining tree was constructed in Mega5, 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.

12In 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.

13We 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.

14The 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 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).

15In 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).

16We 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.

17Although 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).

18When 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.

19Average 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.

20We 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).

21Populations of Tigriopus californicus show strikingly different survivorship following acute thermal stress, 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. For example, even at the interspecific level, Lockwood et al. 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.

22We 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. 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. 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.

23In 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, 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, 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, where northern populations exhibit an up-regulation of proteolytic genes.

24In 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, other studies have shown that gene expression changes in experimental treatments are relevant predictors of organismal physiology, and that there is often a high correlation between gene expression and protein abundance, including the Hsp gene family.

25Numerous species exhibit morphological and life history variation along latitudinal clines. These traits typically involve changes in body size, phenology, or temperature adaptation. Recently, Willett and Kelly et al. 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. 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.

26Our 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. 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, 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, as well as responding to other forms of environmental stress. 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.

27In 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.

28The authors’ have no competing interests to declare.

29SDS 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.