The challenge of separating signatures of local adaptation from those of isolation by distance and colonization history: The case of two white pines

Abstract Accurately detecting signatures of local adaptation using genetic‐environment associations (GEAs) requires controlling for neutral patterns of population structure to reduce the risk of false positives. However, a high degree of collinearity between climatic gradients and neutral population structure can greatly reduce power, and the performance of GEA methods in such case is rarely evaluated in empirical studies. In this study, we attempted to disentangle the effects of local adaptation and isolation by environment (IBE) from those of isolation by distance (IBD) and isolation by colonization from glacial refugia (IBC) using range‐wide samples in two white pine species. For this, SNPs from 168 genes, including 52 candidate genes for growth and phenology, were genotyped in 133 and 61 populations of Pinus strobus and P. monticola, respectively. For P. strobus and using all 153 SNPs, climate (IBE) did not significantly explained among‐population variation when controlling for IBD and IBC in redundancy analyses (RDAs). However, 26 SNPs were significantly associated with climate in single‐locus GEA analyses (Bayenv2 and LFMM), suggesting that local adaptation took place in the presence of high gene flow. For P. monticola, we found no evidence of IBE using RDAs and weaker signatures of local adaptation using GEA and F ST outlier tests, consistent with adaptation via phenotypic plasticity. In both species, the majority of the explained among‐population variation (69 to 96%) could not be partitioned between the effects of IBE, IBD, and IBC. GEA methods can account differently for this confounded variation, and this could explain the small overlap of SNPs detected between Bayenv2 and LFMM. Our study illustrates the inherent difficulty of taking into account neutral structure in natural populations and the importance of sampling designs that maximize climatic variation, while minimizing collinearity between climatic gradients and neutral structure.


| INTRODUCTION
Climate is a major factor affecting the distribution of genetic diversity among natural populations of plants. Tree species generally exhibit moderate to high among-population genetic variation for adaptive traits along climatic gradients (Alberto et al., 2013;Savolainen, Pyhäjärvi, & Knürr, 2007). Despite such evidence of local adaptation from common-garden studies, patterns of population structure observed at nuclear loci are often considered to result from neutral processes affecting the whole genome, including genetic drift, gene flow, and past demographic events (e.g., recent range contractions and expansions). A more recent view is that natural selection can also affect genomewide population divergence if gene flow among ecologically divergent habitats is reduced because of selection acting against nonlocally adapted migrants (Hendry, 2004;Nosil, Vines, & Funk, 2005), or because of other nonadaptive processes (Wang & Bradburd, 2014).
Disentangling IBE from neutral patterns of genetic variation is challenging (Shafer & Wolf, 2013;Wang & Bradburd, 2014). For example, decreasing gene flow with increasing geographic distance due to restricted dispersal (i.e., isolation by distance, IBD; Wright, 1943) can produce patterns similar to IBE when geography is correlated with environmental variation (Meirmans, 2012;Orsini, Vanoverbeke, Swillen, Mergeay, & De Meester, 2013). Postglacial recolonization can also generate allele frequency gradients similar to IBE or IBD as a result of repeated founder events and "allele surfing" along the colonization front (de Lafontaine, Ducousso, Lefèvre, Magnanou, & Petit, 2013) because colonization routes often covary with environmental gradients. Furthermore, postglacial recolonization from different glacial refugia followed by secondary contact can also create genetic barriers (hereafter referred to as isolation by colonization, IBC) that often coincide with environmental clines (e.g., Bierne, Welch, Loire, Bonhomme, & David, 2011;Richardson, Rehfeldt, & Kim, 2009). Hence, because the selective climatic gradients, geography, and postglacial recolonization routes are often spatially correlated in natural populations, it is extremely difficult to separate the relative effects of IBE from those of IBD and IBC. However, disentangling these effects is important to accurately control for neutral population structure (e.g., IBD and IBC) when looking for signatures of local adaptation.
The sampling design also impacts the ability to detect IBE (Wang & Bradburd, 2014), F ST outliers, and GEAs (Lotterhos & Whitlock, 2015;Meirmans, 2015). For GEA and IBE analyses, power can be improved by sampling individuals from as many climatically variable populations as possible across the range of a species, at the cost of sampling fewer individuals per population (De Mita et al., 2013;Poncet et al., 2010;Wang & Bradburd, 2014). In addition, simulations showed that increasing the total number of sampled individuals increased the power of GEA and F ST outlier analyses (Lotterhos & Whitlock, 2015). Statistical methods that take into account uncertainty due to small population sample sizes (e.g., Coop et al., 2010;Foll & Gaggiotti, 2008;Frichot et al., 2013) are well suited for sampling schemes that aim to maximize environmental variation by including a large number of populations in order to improve our ability to detect signatures of local adaptation in natural populations.
Another promising avenue to detect adaptive loci of importance is to compare signatures of adaptation among closely related species or evolutionary lineages using a set of orthologous genes (i.e., genes that descended from a common ancestral gene by speciation; e.g., Grivet et al., 2011;Mosca et al., 2012;Chen et al., 2014;Zhou, Zhang, Liu, Wu, & Savolainen, 2014). Evidence of convergent evolution (Arendt & Reznick, 2008), or the repeated evolution of similar phenotypes from similar genetic mechanisms is increasing (Stern, 2013;Jones et al., 2012;Yeaman et al., 2016), but it is currently limited to a few taxa. Eastern white pine (Pinus strobus, Figure 1) and western white pine (P. monticola) diverged <12 million years ago (Gernandt et al., 2008) and are well suited for studying local adaptation as both species are distributed latitudinally and longitudinally across a wide variety of climates in North America. However, these two species have demographic histories that could complicate the detection of signatures of local adaptation and IBE. Populations of both species cluster into southern and northern genetic groups, likely resulting from range expansion from multiple glacial refugia (Rehfeldt, Hoff, & Steinhoff, 1984;Nadeau et al., 2015; but see Richardson et al., 2009;; who suggested a single refugium).
Differentiation between the phylogeographic groups may also be in part due to adaptation to contrasting climates as the northern and southern groups differ in their adaptive traits (e.g., height growth potential, cold hardiness; Rehfeldt et al., 1984;Richardson et al., 2009;Joyce & Rehfeldt, 2013).
Here, we look for evidence of signatures of local adaptation and IBE using single nucleotide polymorphism (SNP) markers developed from 168 orthologous genes and genotyped on 133 P. strobus and 61 P. monticola populations distributed across their natural ranges. We addressed the following questions: (1) "Can we detect genes showing signatures of local adaptation to climate in each species and in both species?" and (2) "Did local adaptation to climate contribute to the observed population structure (IBE) in P. monticola and P. strobus, or was it mostly driven by neutral processes (i.e., IBD or IBC)?"

| Sampling and SNP dataset
To investigate patterns of adaptation, we used a previously developed dataset (Nadeau et al., 2015), in which 153 (120 genes) and 158 SNPs (127 genes) were genotyped on 831 individuals (133 populations) of P. strobus and 348 individuals (61 populations) of P. monticola ( Figure 2). A selection from samples available in provenance trials and seedbanks (see Nadeau et al., 2015 for details) was made to cover a large range of climatic conditions across the natural distribution of each species. To do so, we performed a principal component analysis (PCA) on seven climatic variables (see "2.2 Climatic data"; Table 1) obtained for all available samples of each species, using the prcomp function in R (R Core Team 2015), and we selected populations that covered a wide range of variation in the first two principal components ( Figure 3). Note that for P. monticola, many provenances from southern Oregon and California were not available because they had died in the Whidbey Island provenance trial (WA, USA) before sampling.
SNP development was conducted in parallel using putative orthologous gene sequences available for both P. strobus and P. monticola (i.e., sequences amplified using the same primers in both species ;Nadeau et al., 2015). Briefly, an initial set of 118 gene sequences from the White Pine Resequencing Project (WHISP, http://dendrome.ucdavis.edu/ wpgp; Eckert et al., 2013), randomly distributed across the genome, was selected. We also included 23 candidate genes for growth, phenology, and cold hardiness in Picea glauca, 24 candidate genes for wood formation in P. glauca, one candidate gene for adaptation to aridity in Pinus taeda, and two gene sequences available from GenBank (see Nadeau  Rigault et al., 2011), yielding a total of 52 candidate genes putatively involved in growth, phenology, and cold hardiness (El Kayal et al., 2011;Holliday, Ralph, White, Bohlmann, & Aitken, 2008;Pelgas, Bousquet, Meirmans, Ritland, & Isabel, 2011). To provide complementary information for the white pine sequences that did not had a significant blast hit on the RefSeq database, we obtained the P. glauca best-ortholog annotations (GCAT 3.3 sequences are complete or near complete) from the Arabidopsis database (TAIR, https://www.arabidopsis.org/index.jsp). This was particularly useful for partial white pine sequences that were mainly composed of intron sequences.
Of 168 orthologous genes, 79 contained SNPs in both species, including 34 orthologous SNPs (i.e., occurring at the same nucleotide position in both species). Sixty-eight Pinus strobus SNPs and 72 P. monticola SNPs occurred at different nucleotide positions within orthologous genes. Forty-one genes (51 SNPs) and 48 genes (52 SNPs) contained SNPs only in P. strobus and P. monticola, respectively. We deduced SNP annotations (i.e., noncoding, synonymous, nonsynonymous) for 71 fully annotated genes from the WHISP dataset . For the other gene sets, the Picea glauca gene catalog was used to deduce coding regions and SNP annotations.

| Climatic data
Climate normals for each population for the 1961-1990 period were obtained using ClimateNA (Wang, Hamann, Spittlehouse, & Carroll, 2016). We selected seven climatic variables that did not covary strongly (r < .80) in at least one of the species (Table 1). In other words, a climatic variable that was highly correlated (r > .80) in one species could still be retained if it was less correlated in the other species (r < .80) to ensure that we did not miss any important climatic variation ( Figure S1, Appendix S2). In addition, to reduce collinearity with the main axes of ancestry, we ensured that the selected variables were not highly correlated (r < .70) with the Q-values from STRUCTURE (K = 2 within each species) obtained from Nadeau et al. (2015). We also included elevation as a climatic surrogate (eighth climatic variable) as it represents many climatic gradients upon which selection  Nadeau et al., 2015). Variation along PC1 (x-axis) and PC2 (y-axis) was used to select samples for genotyping in order to cover a wide range of environmental variation. Genotyped populations are colored according to their genetic group membership as in Figure 2. Available populations that were not genotyped (gray dots) were either not sampled or failed genotyping. Ellipses represent the 95% confidence intervals for each group. Insets show the proportion of variation explained by each PC can act and should not be strongly correlated with patterns of IBD (e.g., two populations on different mountain tops may have similar climates, but each of them is spatially closer to their warmer, lower-lying mountain flanks than to the other mountain top). Reduction in climatic variables to principal components was avoided to make possible direct and easily interpretable comparisons between our study species.

| F ST outlier tests
All analyses were performed separately for each species. We first looked for signatures of selection using F ST (Table S1, Appendix S1). We chose to report results with PO = 1,000 because increasing PO has been shown to reduce the number of false positives without greatly affecting the ability to detect true positives (Lotterhos & Whitlock, 2014). The internal q-value function provided in BayeScan was used to assess significance, and outliers were reported at FDR <5% (q < 0.05).

| Genetic-environment associations
Signatures of local adaptation to climate were investigated using two GEA methods that take into account neutral population structure: Bayenv2 Gunther & Coop, 2013) and LFMM (Frichot et al., 2013). We first ran Bayenv2 using the entire SNP dataset and 100,000 Markov Chain Monte Carlo (MCMC) runs to estimate the covariance matrix ( Figure S2, Appendix S2). We then tested associations between each SNP and each of the eight climatic variables, while including the covariance matrix as a null model, by running Bayenv2 in "test mode" with 100,000 MCMC runs. Bayes factors (BF) were averaged across 10 replicates using 10 independent estimates of the covariance matrix. The average correlation among replicates of the covariance matrix (P. strobus: r = .694; P. monticola: r = .794) and of BFs (P. strobus: r = .863; P. monticola: r = .716) were moderately high. The significance of each tested locus was determined according to Jeffrey's scale of evidence (Jeffrey, 1961): BF > 3, BF > 10, BF > 32, and BF > 100 indicated substantial, strong, very strong, and decisive support for selection, respectively.
The second GEA method used latent factor mixed models (LFMM), as implemented in the software LFMM v.1.4 (Frichot et al., 2013). This method uses a hierarchical Bayesian mixed model based on a variant of PCA, in which neutral population structure is introduced via (k) unobserved or latent factors. We implemented the LFMM method using the default individual-based data specification to avoid potential biases due to unequal population sample sizes (de Villemereuil et al., 2014).
To determine k, we performed a PCA on individual allele frequencies using the LEA package in R . For each species, a Tracy-Widom test indicated that seven principal components significantly explained genetic variation (Table S2, Appendix S1), so we ran LFMM using k = 7 for each species. For each test, 10,000 iterations of the Gibbs sampling algorithm were run, with the first 5,000 iterations discarded as burn-in. Z-scores from 10 independent replicate runs were combined using the Fisher-Stouffer method, and the resulting p-values were adjusted using the genomic inflation factor (λ). For k = 7, the average λ was close to or lower than 1 in each species (P. monticola: λ = 0.92; P. strobus: λ = 1.68) as recommended. A Benjamini and Hochberg (1995) FDR correction of 5% was applied to p-values using the qvalue package in R (Storey, 2002). Figure S3 (Supporting information) shows the effect of the choice of k on the number of SNPs associated with each climatic variable. The overlap in outlier SNPs among analyses using different values of k was generally high, and the smaller overlap in P. monticola was mostly due to a decrease in the number of outlier SNPs with increasing k ( Figure S4, Appendix S2).

| IBE, IBD, and IBC
To test for IBD, we estimated the correlation between a matrix of pairwise Slatkin's linearized F ST (F ST /(1-F ST )) and the matrix of log-transformed geographic distances (Rousset, 1997) calculated using the Geographic Distance Matrix Generator online tool (http:// biodiversityinformatics.amnh.org/open_source/gdmg). We then tested for IBE, while controlling for IBD, using partial Mantel tests.
Climatic distances for the eight climatic variables were computed as the Euclidean distance between pairs of populations using the dist function in R. The correlation between linearized F ST and climatic distance was tested for each climatic variable separately, with geographic distance included as a covariate. The significance of Mantel's r statistics for IBD and IBE was tested using n = 1,000 random permutations using the mantel function in the R ecodist package (Goslee & Urban, 2007). To control for multiple testing, p-values were converted into q-values, and a FDR of 5% was applied based on the Benjamini and Hochberg (1995) criteria using the qvalue package in R.
A series of redundancy analyses (RDAs) were performed to partition the among-population genetic variation into three components: IBE, IBD, and IBC. RDA is a multiple linear regression method performed between a matrix of dependent variables and matrices of independent variables. This type of multivariate analysis is more appropriate than Mantel tests when multiple climatic variables are analyzed to identify ecological drivers of population genetic structure (Orsini et al., 2013).
The dependent matrix contained allele frequencies for each population.
We included three independent matrices: (1) the eight climatic variables (representing IBE); (2) geographic variables (IBD); and (3) a north-south ancestry variable (IBC). For the geography matrix, we used a trend surface analysis (Borcard, Legendre, & Drapeau, 1992) to calculate secondorder polynomials and combinations of the coordinates of sampling locations (x, y, xy, x 2 , y 2 ) to ensure that linear gradients in the data, as well as more complex patterns, were extracted. To prevent overfitting, we used a forward selection procedure with a stringent alpha value of 0.01 (Lee & Mitchell-Olds, 2011). This resulted in the retention of four geographic variables for P. monticola (x, y, xy, y 2 ) and three for P. strobus (x, y, xy). Results were similar if we included only x and y (Table S3, Appendix S1). The north-south ancestry variable was the population Q-values from STRUCTURE, which separated populations into K = 2 glacial lineages (northern and southern) within each species (Nadeau et al., 2015). All three independent matrices were scaled to a mean of zero and a variance of one prior to analyses, but the dependent matrix was left untransformed. Among-population variation in each species was partitioned into exclusive effects of climate, geography, and northsouth ancestry (i.e., constrained by the effects of the remaining two independent matrices), as well as all possible combinations of these three matrices, using the varpart and rda functions of the vegan package in R (Oksanen et al., 2013). Significance of each partition was tested with the anova.cca function of vegan with a step size of 1000, resulting in at least 999 permutations.
We first performed RDAs using all SNPs to determine the main drivers of genomewide population structure. To see whether among-population differentiation at loci under divergent selection could be explained by climate, we performed two additional sets of RDAs using subsets of SNPs  Figure 2).

| F ST outlier and GEA tests
We first looked for signatures of selection using F ST Table 3). Two of these SNPs in P. strobus (M-015, M-016) were located within the same gene and were in moderate linkage disequilibrium (r = .402; Nadeau et al., 2015).
Finally, we looked for orthologous SNPs or genes that were detected as outliers in both species by any of the three methods (BayeScan, Bayenv2, LFMM). Three of the 79 orthologous genes contained outlier SNPs in both species (Table 4). Given the number of genes that contained outlier SNPs in each species separately (P. strobus: 25 genes; P. monticola: 18 genes), the number of genes containing outlier SNPs in both species did not differ from random expectation (p (≥3) = .415; permutation test with 10,000 random draws in R). None of the 34 orthologous SNPs (i.e., occurring at the same nucleotide position) were outliers in both species.

| Mantel tests and RDAs using all SNPs
We investigated the importance of IBE, IBD, and IBC as drivers of genomewide population structure with partial Mantel tests and RDAs.
Mantel tests detected significant IBD in both species as genetic distance increased with geographic distance (Table 5). In P. strobus, partial Mantel tests found that climatic distances for TD, bFFP, and eFFP were significantly correlated with genetic distances when controlled for geographic distance (p < .05), but only TD remained significant after correction for multiple testing (q = 0.036). In P. monticola, only elevational distance significantly explained genetic distance (p < .017), and this effect was marginally significant after correction for multiple testing (q = 0.077).
Using RDAs, we partitioned among-population genetic differentiation into three components: climate (IBE), geography (IBD), and north-south ancestry (Q-values from STRUCTURE) representing recolonization history from northern and southern glacial refugia (IBC). In the uncorrected RDAs, climate, geography, and north-south ancestry each explained significant proportions of the genetic variation in both P. strobus and P. monticola, as measured by the adjusted R 2 ("combined fractions" in Table 6). A series of partial RDAs were performed to decompose their contribution to among-population variation ("individual fractions" in Table 6; displayed as Venn diagrams in Figure 5). A total of 8.4% and 17.6% of the variation in P. strobus and P. monticola, respectively, could be explained by the three components and their various combinations ("Total explained" in Table 6). In P. strobus, north-south ancestry (1.8%, p < .001, constrained by climate and geography) and geography (0.7%, p = .023, constrained by climate and north-south ancestry) explained significant proportions of variation, but climate did not (0.1%, p = .382, constrained by north-south ancestry and geography). Similarly, in P. monticola, significant variation could be attributed exclusively to north-south ancestry (2.1%, p < .001) and to geography (2.5%, p = .006), but not to climate (0%, p = .722). For both species, 69.0 to 73.9% of the explained variation was confounded between the effects of climate, north-south ancestry, and geography ("Total confounded" in Table 6). Finally, a large portion of the variation remained unexplained (P. strobus: 91.9%; P. monticola: 82.4%). This unexplained variation could be due to environmental variables that we did not take into account (e.g., soil composition or biotic interactions), but most of T A B L E 3 Highly supported candidate SNPs, that is, detected by a minimum of two different methods. Variable names in the "Bayenv2" and "LFMM" columns refer to climatic variables that were significantly associated with the SNPs

| RDAs using subsets of SNPs detected by Bayenv2 and LFMM
Finally, we performed RDAs using the subsets of candidate SNPs that showed signatures of local adaptation in Bayenv2 and LFMM analyses.
In P. strobus, a significant proportion of the among-population variation could be attributed exclusively to climate (LFMM: 2.5%, p = .010; marginally significant for Bayenv2: 1.6%, p = .091), but not in P. monticola (Table 6). A greater total proportion of the variation could be explained by climate, geography, north-south ancestry, and their various combinations when using SNPs detected by Bayenv or LFMM ("Total explained" in

T A B L E 5 Mantel and partial Mantel tests in Pinus strobus and P. monticola.
Correlation coefficients (r) between (1) genetic distance (Y) and geographic distance (D); and (2) between genetic distance (Y) and each of the eight climatic variables after controlling for D T A B L E 6 Redundancy analyses (RDAs) to partition among-population genetic variation (F) in Pinus strobus and P. monticola into three components: climate (IBE); geography (IBD); and north-south ancestry (IBC). F~clim.+geog. | anc.   isolation by distance (IBD) or recolonization history from glacial refugia (isolation by colonization, IBC), in shaping among-population genetic differentiation across the distribution of P. monticola and P. strobus. Using three GEA and F ST outlier methods, we detected signatures of local adaptation in P. strobus, but such signatures were weaker in P. monticola.
We found that, for the most part, the among-population genetic differentiation could not be partitioned into exclusive effects of IBE, IBD, and IBC in both species, thus making it difficult to separate signatures of local adaptation from neutral patterns of population structure.

| Signatures of local adaptation and IBE in Pinus strobus and P. monticola
Patterns of IBE occur when selection against nonlocally adapted migrants increase the genetic divergence among populations from different environments. IBE can be detected at neutral loci when divergence at selected loci extends to surrounding loci by "divergence hitchhiking" and eventually to the entire genome via "genome hitchhiking" (Feder, Egan, & Nosil, 2012;Feder & Nosil, 2010). In P. strobus, when using all SNPs, we did not detect strong evidence of IBE using partial Mantel tests (only TD-continentality-was significant; q < 0.05), and climate did not explain among-population variation in a RDA that controlled for IBD and IBC. However, single-locus GEA analyses found a relatively large number of SNPs associated with climate and a significant proportion of the variation at these SNPs could be exclusively attributed to climate in multilocus RDAs. This corresponds to a scenario where gene flow is reduced among ecologically distant populations at loci directly involved in local adaptation or at closely linked loci, while there are no selective constraints on gene flow among environments for the rest of the genome (Barton, 2000;Gavrilets & Vose, 2005;Wu, 2001). This result is not surprising considering the rapid decay of linkage disequilibrium in large outcrossing populations of conifers (Namroud, Guillet-Claude, Mackay, Isabel, & Bousquet, 2010) and the high levels of gene flow across the range of P. strobus (Mehes, Nkongolo, & Michael, 2009;Nadeau et al., 2015), which should uniformize among-population genetic variation at neutral loci. Provenance trial studies have previously found moderate among-population genetic variation for adaptive traits in P. strobus (e.g., height growth, bud phenology, cold hardiness; Li, Beaulieu, Daoust, & Plourde, 1997;Joyce & Sinclair, 2002;Lu, Joyce, & Sinclair, 2003a,b). Interestingly, the climatic variable "degree-days above 5°C" was involved in the greatest number of GEAs using both Bayenv2 and LFMM, and it was also the best climatic predictor of range-wide genetic variation in growth potential and phenology ( 2.5% ** 2013). Thus, some of the SNPs we detected in GEAs may be important for growth potential or phenology, but confirmatory evidence would be needed from common-garden or functional studies.
In a similar study on P. strobus, Rajora, Eckert, and Zinck (2016) did not detect many signatures of local adaptation in single-locus analyses, but detected significant associations with climate using multilocus analyses with a set of 44 candidate SNPs (25 genes). In their discussion, the authors suggested that local adaptation to climate was occurring via covariance in allele frequencies among loci of small effects, rather than via allele frequency changes at a few loci of larger effects (Latta, 1998(Latta, , 2004.  Lotterhos & Whitlock, 2015;Yeaman, 2015). Given the small proportion of the genome surveyed here, it seems unlikely that we captured a significant amount of adaptive covariance among loci, and so we abstain from drawing conclusions about the genetic architecture of local adaptation to climate in P. strobus.
For P. monticola, F ST outlier and GEA tests detected a smaller number of SNPs (of generally lower significance) than in P. strobus.
Moreover, climate did not explain among-population variation in RDAs after controlling for IBD and IBC, even for the subsets of SNPs that were detected by GEA methods. Previous studies showed no or little differentiation in phenotypic traits among populations within the large northern group, and it has been suggested that P. monticola has adapted to a wide variety of climates mostly via phenotypic plasticity (Chuine, Rehfeldt, & Aitken, 2006;Rehfeldt et al., 1984). For example, shoot elongation in this species is initiated later than in most temperate conifers due to a high threshold for forcing temperatures (average 10.2°C), with little genetic variation among populations (Chuine et al., 2006). Delayed shoot elongation would allow P. monticola to avoid late spring frost damage and to survive in a wide range of environments without the need to be locally adapted. Genetic differences for height growth potential and cold hardiness exist between the northern and southern group (Rehfeldt et al., 1984;Richardson et al., 2009).
However, the small sample size for the southern group and the severe corrections for population structure applied by Bayenv2 and LFMM ( Figure S5, Appendix S2) may have prevented us from separating signatures of selection from the neutral structure.

| The role of IBE, IBD, and IBC in shaping population structure
We attempted to determine whether the genomewide population structure of both species (i.e., using all SNPs) was determined by local adaptation to climate (IBE), geography (IBD), or postglacial recolonization from glacial refugia (IBC). For both species, IBD and IBC were significant drivers of population structure, but climate alone was not.
After controlling for IBE and IBD, north-south ancestry (Q-values from STRUCTURE) explained the largest amount of among-population variation in P. strobus and the second largest in P. monticola. This was expected since STRUCTURE looks for the dominant population structure patterns. For both species, populations from the northern and southern genetic groups detected by STRUCTURE may have originated from different glacial refugia, thus representing IBC (Nadeau et al., 2015), although others have suggested a single refugium (Richardson et al., 2009;. A portion of the genetic variation included in the north-south ancestry variable could also be explained by genetic differences for adaptive traits between the northern and southern groups of each species (Joyce & Rehfeldt, 2013;Rehfeldt et al., 1984;Richardson et al., 2009). Results were similar when we did not control for IBC in RDAs: A significant proportion of the variation could be attributed exclusively to IBD, but not to IBE, and the majority of the explained variation was confounded between IBD and IBE (not shown). Thus, we were unable to determine whether local adaptation has contributed to the genetic differentiation between the northern and southern groups in either species.  ulations (DeWoody et al., 2015;Papadopulos et al., 2014). Therefore, disentangling their relative contributions is very challenging, and attributing patterns of genetic variation to a single factor can be an oversimplification of the processes involved.

| Comparisons between GEA methods
Bayenv2 and LFMM control for population structure in different ways. Bayenv2 first estimates a covariance matrix of allele frequencies among populations and then tests for significant genotype-environment correlations using this covariance matrix as a null model.
LFMM estimates genotype-environment correlations while jointly estimating population structure via a number of latent factors (k, related to principal components). Although both methods essentially operate based on the same principles, that is, they test for GEAs after controlling for the portion of variation that is due to neutral population structure, their results differed greatly (P. strobus: 19% overlap; P. monticola: no overlap). The relative performance of Bayenv and LFMM depends on the demographic scenario and sampling design, and a relatively low overlap between the two methods has previously been observed in simulation studies (Lotterhos & Whitlock, 2015;de Villemereuil et al., 2014).
When the selective climatic gradients are highly collinear with neutral patterns of population structure, it becomes harder to separate neutral from selected loci, especially under weak selection (Lotterhos & Whitlock, 2015). For the subsets of SNPs detected by GEA methods in this study, the majority of the explained among-population variation (70.5 to 95.6%) was confounded between the effects of climate, geography, and north-south ancestry, leaving only a small proportion of the variation attributed exclusively to climate ( Figure 5). Depending on the underlying correction for population structure, GEA methods can attribute this confounded variation either to neutral structure (i.e., overcorrection resulting in false negatives) or to variation due to climate (i.e., undercorrection resulting in false positives

| Importance of the sampling strategy
The sampling design is one of the most important aspects to consider when looking for signatures of local adaptation (Meirmans, 2015).
In this study, we selected a large number of populations to cover a wide range of environmental variation across the natural distribu- Simulations showed that for a large total sample size (~900 diploids) there was little benefit in allocating individuals to more or less populations (Lotterhos & Whitlock, 2015). In a study conducted at a similar scale in P. strobus, Rajora et al. (2016) opted for sampling a larger number of individuals per population (22 individuals per population for SNP data) in 50 populations (total of 1100 individuals). They found only two significant SNP-environment correlations of 44 a priori candidate SNPs (4.5%, Spearman's rank correlation tests corrected for latitude and longitude). In comparison, our dataset included candidate genes (i.e., the 52 putative Picea glauca orthologs that were candidates for growth, phenology, and cold hardiness) as well as noncandidate genes.
The higher percentage of GEAs detected in our study (Pinus strobus: Bayenv2: 9.8%; LFMM: 13.7%) may reflect the greater power derived from sampling a wider range of environmental variation. We included a larger number of populations from the southern and western edges of the distribution that experience different temperature and precipitation regimes compared with the rest of the range. Moreover, we included eight populations (104 individuals) from the southern genetic group detected in Nadeau et al. (2015), whereas Rajora et al. (2016) included only one (North Carolina). Therefore, our GEA analyses may have detected genes that are involved in local adaptation (exogenous loci) or in maintaining an endogenous barrier between the southern and northern genetic groups (Bierne et al., 2011). The different GEA methods used are also likely to account for the differences between both studies for the reasons discussed in "Comparisons between GEA methods".
Because the effects of climate are confounded with IBD and IBC in P. strobus and P. monticola, sampling designs that minimize collinearity between environmental gradients and neutral population structure could greatly increase the power to detect signatures of local adaptation and IBE (Lotterhos & Whitlock, 2015;de Villemereuil et al., 2014;Wang & Bradburd, 2014). Sampling pairs of populations from contrasting environments, but which are relatively closely located in order to minimize differences in ancestry, showed higher power in simulations (Lotterhos & Whitlock, 2015). However, this strategy would be difficult to implement in species like P. strobus because climatic gradients and patterns of differentiation for adaptive traits occur at wide geographic scales, mostly along northward and westward postglacial colonization routes (Joyce & Rehfeldt, 2013;Nadeau et al., 2015;. In such case, increasing the total number of sampled individuals (Lotterhos & Whitlock, 2015) from many ecologically different populations (including climate extremes) may be the best strategy. Analyses performed within phylogeographic genetic groups would remove the confounding effect of IBC, but may also miss some important climatic adaptation between groups. If possible, replicated transects along climatic, edaphic or phenotypic gradients that are less correlated with the main axes of neutral population structure could increase the power to detect signatures of local adaptation.

| Highly supported candidate genes
Simulation studies showed that combining results from a number of methods can reduce false discovery rates (de Villemereuil et al., 2014) and detect loci under strong selection (Lotterhos & Whitlock, 2015).
By combining the results from F ST outlier tests and two GEA methods, we identified four and one highly supported candidate genes in P. strobus and P. monticola, respectively (Table 3). Putative orthologs of three of those genes were previously found to be important for growth and phenology in Picea glauca (El Kayal et al., 2011;Pelgas et al., 2011), one of which (CL3539-Contig1_01) was included in this study for this reason (i.e., part of the 23 candidate genes for growth, phenology, and cold hardiness, see "Materials and Methods"). We detected a larger number of highly supported genes in Pinus strobus than in P. monticola. These included a serine-threonine-protein kinase and a galacturonosyltransferase that were both found to be differen-

| Overlap of outlier loci between species
Three orthologous genes showed signatures of selection in both P. strobus and P. monticola (though only detected by a single method).
The discovery of genes that are involved in local adaptation to climate in both species could be expected given their relatively recent divergence (< 12 MYA) and the high degree of synteny among conserved orthologous genes in conifers (Pavy et al., 2012). One of those genes, a flavodoxin family protein (0_7001_01), was associated with the end of the frost-free period (eFFP) and temperature-related variables (DD5, TD) in both species. These climatic variables could be important drivers of local adaptation as phenotypic traits such as timing of budburst, timing of budset (or growth initiation and cessation), and cold hardiness vary significantly among populations of both species (Joyce & Rehfeldt, 2013;Joyce & Sinclair, 2002;Li et al., 1997;Lu et al., 2003a,b;Rehfeldt et al., 1984). The putative ortholog (i.e., gene amplified using the same primers as those used in this study) was also detected as an F ST outlier among environmental groups defined based on DD5 in June and precipitation in December in Larix decidua (Mosca et al., 2013), and associated with spring-fall precipitation and aridity in P. taeda (Eckert et al., 2010). Thus, this gene may have evolved in response to climatic constraints in multiple tree species.
Overall, we found that the number of genes carrying outlier SNPs in both species did not differ from random expectations and that the majority of outliers were species specific. In similar comparisons made between Picea mariana and P. glauca (divergence time >10 MYA), Prunier, Laroche, Beaulieu, and Bousquet (2011) found more adaptive similarities at the gene family level (paralogs) than at the gene level (orthologs). The redundancy of functions among recently duplicated genes in conifers could have offered the possibility for selection to act on paralogous genes in different species (Namroud et al., 2010).
In distantly related P. glauca and Pinus contorta (divergence time ~140 MYA), an exome-wide study detected 47 genes (~10-18% of top candidate genes) with convergent signatures of local adaptation to low temperatures (Yeaman et al., 2016). Paralogous genes in either species were more likely to show strong signatures of convergence than one-to-one orthologs. Yeaman et al. (2016) sidestepped the problem of overadjustment for population structure by using uncorrected genotype-environment and genotype-phenotype correlations in each species separately to identify top candidate genes. Then, they looked for enrichment of signatures of local adaptation between both species under the assumption that genetic drift is unlikely to affect the same genes similarly between distantly related species and give rise to the same false positives. In cases where patterns of adaptation covary with neutral population structure, this method is more powerful to identify convergent loci. However, it could not be used in the current study due to the relatively modest number of loci.

| CONCLUSIONS AND PERSPECTIVES
In this study, we attempted to disentangle signatures of local adaptation and IBE from those of IBD and IBC in two white pine species with different demographic histories (Nadeau et al., 2015). P. strobus could be considered an ideal species in which to look for signatures of local adaptation since it shows moderate among-population genetic variation for adaptive traits, but weak neutral population structure. We found in both species that a large amount of the explained amongpopulation genetic variation was confounded between the effects of climate (IBE), IBD, and IBC, with only a small proportion of the variation attributed exclusively to climate in P. strobus. Such confounding of patterns of local adaptation with neutral population structure is expected to be common in natural landscapes (e.g., Lasky et al., 2012;Lee & Mitchell-Olds, 2011;Sork et al., 2010). Two main reasons can explain these patterns: (1)  GEA methods can account differently for the confounded variation between the effects of climate, IBD, and IBC and, consequently, overcorrect or undercorrect for population structure. Combining results across a number of different methods should be standard practice to detect strong candidate genes (De Mita et al., 2013;Lotterhos & Whitlock, 2014de Villemereuil et al., 2014). The integration of phenotypic and genotypic information from populations growing in common-garden experiments might be the most informative approach to discover loci important for local adaptation (Sork et al., 2013;Yeaman et al., 2016), and it could be used to validate candidate SNPs detected by F ST outlier or GEA methods (e.g., De Kort et al., 2014;Jaramillo-Correa et al., 2015;Yoder et al., 2014). Another promising avenue would be taking advantage of the annual tree rings to establish relationships between annual growth and climatic variation in commongarden experiments over a number of years. A genotype-phenotype association study using variation in growth responses to climate among genotypes is underway to validate the outlier loci found in this study.