The thermal dependence and molecular basis of physiological color change in Takydromus septentrionalis (Lacertidae)

Color change is mediated by synchronous intracellular transport of pigmented organelles in chromatophores and has functional roles in signaling (visual communication), background matching (camouflage), and thermoregulation (Baling et al., 2016; Smith et al., 2016a,b). Previous studies have focused more intensively on the former two aspects than on the latter one. Thermoregulation is important for animals because body temperature affects numerous physiological and behavioral performances (Angilletta, 2001). Animals regulate body temperature by a combination of behavioral and physiological mechanisms. Body temperature is dependent on environmental temperature in ectotherms, and this is especially true for those in thermally uniformed environments where behavioral thermoregulation is constrained (Angilletta, 2001; Muri et al., 2015). Many animals are vulnerable to climate warming primarily because the increased heat loads constrain the activities necessary to maintain a normal life (Bakken, 1992; Huey et al., 2009, 2010; Sinervo et al., 2010; Wang et al., 2017). One way of thermoregulation involves a rapid change in coloration, which alters absorption and reflection of radiant heat (Stuart-Fox and Moussalli, 2009; Smith et al., 2016b). There are two major types of color change: morphological color change occurs over days or months (Aspengren et al., 2009); physiological color change occurs much more rapidly in a few seconds or minutes (Fujii, 1993). Chromatophores originate from the embryonic nerve crest and form units of migrating dermal chromatophores on the outer layer of the epidermis (Bagnara and Hadley, 1969; Fujii, 1993). The density and aggregation of pigments in chromatophores determine skin color (Grether et al., 2004). Morphological color is determined by the number and quality of chromatophores in the dermis, and physiological color by the relative state of pigment dispersion and aggregation within each chromatophore (Ligon and McCartney, 2016).

Physiological color change occurs in response to environmental signals and is mediated by nervous-fibers or the endocrine system (Fujii, 1993; Ligon and McCartney, 2016). In fish and reptiles, the sympathetic nervous system responds to stress and releases catecholamine (Summers and Greenberg, 1994), which binds to adrenoceptors on the plasma membrane and increases the concentration of secondary messengers such as cyclic adenosine monophosphate (cAMP) (Nery and Castrucci, 1997). High concentrations of cAMP activate protein kinase-A, leading to protein phosphorylation (Nery and Castrucci, 1997). Phosphorylated proteins detach from the cytoskeleton and move along cytoskeletal filaments to the cell periphery (Nery and Castrucci, 1997). Microtubule- and microfilament-dependent pigment transportation requires the assistance of motor proteins such as kinesin, dynein and myosin (Rogers and Gelfand, 1988; Ligon and McCartney, 2016). The genes MYO3 and MYO5 encode classical myosin I proteins (Goodson et al., 1996). Noradrenaline, a strong pigment-aggregating hormone, can activate cAMP synthetase adenylate cyclase and stimulate melanosome diffusion in fish and amphibians (Nilsson et al., 2013). Melanocyte-stimulating hormone (MSH) and melanocyte-concentrating hormone (MCH) regulate body color change (Sugimoto, 2002; Mizusawa et al., 2013; Nilsson et al., 2013). The MSH and MCH receptors, MC1R and MCHR, both participate in color change in fish, amphibians, birds and mammals (Nilsson et al., 2013; Bertolesi and McFarlane, 2018). Opn4x in amphibians and the orthologous Opn4m in mammals regulate α-MSH levels though a neuroendocrine circuit (Bertolesi and McFarlane, 2018).

Many lizards have the ability to change body color. For example, bearded dragons (Pogona vitticeps) darken or lighten their color, thereby altering the amount of solar radiation absorbed or reflected (Stuart-Fox and Moussalli, 2009; Nilsson et al., 2013; Smith et al., 2016b). Solar radiation comprises ultraviolet light (UV, 290–400 nm), visible light (Vis, 400–700 nm) and near infrared light (NIR, 700–2600 nm) (Porter, 1967). Most energy exists within the Vis and NIR parts of the spectrum; however, NIR has little to no effect on camouflage because of the insensitivity of the visual system in animals to NIR (Clusella-Trullas et al., 2007). UV-Vis is the most important spectrum component affecting absorption of solar radiation in P. vitticeps (Smith et al., 2016b). As the dorsal surface is able to reflect light within the Vis wavelengths in lizards (Murillo et al., 2020), one may hypothesize that Vis is the most important part of solar radiation for thermoregulation. Here, we quantified dorsal coloration based on the RGB value in adult male northern grass lizards (Takydromus septentrionalis) and then examined the transcriptomic basis of that color change to test the hypothesis. The RGB color model is an additive color model used for numerical color specifications. The RGB value refers to luminance, which is the most important color component affecting absorption of solar radiation; the RGB value can be transformed to YCbCr color spaces, in which Y represents the luminance component (Ahirwal et al., 2007; Smith et al., 2016a).

Temperature affected body color luminance (F_{6, 342}=6.76, P<0.001), which generally decreased with increasing temperature at lower temperatures (from 24°C to Tp) and increased with increasing temperature at higher temperatures (from Tp to 36°C) (Fig. 1a). Luminance at temperatures close to Tp was lowest in all three populations (Fig. 1a). Y values varied among populations (F_{2, 57}=5.668, P<0.01), with the mean value being greatest in the ND population and smallest in the LS population (Fig. 1a; Tukey's test, all P<0.05). The temperature×population interaction was not a significant source of variation in luminance (F_{12, 342}=1.43, P=0.150).

In both TG1 and TG2 and for each population, mean values for luminance differed between lizards measured at 24°C or 36°C and at Tp (paired-sampled t-test, all P<0.001). The differences (ΔY) in Y values differed among the three populations (TG1: F_{2, 57}=5.782, TG2: F_{2, 57}=7.553, both P<0.01). In both TG1 and TG2, mean ΔY values were smaller in the LS population than in the other two populations (Fig. 1b). From May to September, patterns of monthly changes in mean temperature and solar radiation were similar among the three populations (Fig. 1c). Cluster analysis using 19 climate variables showed that the GY and ND populations formed a clade climatically differing from the LS population (Fig. 1d).

We de novo assembled a reference transcriptome because of the lack of genomic information on T. septentrionalis. A total of 487,838,568 raw reads were yielded from nine libraries. After adaptor trimming and quality filtering, 472,500,030 clean reads remained. A total of 230,477 unigenes were assembled with a mean length of 831 bp, and an N50 of 987 bp. Furthermore, ∼82% of the unigenes ranged from 300 to 1000 bp, ∼11% ranged from 1000 to 2000 bp, and ∼7% exceeded 2000 bp.

A total of 100,795 (43.7%) unigenes were annotated according to seven databases. Based on the Nr and Nt databases, 36,135 (15.7%) and 78,697 (34.1%) open reading frames were identified, respectively. In total, 26,458 (11.5%) protein-coding genes were annotated according to the Swiss-Prot database. Because of the lack of genomic information on Takydromus, top-hit species distribution revealed a close kinship between two lizard species, P. vitticeps and Anolis carolinensis, based on Nr annotation (Fig. S1a).

A total of 9748 unigenes were classified into 26 KOG categories; the largest group was ‘general function prediction only’ with 1625 unigenes, followed by ‘signal transduction mechanisms’ with 1568 unigenes (Fig. S1b). A total of 12,458 unigenes were mapped to 232 pathways based on KEGG Automatic Annotation Server. ‘Signal transduction’ was the most represented category, containing 1521 unigenes, followed by ‘endocrine system’ (761 unigenes), ‘immune system’ (738 unigenes), and ‘transport and catabolism’ (706 unigenes) (Fig. S1c). A total of 38,198 unigenes were classified to 1536 GO annotation; 24,042 unigenes were attributed to biological process, including 1049 GO terms; 15,712 unigenes were attributed to cellular component, including 323 terms; and 22,683 unigenes were attributed to molecular function, including 164 terms. The GO terms attributed to the greatest number of genes were ‘cellular process’, ‘function binding’, ‘metabolic process’, and ‘single-organism process’ (Fig. S1d).

PCA including all unigenes revealed differences in gene expression among three temperature treatments, with the first two components explaining 79.3% of the variance (Fig. S2). We found that highly expressed genes were from the COX, HBA, HBE, MYO, and KRT families. Of the genes affecting hair color in mammals (Hoekstra, 2006), 46 homologous unigenes were detected in T. septentrionalis, and 40 unigenes were expressed (mean FPKM>0.3 in least two temperature treatments), with four highly expressed (mean FPKM>15) and 36 lowly expressed (mean FPKM<15) in the skin.

A total of 2562 upregulated and 827 downregulated unigenes were identified when the 24°C and 30°C treatments were compared (24/30); 594 upregulated and 503 downregulated unigenes were detected when the 36°C and 30°C treatments were compared (36/30); 2483 and 794 DEGs were annotated in at least one database, respectively (Fig. 2a,b; Tables S2, S3). The results of qRT-PCR were consistent with those from RNA-seq. Six DEGs (two motor protein-encoding genes, two keratin encoding genes, and two transcription factors) were upregulated at 24°C and 36°C. Expression of the downregulated cAMP-responsive element modulator (CREM) was suppressed at 24°C and 36°C, and the expression of OCA1 was similar among the three treatments (Fig. 2c).

Among the DEGs, 384 were co-upregulated and 135 were co-downregulated at 24°C and 36°C, while four DEGs showed an opposite expression pattern (Fig. 2b). Based on the K-means cluster analysis, the expression of numerous genes followed the same pattern at 24°C and 36°C (Fig. 3, Fig. S3). Based on the annotation of DEGs, some kinase-related (Fig. 3a) and motor protein-related (Fig. 3b) genes, and cytoskeletal genes (Fig. 3c,d), were identified, with most highly expressed at 24°C and 36°C (Fig. 3).

The expression of CREM was lower at 24°C and 36°C than at 30°C (Fig. 2c). Hundreds of kinase-related genes were differentially expressed, with some involved in phosphorylation (Tables S2, S3); the expression of phosphorylation-related kinase genes was thermally sensitive. Compared with the 30°C treatment, 30 unigenes were upregulated and three were downregulated at 24°C, while nine unigenes were upregulated and four were downregulated at 36°C (Fig. 3a). Six phosphorylation-related kinase genes were co-expressed at 24°C and 36°C, with all upregulated.

The GO analysis of DEGs from 24/30 and 36/30 revealed that 23 and 22 enzyme activity GO categories were enriched, respectively (Table S3). Among them, ‘the acid phosphatase activity’ was found in 24/30 and 36/30, with 12 and nine DEGs, respectively (Table 1). Nineteen and 12 DEGs in the 24°C and 36°C treatments were involved in ‘phosphatase activity’. ‘G-protein coupled receptor activity’ was also enriched in 24/30 and 36/30, including 29 and 16 DEGs, respectively, among them 24 and 11 DGEs were upregulated (Table 1).

Hundreds of kinesin-, dynein-, and myosin-related genes were annotated, for instance, moy5, KRT1, and KIF. Compared with the 30°C treatment, 30 motor protein genes (15 myosin-, eight dynein-, and seven kinesin-related genes) were differentially expressed at 24°C, and 17 genes (seven myosin-, four dynein-, and six kinesin-related genes) were differentially expressed at 36°C (Fig. 3b). More than a half of motor protein related DEGs were highly expressed at 24°C and 36°C, with 23 (12 myosin-, six dynein-, and five kinesin-related genes) upregulated at 24°C, 12 (three myosin-, four dynein-, and five kinesin-related genes) upregulated at 36°C, and nine (two myosin-, three dynein-, and four kinesin-related genes) co-upregulated at the two temperatures.

Some GO terms were associated with motor function. More than ten DEGs were contributed to ‘dynein complex’, ‘dynein binding’, and ‘motor activity’ in 24/30 and 36/30. Protein accumulation in the nucleus of cells was highly enriched, 14 and eight DEGs were annotated to ‘protein import into nucleus’, ‘protein localization to nucleus’, and ‘protein targeting to nucleus’ in 24/30 and 36/30 respectively (Table 1).

Numerous tubulin and filament genes were classified, with many of them differentially expressed at different temperatures. Compared with those at 30°C, 14 tubulin-related (nine upregulated and five downregulated) and 32 filament-related (25 upregulated and seven downregulated) genes were detected at 24°C; and at 36°C ten tubulin-related (nine upregulated and one downregulated) and 25 filament-related (23 upregulated and two downregulated) genes were detected (Fig. 3c). Most of these genes were co-upregulated, including seven out of eight co-expressed tubulin-related genes and 17 out of 18 co-expressed filament-related genes (Fig. 3c; Table S2). Several other cytoskeleton related genes were differentially expressed, with most of them upregulated at 24°C and 36°C (Fig. 3d).

GO analysis of DEGs revealed many functions related to the cytoskeleton. Twenty-nine and 25 DEGs were related to ‘intermediate filament’, 27 and 21 DEGs were annotated to ‘microtubule cytoskeleton’, 56 and 41 DEGs were contributed to ‘cytoskeletal part’, and 12 and ten DEGs were related to ‘structural constituent of cytoskeleton’ in the 24°C and 36°C treatments, respectively (Tables 1, Table S3). Six DEGs were involved in ‘microtubule motor activity’ in the 24°C and 36°C treatments (Tables 1, Table S3).

MCHR1-1, MCHR1-2, MC1R and the OPN4x orthologue were identified in our data set. QRT-PCR revealed that the expression of these genes in pituitary glands differed among the three temperature treatments. Compared with the 30°C treatment, the expression of MC1R was upregulated and the expressions of MCHR1-2, pro-MCH, and OPN4x were downregulated at 24°C and 36°C, and the expression of MCHR1-1 was downregulated at 24°C and upregulated at 36°C (Fig. 4).

Several DEGs involved in disease-related pathways (Fig. S4), including ‘African trypanosomiasis’ and ‘amoebiasis’. In addition, at 24°C and 36°C, ‘cytokine-cytokine receptor interaction’, ‘proximal tubule bicarbonate reclamation’, ‘steroid biosynthesis’ and some metabolism-related pathways were classified.

Body color becomes dark when pigments are dispersed throughout the cell, and light when the pigments are aggregated (Aspengren et al., 2009; Nilsson et al., 2013; Smith et al., 2016b). Kirchhoff's law of thermal radiation predicts that absorption and release of radiation energy are positively correlated (Robitaille, 2003). Accordingly, individuals with a dark color absorb or release radiant heat more rapidly than those with a light color (Seebacher, 2000). Our data show that body color is thermally sensitive in T. septentrionalis, and that lizards quickly change body color in response to the thermal environment around them. From the pattern of temperature-dependent color change found in T. septentrionalis we may conclude that pigments are more dispersed at body temperatures around Tp than at other lower or higher levels. Body temperatures vary within a range narrower than that of natural thermal fluctuations experienced by reptiles including T. septentrionalis, as revealed by the fact that their body temperatures are higher than ambient temperatures in the low-temperature environment and lower than ambient temperatures in the high-temperature environment (Wang and Xu, 1987; Seebacher and Shine, 2004). Thus, to reduce heat loss (due to increased release of radiation energy) at low temperatures or to avoid overheating at high temperatures, lizards should maintain body color with a relatively high luminance. Our finding that rapid color change functions as a mechanism of thermoregulation is consistent with earlier studies on P. vitticeps (Cadena and Tattersall, 2009; Smith et al., 2016b). Moreover, regardless of the effect of temperature, the ranges of body color variation are finite in T. septentrionalis from the three populations and more similar between the GY and ND populations, which are climatically more similar than other two population pairs (Fig. 1d).

The genomes have been sequenced and characterized in several species of reptiles not including T. septentrionalis. Transcriptase analysis may provide useful information for understanding genetic information and molecular mechanisms. Currently only a few complete genome sequences of reptiles of the suborder Sauria have been available, including A. carolinensis (Alföldi et al., 2011), Shinisaurus crocodilurus (Gao et al., 2017) and P. vitticeps (Georges et al., 2015). The assembled genome of sequenced lizards contains almost 30,000 protein-coding genes (Alföldi et al., 2011; Georges et al., 2015; Gao et al., 2017). We obtained 230,477 unigenes from RNA-seq, far exceeding the number of genes in the species with known genomic information. Nevertheless, this is common in transcriptome analyses of animals (Brandley et al., 2012), due to complex and abundant alternative splicing (Pan et al., 2008; Black, 2003). In humans, transcripts from 95% of multiexon genes undergo alternative splicing, and there are around 100,000 intermediate-to-high abundance alternative splicing events in major human tissues (Pan et al., 2008). In P. vitticeps, 19,406 protein-coding genes were obtained through genome assembly, while 595,564 contigs were obtained with Trinity for transcriptome analysis (Georges et al., 2015). Therefore, alternative splicing may be common in T. septentrionalis.

OAC1 and P-protein participate in pigment biosynthesis (Brilliant, 2001). In the present study, there was no significant change in the expression of these homologous genes in the three constant-temperature treatments. These results indicated that the body color change we observed in T. septentrionalis was not a morphological but a physiological response. High levels of cAMP facilitate protein phosphorylation and promote pigment dispersal (Nery and Castrucci, 1997; Rodionov et al., 2003). While high concentration G protein induces pigment aggregation, because G protein-coupling system inhibits adenylate cyclase activity, which can transform ATP into cAMP and result in lowered cAMP levels (White et al., 1987; Nery and Castrucci, 1997; Aspengren et al., 2009). In the present study, high expression of G protein-related genes was accompanied by low CREM expression. For example, at 24°C and 36°C, G protein-coupled receptor genes were highly expressed, while CREM expression was lower. This finding suggests that T. septentrionalis shares the same G protein-cAMP signaling pathway with fish (Nery and Castrucci, 1997; Nilsson et al., 2013). Pigment aggregation/dispersion is microtubule-dependent, and the concerted efforts of both microtubule and actin-based motors are required for proper pigment distribution in melanophores (Aspengren et al., 2009; Bertolesi and McFarlane, 2018). It seems likely that expression oscillations of cytoskeletal and motor-related genes at different temperatures drive pigment movement and thereby cause body color change in T. septentrionalis.

MSH and MCH are commonly involved in physiological color change in fish (Suzuki et al., 1995; Sugimoto, 2002; Nilsson et al., 2013). In the barfin flounder, a-MSH-rp exhibits pigment-dispersing activities in vitro, but not in vivo, because of inhibition from the nervous system (Kobayashi et al., 2010; Mizusawa et al., 2013). In our study, however, MSH may influence the pigment-dispersing in vivo, because the expression of MC1R and OPN4x were both changed, making it possible that the function of MSH was influenced by nervous system. In teleost fish, MCH triggers color change in vivo and in vitro (Yamanome et al., 2007; Mizusawa et al., 2013). Moreover, the findings of our study indicate that the expression of MCH-related genes was altered; therefore, we may hypothesize that MCH also affects color change in vivo in reptiles. In teleost fish, MCH results in the aggregation of melanin granules in melanophores (White et al., 1987; Rodionov et al., 2003; Aspengren et al., 2009; Nilsson et al., 2013). While in frogs and lizards melanophores MCH induces pigment dispersion from the center of cells (Bertolesi and McFarlane, 2018). In the present study, the expression of MCHR1-2 and pro-MCH genes was significantly higher at 30°C than at 24°C and 36°C, thus resulting in pigment dispersion at 30°C. Clearly, as in frogs, MCH may have a melanin-dispersing effect in T. septentrionalis. In teleost fish, MCH directly binds to seven G protein-coupled transmembrane receptors, activates G proteins, and then causes pigment to aggregate rapidly (Kawauchi, 2006; Aspengren et al., 2009). However, a human and a mouse melanopsin variant (hOpn4L and mOpn4L) can activate or deactivate G-protein pathways by shifting blue/yellow light (Spoida et al., 2016). In T. septentrionalis, high expression of MCH-related genes was accompanied by lower expression of G protein-related genes. This finding suggests that MCH may deactivate G proteins, enhance cAMP concentration, and cause pigment dispersion in lizards.

In this study, we quantified dorsal coloration based on the RGB values in adult males of T. septentrionalis from three populations and exposed to different thermal conditions to test whether and, if so, how body color changes in response to the thermal environment around them. Our results showed that physiological color change was thermally sensitive in T. septentrionalis and had a functional role in thermoregulation achieved by altering the arrangement of pigments and thereby changing the absorption or reflection of solar radiation. Body color became increasingly light as body temperatures deviated from the level preferred by the lizard, and the thermal dependence of body color change differed among populations. From our analysis of the skin transcriptome in individuals exposed to three constant temperatures (24, 30 and 36°C) we know the following. First, G protein-coupling system inhibits adenylate cyclase activity to transform ATP into cAMP, causing the pigment to aggregate quickly. Second, MCH deactivates the G proteins and thereby initiates pigment dispersion in T. septentrionalis. Third, the molecular mechanism of physiological color change found in T. septentrionalis differs from that reported for teleost fish where the MCH activates the G proteins and thereby causes pigment aggregation. Future work could usefully investigate whether the molecular basis of the thermal dependence of physiological color change observed in T. septentrionalis is generalizable to lizards.

Sixty adult males, 20 from each of three populations in Lishui (LS; 28°26′N, 119°55′E), Ningde (ND; 26°53′N, 119°39′E) and Guiyang (GY; 25°41′N, 112°34′E), were collected by noose from the field in April 2017. Lizards were brought to Nanjing, where ten from the same population were housed in each 1.2×0.6×0.4 m cage containing 50 mm-depth moist soil, grasses and pieces of clay tile. Cages holding lizards from different populations were all placed in a room where temperatures varied from 20−28°C for at least two weeks prior to testing the thermal dependence of body color change, thereby removing or minimizing the effects of potential factors other than temperature. Mealworm larvae (Tenebrio molitor), house crickets (Achetus domestica) and water were provided daily, and thermoregulatory opportunities were provided during daytime hours (08:00 to 17:00 h) using a 100 W heating light. Our experimental procedures complied with the current laws of China for the care and use of experimental animals and were approved by the Animal Research Ethics Committee of Nanjing Normal University (IACUC-200411).

We first used a seven temperatures × three populations repeated-measures factorial design experiment to examine the effects of temperature, population and their interaction on body color. Lizards were maintained at one of seven temperatures in a randomized sequence (24, 26, 28, 30, 32, 34 and 36°C) for 30 min to adjust body temperature at the test level. Lizards were then photographed vertically (always by Kun Guo and from the same distance) with a digital camera (Canon, Japan) using manual settings and an in-built flash unit as the light source. A black scale was included in the picture for calibration.

We then followed the procedures described above to photograph lizards on two thermal gradients (TGs) created in the same room pre-set at different temperatures. TG1 was achieved by moving cages to the room pre-set at 24°C, where a thermal gradient from 24−50°C was created by a 100 W heating light (Yang et al., 2008). TG2 was achieved by moving cages to the room pre-set at 36°C, where a thermal gradient from 12−36°C was created by a 500 g ice cube at one end of each cage. Two photos (P) were taken for each lizard: P1 at 24°C (for the TG1 treatment) or 36°C (for the TG2 treatment); P2 at the preferred body temperature (Tp ≈ 30°C), which could be achieved after a lizard was on a thermal gradient for 30 min (Yang et al., 2008).

We used MATLAB to calculate RGB values, which were transformed to luminance values (Y) using an equation of Y=16+(0.257×R+0.504×G+0.09×B) (Ahirwal et al., 2007). We used the NicheMapR script in R to obtain climate date for the three populations from global monthly hourly solar radiation and shadowless air temperature data for 2013 from the microclimate model (Kearney et al., 2014).

Three GY males (5.63±0.38 g) were exposed to one of three temperatures (24, 30 and 36°C) for 30 min. The dorsal skin of each lizard killed by rapid freezing in liquid nitrogen was then collected and stored at −80°C for later total RNA isolation using the PureLink^® RNA Mini Kit (Invitrogen, USA). RNA quality detection, library preparation, and sequencing were performed according to Zhong et al. (2018). De novo transcriptome assembly of the clean reads was performed using Trinity 2.4.0 (Grabherr et al., 2011) with min_kmer_cov set to 3 by default. Hierarchical clustering was performed by Corset 1.05 to obtain the longest cluster sequence for subsequent analysis (Davidson and Oshlack, 2014). The final unigenes were annotated through NCBI non-redundant protein sequences (NR), euKaryotic Ortholog Groups (KOG), and Swiss-Prot, NCBI non-redundant nucleotide sequences (Nt), Protein family (Pfam), KEGG orthology database (KO) and Gene ontology (GO) as usual (Lan et al., 2018). We used edgeR 3.24.3 to analyze differentially expressed genes (DEGs) with absolute change ≥2 and FDR ≤0.05.

Skin and pituitary glands collected from lizards exposed to three temperatures were ground in liquid nitrogen for RNA extraction. Total RNA was isolated using TRIzol Reagent (Invitrogen, USA). Reverse transcription, qRT-PCR and data analysis were performed according to Zhong et al. (2018). Ubiquitin gene was used as the internal control. Table S1 shows primer sequences.

We used repeated-measures analysis of variance (ANOVA) to examine the effects of temperature, population and their interaction on body color. Tukey's tests were performed to examine whether the differences between temperature treatments or populations were significant, and P values were corrected for multiple tests using the false discovery rate (FDR). Differences in ΔY (the absolute difference in Y between P1 and P2) among populations were tested using analysis of covariance (ANCOVA). For all the unigenes and all samples, their expression values, in terms of fragments per kilobase million, were used in principal components analysis (PCA).

We are grateful to Zhi-Hua Lin, Xiao-Li Fan and Guo-Hua Ding for their assistance.