Reduced signal for polygenic adaptation of height in UK Biobank

1In recent years, there has been great progress in understanding the polygenic basis of a wide variety of complex traits. One significant development has been the advent of ‘polygenic scores’, which aim to predict the additive genetic component of individual phenotypes using a linear combination of allelic contributions to a given trait across many sites. An important application of polygenic scores has been the study of polygenic adaptation—the adaptive change of a phenotype through small allele frequency shifts at many sites that affect the phenotype.

2Thus far, the clearest example of polygenic adaptation in humans has come from analyses of polygenic scores for height in Europe. However, as we will show here, this signal is strongly attenuated or absent using new data from the UK Biobank, calling this example into question.

3Starting in 2012, a series of papers identified multiple lines of evidence suggesting that average polygenic scores for height increase from south-to-north across Europe (Table 1). Analyses from multiple groups have concluded that the steepness of this cline is inconsistent with a neutral model of evolution, suggesting that natural selection drove these differences in allele frequencies and polygenic scores. Significant differences in polygenic scores for height have also been reported among ancient populations, and these are also argued to have been driven by selection. In parallel, developed the Singleton Density Score (SDS)—which compares the distance to the nearest singleton on two alternative allelic backgrounds—to infer recent changes in allele frequencies, and used it to analyze a large sample of British individuals (the UK10K; ). They found a significant covariance of SDS and effect on height, suggesting that natural selection drove a concerted rise in the frequency of height-increasing alleles in the ancestors of modern British individuals during the last 2,000 years.

4All such studies rely on estimates of individual allelic effects on height, as calculated from genome-wide association studies (GWAS). These GWAS estimates are then combined with population-genetic analysis to test for selection. Under a null model of no directional change, we would not expect ‘tall’ alleles to increase (or decrease) in frequency in concert; thus, loosely speaking, a systematic shift in frequency of ‘tall’ alleles in the same direction has been interpreted as evidence for selection.

5While our focus here is on the the distribution of height polygenic scores in Europe, we see this as a case study for understanding the challenges in comparing polygenic scores across populations in general. Compared to other complex traits, height is particularly well-characterized, and the evidence for adaptation of height in Europeans seemed clear. Thus our work highlights a need for caution in this area until these issues are more fully understood.

6Until recently, the largest height GWAS dataset came from the GIANT consortium (253,288 individuals as of 2014; ). This is the primary GWAS underlying most studies of adaptation of height. Additionally, several groups have used other, smaller, datasets to replicate signals found using GIANT. In particular, because it is known that population structure may be a confounder in GWAS, leading to false positive inferences of polygenic adaptation, several groups sought to replicate signals using family-based analyses, which protect against confounding due to stratification.

7The first replication, by, showed that the effect sizes of the top 1,400 GIANT associations (based on an earlier version of GIANT, published by were statistically consistent with effect sizes re-estimated in a smaller sibling-based regression approach using data from the Framingham Heart Study (4819 individuals across 1761 nuclear sibships from ). Sibling-based regression is considered to be immune to confounding by population structure, and so the agreement of effect sizes between studies was taken as validation of the north-south gradient observed when using the GIANT effect sizes.

8The second, partially independent, replication came from, who selected 691 height-associated SNPs on the basis of the GIANT association study, and then computed polygenic scores using effect sizes re-estimated in a cohort of 6307 individuals of Sardinian ancestry. They determined that the average polygenic score of Sardinian individuals was significantly lower than observed for other European populations, consistent with the previously reported north-south gradient of polygenic scores.

9A third replication was performed by, who used a different, larger sibling-based GWAS to identify associations ( sibling pairs from ). We refer to this sibling-based dataset as ‘R15-sibs’. These authors showed that the north-south frequency gradient replicates using SNPs ascertained from the sibling-based GWAS. This replication is stronger than that performed by either, or that by, as the cohort is larger and the SNP ascertainment did not rely on GIANT. As pointed out in the supplementary note of, this two-step procedure—ascertaining with a large but potentially biased GWAS like GIANT, before switching to a less powerful but hopefully unbiased replication GWAS—has the potential to introduce an ascertainment bias, even if the effects are correctly estimated in the replication study (we note that a small fraction of the GIANT samples are contained within the R15-sibs analysis, so the effect sizes are not strictly independent; however, because of the sibling design, any bias due to stratification in GIANT should be absent in R15-sibs). The R15-sibs study was also used by to verify a signal of recent selection in ancestors of the present day British population. found that the signal of selection was fully replicated when using R15-sibs data.

10Lastly, also used LD Score regression to test for height adaptation in the British while controlling for population structure. While LD Score regression is typically used to estimate genetic covariance between two phenotypes, used it to test for a relationship between height effects and a recent increase in frequency (measured by SDS)—and found a strong covariance of the two consistent with selection driving allele frequency change at height loci.

11Here, we reassess these previously reported signals using data from the UK Biobank (UKB) with genotype and phenotype data for nearly 500,000 residents of the United Kingdom. UKB has recently become a key resource for GWAS, thanks to its large sample, the relatively unstructured population (compared to international studies such as GIANT), and the opportunity for researchers to work directly with the genotype data rather than with summary statistics.

12This paper has two aims. First, we will show that previously reported signals of directional selection on height in European populations generally do not replicate when using GWAS effect estimates from the UK Biobank. Similar findings have been obtained independently by other groups working in parallel. Second, we will show that both the GIANT and R15-sibs GWAS are confounded due to stratification along the North-South gradient where signals of selection were previously reported. Signals detected using R15-sibs effect estimates were previously used as a significant source of evidence in favor of adaptation by, as well as in. However, the investigators leading the GWAS have now confirmed that the effect size estimates released from their 2015 study were strongly affected by population structure due to a computational bug. We include an analysis of the R15-sibs GWAS here to document how these spurious signals affected previous inferences, as well as the evidence that indicated the presence of confounding in the data prior to detection of the bug.

13The conclusion that adaptation signals in GIANT were spurious has broader implications for GWAS analysis, as it indicates that standard approaches for population structure correction may not always be sufficient, and that further study is needed to understand their limitations. While we anticipate that current methods are likely adequate for many applications, in particular for identification and broad-scale localization of strong genotype-phenotype associations—they may be insufficient for applications such as phenotypic prediction and the detection of polygenic adaptation as these can be highly prone to the cumulative bias through uncorrected structure. Such analysis should be undertaken with great care.

14We downloaded or generated seven different height GWAS datasets, each relying on different subsets of individuals or using different analysis methods. The bold-faced text give the identifiers by which we will refer to each dataset throughout this paper. These include two previous datasets that show strong evidence for polygenic adaptation, as well as an updated version of the R15-sibs dataset released in response to results in the initial preprint version of this manuscript:

15GIANT: ( = 253 k) 2014 GIANT consortium meta-analysis of 79 separate GWAS for height in individuals of European ancestry, with each study independently controlling for population structure via the inclusion of principal components as covariates.

16R15-sibs: ( = 35 k) Family-based sib-pair analysis of data from European cohorts. The effect sizes associated with this paper were publicly released in 2016.

17R15-sibs-updated: ( = 35 k) In November 2018, while this paper was in the final stages of revision, the authors of the R15-sibs data reported that their earlier data release failed to correct properly for structure confounding. They released this corrected dataset as a replacement.

18We also considered four different GWAS analyses of the UK Biobank data, using different subsets of individuals and different processing pipelines:

19UKB-GB: ( = 337 k) Linear regression controlling for 10 principal components of ancestry (unrelated British ancestry individuals only).

20UKB-Eur: ( = 459 k) All individuals of European ancestry, including relatives. Structure correction was performed using a Linear Mixed Model (LMM) approach, which controls for genetic stratification effects by modeling the genetic background as a random effect with covariance structure given by the kinship matrix. Mild amounts of environmental stratification are subsumed into this term, and therefore controlled implicitly.

21UKB-GB-NoPCs: ( = 337 k) Linear regression without any structure correction—with only genotype, age, sex and sequencing array as covariates (unrelated British ancestry individuals only) [newly calculated by us, see Materials and methods].

22UKB-sibs: ( = 35 k) Family-based sib-pair analysis [newly calculated by us, see Materials and methods].

23To understand the extent to which these different datasets capture a shared signal, we treated each set of summary statistics as if it were derived from a GWAS of a different phenotype and estimated the genetic correlation between them using bivariate LD Score regression. We find that all of these studies show high pairwise genetic correlations, consistent with the view that all of them estimate a largely-similar genetic basis of height (Table 2).

24One well-studied signal of adaptation of height in Europe has been the observation that, among height-associated SNPs, the ‘tall’ alleles tend to be at higher frequencies in northern populations. Equivalently, the average polygenic scores of individuals in northern populations tend to be higher than individuals in southern populations. To evaluate this signal for each dataset, we independently ascertained the SNP with the smallest p-value within each of 1700 approximately independent LD blocks (subject to the constraint that MAF 0.05 within the GBR 1000 Genomes population). We used these loci to calculate average polygenic scores for each of a set of European population samples taken from the 1000 Genomes and Human Origins panels (see Materials and methods for statistical details).

25As expected, we find highly significant latitudinal gradients in both the GIANT and R15-sibs data. However, this signal does not replicate in any of the four UK Biobank datasets (Figure 1, top row), or in the R15-sibs-updated dataset (Figure 1—figure supplement 1).

26We also tested whether the polygenic scores are over-dispersed compared to a neutral model, without requiring any relationship with latitude (the test from ). Here we find a similar pattern: we strongly reject neutrality using both the GIANT and R15-sibs datasets, but see little evidence against neutrality among the UK Biobank datasets, or the R15-sibs-updated dataset (Figure 1). The sole exception is for the UKB-GB dataset, though the rejection of neutrality in this dataset is marginal compared to that observed with GIANT and R15-sibs, and it does not align with latitude.

27While most studies have focused on a latitudinal cline in Europe, a preprint by also recently reported a cline of polygenic scores decreasing from west to east across all of Eurasia. Extending this analysis across all six datasets, we observe similarly inconsistent signals (Figure 1, bottom row). Only the GIANT dataset shows the clear longitudinal signal reported by, though the R15-sibs dataset is again strongly over-dispersed in general, and retains some of the longitudinal signal. Interestingly, we find a weakly significant relationship between longitude and polygenic score in the UKB-Eur dataset (though not in the other UKB datasets), suggesting there may be systematic differences between the results based on British-only and pan-European samples, even when state of the art corrections for population structure are applied.

28We also experimented with a larger number of SNPs using a procedure similar to (Appendix 1). We found that in some cases this led to significant values of when using UKB-GB effect sizes to ascertain SNPs. However, this signal was sensitive to the particular method of ascertainment, and seems to be diffuse (i.e. spread out across all axes of population structure, Appendix 1–figure 6), with part of the signal coming from closely linked SNPs. Thus we conclude that this signal is not robust and may, at least partially, arise from a violation of the assumption of independence among SNPs that underlies our neutral model. We also tested different frequency, effect size and probability-of-association cutoffs to determine which SNPs we included in the computation of the scores, but found none of these cutoffs affected the discrepancy observed between the GIANT and UKB-GB datasets (Appendix 2).

29We next evaluated the Singleton Density Score (SDS) signal of selection for increased height in the British population, previously reported by. SDS estimates recent changes in allele frequencies at each SNP within a population by comparing the distances to the nearest singleton variants linked to each of the focal SNP’s alleles. applied SDS calculated across the UK10K sample to investigate allele frequency changes in the ancestors of modern British individuals. SDS can be polarized according to the sign of a GWAS effect at each SNP–this is denoted trait-SDS, or tSDS. Here, tSDS 0 indicates that a height-increasing allele has risen in frequency in the recent past; tSDS 0 correspondingly indicates a decrease in frequency of the height-increasing allele. A systematic pattern of tSDS 0 is consistent with directional selection for increased height.

30Using both GIANT and R15-sibs, found a genome-wide pattern of positive tSDS, indicating that on average, height-increasing alleles have increased in frequency in the last ∼75 generations. tSDS also showed a steady increase with the significance of a SNP’s association with height. We replicate these trends in Figure 2A,B.

31This tSDS trend is greatly attenuated in all four GWAS versions performed on the UK Biobank sample (Figure 2C–F), as well as the R15-sibs-updated dataset (Figure 2—figure supplement 1). The correlation between UKB-GB GWAS p-value and tSDS is weak (Spearman, block-jackknife ). This correlation is stronger for the UKB-Eur GWAS. Since the UKB-Eur GWAS is not limited to British individuals—but instead includes all European ancestry individuals—this might suggest that residual European population structure continues to confound UKB-Eur effect estimates, despite the use of LMM correction for structure, similar to the longitudinal signal detected above for this same dataset (Figure 1).

32We wondered whether the main reason for the weakened trend in UKB-GB is an overly conservative PC-correction. This could occur if the genetic contribution to height is highly correlated with population structure axes. If this were the case, we would expect the correlation between GWAS p-value and tSDS to still be observed in a UKB GWAS without population structure correction (namely, in UKB-GB-NoPCs). However, we see no evidence for this correlation (block jackknife p = 0.6). Taken together with the UKB-GB-NoPCs polygenic score analysis (Figure 1), the lack of signal in UKB-GB-NoPCs suggests that the main reason that UKB is less confounded by population structure than GIANT is the relatively-homogeneous ancestry of the UKB British sample—rather than differences in GWAS correction procedures.

33Lastly, we examined tSDS at the most significant height-associated SNPs of each UKB GWAS (as before, ascertained in approximately-independent LD blocks). Significant SNPs show a positive average tSDS (Figure 2I,K,L; t-test p 0.05)—with the exception of the UKB-GB-NoPCs GWAS (Figure 2J) in which the average tSDS is not significantly different from zero (t-test p = 0.6).

34We have now shown that signals of polygenic adaptation of height are greatly reduced in the UKB data relative to the GIANT and R15-sibs datasets. To better understand the differences among the datasets, we ascertained 1,652 approximately-independent lead SNPs based on the GIANT p-values to form the basis of comparison between the GIANT and UKB-GB datasets.

35Figure 3A shows the effect sizes of ancestral alleles, as estimated using GIANT (x-axis) and UKB-GB (y-axis). The two datasets are highly correlated, consistent with the strong genetic correlation estimated in Table 2. The fact that the slope is probably reflects, at least in part, the standard winner’s curse effect for SNPs ascertained in one study and replicated in another.

36Importantly however, we also see clear evidence that the differences between the GIANT and UKB-GB effect sizes are correlated with European population structure (Figure 3B). Specifically, for each SNP we plotted the difference in effect size between GIANT and UKB-GB against the difference in allele frequency between northern and southern European samples (specifically, between the British (GBR) and Tuscan (TSI) subsets of 1000 Genomes). These differences have a significant correlation, indicating that alleles that are more frequent in GBR, compared to TSI, tend to have more positive effect sizes in GIANT than in UKB-GB, and vice versa. We also observed a similar signal for frequency differences between TSI and the Han Chinese in Beijing (CHB, Figure 3—figure supplement 1), suggesting that longitudinal patterns observed by were also likely driven by incompletely controlled stratification in GIANT.

37Similar patterns are present in a comparison of the R15-sibs and UKB-GB datasets when ascertaining from R15-sibs p-values (Figure 3, panels C and D; 1,642 SNPs). Here, the correlation between effect size estimates is much lower, likely due to the much smaller sample size of R15-sibs, and therefore elevated winner’s curse effects. However, the correlation between the effect-size difference and the GBR-TSI allele frequency difference remains. In contrast, when SNPs are ascertained on the basis of their UKB-GB p-value, these patterns are considerably weaker in both the GIANT and R15-sibs datasets (Figure 3—figure supplement 2).

38Finally, an unexpected feature of the R15-sibs dataset can be seen in Figure 3D: there is a strong skew for the ancestral allele to be associated with increased height (1,201 out of the 1,642 SNPs ascertained with R15-sibs p-values have positive effect sizes in R15-sibs). This pattern is not present in the R15-sibs-updated dataset (851 out of 1,699 SNPs with positive effects), or any other dataset we analyzed, suggesting that it is likely a result of the failure to control for population structure.

39Together, these observations suggest that while all of the datasets primarily capture real signals of association with height, both the GIANT and R15-sibs effect size estimates suffer from confounding along major axes of variation in Europe and Eurasia. This could drive false positive signals in geographic-based analyses of polygenic adaptation. Furthermore, since SDS measured in Britain correlates with north-south frequency differences, this could also drive false positives for SDS.

40To explore this further, we next turn to an analysis of the datasets based on LD Score regression.

41Another line of evidence from came from LD Score regression. LD Score regression applies the principle that, under a polygenic model, SNPs in regions of stronger LD (quantified by a SNP's ‘LD Score’) should tag more causal variants and therefore have larger squared effect size estimates. Similarly, if two traits share a genetic basis, then the correlation between estimated effect sizes of these traits should increase with LD Score. Meanwhile, confounders such as population structure are argued to affect SNPs of different LD Score equally, and therefore affect the intercept but not the slope of a linear regression to LD Score (we return to this point below; ).

42While LD Score regression is commonly used to estimate the genetic covariance between pairs of phenotypes, used it to test for a relationship between height and SDS. SDS is similar to GWAS effect estimates in that the expected change in frequency of an allele depends on both direct selection it experiences due to its own fitness effect as well as correlated selection due to the effects of those in linkage disequilibrium with it. predicted that the covariance between estimated marginal height effect and SDS should increase with LD Score—and found this to be the case using both GIANT and R15-sibs. This provided further evidence for polygenic adaptation for increased stature in Britain.

43Here, we revisit ’s observations (Figure 4A,B). Both GIANT and R15-sibs exhibit a highly significant LD Score regression slope (scaled GIANT slope,; scaled R15-sibs slope, ), as well as a highly significant intercept (GIANT intercept,; R15-sibs intercept, ). These large intercepts suggest that both GWAS suffer from stratification along an axis of population structure that is correlated with SDS in the British population. In contrast, in LD Score regression with the UKB-GB GWAS, the intercept is not significant, suggesting that UKB-GB is not strongly stratified (or at least, not along an axis that correlates with SDS). The slope is ∼1/3 as large as in GIANT, though still modestly significant (, Figure 4C). There is no significant slope or intercept in R15-sibs-updated (Figure 4—figure supplement 1B), and analyses of other UKB datasets give similar results to those for UKB-GB (Figure 4—figure supplement 2).

44We next extended Field et al.’s LD Score rationale from SDS to test whether SNP effects on height affected allele frequency differentiation between northern and southern Europe. We used 1000 Genomes British (GBR) and Tuscan (TSI) samples as proxies for northern and southern ancestry respectively. To control for the correlation between allele frequencies and LD Score, we normalized the frequency differences to have variance 1 within 1% average minor allele frequency bins. For shorthand, we refer to this measure as [GBR-TSI]. Under a model of selection driving allele frequency differences, we would expect the covariance of [GBR-TSI] and effect sizes to increase with LD Score. To test this, we regressed the product [GBR-TSI] effect size (estimated in previous and UKB GWAS) against LD Score.

45In contrast to SDS, we find that none of the GWAS datasets show a strongly positive slope (Figure 4D–F): the slope is approximately zero in GIANT, weakly positive in R15-sibs, and weakly negative in UKB-GB Results were similar for the other UKB datasets (Figure 4—figure supplement 3), and for R15-sibs-updated (Figure 4—figure supplement 1A). We see extremely strong evidence for positive intercepts in GIANT (p = 4) and R15-sibs (p = 9), but not in UKB-GB, R15-sibs-updated or the other UKB datasets. The large intercepts in GIANT and R15-sibs are consistent with stratification affecting both of these GWAS, as the North-South allele frequency difference is systematically correlated with the effect sizes in these GWAS independently of LD Score (see the Materials and methods for a more technical discussion). However, the relative lack of slope in these analyses suggests that the LD Score signal for SDS must be driven by a component of frequency change that is largely uncorrelated with the [GBR-TSI] axis of variation.

46The original LD Score regression paper noted that in the presence of linked selection, allele frequency differentiation might plausibly increase with LD Score. However, they concluded that this effect was negligible in the examples they considered, and subsequent applications of the LD Score regression approach have generally assumed that the two are independent. We find that in bivariate LD Score analyses of SDS, both the intercept and the slope differ significantly from zero for precisely the same GWAS datasets that show evidence of stratification (GIANT and R15-sibs), while both the slope and intercept are much reduced in all of the UK Biobank datasets. This suggests that the LD Score regression slope may not be as robust to stratification as hoped, prompting us to revisit the assumption of independence.

47We find that the squared allele frequency difference [GBR-TSI]2 is significantly correlated with LD Score (, Figure 5A), as are squared allele frequency contrasts for much lower levels of differentiation (i.e. between self-identified ‘Irish’ and ‘White British’ individuals in the UK Biobank (, Figure 5—figure supplement 1), and SDS2 (, Figure 5B). Strikingly, squared measures of more recent allele frequency change (i.e. SDS2 and the squared Irish vs. White British contrast) are much more tightly correlated with LD Score than that of more diverged populations [GBR-TSI]2, suggesting that the LD Score regression slope may be equally vulnerable to stratification involving closely related populations than for those that are more distantly related. Finally, the product [GBR-TSI] SDS is also correlated with LD Score (Figure 5C), demonstrating that the general signal of greater allele frequency change in regions of stronger LD is also shared between [GBR-TSI] and SDS.

48As noted above, the correlation we observe between allele frequency differentiation and LD Score could be generated by the genome-wide effects of linked selection. While a range of different modes of linked selection likely act in humans, one of the simplest is background selection. Background selection (BGS) on neutral polymorphisms results from the purging of linked, strongly deleterious alleles. Because any neutral allele that is in strong LD with a deleterious mutation will also be purged from the population, the primary effect of BGS is a reduction in the number of chromosomes that contribute descendants in the next generation. The impact of BGS can therefore be approximately thought of as increasing the rate of genetic drift in genomic regions of strong LD relative to regions of weak LD. Therefore, SNPs with larger LD Scores will experience stronger BGS and a higher rate of genetic drift, and this effect could generate a positive relationship between LD Scores and allele frequency differentiation.

49In Appendix 3, we derive a simple model for the effect of BGS on the relationship between allele frequency divergence and LD Scores. Empirically, we find that LD Scores are positively correlated with the strength of background selection (Appendix 3—figure 1), and that our simple model of background selection is capable of explaining much of the relationship between LD Score and allele frequency divergence that we observe in Figure 5 (Appendix 3—figures 1–3). Further, in the presence of BGS, bivariate regression of a stratified GWAS together with a measure of allele frequency differentiation can result in a positive slope, provided that the axis of stratification is correlated with the chosen measure of allele frequency divergence (Appendix 3).

50What conclusions should we draw from our LD Score regression analyses? The significant positive intercepts observed for LD Score regression of both GIANT and R15-sibs with [GBR-TSI] suggest that both datasets suffer from confounding due to stratification along a north-to-south axis within Europe. These observations are consistent with the evidence presented in Figure 3. A positive slope in such analyses was previously interpreted as evidence of positive selection either on height or a close genetic correlate (and presumed to be robust to stratification). However, BGS can, and empirically does, violate the assumptions of LD Score regression in a way that may generate a positive slope. We therefore interpret the positive slopes observed for the LD Score regression signals for GIANT and R15-sibs with SDS as likely resulting from a combination of stratification and BGS. A similar conclusion applies to the positive slope observed for R15-sibs [GBR-TSI]. It is unclear why stratification plus BGS should have elevated the slope for GIANT SDS, but not for GIANT [GBR-TSI]. This may suggest that the apparent SDS selection signal found in GIANT may be driven by an axis of variation that is not strongly correlated with [GBR-TSI]. We view this as an area worthy of further exploration.

51To summarize the key observations, we have reported the following:

52Multiple analyses based on GIANT and R15-sibs indicate strong signals of selection on height.

53However, the same signals of selection are absent or greatly attenuated in UK Biobank data. In some, but not all, analyses of frequency differentiation and SDS we still detect weakly significant signals of polygenic adaptation (Figures 1 and 2).

54The GIANT height GWAS is overall highly correlated with UKB-GB, but differs specifically by having an additional correlation with the main gradient of allele frequency variation across Europe, as modeled by frequency differences between GBR and TSI (Figure 3). LD Score analysis of [GBR-TSI] GIANT effect-size also suggests that GIANT is stratified along this axis (positive intercept in Figure 4D).

55Selection signals in the R15-sibs data are consistent with, and in some cases even stronger, than the corresponding signals in GIANT, but are inconsistent with analyses using UK Biobank data. While correctly implemented sib-based studies are designed to be impervious to population structure, [GBR-TSI] R15-sibs effect size also shows a highly positive intercept in the LD Score analysis presented in Figure 4E. As discussed below, the R15-sibs authors have recently identified a bug in the pipeline that generated this dataset. Analysis of corrected summary statistics does not show such a signal (Figure 4—figure supplement 1).

56LD Score analyses show a much stronger relationship between SDS and GIANT or R15-sibs than between SDS and UKB-GB. LD Score regression is generally considered to be robust to population structure (but see the discussion in ). However, the intensity of background selection increases with LD Score (Figure 5, Appendix 3), and this has likely inflated the LD Score-based signal of selection in GIANT and R15-sibs.

57In principle, it is possible that height in the UKB is confounded in a way that suppresses the signal of height adaptation. Instead, multiple lines of evidence strongly suggest that population-structure confounding in GIANT and R15-sibs is the main driver of the discrepancy with UKB-based analyses.

58The sib design used by Robinson et al offered a strong independent replication of the polygenic adaptation signal, which should have been impervious to population structure concerns. However, our analyses highlight multiple signs of stratification in this study. Robinson et al have now confirmed (as of November 2018) that effect sizes they released from their 2015 study were strongly affected by population stratification due to a bug. Furthermore, they have now stated that the effect sizes that they released publicly in 2016 were not the effect sizes used in their 2015 paper. As part of our own investigation, we have independently confirmed that sib-studies conducted using PLINK v1.90b5 are robust to environmental stratification (Appendix 4–figure 1). Our analyses using the newly rerun effect sizes released by show no consistent signal of selection (Figure 1—figure supplement 1, Figure 2—figure supplement 1. and Figure 4—figure supplement 1), in line with our UKB-based analyses.

59GIANT was conducted as a collaboration among a large number of research groups that provided summary statistics to the overall consortium. While the overall value of this pioneering dataset is not in question, it would not be surprising in retrospect if this GWAS were impacted by residual stratification along major axes of population structure.

60Lastly, we must conclude that the strong signal of LD Score genetic covariance between SDS and both GIANT and R15-sibs is largely spurious. This would imply that the LD Score regression slope is not robust to population structure confounding. Specifically, we demonstrated that background selection—through its correlation with LD Score—can potentially generate a spurious LD Score regression slope.

61Taken together, these observations lead us to conclude that what once appeared an ironclad example of population genetic evidence for polygenic adaptation now lacks any strong support. That said, there is still strong evidence that typical GWAS, including GIANT, do capture genuine signals of genotype-phenotype associations. For example, GWAS datasets regularly show strong functional enrichments of heritability within active chromatin from trait-relevant tissues, and the observation that top SNPs identified in GIANT tend to replicate in UKB-GB (Figure 3A), together with the high genetic correlations among all of the datasets (Table 2), suggest that the vast majority of the signal captured by GIANT is real.

62Nonetheless, we have shown that GIANT effect-size estimates contain a component arising from stratification along a major axis of European population structure (Figure 3B), and one would like to know the extent to which the conclusions from other analyses of GIANT, or other GWAS, may be affected. A complete investigation of this is beyond the scope of this study, and will depend on the nature of the analyses performed. The problem is likely most acute for the analysis of polygenic scores in samples drawn from heterogeneous ancestry. This is because while the bias in detection and effect sizes at any individual locus is small, the systematic nature of biases across many loci compound to significant errors at the level of polygenic scores. This error substantially inflates the proportion of the variance in polygenic scores that is among populations. Individual level prediction efforts therefore suffer dramatically from stratification bias, as even small differences in ancestry will be inadvertently translated into large differences in predicted phenotype. This seems likely to remain a difficult complication even within datasets such as the UK Biobank, though we suspect that meta-analyses such as GIANT, which collate summary statistics from many sources, may be particularly sensitive to structure confounding.

63These issues are apparent even within our UK Biobank results, where we see marked differences between results based on UKB-GB and UKB-Eur (Figure 2C,I vs. D, J and Figure 4—figure supplement 3). The study subjects in the two datasets were largely overlapping, and both were computed using widely-accepted structure-correction methods, suggesting that in the more demanding setting of broad European ancestry variation, the linear mixed model approach did not provide complete protection against stratification. This highlights a need for renewed exploration of the robustness of these methods, especially in the context of polygenic prediction.

64The study of polygenic scores across ancestry and environmental gradients offers a range of promises and pitfalls. Looking forward, we recommend that studies of polygenic adaptation should focus on datasets that minimize population structure (such as subsets of UKB), and where the investigators have access to full genotype data, including family data, so that they can explore sensitivity to different datasets and analysis pipelines.

65Figure 1 and Figure 2 display analyses based on six different GWAS. Two of these GWAS were newly calculated by us using UK Biobank data. Below, we describe the specifics of these two GWAS.

66To preform this GWAS, we used the following plink v. 2.0 with command line as follows:

67plink2 --memory 64G --threads 16 --linear

68--bpfile UKB imputed SNPs bp file

69--keep id list of individuals self-identified as ‘White British’

70--out output file

71--pheno standing height phenotype file (UKB phenotype 50.0.0)

72---covar covariates file

73The covariates file included only the sex, age and sequencing array for each individual id. We filtered all AG or CT SNPs–to prevent the possibility of strand errors. Finally, we excluded SNPs for which SDS was not calculated in.

74We used the estimated kinship coefficient and the proportion of SNPs for which the individuals share no allele (IBS0) provided by the UK Biobank, to call siblings as pairs withand IBS0 —following the conditions used by. We further filtered sibling pairs such that both individuals were ‘White British’, their reported sex matched their inferred sex, were not identified by the UK Biobank as ‘outliers’ based on heterozygosity and missing rate nor had an excessive number relatives in the data, and had height measurements. We standardized height values for each sex based on its mean and standard deviation (SD) values in the sample of 336,810 unrelated British ancestry individuals: mean 175.9 cm and SD 6.7 cm for males, and mean 162.7 cm and SD 6.2 cm for females. We also removed pairs if one of the siblings was more than 5 SD away from the mean. After applying all filters, 19,268 sibling pairs remained, equaling 35,524 individuals in 17,275 families. We performed an association analysis on 10,879,183 biallelic SNPs included in UKB-GB (converting dosages from imputation to genotype calls using no hard calling threshold), using plink v. 1.9 QFAM procedure with the following command:

75plink --bfile UKB hard-called SNPs file

76-out output file

77-qfam mperm = 100000

78The family relationships, as well as the phenotypic values, were encoded in plink FAM files.

79Individuals from the GBR and TSI populations from 1000G Phase 3 (N = 189) were assigned binary phenotype labels and a test was run using plink with a Hardy-Weinberg equilibrium cutoff of 1e-6 (--hwe 1e-6) and missing genotype rate of 0.05 (--geno 0.05), but otherwise with default parameters. Additionally, a firth adjusted logistic regression was run and produced qualitatively similar results (data not shown).

80Unrelated individuals, defined using estimates from KING, who self-identified as White British or White Irish in the UK Biobank were compared with distinct phenotype labels. Logistic regression was run on the genotyped SNP set using plink2 with a Hardy-Weinberg equilibrium cutoff of 1e-6 (--hwe 1e-6) and missing genotype rate of 0.05 (--geno 0.05).

811000 Genomes Phase 3 VCF files were downloaded from the 1000 Genomes website, and VCF files for the Human Origins dataset were downloaded from the ‘Affymetrix Human Origins fully public dataset’ link on the Reich lab website and subsequently imputed to full genomes using the Michigan imputation server. Because the Human Origins panel includes some 1000 Genomes populations, individual IDs were compared between the two datasets, and any duplicates were removed from the Human Origins dataset. Individuals were then clustered into populations based on groupings provided by each data resource, and allele frequencies were calculated using VCFtools version 0.1.15.

82In Figure 1, we employ two separate tests to assess the evidence that the distribution of polygenic scores among populations is driven in part by adaptive divergence. Both are based on a simple null model introduced by, which states that the distribution of polygenic scores under neutrality should be approximately multivariate normal. Here, we give a brief overview of the assumptions and calculations underlying the null model, before describing the two tests used in Figure 1. For a more complete treatment, see.

83Let be the vector of population allele frequencies at SNP, while is the effect size for SNP. Then, population level polygenic scores are given by

84Under neutrality, the distribution of polygenic scores among populations should be approximatelywherewhere is the mean of across populations. The matrix gives the population level co-ancestry among populations. Here, we calculate the matrix directly from the same set of SNPs used to calculate polygenic scores, which is a conservative procedure. Concretely, let

85Then, if is a matrix with the as columns, we have

86Now, based on this null model, we perform two separate neutrality tests. One is a general over-dispersion test (i.e. the ' test’ from ), for which the test statistic is

87For populations, this statistic is expected to have a distribution under the multivariate normal null model (Equation 2). An unusually large value of indicates that the neutral null model is a poor fit, and is therefore taken as evidence in favor of selection.

88We also apply a second, more specific test, to test for evidence of a correlation with a specific geographic axis that is unusually strong compared to the neutral expectation. For any vector, if has a multivariate normal distribution given by Equation 2, thenand thereforeunder the multivariate normal null. This fact can be used to test for an unexpectedly strong association between polygenic scores and a geographic axis by choosing to be the vector of latitudes or longitude across populations.

89To analyze tSDS as a function of GWAS p-value, we first divided SNPs into 5% minor allele frequency bins. We standardized SDS values—subtracted the mean and divided by the standard deviation—within each bin. While SDS values were already standardized in a similar manner by, we re-standardized SDS because the post-filtering composition of SNPs after in each GWAS was variable across GWAS. We then assigned tSDS values to each SNP by polarizing SDS to the tall allele. In other words, we setwhere derived is the derived allele in UK10K (by which SDS was polarized in ), and tall is the height increasing allele in the GWAS. We only used sites for which SDS values are available. Notably, this implicitly means that sites with minor allele frequency lower than 5% in UK10K were filtered out, due to the filtering used in.

90Figure 2 illustrates the correlation between tSDS and GWAS p-value (p-value for the strength of association with height). We assessed the significance of the correlation between the two while accounting for LD between SNPs. To do this, we used a blocked-jackknife approach to estimate the standard error of our Spearman correlation point estimate,. For each GWAS, SNPs were assigned to one of b = 200 contiguous blocks based on concatenated genomic coordinates. tSDS values should not be correlated across such large blocks. For each block, we computed the Spearman correlation in the ’th jackknife sample, —that is the Spearman correlation across all SNPs but the SNPs in block. We then estimated the standard error of the point Spearman estimate by, whereand.is the average of jackknife samples. Finally, we compute a p-value for the null hypothesis.by approximating as Normally distributed under the null with standard deviation, namely.

91Summary statistics for traits were filtered and allele flipped using munge_sumstats.py (a python program provided by ), with the default filters. All regressions were performed using the LD Score Regression package, using the LD Scores derived from the 489 unrelated European individuals in 1000 Genomes Phase III and a modified SNP set that excluded the HLA, LCT, and chromosome eight inversion loci.

92For genetic correlations of traits presented in Table 2, raw summary statistics were used. For other analyses, effect sizes of SNPs within each 1% minor allele frequency bin (as estimated by the 489 Europeans) were normalized to mean 0 and standard deviation 1, and those normalized statistics were used for downstream analyses. The standard two-step regression method from LD Score regression was used, with the default of 200 jackknife bins and a chi-square cutoff of 30, though results with UKB-GB were reasonably robust to a wide range of bin sizes and cutoffs.

93The GWAS generated from the UK Biobank for this paper have been uploaded to Dryad: https://doi.org/10.5061/dryad.mg1rr36

94The study also makes use of various publicly available GWAS datasets:

95The data from the GIANT consortium GWAS, conducted by is available at the GIANT consortium website: https://portals.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium_data_files

96The UK Biobank GWAS of individuals of ‘White British’ ancestry only (UKB-GB), conducted by, is available at:http://www.nealelab.is/blog/2017/7/19/rapid-gwas-of-thousands-of-phenotypes-for-337000-samples-in-the-uk-biobank

97The UK Biobank GWAS of individuals of broadly European ancestral (UKB-Eur), conducted by, is available at: https://data.broadinstitute.org/alkesgroup/UKBB/

98Sibling GWAS data from released in 2016 and now known to have been impacted by a computational bug, (R15-sibs,) as well as the newly rerun 2018 data (R15-sibs-updated,) are both available at: http://cnsgenomics.com/data.html

99Copies of these datasets are independently archived at: https://github.com/jjberg2/height-data.