Simple scaling laws control the genetic architectures of human complex traits

Authors: Yuval B. Simons, Hakhamanesh Mostafavi, Huisheng Zhu, Courtney J. Smith, Jonathan K. Pritchard, Guy Sella

Abstract

Genome-wide association studies have revealed that the genetic architectures of complex traits vary widely, including in terms of the numbers, effect sizes, and allele frequencies of significant hits. However, at present we lack a principled way of understanding the similarities and differences among traits. Here, we describe a probabilistic model that combines the effects of mutation, drift, and stabilizing selection at individual sites with a genome-scale model of phenotypic variation. In this model, the architecture of a trait arises from the distribution of selection coefficients of mutations and from two scaling parameters. We fit this model for 95 highly polygenic quantitative traits of different kinds from the UK Biobank. Notably, we infer that all these traits have fairly similar, though not identical, distributions of selection coefficients. This similarity suggests that differences in architectures of highly polygenic traits arise mainly from the two scaling the mutational target size and heritability per site, which vary by orders of magnitude among traits. When these two scale factors are accounted for, we find that the architectures of all 95 traits are very similar.

Introduction

A central goal of genetics is to understand how genetic variation maps to phenotypic variation. Starting in the late 20th century, there was huge progress toward identifying the genes for Mendelian traits. But most phenotypic variation in humans is genetically complex, and it is only in the last 15 years that genome-wide association studies (GWAS) have started to reveal the genetic basis of variation in a wide array of complex traits [1]. These studies have now identified tens of thousands of robust associations between genetic variants and a wide array of traits and diseases.

One intriguing observation from this work is the striking variation in genetic architecture among complex traits [2–5]. (Here, we use the term architecture to refer to the numbers of causal variants and their joint distribution of allele frequencies and effect sizes.) Traits have been found to vary in all aspects of genetic architecture, the number and magnitude of significant signals found at a given sample size [6]; the fraction of heritability explained by lead GWAS signals [3]; the allele frequency distributions of significant variants [4,7]; the estimated numbers of causal variants [8,9]; and the SNP-based heritability [10,11].

Nonetheless, diverse traits do show important similarities. First, most complex traits are influenced by large numbers of variants with small effects, only a small fraction of which can be confidently detected at current sample sizes [12–14]. Indeed, even relatively “simple” complex traits such as molecular biomarkers are highly polygenic with ∼10^4^ causal variants spread widely across the genome, compared to ∼10^5^ or more variants for traits such as height or BMI [2,9,15–17].

Second, the distributions of effect sizes of causal variants are not fit well using standard modeling assumptions such as normal distributions. Instead, effect sizes typically span several orders of magnitude, much like power-law distributions [18,19].

Third, trait-associated variants are often highly i.e., they influence many traits simultaneously. Many pairs of traits show significant genetic correlations, indicating that allelic impacts are often shared [16,20,21]; moreover, whenever different traits are mediated through overlapping cell types or pathways, we can expect that they will share many of the same regulatory variants even if the directions of effects are uncorrelated [5,14,16,22].

Fourth, selection plays a central role in shaping complex trait architecture. Evolutionary theory and empirical evidence indicate that variants with phenotypic effects would usually be under selection and, in particular, that selection is usually stronger for larger-effect variants [23]. Consistent with this, variants with larger effect sizes tend to be at lower frequencies, suggesting that selection prevents such variants from reaching high frequencies [24–26]. Since heritability depends on both effect sizes and allele frequencies, an important consequence is that the genes that are most important for a trait contribute less to heritability than would be expected in the absence of selection, thus flattening the heritability distribution across genes [17,27,28].

Here we develop a principled approach for understanding similarities and differences in genetic architecture. Specifically, we want to understand how the population genetic processes of mutation, selection, and drift alongside properties of individual traits determine the numbers of variants, as well as the joint distributions of allele frequencies and effect sizes. What features of these processes are shared across traits? And which are different? And, consequently, what features of the genetic architecture are shared or differ among traits?

To answer these questions, we require a model for how population genetic processes shape complex trait architecture. Current models differ primarily in their assumptions about the relationship between selection on alleles (or alternatively their frequencies), and the effect sizes of those alleles on a trait of interest [23]. The heuristic ‘α-model’, developed for estimating SNP heritability, assumes a particular parametric relationship between allele frequencies and effect sizes [25,29,30]. While the α-model is motivated by the observed inverse relationship between effect size and frequency, the precise functional form is arbitrary. In turn, several evolutionary models postulate particular parametric relationships between the strength of selection on alleles and their effect sizes, and then rely on explicit population genetics models to derive the relationship between allele frequencies and effect sizes and other aspects of genetic architecture [31–36]. These models, however, differ in their predictions about architecture, owing to the various ad-hoc parametric relationships they assume.

Simons et al. (2018) introduced an evolutionary model that moves beyond ad-hoc choices by deriving the relationship between selection on alleles and their effects on a trait under an explicit, interpretable, biological model [27]. Motivated by extensive evidence that many quantitative traits are subject to stabilizing selection, where fitness declines with displacement from an optimal trait value [23,37,38], and that genetic variation affecting one trait often affects many others [5,16,22], they modeled selection on alleles that arises from stabilizing selection in a multi-dimensional trait space. They then used an explicit population genetic model to derive the genetic architecture of a focal trait with mutation, genetic drift, and stabilizing selection in a multi-dimensional trait space.

As we will show in the next section, the Simons et al. model can be reframed as a generative (statistical) model for the genetic architecture of a continuous complex trait, which depends on a few biologically interpretable parameters. We next describe how to infer these parameters and test the model fit based on data from GWAS. Applying our inference to 95 highly-polygenic quantitative traits of different kinds from UK Biobank, we show that this model provides an excellent fit to the data. Surprisingly, we find that most variation in architecture among traits is explained by differences in just two scaling the mutational target size and heritability per site.

Results

A population genetic model of complex traits. As a starting point, we assume that phenotypic variation exists in a high-dimensional trait space under stabilizing selection. Here we outline key elements of the model; further details and biological motivation can be found in S1 Note and in [23,27].

We model each person’s phenotype as a point in an n-dimensional trait space, and assume that this dimension is high (n≥10). To model stabilizing selection we assume that there is an optimal phenotype, and that fitness decreases with Euclidean distance from the optimum (Fig 1A).

The phenotypic effect of each variant is represented by a random vector in the n-dimensional trait space, namely b⇀=(b1,b2,...bn), where bi is the additive effect of the derived allele on the i-th trait. We assume that a person’s phenotype arises from their genotype according to the standard additive model in n-dimensions: it is a vector sum over the effects of all variants plus a random vector representing the environmental effects [39].

The model thus far is mathematically similar to Fisher’s Geometric model [40], which Fisher and others used to study adaptive processes [41], but we consider a different question and a different evolutionary setting. We focus on the genetic architecture of a single highly polygenic trait that arises in the balance between mutation, stabilizing selection in the multidimensional trait space, and genetic drift.

Mutation, selection, and drift at individual sites. Each generation, mutation introduces new trait-affecting variants into the population at a rate μ per site, per gamete, per generation. The long-term fate of variants is determined by the combined action of selection and drift.

Under stabilizing selection, at equilibrium, selection holds the population’s phenotypic mean very close to the optimal phenotype, and thus acts against mutations (and against variation in general). The strength of selection, s, acting against a variant is proportional to its squared magnitude in the n-dimensional trait space (Fig 1B):

with VS reflecting the width of the fitness function around the optimum.

Given s we can compute the present day allele frequency distribution, as follows. When traits are subject to stabilizing selection, selection at individual sites is under-dominant, meaning that selection acts against minor alleles, regardless of the direction of effect [27,42,43]. At strongly selected sites, this approximates the standard model of selection against deleterious alleles. Specifically, the expected change in allele frequency at an autosomal site in a single generation, given current derived allele frequency q is

Meanwhile, the variance in the change in allele frequencies, i.e., drift, scales inversely with population size [44]. Hence, the distribution of present-day allele frequencies is the result of a stochastic process including past mutations, selection, and drift – which depends on the history of population sizes. For our analysis here, we computed the distribution of present day allele frequencies under the stabilizing selection model using a demographic model estimated for the British population [45]. As expected, strongly selected variants (large s) tend to be rare, while nearly-neutral variants (small s) can drift to high frequencies (Fig 1C).

The relationship between selection and effect sizes. Next we need to understand how selection in the multi-dimensional trait space relates to the genetic architecture of a single focal trait of interest. Without loss of generality, we focus on the first dimension in the n-dimensional trait space, and to simplify the notation we denote the effect size b1 of a variant on trait 1 simply as b.

Causally, effect sizes determine selection coefficients. However, since we want to describe the co-distribution of frequencies and effect sizes we need to invert this specifically, we need the conditional distribution of the effect size on trait 1 given s. This conditional distribution reflects uncertainty about the projection of b⇀ onto the first dimension if all we know is s (or equivalently ||b⇀||2). Fortunately, when the number of traits is sufficiently large, this conditional distribution is well approximated by a Normal distribution (see [27] and Sect 1.3 in S1 Note),

where c=VS/n. Intuitively, variants under weak selection (small s) tend to have small squared effect sizes (b^2^) and variants under strong selection (large s) tend to have larger squared effect sizes (Fig 1D).

With the distributions for b and q given s, we can now compute the expected per site contribution to phenotypic variance as a function of s, given by E[2b^2^q(1 − q)|s]. Under non-equilibrium demography, the expectation does not have a simple form but it is plotted in Fig B in S1 Note. At sites where selection is weak, b^2^ is small, and these sites contribute little to the additive genetic variance, VA. When selection is strong, b^2^ is large, but selection holds q(1–q) low, and these effects cancel out, so these sites are capped in terms of how much they can contribute to VA [27].

Single-site dynamics and heritability. Moving from single sites to a genome-wide model, let L be the number of sites in the genome at which mutations can affect trait 1 (more precisely, we require mutations to have effects that exceed some predefined small value); we refer to L as the mutational target size. We use f(s) to denote the unknown distribution of selection coefficients of mutations at these L sites. Then the expected total additive genetic variance for trait 1 is given by

Next, we rescale b from the original but arbitrary measurement units into units of standard deviations of the trait we define β=b/VP, where VP is the phenotypic variance. Dividing both sides by VP, and noting that VA/VP is the (narrow sense) heritability h^2^, we can relate heritability to the site-level

This equation expresses the key relationship between heritability (h^2^), mutational target size (L), and the expected contribution to variance per site.

Finally, Eq 3 can be rewritten in terms of β and population genetics parameters (for details see Sect 1.5 in S1 Note):

where k is a constant that depends on f(s) and demography and is approximately 1, and μ is the mutation rate. Crucially, Eq 6 shows that the trait’s heritability per site, h^2^/L, is a fundamental scaling factor that relates selection on alleles to their effects on the trait.

Together, these results provide a generative model for the genetic architecture of a complex trait (Fig 1E). Assuming that the demographic history and mutation rate per site are known in advance, this model is fully specified in terms of three the mutational target size, L; the heritability per site, h^2^/L; and the distribution of selection coefficients, f(s). We now describe how we estimate these from GWAS data.

Inference of model parameters from GWAS data. In principle we would want to perform inference using all causal variants, but this is technically challenging since most causal sites have very small effect sizes; hence there is great uncertainty about which sites are causal and their true effect sizes. As a tractable alternative, we restricted our inference to the independent genome-wide significant hits for each trait. We account for this restriction in the inference by noting that we only observe the subset of sites for which the absolute GWAS z-score exceeds 5.45, corresponding to the conventional significance threshold of p<5×10−8, and the minor allele frequency (MAF) exceeds 1%, corresponding to the imputation threshold (see Sect 2.1 in S1 Note). We performed simulations to illustrate the changes in architecture at GWAS hits as a function of each of the main model components (Fig 2; Sect 5 in S1 Note). As expected from

• When selection is weak, causal variants can drift to high frequencies, and most significant hits are at common variants. Conversely, when selection is strong, there is a greater fraction of rare variants among the significant hits, and an inverse relationship between effect size and allele frequency.
• When traits have high heritability per site h^2^/L, the squared effect sizes and z-scores tend to be larger, there are more genome-wide significant hits and they explain a greater proportion of heritability, compared with traits with low h^2^/L.
• When traits have a large mutational target size L (holding h^2^/L and f(s) constant), there are more causal variants, and more genome-wide significant hits, but the distribution of allele frequencies and effect sizes, and the proportion of heritability explained by hits, are unaffected.

We implemented a maximum likelihood method that estimates the components of our model from the joint distribution of q and |z| across significant hits (Fig 2D–2F; Sect 5 in S1 Note). We fit f(s) using a spline function with four knots (see Sect 4.8 in S1 Note for how the number and position of knots were chosen), thus our full model includes six parameters per four for f(s), as well as h^2^/L, and L.

We tested this method using simulated GWAS data under a variety of parameter values. We find that even with modest numbers of hits (∼100) the method provides accurate estimates, while the estimates are noisy for traits with fewer hits (e.g., panel 2E). It may seem surprising that we can estimate f(s) from relatively few observations, but each variant carries considerable information about the strength of the allele frequency bounds s for that variant from above and the effect size bounds s both from below and above, such that jointly they are highly informative (Fig D in S1 Note). Given these results we analyzed traits with at least 100 hits.

We also observed that the data are less informative at both ends of the range of possible selection GWAS has low power to detect very strongly selected variants (s≳10−2) as their allele frequencies are too low, and low power to detect effectively-neutral variants (s≲10−5) as their effect sizes are too small. We therefore implemented a regularization penalty to constrain f(s) to sensible values at these extremes (Sect 3.8 in S1 Note).

Dataset of 95 quantitative traits from the UK BioBank. We selected traits from the UK Biobank for analysis, as follows (Sect 2.4 in S1 Note). Since our model is most directly applicable to quantitative continuous traits, we restricted our analysis to such traits. We identified independent lead variants for each trait using COJO [46]. Since low-frequency variants are often poorly imputed, we removed hits with MAF<1% (see Sect 6.3 in S1 Note for how we account for this in the inference). We excluded traits for which more than 10% of hits fell within a single LD block, as well as hits in regions of extremely high LD (LD score >500). As noted above, we restricted ourselves to traits with at least 100 independent hits; doing so implies that the traits are among the more polygenic and heritable traits in UKBB. For each trait we recorded the number of hits, and the estimated allele frequency, z-score, and effective sample size for each hit.

Our inference uses the architecture of GWAS hits as a proxy for the architecture of the causal variants they tag. In Sect 7 in S1 Note we tested the validity of this approximation using data analysis and simulations. For example, we started with the genotypes from UKB, picked variants to be causal and assigned them effect sizes in various ways, assigned phenotypes to individuals based on the additive model, and then performed GWAS and ran COJO, allowing us to compare the inferred number of hits and their distribution of frequencies and effect sizes with the known underlying architecture. Our results indicate that under sensible assumptions the GWAS hits picked by COJO faithfully reflect the architecture of the underlying causal variants.

Our criteria resulted in a list of 95 traits, with a range of 100 to 1,760 hits per trait (mean=495). These traits include 40 morphometric traits, of which 26 are related to body weight or adiposity (e.g., BMI, waist circumference and birth weight) as well as 14 others (e.g., height, bone mineral density and hand grip strength). The traits also include 27 blood phenotypes (e.g. platelet traits, lymphocyte count, and hemoglobin measurements), and 12 molecular traits sampled from blood or urine (e.g., IGF-1, triglycerides and calcium levels). Additionally, we have 9 cardiovascular traits, including pulse rate, blood pressure measurements, and pulmonary function traits. Lastly, we include 6 ophthalmologic traits, and 1 behavioral trait (age at first sexual intercourse).

Distributions of trait parameters. We applied our inference to all 95 traits. Fig 3A shows the estimated distributions of selection coefficients that vary in polygenicity and kind, from calcium levels to BMI, and are representative of the range spanned by all 95 traits (see Fig F in S1 Note). Note that although our inference is based on significant hits, f(s) represents the distribution of s among new mutations.

Our point estimates of the distribution (solid lines) vary considerably among traits, but they share important features. Most trait-affecting mutations are under weak selection, with modes ranging between 10^−4^ – 10^−5.5^. In this range, the strength of selection is roughly comparable to genetic drift (∼10−4.5 per generation), consistent with the observation that many GWAS variants are common, and more generally that much of the heritable variation in complex traits arises from common variants [3,4,10,12]. However, the distributions have a substantial tail in the strong selection range (s > 10^−3^), and therefore span multiple orders of magnitude. Since selection coefficients span multiple orders of magnitude, so too should effect sizes. This echoes recent results showing that distributions of effect sizes for complex traits do indeed span multiple orders of magnitude (in contrast to the normal distribution which has often been assumed in statistical genetics models) [18,19].

Although our point estimates of the distribution (solid lines) vary considerably among traits, their confidence envelopes largely overlap (shades). With our inferences being based on the architecture of significant hits, the overlap among confidence envelopes suggests that this architecture is largely insensitive to the differences reflected in our point estimates (also see Sect 5.2 in S1 Note). We thus conjectured that we could build a unified model by assuming a Single Shared Distribution for f(s), which we refer to as the SSD, instead of assuming separate Trait-Specific Distributions (TSDs). The SSD is shown in Fig 3B. We also inferred the SSD for a subset of the traits that were chosen to minimize genetic correlations among traits (15 traits with |rg| < 0.2 for any pair) and found that the estimated SSD is insensitive to correlations among traits in the full dataset (Sect 6.4 in S1 Note). As we will show, the SSD provides a useful approximation for the architecture of individual traits, while greatly cutting down the number of model parameters and highlighting important shared features of trait architecture.

In contrast to f(s), our estimates of heritability h^2^ and of the mutational target size L span a wide range that far exceeds the CIs for individual traits (Fig 3C and S1 Table). For example, among the traits in Fig 3C, BMI (541 hits) has an estimated target size of 220 MB, one of the highest estimated target sizes, in contrast to urea (231 hits) with a target size of 4 MB, one of the lowest estimated target sizes. These results are broadly consistent with expectations from previous studies of polygenicity showing that morphological traits including height and BMI have many more contributing variants than do molecular traits, illustrated here by urea and calcium [2,9]. Our estimates of heritability vary by over one order of magnitude among traits, and are concordant with previous estimates (S1 Table, Fig J in S1 Note; [10]). Our estimates of the ratio h^2^/L vary less than h^2^ and L alone, plausibly because conditioning on traits with ≥100 hits biases us toward traits with greater (and thus a more restricted range of) h^2^/L (see Fig 2B and below).

Quantifying model fit. Next, we assessed the fit of our models to the genetic architecture observed in GWAS, including how the fit is affected by using the SSD approximation (Sect 4 in S1 Note). To do so, we computed a measure of model fit using the predicted distribution of z-scores given the allele frequencies and study size. For each variant i, we computed what we refer to as a residual p-value: Pr(|z|>zi | zi>5.45,qi, Model), where zi and qi are the observed z-score and frequency of SNP i, respectively, and Model indicates the SSD or TSD model. The SSD model fit depends on only one trait-specific parameter, h^2^/L, as f(s) is shared across traits. In contrast, the TSD model fit depends on five parameters per h^2^/L and four parameters to fit f(s). The estimated target size, L, does not affect the fit by this measure and more generally (because given the other model parameters, L is estimated by matching the expected and observed number of hits).

The residual p-value has a simple If we correctly model the distribution of z-scores among significant hits, then the distribution of residual p-values will be uniform between 0 and 1. If the observed z-scores are too small then the residual p-values will skew toward 1, and if the z-scores are too high they will skew toward 0. To avoid overfitting, we split the genome into approximately independent blocks [47], each time inferring the model on 90% of the blocks and computing residual p-values for the held-out 10%.

We first considered the fit of our model for height, the trait with the greatest number of hits in our dataset (Fig 4A). For height, the distribution of z-scores across the 1760 hits is fit essentially perfectly by the SSD model (Fig 4A). In contrast, two simpler heuristic models provide a poor fit to the distribution of z-scores (Fig 4A). First, when we assumed that effect sizes are normally distributed, the resulting residual p-values deviate greatly from uniformity with many extremely small p-values. We also considered a version of the α-model with a normal density of effect sizes conditional on allele frequencies [30]. By fitting the inverse relationship between allele frequencies and effect sizes, the alpha model improves the overall fit, but still has an excess of tiny p-values. In both cases, the underlying normal distribution is too narrow to accommodate the wide variation in observed effect sizes, including many hits close to the significance threshold and a minority of much stronger hits.

For BMI, with 541 hits, our SSD model fits most of the distribution well, but the top hits are larger than expected (Fig 4B). In particular, the residual p-value of one SNP is significant even after Bonferroni correction. This outlier represents the well-known FTO signal that was detected even in very early GWAS studies [48,49].

We find 14 additional outlier SNPs for a variety of other traits. The 14 outliers include both missense and noncoding variants, and are all found near genes previously implicated in the biology of the relevant traits (Table A in S1 Note and S2 Table). Our inferences are insensitive to the exclusion of these outliers (see Sect 4.5 in S1 Note). We hypothesize that these outliers violate our model assumptions in some way that allows them to be common despite having large effects. For example, they might have much smaller pleiotropic effects than most other variants affecting those traits, leading to weaker selection than expected given their effect sizes [50]. Alternatively, they may have been targets of strong positive or balancing selection that allowed them to reach high frequencies despite their large effect sizes.

We next performed goodness-of-fit tests for each trait to determine whether the overall distributions of residual p-values match the expected uniform distribution, using Kolmogorov-Smirnov statistics (Fig 4C; Sect 4 in S1 Note). We find that the α-model fits the data far worse than either the TSD and SSD models; notably, we can reject it for 93 out of 95 traits at a FDR of 0.05 (Sect 4 in S1 Note). In contrast, both the TSD and SSD models appear to fit the data well, with a modest inflation of p-values, perhaps relating to simplifying assumptions of the model and inference (Fig 4C). For most traits, the TSD model fits slightly better than the SSD model, with no trait rejected for the TSD and 26 traits for the SSD at a FDR of 0.05 (Sect 4 in S1 Note), indicating that the TSD correctly identifies some degree of trait specific signal. Nonetheless, the SSD model fits the data quite well and is far more parsimonious, with one trait specific parameter compared to five for the TSD model.

Prediction of allele ages. The previous results show that the model provides a good empirical fit to the data. We next wanted to evaluate whether it can also predict the evolutionary processes underlying the genetic architecture. To this end, we turned to an entirely different type of predictions from our model. Conditional on allele frequency, deleterious alleles tend to be younger than neutral alleles [24,51]; and our parameter estimates can be used to predict the extent of this effect for GWAS hits. We compared our model’s predictions to allele ages estimated from a reconstruction of the ancestral recombination graph using Relate [45] (see Sect 8 in S1 Note for details, including a correction for the estimator of allele ages).

Looking at the allele ages of GWAS hits for all 95 traits as estimated by Relate, we see that they are much younger than frequency-matched, putatively neutral alleles (Fig 5). We estimate that the median age of a GWAS hit variant is 137,000 years, compared to 594,000 years for matched neutral variants. These observations highlight the competing influences of selection and drift on GWAS due to selection GWAS hits are much younger than matched neutral variants, yet at the same time selection is weak enough that most GWAS hits are fairly old, predating the out-of-Africa bottleneck. This finding echoes previous ones showing that many common GWAS hits are shared among African and non-African populations [52,53].

Most importantly, the shift in ages of GWAS hits compared to matched control SNPs is predicted well by the SSD model, indicating that this model captures the right magnitude of selection coefficients. (For results about our ability to distinguish selection coefficients with this analysis see Fig X in S1 Note.) Since the allele ages inferred by Relate are estimated from local haplotype structure, information that isn’t used by our inference, this concordance provides an external validation of the SSD model.

Simple scaling rules control differences in trait architectures. The fit of the SSD model suggests an intriguing that the differences in genetic architecture among traits are primarily due to just two trait-specific the heritability per site (h^2^/L) and the mutational target size (L).

To see why, consider the full genetic architecture of a trait, of all variants regardless of whether they can be detected by GWAS. First, for a given demographic history, the distribution of allele frequencies depends only on the distribution of selection coefficients, f(s). Hence, under the SSD approximation, the allele frequency distribution is shared across traits and can be predicted from our estimate of f(s) (black line, Fig 6A).

Next, while f(s) represents the distribution of selection coefficients among new mutations, strongly selected variants are less likely to reach high frequencies. Hence, the distribution of selection coefficients shifts to be smaller with increasing allele frequencies, as shown in Fig 6B.

Given the distributions of selection coefficients at different allele frequencies from Fig 6B, we can compute the distribution of squared effect sizes as a function of allele frequencies, by integrating Eq 6 over s. Crucially, under the SSD model, these distributions are identical across traits, if the effect sizes on the x-axis are scaled in terms of h^2^/L (black lines, Fig 6C).

How do these distributions affect GWAS hits? Unlike the underlying distributions, the power to detect significant variants in GWAS depends on the actual squared effect sizes β2, not scaled by h^2^/L (and it depends on allele frequency and sample size). Consequently, there is more power to detect variants for traits with higher h^2^/L – this is intuitive, because higher h^2^/L implies that on average each site explains more variance in the trait. This is illustrated in Fig 6C: for traits with high h^2^/L, variants within both light and dark blue regions are genome-wide significant, but for traits with low h^2^/L, only variants within the dark blue regions are detected.

Moreover, even though the underlying distribution of causal variant allele frequencies is shared among traits, the frequency distribution of hits is predicted to vary. For traits with high h^2^/L there is relatively more power to detect low frequency variants than for traits with low h^2^/L (Fig 6A).

Lastly, the second scaling parameter, L, represents the mutational target size. Conditional on h^2^/L, changing L only changes the numbers of causal variants (and numbers of hits) but does not change any of the distributions.

We first tested these predictions for three height, platelet crit (a blood phenotype), and FEV1 (a measure of lung function), where we reduced the sample size to 330,000 so it is identical for all traits (Sect 6.7 in S1 Note). Fig 7A shows height and platelet crit. These two traits differ greatly in their number of hits (1533 vs 648), estimated heritability (77% vs 32%) and mutational target size (30 MB vs 10 MB). However, we estimate that they have very similar values of h^2^/L (∼3 × 10^−8^). Consistent with our prediction, the distributions of z-scores, effect sizes, and MAFs for significant hits are nearly identical for the two traits.

In contrast, height and FEV1 have similar mutational target sizes, but h^2^/L for height is about 6-fold higher than for FEV1 (Fig 7B). As expected, this results in greater power to detect variants associated with height than with FEV1 (40% vs. 14% of h^2^ explained by significant hits). Consequently, genome-wide significant hits for height have larger effect sizes and larger z-scores than for FEV1. We also see slightly fewer low-MAF hits for FEV1 (p = 0.01, KS test), due to the lower power compared to height.

We predicted that after rescaling the z-scores (and effect sizes) for these traits by h2/L, their architectures should become nearly identical. On this scale, the significance cutoff is 5.45/h2/L, which is higher for FEV1 than for height (illustrated by the light and dark blue regions in Fig 6C). We therefore compared the architecture of genome-wide significant hits for both traits using the higher threshold in the scaled units. After doing so, we only have 95 hits for height, but it is apparent that the summary properties for both traits are indeed highly similar (Fig 7C).

We repeated this scaling procedure for the 57 other traits in our dataset whose scaled threshold is below that of FEV1. After doing so, the distributions for all hits are highly similar to the scaled distributions for height and FEV1 (dotted line, Fig 7C). We repeated similar analyses for all traits, and found that after applying these scaling laws the architectures of all 95 traits are nearly identical (Sect 10 in S1 Note).

Discussion

What determines the genetic architecture of complex traits? Does genetic variation in a given trait primarily reflect the idiosyncrasies of its biology–or alternatively, does it reflect processes that are shared among different traits?

Here we describe a principled approach to tackle these questions. Our point of departure is an evolutionary model of genetic architecture based on empirically motivated and interpretable biological assumptions. Given arguments and evidence that many quantitative (continuous) traits are subject to stabilizing selection [23,37,38], and that genetic variation affecting one trait often affects many others [5,20,21], we model selection on alleles that arises from stabilizing selection in multi-dimensional trait space. Otherwise, we assume the standard population genetic model incorporating the effects of mutation, selection, genetic drift and demographic history. This model gives rise to a family of genetic architectures, where the architecture of a given trait is determined by the distribution of selection coefficients at trait-affecting sites as well as two scaling the mutational target size and the heritability per site.

We performed inference allowing all these parameters to vary among traits. We found that the model provides a good fit to the joint distribution of allele frequencies and effect sizes at genome-wide significant hits for the 95 quantitative traits in our dataset. Intriguingly, we also found that the confidence envelopes for the distribution of selection coefficients largely overlap among these traits, suggesting that their genetic architecture can be approximated by a Single Shared Distribution (SSD) of selection coefficients and the two additional scaling parameters per trait. We then estimated the SSD using the data from all traits jointly, showing that it fits the architecture of individual traits well, and validated our estimate of the SSD by showing that it accurately predicts the distribution of allele ages at GWAS hits for all traits.

The fit of the SSD model implies that, aside from the two scaling parameters, the genetic architecture of genome-wide significant hits is similar among all 95 traits. Cross-sections of the estimated shared architecture are visualized in Fig 6. Indeed, as predicted, after we rescale the effect sizes by the estimated heritability per site, we find that the joint distributions of effect sizes and allele frequencies for GWAS hits are similar among traits. Meanwhile the number of hits for a given trait is proportional to estimated target size (Fig 7).

These findings delineate the attributes of genetic architecture that are shared among traits and those that are trait-specific. In doing so, they raise new insights and

How similar are the distributions of selection coefficients underlying heritable variation in different complex traits? Our results indicate that the distributions of selection effects among the 95 traits we examined were similar enough for the architectures of GWAS hits to be insensitive to the differences among them. Previous work has hinted at this for example, work using the α-model reported broadly similar relationships between MAF and effect size among a variety of traits [25,26,36], though it should be noted that there is no straightforward interpretation of α in terms of selection coefficients [36]. Previous work also hinted at the shared features of the distributions that we found, notably about selection on many GWAS hits being weak enough for them to precede the Out-of-Africa bottleneck and be shared among populations [52,53], and about their effect sizes spanning several orders of magnitude [18,19].

We might expect some similarity in the distribution of selection coefficients among complex traits. In a hypothetical extreme in which all functional variation affects all traits, the distribution would necessarily be shared among traits. This logic may well explain similarities among traits whose mutational target sizes encompass much of the functional portion of the genome; for example, we estimate the target size for BMI at ∼7% of the genome compared to ∼8% estimated to be functional [55–57]. It could also apply to traits that are mediated through the same tissues [14]. However, this logic does not explain the similarity among traits whose target sizes are substantially smaller and are primarily mediated through different tissues [58]. The biomarkers in our dataset, for example, have target sizes that are more than an order of magnitude smaller than BMI (e.g., calcium level with a target size of ∼0.1% of the genome), and are mediated through distinct cell types or tissues [9].

Shared features in the distribution of selection coefficients may reflect similarity in the biological systems in which variation affecting complex traits arises, notably in gene-regulatory networks. Heritable variation in complex traits is spread across most of the genome and is enriched in regulatory regions near most genes that are expressed in the tissues that affect these traits [14]. The lead variants typically explain a tiny fraction of the heritability [3], and the most relevant biological pathways are usually only modestly enriched for heritability [9,14]. In other words, most heritable variation is mediated through the regulation of genes and pathways that are not closely connected to the trait’s biology [9,14,59]. Perhaps, the essential logic of gene regulatory networks and their evolution are sufficiently similar across tissues to generate the shared features we found, even if the specific pathways, genes, enhancers, and variants differ.

The similarity among distributions has important qualifications, however. Even among the 95 traits we examined and with our inference restricted to GWAS hits, we still identified statistically significant differences among the distributions of selection effects for different traits (see above and Sect 4 in S1 Note). As we already noted, restricting our inference to GWAS hits with MAF > 1%, delimits our power to infer selection effects at the weak and strong ends of the range of selection coefficients (e.g., for s≲10−6 and ≳10−2). We therefore cannot rule out more substantial differences in the distribution of, e.g., strongly selected, rare variants, even among the traits we examined (see, e.g., [28]). Additionally, the traits we examined were ascertained based on having ≥ 100 hits. Traits with vastly different distributions of selection coefficients may fail to pass this threshold, e.g. traits for which most heritability arises from strongly selected variants.

Limitations and future analyses. These limitations and others can be addressed in the future. For example, whole genome sequencing [60] would enable greater power for strongly selected variants. Meanwhile, methodologies that integrate over the full distribution of causal variants accounting for LD [11,19,36] may allow us to relax the reliance on genome-wide significant hits, thus increasing the power to identify weakly selected variants, more generally increasing the precision of our inference for both the TSD and SSD models, and allowing for the analysis of traits with fewer GWAS hits.

Future increases in sensitivity and applications of our inference to different kinds of traits may also warrant extensions of our evolutionary model. Notably, our current model is restricted to quantitative traits and is not immediately applicable to binary traits including diseases. One could imagine that variants that affect the risk of some complex diseases are selected to minimize that risk and are thus subject to purifying (directional) selection, but recent studies indicate that, at least common variation affecting them is predominated by the kind of pleiotropic stabilizing selection we assumed here [61,62]. With an appropriate adjustment of the model to fit data from case-control GWAS, we may find that complex diseases have similar architectures to those of the quantitative traits that we examined.

Our model and inference rely on the simplifying assumption that all the genetic variation in a given trait is affected by the same degree of pleiotropic selection (this degree was reflected in the dimension of the trait space; see, e.g., Eq 3 and Simons et al. 2018). We now have strong evidence that this assumption is violated in ways that affect the kinds of variants that are identified in GWAS [50]. Future extensions of the model may incorporate variation in the degree of pleiotropic selection, where this variation may also differ among complex traits.

Why do trait-specific scaling parameters vary? Future refinements notwithstanding, our evolutionary model with a shared distribution of selection effects fits the data from all 95 traits in our dataset remarkably well, indicating, as we have also confirmed, that differences in genetic architecture of genome-wide significant hits among highly polygenic traits are largely determined by the two trait-specific scaling the mutational target size and heritability per site.

The mutational target size varies over 2 orders of magnitude among our 95 traits. While the estimates are novel, the variation among them is hardly surprising given the vast differences in the biology of these traits. These traits vary in being affected by few to many tissues and by the number and properties of genes and pathways that directly affect the trait in these tissues [59,63]. Moreover, the pathways that directly affect these traits plausibly differ in how buffered they are against genetic (and environmental) variation and in their modularity, plausibly reflective of the kinds of traits and of selection pressures over longer evolutionary timescales than the turnover time of heritable variation [64,65]. Our estimates of heritability are less variable, but still span an order of magnitude, and vary among different kinds of traits (S1 Table). Both of these observations have been known for almost a century, and yet the question about their causes remains largely open [23,37,38]. Our results indicate that other differences in architecture among traits are dwarfed by the variation caused by these two scaling parameters.

Outlook. Taking a step back, our results highlight that evolutionary thinking is essential to understanding of the findings emerging from human GWASs, and more generally, heritable variation in complex traits. This insight is consistent with long-standing thinking in the field, given that heritable variation in complex traits reflects the outcome of numerous genetic perturbations to the phenotype that passed through the sieve of evolution [23,38]. Specifically, the signal measured in GWASs reflects causal variants’ contribution to heritability, which depends on their effect on the trait under consideration, but also on their minor allele frequency, where the relationship between the two is mediated by natural selection. What is surprising, at least to us, is how far an analysis based on evolutionary modeling can go, in this case, showing that the genetic architectures of many highly polygenic quantitative traits are largely shared. This finding carries many implications about human GWASs and their applications, some of which we plan to explore elsewhere. Alongside other evidence, it also hints at underlying biology that largely remains to be discovered, plausibly relating to properties of gene regulatory networks. We think that a combination of evolutionary reasoning alongside a systems approach to gene regulation would move us closer to answering questions tracing back to the very beginning of the field of genetics about the mapping from genetic variation to phenotypes.

Supporting information