Genomic and phenotypic signatures of climate adaptation in an Anolis lizard

Abstract Integrated knowledge on phenotype, physiology, and genomic adaptations is required to understand the effects of climate on evolution. The functional genomic basis of organismal adaptation to changes in the abiotic environment, its phenotypic consequences, and its possible convergence across vertebrates are still understudied. In this study, we use a comparative approach to verify predicted gene functions for vertebrate thermal adaptation with observed functions underlying repeated genomic adaptations in response to elevation in the lizard Anolis cybotes. We establish a direct link between recurrently evolved phenotypes and functional genomics of altitude‐related climate adaptation in three highland and lowland populations in the Dominican Republic. We show that across vertebrates, genes contained in this interactome are expressed within the brain, the endocrine system, and during development. These results are relevant to elucidate the effect of global climate change across vertebrates and might aid in furthering insight into gene–environment relationships under disturbances to homeostasis.

be difficult to discern from the genes that adapt to other factors of selection along these gradients (e.g., available resources). This problem can be approached through the integration of independent evidence from multiple taxa, or by examining natural replicates. In two such studies, text-and genome mining has been used to reveal that candidate genes for climate adaptation initially identified from diverse lineages of vertebrates such as human and fish are in close functional relationship with each other (Porcelli, Butlin, Gaston, Joly, & Snook, 2015;Wollenberg Valero et al., 2014). Wollenberg Valero et al. (2014) found a network of genes identified from over 150 publications on vertebrate temperature adaptation and stress response to be functionally more closely related than expected by chance. Porcelli et al. (2015) subsequently found evidence for convergent functions of thermal adaptation genes among all Metazoans. Predicted functions align to empirically derived biochemical evidence (e.g., for vertebrate cold and freeze tolerance, Storey & Storey, 2004;Storey, 2006;Storey & Storey, 2017). Overall, these in silico studies indicate that a functional genomic network for thermoregulation and/or thermal adaptation to the environment exists, which might be regulated via similar evolutionarily conserved pathways in different vertebrate lineages, which validates a comparative approach to understanding thermal adaptation. These predictions of a conserved set of functions that adaptively evolve in response to changing abiotic environmental parameters have not been verified in a natural setting yet. Previous studies identified a high level of determinism in the evolution of Caribbean Anolis lizards (Losos, Jackman, Larson, de Queiroz, & Rodriguez-Schettino, 1998;Mahler, Revell, Glor, & Losos, 2010;Wollenberg, Wang, Glor, & Losos, 2013). Communities of Anolis ecology-morphology-behavior specialists (ecomorphs) are filtered by environmental niche parameters associated with elevational gradients (Wollenberg, Veith, & Lötters, 2014) and show evidence for environmental selection of thermal performance (Kolbe, Ehrenberger, Moniz, & Angilletta, 2013;Logan, Cox, & Calsbeek, 2014).
The Hispaniolan species Anolis cybotes ( Figure 1) is a common and widely distributed trunk-ground ecomorph comprised of a cluster of genetically divergent populations including three paraphyletic taxa (A. cybotes, A. armouri, and A. shrevei). Intraspecific genetic divergence in this species is aligned to spatial convergence in perch use, osteomorphological measurements, toe pad scalation, physiology, and behavior (Glor, Kolbe, Powell, Larson, & Losos, 2003;Muñoz et al., 2014;Wollenberg et al., 2013). The repeated nature of these patterns across three elevational gradients inhabited by genetically divergent populations suggests that they evolved under selection, and temperature influencing physiology has been identified as a potential selective agent (Hertz & Huey, 1981;Muñoz et al., 2014;Wollenberg et al., 2013). To date, no comprehensive picture exists about the functional genomic basis of thermal adaptation and its degree of conservation across vertebrates. We hypothesize that A. cybotes populations inhabiting elevational gradients will exhibit signatures of genomic adaptation targeting genes, or gene functions that have previously been predicted to underlie vertebrate thermal adaptation. To test this hypothesis, we use correlative evidence to investigate environment-gene-phenotype interactions across elevational gradients.

| RAD sequencing, genomic population structure, and outlier identification
The A. cybotes species complex (A. cybotes, A. breslini, A. armouri, A. shrevei) is a monophyletic clade with paraphyletic positioning of nominal taxa, in following referred to as A. cybotes pending taxonomic revision. It inhabits a wide range of semi-xeric to mesic habitats from sea level to high altitudes on the island of Hispaniola (Alföldi et al., 2011;Glor et al., 2003;Wollenberg et al., 2013). Specimens of A. cybotes were collected from elevational gradients containing distinct populations identified in a prior publication (Glor et al., 2003), including the lowlands in the eastern portion of the islands, and highlands around the Cordillera Central, Sierra Baoruco, and Sierra de Neyba (see Table S1 for spatial distribution of samples). Genomic DNA was extracted from tail tips or muscle tissue preserved in Ethanol using a commercial extraction kit by Eurofins Genomics (Ebersberg, Germany), and DNA yield was checked with a Qubit Fluorometer. Yield was normalized to 20 ng/μl. Illumina RAD sequencing, a proven method to detect adaptive divergence across populations (Hohenlohe, Bassham, Currey, & Cresko, 2012;Hohenlohe et al., 2010), was performed by Floragenex (Oregon). Genomic DNA was digested with Sbfl, and each sample was barcoded with specific base tags to allow multiplexing. The resulting DNA fragments library was randomly sheared and range-selected for fragments of 200-400 bp, to obtain an average final size of the library of 392-bp fragments, which were then sequenced on an Illumina HiSeq 2000 using single-end 100-bp chemistry. We verified the extent to which Sbfl enzyme digestion could provide an unbiased reduced representation of an Anolis genome by annotating Sbfl restriction sites along the related species Anolis carolinensis' genome using an in silico restriction experiment in CLC Genomics Workbench 8.0.3 (https:// www.qiagenbioinformatics.com). Results indicated that this enzyme can provide a reduced representation of an Anolis genome (total restriction sites = 17,183) without any obvious regional bias throughout the genome ( Figure S1). Quality-trimmed Illumina reads were then aligned F I G U R E 1 The study species, Anolis cybotes. Female from the Barahona peninsula on the island of Hispaniola. Copyright KCWV to a de novo assembled genome of A. cybotes. Average sample coverage was 49.11X ± 24.27 (mean:SD). SNPs with coverage >10 reads were thus genotyped from this dataset using STACKS (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013) and later filtered using VCFtools (Danecek et al., 2011). The filtering included all SNPs genotyped in >95% of the samples and excluded those with more than two alleles and rare allele frequency <1%. Population structure was determined using three alternative methods. First, we used fastSTRUCTURE (Raj, Stephens, & Pritchard, 2014), exploring a range of K values (2-36) with default settings (simple prior) and later refining the search using logistic priors. Second, we used the snmf function in the R-package LEA (Landscape and Ecological Associations, Frichot & Francois, 2015) to calculate estimates of individual admixture coefficients over a range of K values (1-16, parameter restrictions adjusted according to fast-STRUCTURE results). We determined the number of ancestral populations by comparing the cross-entropy values for each K. Finally, we also ran principal component analyses, as implemented in LEA and the Rpackage adegenet (Jombart & Ahmed, 2011), to visualize the variation in allelic frequencies between the samples and identify population groups.
We then used latent factor mixed models (LFMM, Frichot, Schoville, Bouchard, & Francois, 2013) to detect relationships between allele frequencies and environmental values, which is a method that takes into account and corrects for underlying population structure (Frichot et al., 2013;de Villemereuil, Frichot, Bazin, Francois, & Gaggiotti, 2014). To summarize the climatic variation across populations, we ran a principal component analysis on all Bioclim variables (bio1-bio19, from Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), combined with elevation of the sampling localities and extracted five environmental principal components (Table S2). Using the first principal component that had the highest proportion of explained variance, we performed five different runs of the LFMM algorithm using 100,000 iterations and 50% burn-in, setting the number of populations as latent factors (k = 6, see Section 3). We performed five different runs of the LFMM algorithm using 100,000 iterations and 50% burn-in each. The results of the five runs were combined by calculating the median z-scores and re-adjusted p-values (using the Fisher-Stouffer procedure, according to Frichot & Francois, 2015). A list of candidate SNPs associated with the bioclimatic variable was obtained after controlling for false positives with the Benjamini-Hochberg algorithm (Benjamini & Hochberg, 1995) considering the number of tested loci (L = 13,652) and 0.001 false detection ratio (FDR) as implemented in LEA package.
Sequence reads corresponding to the significant outlier SNPs were aligned to the genome of the related green anole, A. carolinensis (Genome assembly AnoCar 2.0), using BLAT in ENSEMBL (V.83) and the tblastx algorithm in the software Blast2GO (Conesa et al., 2005). Fragments with significant hits in the Anolis genome were then searched with RepeatMasker (Smit, Hubley, & Greem, 2013 to exclude transposable elements. In summary, this analysis retained 14 protein-coding outlier genes. In order to visualize covariation of allele frequencies and variation in phenotype, phenotypic measurements of A. cybotes that are known to have produced similar phenotypic adaptations across to similar environments were obtained from the data of Wollenberg et al., 2013. These data included two variables, (1) as a measure of body condition, SVL/field weight in g, and (2) XPC1, the first principal component obtained from a PCA of osteomorphological measurements. This principal component explained 43% of osteomorphological variance and included the quantitative measurements: hindlimb metatarsal length, tibia length, fibula length, femur length, ulna length, radius length (see Wollenberg et al., 2013). The factor loadings of these original variables to the XPC1 principal component variable were negative, so results have to be interpreted accordingly. We correlated body condition values for 40 and XPC1 values of 42 of the sequenced A. cybotes, and frequencies of rare alleles with elevation (in STATISTICA V.13; Dell Inc.).

| Vertebrate interactome and GO network generation
In order to infer whether newly identified protein-coding outlier genes in A. cybotes fall into previously published functional categories for thermal adaptation in vertebrates (Porcelli et al., 2015;de Villemereuil et al., 2014;Wollenberg Valero et al., 2014), we constructed a functional genomic network. As protein-protein interaction data for vertebrates are most complete for humans, we used human interaction data to explore these functional categories based on the assumption of functional similarity of orthologs. Support for this approach comes from the fact that orthologs across taxa are, despite their phylogenetic definition, most commonly identified by conserved function (Altenhoff & Dessimoz, 2009). Corresponding human gene symbols of A. carolinensis gene symbols to which A. cybotes outlier loci could be mapped were subsequently used to construct a functional genomic network in CYTOSCAPE (V.3; Shannon et al., 2003). Here, we used the STRING and UniProt databases to model human gene interactions, using "human" as a taxonomy filter.
The gene coding for Apolipoprotein B, APOB, was found as interaction partner of seven of the 14 protein-coding outlier genes after a first trial run of the analysis and has been published to be associated with cold adaptation (Gracey et al., 2004). We therefore repeated the network with the outlier genes plus APOB. To infer to which extent outlier loci among A. cybotes populations inhabiting different elevations have similar functions as predicted for vertebrate thermal adaptation candidate genes, we then constructed a network (with Kappa 0.3) of gene ontologies and biological functions using ClueGO (V. 3.0) in CYTOSCAPE. This resulted in a list of significantly overrepresented GO terms (after termwide p-value Bonferroni correction). We compared significantly overrepresented functions of the interactome, and functions of outlier genes, to candidate functions for thermal adaptation based on Porcelli et al. (2015), and Wollenberg Valero et al. (2014). Additionally, we annotated functions of A. cybotes outlier loci using information from online databases (e.g., www.genecards.org, OMIM) and literature searches. In order to explore nonrandom characteristics of the A. cybotes interactome, we also compared it to randomizations of itself; 100 randomizations were created using the NetworkRandomizer plugin to CYTOSCAPE (V.1.1.1., Aug. 2016). These randomizations were based on the Barabási-Albert model, which generates scale-free networks of a structure that is frequently encountered in nature (Albert & Barabási, 2002). The initial node degree for the algorithm was obtained from the A. cybotes interactome (m = 4). Network statistics (t tests) were then computed for differences between randomized and observed interactomes.
In order to test whether recovered candidate markers for thermal adaptation are represented with significantly higher proportion in the interactome than expected by chance, we also generated 50 random samples of 15 protein-coding genes drawn from the ENSEMBL A. carolinensis genome annotation using the UCSC genome browser interface (accessed 9/29/2016). Interactomes and basic statistics were computed for each list of 15 markers and the number of candidate markers for thermal adaptation. A. cybotes outlier gene symbols were counted within these networks, and the significance of results was assessed with t tests.

| Gene expression: tissue specificity in different vertebrates
The comparative functional genetic approach in our study led us to question, in which tissues the loci recovered in the lizard outlier + candidate gene network are expressed. Shared orthology of genes may predict similarity in tissue expression across vertebrate lineages. For example, Chan et al. (2009) verified strong evolutionary conservatism of tissue expression across vertebrates (pufferfish, frog, chicken, mouse, and human). On the same notion, we obtained gene expression data for outlier and candidate genes from a set of vertebrates (frog, lizard, zebrafish, chicken, and human). Expression data for human were again more exhaustive than data for other vertebrates, on which we consequently focus our interpretation. Human gene expression data measured on Affymetrix chip U133 were obtained for a list of 31 candidate genes plus 14 outlier genes from BioGPS (Wu et al., 2009) and imported into TABLEAU (V.9.0, Seattle, WA, USA) after data transformation (Isokpehi et al., 2015), and standardization to mean and standard deviation. Probe sets were compared for expression outlier values in different tissues. For this purpose, we categorized the 84 human tissue types for which gene expression is available in BioGPS into 11 tissue categories. We constructed a second expression dataset for the other vertebrates using the NCBI EST database to infer gene expression in 22 different tissue types.

| Anolis cybotes exhibits genomic changes in response to temperature and elevation, independently of population structure
We generated sequences from 87 samples meeting the gDNA quality threshold. Illumina Hiseq2000 sequencing was performed with F I G U R E 2 Sampling localities and genetic structure of Anolis cybotes. Top: Results from phylogenetic PCA and snmf-LEA recovering six genomic populations (red-Matadero de Honduras; yellow-Sierra de Neyba; green, northern part of Sierra Baoruco; cyan-southern part and top of Sierra Baoruco; navy-eastern lowlands; pink-Cordillera Central). Center: Distribution of samples belonging to these populations across Hispaniolan mountains and lowlands. Bottom: Histogram showing the elevational distribution of the six populations and number of sequenced samples a sequence threshold of 2,000,000 reads/sample, and the obtained mean number of reads per sample was 2,502,196. In total, the dataset contained 217,691,060 sequence reads. From this dataset, we geno-  of driest quarter (bio 9), mean temperature of warmest (bio 10), and mean temperature of coldest quarters (bio 11, Table S2).
The LFMM analysis recovered 84 outlier SNPs with convergent bioclimatic/elevational adaptation across the genetic populations. Of these, 33 SNPs were located within transposable elements, and 29 were either located in intergenic regions or could not be annotated.
The remaining 25 SNPs were located in spliced introns, transcribed introns (evidenced by the presence of cDNA in ENSEMBL), and exons of 14 protein-coding genes located on different chromosomes (see information provided in Table 1). Regardless of the annotation of SNPs, equal numbers of SNPs (between one and four) were contained within single RAD fragments. Of the SNPs located within the 14 protein-coding genes, seven SNPs were located in exons of five genes, eight SNPs were located within the introns of five genes, three SNPs were located in Intron/Exon boundaries of two protein-coding genes, and four SNPs were located in Intron/Exon boundaries of transcribed Introns (Table 1). Across all outlier SNP's, rare allele frequencies significantly increased with elevation (r = .287; p = .0186, Figure 3), as did XPC1 representing the inverse of relative bone length (r = .4337, p = .0041), and body condition (r = .52, p = .0007, Figure 3).

| The outlier interactome includes both predicted candidate genes and candidate functions for vertebrate thermal adaptation
We screened the network constructed in CYTOSCAPE (Figure 4) based on human orthologues of A. cybotes outlier loci for candidate genes previously identified in vertebrate thermal adaptation. We re- interactome (Ozgur, Vu, Erkan, & Radev, 2008), but have to be considered relative to other networks. In Wollenberg Valero et al. (2014), the same statistics were estimated for the 48-candidate genes network predicted to be conserved in vertebrate thermal adaptation.
Comparing parameters between these two networks reveals that the clustering coefficient constructed with A. cybotes outlier loci + APOB is lower than the predicted candidate gene network (which had a clustering coefficient of 0.02), which means that there are also discrete functional pathways contained within the A. cybotes outlier + candidate gene network. The network density was lower than in the original candidate gene network (which had a density of 0.005), which hints at a lower degree of common functions in the outlier + candidate gene network than in the candidate network. This is not surprising, as the RAD sequencing technique was not geared toward recovering specific genes from the original interactome and thus was more likely to reveal more distant interactions. Network heterogeneity, in contrast, was higher than in the network of 48 thermal adaptation candidate genes Functional annotation of the outlier loci (Table 1) revealed a diverse array of gene functions, with clinical significance in humans in diseases of the brain, bones, and cancer. We allocated these genes to previously predicted functions of thermal adaptation in Table 2.
Functions retrieved from online resources, and the literature (given as   A-Predicted functions of genes underlying climate adaptation. B-Thermal adaptation candidate gene list. C-Genes verified in Anolis cybotes interactome. In bold: outlier loci; In regular: predicted candidate genes; In parentheses: novel genes with similar function as predicted genes. Table References (

| Tissue-specific gene expression of thermal adaptation genes across vertebrates shows common characteristics
In humans, a majority of candidate + outlier genes in this climate adaptation interactome is expressed within different parts of the brain, the endocrine glands, and the developing fetus ( Figure 6). Specifically, high expression values are found in the pineal gland at day and night, the pituitary gland, and in the thyroid (annotated as Endocrine system in Figure 6). Further high expression values were found in developing tissues, such as fetal brain, fetal liver, and fetal thyroid (Figure 7).
Pooled gene expression data for Xenopus laevis, Xenopus tropicalis, A. carolinensis, Danio rerio, and Gallus gallus showed similar high expression values in the brain (Figure 8), but low expression values in bone. However, available data for tissue-specific expression of these genes in other vertebrates are not exhaustive yet, and information for developmental stages in this database may be imprecise, so that these results need to be interpreted with some caution.

| Functions of Anolis cybotes outlier loci
Our dataset combined samples of A. cybotes from three montane gradients and lowlands with similar bioclimatic niches. These samples belonged to six genomic populations corresponding to clades previously identified using mtDNA (Glor et al., 2003;Wollenberg et al., 2013); 14 protein-coding genes were found to contain outlier SNPs that could be adaptive to high elevation climates in the Dominican Republic. It is also likely that in addition to the coding genes, some of the intergenic SNPs we found to be under selection could also affect the expression of coding genes. These loci could additionally be affected by Linkage, which in turn may be functionally associated with thermal adaptation (see Hoffmann & Rieseberg, 2008). However, lacking both a properly annotated genome of the target species and a reference databases for gene regulatory element functional pathways, it is not yet possible to discern these effects yet and we therefore focus this manuscript on the functional interpretation of obtained intragenic SNPs and functions of their loci. Several outliers were located within genes involved in bone formation and development: SOX9, CALCR, and SPTBN1 are known to directly influence vertebrate skeletal morphogenesis, cartilage formation, and ossification (Estrada et al., 2012;Hittmeier, Grapes, Lensing, Rothschild, & Stahl, 2006). CALCR (the calcitonin receptor) has been associated with osteoporosis in humans (Table 1). It is also a biomarker for inflammation and stress response (Becker, Nylen, White, Muller, & Snider, 2004) and is overexpressed during retinal aging (Chen, Muckersie, Forrester, & Yu, 2010). The outlier locus SIRT6 (Sirtuin-6) has been identified as a biomarker for mammalian aging with DNA repair function (Mostoslavsky et al., 2006) and functions as a gatekeeper for the glucose metabolism under conditions of hypoxia (Zhong & Mostoslavsky, 2010). Locus SOD2 (Superoxide Dismutase 2) is F I G U R E 5 Predicted and recovered thermal adaptation functions. Matches between statistically overrepresented Gene Ontologies of Anolis cybotes outlier interactome presented in Figure 3 that correspond to predicted functions of climate adaptation (see Table 2). Of 120 significantly overrepresented functions in the network, 103 correspond to a predicted function for thermal adaptation. First column shows which genes or combination of genes (candidate genes only, or A. cybotes outlier genes) of the interactome these correspond to. Colors in table represent term p-values after Bonferroni stepdown correction, complete table with GOs is available as  Table S3 normally associated with the mitochondrial oxidative stress response, which has been identified as one possible factor for aging, and is a predicted function for thermal adaptation (Harman, 1956;Lund, Chu, Miller, & Heistad, 2009). Like other genes in this study, it shows high expression levels in the endocrine system, reproductive tract, brain, and during development. SOX9 (SRY-Box9) defects in humans can lead to skeletal malformations in early development and is involved in vertebrate tail regeneration (including Anolis, regulated by TGF-β, Lozito & Tuan, 2015) and temperature-dependent sex determination (Western, Harry, Graves, & Sinclair, 1999).

| Conserved functional genetic pathways for thermal adaptation and tissue specificity
Our sequencing strategy based on RAD sequencing constitutes a re- "Aging and DNA Repair" was a novel, not previously predicted function of outlier genes, and might be related to the predicted "Stress Response"/"Oxidative Stress Response" categories. The functional category "Lipid metabolism" is of known relevance for endotherm vertebrate thermogenesis, but in our study was covered by three A. cybotes outlier loci (SOX9, CALCR, MROH7). While ectotherms do not have BATassociated thermogenesis, these genes are also linked to hypoxia tolerance. This suggests a functional connection between lipid metabolism and hypoxia response. Hypercholesterolemia was shown to increase hypoxia tolerance in rats (Xi et al., 2006), but its link to thermal adaptation of ectotherms is not well studied and warrants further attention.
Independent support for our findings of conserved functional categories for reptile thermal adaptation comes from a recent study of another anole species, A. carolinensis, which recovered a genomic outlier locus for latitudinal climate adaptation with one of our predicted candidate functions "vasodilation, constriction and regulation of blood pressure" (Campbell-Staton, Edwards, & Losos, 2016).
Comparisons of the predicted thermal adaptation network (from Wollenberg Valero et al., 2014) with the A. cybotes outlier network revealed comparable properties of information content and information flow, but less common functions of network nodes. This may be an outcome of the random nature of the RAD sequencing strategy not being targeted to capture those specific genes predicted to be adaptive to thermal selection. Network randomizations revealed that the network structure of the A. cybotes interactome indicates more discrete functional pathways, a lower speed of information flow, and a more centralized structure with a higher tendency to contain hub nodes than the structure of random networks drawn from the same set of genes. This result also agrees with our hypothesis that discrete functional categories mediate thermal adaptation on a genomic level.
Furthermore, comparisons of the A. cybotes interactome with randomly generated networks indicate that the high number of previously predicted genes recovered in this network is unlikely to be a random event.
We found gene orthologs of the A. cybotes interactome in humans to be highly expressed within the brain and endocrine system. Like in mammals, the endocrine system of reptiles and amphibians is involved in various aspects of bio-regulation, such as temperature, oxygen, and humidity homeostasis. Additionally, it regulates molting and metamorphosis. In reptiles, thyroid removal reduces rates of oxygen consumption (Norris & Carr, 2013 Highest expression values across these species are found in digestive tract, brain, and head. Note that tissue-specific gene expression data are incomplete for most species. "Int and Ext Developmental" refer to internal and external "fetal" or "embryonal" tissues which, in absence of more specific information, are summarized under these categories for this manuscript. Synapomorphy-genes expressed in a structure not unambiguously comparable among species, for example, "tadpole" or "beak"

| Relationship between genotype, phenotype, and environment
Among the predicted biological functions of outlier loci, response to temperature and stress have previously been linked to phenotypic modification of A. cybotes: Muñoz et al. (2014) already suggested that A. cybotes is under greater selective pressure to adapt to cold temperatures at high elevations than to hot temperatures at low elevations. Across A. cybotes populations, average frequency of rare alleles of the intragenic SNPs increased with elevation, along with increasing body mass and decreasing length of limbs. Body condition is related to body mass (Green, 2001), which may be determined by water content, bone mass and density, or adipose tissue mass.
Despite the ability of this lizard species to behaviorally compensate for temperature fluctuation, we can show here that its phenotypic divergence is aligned with genomic divergence in response to elevation-related change in climate. We suggest that the observed phenotypic variation may have a functional basis in adaptive evolution of genes involved in bone formation and adipose tissue content.
Abiotic factors other than temperature, such as precipitation and relative humidity, also co-vary with elevation. While temperature is the abiotic factor most likely shaping reptile diversity across global elevational gradients (McCain, 2010), the challenge remains to disentangle the effect of these different abiotic climate variables on organismal evolution.
The fact that A. cybotes shows threefold recurrent similar phenotypic and ecological adaptations in different populations across the Cordillera Central, Sierra Baoruco, and Sierra de Neyba, and that repeated genomic changes arise in genes with similar functions across this dataset, can serve as first evidence from natural populations for the in silico derived hypothesis that functional pathways mediating environmental adaptation across vertebrates may be conserved (Porcelli et al., 2015;Wollenberg Valero et al., 2014).

| CONCLUSION
In conclusion, our results show that genes for vertebrate thermal adaptation can be functionally classified and seem to be conserved through evolution. While the finding of genomic in concert with phenotypic adaptation across several genetically distant populations inhabiting similar environmental gradients is indicative of recurrent adaptive responses to similar environmental selection, it does not necessarily imply convergent evolution. More research is needed to determine whether the observed pattern is an outcome of convergent evolution, of parallel evolution, or the result of recurrent selection on standing variation.
We have identified signatures of selection in several genes across three elevational gradients in a lizard. The gene functions involved suggest a connection between long-term evolutionary adaptation to an environmental selective agent, and short-term stress response in individual organisms, which emerges as an interesting novel perspective for further study.

ACKNOWLEDGMENTS
Miguel Vences (Technical University of Braunschweig, Germany) kindly provided bioinformatics resources. Jonathan Losos (Harvard University) is thanked for providing resources and helpful discussions.
The MCZ Herpetology Collection and Cryogenic Collection is acknowledged for loan of material. We thank Michael Veith and Stefan Lötters (University of Trier, Germany), for providing logistic assistance with sample shipping and quantitation. We thank Raphael D. Isokpehi (Bethune-Cookman University, FL, USA) for helpful discussions.

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
KWV conceived and designed the study and collected samples. AR and KWV performed the laboratory work, analyzed the data and wrote the manuscript. TR, RH, LH, DJ contributed to data analysis and wrote the manuscript. All authors approved of the final version.