A Population Genetic Signal of Polygenic Adaptation

Summary

The process of adaptation is of fundamental importance in evolutionary biology. Within the last few decades, genotyping technologies and new statistical methods have given evolutionary biologists the ability to identify individual regions of the genome that are likely to have been important in this process. When adaptation occurs in traits that are underwritten by many genes, however, the genetic signals left behind are more diffuse, and no individual region of the genome is likely to show strong signatures of selection. Identifying this signature therefore requires a detailed annotation of sites associated with a particular phenotype. Here we develop and implement a suite of statistical methods to integrate this sort of annotation from genome wide association studies with allele frequency data from many populations, providing a powerful way to identify the signal of adaptation in polygenic traits. We apply our methods to test for the impact of selection on human height, skin pigmentation, body mass index, type 2 diabetes risk, and inflammatory bowel disease risk. We find relatively strong signals for height and skin pigmentation, moderate signals for inflammatory bowel disease, and comparatively little evidence for body mass index and type 2 diabetes risk.

Introduction

Population and quantitative genetics were in large part seeded by Fisher's insight [1] that the inheritance and evolution of quantitative characters could be explained by small contributions from many independent Mendelian loci [2]. While still theoretically aligned [3], these two fields have often been divergent in empirical practice. Evolutionary quantitative geneticists have historically focused either on mapping the genetic basis of relatively simple traits [4], or in the absence of any such knowledge, on understanding the evolutionary dynamics of phenotypes in response to selection over relatively short time-scales [5]. Population geneticists, on the other hand, have usually focused on understanding the subtle signals left in genetic data by selection over longer time scales [6]–[8], usually at the expense of a clear relationship between these patterns of genetic diversity and evolution at the phenotypic level.

Recent advances in population genetics have also allowed for the genome-wide identification of individual recent selective events either by identifying unusually large allele frequency differences among populations and environments or by detecting the effects of these events on linked diversity [9]. Such approaches are nonetheless limited because they rely on identifying individual loci that look unusual, and thus are only capable of identifying selection on traits where an individual allele has a large and/or sustained effect on fitness. When selection acts on a phenotype that is underwritten by a large number of loci, the response at any given locus is expected to be modest, and the signal instead manifests as a coordinated shift in allele frequency across many loci, with the phenotype increasing alleles all on average shifting in the same direction [10]–[14]. Because this signal is so weak at the level of the individual locus, it may be impossible to identify against the genome-wide background without a very specific annotation of which sites are the target of selection on a given trait [15], [16].

The advent of well-powered genome wide association studies with large sample sizes [17] has allowed for just this sort of annotation, enabling the mapping of many small effect alleles associated with phenotypic variation down to the scale of linkage disequilibrium in the population. The development and application of these methods in human populations has identified thousands of loci associated with a wide array of traits, largely confirming the polygenic view of phenotypic variation [18].

Although the field of human medical genetics has been the largest and most rapid to puruse such approaches, evolutionary geneticists studying non-human model organisms have also carried out GWAS for a wide array of fitness-associated traits, and the development of further resources is ongoing [19]–[21]. In human populations, the cumulative contribution of these loci to the additive variance so far only explain a fraction of the narrow sense heritability for a given trait (usually less than 15%), a phenomenon known as the missing heritability problem [22], [23]. Nonetheless, these GWAS hits represent a rich source of information about the loci underlying phenotypic variation.

Many investigators have begun to test whether the loci uncovered by these studies tend to be enriched for signals of selection, in the hopes of learning more about how adaptation has shaped phenotypic diversity and disease risk [24]–[27]. The tests applied are generally still predicated on the idea of identifying individual loci that look unusual, such that a positive signal of selection is only observed if some subset of the GWAS loci have experienced strong enough selection to make them individually distinguishable from the genomic background. As noted above, it is unlikely that such a signature will exist, or at least be easy to detect, if adaptation is truly polygenic, and thus many selective events will not be identified by this approach.

Here we develop and implement a general method based on simple quantitative and population genetic principals, using allele frequency data at GWAS loci to test for a signal of selection on the phenotypes they underwrite while accounting for the hierarchical structure among populations induced by shared history and genetic drift. Our work is most closely related to the recent work of Turchin et al [28], Fraser [29] and Corona et al [30], who look for co-ordinated shifts in allele frequencies of GWAS alleles for particular traits. Our approach constitutes an improvement over the methods implemented in these studies as it provides a high powered and theoretically grounded approach to investigate selection in an arbitrary number of populations with an arbitrary relatedness structure.

Using the set of GWAS effect size estimates and genome wide allele frequency data, we estimate the mean genetic value [31], [32] for the trait of interest in a diverse array of human populations. These genetic values may often be poor predictors of the actual phenotypes for reasons we address below and in the Discussion. We therefore make no strong claims about their ability to predict present day observed phenotypes. We instead focus on population genetic modeling of the joint distribution of genetic values, which provides a robust way of investigating how selection may have impacted the underlying loci.

We develop a framework to describe how genetic values covary across populations based on a flexible model of genetic drift and population history. In Figure 1 we show a schematic diagram of our approach to aid the reader. Using this null model, we implement simple test statistics based on transformations of the genetic values that remove this covariance among populations. We judge the significance of the departure from neutrality by comparing to a null distribution of test statistics constructed from well matched sets of control SNPs. Specifically, we test for local adaptation by asking whether the transformed genetic values show excessive correlations with environmental or geographic variables. We also develop and implement a less powerful but more general test, which asks whether the genetic values are over-dispersed among populations compared to our null model of drift. We show that this overdispersion test, which is closely related to [33], [34] and a series of approaches from the population genetics literature [35]–[39], gains considerable power to detect selection over single locus tests by looking for unexpected covariance among loci in the deviation they take from neutral expectations. Lastly, we develop an extension of our model that allows us to identify individual populations or groups of populations whose genetic values deviate from their neutral expectations given the values observed for related populations, and thus have likely been impacted by selection. While we develop these methods in the context of GWAS data, we also relate them to recent methodological developments in the quantitative genetics of measured phenotypes (as opposed to allele frequencies) [40], [41], highlighting the useful connection between these approaches. An implementation of the methods described here in the form of a collection of R scripts is available at https://github.com/jjberg2/PolygenicAdaptationCode.

Results

Estimating Genetic Values with GWAS Data

Consider a trait of interest where loci (e.g. biallelic SNPs) have been identified from a genome-wide association study. We arbitrarily label the phenotype increasing allele and the alternate allele at each locus. These loci have additive effect size estimates , where is the average increase in an individual's phenotype from replacing an allele with an allele at locus . We have allele frequency data for populations at our SNPs, and denote by the observed sample frequency of allele at the locus in the population. From these, we estimate the mean genetic value in the population as

and we take to be the vector containing the mean genetic values for all populations.

A Model of Genetic Value Drift

We are chiefly interested in developing a framework for testing the hypothesis that the joint distribution of is driven by neutral processes alone, with rejection of this hypothesis implying a role for selection. We first describe a general model for the expected joint distribution of estimated genetic values () across populations under neutrality, accounting for genetic drift and shared population history.

A simple approximation to a model of genetic drift is that the current frequency of an allele in a population is normally distributed around some ancestral frequency (). Under a Wright-Fisher model of genetic drift, the variance of this distribution is approximately , where is a property of the population shared by all loci, reflecting the compounded effect of many generations of binomially sampling [42]. Note also that for small values, is approximately equal to the inbreeding coefficient of the present day population relative to the defined ancestral population, and thus has an interpretation as the correlation between two randomly chosen alleles relative to the ancestral population [42].

We can expand this framework to describe the joint distribution of allele frequencies across an arbitrary number of populations for an arbitrary demographic history by assuming that the vector of allele frequencies in populations follows a multivariate normal distribution

where is an by positive definite matrix describing the correlation structure of allele frequencies across populations relative to the mean/ancestral frequency. Note again that for small values it is also approximately the matrix of inbreeding coefficients (on the diagonal) and kinship coefficients (on the off-diagonals) describing relatedness among populations [38], [43]. This flexible model was introduced, to our knowledge, by [44] (see [45] for a review), and has subsequently been used as a computationally tractable model for population history inference [42], [46], and as a null model for signals of selection [38], [39], [47], [48]. So long as the multivariate normal assumption of drift holds reasonably well, this framework can summarize arbitrary population histories, including tree-like structures with substantial gene flow between populations [46], or even those which lack any coherent tree-like component, such as isolation by distance models [49], [50].

Recall that our estimated genetic values are merely a sum of sample allele frequencies weighted by effect size. If the underlying allele frequencies are well explained by the multivariate normal model described above, then the distribution of is a weighted sum of multivariate normals, such that this distribution is itself multivariate normal

where and are respectively the expected genetic value and additive genetic variance of the ancestral (global) population. The covariance matrix describing the distribution of therefore differs from that describing the distribution of frequencies at individual loci only by a scaling factor that can be interpreted as two times the contribution of the associated loci to the additive genetic variance present in a hypothetical population with allele frequencies equal to the grand mean of the sampled populations.

The assumption that the drift of allele frequencies around their shared mean is normally distributed (2) may be problematic if there is substantial drift. However, even if that is the case, the estimated genetic values may still be assumed to follow a multivariate normal distribution by appealing to the central limit theorem, as each estimated genetic value is a sum over many loci. We show in the Results that this assumption often holds in practice.

It is useful here to note that the relationship between the model for drift at the individual locus level, and at the genetic value level, gives an insight into where most of the information and statistical power for our methods will come from. Each locus adds a contribution to the vector of deviations of the genetic values from the global mean. If the allele frequencies are unaffected by selection then the frequency deviation of an allele at locus in population will be uncorrelated in magnitude or sign with both the effect at locus and the allele frequency deviation taken by other unlinked loci. Thus the expected departure of the genetic value of a population from the mean is zero, and the noise around this should be well described by our multivariate normal model.

The tests described below will give positive results when these observations are violated. The effect of selection is to induce a non-independence of the allele frequency deviation () across loci, determined by the sign and magnitude of the effect sizes [10]–[14] and as we demonstrate below, all of our methods rely principally on identifying this non-independence. This observations has important considerations for the false positive profile of our methods. Specifically, false positives will arise only if the GWAS ascertainment procedure induces a correlation between the estimated effect size of an allele () and the deviation that this allele takes across populations . This should not be the case if the GWAS is performed in a single population which is well mixed compared to the populations considered in the test. False positives can occur when a GWAS is performed in a structured population and fails to account for the fact that the phenotype of interest is correlated with ancestry in this population. We address this case in greater depth in the Discussion.

These observation also allows us to exclude certain sources of statistical error as a cause of false positives. For example, simple error in the estimation of , or failing to include all loci affecting a trait cannot cause false positives, because this error has no systematic effect on across loci. Similarly, if the trait of interest truly is neutral, variation in the true effects of an allele across populations or over time or space (which can arise from epistatic interactions among loci, or from gene by environment interactions) will not drive false positives, again because no systematic trends in population deviations will arise. This sort of heterogeneity can, however, reduce statistical power, as well as make straightforward interpretation of positive results difficult, points which we address further below.

Fitting the Model and Standardizing the Estimated Genetic Values

As described above, we obtain the vector by summing allele frequencies across loci while weighting by effect size. We do not get to observe the ancestral genetic value of the sample , so we assume that this is simply equal to the mean genetic value across populations . This assumption costs us a degree of freedom, and so we must work with a vector , which is the vector of estimated genetic values for the first populations, centered at the mean of the (see Methods for details). Note that this procedure will be the norm for the rest of this paper, and thus we will always work with vectors of length that are obtained by subtracting the mean of the vector and dropping the last component. The information about the dropped population is retained in the mean of the length vectors, and thus the choice of which population to drop is arbitrary and does not affect the inference.

To estimate the null covariance structure of the populations we sample a large number K random unlinked SNPs. In our procedure, the SNPs are sampled so as to match certain properties of the GWAS SNPs (the specific matching procedure is described in more depth below and in the Methods section). Setting to be the mean sample allele frequency across populations at the SNP, we standardize the sample allele frequency in population as . We then calculate the sample covariance matrix () of these standardized frequencies, accounting for the rank of the matrix (see Methods). We estimate the scaling factor of this matrix as

We now have an estimated genetic value for each population, and a simple null model describing their expected covariance due to shared population history. Under this multivariate normal framework, we can transform the vector of mean centered genetic values () so as to remove this covariance. First, we note that the Cholesky decomposition of the matrix is

where is a lower triangular matrix, and is its transpose. Informally, this can be thought of as taking the square root of , and so can loosely be thought of as analogous to the standard deviation matrix.

Using this matrix we can transform our estimated genetic values as:

If then , where is the identity matrix. Therefore, under the assumptions of our model, these standardized genetic values should be independent and identically distributed random variates [39].

It is worth spending a moment to consider what this transformation has done to the allele frequencies at the loci underlying the estimated genetic values. As our original genetic values are written as a weighted sum of allele frequencies, our transformed genetic values can be written as a weighted sum of transformed allele frequencies (which have passed through the same transform). We can write

and so we can define the vector of transformed allele frequencies at locus to be

This set of transformed frequencies exist within a set of transformed populations, which by definition have zero covariance with one another under the null, and are related by a star-like population tree with branches of equal length.

As such, we can proceed with simple, straightforward and familiar statistical approaches to test for the impact of spatially varying selection on the estimated genetic values. Below we describe three simple methods for identifying the signature of polygenic adaptation, which arise naturally from this observation.

Environmental Correlations

We first test if the genetic values are unusually correlated with an environmental variable across populations compared to our null model. A significant correlation is consistent with the hypothesis that the populations are locally adapted, via the phenotype, to local conditions that are correlated with the environmental variable. However, the link from correlation to causation must be supported by alternate forms of evidence, and in the lack of such evidence, a positive result from our environmental correlation tests may be consistent with many explanations.

Assume we have a vector , containing measurements of a specific environmental variable of interest in each of the populations. We mean-center this vector and put it through a transform identical to that which we applied to the estimated genetic values in (7). This gives us a vector , which is in the same frame of reference as the transformed genetic values.

There are many possible models to describe the relationship between a trait of interest and a particular environmental variable that may act as a selective agent. We first consider a simple linear model, where we model the distribution of transformed genetic values () as a linear effect of the transformed environmental variables ()

where under our null is a set of normal, independent and identically distributed random variates (i.e. residuals), and can simply be estimated as . We can also calculate the associated squared Pearson correlation coefficient () as a measure of the fraction of variance explained by our variable of choice, as well as the non-parametric Spearman's rank correlation , which is robust to outliers that can mislead the linear model. We note that we could equivalently pose this linear model as a mixed effects model, with a random effect covariance matrix . However, as we know both and , we would not have to estimate any of the random effect parameters, reducing it to a fixed effect model as in (9) [51].

In the Methods (section “The Linear Model at the Individual Locus Level”) we show that the linear environmental model applied to our transformed genetic values has a natural interpretation in terms of the underlying individual loci. Therefore, exploring the environmental correlations of estimated genetic values nicely summarizes information in a sensible way at the underlying loci identified by the GWAS.

In order to assess the significance of these measures, we implement an empirical null hypothesis testing framework, using , , and as test statistics. We sample many sets of SNPs randomly from the genome, again applying a matching procedure discussed below and in the Methods. With each set of SNPs we construct a vector , which represents a single draw from the genome-wide null distribution for a trait with the given ascertainment profile. We then perform an identical set of transformations and analyses on each , thus obtaining an empirical genome-wide null distribution for all test statistics.

Excess Variance Test

As an alternative to testing the hypothesis of an effect by a specific environmental variable, one might simply test whether the estimated genetic values exhibit more variance among populations than expected due to drift. Here we develop a simple test of this hypothesis.

As is composed of independent, identically distributed standard normal random variables, a natural choice of test statistic is given by

This statistic represents a standardized measure of the among population variance in estimated genetic values that is not explained by drift and shared history. It is also worth noting that by comparing the rightmost form in (10) to the multivariate normal likelihood function, we find that is proportional to the negative log likelihood of the estimated genetic values under the neutral null model, and is thus the natural measurement of the model's ability to explain their distribution. Multivariate normal theory predicts that this statistic should follow a distribution with degrees of freedom under the null hypothesis. Nonetheless, we use a similar approach to that described for the linear model, generating the empirical null distribution by resampling SNPs genome-wide. As discussed below, we find that in practice the empirical null distribution tends to be very closely matched by the theoretically predicted distribution.

Values of this statistic that are in the upper tail correspond to an excess of variance among populations. This excess of variance is consistent with the differential action of natural selection on the phenotype among populations (e.g. due to local adaptation). Values in the lower tail correspond a paucity of variance, and thus potentially to widespread stabilizing selection, with many populations selected for the same optimum. In this paper we report upper tail p-values from the empirical null distribution of both for our power simulations and empirical results. A two tailed test would be appropriate in cases where stabilizing selection is also of interest, however such signals are likely to be difficult to spot with GWAS data because the we are missing the large effect, low frequency alleles most likely to reveal a signal of stabilizing selection.

The relationship of to previous tests

Our statistic is closely related to , the phenotypic analog of , which measures the fraction of the genetic variance that is among populations relative to the total genetic variance [33], [34], [52]. is typically estimated in traditional local adaptation studies via careful measurement of phenotypes from related individuals in multiple populations in a common garden setting. If the loci underlying the trait act in a purely additive manner and are experiencing only neutral genetic drift, then [53], [54].

If both quantities are well estimated, and we also assume that there is no hierarchical structure among the populations, then is known to have a distribution under a wide range of models [55]–[57]. This statistic is thus a natural phenotypic extension of Lewontin and Krakauer's based-test (LK test) [35].

To see the close correspondence between and , consider the case of a starlike population tree with branches of equal length (i.e. and ). Under this demographic model, we have

where is an estimated value for obtained from our estimated genetic values. This relationship between and breaks down when some pairs of populations do not have zero covariance in allele frequencies under the null, in which case the distribution of the LK test also breaks down [36], [37]. Bonhomme and colleagues [38] recently proposed an extension to the LK test that accounts for a population tree, thereby recovering the distribution (see also [39], which relaxes the tree-like assumption), and our statistic is a natural extension of this enhanced statistic to the problem of detecting coordinated selection at multiple loci. This test is also nearly identical to that developed independently by Ovaskainen and colleagues for application to direct phenotypic measurements [40].

Writing in terms of allele frequencies

Given that our estimated genetic values are simple linear sums of allele frequencies, it is natural to ask how can be written in terms of these frequencies. Again, restricting ourselves to the case where is diagonal, (i.e. and ), we can express as

which can be rewritten as

The numerator of the first term inside the parentheses is the weighted sum of the variance among populations over all GWAS loci, scaled by the contribution of those loci to the additive genetic variance in the total population. As such this first term is similar to calculated for our GWAS loci, but instead of just averaging the among population and total variances equally across loci in the numerator and denominator, these quantities are weighted by the squared effect size at each locus. This weighting nicely captures the relative importance of different loci to the trait of interest.

The second term in (13) is less familiar; the numerator is the weighted sum of the covariance of allele frequencies between all pairs of GWAS loci, and the denominator is again the contribution of those loci to the additive genetic variance in the total population. This term is thus a measure of the correlation among loci in the deviation they take from the ancestral value, or the across population component of linkage disequilibrium. For a more in depth discussion of this relationship in the context of , see [10]–[14].

As noted above (8), when is non-diagonal, our transformed genetic values can be written as a weighted sum of transformed allele frequencies. Consequently, we can obtain a similar expression to (13) when population structure exists, but now expressed in terms of the covariance of a set of transformed allele frequencies in populations that have no covariance with each other under the null hypothesis. Specifically, when the covariance is non-diagonal we can write:

We refer to the first term in this decomposition as the standardized -like component and the second term as the standardized LD-like component. Under the neutral null hypothesis, the expectation of the second term is equal to zero, as drifting loci are equally likely to covary in either direction. With differential selection among populations, however, we expect loci underlying a trait not only to vary more than we would expect under a neutral model, but also to covary in a consistent way across populations. Models of local adaptation predict that it is this covariance among alleles that is primarily responsible for differentiation at the phenotypic level [10]–[14], and we therefore expect the statistic to offer considerably increased power as compared to measuring average or identifying outliers. We use simulations to demonstrate this fact below, and also demonstrate the perhaps surprising result that for a broad parameter range the standardized LD-like component exhibits almost no loss of power when used as a test statistic.

Identifying Outlier Populations

Having detected a putative signal of selection for a given trait, one may wish to identify individual regions and populations which contribute to the signal. Here we rely on our multivariate normal model of relatedness among populations, along with well understood methods for generating conditional multivariate normal distributions, in order to investigate specific hypotheses about individual populations or groups of populations. Using standard results from multivariate normal theory, we can generate the expected joint conditional distribution of genetic values for an arbitrary set of populations given the observed genetic values in some other set of populations. These conditional distributions allow for a convenient way to ask whether the estimated genetic values observed in certain populations or groups of populations differ significantly from the values we would expect them to take under the neutral model given the values observed in related populations.

Specifically, we exclude a population or set of populations, and then calculate the expected mean and variance of genetic values in these excluded populations given the values observed in the remaining populations, and the covariance matrix relating them. Using this conditional mean and variance, we calculate a Z-score to describe how well fit the estimated genetic values of the excluded populations are by our model of drift, conditional on the values in the remaining populations. In simple terms, the observation of an extreme Z-score for a particular population or group of populations may be seen as evidence that that group has experienced directional selection on the trait of interest (or a correlated one) that was not experienced by the related populations on which we condition the analyses. The approach cannot uniquely determine the target of selection, however. For example, conditioning on populations that have themselves been influenced by directional selection may lead to large Z-scores for the population being tested, even if that population has been evolving neutrally. We refer the reader to the Methods section for a mathematical explication of these approaches.

Datasets

We conducted power simulations and an empirical application of our methods based on the Human Genome Diversity Panel (HGDP) population genomic dataset [58], and a number of GWAS SNP sets. To ensure that we made the fullest possible use of the information in the HGDP data, we took advantage of a genome wide allele frequency dataset of 3 million SNPs imputed from the Phase II HapMap into the 52 populations of the HGDP. These SNPs were imputed as part of the HGDP phasing procedure in [59]; see our Methods section for a recap of the details. We applied our method to test for signals of selection in six human GWAS datasets identifying SNPs associated with height, skin pigmentation, body mass index (BMI), type 2 diabetes (T2D), Crohn's Disease (CD) and Ulcerative Colitis (UC).

Choosing null SNPs

Various components of our procedure involve sampling random sets of SNPs from across the genome. While we control for biases in our test statistics introduced by population structure through our matrix, we are also concerned that subtle ascertainment effects of the GWAS process could lead to biased test statistics, even under neutral conditions. We control for this possibility by sampling null SNPs so as to match the joint distribution of certain properties of the ascertained GWAS SNPs. Specifically, we chose our random SNPs to match the GWAS SNPs in each study in terms of their minor allele frequency (MAF) in the ascertainment population and the imputation status of the allele in our population genomic dataset (i.e. whether the allele was imputed or present in the original HGDP genotyping panel). In addition, we were concerned that GWAS SNPs might be preferentially found close to genes and in low recombination regions, the latter due to better tagging, and as such may be subject to a high rate of drift due to background selection, leading to higher levels of differentiation at these sites [60]. Therefore, in addition to MAF and imputation status, we also matched our random SNPs to an estimate of the background selection environment experienced by the GWAS SNPs, as measured by B value [61], which is a function of both the density of functional sites and recombination rate calibrated to match the reduction in genetic diversity due to background selection. We detail the specifics of the binning scheme for matching the discretized distributions of GWAS and random SNPs in the Methods.

Power Simulations

To assess the power of our methods in comparison to other possible approaches, we conducted a series of power simulations. There are two possible approaches to simulate the effect of selection on large scale allele frequency data of the type for which our methods are designed. The first is to simulate under some approximate model of the evolutionary history (e.g. full forward simulation under the Wright-Fisher model with selection). The second is to perturb real data in such a way that approximates the effect of selection. We choose to pursue the latter, both because it is more computationally tractable, and because it allows us to compare the power of our different approaches for populations with evolutionary histories of the same complexity as the real data we analyze. Each of our simulations will thus consist of sampling 1000 sets of SNPs matched to the height dataset (in much the same way we sample SNPs to construct the null distributions of our test statistics), and then adding slight shifts in frequency in various ways to mimic the effect of selection.

Below we first describe the set of alternative statistics to which we compare our methods. We then describe the manner in which we add perturbations to mimic selection, and lastly describe a number of variations on this theme which we pursued in order to better demonstrate how the power of our statistics changes as we vary parameters of the trait of interest, evolutionary process, or the ascertainment.

Statistics tested

For our first set of simulation experiments we compared two of our statistics, ( and ) against their naive counterparts, which are not adjusted for population structure (naive and ). We also include the adjusted -like and LD-like components of as their behavior over certain parameter ranges is particularly illuminating. For , , and it's components, we count a given simulation as producing a positive result if the statistic lies in the upper 5% tail of the null distribution, whereas for the environmental correlation statistics ( and naive ) we use a two-tailed 5% test. We also compared our tests to a single locus enrichment test, where we tested for an enrichment in the number of SNPs that individually show a correlation with the environmental variable. We considered this test to produce a positive result if the number of individual loci in the 5% tail of the null distribution for individual locus was itself in the 5% tail using a binomial test. We do not include our alternative linear model statistics and in these plots for the sake of figure legibility, but they generally had very similar power to that of . While slightly more powerful versions of the enrichment test that better account for sampling noise are available [47], note that our tests could be extended similarly as well, so the comparison is fair.

Simulating selection

We base our initial power simulations on empirical data altered to have an increasing effect of directional selection along a latitudinal gradient. In order to mimic the effect of selection, we generate a new set of allele frequencies () by taking the original frequency () and adding a small shift according to

where is the effect size assigned to locus , and is the mean centered absolute latitude of the population. We use 1000 simulations at to form null distribution for each of our test statistics, and from this established the significance level. We then increment and give the power of each statistic as the fraction of simulations whose test statistic falls beyond this cutoff. While this approach to simulating selection is obviously naive to the way selection actually operates, it captures many of the important effects on the loci underlying a given trait. Namely, loci will have greater shifts if they experience extreme environments, have large effects on the phenotype, or are at intermediate frequencies. Because we add these shifts to allele frequencies sampled from real, putatively neutral loci, the effect of drift on their joint distribution is already present, and thus does not need to be simulated. The results of these simulations are shown in Figure 2A.

Our population structure adjusted statistics clearly outperform tests that do not account for structure, as well as the single locus outlier based test. Particularly noteworthy is the fact that the power of a test relying on and that using only the LD-like component are essentially identical over the entire range of simulation, while the -like component achieves only about power by the point at which the former statistics have reached 100%. This reinforces the observation from previous studies of that for polygenic traits, nearly all of the differentiation at the trait level arises as a consequence of across population covariance among the underlying loci, and not as a result of substantial differentiation at the loci themselves [14]. While our environment-genetic value correlation tests considerably outperform , this is somewhat artificial as it assumes that we know the environmental variable responsible for our allele frequency shift. In reality, the power of the environmental variable test will depend on the investigator's ability to accurately identify the causal variable (or one closely correlated with it) in the particular system under study, and thus in some cases may have have higher power in practice. Panels A and B from Figure 2 with SNPs matched to each of the other traits we investigate can be found in Figures S1–S5.

Pleiotropy and correlated selection

We next considered the fact that many of the loci uncovered by GWAS are may be relatively pleiotropic, and thus may simultaneously respond to selection on multiple different traits. To explore how our methods perform in the presence of undetected pleiotropy, we consider the realization that from the perspective of allele frequency change there is only one effect that matters, and that is the effect on fitness. We therefore chose a simple and general approach to capture a flavor of this situation. We simulate the effect of selection as above (15), but give each locus an effect on fitness () that may be only partially correlated with the observed effect sizes for the trait of interest (with the unaccounted for effect on fitness coming via pleiotropic relationships to any number of unaccounted for phenotypes). For simplicity we assume that and have a bivariate normal distribution around zero with equal variance and correlation parameter . We then simulate from its conditional distribution given (i.e. ). For each SNP in (15) we replaced by its effect on the unobserved phenotype, but then perform our tests using the measured for the trait of interest. Here can be thought of as the genetic correlation between our phenotype and fitness if this simple multivariate form held true for all of the loci contributing to the trait. The extremes of and respectively represent the cases where selection acts only on the focal trait and that were all the underlying loci are affected by selection, but not due to their relationship with the focal trait. These simulations can also informally be seen as modeling the case where the GWAS estimated effect sizes are imperfectly correlated with the true effect sizes that selection sees, for example due to measurement error in the GWAS.

In Figure 2B we hold the value of constant at 0.14 and vary the genetic correlation from one down to zero. Predictably, our GWAS genetic value based statistics lose power as the the focal trait becomes less correlated with fitness but do retain reasonable power out to quite low genetic correlations (e.g. our out performs the single locus metrics until ). In contrast, counting the number of SNPs that are significantly correlated with a given environmental variable remains equally powerful across all genetic correlations. This is because the single locus environmental correlation tests treat each locus separately with no regards to whether there is agreement across alleles with the same direction of effect size. This may be a desirable property of the environmental outliers enrichment approach, as it does not rely on a close relationship between the effect sizes and the way that selection acts on the loci. On the other hand, this is also problematic, as such tests may often be detecting selection on only very weakly pleiotropically related phenotypes. Our approaches, however, are more suited to determining whether the genetic basis of a trait of interest, or one that is genetically correlated it, has been affected by differentiating selection.

Ascertainment and genetic architecture

We next investigated the relationship between statistical power, the number of loci associated with the trait, and the amount of variance explained by those loci. Our simulations were motivated by the fact that the number of loci identified by a given GWAS, and the fraction of variance explained by those loci, will depend on both the design of the study (e.g. sample size) and the genetic architecture of the trait. To illustrate the impact these factors have on the power of our methods, we performed two experiments in which we again held constant at 0.14. In the first, for each of the 1000 sets of 161 loci chosen above to mimic the height data ascertainment, we randomly sampled loci, without regard to effect sizes, and recalculated the null distribution and power for these reduced sets, allowing to range from 2 to 161. The results of these simulations are shown in Figure 2C. This corresponds to imagining that fewer loci had been ascertained by the initial GWAS, and estimating the power our methods would have with this reduced set of loci. As we down sample our loci without regard to effect sizes, the horizontal axis of Figure 2C is proportional to the phenotypic variance explained, e.g. the simulations in which only 80 loci are subsampled correspond to having a dataset which explains only 50% of the variance explained in those for which all 161 were used.

The second experiment is nearly identical to the first, except that before adding an effect of selection to the subsampled loci, we linearly rescale the effect sizes such that is held constant at the value calculated for the full set of 161 loci. The results of these simulations are shown in Figure 2D. These simulations correspond to imagining that we have explained an equivalent amount of phenotypic variance, but the number of loci over which this variation is partitioned varies.

Our results (Figure 2C and 2D) demonstrate that even if only a small number of loci associated with the phenotype have been identified, our tests offer higher power than single locus-based tests. Moreover, for statistics that appropriately deal with both covariance among loci and among populations ( and ), power is generally a constant function of variance explained by the underlying loci, regardless of the number of loci over which it is partitioned. Notably, most the power of comes from the LD-like component, especially when the number of loci is large. Statistics that rely on an average of single locus metrics (the -like component of ), and those that rely on outliers ( enrichment) all lose power as the the variance explained is partitioned over more loci, as the effect of selection at each locus is weaker. Somewhat surprisingly, the versions of our tests that fail to adequately control for population structure (naive and ) also lose power as the phenotypic variance is spread among more loci. We believe this reflects the fact that they are being systematically mislead by LD among SNPs due to population structure, a problem which is compounded as more loci are included in the test. Overall these results suggest that accounting for population structure and using the LD between like effect alleles is key to detecting selection on polygenic phenotypes.

Localizing signatures of selection

Lastly, we investigated the power of our conditional Z-scores to identify signals of selection that are specific to particular populations or geographic regions, and contrast this with the power of the global statistic to detect the same signal. We again perform two experiments. In the first, we choose a single population whose allele frequencies to perturb, and leave all other populations unchanged. In other words, an effect of selection is mimicked according to (15), but with set equal to one for a single population, and zero for all others. We then increment to see how power changes as the effect of selection becomes more pronounced. In Figure 2E we display the results of these simulations for five populations chosen to capture the range of power profiles for the populations we consider in our empirical applications. In the last experiment, we chose a group of populations to which to apply the allele frequency shift, again consistent with (15), but now with set equal to 1 for all populations in an entire region, and zero elsewhere. In Figure 2F, we show the results of these simulations, with each of the seven geographic/genetic clusters identified by Rosenberg et al (2002) [62], chosen in turn as the affected region.

These simulations demonstrate that the conditional Z test can detect subtler frequency shifts than the global test, provided one knows which population(s) to test a priori. They also show how unusual frequency patterns indicative of selection are easier to detect in populations for which the dataset contains closely related populations that are unaffected (e.g. compare the Han and Italian to the San and Karitiana at the individual population level, or Europe, the Middle East and Central Asia to Africa, America, and Oceania at the regional level). Lastly, note that the horizontal axes in Figure 2E and 2F are equivalent in the sense that for a given value of , alleles in (say) the Italian population have been shifted by the same amount in the Italian specific simulations in Figure 2E as in the Europe-wide simulation in Figure 2F, indicating that the HGDP dataset, power is similar in efforts to detect local, population specific events, as well as broader scale, regional level events.

Empirical Applications

We estimated genetic values for each of six traits from the subset of GWAS SNPs that were present in the HGDP dataset, as described above. We discuss the analysis of each dataset in detail below, and address general points first. For each dataset, we constructed the covariance matrix from a sample of approximately appropriately matched SNPs, and the null distributions of our test statistics from a sample of sets of null genetic values, which were also constructed according to a similar matching procedure (as described in the Methods).

In an effort to be descriptive and unbiased in our decisions about which environmental variables to test, we tested each trait for an effect of the major climate variables considered by Hancock et al (2008) [63] in their analysis of adaptation to climate at the level of individual SNPs. We followed their general procedure by running principal components (PC) analysis for both seasons on a matrix containing six major climate variables, as well as latitude and longitude (following Hancock et al's rationale that these two geographic variables may capture certain elements of the long term climatic environment experienced by human populations). The percent of the variance explained by these PCs and their weighting (eigenvectors) of the different environmental variables are given in Table 1. We view these analyses largely as a descriptive data exploration enterprise across a relatively small number of phenotypes and distinct environmental variables, and do not impose a multiple testing penalty against our significance measures. A multiple testing penalization or false discovery rate approach may be needed when testing a large number phenotypes and/or environmental variables.

We also applied our test to identify traits whose underlying loci showed consistent patterns of unusual differentiation across populations, with results reported in Table 2. In Figure 3 we show for each GWAS set the observed value of and its empirical null distribution calculated using SNPs matched to the GWAS loci as described above. We also plot the expected null distribution of the statistic (). The expected null distribution closely matches the empirical distribution in all cases, suggesting that our multivariate normal framework provides a good null model for the data (although we will use the empirical null distribution to obtain measures of statistical significance).

For each GWAS SNP set we also separate our statistic into its -like and LD-like terms, as described in (14). In Figure 4 we plot the null distributions of these two components for the height dataset as histograms, with the observed value marked by red arrows (Figures S6–S10 give these plots for the other five traits we examined). In accordance with the expectation from our power simulations, the signal of selection on height is driven entirely by covariance among loci in their deviations from neutrality, and not by the deviations themselves being unusually large.

Lastly, we pursue a number or regionally restricted analyses. For each trait and for each of the seven geographic/genetic clusters described by Rosenberg et al (2002) [62], we compute a region specific statistic to get a sense for the extent to which global signals we detect can be explained by variation among populations with these regions, and to highlight particular populations and traits which may merit further examination as more association data becomes available. The results are reported in Table 3. We also apply our conditional Z-score approach at two levels of population structure: first at the level of Rosenberg's geographic/genetic clusters, testing each cluster in turn for how differentiated it is from the rest of the world, and second at the level of individual populations. The regional level Z-scores are useful for identifying signals of selection acting over broad regional scale or on deeper evolutionary timescales, while the population specific Z-scores are useful for identifying very recent selection that has only impacted a single population. We generally employ these regional statistics as a heuristic tool to localize signatures of selection uncovered in global analyses, or in cases where there is no globally interesting signal, to highlight populations or regions which may merit further examination as more association data becomes available. The result of these analyses are depicted in Figure 5, as well as Tables S3–S14.

Height

We first analyzed the set of 180 height associated loci identified by Lango Allen and colleagues [64], which explain about 10.5% of the total variance for height in the mapping population, or about 15% of heritability [65]. This dataset is an ideal first test for our methods because it contains the largest number of associations identified for a single phenotype to date, maximizing our power gain over single locus methods (Figure 2). In addition, Turchin and colleagues [28] have already identified a signal of pervasive weak selection at these same loci in European populations, and thus we should expect our methods to replicate this observation.

Of the 180 loci identified by Lango Allen and colleagues, 161 were present in our HGDP dataset. We used these 161 loci in conjunction with the allele frequency data from the HGDP dataset to estimate genetic values for height in the 52 HGDP populations. Although the genetic values are correlated with the observed heights in these populations, they are unsurprisingly imperfect predictions (see Figure S11 and Table S1, which compares our estimated genetic values to observed sex-average heights for the subset of HGDP populations with a close proxy in the dataset of Gustafsson and Lindenfors (2009) [66]).

We find a signal of excessive correlation with winter PC2 (Figure 6 and Table 2), but find no strong correlations with any other climatic variables. Our test also strongly rejects the neutral hypothesis, suggesting that our estimated genetic values are overly dispersed compared to the null model of neutral genetic drift and shared population history (Figure 3 and Table 2). These results are consistent with with directional selection acting in concert on alleles influencing height to drive differentiation among populations at the level of the phenotype.

We followed up on these results by conducting regional level analyses, which indicate that our signal of excess variance arises primarily from extreme differentiation among populations within Europe (Table 3). Analyses using the conditional multivariate normal model indicate that this signal is driven largely by divergence between the French and Sardinian populations, in line with Turchin et al's (2012) previous observation of a North-South gradient of height associated loci in Europe. We also find weaker signals of over-dispersion in other regions, but the globally significant statistic can be erased by removing either the French or the Sardinian population from the analysis, suggesting that the signal is primarily driven by differentiation among those two populations.

Skin pigmentation

We next analyzed data from a recent GWAS for skin pigmentation in an African-European admixed population of Cape Verdeans [67], which identified four loci of major effect that explain approximately 35% of the variance in skin pigmentation in that population after controlling for admixture proportion. Beleza et al (2013) report effect sizes in units of modified melanin (MM) index, which is calculated as , i.e. a higher MM index corresponds to darker skin, and a lower value to lighter skin.

We used these four loci to calculate a genetic skin pigmentation score in each of the HGDP populations. As expected, we identified a strong signal of excess variance among populations, as well as a strong correlation with latitude (Figure 6 and Table 2), again consistent with directional selection having acted on the phenotype of skin pigmentation to drive divergence among populations. Note, however, that this signal was driven entirely by the fact that populations of western Eurasian descent have a lower genetic skin pigmentation score than populations of African descent. Using only the markers from [67], light skinned populations in East Asian and the Americas have a genetic skin pigmentation score that is almost as high (dark) as that of most African populations, an effect that is clearly visible when we plot the measured skin pigmentation and skin reflectance of HGDP populations [68], [69] against their genetic values (see Figures S12 and S13). The correlation with latitude is thus weaker than one might expect, given the known phenotypic distribution of skin pigmentation in human populations [68], [70]. To illustrate this point further, we re-ran the analysis on a subsample of the HGDP consisting of populations from Europe, the Middle East, Central Asia, and Africa. In this subsample, the correlation with latitude, and signal of excess variance, was notably stronger (, ; , ).

This poor fit to observed skin pigmentation is due to the fact that we have failed to capture all of the loci that contribute to variation in skin pigmentation across the range of populations sampled, likely due to the partial convergent evolution of light skin pigmentation in Western and Eastern Eurasian populations [71]. Including other loci putatively involved in skin pigmentation (OCA2 and KITLG) [72], [73] decreases the estimated genetic pigmentation score of the other Eurasian populations (Figures S12 and S13 and Table S2), but we do not include these in our main analyses as they differ in ascertainment (and the role of KITLG in human pigmentation variation has been contested by [67]).

Within Africa, the San population has a decidedly lower genetic skin pigmentation score than any other HGDP African population. This is potentially in accordance with the observation that the San are more lightly pigmented than other African populations represented by the HGDP [68] and the observation that other putative light skin pigmentation alleles have higher frequency in the San than other African populations [71]. Although there is still much work to be done on the genetic basis of skin pigment variation within Africa, in this dataset a regional analysis of the six African populations alone identifies a marginally significant correlation with latitude (, ), and a signal of excess variance among populations (, ), suggesting a possible role for selection in the shaping of modern pigmentary variation within the continent of Africa.

Body mass index

We next investigate two traits related to metabolic phenotypes (BMI and Type 2 diabetes), as there is a long history of adaptive hypotheses put forward to explain phenotypic variation among populations, with little conclusive evidence emerging thus far. We first focus on the set of 32 BMI associated SNPs identified by Speliotes and colleagues [74] in their Table 1, which explain approximately 1.45% of the total variance for BMI, or about 2−4% of the additive genetic variance. Of these 32 associated SNPs, 28 were present in the HGDP dataset, which we used to calculate a genetic BMI score for each HGDP population. We identified no significant signal of selection at the global level (Table 2).

Our regional level analysis indicated that the mean genetic BMI score is significantly lower that expected in East Asia (; see also Figure 5 and Table S7), while marginal statistics identify excess intraregional variation within East Asia and the Americas (Table 3). While these results are intriguing, given the small fraction of the additive genetic variance explained by the ascertained SNPs and the lack of a globally significant signal or a clear ecological pattern or explanation, it is difficult to draw strong conclusions from them. For this reason BMI and other related traits will warrant reexamination as more association results arise and methods for analyzing association results from multiple correlated traits are developed.

Type 2 diabetes

We next investigated the 65 loci reported by Morris and colleagues [75] as associated with T2D, which explain of the total variance for T2D susceptibility, or about 8–9% of the additive genetic variance. Of these 65 SNPs, 61 were present in the HGDP dataset. We used effect sizes from the stage 1 meta-analysis, and where a range of allele frequencies are reported (due to differing sample frequencies among cohorts), we simply used the average. Where multiple SNPs were reported per locus we used the lead SNP from the combined meta-analysis. Also note that Morris and colleagues report effects in terms odds ratios (OR), which can be converted into additive effects by taking the logarithm (the same is true of the IBD data from [26], analyzed below).

The distribution of genetic T2D risk scores showed no significant correlations with any of the five eco-geographic axes we tested, and was in fact fairly underdispersed worldwide relative to the null expectation due to population structure (Table 2), suggesting we have little to no evidence that differential selection has influenced the distribution of T2D risk across human populations.

We note that our regional level analyses do find that European populations have a significantly lower T2D risk score than expected due to drift (, ), while Middle Eastern populations have significantly higher risk score than expected (, ). However, we are skeptical that this represents a meaningful signal of selection for two reasons. The first is that we have probed the data quite deeply despite seeing no evidence for adaptive differentiation at the global level. Second, expanding to the next closest region, an analysis in which we treated regional membership as the linear predictor was unable to find significant differentiation between Central Asia and Europe () or between Central Asia and the Middle East (). Our results are therefore consistent at most with a recent, but fairly weak selective event which influenced only European and Middle Eastern populations, but we do not feel our results count as strong evidence for this hypothesis.

A number of investigators have claimed that individual European GWAS loci for Type 2 Diabetes show signals of selection [63], [76]–[78], a fact that may be seen as support for the idea that genetic variation for T2D risk has been shaped by local adaptation, potentially consistent with a variation on the thrifty genotype hypothesis [79]. However, our result suggest that local adaptation has not had a large role in shaping the present day world-wide distribution of T2D susceptibility alleles (as mapped to date in Europe). One explanation of this discrepancy is that it is biologically unrealistic that the phenotype of T2D susceptibility would exhibit strong adaptive differentiation. Rather, local adaptation may have shaped some pleiotropically related phenotype (which shares only some of the loci involved). However, as seen in Figure 2, our methods have better power than single locus statistics so long as there is a reasonable correlation () between the focal phenotype and the one under selection. As such, the intersection of our results with previous studies support the idea that local adaptation has had little direct influence on the genetic basis of T2D or closely correlated phenotypes, but that a handful of individual SNPs associated with T2D may have experienced adaptive differentiation as a result of their function in some other phenotype.

IBD

Finally, we analyzed the set of associations reported for Crohn's Disease (CD) and Ulcerative Colitis (UC) [26]. Because CD and UC are closely connected phenotypes that share much of their genetic etiology, Jostins and colleagues used a likelihood ratio test of four different models (CD only, UC only, both CD and UC with equal effects on each, both CD and UC with independent effects) to distinguish which SNPs where associated with either or both phenotypes, and to assign effect sizes to SNPs (see their supplementary methods section 1d). We take these classifications at face value, resulting in two partially overlapping lists of 140 and 135 SNPs associated with CD and UC, which explain and of disease susceptibility variance respectively. Of these, there are 95 SNPs for CD and 89 SNPs for UC were present in our HGDP dataset, and these remaining SNPs on which our analyses are based explain 9% and 5.1% of the total variance. For now, we treat these sets of loci independently, and leave the development of methods that appropriately deal with correlated traits for future work.

We used these sets of SNPs to calculate genetic risk scores for CD and UC across the 52 HGDP populations. Both CD and UC showed strong negative correlations with summer PC2 (Figure 7), while CD also showed a significant correlation with winter PC1, and a marginally significant correlation with summer PC1 (Table 2).

We did not observe any significant statistics for either trait, either at the global or the regional level, suggesting that our environmental correlation signals most likely arise from subtle differences between regions, as opposed to divergence among closely related populations. Indeed, we find moderate signals of regional level divergence in Europe (UC: ), Central Asia (CD: ), and East Asia (CD: and UC: ; see also Figure 5 and Tables S13 and S14).

Discussion

In this paper we have developed a powerful framework for identifying the influence of local adaptation on the genetic loci underlying variation in polygenic phenotypes. Below we discuss two major issues related to the application of such methods, namely the effect of the GWAS ascertainment scheme on our inference, and the interpretation of positive results.

Ascertainment and Population Structure

Among the most significant potential pitfalls of our analysis (and the most likely cause of a false positive) is the fact that the loci used to test for the effect of selection on a given phenotype have been obtained through a GWAS ascertainment procedure, which can introduce false signals of selection if potential confounds are not properly controlled. We condition on simple features of the ascertainment process via our allele matching procedure, but deeper issues may arise from artifactual associations that result from the effects of population structure in the GWAS ascertainment panel. Given the importance of addressing this issue to the broader GWAS community, a range of well developed methods exist for doing GWAS in structured populations, and we refer the reader to the existing literature for a full discussion [80]–[86]. Here, we focus on two related issues. First, the propensity of population structure in the GWAS ascertainment panel to generate false positives in our selection analysis, and second, the difficulties introduced by the sophisticated statistical approaches employed to deal with this issue when GWAS are done in strongly structured populations.

The problem of population structure arises generally when there is a correlation in the ascertainment panel between phenotype and ancestry such that SNPs that are ancestry informative will appear to be associated with the trait, even when no causal relationship exists [81]. This phenomenon can occur regardless of whether the correlation between ancestry and phenotype is caused by genetic or environmental effects. To make matters worse, multiple false positive associations will tend to line up with same axis of population structure. If the populations being tested with our methods lie at least partially along the same axis of structure present in the GWAS ascertainment panel, then the ascertainment process will serve to generate the very signal of positive covariance among like effect alleles that our methods rely on to detect the signal of selection.

The primary takeaway from this observation is that the more diverse the array of individuals sampled for a given GWAS are with respect to ancestry, the greater the possibility that failing to control for population structure will generate false associations (or bias effect sizes) and hence false positives for our method.

What bearing do these complications have on our empirical results? The GWAS datasets we used can be divided into those conducted within populations of European descent and the skin pigmentation dataset (which used an admixed population). We will first discuss our analysis of the former.

The European GWAS loci we used were found in relatively homogeneous populations, in studies with rigorous standards for replication and control for population structure. Therefore, we are reasonably confident that these loci are true positives. Couple this with the fact that they were ascertained in populations that are fairly homogenous relative to the global scale of our analyses, and it is unlikely that population structure in the ascertainment panels is driving our positive signals. One might worry that we could still generate false signals by including European populations in our analysis, however many of the signals we see are driven by patterns outside of Europe (where the influence of structure within Europe should be much lessened). For height, where we do see a strong signal from within Europe, we use a set of loci that have been independently verified using a family based design that is immune to the effects of population structure [28].

We further note that for a number of GWAS datasets, including some of those analyzed here, studies of non-European populations have replicated many of the loci identified in European populations [87]–[93], and for many diseases, the failure of some SNPs to replicate, as well as discrepancies in effect size estimate, are likely due to simple considerations of statistical power and differences in patterns of LD across populations [94], [95]. This suggests that, at least for GWAS done in relatively homogenous human populations, structure is unlikely to be a major confounding factor.

The issue of population structure may be more profound for our style of approach when GWAS are conducted using individuals from more strongly structured populations. In some cases it is desirable to conduct GWAS in such populations as locally adaptive alleles will be present at intermediate frequencies in these broader samples, whereas they may be nearly fixed in more homogeneous samples. A range of methods have been developed to adjust for population structure in these setting [96]–[98]. While generally effective in their goal, these methods present their own issues for our selection analysis. Consider the extreme case, such as that of Atwell et al (2010) [19], who carried out a GWAS in Arabidopsis thaliana for 107 phenotypes across an array of 183 inbred lines of diverse geographical and ecological origin. Atwell and colleagues used the genome-wide mixed model program EMMA [83], [96], [97] to control for the complex structure present in their ascertainment panel. This practice helps ensure that many of the identified associations are likely to be real, but also means that the loci found are likely to have unusual frequencies patterns across the species range. This follows from the fact that the loci identified as associated with the trait must stand out as being correlated with the trait in a way not predicted by the individual kinship matrix (as used by EMMA and other mixed model approaches). Our approach is predicated on the fact that we can use genome-wide patterns of kinship to adjust for population structure, but this correction is exactly the null model that loci significantly associated with phenotypes by mixed models have overcome. For this reason, both the theoretical distribution of the statistic, as well as the empirical null distributions we construct from resampling, may be inappropriate.

The Cape Verde skin pigmentation data we used may qualify as this second type of study. The Cape Verde population is an admixed population of African/European descent, and has substantial inter-individual variation in admixture proportion. Due to its admixed nature, the population segregates alleles which would not be at intermediate frequency in either parental population, making it an ideal mapping population.

Despite the considerable population structure, the fact that recombination continues to mix genotypes in this population means that much of the LD due to the African/European population structure has been broken up (and the remaining LD is well predicted by an individual's genome-wide admixture coefficient). Population structure seems to have been well controlled for in this study, and a number of the loci have been replicated in independent admixed populations. While we think it unlikely that the four loci we use are false associations, they could in principle suffer from the structured ascertainment issues described above, so it is unclear that the null distributions we use are strictly appropriate. That said, provided that Beleza and colleagues have appropriately controlled for population structure, under neutrality there would be no reason to expect that the correlation among the loci should be strongly positive with respect to the sign of their effect on the phenotype, and thus the pattern observed is at least consistent with a history of selection, especially in light of the multiple alternative lines of evidence for adaptation on the basis of skin pigmentation [68]–[70], [99]–[101].

Further work is needed to determine how best to modify the tests proposed herein to deal with GWAS performed in structured populations.

Complications of Interpretation

Our understanding of the genetic basis of variation in complex traits remains very incomplete, and as such the results of these analyses must be interpreted with caution. That said, because our methods are based simply on the rejection of a robust, neutral null model, an incomplete knowledge of the genetic basis of a given trait should only lead to a loss of statistical power, and not to a high false positive rate.

For all traits analyzed here except for skin pigmentation, the within population variance for genetic value is considerably larger than the variance between populations. This suggests that much of what we find is relatively subtle adaptation even on the level of the phenotype, and emphasizes the fact that for most genetic and phenotypic variation in humans, the majority of the variance is within populations rather than between populations (see Figures S14–S19). In many cases, the influence of the environment likely plays a stronger role in the differences between populations for true phenotypes than the subtle differences we find here (as demonstrated by the rapid change in T2D incidence with changing diet, e.g. [102]). That said, an understanding of how adaptation has shaped the genetic basis of a wide variety of phenotypes is clearly of interest, even if environmental differences dominate as the cause of present day population differences, as it informs our understanding of the biology and evolutionary history of these traits.

The larger conceptual issues relate to the interpretation of our positive findings, which we detail below. A number of these issues are inherent to the conceptual interpretation of evidence for local adaptation [103].

Effect size heterogeneity and misestimation

In all of our analyses, we have simply extrapolated GWAS effect sizes measured in one population and one environment to the entire panel of HGDP populations. It is therefore prudent to consider the validity of this assumption, as well as the implications for our analyses when it is violated. Aside from simple measurement error, there are two possible reasons that estimated effect sizes from GWAS may not reflect the true effect sizes.

The first is that most GWAS hits likely identify tag SNPs that are in strong LD with causal sites that are physically nearby on the chromosome, rather than actual causal sites themselves [94], [95]. This acts to reduce the estimated effect size in the GWAS sample. More importantly for the interpretation of our signals, patterns of LD between tag SNPs and causal sites will change over evolutionary time, and so a tag SNP's allele frequencies will be an imperfect measure of the differentiation of the causal SNP over the sampled populations. This should lead to a reduction in our power to detect the effect of selection in much the same way that power is reduced when selection acts on a trait that is genetically correlated with the trait of interest (Figure 1B). This effect will be especially pronounced when the populations under study have a shorter scale of LD than the populations in which the effect have been mapped (e.g. when applying effect sizes estimated in Europe to populations of African descent). In the case that selection has not affected the trait of interest, the effect sizes have no association whatsoever with the distribution of allele frequencies across populations unless such an association is induced by the ascertainment process, as described above. Therefore, changes in the patterns of LD between identified tag SNPs and causal sites will not lead to an excess of false positives if the loci under study have not been subject to spatially varying selection pressures.

The second is that the actual value of the additive effect at a causal site may change across environments and genetic backgrounds due to genotype-by-genotype (i.e. functional epistasis) and genotype-by-environment interactions. Although the response at a given locus due to selection depends only the additive effect of the allele in that generation, the additive effect itself is a function of the environment and the frequencies of all interacting loci. As all of these can change considerably during the course of evolution, the effects estimated in one population may not apply in other populations, either in the present day, or over history of the populations [1], [104]. We first wish to stress that, as above, because our tests rely on rejection of a null model of drift, differences in additive effects among populations or over time will not lead to an excess of false positives, provided that the trait is truly neutral. Such interactions can, however, considerably complicate the interpretation of positive results. For example, different sets of alleles could be locally selected to maintain a constant phenotype across populations due to gene-by-environment interactions. Such a scenario could lead to a signal of local adaptation on a genetic level but no change in the phenotype across populations, a phenomenon known as countergradient variation [105].

It will be very difficult to know how reasonable it is to extrapolate effect sizes among populations without repeating measurements in different populations and different environments. Perhaps surprisingly, the existing evidence suggests that for a variety of highly polygenic traits, effects sizes and directions may be relatively consistent across human populations [87]–[95]. There is no particular reason to believe that this will hold as a general rule across traits or across species, and thus addressing this issue will require a great deal more functional genetic work and population genetic method development, a topic which we discuss briefly below in Future Directions.

Missing variants

As the majority of GWAS studies are performed in a single population they will often miss variants contributing to phenotypic variation. This can occur due to GxG or GxE interactions as outlined above, but also simply because those variants are absent (or at low frequency) due to drift or selection among the populations. Such cases will not create a false signal of selection if only drift is involved, however, they do complicate the interpretation of positive signals. A particularly dramatic example of this is offered by our analysis of skin pigmentation associated loci, whose frequencies are clearly shaped by adaptation. The alleles found by a GWAS in the Cape Verde population completely fail to predict the skin pigmentation of East Asians and Native Americans. This reflects the fact that a number of the alleles responsible for light skin pigmentation in those populations are not variable in Cape Verde due to the partially convergent adaptive evolution of light skin pigmentation [71]. As a result, when we take the Eurasian HGDP populations we see a significant correlation between genetic skin pigmentation score and longitude (), despite the fact that no such phenotypic correlation exists. While the wrong interpretation is easy to avoid here because we have a good understanding of the true phenotypic distribution, for the majority of GWAS studies such complications will be subtler and so care will have to be taken in the interpretation of positive results.

Loss of constraint and mutational pressure

One further complication in the interpretation of our results is in how loss of constraint may play a role in driving apparent signals of local adaptation. Traits evolving under uniform stabilizing selection across all populations should be less variable than predicted by our covariance model of drift, due to negative covariances among loci, and so should be underrepresented in the extreme tails of our environmental correlation statistics and the upper tail of . As such, loss of constraint (i.e. weaker stabilizing selection in some populations than others), should not on its own create a signal of local adaptation. While the loci underpinning the phenotype can be subject to more drift in those populations, there is no systemic bias in the direction of this drift. Loss of constraint, therefore, will not tend to create significant environmental correlations or systematic covariance between alleles of like effect.

An issue may arise, however, when loss of constraint is paired with biased mutational input (i.e. new mutations are more likely to push the phenotype in one direction than another [106]) or asymmetric loss of constraint (selection is relaxed on one tail of the phenotypic distribution). Under these two scenarios, alleles that (say) increase the phenotype would tend to drift up in frequency in the populations with loss of constraint, creating systematic trends and positive covariance among like effect alleles at different loci, and resulting in a positive signal under our framework. While one would be mistaken to assume that the signal was necessarily that of recent positive directional selection, these scenarios do still imply that selection pressures on the genetic basis of the phenotype vary across space. Positive tests under our methods are thus fairly robust in being signals of differential selection among populations, but are themselves agnostic about the specific processes involved. Further work is needed to establish whether these scenarios can be distinguished from recent directional selection based on only allele frequencies and effect sizes, and as always, claims of recent adaptation should be supported by multiple lines of evidence beyond those provided by population genomics alone.

Future directions

In this article we have focused on methods development and so have not fully explored the range of populations and phenotypes to which our methods could be applied. Of particular interest is the possibility of applying these methods to GWAS performed in other species where the ecological determinants of local adaptation are better understood [19], [20].

One substantial difficulty with our approach, particularly in its application to other organisms, is that genome-wide association studies of highly polygenic phenotypes require very large sample sizes to map even a fraction of the total genetic variance. One promising way to partially sidestep this issue is by applying methods recently developed in animal and plant breeding. In these genomic prediction/selection approaches, one does not attempt to map individual markers, but instead concentrates on predicting an individual's genetic value for a given phenotype using all markers simultaneously [107]–[109]. This is accomplished by fitting simple linear models to genome-wide genotyping data, in principle allowing common SNPs to tag the majority of causal sites throughout the genome without attempting to explicitly identify them [110]. These methods have been applied to a range of species, including humans [111]–[117], demonstrating that these predictions can potentially explain a relatively high fraction of the additive genetic variance within a population (and hence much of the total genetic variance). As these predictions are linear functions of genotypes, and hence allele frequencies, we might be able to predict the genetic values of sets of closely related populations for phenotypes of interest and apply very similar methods to those developed here. Such an approach may allow for substantial gains in power, as it would greatly increase the fraction of the genetic variance used in the analyses. However, if the only goal is to establish evidence for local adaptation in a given phenotype, then because measurements of true phenotypes inherently include all of the underlying loci, the optimal approach is to perform a common garden experiment and employ statistical methods such as those developed by Ovaskainen and colleagues [40], [41], [118], assuming such experiments can be done.

As discussed in various places above, it is unlikely that all of the loci underpinning the genetic basis of a trait will have been subject to the same selection pressures, due to their differing roles in the trait and their pleiotropic effects. One potential avenue of future investigation is whether, given a large set of loci involved in a trait, we can identify sets of loci in particular pathways or with a particular set of functional attributes that drive the signal of selection on the additive genetic basis of a trait.

Another promising extension of our approach is to deal explicitly with multiple correlated phenotypes. With the increasing number of GWAS efforts both empirical and methodological work are beginning to focus on understanding the shared genetic basis of various phenotypes [26], [119]. This raises the possibility that we may be able to disentangle the genetic basis of which phenotypes are more direct targets of selection, and which are responding to correlated selection on these direct targets (for progress along these lines using , see [40], [120]–[122]). Such tools may also offer a way of incorporating GxE interactions, as multiple GWAS for the same trait in different environments can be treated as correlated traits [123].

As association data for a greater variety of populations, species, and traits becomes available, we view the methods described out here as a productive way forward in developing a quantitative framework to explore the genetic and phenotypic basis of local adaptation.

Materials and Methods

Mean Centering and Covariance Matrix Estimation

Written in matrix notation, the procedure of mean centering the estimated genetic values and dropping one population from the analysis can be expressed as

where is an by matrix with on the main diagonal, and elsewhere.

In order to calculate the corresponding expected neutral covariance structure about this mean, we use the following procedure. Let be an by matrix, where each column is a vector of allele frequencies across the populations at a particular SNP, randomly sampled from the genome according to the matching procedure described below. Let and be the mean allele frequency in columns and of respectively, and let be a matrix such that . With these data, we can estimate as

This transformation performs the operation of centering the matrix at the mean value, and rooting the analysis with one population by dropping it from the covariance matrix (the same one we dropped from the vector of estimated genetic values), resulting in a covariance matrix describing the relationship of the remaining populations. This procedure thus escapes the singularity introduced by centering the matrix at the observed mean of the sample.

As we do not get to observe the population allele frequencies, the entries of are the sample frequencies at the randomly chosen loci, and thus the covariance matrix also includes the effect of finite sample size. Because the noise introduced by the sampling of individuals is uncorrelated across populations (in contrast to that introduced by drift and shared history), the primary effect is to inflate the diagonal entries of the matrix by a factor of , where is the number of chromosomes sampled in population (see the supplementary material of [46] for discussion). This means that our population structure adjusted statistics also approximately control for differences in sample size.

Standardized environmental variable

Given a vector of environmental variable measurements for each population, we apply both the and Cholesky tranformation as for the estimated genetic values

This provides us with a set of adjusted observations for the environmental variable which can be compared to the transformed genetic values for inference. This step is necessary as we have rotated the frame of reference of the estimated genetic values, and so we must do the same for the environmental variables to keep them both in a consistent reference frame.

Identifying Outliers with Conditional MVN Distributions

As described in the Results, we can use our multivariate normal model of relatedness to obtain the expected distribution of genetic values for an arbitrary set of populations, conditional on the observed values in some other arbitrary set.

We first partition our populations into two groups, those for which we want to obtain the expected distribution of genetic values (group 1), and those on which we condition in order to obtain this distribution (group 2). We then re–estimate the covariance matrix such that it is centered on the mean of group 2. This step is necessary because the amount of divergence between the populations in group 1 and the mean of group 2 will always be greater than the amount of divergence from the global mean, even under the neutral model, and our covariance matrix needs to reflect this fact in order to make accurate predictions. We can obtain this re-parameterized matrix as follows. If is the total number of populations in the sample, then let be the number of populations in group one, and let be the number of populations in group 2. We then define a new matrix such that the columns corresponding the populations in group one have 1 on the diagonal, and 0 elsewhere, while the columns corresponding to group two have on the diagonal, and elsewhere. We can then re–estimate a covariance matrix that is centered at the mean of the populations in group 2. Recalling our matrices and from (17), this matrix is calculated as

where we write to indicate that it is a covariance matrix that has been re-centered on the mean of group two.

Once we have calculated this re–centered covariance matrix, we can use well known results from multivariate normal theory to obtain the expected joint distribution of the genetic values for group one, conditional on the values observed in group two.

We partition our vector of genetic values and the re–centered covariance matrix such that

and

where and are vectors of genetic values in group 1 and 2 respectively, and , and are the marginal covariance matrices of populations within group 1, within group 2, and across the two groups, respectively. Letting (i.e. the sum of the elements of ), we wish to obtain the distribution

where and give the expected means and covariance structure of the populations in group 1, conditional on the values observed in group 2. These can be calculated as

and

where the one vectors in line (23) are of length and respectively.

This distribution is itself multivariate normal, and as such this framework is extremely flexible, as it allows us to obtain the expected joint distribution for arbitrary sets of populations (e.g. geographic regions or continents), or for each individual population. Further,

and

where denotes the elements of . In words, the conditional expectation of the mean estimated genetic value across group 1 is equal to the mean of the conditional expectations, and its variance is equal to the mean value of the elements of the conditional covariance matrix. As such we can easily calculate a Z score (and corresponding p value) for group one as a whole as

This Z score is a normal random variable with mean zero, variance one under the null hypothesis, and thus measures the divergence of the genetic values between the two populations relative to the null expectation under drift. Note that the observation of a significant Z score in a given population or region cannot necessarily be taken as evidence that selection has acted in that population or region, as selection in some of the populations on which we condition (especially the closely related ones) could be responsible for such a signal. As such, caution is warranted when interpreting the output of these sort of analyses, and is best done in the context of more explicit information about the demographic history, geography, and ecology of the populations.

The Linear Model at the Individual Locus Level

As with our excess variance test, explored in the main text, it is natural to ask how our environmental correlation tests can be written in terms of allele frequencies at individual loci.

As noted in (8), we can obtain for each underlying locus a set of transformed allele frequencies, which have passed through the same transformation as the estimated genetic values. We assume that each locus has a regression coefficient

where is shared across all loci so that

where the are independent and identically distributed residuals. We can find the maximum likelihood estimate by treating as the linear predictor, and taking the regression of the combined vector , across all populations and loci, on the combined vector . As such

we can decompose this into a sum across loci such that

As noted in (8), our transformed genetic values can be written as

and so the estimated slope () of our regression () is

Comparing these equations, the mean regression coefficient at the individual loci (31) and the regression coefficient of the estimated genetic values (33) are proportional to each other via a constant that is given by one over two times the sum of the effect sizes squared (i.e. ). Our test based on estimating the regression of genetic values on the environmental variable is thus mathematically equivalent to an approach in which we assume that the regression coefficients of individual loci on the environmental variable are proportional to one another via a constant that is a function of the effect sizes. Such a relationship can also be demonstrated for the correlation coefficient () calculated at the genetic value level and at the individual locus level (this is not necessarily true for the rank correlation ), however the algebra is more complicated, and thus we do not show it here.

This is in contrast to the enrichment statistic we compute for the power simulations, in which we assume that the correlations of individual loci with the environmental variable are independent of one another, and then perform a test for whether more loci individually show strong correlations with the environmental variable than we would expect by chance.

HGDP Data and Imputation

We used imputed allele frequency data in the HGDP, where the imputation was performed as part of the phasing procedure of [59], as per the recommendations of [124]. We briefly recap their procedure here:

Phasing and imputation were done using fastPHASE [125], with the settings that allow variation in the switch rate between subpopulations. The populations were grouped into subpopulations corresponding to the clusters identified in [62]. Haplotypes from the HapMap YRI and CEU populations were included as known, as they were phased in trios and are highly accurate. HapMap JPT and CHB genotypes were also included to help with the phasing.

Choosing Null SNPs

Various components of our procedure involve sampling random sets of SNPs from across the genome. While we control for biases in our test statistics introduced by population structure through our matrix, we are also concerned that subtle ascertainment effects of the GWAS process could lead to biased test statistics, even under neutral conditions. We control for this possibility by sampling null SNPs so as to match the joint distribution of certain properties of the ascertained GWAS SNPs. Specifically, we were concerned that the minor allele frequency (MAF) in the ascertainment population, the imputation status of the allele in the HGDP datasets, and the background selection environment experienced at a given locus, as measured by B value [61], might influence the distribution of allele frequencies across populations in ways that we could not predict.

We partitioned SNPs into a three way contingency table, with 25 bins for MAF (i.e. a bin size of 0.02), 2 bins for imputation (either imputed or not), and 10 bins for B value (B values range from 0 to 1, and thus our bin size was 0.1). For each set of null genetic values, we sampled one null SNP from the same cell in the contingency table as each of the GWAS SNPs, and assigned this null SNP the effect size associated with the GWAS SNP it was sampled to match. While we do not assign effect sizes to sampled SNPs used to estimate the covariance matrix (instead simply scaling by a weighted sum of squared effect sizes, which is mathematically equivalent under our assumption that all SNPs have the same covariance matrix), we follow the same sampling procedure to ensure that describes the expected covariance structure of the GWAS SNPs.

For the skin pigmentation GWAS [67] we do not have a good proxy present in the HGDP population, as the Cape Verdeans are an admixed population. Cape Verdeans are admixed with African ancestry, and European ancestry in the sample obtained by [67] (Beleza, pers. comm., April 8, 2013). As such, we estimated genome wide allele frequencies in Cape Verde by taking a weighted mean of the frequencies in the French and Yoruban populations of the HGDP, such that . We then used these estimated frequencies to assign SNPs to frequency bins.

[67] also used an admixture mapping strategy to map the genetic basis of skin pigmentation. However, if they had only mapped these loci in an admixture mapping setting we would have to condition our null model on having strong enough allele frequency differentiation between Africans and Europeans at the functional loci for admixture mapping to have power [126]. The fact that [67] mapped these loci in a GWAS framework allows us to simply reproduce the strategy, and we ignore the results of the admixture mapping study (although we note that the loci and effect sizes estimated were similar). This highlights the need for a reasonably well defined ascertainment population for our approach, a point which we comment further on in the Discussion.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20

Attachment 21

Attachment 22

Attachment 23

Attachment 24

Attachment 25

Attachment 26

Attachment 27

Attachment 28

Attachment 29

Attachment 30

Attachment 31

Attachment 32

Attachment 33

Attachment 34

Attachment 35