Locally adapted populations of a copepod can evolve different gene expression patterns under the same environmental pressures

Abstract As populations diverge in allopatry, but under similar thermal conditions, do similar thermal performance phenotypes evolve by maintaining similar gene expression patterns, or does genetic divergence lead to divergent patterns of gene expression between these populations? We used genetically divergent populations of the copepod Tigriopus californicus, whose performance at different thermal conditions is well characterized, to investigate transcriptome‐wide expression responses under two different thermal regimes: (1) a nonvariable temperature regime and (2) a regime with variable temperature. Our results show the expression profiles of the response to these regimes differed substantially among populations, even for populations that are geographically close. This pattern was accentuated when populations were raised in the variable temperature environment. Less heat‐tolerant populations mounted strong but divergent responses to the different thermal regimes, with a large heat‐shock response observed in one population, and an apparent reduction in the expression of genes involved in basic cellular processes in the other. Our results suggest that as populations diverge in allopatry, they may evolve starkly different responses to changes in temperature, at the gene expression level, while maintaining similar thermal performance phenotypes.

these genes are similar across studies, there is a lot of variation in the set of genes that respond to thermal stress. This is in part because of differences in study design, but also because these divergent taxa have evolved different ways to deal with thermal changes. Therefore, we were interested in determining how much the transcriptome-wide response varies when divergent populations of the same species are compared under different thermal regimes. As populations diverge in similar environments but in allopatry, does selection favor similar genetic responses across the populations? Or do novel mutations that appear in each population contribute to increased divergence in their response to thermal changes?
The copepod Tigriopus californicus is an ideal species in which to test this, as it gives us the opportunity to compare a number of locally adapted populations across a broad latitudinal spectrum differing in thermal tolerance. Populations of this species inhabit splash pools on the Pacific coast of North America and can be exposed to large changes in temperature on a daily basis. Gene flow is extremely low, which has led to the formation of many genetically divergent populations, even for populations that are only a few kilometers apart (Burton, 1997;Willett & Ladner, 2009). Within California, northern and southern populations fall out into two clades, with the northern clade encompassing populations north to Alaska (Burton, 1998;Edmands, 2001;Willett & Ladner, 2009;Figure 1a). Populations from the two clades show significant differences in local adaptation, particularly for thermal tolerance where upper thermal tolerance increases as latitude decreases, with lethal temperatures ranging from 35° in the northern clade to 38° in the southern clade (Kelly, Sanford, & Grosberg, 2012;Pereira, Barreto, & Burton, 2014;Tangwancharoen & Burton, 2014;Willett, 2010). In a study of acute thermal stress (35°), Schoville et al. (2012) found that T. californicus from a southern clade population (San Diego, CA) showed much greater upregulation of genes that are known to respond to heat stress, such as heat-shock proteins, than did individuals from a northern clade population (Santa Cruz, CA), suggesting that their ability to upregulate these genes may be at least in part responsible for their higher heat tolerance. The northern population differentially expressed a much higher number of genes overall (both up-and downregulated), but to much lower levels of fold change.
Differences in thermal adaptation in this species can also be observed as a shift in the thermal performance curve (TPC) of populations from south to north (Hong & Shurin, 2015). The width of the TPC remains the same, but the fitness peak shifts from warmer to cooler temperatures as you move north. In agreement with this, Willett (2010) showed that at higher nonlethal daily variable temperatures (20°-28°), two different southern populations outcompete two other northern ones, and the opposite is true for colder variable temperatures (16°-25°) or lower nonvariable temperatures (16°). The populations have roughly equivalent fitness at 20° (Willett, 2010). Both southern and northern populations commonly experience temperatures to and beyond 28° in nature, and when in the high variable temperature environment by themselves, all populations are able to develop and have offspring normally (Willett, 2010). Therefore, it seems northern populations are unable to survive in these high temperatures because of low competitiveness compared to the more heat-tolerant southern populations.
In contrast to numerous studies carried out at high, nearly lethal stressful temperatures, fewer studies of genetic response have been carried out under moderate and variable temperatures that many organisms may encounter on a regular basis (but see Barshis et al., 2013;Kenkel, Meyer, & Matz, 2013;Franssen et al., 2014;Dayan et al., 2015;Kenkel & Matz, 2016). Here, we look at transcriptome-wide gene expression under such moderately stressful temperatures in different populations of T. californicus. This study has two aims: (1) to determine how similar the transcriptome-wide response is between relatively geographically proximate populations with similar thermal phenotypes (the two northern and two southern populations used in Willett, 2010) and (2) to elucidate the molecular responses that underlie the ability of copepods from these southern populations to outcompete the northern ones at higher, variable temperatures but not at lower temperatures (Willett, 2010). We use two different thermal regimes, with a total of four treatments, and compare the transcriptome-wide response between the four populations from Willett (2010). By comparing the response from each population, we get a picture of how each population has adapted to deal with these temperature changes; while comparisons between treatments within populations show us the plastic response, these populations are able to mount given changes in temperature. The combined results show that the transcriptome-wide response from each population differs greatly even between populations that are only 8 km apart (the two southern populations) and is even greater between the northern populations. The difference in gene expression changes also increased when organisms were raised in the variable temperature regime, compared to a nonvariable regime followed by a one-time moderate heat stress. Our results suggest that the less heat-tolerant northern populations are potentially outcompeted by southern ones at least in part due to the costs associated with changing the expression of a much higher number of genes when in the moderately high variable temperature regime. Aquarium Systems) and consumed both commercial fish food and natural algae growth. Cultures were kept in incubators with 12-hr light:dark cycle at 20° for two to three generations before beginning the experiment. Copepods were exposed to two thermal regimes: (1) a nonvariable 20° and (2) a variable environment with 12 hr at 20° and 12 hr at 28° each day. Although the variable environment is more similar to the natural environment where daily temperature fluctuations are common, aspects such as rate of temperature change may be less representative (but are identical to those used in the Willett (2010) study). For both regimes for all populations, 25 gravid females were transferred to petri dishes (four dishes per population per treatment) and randomly assigned to one of the thermal regimes. Mothers were removed when offspring reached the copepodid stage. Therefore, copepods were born in their respective thermal regimes and stayed in this regime until they reached the adult life stage (copepodid stage 6) at which point RNA was extract for the different treatments.

| Copepod collection, culture, and thermal treatments
Differences in gene expression for each population, between thermal regimes, should be due to plasticity differences and not adaptation, as F I G U R E 1 Phylogeography of populations and experimental design. (a) Phylogeny and sampling locations for the four populations in this study. Branch lengths reflect average genetic divergence. Numbers below branches are bootstrap support. Percentages in parentheses are averages of uncorrected divergence based on 11,560 nuclear loci (figure modified from Pereira et al., 2016); (b) populations were exposed to two thermal regimes (nonvariable and variable). Gene expression was assayed from both regimes at 20° at the end of the 20° portion of the variable regime. Copepods from the nonvariable regime were moved to 28° at this point. RNA was isolated from both regimes at 28°, after two hours at this temperature. Dashed arrows indicate pairwise comparisons that were made to calculate relative gene expression between treatments. In red: nonvariable at 20° (20NV) compared to stress at 28° (28ST); in blue: variable at 20° (20V) compared to variable at 28° (28V); in dark gray: 20NV compared to 20V the copepods were in the different regimes only for the one generation. On the day before RNA isolation, 100 adult copepods (50 males and 50 females) were transferred to a new petri dish, without additional food. Each petri dish with 100 individuals was a replicate for one of the treatments. For each treatment, all populations had two dishes with 100 copepods and a separate mRNA library was created from the pool of 100 copepods for each population/treatment replicate (Fig. S1). As shown in Figure 1b, there were two sampling points for each thermal regime: (1) At the end of the 20° period in the variable regime, RNA was extracted from copepods from both regimes (20NV and 20V); at this point, copepods from the nonvariable regime were transferred to 28°; (2) 2 hr after the temperature shift to 28° on the variable regime, RNA was extracted from copepods at 28° from both regimes (28ST and 28V) (Figures 1b; S1.).

| RNA extraction and Illumina sequencing
For each treatment at the appropriate sampling time, copepods from each dish (~100 individuals) were collected onto filter paper, rinsed with filtered seawater, and immediately transferred to 250 μl of TRI Reagent (Sigma). RNA isolation was performed following the manu-

| Transcriptome assembly and annotation
An annotated reference genome for the SD-S population has been published online (https://i5k.nal.usda.gov/Tigriopus_californicus), and it was used to build transcriptome references for all populations in this study as follows. Annotated genes from the SD-S reference genome were extracted, and alternative splice variants were removed, keeping only the longest splice variant of a gene. To avoid a bias in mapping of reads from SD-S to its own high-quality reference with reads from the other populations likely to have some SNP differences, we built references from each of the four populations by mapping them to the SD-S genes. Trimmed reads from all treatments were combined for each single population and mapped to the SD-S genes using BWA MEM using default parameters (Li & Durbin, 2009). Following mapping, reads with low mapping quality (MAPQ < 20, likely the result of incorrect alignment) were removed, and the consensus sequences for each population were extracted using SAMtools and BCFtools (Li, 2011;Li et al., 2009). We assessed the quality of the references by calculating the total percentage "N" in the reference transcript set of each population. We also compared the sizes of the orthologous contigs and the percentage "N" in each contig for all transcripts that occurred in all four populations. As our gene expression analysis only compared within-population differences in treatment reads mapped to their own reference, no further action was taken at this stage of the analysis to remove short transcripts with a high percentage of "N," as these were removed during gene expression analysis.
The current SD-S genome assembly is thought to include ~80% of the entire genome, and certain genes are known to not be in this reference. To avoid excluding these genes from our analysis, we used de novo transcriptome assemblies from Pereira et al. (2016) to complement the pool of genes in our analysis as the populations used in this study are the same as the ones used in Pereira et al. (2016). To avoid confusion, we will refer to the set of genes relating to the reference genome as genomic transcripts (GTs), and from the de novo transcriptome assembly as de novo transcripts (DNTs). Transcripts originating from the different references can be differentiated based on their ID, where GT IDs begin with "TCALIF," while DNT IDs begin with either "Contig" or "comp." BLAST2GO (Conesa et al., 2005) was used to annotate transcripts from the references. We used the SD-S reference for annotation as it had slightly higher quality than the others. BLASTX searches were performed against the "nr" NCBI protein database, retaining hits with E ≤ 10 −3 . Gene Ontology (GO) terms (Ashburner et al. 2000) were retrieved for contigs with positive BLAST hits, with an E ≤ 10 −3 .

| Mapping and identification of differentially expressed genes
Trimmed reads for each treatment were mapped to their respective population's GT references in CLC GW (mapping parameters: similarity fraction = 0.99; length fraction = 0.8; mismatch cost = 2; insertion cost = 3; deletion cost = 3). Unmapped reads were retained and mapped to the DNT references using the same mapping parameters. Only unique mapped reads from the two sets of mapping files were used in further analyses. Only orthologous genes that occurred in all four populations were considered. Orthologous genes from the GT reference already had the same ID, while those from the DNT reference were extracted from the orthologous list in Pereira et al. (2016) and were given the SD-S ID for ease of comparison in later steps.
Differential expression was determined using the Bioconductor package edgeR (Robinson, McCarthy, & Smyth, 2010), in pairwise comparisons between treatments within each population. As the pairwise comparisons were not simply "control versus treatment," separate files for each pairwise comparison for each population were created and analyzed separately. Genes with very low levels of expression in both treatments were removed by retaining only those that accumulated at least two counts per million in at least two of four samples, allowing for both replicates of a treatment to have zero mapped reads, if both replicates in the other treatments had enough reads mapped to them. This also removed all genes that had zero mapped reads in the DNT reference, because the appropriate reads for these transcripts had already mapped to the GT reference.
Compositional differences between the libraries were normalized using the trimmed mean of M-value method (Robinson & Oshlack, 2010). For each population, all treatments were examined for possible batch effects between the replicates, using a multidimensional scaling (MDS) plot to visualize the level of similarity between each RNA sample for each population (Fig. S2). To account for batch effects between the replicates, a negative binomial generalized linear model (GLM) was fit, where "sequencing lane" was used as a blocking factor.
Likelihood ratio tests were used to determine genewise expression differences between two treatments, followed by a false discovery rate (FDR) correction of p values set to 5% (Benjamini & Hochberg, 1995). Genes that were detected as differentially expressed (DE), but had 0 reads mapped to both replicates in one of the treatments, were considered DE but were excluded in comparisons of the magnitude of their relative expression, as they showed abnormally high levels of fold change. Gene expression data are available through the NCBI Gene Expression Omnibus under Accession No. GSE80737.

| Pairwise comparisons of relative gene expression
Three pairwise comparisons between treatments were performed in EdgeR, each aimed at answering a specific question ( Figure 1b). To determine how populations respond to potentially low levels of thermal stress, relative expression between the nonvariable 20° (20NV) and stress 28° (28ST) treatments was calculated. To determine how populations respond to daily fluctuations in temperature, relative expression between the variable 20° (20V) and variable 28° (28V) was calculated. Differences between these two comparisons (28ST vs. 20NV vs. 28V vs. 20V) should then reflect the effects of phenotypic plasticity responses to exposure to 28°. To determine whether genes were being differentially expressed between variable and nonvariable regimes, relative expression between 20NV and 20V was calculated.

| Gene ontology enrichment analysis
Enrichment of Gene Ontology (GO) terms was assessed in Blast2GO using Fisher's exact tests for every GO term that appeared in a subset of genes, compared to all genes used in the analysis for each population (reference set). p-values were adjusted for multiple comparisons using a FDR set at 5%. GO enrichment analysis was performed for up-and downregulated genes separately for each pairwise comparison in each population. In many cases, redundant GO terms were significantly enriched in a dataset, in which case only one of these terms was included in our results using the following criteria: (1) If the terms associated with the same transcript belonged to different GO categories (i.e., biological process (P), molecular function (F), or cellular component (C)), they were all kept; (2) the most specific term was kept if it included approximately the same number of genes as the more broad term (e.g., response to stress (more specific) had 13 genes, response to stimulus (broader) also had 13 genes, response to stress was kept); (3) when multiple more specific terms combined included most of the genes in their common broader term, the specific ones were kept.

| Illumina sequencing and RNA-seq mapping
RNA sequencing yielded ~7.6-29 million reads per sample after trimming (Table S1) The MDS plot showed that all RNA samples clustered strongly by population (Fig. S2a). Within each population, treatments and replicates showed some separation based on the sequencing lane they were in, but there was also a lack of consistent grouping between replicates of the same treatment ( Fig. S2b-e). A likely explanation for this is the fact that our treatment temperatures were not extreme, 28° is only moderately stressful to these populations, and are not expected to cause large consistent changes in expression for the large majority of the transcriptome. RNA was extracted from a pool of 100 individuals for each treatment, and only a small percentage of all the transcriptome is responding in the same way across all tissues of all individuals. However, this should not be seen as negative outcome; as suggested by Sorensen et al. (2005), this averaging across tissues is an advantage as we are able to identify general responses taking place in several tissues and across several individuals. This should increase our power to detect genes that are globally important for thermal regulation.

| Response to moderate thermal stress in the two thermal regimes
When the different populations experience 28°, we expect northern populations to be more stressed as they have lower heat tolerance and their thermal performance curve is shifted to colder temperatures. Therefore, in both thermal regimes [nonvariable (NV) and variable (V)], we hypothesized that the northern populations will have to mount a stronger heat-shock response, which should translate into a larger number of differentially expressed (DE) genes, as well as higher magnitude of change for these genes. We expect, especially in northern populations, upregulated genes to be involved in the heatshock response pathway (such as HSP genes), but not necessarily upregulation of genes associated with proteolysis or apoptosis as the high-temperature treatments are well below the populations' lethal temperatures. We also expect to see the downregulation of genes associated with normal cell maintenance as the heat-shock response is known to favor the synthesis of HSPs while suppressing the expression of other genes (Tomanek, 2010). Southern and northern populations are expected to share more DE genes when compared to the population in the same region, as opposed to comparisons between regions. Our results, however, did not support (or only partially supported) our expectations. Below we detail the results for each thermal regime comparison.

| Response to moderate thermal stress after nonvariable temperature regime (28ST vs. 20NV)
For the northern populations, BB-N differentially expressed a higher number of genes than the two southern populations (204 up (Table S2).
In all four populations, GO enrichment analysis found that terms associated with response to stress, unfolded protein binding, and protein folding (except for BB-N) were enriched in DE genes that were upregulated (Table 1). Heat-shock protein genes contributed to the enrichment of these terms and were the most common class of upregulated genes (16 in SD-S, 13 in BR-S, 14 in SC-N, and 15 in BB-N). GO terms related to protein targeting and cell motility were also enriched in SD-S, while structural molecule activity and extracellular region terms were enriched in downregulated genes in BB-S. In this last case, these terms are overrepresented in large part due to the large number of cuticle proteins that are DE (especially in BB-N; Table 1).

| Response to moderate thermal stress in variable versus nonvariable temperature regimes (28ST vs. 20NV compared to 28V vs. 20V)
Next, we looked at the effects that raising individuals in a variable temperature regime (20°-28°) would have compared to the effects of raising them at a nonvariable 20° and exposing them to 28° as adults.
Between 28V and 20V, BB-N differentially expressed a larger number of genes than the southern populations (107 up, 384 down versus 94 F I G U R E 3 Percentage of genes with same versus opposite fold change direction in pairwise comparisons. (a) Stress at 28° (28ST) versus nonvariable 20° (20NV); (b) variable 28° (28V) versus variable 20° (20V). Fold change direction (but not magnitude) was determined for genes that were differentially expressed in one or both populations for all six pairwise comparisons. In both (a) and (b), within-region comparisons do not have a higher percentage of genes changing in the same direction than the betweenregion comparisons. Percentage of genes changing in the same direction is lower in the variable temperature regime for all comparisons. "Same (DE in both)" are genes that are differentially expressed in both populations, and fold change is in the same direction. "Same (DE in one)" are genes that are differentially expressed in only one of the two populations, but fold change is in the same direction in both. "Opposite" are genes with fold change in opposite direction in each of the populations T A B L E 1 Enriched gene ontology (GO) terms for differentially expressed genes in the three treatment comparisons  (Figure 2b and Appendix S1). Only shared 10 genes were DE in all four populations in this comparison (Table S2) (Figure 4; Appendix S1). Even though SC-N had a much smaller number of DE genes than any of the other populations, the magnitude of the potentially plastic gene expression change was high for several genes. For example, fold change between three HSP 70 genes (TCALIF_06728, TCALIF_04517, comp38417) decreased from 54.94, 24.16, and 22.90, to 14.18, 12.69, and 12.55, respectively (Appendix S1). However, even though the fold change decreased in the variable regime, fold change of these HSP genes is the highest  (Table 1). As mentioned before, cuticle protein genes contribute to this enrichment.

| Differential expression between thermal regimes (20V versus 20NV)
Lastly  Gray lines separate the four quadrants for up-and downregulated genes in the two comparisons were enriched in downregulated genes in these populations. In BR-S, extracellular region, carbohydrate metabolic process, and hydrolase activity terms were enriched in upregulated genes, and peptidase activity terms are enriched in downregulated genes. In SC-N, eight terms were enriched for upregulated genes, while nine were enriched for downregulated genes (Table 1). Several of these enriched GO terms in downregulated genes indicate that SC-N has to downregulate a number of genes that are involved in basic cell maintenance processes.
In BB-N, genes associated with structural molecule activity are overrepresented in upregulated genes, while three catalytic activity terms, as well as, extracellular region, and carbohydrate metabolic process genes are enriched in downregulated genes (Table 1).
One reason for genes to be upregulated at 20V compared to 20NV could be in "anticipation" of the higher temperatures individuals in the variable regime experience daily compared to those in the nonvariable regime (frontloading). In this study, these frontloaded genes would be genes that were upregulated in 28ST versus 20NV as well as in 20V versus 20NV. SD-S had 18 frontloaded genes (including six HSPs and one cuticle protein), BR-S had 10 frontloaded genes (including one small HSP and two cuticle proteins), SC-N had four frontloaded genes (three HSPs), and BB-N had 84 frontloaded genes (including eight HSPs and seven cuticle proteins) (Appendix S1).

| DISCUSSION
We used RNA-seq to determine transcriptome-wide patterns of gene expression for locally adapted populations of the copepod T. californicus in two thermal regimes. The higher temperature experienced by these populations (28°) should present a moderate stress, but all of the populations encounter this temperature in nature (Kelly et al., 2012), and it is well below their acute lethal temperature (Kelly et al., 2012;Pereira et al., 2014;Tangwancharoen & Burton, 2014;Willett, 2010).
Unlike studies that expose organisms to an acute thermal stress, the changes in gene expression we expect to observe are not the maximum response these organisms can mount, but the level of their response (both as the number of DE genes, as well as the magnitude of the expression change) should indicate the level of thermal stress each population is experiencing. Our results can also give us insights into the differences in fitness that have been observed between these populations under these specific temperature regimes (Willett, 2010).

| Transcriptome-wide response to moderate thermal stress differs between populations
The most striking result in this study was the level of differentiation in gene expression observed between the populations when they were exposed to a moderate heat stress (28°). This was particularly surprising for the two southern populations (SD-S and BR-S), which are separated by only 8 km and should share roughly similar thermal environments overall. There are, however, aspects of the biology of these copepods that can explain our results. There is extremely low gene flow between these two populations (Burton, 1997;Willett & Ladner, 2009), and shared polymorphism is also very low. Pereira et al. (2016) showed that shared polymorphism rapidly decreases as divergence increases between populations of this species. Mitochondrial DNA divergence is ~10% between SD-S and BR-S, and it can be >20% between the southern and northern populations used in this study (Burton, 1998;Willett & Ladner, 2009). Our results may also help partially explain the observed transgressive segregation that has been shown in hybrids between the SD-S and BR-S populations (Pereira et al., 2014). Late-generation hybrids (F9) between these two populations have higher thermal tolerance than either parental population, and some hybrid lines were even able to survive temperatures that are lethal to both parental populations. The results presented would suggest that this increased thermal tolerance of hybrids could be due to these two populations having evolved different ways to deal with increases in temperatures. When the two genomes are combined in hybrids, complementary gene action of factors that are associated with higher heat tolerance could then lead to transgressive, higher tolerance in these hybrids; however, this remains to be tested.
When we consider the two northern populations, where differences in gene expression profiles were even more dramatic, a combination of differences in abiotic selective pressures and population history is likely to play a more important role than between the southern populations. These northern populations are approximately 240 km apart and are twice as genetically divergent as the southern populations ( Figure 1a; Pereira et al., 2016). It is important to remember, however, that these two populations have similar upper thermal limits (Willett, 2010), have similar thermal performance curves (Hong & Shurin, 2015), and are both outcompeted by the southern populations when raised in variable 20°-28° environment (Willett, 2010).
Therefore, these thermal phenotypes might have different genetic underpinnings between these populations. While it is not uncommon for species or populations that differ in thermal tolerance to show different sets of DE genes with small overlap between them ( Barshis et al., 2013;Dayan et al., 2015;Franssen et al., 2014;Gleason & Burton, 2015;Narum & Campbell, 2015;Wang et al., 2014), our results are unique as we observe large differences between populations of the same species with similar thermal tolerance phenotypes. Willett (2010) (Table 1; Appendix S1). Therefore, it is possible that SC-N may have lower fitness at these variable temperatures not solely due to their level of heat-shock response when exposed to moderately stressful temperatures on a regular basis, as appears to be the case for BB-N, but also potentially because of how much it has to change its metabolic framework as a whole when it lives in these variable temperatures.

| General patterns of gene expression across all populations
It is interesting that even at this moderately high temperature (28°), all populations display some level of heat-shock response, upregulating a number of HSP genes and other genes associated with protein biding and refolding (Table 1, Appendix S1). This suggests that all populations, but especially the southern ones, often have to elicit a heat-shock response during the warmer months in nature (Kelly et al., 2012). The production of HSPs, and their dependence on ATP for function (HSP 70 and 90), can add considerable ATP demand to the cell and negatively affect the organism (Clare & Saibil, 2013;Feder et al., 1992;King & MacRae, 2015;Tomanek, 2010). Therefore, even though this heat-shock response is likely an adaptation that allows these organisms to persist when exposed to these temperatures, increases in mean temperature due to climate change may have a strong negative impact in this species. The southern populations already experience these temperatures (or higher) often in nature, and as is the case for many other intertidal species (Tomanek, 2010), increases in their ambient temperature would mean they have to mount this heatshock response for a higher proportion of time, possibly leading to decreased fitness.
As was observed in a previous RNA-seq study in T. californicus (Schoville et al., 2012), a number of genes that are annotated as cuticle protein genes are DE in all populations and are overrepresented in at least one of the treatments comparisons, although the direction of the expression change differs between the populations (Table 1).
While we do not know the exact function these genes/proteins have in response to thermal stress, cuticle protein genes have been observed to be DE in studies of thermal adaptation in Drosophila (Zhao et al., 2015), and a gene homologous to a DE T. californicus cuticle protein gene (Contig_59_58) in Anopheles gambiae (agap006369-pa) is annotated with a GO term associated with stress response. This gene occurs in A. gambiae within an inversion that contains a large cluster of cuticle protein genes, as well as three hsp83 genes (White et al., 2007); however, it is still unknown the role these cuticle protein genes play in heat or desiccation resistance in this mosquito. While we do not know the function these genes are serving, it is possible that cuticleassociated proteins are part of the thermal response in arthropods.
Studies in Chlorostoma snails, a species of Acropora corals, and two species of seagrass show that more thermally tolerant populations or species in these groups have higher constitutive levels of HSP gene expression (frontloading), which may enable them to more readily respond to thermal stress (Barshis et al., 2013;Dong et al., 2010;Franssen et al., 2014;Gleason & Burton, 2015;Tomanek & Somero, 1999). The same may be expected in individuals that were raised in a variable environment compared to a nonvariable one, where some genes that respond to thermal stress in the nonvariable temperature environment (28ST versus 20NV) would also be upregulated at 20° in the variable environment compared to 20° in the nonvariable environment (20NV-20V). Frontloading genes can be seen as a form of plastic response to higher daily temperatures, and while all populations have at least some frontloaded genes, BB-N has by far the most, including not only several HSP genes, but also several cuticle protein and lectin genes (Appendix S1). SC-N has a small number of upregulated genes at 20NV-28ST, which limits the number of frontloaded genes it can have in this case, but of the four genes that are frontloaded, three are HSP genes (Appendix S1). The two southern populations have intermediate numbers of frontloaded genes (18 for SD-S and 10 for BR-S), again suggesting that the two northern populations mount very different molecular responses to the variable temperature environment.
The present study determined a small set of candidate genes that are likely to be the product of adaptive plastic response to the different temperature treatments. These are the genes that are DE in all populations in the different treatment comparisons: 25 DE genes for 28ST versus 20NV, and 10 for 28V versus 20V (seven common to both comparisons; Table S2). The large majority of these genes are associated with heat-shock response, with a strong overrepresentation of HSP genes (Table S2). One of these genes codes for hsp beta-1, a gene on which Barreto, Schoville, and Burton (2015) performed RNAi using T. californicus, and showed that individual copepods that had this gene knocked down had lower thermal tolerance compared to control individuals. In the present study, this gene is not only upregulated in all four populations at 28° in both thermal regimes, but it also shows strong signs of plasticity with a dampening of gene expression in comparisons of nonvariable and variable temperature environment in all but one population; SC-N has very similar levels of fold change between 28° and 20° in both thermal regimes. While we do not have direct evidence that the other genes in this group are the result of adaptive plastic response, they are strong candidates for further work into their function in thermal stress response.

| CONCLUSION
The present study highlights some key ways that local adaptation can impact the manner in which an organism deals with temperature changes at the level of gene expression. First, it is clear that locally adapted populations of the same species can display different responses at the molecular level despite showing similar thermal performance phenotypes. This was particularly striking between SD-S and BR-S, which are only 8 km apart and yet show substantial differences in their responses to both temperature variability and moderate heat stress. Therefore, even for closely related populations, the molecular mechanisms they use to deal with temperature changes can be different. Second, as seen in the northern populations, the molecular response to changes in thermal variability (i.e., 20NV vs. 20V comparison) may be drastically different, even though their upper thermal limit is very similar. Therefore, we may be underestimating the amount of variation to stress response at the molecular, and potentially physiological, level in species with very segregated populations, by assuming that populations with similar thermal tolerances respond to changes in temperature in the same way.