Ecology and genomic background shape the probability of parallel adaptation to climate


 Evolution can repeat itself, resulting in parallel adaptations in independent lineages occupying similar environments. Moreover, parallel evolution sometimes, but not always, uses the same genes. Two main hypotheses have been put forth to explain the probability and extent of parallel evolution. First, parallel evolution is more likely when shared ecologies result in similar patterns of natural selection in different taxa. Second, parallelism is more likely when genomes are similar, because of shared standing variation and similar mutational effects in closely related genomes. Here we combine ecological, genomic, experimental, and phenotypic data with randomization tests and Bayesian modeling to quantify the degree of parallelism and study its relationship with ecology and genetics. Our results show that the probability of parallel adaptation to climate among species of Timema stick insects is shaped collectively by shared ecology and genomic background. Specifically, the probability of genetic parallelism decays with divergence in climatic (i.e., ecological) conditions and genomic similarity. Moreover, we find that climate-associated loci are likely subject to selection in a field experiment, overlap with genetic regions associated with cuticular hydrocarbon traits, and are not strongly shaped by introgression between species. Our findings shed light on when evolution is most expected to repeat itself.


INTRODUCTION 39
To what extent is evolution predictable and repeatable? This was the question posed by Stephen 40 J. Gould's famous thought experiment on whether repeatedly 'replaying the tape of life' would 41 yield similar evolutionary outcomes (Gould 1990) . Gould  forecasting organismal responses to natural and human-induced environmental change 51 (Waldvogel et al. 2020), the planning of plant and animal breeding programs, and the design of 52 medicines and strategies to combat the spread of disease (Lieberman et al. 2011). 53 It is now known that evolution can repeat itself but does not always do so (Grant and Grant 2002; 54 Bolnick et al. 2018). Moreover, parallelism has been documented at the genetic level, with 55 striking cases of parallel evolution at single genes of major effect, at both the within-and among-56 species level. For example, the Ectodysplasin gene controlling body armor has repeatedly been 57 used by numerous populations of stickleback fish during freshwater adaptation (Colosimo et al. 58 2005). Likewise, the Agouti and Mc1R genes control coloration in diverse organisms ( shown that parallel phenotypic evolution often involves the same genomic region (Elmer and 61 In turn, we use a field experiment and genetic mapping to bolster the evidence that climate-122 associated loci are likely subject to selection and to identify some of the traits involved. Finally, 123 we conduct genomic analyses that test the role of evolutionary history, specifically gene flow 124 and introgression, in observed parallelism. The collective results yield a comprehensive 125 evaluation of parallel evolution at the genome-wide level, in the context of an environmental 126 pressure of high current interest (i.e., climate), and in a system where comparison can be made to 127 parallelism seen at a single, major locus (i.e., Mel-Stripe). suggesting a diversity of mechanisms for repeated evolution. 134

RESULTS 135
Climatic variation within-and among-species. We studied eight Timema species across 53 136 geographic localities (n = 1420 individuals). Due to high correlations among the 22 studied 137 climate variables, we performed an ordination analysis using principal component analysis 138 (PCA) of the 22 climate variables for all populations included in the study (see Figure 2A for 139 species range map). This revealed that most of the variation in climate variables was explained 140 by the first three principal components (PC) (Total = 92.2%, PC1 = 51.7%, PC2 = 24.4% and 141 PC3 = 16.1%), which we hereafter focus on and refer to as PC1, PC2, and PC3 (Table S2 for PC 142 loadings, Figure S1). 143 The PC loadings allowed interpretation of the main factors contributing to each PC (variables 144 with PC loadings > 0.25 were considered top variables). For PC1, the top four variables were (i) 145 elevation (Ele), (ii) precipitation of warmest quarter (BIO18), (iii) precipitation of driest quarter 146 (BIO17), and (iv) precipitation of driest month (BIO14) (Figure S1A, Figure S1C, Table S2). 147 Therefore, PC1 is a general axis of elevation and precipitation variation, with high values 148 representing wet localities at high elevation (and conversely, low values representing low-149 elevation dry sites). For PC2, the top four explanatory variables were (i) mean temperature of 150 warmest quarter (BIO10), (ii) maximum temperature of warmest quarter (BIO5), (iii) mean 151 temperature of driest quarter (BIO9), and (iv) annual mean temperature (BIO1) ( Figure S1A, 152 Figure S1B, Table S2). Therefore, PC2 is a general axis of temperature variation, with high 153 loadings representing localities experiencing high temperatures. Lastly, for PC3, the top four 154 explanatory variables were: (i) precipitation of wettest quarter (BIO16), (ii) precipitation of 155 wettest month (BIO13), (iii) precipitation of coldest quarter (BIO19), and (iv) mean temperature 156 of driest month (BIO9) (Figure S1B, Figure S1C, Table S2). Therefore, PC3 is an axis of 157 contrasting variation in precipitation and temperature, with high values representing localities 158 closer to the coast experiencing greater temperature and precipitation fluctuations. 159 Overall, precipitation variables consistently loaded strongly on the first 3 PCs, followed by 160 temperature and elevation. One way ANOVA revealed significant between-species variation for 161 all three PCs (PC1: F value 104.5, P-value = < 0.0001; PC2: F value 104.5, P-value = < 0.001; 162 PC2: F value 6.803, P-value = < 0.0001; PC3: F value 28.07, P-value = < 0.0001). We also 163 detected clear within-species variation (range of median PC scores values across the eight 164 species were -3 and 5.8 for PC1, -2.5 to 6.5 for PC2, and -1.6 to 3.5 for PC3; Figure 2C-D). We 165 draw attention particularly now to PC3, as several of the most compelling results that follow 166 concern this PC. We describe a summary of our analyses in Figure 1. 167

Identifying climate-associated genomic regions within species 168
We first identified the genomic regions most strongly associated with climatic variation within 169 each of the eight species. To do so, we analyzed single nucleotide polymorphisms (SNPs) 170 obtained through previous genotyping-by-sequencing of natural populations (Riesch et al. 2017 quantile ( climate-associated SNP windows" hereafter). By using this approach, we incorporated 187 multiple loci spread across the genome in our analyses and avoided focusing on a single region 188 of the genome. Additionally, we expected this range of quantile to be enriched for loci truly 189 involved in climate adaptation. Finally, this approach allowed us to make comparisons among 190 species in downstream analysis. In all species, the top climate-associated SNP windows were 191 widely distributed across the genome and found on all 13 linkage groups (LGs) (Figure 3,Figure 192 S2, Figure S3). 193

Parallel evolution of climate-associated genomic regions across species 194
We next quantified the extent to which climate-associated SNP windows were parallel (i.e., the 195 same) across the eight species of Timema that we studied. We did so by testing if such windows 196 exhibit excess overlap across species relative to that expected by chance. We used two related 197 approaches. First, we asked if climate-associated SNP windows exhibited excess overlap 198 between each pair of species ( pairwise-comparison", Figure 1B). Second, we asked if climate-199 associated SNP windows exhibited excess overlap among multiple species i.e., if the same SNP 200 windows show association with climate PCs between 3, 4, 5, 6, 7 or 8 species ( multi-species 201 comparison", Figure 1B). To do so, we conducted randomisation tests to quantify excess overlap 202 of windows relative to expectations assuming that SNP-climate associations were independent 203 across species, for pairwise as well as multi-species comparisons ( Figure 1A). As an example, an 204 x-fold enrichment of 2.0 in the parallelism analyses would indicate that the evidence for overlap 205 of climate-associated SNP windows for a given comparison was two times higher than expected 206 by chance (based on the mean of the null). For this, we focused on windows with the greatest 207 (top 10%) climate association in nature (for all three climate PCs.) 208 These analyses revealed compelling evidence for parallel climate-associated SNP windows 209 across species. For pairwise comparisons, we detected significant excess overlap in climate-210 associated SNP windows for all three climate PCs for at least some species pairs ( Figure 4A, 211 Figure S4A, and Figure S5A). For example, we observed excess overlap between eight species 212 pairs for PC1 ( Figure S4B), one species pair for PC2 ( Figure S5B), and five species pairs for 213 PC3 ( Figure 4B). The quantitative degree of excess overlap was somewhat variable. For 214 example, for PC1, the overlap of climate-associated SNP windows among the seven species pairs 215 that showed a significant excess of overlap ranged from an x-fold enrichment of ~1.3x for T. 216 knulli and T. poppensis (observed = 36, expected = 27, x-fold enrichment = 1.31, P-value = 0.04) 217 to a maximum x-fold enrichment of ~2x for T. podura and T. chumash (observed = 23, expected 218 = 11, x-fold enrichment = 2.20, P-value = 0.01) ( Figure S4A-B). For PC2, there was only one 219 species pair, T. bartmani and T.californicum, which showed significant excess overlap of SNP 220 windows with largest median Bayes factors (observed = 29, expected = 19, x-fold enrichment = 221 1.48, P-value <.01; Figure S5B). Lastly, for PC3, the overlap of climate-associated SNP 222 windows among the five species pairs that showed a significant excess of overlap ranged from an 223 x-fold enrichment of ~1.2x for T. cristinae and T. landelsensis (observed = 49, expected = 38, x-224 fold enrichment = 1.29, P-value = 0.03) to ~1.5x for T. bartmani and T. cristinae (observed = 27, 225 expected = 17, x-fold enrichment = 1.53, P-value = 0.01) ( Figure 4A-B). 226 Parallelism was also supported in analyses examining multi-species (i.e., beyond pairwise) 227 comparisons. For PC1, there was significant excess overlap of SNP windows with largest median 228 Bayes factors between three or more species ~2x more than expected by chance (observed = 60, 229 expected = 26.77, x-fold enrichment = 2.25, P-value < 0.01; Figure S4C), for four or more 230 species about ~3x more than expected by chance (observed = 4, expected = 1.03, x-fold 231 enrichment = 3.87, P-value 0.02; Figure S4C). For PC2, there was significant excess overlap of 232 SNP windows with largest median Bayes factors between three or more species about ~1.5x 233 more than expected by chance (observed = 42, expected = 26.41, x-fold enrichment = 1.59, P-234 value <.01; Figure S5C), for four or more species about ~4x more than expected by chance 235 (observed = 5, expected = 1.19, x-fold enrichment = 4.17, P-value = 0.007; Figure S5C). Lastly, 236 for PC3 there was significant excess overlap of climate-associated SNP windows between three 237 or more species ~1.6x more than expected by chance (observed = 43, expected = 26, x-fold 238 enrichment = 1.63, P-value < 0.01; Figure 4C). This effect was even stronger for shared 239 windows among four or more species (observed = 5, expected = 1.10, x-fold enrichment = 4.53, 240 P-value = 0.006; Figure 5C). Similar results were observed for PC1 and PC2 ( Figure S4C and 241 Figure S5C). After finding this evidence for parallelism, we turned to testing our core two 242 hypotheses, by quantifying the extent to which parallelism decays with ecological divergence, 243 genomic divergence, or both. 244

Parallelism declines predictably with ecological and genomic divergence between species 245
We next tested the shared ecology and shared genetics hypotheses for variation in the degree of 246 parallelism. Shared ecology would cause a higher degree of parallelism due to similar selective 247 pressures from similar environments. One the other hand, shared genetics would a higher degree 248 of parallelism due to higher probability of gene reuse. Here, we quantified parallelism as the 249 degree of excess overlap of climate-associated SNP windows relative to null expectations, for 250 pairwise comparisons ( Figure 2B). We estimated ecological similarity between pairs of species 251 using climatic data and genetic similarity based on a previously published phylogeny. We then 252 fit Bayesian linear mixed models to explicitly compare models where the degree of parallelism is 253 determined by ecological similarity, genetic similarity, and both. Notably, this mixed model 254 approach accounts for the non-independence of pairwise distances (Gompert et al. 2014 for 255 details). Specifically, for each climatic PC variable, we modeled parallelism as the x-fold excess 256 in shared top climate-associated SNP windows (as described in the preceding section) as a 257 function of climatic distance which was calculated as the average difference in climate PC scores 258 between a given pair of species (hereafter referred to as ecology indicating divergence in 259 ecology"), phylogenetic distance (hereafter referred to as genes indicating divergence in 260 genetics"), or both. We tested the effect of both ecology and genetics on parallelism ( Figure 5B, 261 Figure S6B, Figure S7B). The fit of the different models with or without ecology or genetics was 262 compared using deviance information criterion (DIC), which is a metric of predictive 263 performance. 264 Our analyses revealed evidence for effects of both ecology and genes on the probability of 265 parallelism, with results that varied among the climate PCs ( Figure 5C-D for PC3, Figure S6C-D 266 for PC1, Figure S7C-D for PC2). For PC3, the full model (ecology and genes) was the best 267 model, with similar, negative effects on parallelism of divergence in ecology (standardized beta 268 0.87 to -0.21; Figure 5E; Table S3). For PC1, the genes-only model was the best model 270 (standardized beta = -0.55, 95% CI = -0.8 to -0.25; Figure S6E, Table S3). The second-best 271 model was the full model, but this included a positive rather than negative effect of ecological 272 distance on parallelism. Lastly, for PC2 the best model was a null model of no effect of genes or 273 ecology on parallelism ( Figure S7E, Table S9). The results thus provide variable support for both 274 the shared ecology and shared genetics hypotheses, dependent on climate PC and strongest for 275 PC3. Having tested these hypotheses, we next tested for additional evidence, beyond parallelism, 276 that the climate-associated SNP windows have been affected by natural selection. 277 Climate-associated regions exhibit elevation-dependent allele-frequency change in a field 278 experiment and co-vary with CHCs in T. cristinae 279 To potentially bolster the evidence that climate-associated SNP windows are enriched for genetic 280 variants experiencing natural selection, we next tested if such windows exhibited exceptional 281 patterns of allele-frequency change in a transplant-and-sequence field experiment and if they 282 overlap with regions associated with phenotypic variation in genetic mapping analyses. Here, an 283 x-fold enrichment of 2.0 in the CHC analysis would indicate that the evidence for SNP-CHC 284 associations in climate-associated SNP windows was two times higher than expected by chance 285 (based on the mean of the null, see below for details). 286 The transplant experiment measured 8-day survival and associated genome-wide allele frequency 287 change during this period in 500 T. cristinae transplanted to 10 experimental bushes comprising 288 two host plants and all occurring at a gradient of higher elevations than the source population for 289 the experiment. A previous analysis of this experiment documented evidence of selection 290 associated with elevation. Here, as a metric of possible elevation (environment)-dependent 291 selection, we calculated the Pearson correlation between transplant elevation and allele 292 frequency change caused by mortality during the experiment for each SNP. In this current 293 analysis, we found that 100 kb windows exhibiting patterns of allele frequency change most 294 strongly associated with elevation in the transplant experiment coincided modestly with climate-295 associated SNP windows, but more than expected by chance. Specifically, when focusing on the 296 windows with the greatest (top 10%) correlation between change and elevation in the experiment 297 and with the greatest (top 10%) climate association in nature, windows associated with all three 298 climate PCs corresponded with those where change was most strongly associated with elevation 299 ~1.2-1.3 times more than expected under the null hypothesis of independence (constrained 300 randomization test controlling for SNP density within windows based on 1000 randomizations; 301 PC1: observed = 108 shared windows, P = 0.005; PC2: observed = 101 shared windows, P = 302 0.015; PC3: observed = 105, P = 0.021 windows) ( Figure 6). Similar patterns were observed 303 where more extreme top percentiles were considered, and when using an unconstrained 304 randomization test (Table S4). These patterns are consistent with the hypothesis that multiple 305 genetic variants in these windows are subject to selection in nature. 306 For the CHC analyses, we considered three compound classes -pentacosanes, heptacosanes, and 307 nonacosanes -in males and in females (i.e., six CHC traits total). We found evidence of heritable 308 variation for each compound in both male and female T. cristinae, with 50.8% (male 309 nonacosanes) to 89.7% (female pentacosanes) of the variability in these traits explained by ~176 310 thousand sequenced SNPs in a mapping population (these values denote Bayesian point 311 estimates; these results are based on 602 T. cristinae from a single population, FHA) (see Table  312 S5 for details). We summarized the evidence that each 100 kb window included CHC-associated 313 SNPs by computing the mean posterior probability of association (i.e., the mean probability of a 314 non-zero genotype-phenotype association, also known as the posterior inclusion probability or 315 PIP) across SNPs in 100 kb windows (i.e., the same 100 kb windows used for summarizing SNP-316 climate associations). Then, based on a randomization test, we found that for some CHC traits, 317 the average posterior inclusion probability for SNPs in the top climate-associated SNP windows 318 in T. cristinae was marginally but significantly greater than expected by chance. Specifically, the 319 average probability of SNPs being associated with female pentacosanes was ~1.05 times higher 320 than expected by chance for both the top 10% of PC2 and PC3 climate-associated SNP windows 321 (P-value = 0.009 for PC2 and P-value = 0.010 for PC3 based on 1000 permutations; Figure 5B, 322 Table S7 and S8). We also detected a marginally non-significant increase in the average 323 posterior inclusion probability for SNP associations with female nonacosanes in the top 10% of 324 PC3 climate-associated SNP windows (x-fold increase in mean inclusion probability = 1.03, P-325 value = 0.051, 1000 permutations, Table S8). These results from CHCs support the hypothesis 326 that at least a subset of the top climate-associated SNP windows are associated with traits 327 involved in climatic adaptation in Timema. Thus, together with the results presented in the 328 previous paragraph, these results suggest a polygenic basis for climatic adaptation in T. cristinae 329 with at least a modest correspondence between our top climate-associated windows and the 330 actual loci involved in climate adaptation. With this evidence, we next turn to additional analyses 331 concerning the evolutionary history and mechanisms of parallelism, namely the potential role of 332 introgression between species in promoting parallelism. 333

Introgression between species does not contribute strongly to parallel evolution 334
We conducted two analyses, focused on different time scales, to ask if introgression and gene 335 flow between species promotes gene sharing and thus climate-associated parallel evolution. First, 336 we identified historical patterns of introgression using a population tree-based approach. Second, 337 we identified contemporary patterns of gene flow using an admixture model. Both these analyses 338 helped us to assess the degree of genetic independence in adaptation to climate within each 339 species. 340 To identify historical patterns of introgression, we used TREEMIX to generate a population tree 341 for all populations and species while allowing for historical admixture or gene flow among intra-342 specific or inter-specific populations. For this analysis, we realigned GBS sequence data for all 343 1420 individuals included in this study to the T. cristinae genome. We then called and filtered 344 single nucleotide polymorphisms (SNPs) to identify final set of 8787 SNPs for this analysis. Our 345 results from TREEMIX yielded a population graph or bifurcating tree depicting relationships 346 between focal localities in this study. The best bifurcating tree explained 99.6% of the variation 347 in the population allele-frequency covariances. In this tree, Timema populations formed eight 348 major clades that grouped populations by species ( Figure 5A). Adding migration edges to the 349 tree increased the variance explained by a negligible extent (Table S9) We used GBS data from 1420 individuals across eight species, combined with other forms of 366 data, to show that adaptation to climate is occurs in parallel across species but as a function of 367 the ecological and genomic divergence between species. Our results inform three fundamental 368 issues in biology, namely the repeatability of evolution, the effect of ecology and genetics on 369 parallelism, and the processes promoting parallelism. We treat these in turn. In summary, our study shows how local adaptation even among species with minimal gene flow 400 can occur and consequently be crucial for predicting evolution in response to rapidly changing 401 environments and climate. Furthermore, our results bolster evidence for selection beyond a 402 correlational genome scan because we found that the genomic regions which underlie parallelism 403 also showed marked allele-frequency change in an experiment and were associated with 404 ecologically relevant CHC traits. Thus, together these results suggest that allele reuse through 405 standing genetic variation, new mutations, and selection can be powerful drivers of local 406 adaptation. 407

METHODS 408
Below we describe details of all our methods and analyses, and we provide a graphic summary in 409  (Table S1). GBS data for 417 this study has been previously published in a study of speciation continuum in Timema (Riesch 418 et al. 2017) (also see data availability for more details). The genomic data in the transplant 419 experiment and used for genetic mapping of cuticular hydrocarbons is independent from these 420 data and is described in detail below. bases with base quality 13, and ignored insertion-deletion polymorphisms. We then set the prior 432 on single nucleotide polymorphisms (SNPs) to 0.001 (-P) and called SNPs when the posterior 433 probability that the nucleotide was invariant was 0.01 (-p). We then performed two rounds of 434 filtering to retain final sets of SNPs. In the first round, we filtered the initial set of SNPs to retain 435 only those with sequence data for at least 80% of the individuals, a mean sequence depth of two 436 per individual, at least four reads of the alternative allele, a minimum quality score of 30, a 437 minimum (overall) minor allele frequency of at least 5%, and no more than 0.01% of the reads in 438 the reverse orientation. In the second round of filtering, we removed SNPs with excessive 439 coverage (2 standard deviations above the mean) or that were tightly clustered ( Figure 2B, Table S2, Figure S1), we used these three PCs to study genomic 460 associations with climate in all further analyses. Factors for 100 kilobase (kb) non-overlapping SNP windows (i.e., the same windows were used 474 in every species, facilitating comparisons among them). Our downstream analyses described 475 below focus on these windows. We delimited climate-associated SNP windows as those with 476 greatest association with the three climate PCs, specifically as the windows in the top 10% 477 quantile. We refer to such windows as "climate-associated SNP windows" hereafter. 478

Quantifying parallel genomic associations with climate across species 479
We quantified parallel genomic associations with climate across species (using the results 480 described above from BayPass) and used randomization tests to measure the extent to which the 481 observed parallelism exceeded that expected by chance. We report this excess as 'x-fold' 482 enrichments, relative to null expectations, also reporting associated P-values for statistical 483 significance. 484 We first tested for excess overlap of climate-associated SNP windows between pairs of species 485 ( pairwise comparisons"), for each of the three climate PCs. To do this, we used randomization 486 tests (10,000 randomizations per test) to generate null expectations for the proportion of top 487 climate-associated SNP windows shared between a given pair of species and tested whether this 488 was significantly more than expected by chance (x-fold enrichments and P-values). As an 489 example, an x-fold enrichment of 2.0 would indicate that twice as many climate-associated SNP 490 windows showed overlap between a species pair than was expected by chance (based on the 491 mean of the null). We then quantified overlap in climate-associated SNP windows between 492 multiple species ( multi-species comparisons") i.e., we tested if similar SNP windows show 493 association with climate PCs between 2, 3, 4, 5, 6 or 8 species. We did this by using similar 494 randomization tests as for pairwise-comparisons, to generate x-fold enrichments and P-values. 495

Testing the shared ecology and shared genetics hypotheses 496
We tested the contribution of shared ecology versus shared genetics to the observed degree of 497 parallelism. We expect both shared ecology and genetics to have an effect on the probability of 498 parallelism. To do so, we fit Bayesian linear mixed models to explicitly compare models where 499 parallelism is determined by ecological similarity, genetic similarity, or both. This Bayesian 500 regression analysis is based on the mixed model framework proposed by (Clarke, Rothery, and 501 Raybould 2002) and extended by (Gompert, Lucas, et al. 2014). Our method accounts for the 502 correlated error structure inherent in pairwise covariates and response variables (e.g., ecological 503 or genetic distances). In this analysis, our response variable was the x-fold excess in shared top 504 climate-associated SNP windows for a given PC (we did analyses separately for each climate 505 PC). Our independent variables were ecological and genetic distances, estimated as follows. 506 Climatic (i.e., ecological) distance was calculated as pairwise absolute mean difference of PC 507 scores of each species. We calculated genetic (i.e., phylogenetic) distances based on the 508 previously published phylogeny described in Riesch et al. 2017. Briefly, we used the data from 509 this previous phylogeny constructed using Bayesian phylogenetic inference with BEAST 510 (version 2.1.387) for 11 Timema species based on GBS data of curated dataset of 19,556 single-511 nucleotide variants. For our current study, we used pairwise phylogenetic distances for the eight 512 Timema species as our metrics of genetic distances for this analysis. All variables were 513 standardized (given mean 0 and standard deviation of 1) before analysis. 514 We then considered four alternative models: (i) a null model without covariates, (ii) a model 515 including only phylogenetic distance, (iii) a model with only climatic distance, and (iv) a model 516 with both climate and phylogenetic distance. We fit the models in R using the rjags (version 4.8) 517 interface with Jags (version 4.3.0). We used minimally informative priors for the regression 518 coefficients (i.e., normal with μ = 0 and precision τ = 0.001) and for the population random 519 effects and residual errors, all gamma (1, 0.01). Deviance information criterion was used for 520 model comparison. Parameter estimates and DIC estimates were obtained via MCMC. For each 521 analysis and model, we ran three chains each comprising 10,000 sampling iterations, a 2000-522 iteration burn-in, and a thinning interval of 5. 523

Overlap of climate-associated SNP windows with genetic regions showing elevation-524 associated change in a field experiment 525
We quantified overlap between climate-associated SNP windows and windows that exhibited 526 elevation-dependent allele-frequency change in a previously published release-recapture field 527 experiment. We then tested if this overlap was greater than expected by chance. Survivors were recaptured eight days later. Whole-genome sequence data, which we analyze 533 here, was obtained from 491 of the 500 stick insects (Nosil et al. 2018). 534 For the current study, we estimated allele frequencies in the released and recaptured stick insects 535 at the 6,175,495 bi-allelic SNPs identified by Nosil et al. 2018. This was done using an 536 expectation-maximization (EM) algorithm as implemented in the program estpEM (version 0.1) 537 with tolerance of 0.001 and a maximum of 50 EM iterations (Soria-Carrasco et al. 2014). We 538 then used these estimates to compute allele-frequency change between the start and end of the 539 experiment. Then, for each SNP we calculated the Pearson correlation between allele frequency 540 change and the elevation at each of the ten transplant sites. Finally, we determined the average 541 correlation between change and elevation for the 100 kb windows across the genome. Windows 542 with fewer than four SNPs were ignored. These steps were done using R (version 3.4). 543 We then calculated the number of 100kb windows that were among the top 10% for both 544 elevation-dependent change during the experiment (highest average absolute correlation) and for 545 climate-association (highest average Bayes factor for each climate PC). We used a constrained 546 randomization procedure to generate null expectations for such concordance between change and 547 climate-association windows, using a separate randomization for each PC. Specifically, we 548 randomized mean change metrics across windows, but only among windows with similar SNP 549 densities (10 equally sized bins were used for this). This was done because we observed a 550 positive correlation between SNP density and mean change-elevation correlations per window 551 (Pearson R = 0.069, 95% CI = 0.047-0.091, P < 0.001), and we wanted to control for this. Null 552 distributions and P-values were based on 1000 randomizations and are reported for each climate 553

PC. 554
Overlap of climate-associated SNP windows with genetic regions associated with cuticular 555 hydrocarbon variation 556 Our next set of analyses concern climate-associated SNP windows and genomic regions 557 associated with cuticular hydrocarbon (CHC) variation. The logic here is that CHCs tend to play 558 a role in desiccation tolerance and climatic adaptation in insects (e.g., (Rajpurohit et al. 2017)), 559 such that genetic regions associated with climate versus CHCs might overlap. We thus 560 specifically quantified the extent to which climate-associated SNP windows overlapped with 561 windows harbouring SNPs associated with CHCs, and whether this overlap was greater than 562 expected by chance. The CHC data were originally described and analyzed by Riesch et al. 563 (2017). Specifically, for each insect we had quantified the proportional abundance of 26 different 564 mono-and di-methylated CHCs, which comprised eight pentacosanes, eight heptacosanes and 565 ten nonacosanes, and then applied log-contrasts. For the current dataset, we used those values to 566 calculate the proportional abundance of the sum of all pentacosanes, the sum of all heptacosanes 567 and the sum of all nonacosanes (henceforth: pentacosanes, heptacosanes and nonacosanes). 568 Therefore, the six CHC traits considered were pentacosanes, heptacosanes, and nonacosanes in 569 males and females (i.e., three molecule types in each of two sexes). 570 Here, we first re-aligned the GBS data from Riesch et al. (2017)  for CHC data was also collected. These data were aligned to the genome using the BWA ALN 574 algorithm (version 0.7.17-r1188) (Li 2013). We allowed for 5 miss-matches total, and not more 575 than 2 miss-matches in the first 20 bp. Only reads with a mapping quality greater than 10 were 576 retained. We then compressed, sorted and indexed the alignments with SAMTOOLS and 577 BCFTOOLS (version 1.2) (Li et al. 2009;Danecek et al. 2021). Next, we used SAMTOOLS and 578 BCFTOOLS to identify SNPs and calculate genotype likelihoods. For this, we used the 579 recommended mapping quality adjustment (-C 50), only considered alignments with mapping 580 qualities of 20 or more and SNPs with base qualities of 30 or more, and only called variants 581 when the posterior probability that locus was invariant was less than 0.01 given a prior mutation 582 rate parameter of 0.001. We then used custom perl scripts to filter out variants with a mean 583 coverage of less than 2x, fewer than 10 non-reference reads total, mapping quality less than 30, 584 minor allele frequency less than ~0.005, more than 1% of reads in the reverse orientation (with 585 our GBS method, all reads should have the same orientation), missing data (no reads) for more 586 than 20% of individuals, SNPs with more than two alleles, and SNPs with coverage exceeding 587 three standard deviations above the mean. Finally, we obtained Bayesian point estimates 588 (posterior means) of genotypes for each locus and individual based on the genotype likelihoods 589 and used the estimated allele frequencies to parameterize a binomial prior. 590 We then conducted genetic mapping of CHC variation using a polygenic genome-wide 591 association (GWA) mapping approach, that controls for linkage disequilibrium among SNPs and 592 background population structure as detailed below. We specifically fit Bayesian sparse linear 593 mixed models (BSLMMs) to determine the contribution of additive genetic variation (as 594 captured by our collective SNP data set) to each of six CHC traits, and to determine the 595 probability of association (posterior inclusion probability, PIP) of each individual SNP with each 596 trait (this PIP value is computed from, i.e., equal to, the proportion of MCMC samples that 597 included each SNP in the polygenic regression model). We fit this model using gemma (version 598 0.95a) (Zhou, Carbonetto, and Stephens 2013), a polygenic GWA mapping method that fits a 599 single model with all SNPs while accounting for uncertainty and redundancy in genotype-600 phenotype associations, for example by controlling for linkage disequilibrium among SNPs, and 601 background polygenic effects. The latter is inferred based on a kinship matrix derived from the 602 collective SNPs, which also serves to control for population structure when estimating effects for 603 individual SNPs. Models were fit using MCMC, with each mapping exercise involving 10 604 independent chains each comprising 1 million sampling iterations and a 200,000-iteration burn-605 in. 606 Based on these analyses, we then computed the mean PIP (i.e., probability of a genotype-607 phenotype association) across all SNPs in 100 kb windows for each of the six CHC traits. Then, 608 we asked whether the average association with CHCs (averaged over windows) was higher for 609 the climate-associated SNP windows than expected by chance. Randomizations (1000) were 610 used to generate a null distribution. Specifically, mean posterior probabilities for SNP-CHC 611 associations were permuted across 1000 kb windows and the number windows in the top 10% 612 for climate association and (permuted) CHC posterior inclusion probabilities was determined. 613 Note that we conducted this test independently for each of the six CHC traits and each of the 614 three climate PCs. We then examined the combination of these results to assess the total 615 evidence that SNP windows associated with climate adaptation are enriched for those regions of 616 the genome possibly affecting CHC variation. 617

Testing for introgression and quantifying population structure 618
We quantified both historical and contemporary gene flow patterns, respectively as follows. For the analysis in our study here, we re-aligned the GBS sequences for 1420 individuals (across 624 53 populations) included in this study to the T. cristinae genome (draft version 0.3). We did this 625 by using the MEM algorithm from BWA (version 0.7.17-r1188). We ran BWA MEM with a 626 minimum seed length of 15, internal seeds of longer than 20 bp, and only output alignments with 627 a quality score >= 30. We then used SAMTOOLS (version 1.6) to compress, sort and index the 628 alignments (Li et al. 2009). We then identified SNPs using SAMTOOLS and BCFTOOLS 629 (version 1.6). For variant calling, we used a mapping quality of 50, skipped alignments with 630 mapping quality lower than < 20, skipped bases with base quality <15, and ignored insertion-631 deletion polymorphisms. We set the prior on SNPs to 0.001 and called SNPs when the posterior 632 probability that the nucleotide was invariant was <=0.01. After we got the initial set of variants, 633 we filtered them to retain only those SNPs with sequence data for at least 80% individuals, a 634 mean sequence depth of two per individual, at least 4 reads of the alternative allele, a minimum 635 quality score of 30, a minimum overall) minor allele frequency of at least 0.005, and no more 636 than 1% of the reads in the reverse orientation (this is an expectation for our GBS method). We 637 further removed SNPs with excessive coverage (3 standard deviations above the mean) or that 638 were tightly clustered (within 3 bp of each other), as these could be poor alignments (e.g., reads 639 from multiple paralogs mapping to the same region of the genome). This left us with 8787 SNPs 640 for this analysis. We used these variants to run TREEMIX to fit trees allowing 0-9 admixture

914
Here red dots denote climate-associated SNP windows for each species. We then quantify overlap in these windows 915 between the two species ( N").