Polygenic adaptation on height is overestimated due to uncorrected stratification in genome-wide association studies

1Most human complex traits are highly polygenic. For example, height has been estimated to be modulated by as much as 4% of human allelic variation. Polygenic traits are expected to evolve differently from monogenic ones, through slight but coordinated shifts in the frequencies of a large numbers of alleles, each with mostly small effect. In recent years, multiple methods have sought to detect selection on polygenic traits by evaluating whether shifts in the frequency of trait-associated alleles are correlated with the signed effects of the alleles estimated by genome-wide association studies (GWAS).

2Here we focus on a series of recent studies—some involving co-authors of the present manuscript—that have reported evidence of polygenic adaptation at alleles associated with height in Europeans. One set of studies observed that height-increasing alleles are systematically elevated in frequency in northern compared to southern European populations, a result that has subsequently been extended to ancient DNA. Another study using a very different methodology (singleton density scores, SDS) found that height-increasing alleles have systematically more recent coalescence times in the United Kingdom (UK) consistent with selection for increased height in the last few thousand years. In the present work, we assess polygenic adaptation on human height as a particular case of the effects that uncorrected population structure in GWAS can have on studies of complex traits.

3Most of these previous studies have been based on SNP associations and effect sizes (summary statistics) reported by the GIANT Consortium, which most recently combined 79 individual GWAS through meta-analysis, including a total of 253,288 individuals. Here, we show that the selection effects described in these studies are severely attenuated and in some cases no longer significant when using summary statistics derived from the UK Biobank, an independent and larger study that includes 336,474 genetically unrelated individuals who derive their recent ancestry almost entirely from the British Isles (identified as ‘white British ancestry’ by the UK Biobank) (Supplementary file 1). The UK Biobank analysis is based on a single cohort drawn from a relatively homogeneous population enabling better control of population stratification. Both datasets have high concordance even for low P value SNPs which do not reach genome-wide significance (Figure 1—figure supplement 1; genetic correlation between the two height studies is 0.94 [se = 0.0078]). Despite this concordance, we observe that small but systematic biases lead to the two datasets yielding qualitatively different conclusions with respect to signals of polygenic adaptation.

4To study population level differences among ancient and present-day European samples, we began by estimating ‘polygenic height scores’ as sums of allele frequencies at independent SNPs weighted by their effect sizes from GIANT. We used a set of different significance thresholds and strategies to correct for linkage disequilibrium as employed by previous studies, and replicated their signals for significant differences in genetic height across populations (Figure 1a, Figure 1—figure supplement 2). We then repeated the analysis using summary statistics from a GWAS for height in the UK Biobank restricting to individuals of British Isles ancestry (hereafter referred to as the ‘white British' (WB) subset) and correcting for population stratification based on the first ten principal components (UK Biobank [UKB]; also referred to as ‘UKB Neale’ in the supplementary figures). This analysis resulted in a dramatic attenuation of differences in polygenic height scores (Figure 1a, Figure 1—figure supplements 2–4). The differences between ancient European populations also greatly attenuated (Figure 1a, Figure 1—figure supplement 5). Strikingly, the ordering of the scores for populations also changed depending on which GWAS was used to estimate genetic height both within Europe (Figure 1a, Figure 1—figure supplements 2–5) and globally (Figure 1—figure supplement 6), consistent with reports from a recent simulation study. The height scores were qualitatively similar only when we restricted to independent genome-wide significant SNPs in GIANT and the UK Biobank (p<5×10−8) (Figure 1—figure supplement 2b). This replicates the originally reported significant north-south difference in the allele frequency of the height-increasing allele or in genetic height across Europe, as well as the finding of greater genetic height in ancient European steppe pastoralists than in ancient European farmers, although the signals are attenuated even here. Our observations suggest that tests of polygenic adaptation based on genome-wide significant SNPs are relatively consistent across different GWAS (Figure 1—figure supplement 2b) and that our concern is primarily directed towards the use of sub-significant SNPs in polygenic scores (Figure 1a, Figure 1—figure supplement 2a).

5Next, we assessed if an independent measure, the ‘singleton density score (SDS)', which uses a coalescent approach to infer adaptation within a population, is equally as susceptible to biases in GWAS. SDS can be combined with GWAS effect size estimates to infer polygenic adaptation on complex traits (generating a ‘tSDS score’ by aligning the SDS sign to the trait-increasing allele). A tSDS score larger than zero for height-increasing alleles implies that these alleles have been increasing in frequency in a population over time due to natural selection. We replicate the original finding that SDS scores of the height-increasing allele computed in the UK population (using the UK10K dataset) increase with stronger association of the alleles to height as inferred by GIANT across the entire P value spectrum (Spearman’s ρ = 0.078, p=1.55×10−65, Figure 1b). However, we observed that this signal of polygenic adaptation in the UK, measured using a Spearman correlation across all GWAS SNPs, disappeared when we used the UK Biobank height effect size estimates (ρ = 0.009, p=0.077, Figure 1b). These observations suggest that concerns about sub-significant SNPs should not only be directed towards population-level differences using polygenic scores but also to analyses of adaptation within a single population.

6We propose that the qualitative difference between the polygenic adaptation signals in GIANT and the UK Biobank is due to the cumulative effect of subtle biases in each of the SNPs estimated in GIANT. This bias can arise due to incomplete control of the population structure in GWAS. For example, if height were differentiated along a north-south axis because of differences in environment, any variant that is differentiated in frequency along the same axis would have an artificially large effect size estimated in the GWAS. Population structure is substantially less well controlled for in the GIANT study than in the UK Biobank study. This is both because the GIANT study population is more heterogeneous than that in the UK Biobank, and because population structure in the GIANT meta-analysis may not have been well controlled in some component cohorts due to their relatively small sizes (i.e., the ability to detect and correct population structure is dependent on sample size. The GIANT meta-analysis also found that such stratification effects worsen as SNPs below genome-wide significance are used to estimate height scores, consistent with our finding that the differences in genetic height among populations increase when including these SNPs.

7We obtained direct confirmation that population structure is more correlated with effect size estimates in GIANT than to those in the UK Biobank. Figure 2a shows that the effect sizes estimated in GIANT, in contrast to those in the UKB, are highly correlated with the SNP loadings of several principal components of population structure (PC loadings). We also find that the UK Biobank estimates including individuals of diverse ancestry and not correcting for population structure (UKB all no PCs) show the same stratification effects as GIANT (Figure 2—figure supplements 1–3). Further, in line with our intuition regarding the effects of residual stratification on GWAS effect size estimates, we find that alleles that are more common in the Great Britain population (1000 genomes GBR) than in the Tuscan population from Italy (1000 genomes TSI) tend to be preferentially estimated as height-increasing according to the GIANT study but not according to the UKB study (Figure 2c, Figure 2—figure supplements 2–3).

8We analyzed previously released family-based effect size estimates based on an approach of (NG2015 sibs). Surprisingly, we found that while these summary statistics produced significant polygenic adaptation signals, they were also correlated with PC loadings as well as with GBR-TSI allele frequency differences (Figure 2—figure supplements 1–3). This suggests that these estimates are also affected by population structure despite being computed within families and, therefore, in principle, robust to structure. Our own family-based estimates in the UK Biobank (UKB sibs all, UKB sibs WB) appear unconfounded and do not produce significant adaptation signals across the spectrum of associated SNPs (Figure 2—figure supplements 1–3). The residual structure in the original NG2015 sibs dataset is likely to reflect a technical artifact (personal communication from Peter Visscher, and note on their website). Berg and colleagues show that the updated NG2015 sibs summary statistics (posted in the public domain in November 2018 during the revision of this manuscript) do not show significant signals of polygenic adaptation using either polygenic score differences in Europe or the tSDS metric in the UK.

9We also note that the white British subset of the UKB data is not completely free of population stratification (as shown previously), although the magnitude of the potential confounding is much smaller than in the Continental European population (Figure 2—figure supplements 1–2). Interestingly, the north-south genetic cline in the UK tracks the height gradient in the opposite direction than in Continental Europe (Figure 2—figure supplements 2 and 4), and after correcting with principal components, we do not observe any evidence of residual stratification in comparison with the 1000 genomes data (Figure 2a,c). However, we cannot exclude the possibility of uncorrected population stratification, even in the UK Biobank, along axes not captured by the principal components of the 1000 genomes project data. For example, even for genome-wide significant SNPs (Figure 1—figure supplement 2b), polygenic scores for both modern and ancient individuals change when UKB summary statistics (WB ancestry controlling for 10 PCs) are used instead of GIANT. This shift, for example, for the ancient European hunter-gatherer polygenic score is troubling as different European populations are shown to have variable amounts of genetic ancestry from ancient ‘hunter-gatherer’ vs. ‘early farmer’ vs. ‘steppe ancestry’ populations, and could reflect residual stratification in the UKB GWAS not captured by the 1000 genomes PCs.

10We proceeded to investigate the effects of uncontrolled population stratification in GWAS discussed above on a coalescent approach such as tSDS that relies on singleton density. In principle, this approach is robust to the type of population stratification that affects the allele-frequency based tests. However, there is a north-south cline in singleton density in Europe due to lower genetic diversity in northern than in southern Europeans, leading to singleton density being lower in northern than in southern regions. As a consequence, SDS tends to be higher (corresponding to fewer singletons) in alleles more common in GBR than in TSI (Figure 2d). This cline in singleton density coincidentally parallels the phenotypic cline in height and the major axis of genome-wide genetic variation. Therefore, when we perform the tSDS test using GIANT, we find a higher SDS around the inferred height-increasing alleles, which tend, due to the uncontrolled population stratification in GIANT, to be at high frequency in northern Europe (Figure 2c). This effect does not appear when we use UK Biobank summary statistics because of the much lower level of population stratification and more modest variation in height. We find that SDS is not only correlated with GBR-TSI allele frequency differences, but with several principal component loadings across all SNPs (Figure 2b), and that these SDS-PC correlations often coincide with correlations between GIANT-estimated effect sizes and PC loadings (Figure 2a). We further find that the tSDS signal which is observed across the whole range of P values in some GWAS summary statistics can be mimicked by replacing SDS with GBR-TSI allele frequency differences (Figure 3a and c, Figure 3—figure supplements 1–4), suggesting that the tSDS signal at non-significant SNPs may be driven in part by residual population stratification.

11For polygenic adaptation within a population, a small but significant tSDS signal is observed in the UK when we restrict to genome-wide significant SNPs (p<5×10−8). This effect persists when using UK Biobank family-based estimates (UKB sibs WB) for genome-wide significant SNPs (Figure 3b), and is not driven by allele frequency differences between GBR and TSI (Figure 3d), suggesting an attenuated signal of polygenic adaptation in the UK that is driven by a much smaller number of SNPs than previously thought. Indeed, under most genetic architectures, a tSDS signal which is driven by natural selection is not expected to lead to an almost linear increase over the whole P value range in a well-powered GWAS. Instead, we would expect to see a greater difference between highly significant SNPs and non-significant SNPs, similar to the pattern observed in the UK Biobank (Figure 3a).

12For population-level differences in height, we assessed whether any remaining variation in height polygenic scores among populations is driven by polygenic adaptation by testing against a null model of genetic drift. We re-computed polygenic height scores in the POPRES dataset to increase power for this analysis as it has larger sample sizes of northern and southern Europeans than the 1000 Genomes project. We computed height scores using independent SNPs that are 1) genome-wide significant in the UK Biobank (‘gw-sig’, p<5×10−8) and 2) sub-significantly associated with height (‘sub-sig’, p<0.01) in different GWAS datasets. For each of these, we tested if population differences were significant due to an overall overdispersion (PQx), and if they were significant along a north-south cline (Plat) (Figure 4, Figure 4—figure supplements 1–2). Both gw-sig and sub-sig SNP-based scores computed using GIANT effect sizes showed significant overdispersion of height scores overall and along a latitude cline, consistent with previous results (Figure 4, Figure 4—figure supplements 1–2). However, the signal attenuated dramatically between sub-sig (Qx = 1100, PQx = 1×10−220) and gw-sig (Qx = 48, PQx = 2×10−4) height scores. In comparison, scores that were computed using the UK Biobank (UKB) effect sizes showed substantially attenuated differences using both sub-sig (Qx = 64, PQx = 5×10−7) and gw-sig (Qx = 33, PQx = 0.02) SNPs, and a smaller difference between the two scores. This suggests that the attenuation of the signal in GIANT is not only driven by a loss of power when using fewer gw-sig SNPs, but also reflects a decrease in stratification effects. The overdispersion signal disappeared entirely when the UK Biobank family based effect sizes were used (Figure 4, Figure 4—figure supplements 1–2). Moreover, Qx P values based on randomly ascertained SNPs and UK Biobank summary statistics are not uniformly distributed as would be expected if the theoretical null model is valid and if population structure is absent (Figure 4—figure supplement 3). The possibility of residual stratification effects even in the UK Biobank is also supported by a recent study. Therefore, we remain cautious about interpreting any residual signals as ‘real’ signals of polygenic adaptation.

13We have shown, by conducting a detailed analysis of human height, that estimates of population differences in polygenic scores are reduced when using the UK Biobank GWAS data relative to claims of previous studies that used GWAS meta-analyses such as GIANT. We find some evidence for population-level differences in genetic height, but it can only be robustly seen at highly significant SNPs, because any signal at less significant P values is dominated by the effect of residual population stratification. Even genome-wide significant SNPs in these analyses may be subtly affected by population structure, leading to continued overestimation of the effect. Thus, it is difficult to arrive at any quantitative conclusion regarding the proportion of the population differences that are due to statistical biases vs. population stratification of genetic height. Further, estimates of the number of independent genetic loci contributing to complex trait variation are sensitive to and likely confounded by residual population stratification.

14We conclude that while effect estimates are highly concordant between GIANT and the UK Biobank when measured individually (Supplementary file 5–7, Figure 1—figure supplement 1), they are also influenced by residual population stratification that can mislead comparisons of complex traits across populations and inferences about polygenic adaptation. Although these biases are subtle, in the context of tests for polygenic adaptation, which are driven by small systematic shifts in allele frequency, they can create highly significant artificial signals especially when SNPs that are not genome-wide significant are used to estimate genetic height. Our results do not question the reliability of the genome-wide significant associations discovered in the GIANT cohort. However, we urge caution in the interpretation of signals of polygenic adaptation or between-population differences that are based on large number of sub-significant SNPs–particularly when using effect sizes derived from meta-analysis of heterogeneous cohorts which may be unable to fully control for population structure.

15Our results have implications in other areas of human genetics research. For example, there is growing interest in polygenic scores that predict complex phenotypes from the aggregate effects of all allelic variants. The observation that individuals with extreme values of polygenic scores exhibit many-fold elevated risk of common diseases raises hopes for their potential clinical utility, and use for sociogenomics applications. It is already clear that polygenic scores derived from European populations do not translate across populations on a global scale. Our analysis further suggests that subtle population structure, especially in GWAS that are meta-analyses of independent cohorts, could be an additional source of error in polygenic scores and affect their applicability even within populations. We also note that other factors such as gene by environment interactions can be an alternative confounding factor for GWAS effect sizes and polygenic scores.

16We analyzed height using publicly available summary statistics that were obtained either by meta-analysis of multiple GWAS or by a GWAS performed on a single large population. We used results from the GIANT Consortium (N = 253,288) and a GWAS performed on individuals of the UK Biobank (‘UKB Neale’ or simply ‘UK Biobank (UKB)', N = 336,474) who derive their ancestry almost entirely from the British Isles (identified as ‘white British ancestry (WB)’ by the UK Biobank). The Neale lab’s GWAS uses a linear model with sex and 10 principal components as covariates. We also used an independent GWAS that included all UK Biobank European samples, allowing related individuals as well as population structure (‘UKB Loh’, N = 459,327). Loh et al.’s GWAS uses a BOLT-LMM Bayesian mixed model. Association signals from the three studies are generally correlated for SNPs that are genome-wide significant in GIANT (see ).

17We also used previously published family-based effect size estimates (‘NG2015 sibs’) as well as a number of test summary statistics on the UK Biobank that we generated to study the effects of population stratification. These are: ‘UKB Neale new’ (Similar to UKB Neale, with less stringent ancestry definition and 20 PCs calculated within sample), ‘UKB all no PCs’ (All UK Biobank samples included in the GWAS without correction by principal components), ‘UKB all 10 PCs’ (All UK Biobank samples included in the GWAS with correction by 10 principal components), ‘UK WB no PCs’ (Only ‘white British ancestry’ samples included in the GWAS without correction by principal components), ‘UKB WB 10 PCs’ (Only ‘white British ancestry’ samples included in the GWAS with correction by 10 principal components), ‘UKB sibs all’ (All UK Biobank siblings included in the GWAS), ‘UKB sibs WB’ (Only UK Biobank ‘white British ancestry’ siblings included in the GWAS) (Please see Supplementary file 1 for sample sizes and other details).

18We analyzed ancient and modern populations for which genotype data are publicly available. For ancient samples, we computed scores after dividing populations into three previously described broad ancestry labels (HG = Hunter Gatherer (n = 162 individuals), EF = Early Farmer (n = 485 individuals), and SP = Steppe Ancestry (n = 465 individuals)). For modern samples available through the 1000 genomes phase three release, we computed scores in two populations each from Northern Europe (GBR, CEU), Southern Europe (IBS, TSI), Africa (YRI, LWK), South Asia (PJL, BEB) and East Asia (CHB, JPT) (Figure 1a). In total, we analyzed 1112 ancient individuals, and 1005 modern individuals from 10 different populations in the 1000 genomes project (Supplementary file 2). We used the allele frequency differences between the GBR and TSI populations for a number of analyses to study population stratification (Figures 2–3). We also analyzed 19 European populations from the POPRES dataset with at least 10 samples per population (Figure 4—figure supplement 4).

19All ancient samples had ‘pseudo-haploid’ genotype calls at 1240k sites generated by selecting a single sequence randomly for each individual at each SNP. Thus, there is only a single allele from each individual at each site, but adjacent alleles might come from either of the two haplotypes of the individual. We also re-computed scores in present-day 1000 genomes individuals using only pseudo-haploid calls at 1240 k sites to allow for a fair comparison between ancient and modern samples (Figure 1—figure supplement 6).

20The polygenic scores, confidence intervals and test statistics (against the null model of genetic drift) were computed based on the methodology developed in references and. We computed the polygenic score (Z) for a trait in a population by taking the sum of allele frequencies in that population across all L sites associated with the trait, weighting each allele’s frequency by its effect on the trait.

21Al polygenic scores are plotted in centered standardized form,

22where,, and is the mean allele frequency across all populations analyzed. Source code repositories for the polygenic score analysis and computing scripts and source data for all the main figures have been made available at https://github.com/msohail88/polygenic_selection (; copy archived at https://github.com/elifesciences-publications/polygenic_selection) and https://github.com/uqrmaie1/sohail_maier_2019 (; copy archived at https://github.com/elifesciences-publications/sohail_maier_2019).

23Polygenic scores were computed using independent GWAS SNPs associated with height in three main ways: (1) The genome was divided into ~1700 non-overlapping linkage disequilibrium (LD) blocks (using the approximately independent linkage disequilibrium blocks in the EUR population computed in ), and the SNP with the lowest P value within each block was picked to give a set of ~1700 independent SNPs for each height GWAS used (all SNPs for which effect sizes are available were considered) similar to the analysis in. In (2) and (3), Plink’s clumping procedure was used to make independent ‘clumps’ of SNPs for each GWAS at different P value thresholds. This procedure selects SNPs below a given P value threshold as index SNPs to start clumps around, and then reduces all SNPs below a given P value threshold that are in LD with these index SNPs (above an r2 threshold, 0.1) and within a physical distance of them (1 Mb) into clumps with them. Clumps are preferentially formed around index SNPs with the lowest P value in a greedy manner. The index SNP from each clump is then picked for further polygenic score analyses. The algorithm is also greedy such that each SNP will only appear in one clump if at all. We clumped each GWAS to obtain (2) a set of independent sub-significant SNPs associated with height (p<0.01) similarly to, and (3) a set of independent genome-wide significant SNPs associated with height (p<5×10−8). The 1000 genomes phase three dataset was used as the reference panel for computing LD for the clumping procedure.

24The estimated effect sizes for these three sets of SNPs from each GWAS was used to compute scores. Only autosomal SNPs were used for all analyses to avoid creating artificial mean differences between populations with different numbers of males and females.

25The 95% credible intervals were constructed by assuming that the posterior of the underlying population allele frequency is independent across loci and populations and follows a beta distribution. We updated a Uniform prior distribution with allele counts from ancient and modern populations to obtain the posterior distribution at each locus in each population. We estimated the variance of the polygenic score using the variance of the posterior distribution at each locus, and computed the width of 95% credible intervals as for each population.

26The Qx test statistic measures the degree of overdispersion of the mean population polygenic score compared to a null model of genetic drift. It assumes that the vector of mean centered mean population polygenic score follows a multivariate normal distribution: Z ~ MVN(0, 2 VA F), where VA is the additive genetic variance of the ancestral population and F is a square matrix describing the population structure. This is equivalent to the univariate case of the test statistic used in. The latitude test statistic assumes that Y’Z ~ N(0, 2 VA Y’FY), where Y is a mean centered vector of latitudes for each population.

27The Singleton Density Score (SDS) method identifies signatures of recent positive selection based on a maximum likelihood estimate of the log-ratio of the mean tip-branch length of the derived vs. the ancestral allele at a given SNP. The tip-branch lengths are inferred from the average distance of each allele to the nearest singleton SNP across all individuals in a sequencing panel. When the sign of the SDS scores is aligned with the trait-increasing or trait-decreasing allele in the effect estimates of a GWAS, the Spearman correlation between the resulting tSDS scores and the GWAS P values has been proposed as an estimate of recent positive selection on polygenic traits.

28Here, we applied the tSDS method using pre-computed Singleton Density Scores for 4,451,435 autosomal SNPs obtained from 3195 individuals from the UK10K project for SNPs associated with height in GIANT and the UK biobank (Figure 1b) and in different summary statistics (Figure 3). After normalizing SDS scores in each 1% allele frequency bin to mean zero and unit variance, excluding SNPs from the MHC region on chromosome six and aligning the sign of the SDS scores to the height increasing alleles (resulting in tSDS scores), we computed the Spearman correlation coefficient between the tSDS score and the GWAS P value. The tSDS Spearman correlation standard errors and P values were computed using a block-jackknife approach, where each block of 1% of all SNPs ordered by genomic location was left out and the Spearman correlation coefficient was computed on the remaining SNPs. We also compared the tSDS score distributions for only genome-wide significant SNPs (Figure 3b).

29To compute SNP loadings of the principal components of population structure (PC loadings) in the 1000 genomes data (Figure 2), we first computed PC scores for each individual. We used SNPs that had matching alleles in 1000 genomes, GIANT and UK Biobank, that had minor allele frequency >5% in 1000 genomes, and that were not located in the MHC locus, the chromosome eight inversion region, or regions of long LD. After LD pruning to SNPs with r2 <0.2 relative to each other, PCA was performed in PLINK on the 187,160 remaining SNPs. In order to get SNP PC loadings for more SNPs than those that were used to compute PC scores, we performed linear regressions of the PC scores on the genotype allele count of each SNP (after controlling for sex) and used the resulting regression coefficients as the SNP PC loading estimates. The 1000 genomes phase one dataset was used to compute the PC loadings.