Evidence of widespread selection on standing variation in Europe at height-associated SNPs

Strong signatures of positive selection at newly arising genetic variants are well-documented in humans1–8, but this form of selection may not be widespread in recent human evolution9. Because many human traits are highly polygenic and partly determined by common, ancient genetic variation, an alternative model for rapid genetic adaptation has been proposed: weak selection acting on many pre-existing (standing) genetic variants, or polygenic adaptation10–12. By studying height, a classic polygenic trait, we demonstrate the first human signature of widespread selection on standing variation. We show that frequencies of alleles associated with increased height, both at known loci and genome-wide, are systematically elevated in Northern Europeans compared with Southern Europeans (p<4.3×10−4). This pattern mirrors intra-European height differences and is not confounded by ancestry or other ascertainment biases. The systematic frequency differences are consistent with the presence of widespread weak selection (selection coefficients ~10−3–10−5 per allele) rather than genetic drift alone (p<10−15).


Supplementary Tables
Supplementary Table 1 LRT results for forward genetic drift simulations given different combinations of input parameters s (the actual selective coefficient is w = s*β, where β is the estimated effect size on height) and N e (the effective population size). Supplementary Table 3A shows analysis using only the top ~1,400 SNPs of the simulated data, and Supplementary Analyses were carried out as described in Table 1 Analyses were carried out as described in Supplementary Table 4, but using individuals from the UK and Italy from the POPRES dataset to estimate Northern-and Southern-European allele frequencies. Analyses were carried out as described in Supplementary Table 4, but using individuals from the UK and Portugal from the POPRES dataset to estimate Northern-and Southern-European allele frequencies.

Supplementary
Supplementary Analyses were carried out as described in Supplementary Table 4, but using individuals from Switzerland (French) and Italy from the POPRES dataset to estimate Central-and Southern-European allele frequencies.
Supplementary Analyses were carried out as described in Supplementary Table 4, but using individuals from Switzerland (French) and Portugal from the POPRES dataset to estimate Central-and Southern-European allele frequencies. Analyses were carried out as described in Supplementary Table 4, but a genome-wide data set of ~56,000 SNPs was used. In each case, the model incorporating proportional selection effects (w = s*β) showed the best fit to the AF data. Northern-and Southern-European allele frequencies are estimated from the MIGen Harps (NEur) and Regicor (SEur) datasets. Analyses were carried out as described in Supplementary Table 11, but the Northern-and Southern-European allele frequencies were estimated using UK and Italian individuals from the POPRES dataset.

Supplementary Table 13. Sign test and mean NEur -SEur AF difference tests for directly genotyped height-associated SNPs or height-associated SNPs with proxies
Supplementary

Supplementary Figures and Legends
Supplementary Figure  NEur -SEur AF difference for sets of 500 independent (r 2 < 0.1) SNPs across the genome using POPRES data. SNPs were sorted by GIANT height association p-value. Shown in red is the curve of best fit, in purple the genome-wide mean NEur -SEur AF difference, and in blue the expected mean difference (y=0). NEur, Northern European. SEur, Southern European. AF, allele frequency. 1,000 simulations were conducted using the ~56,000 SNP MIGen dataset where for each simulation, GIANT effect size directions were randomly flipped so that ~50% were in the opposite direction than from the original data. For each of these simulations, we compared the fit to data using models incorporating selection (w = svs. drift alone using the likelihood ratio test (LRT) as described in Table 1. Over 1,000 simulations, the maximized s parameters were normally distributed around 0 (Supplementary Figure 6), and the distribution of LRT p-values were uniform, consistent with accepting the null hypothesis of drift alone (i.e. no selection). The false discovery rate at p < 0.05 within these 1,000 simulations was 4.7%. We performed within-family association analyses of height and the Northern-predominant alleles across the genome, ranking SNPs by the SEur AF-adjusted p-values (see Supplementary Methods), and using AF data from MIGen. Analyses were conducted as described in Figure 3, but using the MIGen dataset to define the Northern-predominant allele, and ranked by SEur AFadjusted GIANT height p-values. a, Regression coefficients are plotted on the y-axis for groups of 500 independent SNPs on the x-axis, with the curve of best fit shown in red, the curve of best fit for the GIANT effect sizes in purple, and y=0 in blue. b, Running averages of the regression coefficients (red and black filled circles) are plotted on the y-axis and compared against the running averages of regression coefficients from 1,000 analyses where phenotypes were Given a set of independent SNPs, we calculate the predicted population phenotypic differentiation by multiplying the effect size (in units of standard deviations; 1 standard deviation is approximately 6.5 cm, or 2.5 inches) of each height-increasing allele by the difference in allele frequency between NEur and SEur, and summing this over all the SNPs in the set. Based on ~1,400 SNPs with the best evidence of association with height from the GIANT dataset (see Supplementary Methods), we compared the population phenotypic differentiation predicted to be accounted for by these height-associated SNPs (red solid line) to that produced by 10,000 sets of randomly-drawn, frequency-matched SNPs presumed to be under drift (histogram, mean indicated by blue dashed line). The matched SNPs have the same ancestral allele frequency, and absolute magnitude of effect sizes as the height-associated SNPs, and their allele frequencies are drawn from the actual NEur and SEur AF distributions, but they have no bias in the direction of effect with respect to the NEur -SEur AF. The height-associated SNPs produced a predicted phenotypic differentiation significantly greater than that produced by neutrally evolving variants assumed to have the same genetic architecture (p = 3.9x10 -17 , based on the mean and standard deviation for the 10,000 matched sets). This suggests that the observations presented in this study are not simply the result of ascertaining a trait/population pair from the extreme end of a plausibly-sized set of neutrally evolving traits, but rather reflect the effects of selection. NEur, Northern European. SEur, Southern European. AF, allele frequency.

Supplementary
Supplementary Figure 10 Supplementary Figure 10. Comparisons of exact and linearly approximated AF differences. For selection coefficients, s, of 10 -6 and 10 -5 , exact AF differences were calculated at starting AFs from .01 to .99 over generation times of 100, 500, 1000, 1500 and 2000. For s = 10 -4 generation times of 50, 100, 200 and 400 were used, and for s = 10 -3 generation times of 20 and 50 were used. Linear approximations of these same AF differences were then also calculated as described in section 2.5. The results from these two approaches were then compared against one another, with exactly calculated AF differences on the x-axis and linear approximated AF differences on the y-axis. For all selection coefficients tested, the maximum error was 0.12% (mean error: 0.0152%). AF, allele frequency.

GIANT height association dataset
The GIANT (Genetic Investigation of ANthrometrics Traits) Consortium was a study that meta-analyzed the genome-wide association results from across 46 different cohorts of European descent containing ~129,000 individuals after quality control filters 6 . The cohorts in GIANT were largely European-American and Northern European (predominantly UK, Iceland, Finland, Sweden, Netherlands and Northern Germany) but also included cohorts from Southern Germany, Switzerland, Tyrolia, Croatia, Sardinia, and Spain. For adult human height, the GIANT consortium identified 180 loci that surpassed genome-wide significance 6 . The genome-wide association dataset included the average allele frequency across the 46 cohorts and the evidence of height association for 2,834,209 SNPs. We focused on the subset of the SNPs that were directly genotyped in either the MIGen study 7 or the POPRES dataset 8 (see below). All alleles are reported on the positive strand using HapMap phase 2 CEU panel 9 as reference.
For GIANT authors and affiliations, see section 5.

European allele frequency datasets
Two datasets were used to estimate Northern-and Southern-European allele frequencies.
The first dataset consisted of cohorts from the MIGen (Myocardial Infarction Genetics) well as Swedish individuals from the MIGen Malmo dataset. These cohorts were subjected to the same quality control measures as described above, and were down-sampled as needed in order to make the appropriate comparison in our analysis.

Framingham Heart Study
The Framingham Heart Study is a multi-generational family-based cohort based in

Modeling genetic drift and selection
To model the effect of drift and selection on the observed Northern-to Southern-European frequency differences, we first estimated the expected amount of allele frequency differences that could be attributed to selection using the following equation (see section 4 for derivation): where p is the ancestral allele frequency (estimated as the average of Northern-and Southern-European allele frequency), T is the number of generations since the two populations have split, and w is the selective pressure experienced by the population under different models of ongoing selection.
Note that the above equation for changes of allele frequency is only an approximation, as the changes in allele frequency are also a function of the allele frequency itself. However, this effect is negligible when the changes in allele frequency are very small between generations, as is the case here, since (1+x) T ≈ 1+Tx when Tx<<1. To confirm the validity of this assumption for the ranges of selection coefficients that we estimated, we calculated the exact changes in allele frequency expected for selection coefficients of 10 -5 and 10 -6 , for starting allele frequencies between .01 to .99 at generation times of 100, 500, 1000, 1500 and 2000. For a selection coefficient of 10 -4 , generation times of 50, 100, 200 and 400 were used, and for a selection coefficient of 10 -3 , generation times of 20 and 50 were used. At each generation, we calculate the expected allele frequency for the next generation conditional on the selection coefficient and the allele frequency of the current generation. We then compared the exact calculations with that obtained from the linear approximation (calculated above; Supplementary Figure 10). For all selection coefficients tested, the maximum error ((exact allele frequency -approximated allele frequent) / exact allele frequency) was 0.12% (mean error: 0.0152%).
Based on these results, we conclude that our approximation is sufficiently accurate for very small changes in allele frequency. This simplifying (and computationally less intensive) approximation was used to produce the LRT results in Supplemental Tables 3-12. Note that results presented in Table 1 were obtained using the exact calculation.
We examined two different models of selection: 1. constant selection, where w was set to a constant selective parameter, s; 2. selection modulated by effects on height, where w was set to the product of a constant selective parameter and the SNP's effect on height (obtained from GIANT dataset, in units of s.d.), s * β. We then determine ΔAF Sel over a range of selective parameters (0 to 0.1) across multiple generational times since divergence of the two populations (20 to 2,000). Though NEur and SEur populations are unlikely to have diverged as recently as 20 generations ago (~400-500 years), these values of T were included to account for the likely migration between NEur and SEur since divergence, which would decrease the effective T and hence increase the estimate of s. We performed a grid-like search over both parameter spaces: selection parameters were stepped through by increments of 1x10 -4 from 0.01 to 0.1, increments of 1x10 -5 from 1x10 -5 to 0.01 and increments of 1x10 -6 from 0 to 1x10 -5 , and generational times were stepped through by increments of 100 generations from 100 to 2000. For each combination of parameters, we subtracted the change in allele frequency attributable to selection from the observed Northern-to Southern-European allele frequency differences, and determined the likelihood of observing the remaining allele frequency differences due to drift as described above. The parameters that maximized the likelihood of each model for each generational time since divergence were determined. The maximum likelihood parameter for all three models incorporating selection were compared against one another as well as the model with drift-alone, in order to determine the model that best fit our observed data using the likelihood ratio test (LRT). The likelihood ratio statistic (twice the difference in negative log likelihoods) was compared with a χ 2 distribution with the degrees of freedom equal to the difference in the number of free parameters. We also explored a third model, w = s 1 + (s 2 * β). However, in this case the maximum likelihood estimate for s 1 is zero across all generational times tested, so this class of models was not considered further.
We performed a number of simulations to test the behavior of our model. Specifically, we tested the performance of the LRT in data sets where we randomly reversed half of the signs of the effects on height, thereby creating demographically-matched data sets in which there is no effect of selection (described in section 3.4). We also performed forward simulations of populations under drift with or without selection, to test whether our LRT procedure could accurately estimate the selection parameters, and whether the method had sufficient power to detect widespread weak selection with coefficients in the ranges estimated by our data (described in section 3.5). Finally, in our modeling, we assumed a simplistic process of concurrent selection at many sites, which does not account for the more complex effects of multiple selected loci that ensue as allele frequencies begin to change substantially 14 . However, because the observed frequency differences are typically small (mean frequency difference = ~0.01 even for the most strongly associated SNPs), we can still use our models to roughly estimate the typical selective pressures on height-associated variants that would be consistent with our data. We considered whether the results of our analysis would remain the same if the GIANT consortium had been composed entirely of Southern-European populations, and the SNPs had the same effects on height in Northern Europeans and Southern Europeans. To do this, we adjusted the height-association p-value for each SNP in the GIANT dataset by assuming that the allele frequency in GIANT was equal to the estimated Southern-European allele frequency, but effect sizes remained the same. Specifically, for each SNP in the GIANT dataset, we first converted the height-association p-values into quantiles of a standard normal distribution. As the quantiles (or z-statistics) are functions of the effect size (unchanged) and the square root of the heterozygosity in the population, we adjust the GIANT z-statistics by (pq SEur /pq GIANT ) 1/2 , where p and q are the allele frequencies of the major and minor allele for that SNP. The adjusted zstatistics was then converted into Southern-European frequency-adjusted p-values to obtain the hypothetical evidence of association for each SNP, under the scenario where GIANT had studied predominantly Southern-European individuals.
To recreate a list of height loci that would surpass genome-wide significance in this hypothetical scenario, we then followed the analytical steps described by GIANT 6 . Briefly, we clumped our new dataset with an r 2 ≥ 0.1 using the most significantly associated SNP as our lead SNP in each region, and then pruned remaining SNPs within 0.5 Mb of one another, preferentially keeping the SNPs with best association to height within each clump. We identified SNPs with Southern-European frequency-adjusted p-values lower than that of the least significantly associated height loci among the published 180 loci in the original GIANT dataset.
This process generated a list of 185 SNPs (using frequency data from MIGen) and 156 SNPs (using frequency data from POPRES) that in theory would have passed genome-wide significance threshold had the GIANT consortium used predominantly individuals from Southern Europe.
We repeated all of our analyses using the new set of height SNPs and also using a new pvalue-ordered genome-wide set of SNPs based on Southern-European frequency-adjusted pvalues (Supplementary Table 1 and Supplementary Figures 4-5). In general, we observed the same pattern of height-increasing allele being more prevalent in Northern Europeans than in Southern Europeans, suggesting that the variant ascertainment due to the largely Northern-European ancestry of the GIANT consortium did not confound our results.

HapMap SNP ascertainment bias
We also considered the effects of identifying height-associated SNPs from the set of SNPs in HapMap, which were ascertained in part by sequencing in Northern-but not Southern-European samples. We used data from the 1000Genomes Project 15 to simulate two different SNP ascertainment schemes. We catalogued all SNPs present in the low-pass whole genome sequence data in either four unrelated CEU samples from 1000Genomes or in two CEU, 1 YRI and 1 CHB sample. We then drew SNPs from this pool, assembling 1,000 sets of SNPs matched to the ~1,400 height-associated SNPs by their Northern-European minor allele frequency in the 1000Genomes GBR panel. We also assigned the "height-increasing" allele of each matched SNP to be either the GBR major or minor allele, depending on whether the actual heightassociated SNP to which it was matched was the Northern-European major or minor allele. We then calculated Northern-to Southern-European allele frequency differences for the actual set of height-increasing alleles and for the 1000 matched sets of "height-increasing" alleles. For neither ascertainment scheme was the frequency difference of any of the matched sets of "heightincreasing" alleles greater than the frequency difference for the actual height-increasing alleles.
We also calculated the mean and standard deviation of the Northern-to Southern-European allele frequency difference data for the matched sets of SNPs, and used these to calculate a Z-score for the value of frequency difference for the actual height-increasing alleles. The Z-scores were 6.04 and 6.05 for the two ascertainment schemes, corresponding to p values of 7.5x10 -10 and 7.3x10 -10 .

Phenotype ascertainment bias
Another potential bias to our analysis is our ascertainment on a phenotype known to be differentiated between Northern and Southern Europe. As had been discussed by Orr 16 , because we selected a phenotype known to be differentiated between two populations, it might not be surprising that we observed more height-increasing alleles in the taller population, when compared to expected parity. We addressed this potential ascertainment bias by demonstrating that the population height differences between Northern and Southern Europeans are extremely unlikely to be obtained due to just the effects of genetic drift alone, given the known genetic architecture of height.
We performed 10,000 neutral trait simulations starting with the set of 1,437 most strongly height-associated SNPs described in Online Methods. For each simulation, we randomly selected a set of 1,437 variants in the genome, each matched to the corresponding heightassociated SNP by estimated ancestral allele frequency (the average of Northern-and Southern-European allele frequency estimates). Each of these randomly drawn frequency-matched variants was assigned the height effect size from the height-associated variant to which the randomlydrawn variant was matched. We then calculated the level of population phenotypic differentiation that would be predicted to arise from drift, based on the Northern-to Southern-European allele frequency difference of the randomly-drawn SNPs and their assigned effect sizes. We calculate this as: Phenotypic differentiation = Σ β * (NEur-SEur AF). The distribution of phenotypic differentiation over 10,000 sets of matched SNPs is then compared to the phenotypic differentiation calculated using the actual height SNPs. Based on the distribution of phenotypic differentiation for randomly-matched SNPs, we observed that the actual height SNPs are ~8.3 standard deviations away from the mean (Supplementary Figure 9), corresponding to a probability of 10 -17 of observing this strong a phenotypic differentiation between Northern and Southern Europe in a neutrally evolving trait with the genetic architecture of height. In other words, our observed result in height would still be statistically significant after correction for testing multiple phenotypes even had we ascertained height as the most differentiated phenotype among ~10 16 neutrally evolving phenotypes.

Randomized GIANT effect size direction LRT simulations
To test whether our LRT procedure has the appropriate statistical properties, we first estimated the type-1 error rate of our procedure. We started with the actual Northern-and Southern-European allele frequency data and GIANT effect sizes for the set of 1,437 SNPs described in Online Methods, and generated 1,000 matched data sets where there is no effect of selection on the height SNPs. Specifically, for each matched data set, the directions of the GIANT effect sizes for a random ~50% of the analyzed SNPs were reversed, thereby eliminating selection, but retaining all of the other features of our observed data set. We then performed our analyses, comparing a model of w = s*β vs. a model of drift alone, and estimated the maximum likelihood estimates of s and calculated the corresponding LRT p-values. If our method were well-behaved, we expect that only ~5% of the simulations should have p-values < 0.05, and the estimated selection parameters should be normally distributed around zero, with magnitudes much less than we observed in the real data. In 1,000 matched data sets, the false discovery rate was ~4.7% (see Supplementary Figure 5 for QQ-plot for p-values), the estimated selection coefficients were normally distributed around zero, and were typically at least an order of magnitude lower than those estimated from the actual data (Supplementary Figure 6). Thus, our approach has the appropriate type-1 error rate on data sets directly comparable to our real data, and, under the null, produces an unbiased set of estimated selective coefficients.

Forward genetic drift simulations
To test whether our LRT procedure could accurately estimate selection parameters in the ranges that we observed, and whether our method had sufficient power to detect widespread weak selection with coefficients in the ranges estimated by our data, we conducted forward simulations of drift and selection. We used the average allele frequencies of Northern and Southern Europeans for ~56,000 independent SNPs genome-wide, representing an ancestral population, as the starting point for our simulations. The ancestral population was divided into three equal daughter populations of size N e. The SNPs in the first of these simulated ancestral populations underwent drift with no selection. The SNPs in the second population underwent both drift and selection (with selective coefficients proportional to effect size, w = s*β). In the third population, only the 180 height SNPs underwent drift and selection, in order to calibrate the number of generations for which selection should act in the second population. We allowed the simulations to proceed until the 180 height SNPs in the third population reached the 1.2% mean frequency difference observed in the original data. Once this mean frequency difference was met, selection was "turned off" in the second population. We then, if necessary, allowed drift to continue in each of the first two populations until the genome-wide F ST between these two simulated data sets matched the F ST for the observed data. The allele frequency data from these first two simulated populations were used in the analyses. By this procedure, each pair of simulated data sets is matched to the actual Northern-and Southern-European data by the allele frequency spectrum, by the average differentiation of the known height-associated SNPs and also by the overall F ST between NEur and SEur.
We varied the input parameters over a range of s and N e. Five replicate simulations were performed for each combination of input parameters. For each replicate, we then determined if a model incorporating drift plus selection or a model with drift alone better fit the data by calculating the LRT statistic, using either the top ~1,400 SNPs or the larger genome-wide data set of ~56,000 SNPs. We estimated the selective parameter s using the likelihood ratio approach as described in section 2. For each combination of input parameters, we report the results for the median of these five replicates (Supplementary Table 2). For all combinations of parameters tested, the LRT p values were strongly significant (and comparable to those observed in the actual data) and estimated selective coefficients were extremely close to the input parameters used, thereby confirming both the validity and the power of our approach. 4 Derivation of the allele frequency difference expected due to weak selection To calculate the expected allele frequency difference between two populations after selection for T generations, we begin with the frequency of an allele after one generation of selection, namely (1) where p' is the frequency after one generation of selection, p is the frequency before selection, q is 1p, and w 11 , w 12 , and w 22 , are the selection pressures experienced by individuals homozygous for allele 1, heterozygous for the two alleles, and homozygous for allele 2, respectively. We assumed that allele 1 is under selection, thus the relative fitness for w 11 , w 12 , and w 22 are 1+2w, 1+w, and 1, respectively, for some constant w (ref 17).
If the total effect of selection is small (as suggested by our small observed AF differences at height-associated SNPs), the difference in frequency after T generations between two populations, one of which is under selection, is approximately We can replace p' with equation 1 to get: