Cold adaptation drives population genomic divergence in the ecological specialist, Drosophila montana

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Cold adaptation drives population genomic divergence in the ecological specialist, Drosophila montana © 2021 the Authors Published version Wiberg, R. A. W.; Tyukmaeva, V.; Hoikkala, A.; Ritchie, M. G.; Kankare, M. Wiberg, R. A. W., Tyukmaeva, V., Hoikkala, A., Ritchie, M. G., & Kankare, M. (2021). Cold adaptation drives population genomic divergence in the ecological specialist, Drosophila montana. Molecular Ecology, 30(15), 3783-3796. https://doi.org/10.1111/mec.16003 2021 Molecular Ecology. 2021;30:3783–3796. | 3783wileyonlinelibrary.com/journal/mec Received: 22 April 2020 | Revised: 10 May 2021 | Accepted: 20 May 2021 DOI: 10.1111/mec.16003 ORIGINAL ARTICLE Cold adaptation drives population genomic divergence in the ecological specialist, Drosophila montana R. A. W. Wiberg1 | V. Tyukmaeva2 | A. Hoikkala2 | M. G. Ritchie1 | M. Kankare2 This is an open access article under the terms of the Creat ive Commo ns Attri bution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Molecular Ecology published by John Wiley & Sons Ltd. 1Centre for Biological Diversity, School of Biology, University of St Andrews, St Andrews, UK 2Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland Correspondence R. A. W. Wiberg, Centre for Biological Diversity, School of Biology, University of St Andrews, St Andrews, KY16 9TH, Scotland, UK. Email: raw.wibe[email protected] Present address R. A. W. Wiberg, Department of Environmental Sciences, Zoological Institute, University of Basel, Basel, Switzerland V. Tyukmaeva, Centre d'Ecologie, Fonctionelle et Evolutive, CNRS, Montpellier, France Funding information Ella ja Georg Ehrnroothin Säätiö; Academy of Finland, Grant/Award Number: 267244, 268214 and 322980; Natural Environment Research Council, Grant/Award Number: NE/L501852/1 and NE/P000592/1 Abstract Detecting signatures of ecological adaptation in comparative genomics is challenging, but analysing population samples with characterised geographic distributions, such as clinal variation, can help identify genes showing covariation with important ecological variation. Here, we analysed patterns of geographic variation in the coldadapted species Drosophila montana across phenotypes, genotypes and environmental conditions and tested for signatures of cold adaptation in population genomic divergence. We first derived the climatic variables associated with the geographic distribution of 24 populations across two continents to trace the scale of environmental variation experienced by the species, and measured variation in the cold tolerance of the flies of six populations from different geographic contexts. We then performed pooled whole genome sequencing of these six populations, and used Bayesian methods to identify SNPs where genetic differentiation is associated with both climatic variables and the population phenotypic measurements, while controlling for effects of demography and population structure. The top candidate SNPs were enriched on the X and fourth chromosomes, and they also lay near genes implicated in other studies of cold tolerance and population divergence in this species and its close relatives. We conclude that ecological adaptation has contributed to the divergence of D. montana populations throughout the genome and in particular on the X and fourth chromosomes, which also showed highest interpopulation FST. This study demonstrates that ecological selection can drive genomic divergence at different scales, from candidate genes to chromosomewide effects. KEYWORDS chill coma recovery time, cline populations, cold tolerance, CTmin, D. montana, environmental adaptation, genomic divergence 3784 | WIBERG Et al. 1 | INTRODUCTION The geographic structure of a species is a result of its phylogeographic history, influenced by past and present dispersal, population demography, and selection. Obtaining genomewide data on genetic polymorphisms across multiple populations of a species is becoming relatively easy, but interpreting the patterns of geographic variation in such data and identifying genes which vary primarily due to selection remains challenging. Often simple”outlier” approaches using genome scans which measure genetic differentiation such as FST or Dxy are adopted, but results are difficult to interpret due to confounding effects of selection, drift and population structure, or genomic features such as inversions and other causes of variation in recombination rate (Cruickshank & Hahn, 2014; Noor & Bennett, 2009; Ravinet et al., 2017; Wolf & Ellegren, 2016). If environmental data are available, we can use associations with such factors to help identify loci where genetic differentiation covaries with this environmental variation. Some genome scan methods can incorporate environmental variation and simultaneously fit effects for covariance with environmental factors, while controlling for effects of population demography (Foll & Gaggiotti, 2008; de Villemereuil & Gaggiotti, 2015). This approach has successfully identified genetic variation associated with altitude in humans, among other examples (Foll et al., 2014; Gautier, 2015; de Villemereuil & Gaggiotti, 2015) and has become a useful approach to investigate the ecological adaptations underlying population divergence. Clinal patterns of variation in phenotypes or gene frequencies have a long history of being used to infer selection along ecotones, and analyses of cline shape can sometimes identify loci under direct selection from others showing clinal variation for other reasons, such as phylogeographic history (Barton & Gale, 1993). Such studies can be very powerful, especially when independent parallel clines are available. For example, Kolaczkowski et al., (2011) sampled isofemale lines from extremes of a cline in Australian populations of D. melanogaster, and found many genes implicated in clinally varying phenotypes to show highest differentiation. Also, Bergland et al., (2014), and Kapun et al., (2016) sampled North American clines in D. melanogaster and D. simulans over several years to uncover clinal variation at a genome scale. Bergland et al., (2014) also found consistent fluctuations in allele frequencies for a population sampled over several seasons, indicating a regular response to seasonally varying selection pressures. On the other hand, Machado et al., (2015) concluded that migration and gene flow play a greater role than adaptation in the overall clinality of genomic variants in D. simulans than D. melanogaster. While the two species share a significant proportion of the genes showing clinal variation, their differences in overwintering ability, migration and population bottlenecks probably act as additional drivers of differences in patterns of variation between them (Machado et al., 2015). Similar studies of clinal variation in phenotypes and allele frequencies have also been carried out in other insects (Paolucci et al., 2016), plants (Bradbury et al., 2013; Chen et al., 2012), mammals (Carneiro et al., 2013; Hoekstra et al., 2004), fish (Vines et al., 2016), and other organisms (Endler, 1973, 1977; Takahashi, 2015). However, more studies are needed that simultaneously compare patterns of variation in allele frequencies and potentially causal environmental variation or ecologically important traits while controlling for population structure, as such studies are necessary to determine if adaptation to climate is directly driving patterns of genetic differentiation. Here, we investigated geographic variation at both the phenotypic and genetic levels in Drosophila montana samples from two continents. This species has spread around the northern hemisphere (Throckmorton, 1969) and both mtDNA and microsatellite data have revealed genetically distinct Finnish and North American populations (Mirol et al., 2007). Moreover, more recent modelling of genomewide SNP frequencies suggest that the FinnishNorth American split happened around 1.75 Mya and that of the North American populations shortly after that (Garlovsky et al., 2020). It is one of the most coldtolerant Drosophila species (Kellermann et al., 2012) and the basic cold tolerance of D. montana flies can increase towards the cold seasons through two mechanisms, photoperiodic reproductive diapause (Vesala & Hoikkala, 2011) and coldacclimation induced by a decrease in day length and/or temperature (Kauranen et al., 2019; Vesala, Salminen, Kostal, et al., 2012; Vesala, Salminen, Laiho, et al., 2012). D. montana populations have been found to show clinal variation in the critical day length required for diapause induction (CDL; Lankinen et al., 2013; Tyukmaeva et al., 2011). There is also a correlation between CDL and latitudinally covarying climatic factors such as the mean temperature of the coldest month (Tyukmaeva et al., 2020). In addition, D. montana populations from different geographic regions show variation in their courtship cues and mate choice (Klappert et al., 2007; Routtu et al., 2007), which has led to partial reproductive isolation between some distant populations (Jennings et al., 2011, 2014). At the genetic level, differential gene expression studies have identified candidate genes underlying diapause (Kankare et al., 2010, 2016), perception of day length (Parker et al., 2016), and cold acclimation (Parker et al., 2015). Furthermore, a quasinatural selection experiment for shorter CDL, accompanied by a decrease in coldtolerance, induced widespread changes in loci with potential roles with these traits (Kauranen et al., 2019). Finally, population genomic analyses have identified several outlier loci when examining differentiation between North American and European populations (Parker et al., 2018). All this makes D. montana an interesting example of nascent speciation, potentially influenced by ecological adaptation. Here, we sought to ask to what extent patterns in the genomic divergence of D. montana populations across two continents are correlated with climatic variation and phenotypic responses to cold adaptation. We performed pooled wholegenome sequencing (poolseq) on six different populations and used Bayesian methods to examine the association between genomic differentiation between populations and environmental variables across both continents. We also phenotyped populations for two different cold tolerance measures; critical thermal minimum (CTmin), and chill coma recovery time (CCRT), and investigated the associations between them and the genetic and climatic data. Ultimately, we | 3785 WIBERG Et al. asked if the genomic loci showing an association between genetic and environmental differentiation also showed association with population differentiation in cold tolerance phenotypes, and examined the possible overlap between the set of genes close to candidate SNPs with sets of candidate genes from previous studies of cold adaptation in D. montana. If population differentiation is driven by ecological selection then we would predict the extreme cold adaptation of D. montana to have left a signature of genomic divergence associated with environmental and phenotypic differentiation across these loci. 2 | MATERIALS AND METHODS 2.1 | Sample collections and DNA extraction We collected samples of 49– 50 wildcaught flies from six D. montana populations from a range of latitudes from 66°N to 38°N in the spring of 2013 or 2014. Four of these populations represented a range of latitudes in North America (N.A.), and two populations were from a range in Finland (Figure 1; Table 1). Samples of wildcaught flies from the six populations were stored in ethanol (the male/female ratio varied across samples; Table 1) and DNA of individual flies was extracted using CTAB solution and phenolchloroformisoamylalcohol purifications in 2016. Genomic DNA was extracted from individual flies and quantified using Qubit (Thermo Fisher Scientific), and an equal amount of DNA from each individual (50 ng) was pooled into the final sample. Sequencing was performed at the Finnish Functional Genomics Centre in Turku, Finland (www.btk.fi/ funct ional - genomics) on the Illumina HiSeq3000 platform (pairedend reads, read length = 150 bp, estimated coverage ~121x). 2.2 | Phenotyping We measured the critical thermal minimum (CTmin) and chill coma recovery time (CCRT) of flies from six populations. Fly samples for these tests were collected for five populations (Seward, Terrace, Ashford, Crested Butte, and Korpilahti) from massbred populations that have been maintained in the laboratory since 2013– 2014. For the Oulanka population, flies were collected from three isofemale strains (established in 2014), because the massbred population had been contaminated by another species. The massbred populations were originally established from F4 progenies of 20 isofemale strains each (400 flies) and have been maintained in constant temperature (19 ± 1℃) and light (24 h of light, to prevent the flies from entering diapause) regimes for about 20– 25 generations before the experiment. All flies were supplied with fresh malt medium in halfpint bottles every week (Lakovaara, 1969). The newly emerged flies were collected using light CO2 anaesthesia within 24 h after emergence, separated by sex and placed in maltvials in the same conditions until sexual maturity (20– 21 days old) and were then used in CTmin and CCRT tests. The same individual flies were first assayed for CTmin and then for CCRT, the flies were not anaesthetised before these tests. We assayed a total of 328 females and 302 males for CTmin and CCRT. These assays were done in batches of between 22 and 30 flies, and split evenly by sex (21 batches in total). Between 32 and 46 (mean 39) flies per population per sex were tested (for the FIGURE 1 Maps of all D. montana populations in this study. Panels show population locations from (a) Finland and (b) North America showing the locations of all populations sampled. Labelled, blue circles give the locations of populations sampled for phenotyping and sequencing (a) (b) 3786 | WIBERG Et al. Oulanka isofemale strains, between 15 and 39 flies per strain per sex were used). CTmin tests are based on detecting the temperature (CTmin) at which flies lose neuromuscular function and enter reversible state of chill coma (Andersen et al., 2015). In these tests, the flies were placed into tubes sealed with parafilm and submerged into a 30% glycolwater mixture in Julabo F32HL chamber. The temperature was decreased at the rate of 0.5℃ per min (from 19℃ t o – 6 ℃) and CTmin was determined by eye, as the temperature at which a fly was unable to stand on its legs. Immediately after the CTmin test, the temperature was set to – 6℃ and the flies were left in this temperature for 16 h. Vials were then quickly taken out of the glycolwater bath and the flies’ CCRT was determined as the time required for the flies to recover from chill coma and stand on their legs. The ambient room temperatures was recorded during trials but in an initial analysis including source population and room temperature there was no statistically significant effect of room temperature on CCRT (F5,569 = 0.72, p = .4), this variable was therefore left out of further analyses. To investigate population differences in CTmin and CCRT phenotypes we fit general additive mixed models (GAMMs) in r (v. 3.6.3; R Development Core Team, 2020) using the “mgcv” r package (Wood, 2004). In simple linear models including population, sex, and experimental batch as a fixed effects, experimental batch had an effect on CTmin (F20,533 = 2.8, p < .01), although compared to population (F5,533 = 12, p < .01), and sex (F1,533 = 11.5, p < .01) this effect was not strong. Experimental batch had a stronger effect on CCRT (F20,533 = 2.83, p < .01), compared to population (F5,533 = 4.5, p < .01) and sex (F1,533 = 1.5, p < .22). We therefore included it as a random effect in GAMM analyses. The full GAMMs included altitude, latitude, and sex as fixed effects and experimental batch as a random effect. We used a cubic regression spline as the basis smoothing function for both altitude and latitude. The raw data for all phenotyping are given in Table S1. The full models are shown in the Supporting Information (Supplemnetary Material). 2.3 | Bioclimatic variables and population geography We obtained representative climate data from the WorldClim database (Hijmans et al., 2005) for each D. montana population sampled for the poolseq (see above), as well as for 18 additional populations of this species (Table S1; Tyukmaeva et al., 2020). We downloaded the climate data and extracted the values corresponding to population coordinates using the R package “raster” (version 2.58; Hijmans et al., 2016). In total this amounts to 55 bioclimatic variables for each population (Table S1). To reduce the number of variables in the data set a principle components analysis (PCA) was performed using the “PCA()” function from the R package “FactoMineR” (version 1.28; Lê et al., 2008). Principle components were kept for further analysis if their eigenvalues were >1. PCA scores for each population were ztransformed using the “scale()” function in base R. Additionally, CTmin and CCRT were summarised to a mean value for each population. In total, this gives four “environmental” variables measured for each population (PC1, PC2, CTmin, and CCRT). 2.4 | Mapping, SNP calling and genomic analysis Quality of raw reads was checked with fastqc (v. 0.11.5; Andrews, 2015) and reads were trimmed using trimmomatic (v. 0.32; Bolger et al., 2014; see Supporting Information (Supplementary Material) for full trimming parameters). Trimmed reads were mapped to the D. montana reference genome (Parker et al., 2018) using bwa mem (v. 0.7.7; Li, 2013) with the default options but keeping only alignments with a mapping quality of >20 following best practice guidelines for poolseq (Schlötterer et al., 2014). Duplicate alignments were removed with samtools rmdup (v 1.3.1; Li et al., 2009) and regions around indels were realigned using picard (version 1.118; Broad Institute), gatk (v. 3.22; McKenna et al., 2010) and samtools. Separate.bam files for each of the sequenced populations were finally merged using bamtools (v. 2.4.0; Barnett et al., 2011). Over 80% of reads were properly mapped and retained in all samples. The mean coverage for Seward samples was nearly twice that of the other samples (Figure S1 and Figure S2). To remove the potential for this difference to cause artefacts in downstream analyses, the.bam files for Seward were downsampled to contain 94.1 million reads (the average across the five remaining populations). Median empirical coverage was between ~62 and 88x (Table S2; Figure S3) and much less variable among the populations, allowing common maximum and minimum thresholds to be set based on the aggregate distribution. Allele counts for each population at each genomic position were obtained with samtools mpileup (version 1.3.1; Li et al., TABLE 1 The sources of genomic samples (coordinates and the name of the nearest town), altitude of the sampling site, the year in which sampling was performed, and the number of males and females sampled (M/F) for each pool Source Sampling site Year M/F USA, Alaska Seward 60°9’N; 149°27′W Altitude 35 m 2013 30/20 Canada, British Columbia Terrace 54°27′N; 128°34′W Altitude 217 m 2014 22/27 USA, Washington Ashford, 46°45′N; 121°57′W Altitude 573 m 2013 16/34 USA, Colorado Crested Butte 38°54′N; 106°57′W Altitude 2,900 m 2013 36/13 Finland Oulanka 66°22′N; 29°20′E Altitude 337 m 2013 25/25 Finland Korpilahti 62°20′N; 25°34′E Altitude 133 m 2013 27/23 | 3787 WIBERG Et al. 2009) using options to skip indel calling as well as ignoring reads with a mapping quality <20 and sites with a base quality <15. This was followed by the heuristic SNP calling software PoolSNP using a minimum count of 5 to call an allele, and a minimum coverage of 37 and a maximum coverage <95th percentile of the scaffoldwide coverage distribution to call a SNP (Kapun et al., 2020). Even if all these filters were passed, an allele was not considered if its frequency was <0.001. Finally, we only considered SNPs on scaffolds >10 kb in length. The final set consisted of 2,190,511 biallelic SNPs that could be placed on scaffolds ordered along the chromosomes and were used in downstream analyses. To test for an association between the four environmental variables and genetic differentiation we used bayescenv (v. 1.1; Foll & Gaggiotti, 2008; de Villemereuil & Gaggiotti, 2015). bayescenv fits a model of FST to population differentiation for each locus, incorporating environmental differentiation as a predictor (included as the parameter g) while fitting two locusspecific effects, one for environmental selection, and the other for other processes (demography, other types of selection). It therefore controls for confounding effects of population structure/relatedness in testing for an association with environmental variables. bayescenv was run with five pilot runs of 1000 iterations each, followed by a main chain of 4000 iterations of which 2000 were discarded as burnin. Four MCMC chains were run for each analysis to evaluate convergence of the chains to common parameter estimates. Because of the unbalanced number of males (and the resulting variation in ratios of X:Y chromosomes) in the pools, bayescenv analyses were performed separately on SNPs that could be assigned to the autosomal linkage groups (chromosomes) and the X chromosome. Raw count data were used for the autosomal data. For X linked SNPs, allele count data were scaled to the known number of X chromosomes in the pool using neff, the effective sample size taking into account the multiple rounds of binomial sampling inherent to a poolseq design (Feder et al., 2012; Kolaczkowski et al., 2011). Neff scales the allele counts at each SNP downwards based on the known number of chromosomes in the pool (see Table 1). Chains were assessed for convergence with the “coda” r package (v. 0.191; Plummer et al., 2006). Convergence was reached across the four chains for most analyses and parameters (potential scale reduction factors (PSRFs) of ~1 in a GelmanRubin diagnostic test; Figures S4– S7), except for analyses of autosomal SNPs and PC2 as the environmental variable which showed mild convergence problems (PSRF = 1.71), although parameter estimates agreed well with the other chains. Thus, this first chain (black line in Figures S4– S7) was discounted for all parameters and only estimates from the remaining three chains were used. The union of significant SNPs (using qvalues for the g parameter, describing the association between genetic differentiation at a locus and environmental differentiation, to control the FDR at 0.05, i.e., SNPs with qvalues <0.05) across these chains were taken as the final candidate SNPs. Finally, population genetic statistics (π, and Tajima's D) were computed in windows of 10 kb with a step size of 5 kb using methods implemented for poolseq data (Kapun et al., 2020). Windows contained a mean of 488.2 ± 1.2 SNPs (mean ± SE). These statistics were only computed for scaffolds with a length >50 kb. SNPwise FST was computed for each population with the r package “poolfstat” (v. 1.1.1; Hivert et al., 2018) by first computing all pairwise values, and then deriving population specific FST values by averaging across all pairwise values where a population was included. We computed 95% confidence intervals for mean SNPwise FST on each chromosome from the distribution of mean FST values across 100 bootstrap samples of SNPs. See the Supplementary Material (Supporting Information) for pseudocode commands of the key pipeline steps. We performed GO term enrichment analyses with david (v. 6.8; Huang et al., 2009a, 2009b) and gowinda (v. 1.12; Kofler & Schlötterer, 2012), which accounts for differences in gene length. For david, since the D. montana annotation contains information about orthologs in D. virilis, we extracted all genes within 10 kb upstream or downstream of a candidate SNP with an ortholog in D. virilis and submitted them to david. For analyses with gowinda, because the genesets need to be given manually, we obtained genesets from FuncAssociate3 (Berriz et al., 2009). As above, we considered SNPs within 10 kb of a gene (- - genedefinition updownstream 10,000) and performed 1 million simulations to obtain the empirical pvalue distribution. Regulatory elements, such as enhancers, and transcription factor binding sites, can occur up to 1 Mb upor downstream from a target gene in other species (Chan et al., 2010; Maston et al., 2006; Pennacchio et al., 2013; Werner et al., 2010) but generally lie within 2 kb of a gene region in D. melanogaster (Arnosti, 2003), thus 10 kb represents a compromise. 3 | RESULTS 3.1 | Cold tolerance measures, bioclimatic variables and population geography Across individuals, there was no evidence of an association between minimum critical temperature (CTmin; pooled across sexes) and the chill coma recovery time (CCRT; cor = 0.08, p = .07). Neither was there evidence for an association between these traits across populations (cor = 0.22 p = .68). CTmin was significantly different between sexes, with females being more cold tolerant (t = 3.14, p = .002). While latitude appears to have a nonlinear effect on CTmin (Figure 2a and Figure S8A), the effective degrees of freedom (edf = 1) of the partial effect of latitude suggested a largely linear effect after accounting for altitude (F = 40.4, p < .001). Altitude had a more complex effect (edf = 1.262, F = 24.0, p < .001), but with a strong altitude outlier in Crested Butte (Table 1), this result should be taken with some caution. Moreover, the overall adjusted Rsquared was low at only 0.09. For CCRT, only latitude had a significant effect (edf = 1, F = 17.2, p < .001; Figure 2b). Taken together, CTmin was lower for populations at higher latitudes as well as for populations at higher altitudes (Figure 2a and Figure S8A). As one would expect, CTmin showed on average lower values (Figure 2a and Figure S8A), and CCRTs were 3788 | WIBERG Et al. shorter (Figure 2b and Figure S8B), at higher latitudes, meaning that more northern and higher altitude populations show higher cold tolerance. We performed a principal component analysis (PCA) of the WorldClim climate data for a total of 24 D. montana populations, from which we had collected samples and where climate data were available. We wanted to examine whether the populations chosen for the cold tolerance tests and genome sequencing were representative of the range of environmental variation experienced by the species. These analyses identified four principal components (PCs) that together explained about 98% of the variation (Figure 3a and Figure S9). The first two PCs separated the populations roughly by a measure of distance from the coast (PC1) and then by latitude and altitude (PC2). PC1 explained ~54% of the variation (Figure 3a and b) and loaded heavily on climate and biological variables associated with precipitation and temperature such as “mean temperature of coldest FIGURE 2 Population phenotypes. (a) and (b) show the variation in CTmin and CCRT across populations and latitude, respectively. Solid and dashed lines show the predicted values from the partial effect of latitude from the best model (see Results) for males and females, respectively. In (b) although the points are plotted separately for males and females, the best model only included latitude as a covariate. In both panels the points and error bars give the mean and standard errors, respectively. Figure S8 shows the full data for each population (a) (b) FIGURE 3 Principle component analysis (PCA). (a) The distribution of all populations along the two first PC axes. The loadings of variables on each axis can be found in Table S3. Blue circles give the populations that have been poolsequenced for this study, red squares and triangles give the other Finnish and North American populations, respectively. (b) and (c) give PCs 1 and 2 as a function of latitude and altitude, respectively. The legend for altitude in both (b) and (c) is given in (c) (a) (b) (c) | 3789 WIBERG Et al. quarter”, “precipitation of wettest month”, “annual precipitation” (see Figure S9 and Table S3). Meanwhile, PC2 explained ~23% of the variation (Figure 3a and c) and loaded heavily on biological variables that are associated with latitudinal clinality, e.g., “mean diurnal [temperature] range,” and “isothermality” which is the diurnal range divided by the mean “annual [temperature] range”. The remaining PCs (PC3 and PC4) explained about 11.5% and 5% of the variation, respectively and did not capture as much of the climatic variation, and were therefore not considered further. Latitude (Spearman's rank correlation: rho = – 0.50, p = .01) but not altitude (rho = – 0.19, p = .39) correlated significantly with PC1. However, both altitude (rho = – 0.50, p = 0.02) and latitude (rho = – 0.59, p = .003) correlated with PC2 (see also Figure 3c). Importantly, these patterns were fairly robust also when performing PCA using only the six populations for which genomic data were collected. Loadings on PC1 were highly comparable, and although correlations with latitude and altitude did not achieve statistical significance, they were similar in magnitude and direction (latitude vs. PC1: rho = 0.43, altitude vs. PC1: rho = – 0.14, p = .8). For PC2 the correlations were also not statistically significant and differed for latitude both in magnitude and direction (rho = 0.6, p = .2), while for altitude, the correlation remained similar and marginally nonsignificant (rho = – 0.82, p = .06). Moreover, while the range of the sequenced populations covered ~70% of the range of all populations for PC1, they covered only ~50% of the range for PC2 (Figure 3a). This suggests that environmental variation across the six populations selected for sequencing reflects that experienced by all 24 populations at least for the PC1 axis. Therefore, the relationship between environmental variables and genetic differentiation in the samples selected for poolseq is likely to reflect true patterns across populations of D. montana, though power might be somewhat reduced for PC2. For all subsequent analyses, we used the results of the PCA using all 24 populations, and to examine the association between climate and phenotype, we compared these across the six populations. CTmin was positively correlated with PC1 (Pearson's correlation coefficient (cor) = 0.94, p < .01) and had a marginally nonsignificant association with PC2 (cor = 0.75, p = .09). However, CCRT showed no relationship with either PC1 (cor = 0.09, p = .87) or PC2 (cor = – 0.46, p = .36) although the small sample sizes (N = 6 in all cases) make reliable conclusions difficult (see Figure S10). 3.2 | Genomic analyses The number of SNPs with a significant (qvalue <0.05) association between populationbased FST, the two cold tolerance measures (CCRT and CTmin), and the two PCs of the bioclimatic data varied from 312 (chromosome 3, CCRT) to 2480 (chromosome 4, PC2) across the chromosomes (Figure 4). Using PC1 as an environmental variable with bayescenv gave a total of 2976 and 1528 SNPs with a qvalue <0.05 on the autosomes and on the X chromosome, respectively. Interestingly, the distribution across the chromosomes was not random. Using the distribution of all SNPs to obtain expected counts, there was a significant deviation from expectation (Χ2 = 2906.4, d.f. = 4, p < .001). There were many more SNPs than expected on chromosome 4 (1432 vs. 954) and on the X chromosome (1528 vs. 526; Figure 4). Results were similar for PC2 with 6607 and 1861 SNPs with a qvalue <0.05 on the autosomes and X chromosomes, respectively. Again, there was a significant deviation from the expected distribution of SNPs across the chromosomes (Χ2 = 1681.9, d.f. = 4, p < .001) with an overrepresentation on the fourth (2480 vs. 1794) and the X chromosomes (1861 vs. 989; Figure 4). We also used average CTmin values per population in similar analyses and found a total of 2668 and 1272 SNPs with a qvalue <0.05 on the autosomes and X chromosomes, respectively. The pattern of significant deviations from expected distributions (Χ2 = 2,526.2, d.f. = 4, p < .001) was also due to an excess on the fourth (1383 vs. 835) and the X chromosomes (1272 vs. 460; Figure 4). Similar results were found for CCRT with a total of 2240 and 1228 SNPs with a qvalue <0.05, respectively. Once again, there was a significant deviation from the expected distribution of SNPs (Χ2 = 2825.6, d.f. = 3, p <.001) with and excess on the fourth (1252 vs. 735) and the X chromosomes (1228 vs. 405; Figure 4). Manhattan plots of the distribution of SNPs across chromosomes are given in Figures S11– S14. To more closely examine the loci implicated in the four BayScEnv analyses, we identified genes within 10 kb of, or containing the candidate SNPs (Table S4). Overall, there is quite a large overlap among these genes with ~39% (1102 in total) of them being shared by all the four analyses (Table S4, Figure S15). This is far more than would be expected (mean ± SE: 3.6 ± 0.06) by randomly resampling (without replacement) a number of genes equivalent to that associated with top SNPs for each BayScEnv analysis from the D. montana annotation, then computing the fourway intersection. This large observed overlap therefore probably reflects a genuine overlap among the underlying top SNPs (data not shown). Some (147, ~13%) of these common genes are novel to D. montana (i.e., not annotated in D. virilis or in other Drosophila spp.) and therefore have no annotation, but 955 have an identifiable D. virilis ortholog (Table S4). The second largest set are those genes unique to the analysis of PC2 as an environmental covariate (Figure S15). We performed functional enrichment analyses with david keeping all annotation clusters with an enrichment score >1.3 (corresponding to an average corrected pvalue of .05). This revealed several common categories of genes associated with the climatic variables and population phenotypes (Table S5). For example, terms associated with membrane and transmembrane structures, immunoglobulins, HAD hydrolase and nucleotide binding were enriched in most of the variables (Table S5). Interestingly, there were also several gene ontology categories that were only enriched in one of the variables, such as glycoside and ATPase hydrolase in CCRT and ion channels and transport, as well as metal binding in PC1 (Table S5). gowinda analyses revealed significant enrichment of GO terms (after accounting for variation in gene lengths and correcting for multiple testing) only for genes near SNPs associated with variation in CCRT across populations. Interestingly, the term “carbohydrate derivative binding” was identified among the most enriched terms (Table S6), which agrees closely with some terms identified using david for the same variable (Table 3790 | WIBERG Et al. S5). Similarly, for PC2, “nucleotide binding” was among the top scoring enriched terms in both gowinda and david analyses, though this was not significant after controlling for multiple testing in gowinda (Table S5, Table S6). We then compared loci near candidate SNPs with genes implicated in previous studies of climatic adaptation in D. montana including gene expression studies of traits connected to diapause and coldtolerance (Kankare et al., 2010, 2016; Parker et al., 2015, 2016). Additionally, several candidate genes have been identified near the most significantly differentiated SNPs among D. montana populations from Oulanka (Finland), and from North American populations in Colorado and Vancouver (Parker et al., 2018). Finally, quasinatural selection experiments identified several genes within 10 kb of SNPs responding to selection for a shorter CDL for diapause induction (Kauranen et al., 2019). We tested for an overlap between the total set of genes within 10kb of outlier SNPs from all of the BayScEnv runs (N = 2694) and the candidate gene sets identified in earlier studies (see Table S7 for the gene sets and studies used). We computed a bootstrap distribution of overlaps by sampling 2694 random genes from the D. montana annotation (N = 13,683). For each of the gene sets from previous studies this was done 100 times and the distribution compared to the empirical overlap. Results are given in Table 2. In all the cases the empirical overlap was greater than expected by chance with empirical pvalues <.05 (ranging from <.0001 to .01; Table 2). The only gene that was found in all five of the previous studies used and in the comparison here is called sidestep II (sideII; Table S8). Unfortunately, there is no information available about the biological processes or molecular functions connected to it. Moreover, from 44 other genes that were common to four of our previous studies and this study (Table S8) most (27) have an ortholog in D. melanogaster. These genes have molecular functions such as transmembrane signaling or transporter, acetylcholinesterase, ATP binding, protein serine/threonine kinase, carboxylic ester hydrolase, or Rho guanylnucleotide exchange factor activity (Thurmond et al., 2019). Many of the genes are also connected to metal ion, nucleid acid or zinc ion binding (Table S8). After identifying information on molecular or biological function and Interpro domains, eventually only five genes remained for which there was no information available (Table S8). Examination of population genetic parameters identified the Crested Butte population as anomalous. The distribution of Tajima's D is centred close to zero in most populations, being slightly more negative in North American populations (Figure S16A). However, Crested Butte is an outlier with a greatly reduced genomewide Tajima's D (Figure S16A). Similarly, diversity (π) is also lower in this population than in other populations. There is no overall relationship between latitude and π (Spearman's rho = 0.14, N = 6, p = .8; Figure S16B) but there is a strong correlation between latitude and Tajima's D which is influenced by this population (with Crested Butte: rho = 0.88, N = 6, p = .03; without Crested Butte: rho = 0.8, N = 5, p = .13). Although Crested Butte occurs at a much higher altitude (>2800 m) than other populations neither Tajima's D nor π correlated significantly with altitude (Tajima's D: rho = – 0.6, N = 6, p = .24, π: rho = – 0.6, N = 6, p = .24). Furthermore, FST was similar across all populations and chromosomes with the exception of Crested Butte which remained an outlier with high FST (Figure 5). Bootstrapped 95% confidence intervals are presented to give a guide to statistical significance of differences among populations and confirm that Crested Butte has high FST (Figure 5). Finally, FST was always highest FIGURE 4 The observed counts of candidate SNPs (i.e., SNPs with a qvalue <0.05) across chromosomes for each environmental variable (see Section 2). The total number of SNPs, and the proportion of all SNPs, on each chromosome are given below each set of bars. The expected counts on each chromosome, obtained from the proportions of all SNPs across chromosomes, are shown as points aligned with each bar Reference N genes N overlaps Mean (95% CI) resampled N overlaps Empirical pvalue Kankare et al., (2010) 14 82.77 (0, 6) .0004 Parker et al., (2015) 42 13 6.28 (2, 11) .001 Kankare et al., (2016) 3,929 946 773.14 (732, 814) <.0001 Parker et al., (2016) 130 27 18.72 (11, 26) .01 Parker et al., (2018) 2,114 629 414.7 (382, 448) <.0001 Kauranen et al., (2019) 1,751 402 344.64 (314, 375) <.0001 Shown are the citations for the original study, the number of genes identified from the original study, the number of overlapping genes with the candidate gene set from this study, and the mean and 95% CI from 10,000 resampled sets of candidate genes. TABLE 2 The overlap of the union of genes within 10 kb of SNPs associated with PC1, PC2, CTmin, and CCRT and previous candidate gene sets (see Table S7)