Landscapes and demographies

All simulations consisted of 100 demes arranged on a 10 x 10 landscape grid. The 15 levels of landscape-demography were broadly divided into three landscape categories: (i) a stepping stone landscape with latitudinal and longitudinal selective clines (Stepping-Stone Clines, the most commonly simulated scenario in testing methods) (Coop et al. 2010; Frichot et al. 2013; Günther & Coop 2013; de Villemereuil et al. 2014; Lotterhos & Whitlock 2015; Rellstab et al. 2015; Gautier 2015; Forester et al. 2016, 2018), (ii) a stepping stone landscape with one latitudinal cline and one non-linear longitudinal mountain range (Stepping-Stone Mountain, which left the potential for unique architectures to arise to the same selective pressure at different geographic locations), and (iii) an estuary landscape with a latitudinal and longitudinal selective clines (Estuary Clines, which simulated repeated independent bouts of adaptation analogous to oysters or sticklebacks that repeatedly colonize and adapt to isolated freshwater environments connected by gene flow in the marine environment). For simplicity, I refer to the latitudinal environment as Temperature and the longitudinal environment as Env2.

In summary, Stepping-Stone Mountain had a different environmental pattern than Stepping-Stone Clines but the same demography, while Estuary Clines had the same environmental pattern as Stepping-Stone Clines but different demography. The demographic parameters were chosen such that different landscapes achieved similar levels of neutral genetic differentiation and local adaptation. Within each of the three landscapes, 5 demographies were simulated that described the migration rates and effective population sizes on the landscape (see Supplemental Methods).

See the Supplemental Methods for a description of the multivariate continuous space simulations with 6 traits.

Genetic map

The Wright-Fisher simulations were based on a previously published quantitative genetic model and a genetic map (Lotterhos 2019). The genome consisted of 20 linkage groups each with 50,000 sites. The scaled recombination rate (N metapop r = 0.01) gave a resolution of 0.001 cM between proximate bases and a total length of 50 cM for each linkage group. This resolution mimicked a SNP chip, in which SNPs were collected across a large genetic map (Lotterhos 2019). The population-scaled neutral mutation rate was (N metapop μ= 0.001). QTNs could evolve on the first 10 linkage groups, while on the second 10 linkage groups only neutral loci could evolve.

Genetic Architecture and Stabilizing Selection

Mutation: Quantitative trait nucleotides (QTNs) contributed additively to the optimal phenotype for each individual without dominance. Three genic levels were simulated: oligogenic (few loci of large effect on the trait), moderately polygenic (dozens to hundreds of loci with intermediate effects), and highly polygenic (hundreds of loci with small effects). For QTN mutations under 1 trait or 2 traits without pleiotropy, the univariate effect size of a new QTN mutation was drawn from a normal distribution with a mean of 0 and standard deviation sigma QTN. For QTN mutations under 2 traits with pleiotropy, the bivariate effect size was drawn from a multivariate normal distribution with a standard deviation of sigma QTN for both traits and no covariance, which gave flexibility for mutations to evolve with effects on one or both traits. Thus, the distribution of effect sizes and linkage relationships among QTNs was allowed to evolve.

Pleiotropy: Within each genic category were 5 levels of pleiotropy and selection: (i) 1 temperature trait (which adapted to the latitudinal cline), (ii) 2 traits without pleiotropy and equal strengths of selection on both traits, (iii) 2 traits without pleiotropy and with weaker selection on the temperature trait, (iv) 2 traits with pleiotropy (QTNs could evolve effects on one or both traits) and equal strengths of selection on both traits, and (v) 2 traits with pleiotropy and with weaker selection on the latitudinal temperature trait.

Selection: The trait was subject to spatially heterogeneous stabilizing selection with the optimum for each location in space given by the environment. For each individual in each generation, the fitness was determined by a Gaussian function given the difference between the individual’s phenotype and the optimum at that location.

See trait equation details in the supplemental document: equation_details_genetic_architecture_and_stabilizing_selection.pdf

For information on burn-in, adding neutral loci with tree sequencing, filtering, and sampling, see the Supplemental Methods.

Quantifying the degree of local adaptation, divergence, and structure

For each replicate, the degree of local adaptation was measured as (i) the difference between population fitness in sympatry and allopatry following (Blanquart et al. 2013) and (ii) the correlation between the phenotype and environmental cline for each trait. Overall divergence (genetic differentiation) was calculated as Weir and Cockerham’s FST (Weir & Cockerham 1984) in OutFLANK (Whitlock & Lotterhos 2015). Population structure was estimated with a principal component analysis on the genotype matrix.

Quantifying trait and allelic clines

The degree of a trait cline was measured as Kendall's tau rank correlation coefficient between individual trait values and deme environment. The degree of an allele frequency cline was measured as Kendall's tau rank correlation coefficient (Kendall 1938, 1945) between deme allele frequency and deme environment, with significance being determined after Bonferroni correction based on the number of SNPs in the data. QTNs that were significant by this criteria were deemed “clinal QTNs.” The proportion of clinal QTNs excluded minor alleles with frequency < 0.01.

GEA performance

Latent-factor mixed models (LFMM) assess the linear relationship between genotype and environment while controlling for structure as latent factors. LFMM was implemented using the function lfmm2 in the R package LEA v.4.0.3 (Frichot et al. 2013; Frichot & François 2015; Caye et al. 2019). Redundancy analysis (RDA) and the partial RDA (pRDA) including a structure correction (conditional on the first 2 PC axes) were implemented using the ‘rda’ function in the R package vegan (Dixon 2003). See Supplemental Methods for details of implementation and correction for false discovery rate.

The performance of the association metrics were summarized as: (i) false discovery rate (FDR, proportion of outliers that are neutral, lower is better), (ii) true positive rate (TPR, proportion of QTNs that are significant outliers, higher is better), and (iii) the area under the precision-recall curve (AUC-PR, higher is better) (Lotterhos et al. 2022). In order to provide the most optimistic estimate of a method’s performance, the performance statistics were calculated by only including truly neutral loci unaffected by selection on linkage groups 11-20 and the QTNs.

Importance of clinal QTNs to local adaptation

First framework: A linear model was used to conduct the GWAS with individual trait value as the response variable and SNP genotype, PC1, and PC2 as explanatory variables. The proportion of GWAS hits that also showed clines with the environmental variable was compared to the known number of clinal QTNs. Second framework: The proportion of additive genetic variance (VA) for each QTN was calculated as the additive genetic variance for the focal QTN standardized by the total additive genetic variance following (Lotterhos 2019). The proportion of additive genetic variance (VA) explained by clinal QTNs was compared to a null expectation equal to the proportion of QTNs that were clinal. Third framework: I estimated the proportion of local explained by different subsets of QTNs: (i) QTNs with MAF > 0.01, (ii) clinal QTNs, and (iii) clinal QTNs inferred from latent factor mixed models that include a structure correction. For (ii) and (iii), a GEA model was performed for each environment, and then outlier QTNs were combined into a focal QTN set that was used for the local adaptation prediction. For each focal subset of QTNs, the counts of the derived allele were multiplied by the QTN effect size, summed to get a phenotype, and that phenotype was used in an in silico reciprocal transplant using the known phenotype-fitness function to estimate the degree of local adaptation. This estimate was then divided by the total degree of local adaptation (using all QTNs including those below the MAF threshold) to get an estimate of the proportion of local adaptation explained by that focal subset.