Selection on adaptive and maladaptive gene expression plasticity during thermal adaptation to urban heat islands

Phenotypic plasticity enables a single genotype to produce multiple phenotypes in response to environmental variation. Plasticity may play a critical role in the colonization of novel environments, but its role in adaptive evolution is controversial. Here we suggest that rapid parallel regulatory adaptation of Anolis lizards to urban heat islands is due primarily to selection for reduced and/or reversed heat-induced plasticity that is maladaptive in urban thermal conditions. We identify evidence for polygenic selection across genes of the skeletal muscle transcriptome associated with heat tolerance. Forest lizards raised in common garden conditions exhibit heat-induced changes in expression of these genes that largely correlate with decreased heat tolerance, consistent with maladaptive regulatory response to high-temperature environments. In contrast, urban lizards display reduced gene expression plasticity after heat challenge in common garden and a significant increase in gene expression change that is congruent with greater heat tolerance, a putatively adaptive state in warmer urban environments. Genes displaying maladaptive heat-induced plasticity repeatedly show greater genetic divergence between urban and forest habitats than those displaying adaptive plasticity. These results highlight the role of selection against maladaptive regulatory plasticity during rapid adaptive modification of complex systems in the wild.

U nderstanding the mechanisms that generate novel phenotypic diversity is a central goal of evolutionary biology. A long-standing gap in this understanding centers around the interplay between environmentally induced phenotypic variation (phenotypic plasticity) and genetically based evolutionary change. The modern history of biology has seen plasticity viewed alternatively as constraining evolution by buffering genotypes against the full brunt of natural selection [1][2][3][4] , as inconsequential for evolutionary change as a source of non-heritable phenotypic variation [5][6][7] , and as facilitating colonization of new environments -in particular during the incipient stages of adaptation 4 . The circumstances under which plasticity assumes these various roles 8,9 and its role in adaptive lineage divergence 10-12 remain controversial.
The adaptive value of phenotypic plasticity depends on the heritability of environmental response and its direction of effect relative to the phenotypic optimum (Fig. 1). Adaptive plasticity moves individuals towards the phenotypic optimum in a given environment, enabling population persistence 6 and allowing selection to act upon genetic variation underpinning the inducible response. This effect moves the population mean closer to the local adaptive peak 4,[13][14][15] (Fig. 1). The evolution of adaptive plasticity has been demonstrated in a wide range of developmental 16 , morphological 15 , physiological 17 , life history 13 , and behavioral phenotypes 18 .
Conversely, maladaptive plasticity results from an ancestral response that is inappropriate in a novel environment, reducing fitness [19][20][21] . Selection on genetic variation underlying maladaptive plasticity can minimize its magnitude and/or reverse its direction entirely 4 (Fig. 1). For instance, ancestral physiological responses to low oxygen have evolved in lowland species to cope with acute endogenous hypoxia due to anemia, blood loss, or low tissue oxygenation 22 . However, at high altitudes, where low oxygen is exogenous and chronic, these ancestral responses result in maladaptive pathologies including pulmonary hypertension and an overproduction of red blood cells (polycythemia) 22 . As a result, many species native to high-altitude environments have evolved a reduction or complete loss (canalization) of this ancestral plasticity 23-27 . Acclimation and evolutionary adaptation to environmental change often involves complex coordinated biological responses, such as the co-regulation of genes underpinning physiological function. Reinforcement of adaptive plasticity and reduction/ reversal of maladaptive plasticity at the regulatory and genomic levels likely occur in tandem to drive populations toward new adaptive peaks 28 . However, the relative importance of adaptive and maladaptive plasticity as substrates for selection during the incipient stages of lineage divergence remains largely unexplored.
Although identifying the initial mechanisms that drive adaptive divergence is critical for understanding the evolution of complex systems in novel environments, studying such phenomena becomes more difficult as the time since lineage divergence proceeds, due to the accumulation of mutations, complex demographic histories, and subsequent selective events. Examples of contemporary evolution, such as urban adaptation, provide "natural experiments" that offer unique insights into the mechanisms that underpin the evolution of complex environmentally-responsive traits. Such examples provide means to quantitatively explore evolutionary dynamics during the initial stages of lineage divergence and local adaptation in the wild 29,30 . Further, disentangling the interactions among selection, adaptation, and plasticity has become increasingly important for understanding the evolutionary implications of human-induced environmental alteration 30 . Rapid habitat modification and environmental shifts due to direct and indirect human influence may induce evolution of physiological, behavioral, and/or developmental plasticity [31][32][33][34] . Ancestral plasticity may therefore be a common target of natural selection if it results in greater population resilience to climate change and anthropogenic habitat alteration 32,33,35 .
One such anthropogenic habitat alteration is the increase in ambient temperature in urban areas due to increased impervious surface cover and heat production, collectively termed the "urban heat island effect" 36 . For ectothermic species in particular, elevated Fig. 1 Evolution of plasticity in novel environments. Phenotypic plasticity can influence the evolutionary outcome for a population colonizing a novel environment. Ancestral plasticity may move the initial phenotype of a population (x) in any number of directions with respect to the local optimum of the novel environment (adaptive peak, y). These variable responses are adaptive when they move phenotypic values directly into (A) or close to-but outside -the peak y. If sufficient genetic variation is exposed, natural selection can reinforce adaptive plasticity and move phenotypic means towards the new peak (B; dotted lines). However, natural selection is not expected when plasticity puts individuals directly onto the new peak (A). In the case of maladaptive plasticity (C), inappropriate responses to the novel environment move individuals away from the peak, reducing fitness. Selection should reduce/reverse the reaction norm and restore the phenotypic mean back to the original ancestral value. In all cases, the strength of selection increases with distance from the peak (dotted lines in B, C). Based on Fig. 2 in Ghalambor et al. 4 .
urban temperatures have profound impacts on performance [37][38][39][40][41][42] , reproduction, development [43][44][45] , and survival [46][47][48][49] . A growing body of evidence supports resilience to high temperatures as a key facet of urban adaptation for terrestrial ectotherms [37][38][39]41,42,50 . Furthermore, ancestral plasticity of heat stress responses aid in short-term survival of individuals exposed to periodic acute temperature spikes and may enable initial population persistence in novel urban environments 51 . For instance, behavioral thermoregulation can allow individuals to achieve optimal body temperatures by exploiting thermal heterogeneity, thereby minimizing selection on thermal physiology 52 (Bogert Effect 53 ). However, in extreme environments, temperatures may fall outside the ancestral range of variability or provide inadequate thermal heterogeneity for behavioral thermoregulation. As a result, these novel extreme thermal environments may necessitate physiological adaptation. In such cases, plastic physiological response to short-term thermal variation can result in tradeoffs with other traits critical for survival and reproduction 54,55 , including diminished growth rates, body condition, clutch size, and locomotor performance [56][57][58][59] . Colonization of urban heat islands, where heat stress is likely to be common, chronic, and extreme, may exacerbate the maladaptive consequences of plastic physiological response. Consequently, populations colonizing urban environments may face the multifaceted physiological conundrum of employing heat-induced plasticity that contributes to greater physiological resilience at high temperatures, while mitigating associated maladaptive consequences.
The Puerto Rican crested anole (Anolis cristatellus) is emerging as a powerful model to understand the interaction between natural selection and plasticity during the early stages of adaptation to urban heat islands. Urban habitat modification has been ongoing across Puerto Rico since European colonization five centuries ago and A. cristatellus has occupied urban habitats for at least 35 generations (Supplementary Table 1). Campbell-Staton et al. 60 finds that wild-caught lizards in cities across the island display greater heat tolerance when compared to closely related forest counterparts ( Fig. 2A, B). Urban and forest lineages do not display this physiological divergence when born and raised in common garden environments (Fig. 2C), highlighting plasticity of thermal response to native environments as the primary source of observed physiological divergence in the wild.
These findings set up two distinct hypotheses regarding the role of thermal plasticity in facilitating colonization of urban heat islands: first, ancestral thermal plasticity exhibited by forest populations may be completely sufficient to obtain optimal performance at novel urban temperatures (perfect plasticity hypothesis). Alternatively, plasticity of forest populations may be insufficient to reach the new thermal optimum in urban environments and require the additional action of selection modifying genetic variation underlying thermal plasticity to move the colonizing populations closer to the adaptive peak of the novel environment (evolved plasticity hypothesis). The major distinction between these alternative hypotheses is the action of selection on the genetic underpinnings of thermal plasticity. Campbell-Staton et al. 60 reports parallel patterns of elevated genetic divergence associated with large suites of temperature-responsive genes between urban and forest habitats, and identifies a single nonsynonymous polymorphism that displays genotype-byenvironment interactions associated with elevated thermal tolerance, but only within urban heat islands 60 . The persistent and repeated signatures of genomic selection associated with thermal plasticity observed in the wild are congruent with the evolved plasticity hypothesis as a significant contributor to differences in thermal resilience observed between environments 60 . However, it remains unclear how selection acts on standing variation underlying this regulatory plasticity to facilitate rapid alterations of thermal performance.
In this study, we build upon the findings of Campbell-Staton et al. 60 to examine the relative roles of adaptive and maladaptive gene expression plasticity as targets of selection in urban habitats. First, using the gene expression data from Campbell-Staton et al. 60 , we characterize the regulatory underpinnings of temperature-dependent performance in A. cristatellus to identify a set of candidate genes within the skeletal muscle transcriptome of the hind limb displaying significant correlations with heat tolerance (CT MAX , Fig. 3A, B and Supplementary Fig. S2); this candidate set serves as a proxy for transcriptional drivers of thermal performance, representing the combined contributions of phenotypic plasticity (developmental plasticity and phenotypic flexibility) and evolved divergence (Fig. 4). Next, we explore evolved divergence in plasticity between urban and forest lineages by assessing the relative contributions of adaptive and maladaptive plasticity to candidate gene expression, using lizards born and raised under common laboratory conditions. Finally, we compare genetic divergence and mechanisms of selection associated with adaptive and maladaptive gene expression plasticity in four pairs of urban-forest populations. We do this by measuring the relative magnitude of genetic divergence between habitat types for genes associated with each plasticity category (adaptive or maladaptive, Fig. 7) and identify putative molecular targets of selection in each case (Fig. 8). Using these sequential and complementary analyses, we suggest that rapid adaptation of Anolis lizards to urban heat islands is facilitated primarily by selection for reduced and/or reversed heat-induced plasticity that is maladaptive in urban thermal environments.

Results and discussion
Identifying regulatory underpinnings of temperaturedependent performance. We present novel and extended analyses of the data collected by Campbell-Staton et al. 60 , to further investigate transcriptome-wide regulatory evolution in hind limb skeletal muscle, an essential component of locomotion 61,62 , displaying repeated patterns of temperature-dependent plasticity and signatures of selection in urban heat islands. In brief, Campbell-Staton et al. 60 measured CT MAX for a panel of wildcaught lizards from four urban-forest pairs of A. cristatellus from the municipalities of Aguadilla, Arecibo, Mayagüez, and San Juan, which represent independent urban colonization events across the island of Puerto Rico (n = 114 60 ). Critical thermal maximum (CT MAX ) is a measure of physiological response to acute temperature change 63 that estimates the temperature at which an organism experiences systemic dysfunction-an inability to coordinate locomotor performance 64 . In addition, CT MAX was measured for urban and forest lizards from Mayagüez born and raised in common laboratory conditions (n = 16). Subsequent to CT MAX trials, wild-caught individuals were exposed to one of three 2 h temperature treatments at average field-measured body temperatures from forest (day treatment: 25.4 ± 0.614°C; night treatment: 15.49 ± 1.82°C) or urban (day treatment: 32.09 ± 3.1°C) habitats. Common garden individuals were randomly split between the forest day and urban day treatments only. After the acclimation period, skeletal muscle transcriptomes were collected for RNA sequencing (RNAseq) analyses (see "Methods," and Supplementary Fig. 2).
To categorize the regulatory architecture of the Anolis skeletal muscle transcriptome, we then analyzed the resultant expression data for all individuals (wild-caught n = 114, common garden = 16, Supplementary Fig. 1) using weighted gene correlation network analysis (WGCNA 65 ). Weighted correlation network analysis is a systems biology method for identifying groups of highly co-expressed genes (modules), summarizing module-level expression, identifying intramodular genes of importance (hub genes), and relating expression patterns to trait variation 65 . With this approach, we identified genes across the skeletal muscle transcriptome (total gene n = 11,654) whose expression was significantly associated with CT MAX in A. cristatellus using the gene significance function (GS 65 ). To account for modular structure within the transcriptome, we separated genes into seven gene expression modules, which define the entire regulatory architecture identified by WGCNA 60 (Fig. 3A, B). Across these modules, we identified a total of 632 genes with expression values significantly associated with CT MAX , after correcting for multiple testing (GS.q < 0.05) (Fig. 3B). This set of CT MAX -associated genes are considered putative candidate genes for thermal performance or all subsequent analyses.
Functional enrichment analysis of these candidate genes revealed enrichment of biological processes associated with cellular heat shock response. Three interrelated biological processes (GO:0044089-positive regulation of cellular component biogenesis, p = 0.036; GO:0045898-regulation of RNA polymerase II transcriptional preinitiation complex assembly, p = 0.0004; GO:0044265-cellular macromolecule catabolic process, p = 0.004) share a core set of genes associated with the RNAII transcription preinitiation complex (Supplementary Table 2), which plays a critical role in regulating cellular response to heat shock by initiating transcription of genes that produce heat shock proteins 66 Lines connecting treatments highlight the magnitude of mean CT MAX difference between common garden and native environments. No significant differences in CT MAX are observed between urban and forest lineages under common garden conditions, supporting plastic response to divergent native thermal environments as the primary source of observed in situ differences. Ancestral variation in this plastic response may be the target of selection in the novel thermal environments of urban heat islands.
suite of genes involved in DNA replication and repair, which are critical for maintenance of DNA integrity during and after cellular heat shock 68 . These genes include several that encode kinesin, condensin, kinetochore, RNA helicase, and type II DNA topoisomerase proteins (Supplementary Table 2). Together, these functional associations highlight potential cellular mechanisms underpinning individual variation in heat tolerance and potential targets of selection within urban heat islands. We combined all individuals from common garden and wildcaught groups for initial candidate gene identification, in order to recover the most robust set of thermal tolerance-associated candidate genes in the species. Consequently, our candidate gene set captures expression-phenotype correlations due to the combined effects of developmental thermal environment (developmental plasticity), temperature-dependent expression in adult individuals (phenotypic flexibility), and interactions therein. Capturing this variation in potential regulatory contributions was deemed critical for subsequent analyses, as selection on any of these components may play an important role in adaptive temperature-dependent gene expression divergence between lineages.
We next validated the ability of these candidate genes to recapitulate expression-phenotype correlations independently in wild and common garden environments. We re-ran WGCNA twice more as outlined above, separating individuals from the common garden and wild-caught groups. We then recalculated GS for each candidate gene and compared them across both data sets (Fig. 3C). There was a significant correlation between expression-phenotype values of the wild-caught and common garden data (linear model, adjusted R 2 : 0.601; p « 0.001). Relative rank of GS was also highly conserved across groups (Spearman's rank correlation ρ: 0.631, p « 0.001). Among genes that displayed negative correlations with CT MAX (negative regulators), we observed no significant differences in GS between wild and common garden conditions (paired t-test, t = 0.81, df = 295, p = 0.42). Genes positively correlated with CT MAX (positive regulators) displayed significantly higher GS in common garden than in the wild (paired t-test, t = 5.58, df = 335, p « 0.001), suggesting a greater ability to predict expression-phenotype relationships in common garden. Higher variance was observed in the common garden data set, resulting in lower explanatory  60 . Bars represent the number of genes within each module. Inset is a representation of connectivity among co-expression modules. Circle sizes represent the relative size of each module (circle area is proportional to log-transformed number of genes within each module) and line thickness represents the relative strength of connectivity between modules (absolute value of pairwise Pearson's correlation among eigengene values). B Violin plot of gene significance scores (GS) for CT MAX (strength of association between gene expression and CT MAX ). Dots indicate individual genes assigned to each module. Colored dots within each module indicate genes with significant phenotypic correlations after multiple testing correction (GS.q < 0.05). This subset of genes made up the focal data set of the current study. C To validate the candidate gene set, two additional network analyses were conducted for wildcaught (n = 114) and common garden lizards (n = 16), separately. Gene significance (GS) scores were calculated for candidate genes independently for each group. Genes identified as negatively correlated with CT MAX (negative regulators) in the full data set are indicated in blue; genes positively correlated with CT MAX in the full data set (positive regulators) are indicated in red. Black dots in the center of each violin plot are GS scores for each individual gene. Gray lines present the directionality of GS score change for each gene between common garden and wild-caught animals. There was a significant expression-phenotype correlation between wild-caught and common garden data (linear model, R 2 : 0.601; p « 0.001). Relative rank of gene significance was highly conserved across sets (Spearman's rank correlation ρ: 0.75, p « 0.001). Negative regulators displayed no bias in gene significance between groups (paired t-test, p = 0.26). Positive regulators displayed higher gene significance values in the common garden data set than the wild-caught data set (paired t-test, p « 0.001). The gene-gene correlations between the wild-caught and common garden data sets was small but statistically significant in each case (positive regulators: R 2 : 0.04, p < 0.0001; negative regulators R 2 : 0.03, p = 0.003). power within each regulatory category. However, the correlation between the wild-caught and common garden data sets remained significant in each case (positive regulators: R 2 : 0.04, p < 0.0002; negative regulators R 2 : 0.03, p = 0.0003). In summation, our candidate gene set displays similar expression-phenotype correlations and relative rank correlations between common garden and wild-caught individuals, supporting the gene set as suitable for analysis in both native and common garden environments.
To directly assess evolutionary contributions to expression variation in these candidate genes, we searched for genomic signatures of selection in each population using all individuals (wild-caught = 114, common garden = 16). If the plasticity producing observed differences in thermal tolerance in the wild habitats, suggesting ongoing selection on thermal physiology island-wide. In addition, increased genetic divergence between urban and forest habitats (Column 3) in CT MAX -associated SNPs (red) with respect to background (gray) expectations (Arecibo, Mayagüez, and San Juan: Welch t-test; p < 0.05 (asterisks)) supports selection on heat tolerance specific to urban heat islands. LD plots are represented in 1 kb distance bins. Dots and error bars represent means ± 1 SE. All 130 animals for which we had sequence data were used for this analysis.
are the target of selection in urban habitats, we would expect to see significant signatures of selection at the sequence level specific to our candidate loci. We identified 2,161,083 single-nucleotide polymorphisms (SNPs) transcriptome-wide using RNAseq reads from lineages originating in all four municipalities studied in Campbell-Staton et al. 60 (Aguadilla: forest n = 11, urban n = 16; Arecibo: forest n = 18, urban n = 16; Mayaguëz: forest n = 19, urban n = 18; San Juan: forest n = 11, urban n = 21). Using all polymorphic loci within candidate genes as the focal candidate SNP set, we tested for significant deviations from neutrality by direct comparison with the null expectation estimated by background variation 69 across all SNPs within genes that had no significant association with CT MAX (Fig. 4). We calculated pairwise linkage disequilibrium (LD; R 2 within non-overlapping 10 kb windows) and tested for deviations from Hardy-Weinberg (HW) equilibrium across all polymorphic sites of the transcriptome. We then compared patterns of LD (binned in 1 kb increments) between candidate and background SNP sets. As LD and HW equilibrium can be influenced by both selection and demography, we restricted our analyses to comparisons within each site such that the background and candidate SNP sets shared a common demographic history. If our candidate genes have been a common target of selection, we would expect to observe a significant increase in LD at polymorphic sites proximate to these genes when compared to the transcriptomewide background. Consistent with this expectation, we found that SNPs within CT MAX -associated genes displayed significantly elevated LD (linear mixed-effects model: p < 0.05, Fig. 4) across all sampled localities. Significant deviations from HW equilibrium were also more prevalent among SNPs within CT MAXassociated genes compared to the background SNP set (equality of proportions test-Arecibo: χ 2 = 77.721, df = 1, p-value << 0.001; Aguadilla: χ 2 = 28.478, df = 1, p-value << 0.001; Mayagüez: χ 2 = 37.71, df = 1, p-value << 0.001; San Juan: χ 2 = 33.383, df = 1, p-value = << 0.001; Supplementary Table 3).
These genomic patterns suggest that the regulatory underpinnings of heat tolerance are broadly under selection across the range of A. cristatellus, supporting previous work demonstrating the integral role of thermal biology in habitat specialization among Anolis occupying divergent macro-and microclimatic niches [70][71][72] . If, in addition to this broader selection on thermal physiology, populations are under significant contemporary selection within urban heat islands, we would also expect to observe increased genetic divergence (F ST ) across urban-forest boundaries when comparing SNPs within candidate loci to null expectations from the transcriptome-wide background. We found this to be the case in three of our four urban-forest comparisons (Arecibo, Mayagüez, and San Juan: Welch two-sample t-test, p << 0.001, Fig. 4), suggesting that CT MAX -associated genes are repeated targets of selection between urban and forest habitats across the island.
Evolved divergence in gene expression plasticity between urban and forest lineages tested with a common garden experiment. To categorize heat-induced plasticity of a given gene as adaptive or maladaptive, we made the assumption that greater heat tolerance is adaptive at higher temperatures 51 , whereas lower heat tolerance is maladaptive, resulting in reduced fitness 73,74 . Using the GS scores from the WGCNA analyses outlined above, we categorized the direction of expression-phenotype correlation as positive or negative for each CT MAX -associated candidate gene. We then characterized differences in expression observed for each common garden group exposed for 2 h, to either average forest (25°C; forest n = 4, urban n = 5) or average urban (32°C, forest n = 3, urban n = 4) body temperatures, as congruent (adaptive) or incongruent (maladaptive) with increased heat tolerance. Cases in which the direction of common garden plasticity matched the expectation for greater CT MAX were deemed putatively adaptive. Those that opposed this expectation were deemed putatively maladaptive. For an outline of the categorization workflow, see Fig. 5. Forest lineages of A. cristatellus predate anthropogenic habitat modification. Therefore, we used the common garden forest group (n = 7; 25°C n = 3, 32°C = 4) as a proxy for the ancestral state of gene expression for comparison against the derived (urban) gene expression state (n = 9, 25°C n = 5, 32°C = 4).
We found that urban and forest groups from Mayagüez diverged significantly in heat-induced gene expression plasticity under common garden conditions, supporting evolved divergence between lineages. Forest lineages displayed a disproportionately maladaptive response to heat stress (82.3% of genes, n = 520, binomial test: p << 0.001, Fig. 6A, B), with the majority of positive regulators (87.2%) being downregulated (n = 293, p << 0.001) and the majority of negative regulators (76.7%) being upregulated (n = 227, p << 0.001) in response to increased temperature. These results indicate ancestral heat-induced plasticity of gene expression is predominantly associated with lower CT MAX and may therefore be maladaptive for physiological performance in urban heat islands. In contrast, urban lizards from Mayagüez displayed significantly less maladaptive plasticity (forest = 82.3%, urban = 54.3%, exact binomial p < 0.001) and significantly more adaptive plasticity (forest = 17.7%, urban = 45.7%, exact binomial p < 0.001) across candidate genes when exposed to the urban temperature treatment (Fig. 6B).
The observed divergence in proportions of adaptive and maladaptive responses in common garden lineages suggests that selection in urban habitats has acted primarily to reduce and/or reverse maladaptive gene expression plasticity. Under such a scenario, we would expect to observe an overall reduction of heatinduced gene expression plasticity in derived urban lineages 75 . As a result, the direction of evolutionary divergence should be negatively correlated with the direction of ancestral plasticity 75 . To test this hypothesis, we estimated ancestral plasticity of a given gene as the log-fold change of expression in common garden forest lizards across the two temperature treatments (25°C n = 3, 32°C = 4) and evolved divergence as the log-fold change in expression between the forest and urban common garden groups at 32°C (urban n = 4, forest n = 4). We then quantified the direction and strength of the correlation between ancestral plasticity and evolved divergence across the entire CT MAXassociated candidate gene set for the Mayagüez populations.
Consistent with expectations of predominant selection against maladaptive plasticity, ancestral plasticity and evolved divergence displayed a strong negative correlation (Spearman's ρ = −0.61; 82.3% of genes, n = 520, equality of proportions test, p < 0.001, Fig. 6A). To account for potential statistical artifact due to regression towards the mean 76 , we conducted a randomization test to estimate a null distribution of correlation coefficients and found that the strength of our observed correlation was significantly more negative than expected as a result of artifact alone (empirical p < 0.05, Fig. 6A inset). If selection has acted primarily to reduce the magnitude of maladaptive plasticity, we would also expect urban lineages to display a significant reduction in the magnitude of heat-induced gene expression response when compared to forest counterparts. Therefore, we directly compared the magnitude of heat-induced plasticity between common garden forest and urban groups at 25°C vs. 32°C (forest: 25°C n = 3, 32°C = 4; urban: 25°C n = 5, 32°C = 4) and found that urban lizards displayed a blunted heat-induced gene expression response (lower log-fold change in expression) compared to forest lizards (Welch two-sample t-test, p < 0.001, Fig. 6C). Taken Theoretical workflow for categorizing the ancestral plastic response of a given gene to increased urban temperatures as adaptive or maladaptive. First, the direction of gene expression correlation (measured from a panel of 130 lizards) was used to identify a given gene as a positive regulator (positively correlated with heat tolerance, CT MAX ) or a negative regulator (negatively correlated with CT MAX ). Next, the gene was categorized as upregulated (increasing in expression) or downregulated (decreasing in expression) in forest lineages born and raised in common conditions (a proxy for the ancestral condition) in response to increased temperature (25°C-32°C). Theoretical dots and bars represent mean ± 1 SE. Positive regulators displaying increased expression at high temperature and negative regulators displaying decreased expression at high temperature in common garden lizards were categorized as displaying adaptive plasticity. Decreased expression of positive regulators at high temperature and increased expression of negative regulators at high temperature in common garden individuals were deemed maladaptive plasticity. Evolved gene expression divergence between urban and forest lineages measured through a common garden experiment. Lizards from urban and forest lineages from Mayagüez were born and raised in a common garden laboratory setting to study patterns of evolved divergence in gene expression between habitats types. A The relationship between ancestral (forest) plasticity of gene expression (25°C vs. 32°C) and evolved divergence in gene expression under heat challenge (forest vs urban at 32°C). Each dot represents a single gene from the focal set of CT MAX -associated genes (n = 632). Dots displayed in black represent genes with a positive correlation between ancestral plasticity and evolved divergence. Those displayed in white represent genes with a negative association between ancestral plasticity and evolved divergence (Spearman's ρ = −0.61). Among genes associated with heat tolerance, evolved divergence reverses the direction of ancestral plasticity. Inset: the observed correlation coefficient (red line) is more negative than the lower 95th percentile of a randomized null distribution (black line). B Bar graph representing the number of genes displaying adaptive vs. maladaptive plasticity in common garden lineages from forest (green) and urban (gray) habitats. Figure 5 outlines the methodology of plasticity categorization. Among genes associated with heat tolerance, urban lineages possess a significantly larger proportion displaying adaptive plasticity and a significantly lower proportion that displaying a maladaptive plasticity in response to increased temperature. Asterisks represent degree of significance (equality of proportions test within habitat types, exact binomial test between habitat types: *p < 0.05, ***p < 0.001). C Violin plot displaying heat-induced regulatory plasticity (logfold change in expression from 25°C to 32°C) in forest (green) and urban (gray) lineages. Individual genes are represented by gray dots. Urban lineages display an overall lower magnitude of gene expression plasticity (Welch two-sample t-test: p < 0.001), congruent with evolutionary attenuation of maladaptive plasticity. Sixteen animals from common garden conditions were used for these analyses (25°C; forest n = 4, urban n = 5) or average urban (32°C, forest n = 3, urban n = 4).
together, the differences observed between forest and urban groups from Mayagüez in common garden support evolutionary divergence, driven primarily by the reduction and/or reversal of maladaptive plasticity present in ancestral forest lineages.
Comparing genetic divergence and mechanisms of selection for adaptive and maladaptive gene expression plasticity. Results of the common garden experiment suggest that temperaturedependent selection against maladaptive plasticity has played a predominant role in evolutionary gene expression divergence between urban and forest populations of A. cristatellus. Under this hypothesis, we would expect genetic divergence among genes displaying maladaptive plasticity to exceed genetic divergence associated adaptive plasticity and neutral rates of divergence observed broadly across the rest of the transcriptome. Furthermore, differences in the rate of genetic divergence observed across plasticity categories may suggest different targets of selection acting to drive local adaptation in urban environments. Therefore, we directly compared patterns of genetic divergence between urban and forest habitats to test the prediction that loci underpinning maladaptive plasticity diverge more rapidly between novel (urban) and ancestral (forest) lineages compared to those underpinning adaptive plasticity. To test this prediction, we divided CT MAX -associated candidate SNPs into those occurring within the putatively adaptive and maladaptive plasticity categories prescribed above. Significantly elevated genetic divergence (F ST ) and enrichment of F ST outliers observed within the candidate SNP set beyond the background expectation would provide evidence for selection acting on this group of interacting phenotype-associated genes (polygenic selection) 77,78 .
We first quantified genetic divergence associated with plasticity within each pair of urban-forest populations by testing for statistical enrichment of F ST outliers within the adaptive and maladaptive plasticity classes. These comparisons of genetic divergence were conducted independently for each of the four paired urban-forest groups using sequence data from all individuals (wild-caught n = 114, common garden n = 16) representing each municipality. We used a p < 0.05 threshold based on the empirical distribution of genetic divergence for the transcriptome-wide background. We found that genes displaying adaptive plasticity were significantly enriched for F ST outliers in three of the four urban-forest comparisons (Arecibo, Mayagüez, and San Juan), whereas genes displaying maladaptive plasticity were significantly enriched for outliers in all four comparisons (Fig. 7). These results support the action of selection on both adaptive and maladaptive gene expression plasticity.
We then assessed the relative strength of selection on adaptive and maladaptive plasticity by comparing mean F ST between each urban-forest pair (wild-caught n = 114, common garden n = 16). We found that genes associated with maladaptive plasticity had greater genetic divergence across urban-forest boundaries than those displaying adaptive plasticity in all four municipalities (Welch two-sample t-test, p < 0.01) (Fig. 7). Although both adaptive and maladaptive plasticity appear to be common targets of selection in urban heat islands, these data suggest that selection operates more intensely on maladaptive plasticity during the initial stages of genetic divergence across urban-forest boundaries.
Finally, we examined potential genetic mechanisms underpinning adaptive and maladaptive plasticity by evaluating the proportion of noncoding and coding polymorphisms associated with genetic divergence between forest and urban lineages (wildcaught n = 114, common garden n = 16). If selection is acting primarily on functional variation associated with nonsynonymous polymorphisms within coding regions, we would expect to observe a significant excess of nonsynonymous variants within our CT MAX -associated candidate SNP set when compared to the transcriptomic background. Alternatively, if selection is being driven primarily by genetic variation in cis-regulatory elements, we may expect a significant excess of noncoding polymorphisms.
We found that genetic variation within genes putatively displaying adaptive plasticity contained a significant excess of nonsynonymous polymorphisms when compared to the background transcriptome (equality of proportions test: χ 2 = 7.50; df = 1; p = 0.006, Fig. 8A). In contrast, genes displaying putatively maladaptive plasticity displayed no such enrichment (χ 2 = 0.83; df = 1; p = 0.36, Fig. 8A). These results suggest that function-altering coding variation is the primary target of selection for adaptive plasticity in urban heat islands. The lack of enrichment of nonsynonymous variation within genes displaying maladaptive plasticity suggests, alternatively, that cisregulatory variation may be the primary driver of selection against maladaptive plasticity.
To directly test this latter hypothesis, we identified SNPs within noncoding regions of the skeletal muscle transcriptome (putative  Fig. 7 Repeated increases in urban-forest genetic divergence associated with maladaptive plasticity. Comparisons of genetic divergence (F ST ) among heat tolerance-associated genes in four municipalities, representing independent urban colonization events across Puerto Rico 58 . Bar graphs display the proportion single-nucleotide polymorphisms (SNPs) within the adaptive (black) and maladaptive (white) gene sets that occur beyond the p = 0.05 empirical distribution of genetic divergence (F ST ), based on the transcriptome-wide background (gray). Asterisks indicate a significant enrichment of outliers within a gene class (*p < 0.05, **p < 0.01, ***p < 0.001, two-sample equality of proportions test with continuity correction). No outlier gene SNPs were found in the adaptive plasticity gene set within Aguadilla. Dot plots display mean ± 1 SE of genetic divergence (F ST ) for SNPs within the adaptive (black) and maladaptive (white) plasticity gene classes. Asterisks indicate significant differences in mean divergence (*p < 0.05, **p < 0.01, ***p < 0.001, Welch two-sample t-test). All 130 animals for which we had sequence data were used in this analysis.
cis-regulatory elements). We found that noncoding variants were overrepresented among genes displaying putatively maladaptive plasticity (χ 2 = 96.646, df = 1, p << 0.001, Fig. 8B) when compared to the transcriptome-wide background. However, noncoding SNPs associated with adaptive plasticity display no significant enrichment over transcriptome-wide proportions (χ 2 = 1.926, df = 1, p = 0.165, Fig. 8B). Taken together, these data suggest that adaptive and maladaptive plasticity may not only diverge at different rates across urban-forest boundaries, but that selection on adaptive and maladaptive plasticity may have occurred, at least in part, through different mechanisms in urban heat islands -with selection on beneficial coding mutations driving evolution favoring adaptive plasticity and selection against deleterious cisregulatory variants promoting evolution to minimize and/or reverse maladaptive plasticity.
Considerations and limitations of the methodological approach. Previous work has highlighted the complications of interpreting the relationship between ancestral plasticity and evolved divergence in transcriptomic data, as negative correlations emerge due to statistical artifact (i.e., regression towards the mean 79 ). We have taken several measures in this study to account for the potential influence of this phenomenon. First, regression towards the mean is most likely to result when the most extreme cases in one comparison (i.e., differentially expressed genes between two groups) are reassessed in a subsequent comparison 79 . The identification of candidate genes in this study is not based on differential expression with regard to either ancestral plasticity nor evolved divergence, but on the strength of regulatory association with an environmentally relevant phenotype (CT MAX ) with multiple lines of support as a target of selection in urban habitats 60 . In addition, randomization of the expression data reveal that the magnitude of correlation between ancestral plasticity and evolved divergence observed here is significantly greater than expected by chance 75,80 (Fig. 6a inset). Most importantly, such statistical artifact arising from neutral variation in regulatory data is not expected to produce patterns of increased genetic divergence between lineages. The parallel signatures of increased genetic divergence among independently derived urban lineages further supports the role of natural selection in the evolution of the observed gene expression plasticity.
The major goal of this study is to understand the evolutionary mechanisms driving evolved divergence in gene expression associated with the observed difference in temperaturedependent performance between urban and forest anoles. Towards this goal, we test several independent but interrelated hypotheses to understand evolution of gene expression associated with differences in thermal tolerance between forest and urban populations in the wild. Our approach is dependent on isolating and focusing on gene expression patterns associated with thermal performance to assess patterns of evolved gene expression divergence that are relevant to the ecological context. These gene expression and sequence level analyses have provided multifaceted support for a significant evolutionary contribution to the observed difference in thermal performance observed in the wild.
However, our candidate gene approach likely only represents a partial picture of adaptive modification of a more complex, temperature-dependent phenotype (CT MAX ). In order to gain a more complete understanding of the higher-order implications of lineage-specific divergence in expression among candidate genes and their underlying mechanisms, more extensive common garden study of whole organism thermal performance is needed 81 . Utilizing parallel changes across multiple urban-forest pairs, controlled breeding and experimental designs that isolate the impacts of developmental environment (developmental plasticity) and acute temperature change (phenotypic flexibility) on patterns of whole organism performance and its regulatory underpinnings would further contextualize the results of this study and provide greater insights regarding the importance of phenotypic plasticity as a means of rapid adaptation during the incipient stages of divergence. Nevertheless, the approaches presented herein demonstrate the utility of integrating data from field-based studies, common garden experiments, and multiomics data, to reveal detailed mechanisms driving adaptive divergence between lineages. In this case, the results of our study provide fundamental insights into the relative roles of adaptive and maladaptive plasticity in driving the rapid evolution of complex regulatory systems in novel environments over anthropogenic timeframes.
Understanding the interactions between phenotypic plasticity and evolution remains an area of active debate within evolutionary biology 42,82 . This is largely due to a general lack of data regarding ancestral variation in environmentally induced plastic responses and limited information regarding their Evidence for different targets of selection on adaptive and maladaptive plasticity. A Single-nucleotide polymorphisms (SNPs)-associated adaptive plasticity display a significant enrichment of nonsynonymous polymorphisms, whereas those associated with maladaptive plasticity do not. B Singlenucleotide polymorphisms-associated maladaptive plasticity display a significant enrichment of noncoding polymorphisms), whereas those associated with adaptive plasticity do not. Dots and bars represent the proportion of SNPs among genes associated with adaptive (black) and maladaptive (white) plasticity, and 95% confidence intervals, respectively. Horizontal dotted lines represent expected proportions estimated from transcriptome-wide background (SNPs among genes unassociated with heat tolerance). All 130 animals for which we had sequence data were used in this analysis.
adaptive value or subsequent evolutionary impacts in the face of novel environmental challenges. Among studies that have focused on the evolution of plasticity, most have focused explicitly on adaptive plasticity. Although evolutionary change that reduces non-adaptive plasticity has been suggested as a common form of adaptation in the wild 83 , it has received relatively little attention 84 . Our findings support the role of phenotypic plasticity in adaptation to novel environments and highlight reduction/ reversal of maladaptive plasticity as the predominant target of natural selection in urban heat islands. These findings are generally consistent with an emerging pattern: although adaptive plasticity facilitates survival and population persistence at early stages of colonization, genetic changes that minimize and/or reverse, rather than reinforce, phenotypic plasticity appear to be stronger and more frequent targets of selection 28,[85][86][87] .
The common garden experimental design conducted in this study uses non-sibling individuals spread across temperature treatments. As a result, the plasticity reported herein represents plasticity of gene expression at the population level and does not enable us to quantify individual reaction norms or account for background genetic variation among individuals, which we acknowledge as a limitation of the study. However, populationlevel phenotypic differences observed in the wild are the focus of our investigation and despite high degree of gene flow between urban and forest habitats within each municipality, the populations are genetically differentiated. As such, the common garden experimental design does allow for interrogation of environmentally induced gene expression changes underpinning evolved divergence between habitat types.
The predominance of temperature-induced maladaptive plasticity observed in ancestral forest lineages suggests that thermal stress is a major selective pressure during the initial stages of urban colonization. The observed parallel patterns of genetic divergence are consistent with previous work on this system 60 and suggest that behavioral thermoregulation is insufficient to shield these tropical urban populations from selection on temperature-dependent physiological performance 88,89 . The enrichment of gene ontology categories associated with heat shock protein transcription and DNA replication suggest that mediation of proteotoxicity due to the aggregation of misfolded proteins plays an important role in determining individual variation in ectothermic heat tolerance in urban heat islands. Together, this group of genes displays both broad signatures of polygenic selection across the island (consistent with thermal heterogeneity, and micro-and macroclimatic habitat specialization) and repeated divergence between urban and forest habitats. The maladaptive plasticity exhibited by forest lineages appears to be a primary target of natural selection in this latter context, resulting in a reduction and reversal 28 of maladaptive heatinduced responses and a significant shift towards adaptive plasticity in urban lineages. Furthermore, selection on adaptive and maladaptive plasticity appear to occur via different genetic mechanisms-nonsynonymous and cis-regulatory variation, respectively.
The candidate gene approach that forms the foundation of this study uses the correlations between gene expression and CT MAX to infer a group of regulatory drivers of whole organism performance. Although multiple lines of evidence presented here support these genes as direct targets of selection associated with thermal performance, further common garden study will be needed to directly quantify the evolution of thermal plasticity at the whole organism level. Indeed, the results of our study emphasize the necessity of simultaneous and integrative investigations of adaptive and non-adaptive plasticity across multiple levels of biological hierarchy to fully understand evolution and adaptation of complex traits in novel habitats. This integration of methods may be critical for predicting organismal response to anticipated anthropogenic habitat alteration and changes in future climate.

Methods
Raw data collected for this study are the same as those in Campbell-Staton et al. 60 .
All animal experiments were approved by the University of Massachusetts Boston under Institutional Animal Care and Use Committee (IACUC) protocol number 2012001. We summarize the data collection details here followed by the analytical methods unique to this study.
Source populations. Gene expression analyses for this study were performed on two non-overlapping and unrelated groups of lizards: 114 wild individuals ("wildcaught" group) collected in 2016 and 16 captive, reared individuals ("common garden" group). Details of the experimental design and sample sizes are provided in Supplementary Data File 1.
The 114 wild-caught adult male A. cristatellus used in the study were captured using standard methods (hand capture and floss lasso) as encountered from paired urban and forest sites in four municipalities in Puerto Rico: Aguadilla, Arecibo, Mayagüez, and San Juan ( Fig. 2A). Paired sites were within 10 km of each other. We placed lizards in cloth bags for transport to laboratory conditions where they were individually housed for the duration of the experiments. Lizards were allowed to acclimate to room temperature and humidity for 24 h prior to the initiation of thermal tolerance and subsequent acclimation experiments and transcriptome sequencing.
We established a captive common garden colony of A. cristatellus as part of a separate study 90 . In August 2013, we captured male-female pairs of adult A. cristatellus from the same urban site and a forest site sampled in 2016 in Mayagüez, Puerto Rico. We transported lizards to the animal care facility at the University of Massachusetts, Boston, where we paired individuals within their respective populations for breeding. We collected eggs from August 2013 to April 2014, placing them in hydrated vermiculate in an incubator (28°C) until hatching. Hatchling lizards were immediately moved to individual cages in the animal care facility. All lizards were kept at a 13 h daylight schedule with a constant day/night temperature of 26.5°C, daily misting, and vitamin dusted crickets provided every 3 days. All cages were rearranged regularly in the animal facility to minimize any effects of localized thermal variation. In 2016, 16 individuals were chosen from unique families for inclusion in the present study. Common garden lizards selected for the present study were F1 generation individuals 3-3.5 years of age, reared under identical conditions, and no full or half siblings were included in any of the experiments presented herein.
Thermal tolerance and acclimation. Full methods for all thermal tolerance trials are presented in ref. 60 , but we present relevant methods for this study here. We estimated CT MAX by heating lizards at a rate of 1°C min −1 , starting at room temperature, while monitoring body temperature with a thermocouple inserted into the cloaca. After initiation of thermal stress responses (gaping, panting, and lethargy), we tested the righting response by placing the lizard onto its back. The temperature at which the lizard was no longer able to right itself within 30 s was recorded as its CT MAX , at which point the lizard was immediately cooled to room temperature in a water bath. CT MAX estimates generally ranged between 35°C and 43°C. Thermal trials were conducted on all 130 lizards of this study (114 wildcaught, 16 common garden, see Supplementary Fig. 4 and Supplementary Data File 1).
Lizards were allowed a 24 h recovery period following CT MAX trials. We then subjected wild-caught lizards to one of three thermal acclimation conditions for 2 h, approximating nocturnal forest (15°C ± 1.82°C), diurnal forest (25.4°C ± 0.614°C), or diurnal urban (32.09°C ± 3.1°C) conditions. Common garden animals were only subjected to the diurnal forest or diurnal urban acclimations. Following the acclimation period, lizards were anesthetized with 5% aerial isoflurane and killed via cervical dislocation. We immediately collected skeletal muscle from hind limbs and stored tissues in RNAlater at −80°C. Thermal acclimation treatments were conducted for 130 lizards (114 wild-caught, 16 common garden, see Supplementary Fig. 4 and Supplementary Data File 1).
Transcriptome sequencing and filtering. Total RNA was extracted from 130 tissue samples (114 wild-caught, 16 common garden) using Qiagen RNeasy Fibrous Tissue Kits. Messenger RNA libraries were prepared and sequenced for 100 base pair (bp) single-end reads on the Illumina Hiseq 4000, resulting in an average of 22 M ± 2.6 reads. Trimmomatic 91 was used to assess read quality in 4 bp sliding windows and sequences were trimmed when Phred33 quality scores fell below an average of 15 within a window. Sequences were then mapped to the Anolis carolinensis genome 92 using Tophat2 93 .
Statistical analyses. WGCNA 65 was used to obtain a correlation-based de novo regulatory architecture for the skeletal muscle transcriptome using all samples for which transcriptome data were available (n = 130; 114 wild-caught, 16 common garden, Supplementary Fig. 1). Regulatory modules were identified as branches of the resulting cluster tree via the dynamic tree-cutting method and highly correlated modules (R 2 = 0.75) were merged. The GS function in WGCNA was used to identify genes associated with individual variation in CT MAX , correcting for multiple testing within each identified regulatory module. An identical procedure was used to validate the candidate gene set, using only individuals from common garden ( Supplementary Fig. 2) and wild-caught ( Supplementary Fig. 3) animals, respectively.
Gene Ontology enrichment was performed us the gProfileR package in R 94 , using an unordered query analysis with a false discovery rate correction. The list of all genes represented within the skeletal muscle transcriptome was used as a customized background for analysis. All calculations of log-fold change in expression were obtained using the edgeR software package in R 95 . Filtered RNAseq reads were mapped to the A. carolinensis reference genome with STAR v.2.5.b2 96 using a single-pass approach with Ensembl v.90 annotations, retaining reads that were uniquely mapped with fewer than ten mismatches.
For sequence level analyses in this study, we used all individuals from Campbell-Staton et al. 60 for which RNAseq data were available (n = 130; 114 wildcaught, 16 common garden). Variants were called with the GATK HaplotypeCaller and GenotypeGVCF v.3.5 97 using the recommended RNAseq parameters (-stand_call_conf 20 -stand_emit_conf 20 -dontUseSoftClippedBases) and filtered for binary SNPs with quality ≥ 20, minor allele frequency ≥ 5%, and <20% missing data. Genetic divergence estimates were calculated using the analysis of molecular variance (AMOVA) method with a 5 Mb sigma parameter for smoothing within the Stacks software package (version 1.47) 98 . LD analyses were analyzed for 1 kb bins for each collection site using linear mixed-effects models. Pairwise LD was used as the response variable and gene set (a priori vs. background) was assigned as the predictor variable. Both SNP positions in each comparison were modeled separately as random effects. Differences in genetic divergence (F ST ) were analyzed by Welch two-sample t-test. HW deviations were calculated with vcftools 99 and differences in proportions were estimated by two-sample test for equality of proportions with continuity correction. Enrichment of nonsynonymous and noncoding SNPs among groups was assessed by two-sample equality of proportions test with continuity correct, using proportions obtained from the background data set as the null hypothesis. All analyses were conducted with the R statistical software package 100 .
Gene expression data were analyzed using the R package edgeR 95 . Raw read counts were first normalized by total library size and log-transformed using the edgeR functions calcNormFactors and cpm, respectively. The correlation between evolved divergence (average log-fold change in expression between forest vs. urban lizards at 32°C) and ancestral plasticity (log-fold change between forest lizards acclimated to 25°C vs. 32°C) were examined with Spearman's rank correlation coefficient (ρ) in cor.test in R. To account for potential statistical artifact in the regression, a null distribution of correlation coefficients was produced using a randomization and resampling prodecure: normalized gene expression data among the candidate gene set was randomized 1000 times using the randomizeMatrix function in the R package picante 101 , after which we recalculated average log-fold change in expression for ancestral plasticity, evolved divergence, and Spearman's ρ. The upper and lower 95% limits of the null distribution were −0.48 and −0.59, respectively (Fig. 3A inset).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The sequence and metadata data used in this study are available in the NCBI database under accession code PRJNA592594. The data can be directly accessed at. Source data are provided with this paper.

Code availability
All statistical analyses have been conducted according to methods outlined herein, implemented in the associated referenced software packages. Code for reproduction of the figures herein are provided in the supplemental materials.