Linkage Disequilibrium and Heterozygosity Modulate the Genetic Architecture of Human Complex Phenotypes

We propose an extended Gaussian mixture model for the distribution of causal effects of common single nucleotide polymorphisms (SNPs) for human complex phenotypes that depends on linkage disequilibrium (LD) and heterozygosity (H), while also allowing for independent components for small and large effects. Using a precise methodology showing how genome-wide association studies (GWAS) summary statistics (z-scores) arise through LD with underlying causal SNPs, we applied the model to GWAS of multiple human phenotypes. Our findings indicated that causal effects are distributed with dependence on total LD and H, whereby SNPs with lower total LD and H are more likely to be causal with larger effects; this dependence is consistent with models of the influence of negative pressure from natural selection. Compared with the basic Gaussian mixture model it is built on, the extended model – primarily through quantification of selection pressure – reproduces with greater accuracy the empirical distributions of z-scores, thus providing better estimates of genetic quantities, such as polygenicity and heritability, that arise from the distribution of causal effects.


INTRODUCTION
There is currently great interest in the distribution of causal effects among trait-associated single nucleotide polymorphisms (SNPs), and recent analyses of genome-wide association studies (GWAS) have begun uncovering deeper layers of complexity in the genetic architecture of complex human traits and disorders (Gazal et al., 2017;Zeng et al., 2018;Zhang et al., 2018;Frei et al., 2019;Holland et al., 2020). This research is facilitated by using new analytic approaches to interrogate structural features in the genome and their relationship to phenotypic expression. Some of these analyses take into account the fact that different classes of SNPs have different characteristics and play a multitude of roles (Schork et al., 2013;Finucane et al., 2015;. Along with different causal roles for SNPs, which in itself would suggest differences in distributions of effect-sizes for different categories of causal SNPs, the effects of minor allele frequency (MAF) of the causal SNPs and their total correlation with neighboring SNPs are providing new insights into the action of selection on the genetic architecture of complex traits (Gazal et al., 2017;Wray et al., 2018;Zhang et al., 2018).
Any given mutation is likely to be neutral or deleterious to fitness (Fay et al., 2001). Natural selection partly determines how the prevalence of a variant develops over time in a population, and evidence for its action can be found in the relationship between effect size and MAF in complex traits and common diseases (Zeng et al., 2018). Negative selection acts predominantly to keep variants deleterious to fitness at low frequency, or ultimately remove them. The larger the effect of a deleterious variant the more efficient negative selection will be, suggesting that the lower the MAF the larger the effect size -and this is expected under evolutionary models (Pritchard and Cox, 2002;Eyre-Walker, 2010), and consistent with empirical findings (Park et al., 2011) and recent analyses based on genome-wide associations (Zeng et al., 2018;Schoech et al., 2019;O'Connor et al., 2019).
The effect of linkage disequilibrium (LD) has also been studied, suggesting that SNPs with low "levels of LD" (LLD) explain more heritability, which is again consistent with the action of negative selection (Gazal et al., 2017). One unexplored issue is how the prior probability of a SNP being causal depends on its LD score (or related measures). Due to the complexity of genetic forces acting on alleles, it is not clear what form any such dependency might take. However, explicitly modeling any such role promises to yield a closer match between empirical distributions of GWAS summary statistics and model predictions, and thereby can provide more accurate estimates of quantities of interest like polygenicity and heritability.
In previous work (Holland et al., 2020), building on earlier reports of others (e.g., George and McCulloch (1993); Erbe et al. (2012); Zhou et al. (2013)), we presented a basic Gaussian Gaussian mixture model to describe the distribution of underlying causal SNP effects (the per unit allelic effect, β, estimated from simple linear regression). Due to extensive and complex patterns of LD among SNPs, many non-causal SNPs will exhibit strong association with phenotypes, resulting in a far more complicated distribution for the summary z-scores. The basic model for the distribution of the causal β's is a mixture of non-null and null normal distributions, the latter (denoted N (0, 0)) being just a delta function (or point mass at zero): where π 1 is the polygenicity, i.e., the proportion of SNPs that are causal (equivalently, the prior probability that any particular SNP is causal), and σ 2 β is the discoverability (the expectation of the square of the effect size), which was taken to be a constant across all causal SNPs. The distribution of z-scores arising from this shows strong heterozygosity (H=MAF×(1-MAF)) and total LD (TLD) de-pendence; TLD, defined in the subsection 1.3, is similar to LD score  but takes into account more neighboring SNPs.
The recent work by others focusing on selection pressure using a single causal Gaussian (Zeng et al., 2018;Schoech et al., 2019) indicated that it is important to take SNP heterozygosity into account as a component in discoverability (initially explored by examining four fixed relationships in (Speed et al., 2012); see also (Lee et al., 2013)). It was also recently shown that an additional Gaussian distribution for the β's might be appropriate if large and small effects are distributed differently (Zhang et al., 2018). These approaches, however, were not combined -a single causal Gaussian incorporating heterozygosity, versus two causal Gaussians with no heterozygosity dependence. An extra complexity is that measures related to total LD, such as LLD (Gazal et al., 2017), can be expected to play an important role in the distribution of causal effects. It is unclear how all these factors impact each other. A final and very important matter is that many analyses in the literature by other groups were conducted in the "infinitesimal" model framework Gazal et al., 2017;Schoech et al., 2019;Speed and Balding, 2019), where all SNPs are causal, though there is compelling evidence that only a very small fraction of SNPs are in fact causal for any given phenotype (Zhu and Stephens, 2017;Zeng et al., 2018;Zhang et al., 2018;Frei et al., 2019;Holland et al., 2020).
In the current work, we sought to extend our earlier work to incorporate multiple Gaussians, while taking into account TLD as a factor in polygenicity and selection effects reflected in heterozygosity, in modeling the distribution of causal β's. With a wide range of model parameter values across real phenotypes, the specificity of the individual parameters for a given phenotype make them more narrowly defining of the distribution of summary statistics for that phenotype.

Materials and Methods
1.1. The model: an extension of prior work The methodology calls for using an extensive reference panel that likely involves all possible causal SNPs with MAF greater than some threshold (e.g., 1%), and regard the z-scores for typed or imputed SNPs -a subset of the reference SNPs -as arising, directly or through LD, from the underlying causal SNPs in the reference panel.
Since the single causal Gaussian, Eq. 1, has provided an appropriate starting point for many phenotypes, it is reasonable to build from it. With additional terms included, if it turns out that this original term is not needed, the fitting procedure, if implemented correctly, should eliminate it. Also, anticipating extra terms in the distribution of causal β's, we introduce a slight change in labeling the Gaussian variance (σ 2 β → σ 2 b ), and write the distributions for the causal component only -it being understood that the full distribution will include the last term on the right side of Eq. 1 for the prior probability of being null.
Given that for some phenotypes there is strong evidence that rarer SNPs have larger effects, we next include a term that reflects this: a Gaussian whose variance is proportional to H S , where H is the SNP's heterozygosity and S is a new parameter which, if negative, will reflect the noted behavior (Zeng et al., 2018). With the addition of the new term, the total prior probability for the SNP to be causal is still given by π 1 . Thus, extending Eq. 1, we get: where p c (0 p c 1) is the prior probability that the SNP's causal component comes from the "c" Gaussian (with variance σ 2 c H S ), and p b ≡ 1−p c is the prior probability that the SNP's causal component comes from the "b" Gaussian (with variance σ 2 b ). This extension introduces an extra three parameters p c , σ c , and S, assumed for the moment to be the same for all SNPs. Ignoring inflation and implementation details like choice of reference panel and parameter estimation scheme, setting p c ≡ 1 recovers the model distribution assumed in (Zeng et al., 2018); additionally setting π 1 ≡ 1 recovers the primary model distribution assumed in ; and additionally setting S ≡ −1 recovers the model assumed in , while instead setting S ≡ −0.25 partially recovers the "recommended" LD-adjusted kinships (LDAK) model distribution in (Speed et al., 2017;Speed and Balding, 2019). For a further discussion of LDAK  Table G1 for parameter values). These functions can be summarized by three quantities: the maximum value, p c1 , which occurs at L = 1; the total LD value, L = mc, where pc(mc) = p c1 /2, give by the gray dashed lines in the figure; and the total LD width of the transition region, wc, defined as the distance between where pc(L) falls to 95% and 5% of p c1 given by the flanking red dashed lines in the figure. Numerical values of p c1 , mc, and wc are given in Table 1  Within the infinitesimal models (Speed et al., 2012;Finucane et al., 2015;Speed et al., 2017;Gazal et al., 2017;Schoech et al., 2019;Speed and Balding, 2019), it is not clear the degree to which explicit LD-dependence in the variance of the causal effect size merely takes into account the effect on z-scores due to LD with causal SNPs, and how much it models any true effect of LD on underlying causal effect size. Also, such models preclude examining if the TLD of a variant has any bearing on whether the variant is causal. In contrast is the "BayesS" model (Zeng et al., 2018), a causal mixture model (i.e., not infinitesimal), using individual genotype data and a reference panel of ∼484k non-imputed SNPs, that examines the effects of heterozygosity on effect size. The model we present here is in some respects an extension of that, but based on summary statistics, adding TLD dependence and an additional Gaussian, and we fit the model from a reference panel of 11 million SNPs using the exact procedure -convolutionto relate posited underlying distributions of causal effects to empirical distributions of z-scores.
Along with heterozygosity, SNP effect size might independently depend on TLD (for which we use the variable L in equations below), and in principle this could be explored in a manner similar to how heterozygosity is incorporated in the "c" causal effect Gaussian variance (e.g., with an extra factor L T , say, scaling σ 2 c , where T would be a new parameter). However, TLD and heterozygosity are often related (given TLD, the expected heterozygosity shows a distinct well-defined pattern for SNPs with TLD<200, i.e., about 80% of SNPs -see Supplemental Material, Figure  S6), and independent contributions might be difficult to disentangle. Instead, here we explore a separate mathematical role for TLD.
There is no obvious a priori reason why the probability (in a Bayesian sense) of a SNP's being causally associated with any particular trait should be independent of the SNP's TLD. Indeed, with the complex interaction of multiple genetic forces such as mutation, genetic drift, and selection, the net relationship between TLD -through mechanisms like background selection -and causal association with a particular phenotype is not clear. The results of Gazal et al. (2017), however, indicate that SNPs with low LLD have significantly larger per-SNP heritability. Note that in Eq. 2, for σ 2 c H S > σ 2 b the "c" Gaussian will describe larger effect. In this case, the LLD dependence suggests modulating the "c" Gaussian such that p c is larger for smaller TLD. Whatever the relationship, however, the more accurately it is incorporated in a model for the distribution of effect sizes should lead to more accurate reproduction of the distribution of empirical summary statistics and estimation of quantities of interest like polygenicity, heritability, and selection effects.
As heterozygosity decreases, SNPs will continue to have a range of total LD (see Supplemental Material, Figure S7 for the number density of SNPs with respect to heterozygosity and TLD). We explore here the possibility that the prior probability of being causal with large effect decreases with TLD. If the "c" Gaussian is capturing larger effects from rarer SNPs, reflecting selection pressure, we inquire if the prior probability for a causal SNP's contribution from the "c" Gaussian is TLD-mediated. Specifically, instead of treating p c as a constant, we explore the possibility that it is larger for SNPs with lower TLD; in the event that this probability is in fact constant, or increasing, with respect to TLD, we would at least not expect to find it decreasing. This can be accomplished by means of a generalized sigmoidal function that will have a maximum at very low TLD, might maintain that maximum for all SNPs (equivalently, p c is a constant), or decrease in amplitude slowly or rapidly, possibly to 0, for SNPs with higher TLD. With the variable L denoting the TLD of a SNP, such a function of TLD can be characterized by three parameters: its amplitude (at L = 1), the TLD at the midpoint of the sigmoidal transition, and the width of the sigmoidal transition (over a wide or narrow range of TLD). We use the following general form with three parameters, y max , x mid , and x width : defined in the range −∞ < x < ∞, for which y(x) is monotonically decreasing and bounded 0 < y(x) < y max , with 0 y max 1; −∞ < x mid < ∞ locates the mid point of the overall sigmoidal transition (y(x mid ) = y max /2), and 0 < x width < ∞ controls its width (y(x) smoothly changing from a Heaviside step function at x mid as x width → 0 to a constant function as x width → ∞). Examples are shown in Figure 1 (scaled by π 1 , giving the physically interesting total prior probability of a SNP being causal with respect to the selection "c" Gaussian, as a function of the SNP's TLD); parameter values are in 3 Table G1. Mathematically, the curve can continue into the "negative TLD" range, revealing a familiar full sigmoidal shape; since we are interested in the range 1 x max(TLD), below we report the actual mid point (denoted m c ) and width (w c , defined below) of the transition that occurs in this range. Then where p c (L) is the sigmoidal function (0 p c (L) 1 for all L) given by y(x) in Eq. 3 for L = x 1, which numerically can be found by fitting for its three characteristic parameters.
As a final possible extension, we add an extra term -a "d" Gaussian -to describe larger effects not well captured by the "b" and "c" Gaussians. This gives finally: where σ 2 d is a new parameter, p d (L) is another general sigmoid function (0 p d (L) 1 for all L) where now there is the added constraint 0 p c (L) + p d (L) 1, and the prior probability for the "b" Gaussian becomes Depending on the phenotype and the GWAS sample size, it might not be feasible, or meaningful, to implement the full model. In particular, for low sample size and/or low discoverability, the "b" Gaussian is all that can be estimated, but in most cases both the "b"and "c" Gaussians can be estimated, and β will be well characterized by Eq. 4. We refer to the model given by Eq. 1 as model B; models C and D are given by Eqs. 4 and 5, respectively.
As we described in our previous work (Holland et al., 2020), a z-score is given by a sum of random variables, so the a posteriori pdf (given the SNP's heterozygosity and LD structure, and the phenotype's model parameters) for such a composite random variable is given by the convolution of the pdfs for the component random variables. This facilitates an essentially exact calculation of the z-score's a posteriori pdf arising from the underlying model of causal effects as specified by Eq. 1, 4, or 5.
For our reference panel we used the 1000 Genomes phase 3 data set for 503 subjects/samples of European ancestry Sveinbjornsson et al., 2016). In ref. (Holland et al., 2020) we describe how we set up the general framework for analysis, including the use of the reference panel and dividing the reference SNPs into a sufficiently fine 10×10 heterozygosity×TLD grid to facilitate computations.
In our earlier work on the "b" model we gave an expression, denoted G(k), for the Fourier transform of the genetic contribution to a z-score, where k is the running Fourier parameter. The extra complexity in the "c" and "d" models here requires a modification only in this term, which we describe below in the 1.4 subsection. In addition to the parameters presented above, we also include an inflation parameter σ 0 : if z u denotes uninflated GWAS z-scores and z denotes actual GWAS z-scores, then σ 0 is defined by z = σ 0 z u (Devlin and Roeder, 1999). The optimal model parameters for a particular phenotype are estimated by minimizing the negative of the log likelihood of the data (z-scores) as a function of the parameters. This is done with the "b" model as before, and then proceeding iteratively with the more complex models, continually re-estimating the new values of the earlier parameters that maximize the likelihood when a new parameter is introduced, with extensive single and multiple parameter searches (to avoid being trapped in local minima) until all parameters simultaneously are at a convex minimum. This involved many explicit and repeated coarse-and finegrained linear (for each individual parameter) and grid (for multiple parameters simultaneously) searches as well as many Nelder-Mead multidimensional unconstrained nonlinear minimizations. There is no artificial constraining on the parameters; for example, no prior assumption is made about the relative sizes of the causal σ's, or the parameter values of the TLD-dependent prior probabilities (the full range of amplitudes, transition widths, and location of the mid-point of the transitions are searched).
The total SNP heritability is given by the sum of heritability contributions of each SNP in the reference panel from each of the relevant Gaussians. In the Bayesian approach, we do not know the value of the causal effect of any particular SNP, but we assume it comes from a distribution, β(H, L), which characterize our ignorance of it. For a specific Gaussian ("b", "c", or "d") in our model, the contribution of the SNP to heritability is given by the prior probability that the SNP is causal with respect to that Gaussian, times the expected value of the square of the effect size, E(β 2 ), times H. For the "c" Gaussian, for example, the prior probability that the SNP is causal is π 1 p c (L), and E(β 2 ) is just the variance, σ 2 c H S . Thus, the contribution of this SNP to the overall heritability associated with the "c" Gaussian, h 2 c , is π 1 p c (L)Hσ 2 c H S . Below we report the sums over all such contributions. The number of causal SNPs associated with the "c" Gaussian, n c , is given by summing π 1 p c (L) for each reference panel SNP, and similarly for the other Gaussians. All heritabilities and discoverabilities are, as before, corrected with respect to the inflation parameter, i.e., divided by σ 2 0 . 4  Figure 2: QQ plots of (pruned) z-scores for qualitative phenotypes (dark blue, 95% confidence interval in light blue) with model prediction (yellow). See Supplemental Material, Figures S15 to S22. The value given for p c1 is the amplitude of the full pc(L) function, which occurs at L = 1; the values (mc, wc) in parentheses following it are the total LD (mc) where the function falls to half its amplitude (the middle gray dashed lines in Figure 1 are examples), and the total LD width (wc) of the transition region (distance between flanking red dashed lines in Figure 1). Similarly for p d1 (m d , w d ), where given. h 2 b , h 2 c , and h 2 d are the heritabilities associated with the "b", "c", and "d" Gaussians, respectively. h 2 is the total SNP heritability, reexpressed as h 2 l on the liability scale for binary phenotypes. Parameter values are also given in Table 1 and heritabilities are also in Table 3; numbers of causal SNPs are in Table 2. Reading the plots: on the vertical axis, choose a p-value threshold for typed or imputed SNPs (SNPs with z-scores; more extreme values are further from the origin), then the horizontal axis gives the proportion, q, of typed SNPs exceeding that threshold (higher proportions are closer to the origin). See also Supplemental Material, Figure S1, where the y-axis is restricted to 0 −log 10 (p) 10. The All code used in the analyses, including simulations, is publicly available on GitHub (Holland, 2019b,a).  (Willer et al., 2013). Most participants were of European ancestry. (A spreadsheet giving data sources is provided in the Supplemental Material.) In the tables, we also report results for body mass index (GIANT 2015) (N = 233,554) (Locke et al., 2015), and height (UKB-GIANT 2018) .

Data Preparation
In Figure 2 and Supplemental Material, Figure S1 we report effective sample sizes, N ef f , for the case-control GWASs. This is defined as N ef f = 4/(1/N cases +1/N controls ), so that when N cases = N controls , N ef f = N cases +N controls = N , the total sample size, allowing for a straightforward comparison with quantitative traits.
Confidence intervals for parameters were estimated using the inverse of the observed Fisher information matrix (FIM). The full FIM was estimated for up to eight parameters used in model C, and for the remaining parameters that extend the analysis to model D the confidence intervals were approximated ignoring off-diagonal elements.

Intelligence Education
Height ( Figure S2, where the y-axis is restricted to 0 −log 10 (p) 10. For (A) BMI, see also Supplemental Material, Figure S37.
Additionally, the w d parameter was treated as fixed quantity, the lowest value allowing for a smooth transition of the p d (L) function to 0 (see Supplemental Material, Figure S5; for CD, UC, and TC, however, the function p d (L) was a constant (=p d (1)). For the derived quantities h 2 and n causal , which depend on multiple parameters, the covariances among the parameters, given by the off-diagonal elements of the inverse of the FIM, were incorporated. Numerical values are in 3 Tables E1-E3, and 3 Tables F1-F2.
In order to carry out realistic simulations (i.e., with realistic heterozygosity and LD structures for SNPs), we used HAPGEN2 (Li and Stephens, 2003;Spencer et al., 2009;Su et al., 2011) to generate genotypes for 10 5 samples; we calculated SNP MAF and LD structure from 1000 simulated samples.

Total Linkage Disequilibrium
Sequentially moving through each chromosome in contiguous blocks of 5,000 SNPs in the reference panel, for each SNP in the block we calculated its Pearson r 2 correlation coefficients (that arise from linkage disequilibrium, LD) with all SNPs in the central block itself and with all SNPs in the pair of flanking blocks of size up to 25,000 each. For each SNP we calculated its total linkage disequilibrium (TLD), given by the sum of LD r 2 's thresholded such that if r 2 < r 2 min we set that r 2 to zero (r 2 min = 0.05). The fixed window size corresponds on average to a window of ±8 centimorgans. This is deliberately larger than the 1-centimorgan window used to define LD Score , because the latter appears to exclude a noticeable part of the LD structure. In applying the model to summary statistics, we calculated histograms of TLD (using 100 bins) and ignoring SNPs whose TLD was so large that their frequency was less than a hundredth of the respective histogram peak; typically this amounted to restricting to SNPs for which TLD 600. We also ignored summary statistics of SNPs for which MAF 0.01.
Since we are estimating a dozen or fewer parameters from millions of data points, it is reasonable to coarsegrain the data. Knowing the TLD and heterozygosity (H) of each SNP, we divided the full set of GWAS SNPs into a H×TLD coarse-grained grid; we found that 10×10 is more than sufficient for converged results.

Model PDF
When implementing the discrete Fourier transform (DFT) to calculate model a posteriori probabilities for z-score outcomes for a single SNP, we discretize the range of possible z-scores into the ordered set of n (equal to a power of 2) values z 1 , . . . , z n with equal spacing between neighbors given by ∆z (z n = −z 1 − ∆z, and z n/2+1 = 0). Taking z 1 = −38 allows for the minimum p-values of 5.8 × 10 −316 (near the numerical limit); with n = 2 10 , ∆z = 0.0742. Given ∆z, the Nyquist critical frequency is f c = 1 2∆z , so we consider the Fourier transform function for the z-score pdf at n discrete values k 1 , ..., k n , with equal spacing between neighbors given by ∆k, where k 1 = −f c (k n = −k 1 − ∆k, and k n/2+1 = 0; the DFT pair ∆z and ∆k are related by ∆z∆k = 1/n).
In our earlier work on the "b" model (Holland et al., 2020) we gave an expression, denoted G(k j ), for the Fourier transform of the genetic contribution to a z-score, where k j is the discrete Fourier variable described above. In constructing the a posteriori pdf for a z-score, the extra complexity in the "c" and "d" models presented in the current work requires a modification only in this term.
Let H denote the LD and heterozygosity structure of a particular SNP, a shorthand for the set of values {n w , H w , L w : w = 1, . . . , w max } that characterize the SNP, and let M denote the set of model parameters for whichever model , denote the vector of products of Fourier transform values, and let F −1 denote the inverse Fourier transform operator. Then for the SNP in question, the vector of pdf values, pdf z , for the uniformly discretized possible z-score outcomes z 1 , . . . , z n described above, i.e., pdf z = (f 1 , ..., f n ) where f i ≡ pdf(z i |H, M, N ), is Thus, the i th element pdf zi = f i is the a posteriori probability of obtaining a z-score value z i for the SNP, given the SNP's LD and heterozygosity structure, the model parameters, and the sample size.

Data Availability
File S2 contains detailed descriptions of all GWAS data used, all of which is publicly available.

Phenotypes
Summary QQ plots for pruned z-scores are shown in Figure 2 for seven binary phenotypes (for AD we separate out chromosome 19, which contains the APOE gene), and Figure 3 for seven quantitative phenotypes (including two separate GWAS for height), with model parameter values in Table 1. An example of the breakdowns of a summary plot with respect to a 4×4 grid of heterozygosity×TLD (each grid a subset of a 10×10 grid) is in Figure 4 Table 2: Numbers of causal SNPs. n causal is the total number of causal SNPs (from the 11 million in the reference panel); n b , nc, and n d are the numbers associated with the "b", "c", and "d" Gaussians, respectively. 95% confidence intervals are in 3, Table F1.
Material, Figures S15-S29. For each phenotype, model selection (B, C, or D) was performed by testing the Bayesian information criterion (BIC) -see 3 Table D1. For comparison, all QQ figures include the basic (B) model in green; the extended model C or D, consistently demonstrating improved fits, in yellow; and the data in blue. The distributions of z-scores for different phenotypes are quite varied. Nevertheless, for most phenotypes analyzed here, we find evidence for larger and smaller effects being distributed differently, with strong dependence on total LD, L, and heterozygosity, H.
For σ 2 b ≪ σ 2 c H S , and so focusing on the "c" Gaussian (and "d" Gaussian if applicable for σ 2 b ≪ σ 2 d ), our model estimates an effective polygenicity as a one-dimensional function of L. We find that polygenicity is dominated by SNPs with low L. However, the degree of restriction varies widely across phenotypes, depending on the shapes and sizes of p c (L) and p d (L) in Eq. 5, the prior probabilities that a causal SNP belongs to the "c" and "d" Gaussians. These prior probabilities, multiplied by π 1 , are shown in Figure 1 and Supplemental Material, Figures S3-S5. Taking into account the underlying distribution of reference SNPs with respect to heterozygosity, these distributions lead to a varied pattern across phenotypes of the expected number of causal SNPs in equally-spaced elements in a two-dimensional H×TLD grid, as shown for height (2014) in Figure 5 C, and for all phenotypes in Supplemental Material, Figures S11-S14 (third columns). Further, for any given phenotype, the effect sizes of causal variants come from distributions whose variances can be widely differ-ent -by up to two orders of magnitude. Thus, given the prior probabilities (p b , p c , and p d ) by which these distributions are modulated as a function of L, we are able to estimate the expected effect size per causal-SNP, E(β 2 ), in each H×TLD grid element, as shown in Figure 5 D and Supplemental Material, Figures S11-S14 (fourth columns). In general, SNPs with lower L have larger E(β 2 ). However, the selection parameter S in the "c" Gaussian has a large impact on E(β 2 ) as a function of H (see Supplemental Material, Figure S9). As a result, for most phenotypes, we find that the effect sizes for low MAF causal SNPs (H < 0.05) are several times larger than for more common causal SNPs (H > 0.1). We find that heritability per causal-SNP is larger for lower L, a general pattern that follows from at least one of the prior probabilities, p c (L) or p d (L), being non-constant. However, because heritability per causal-SNP is proportional to H, we find that, even with negative selection parameter, S (and thus larger E(β 2 ) for lower H), the heritability per causal-SNP is largest for the most common causal SNPs (H > 0.45).
SNP heritability estimates from the extended model, as shown in Table 3, are uniformly larger than for the exclusive basic model. With the exception of BMI, generally a major portion of SNP heritability was found to be associated with the "selection" component, i.e., the "c" Gaussian, with a smaller contribution from the "d" Gaussian -see the right-most column in Table 3. In 3 Table B2 is a comparison of basic model estimates of the number of  , h 2 c , and h 2 d are the heritabilities associated with the "b", "c", and "d" Gaussians, respectively. The last column, labeled %c, gives the percentage of SNP heritability that comes from the "c" Gaussian. 95% confidence intervals are in 3, Table F2. A comparison of the total heritabilities with estimates from our basic model alone (Holland et al., 2020), and with estimates from other models (Zhang et al., 2018;Zeng et al., 2018;Schoech et al., 2019) are given in 3, Table B1.
causal SNPs and heritability with the corresponding net contributions n c + n d and h 2 c + h 2 d from the large effects components of the extended model. The basic model can be seen as a good approxiation to the large effects components of the extended model. Correspondingly, the "b" Gaussian as a components of the extended model now represents a relatively large number of much weaker effects.

Simulations
To test the specificity of the model for each real phenotype, we constructed simulations where, in each case, the true causal β's (a single vector instantiation) for all reference panel SNPs were drawn from the overall distribution defined by the real phenotype's parameters (thus being the "true" simulation parameters). We set up simulated phenotypes for 100,000 samples by adding noise to the genetic component of the simulated phenotype, and performed a GWAS to calculate z-scores. We then sought to determine whether the true parameters, and the component heritabilities, could reasonably be estimated by our model. In Figures 6 and 7 we show the results for the simulated case-control and quantitative phenotypes, respectively. Overall heritabilities were generally faithful to the true values (the values estimated for the real phenotypes) -see 3, Table C2 -though for Crohn's disease the simu- lated value was overestimated due to the h d component.
Note that for the case-control simulated phenotypes, the heritabilities on the observed scale, denotedĥ 2 in Figure  6, should be compared with the corresponding values in Figure 2, not withĥ 2 l , which denotes heritability on the liability scale, i.e., adjusted for population prevalence; note also that for case-control phenotypes, we implicitly assume the same proportion P of cases in the real and simulated GWAS (for the basic model, assuming the liability-scale model is correct, one can easily see from the definition of h 2 l that discoverability σ 2 β ∝ P (1 − P ); this carries over to σ b , σ c , and σ d used here). Polygenicities and discoverabilities were also generally faithfully reproduced. However, for ALS restricted to chromosome 9, and BMI, the selection parameter was incorrectly estimated, owing to the weak signal in these GWAS (e.g., for BMI, π 1 p c (1) ≃ 8 × 10 −6 , Fig. S4 (A), and only around 5% of SNP heritability was found to be associated with the "c" Gaussian, Table 3) and very low polygenicity (small number of causal SNPs) for the "c" Gaussian. Given the wide variety and even extreme ranges in parameters and heritability components across diverse simulated phenotypes, the multiple simulated examples provide checks with respect to eachother for correctly discriminating phenotypes by means of their model parameter estimates: the results for individual cases are remarkably faithful to the respective true values, demonstrating the utility of the model in distinguishing different phenotypes.

Discussion
We propose an extended Gaussian mixture model for the distribution of underlying SNP-level causal genetic effects in human complex phenotypes, allowing for the phenotypespecific distribution to be modulated by heterozygosity, H, and total LD, L, of the causal SNPs, and also allowing for independent distributions for large and small effects. The GWAS z-scores for the typed or imputed SNPs, in addition to having a random environmental and error contribution, arise through LD with the causal SNPs. Thus, taking the detailed LD and heterozygosity structure of the population into account by using a reference panel, we are able to model the distribution of z-scores and test the applicability of our model to human complex phenotypes.
Complex phenotypes are emergent phenomena arising from random mutations and selection pressure. Underlying causal variants come from multiple functional categories (Schork et al., 2013), and heritability is known to be enriched for some functional categories Gazal et al., 2017;. Thus, it is likely that different variants will experience different evolutionary pressure either due to fitness directly or to pleiotropy with fitness related traits.
Here, we find evidence for markedly different genetic architectures across diverse complex phenotypes, where the effective polygenicity (or, equivalently, the prior probability that a SNP is causal with large effect) is a function of SNP total LD (L), and discoverability is multi-component and MAF dependent.
In contrast to previous work modeling the distribution of causal effects that took total LD and multiple functional annotation categories into account while implicitly assuming a polygenicity of 1 (Gazal et al., 2017), or took MAF into account while ignoring total LD dependence and different distributions for large and small effects (Zeng et al., 2018), or took independent distributions for large and small effects into account (which is related to incorporating multiple functional annotation categories) while ignoring total LD and MAF dependence, here we combine all these issues in a unified way, using an extensive underlying reference panel of ∼11 million SNPs and an exact methodology using Fourier transforms to relate summary GWAS statistics to the posited underlying distribution of causal effects. We show that the distributions of all sets of phenotypic z-scores, including extreme values that are well within genome-wide significance, are accurately reproduced by the model, both at overall summary level and when broken down with respect to a 10×10 H×TLD grid -even though the various phenotypic polygenicities and per causal-SNP heritabilities range over orders of magnitude. Improvement with respect to the basic model can be seen in the summary QQ plots Figs. 2 and 3, and also the H×TLD grids of subplots, such as Figure 4 for HDL (see also Supplemental Material, Figs. S15-S29). Compared with the basic mixture model, the extended model primarily improves the fits to data for GWAS SNPs with low total LD, in several cases for low heterozygosity as Empirical -log 1 0 (q)ˆh 2 =0.05 Model Null Empirical -log 1 0 (q)ˆh 2 =0.7 2 Empirical -log 1 0 (q)ˆh Empirical -log 1 0 (q)ˆh Simulation Coronary Artery Disease Model Null 0 2 4 6 A B C D E F G H Figure 6: QQ plots of (pruned) z-scores for simulated qualitative phenotypes (dark blue, 95% confidence interval in light blue) with model prediction (yellow). See Figure 2. The value given for p c1 is the amplitude of the full pc(L) function, which occurs at L = 1; the values (mc, wc) in parentheses following it are the total LD (mc) where the function falls to half its amplitude (the middle gray dashed lines in Figure  1 are examples), and the total LD width (wc) of the transition region (distance between flanking red dashed lines in Figure 1). Similarly for p d1 (m d , w d ), where given. h 2 b , h 2 c , and h 2 d are the heritabilities associated with the "b", "c", and "d" Gaussians, respectively. h 2 is the total SNP heritability, reexpressed as h 2 l on the liability scale for binary phenotypes. Reading the plots: on the vertical axis, choose a p-value threshold for typed SNPs (SNPs with z-scores; more extreme values are further from the origin), then the horizontal axis gives the proportion, q, of typed SNPs exceeding that threshold (higher proportions are closer to the origin). See 3 Tables C1 and C2 for a comparison of numerical values between model estimates for real phenotypes and Hapgen-based simulations where the underlying distributions of simulation causal effects were given based on the real phenotype model parameters (with σ 0 = 1). Empirical -log 1 0 (q)ˆh Empirical -log 1 0 (q)ˆh well, where very strong signals are evident. This improved model fit can be traced to the structure of the prior probability p c (L) and the variance for effect sizes, σ 2 c H S , features absent in the basic model. But even when model fits are similar, the extended model fit taking selection pressure into account results in higher likelihood and lower Bayesian information criterion.
Negative selection can be expected to result in an increasing number of effects with increasing effect size at lower heterozygosity. It was found in (Gazal et al., 2017) -which, in addition to analyzing total LD, modeled allele age and recombination rates -that common variants associated with complex traits are weakly deleterious to fitness, in line with an earlier model result that most of the variance in fitness comes from rare variants that have a large effect on the trait in question (Eyre-Walker, 2010). Thus, larger per-allele effect sizes for less common variants is consistent with the action of negative selection. Further, based on a model equivalent to Eq. 2 with p c ≡ 1, it was argued in (Zeng et al., 2018), using forward simulations and a commonly used demographic model (Gravel et al., 2011), that negative values for the selection parameter, S, which leads to larger effects for rarer variants, is a signature of negative selection. Similar results were found for the related infinitesimal model (p c ≡ 1 and π 1 ≡ 1) in .
We find negative selection parameter values for most traits, which is broadly in agreement with (Zeng et al., 2018) and  with the exception of BMI, which we find can be modeled with two Gaussians with no or weakly positive (S 0) heterozygosity dependence, though it should be noted that the polygenicity for the larger-effects Gaussian (the "c" Gaussian with the S parameter) is very low, amounting to an estimate of < 100 common causal SNPs of large effect.
A similar situation (S ≃ 0) obtains with ALS restricted to chromosome 9. Here, the sample size is relatively low (12,577 cases), which contributes to the signal being weak, and we estimate only 7 common causal SNPs associated with the "c" Gaussian.
From twin studies, the heritability of sporadic ALS has been estimated as 0.61 (0.38 to 0.78) (Al-Chalabi et al., 2010); a clinically ascertained case series estimated the heritability to be between 0.40 and 0.45 (Wingo et al., 2011). Mutations of the chromosome 9 open reading frame 72 gene C9orf72 are implicated in both familial and sporadic ALS (Van Blitterswijk et al., 2012), but other chromosomes are also known to be involved . However, Supplemental Material, Figure S8, showing QQ plots accounting for all GWAS SNPs, implies that all the GWAS signal comes from chromosome 9. It should be noted that the mixture model, which is primarily focused on characterizing a distribution of undiscovered effect sizes, is designed to capture a broad-based polygenic contribution, i.e., not dominated by singular variants or effects from very large LD blocks. In implementing the model, we limit to SNPs with total LD<600 and randomly prune GWAS SNPs in LD blocks. From Figure 2 B and Figure  S8, however, the few excluded GWAS SNPs clearly can not account for missing the preponderance of heritability. It is likely that the GWAS data lacks rare variants pertinent to ALS, and in addition is underpowered.
Our BMI results are not in agreement with the results of others. For the UKB (2015) (Locke et al., 2015), Zhang et al Zhang et al. (2018) report h 2 = 0.20 and n causal = 18k. It is not clear why our BMI results are in such disagreement with these, except to note that we have a two-component (versus single component) causal Gaussian and the majority of the heritability comes from the small-effects "b" Gaussian. For height, our selection parameter S = −0.46 (2014 GWAS) is in reasonable agreement with S = −0.422 reported in (Zeng et al., 2018), and with α = −0.45 reported in .
Compared with our basic model results (Holland et al., 2020), we generally have found that using the extended model presented here, heritabilities show marked increases, with large contributions form the selection ("c") Gaussian -see 3, Table B1. For example, our basic-model heritability estimate for HDL was 7%, in close agreement with an earlier LDSC estimate (Speed and Balding, 2019) and the M2 model of Zhang et al. (2018) (9%), though the heritability from discovered loci was known to be 12% (Willer et al., 2013); our extended-model estimate is 19%, and the fit to the data is a dramatic improvement over the basic model, as can be seen in Figure 4 (it should be borne in mind that mixture models are designed to capture broad-based polygenic contributions to heritability, not outlier effects). However, due to the small effect sizes associated with the "b" Gaussian as a component of the extended model, for some phenotypes the overall number of causal SNPs can be considerably larger than previously estimated. The concept of effective number of causal SNPs, M e , in O'Connor et al. (2019), roughly corresponds to n c + n d here, though the relationship is not precise (see also 3); e.g., for schizophrenia, M e = 15k versus n c + n d = 26k. Supplemental Material, Figures S11-S14 capture the breakdown in SNP contributions to heritability as a function of heterozygosity and total linkage disequilibrium.
From Eq. 4, β is a weighted sum of contributions from Gaussians of different variance. Since we find σ b < σ c and S < 0 -which increasingly magnifies the difference in variance of the two Gaussians as H gets smaller -we find larger E(β 2 ) for rarer variants. Also, as p c (L) increases (from 0) with decreasing L, for a given H we also find that as L decreases, per causal-SNP heritability and E(β 2 ) increase, consistent with Gazal et al. (2017). These patterns can be seen in Supplemental Material, Figures S11-S14 (second and fourth columns). The per causal-SNP contribution to heritability (second columns) is found to be more smoothly varying across common and low-frequency variants than E(β 2 ), in broad agreement with O'Connor et al. (2019). It was also found in Gazal et al. (2017) that more recent common alleles have both lower LLD and larger per-SNP heritability (all SNPs causal); since selection has had less time to remove recent deleterious alleles, larger per-SNP heritability from SNPs with lower LLD was indicative of negative selection.
Generally, we find evidence for the existence of genetic architectures where the per causal-SNP heritability is larger for more common SNPs with lower total LD. But the trend is not uniform across phenotypes -see Supplemental Material, Figures S11-S14, second columns. The observed-scale heritability estimates given by h 2 b correspond to effects not experiencing much selection pressure. The new final values of h 2 presented here result from a model that, compared with the basic Gaussian mixture model it is an extension of, gives better fits between data and model prediction of the summary and detailed QQ plots, and thus constitute more accurate estimates of SNP heritability. For 2010 and 2014 height GWASs, we obtain very good consistency for the model parameters and therefore heritability, despite considerable difference in apparent inflation. The 2018 height GWAS  has a much larger sample size (almost three quarters of a million); the slightly different parameter estimates for it might arise due to population structure not fully captured by our model (Sohail et al., 2019;Berg et al., 2019). Our h 2 estimates for height, however, remain consistently lower than other reported results (e.g., h 2 = 0.33 in (Zhang et al., 2018), h 2 = 0.527 in (Zeng et al., 2018), h 2 = 0.61 in ). For educational attainment, our heritability estimate agrees with (Zeng et al., 2018) (h 2 = 0.182), despite our using a sample more than twice as large. The difference in selection parameter value, S = −0.335 in Zeng et al. (2018) versus S = −0.44 here, might partially be explained by the model differences (Zeng et al use one causal Gaussian with MAF-dependence but no LD dependence). A comparison of the total heritabilities reported here with estimates from the work of others (Zhang et al., 2018;Zeng et al., 2018;Schoech et al., 2019), in addition to our basic model alone (Holland et al., 2020), are given in 3, Table B1.
For most traits, we find strong evidence that causal SNPs with low heterozygosity have larger effect sizes (S < 0 in Table 1; the effect of this as an amplifier of σ 2 c in Eq. 5 is illustrated in Supplemental Material, Figure S9) -see also Supplemental Material, Figures S11-S14, fourth columns. Thus, negative selection seems to play an important role in most phenotypes-genotypes. This is also indicated by the extent of the region of finite probability for variants of large effect sizes (which are enhanced by having S −0.4) being relatively rare (low H), which will be greater for larger p c1 , the amplitude of the prior probability for the "c" Gaussian (see Supplemental Material, Figures S3 and S4).
The "b" Gaussian in Eq. 4 or 5 does not involve a selection parameter: effect size variance is independent of MAF. Thus, causal SNPs associated with this Gaussian are likely undergoing neutral (or very weakly negative) selection. It should be noted that in all traits examined here, whether or not there is evidence of negative selection (S < 0), the effect size variance of the "b" Gaussian is many times smaller -sometimes by more than an order of magnitude -than that for the "c" Gaussian. Thus, it appears there are many causal variants of weak effect undergoing neutral (or very weakly negative) selection. For the nine phenotypes where the "d" Gaussian could be im-plemented, its variance parameter was several times larger than that of the "c" Gaussian. However, the amplitude of the prior probability for the "d" Gaussian, p d1 , was generally much smaller than the amplitudes of the prior probabilities for the "b" or "c" Gaussians, which translated into a relatively small number of causal variants with very large effect associated with this Gaussian. (Due to lack of power, in three instances -CD, UC, and TC -p d (L) was treated as a constant, i.e., independent of L.) Interestingly, intelligence had the highest number of causal SNPs associated with this Gaussian, while the extent of total LD for associated SNPs was also liberal (m d = 561; see also Supplemental Material, Figure S5). It is possible that some of these SNPs are undergoing positive selection, but we did not find direct evidence of that.
A limitation of this work is that we do not take SNP functional categories into account. Inclusion of SNP annotation -e.g., by having a set of model parameters for each of many functional categories, and subdividing each SNP's total LD into contributions from these categories -is important for deriving more biologically-informed interpretations of genetic effects. However, for the summary quantities estimated here, annotation is not expected to have a large impact; indeed, O'Connor et al. (2019) conclude that the accuracy of S-LD4M is not contingent on modeling annotation. Another limitation, a feature of many large GWAS, is that rare and disease-specific SNPs are not included. Thus, our analysis strictly applies only to the spectrum of relatively common SNPs. Indeed, our results point to the importance of rare variants in order to more comprehensively study the evolutionary architecture of complex phenotypes.
We find a diversity of genetic architectures across multiple human complex phenotypes where SNP total LD plays an important role in effect size distribution. In general, lower total LD SNPs are more likely to be causal with larger effects. Furthermore, for most phenotypes, while taking total LD into account, causal SNPs with lower MAF have larger effect sizes. These phenomena are consistent with models of the action of negative selection. Additionally, for all phenotypes, we find evidence of neutral selection operating on SNPs with relatively weak effect. We did not find direct evidence of positive selection. Compared with the basic Gaussian mixture model, which did not take heterozygosity or total LD into account in the distribution of effect sizes, the extended model consistently provided a much better fit to the distribution of GWAS summary statistics, thus providing more accurate estimates of genetic quantities of interest. Future work will explore SNP functional annotation categories and their differential roles in human complex phenotypes.

Funding
Research Council of Norway (262656,248984,248778,223273) and KG Jebsen Stiftelsen; NIH Grant Number U24DA041123.  Table B1: Comparison of our current extended model heritability estimates (h 2 (l) , see Table 3) with estimates from: our earlier basic model (dressed with a tilde,h 2 (l) ) (Holland et al., 2020); the two-and three-component Gaussian models in (Zhang et al., 2018), denoted with subscript M 2 and M 3, respectively; the single Gaussian model with selection parameter in (Zeng et al., 2018), subscripted M S; and the single Gaussian infinitesimal model (polygenicity=1) with selection parameter in , subscripted α. M 2 and M 3 use the HapMap3 reference panel (Consortium et al., 2010) with 1.07 million common SNPs (MAF 0.05); M S uses an Affymetrix panel with 483,634 SNPs (MAF> 0.01) on UK Biobank data; our results are based on a 1000 Genomes Phase3 reference panel with 11 million SNPs (MAF 0.002). SE denotes standard error. h 2 (l) is our heritability estimate (see Table 3); those obtained from M 2, M 3, and M S are labeled h 2 M 2 , h 2 M 3 , and h 2 M S , respectively. For quantitative phenotypes, values are on the observed scale; for binary phenotypes, values are on the liability scale, using the same population prevalence as used for Table 3. It is important to note that our heritability estimates are corrected for inflation, by dividing by the inflation parameter σ 2 0 ; this is not done for M 2, M 3 or M S. * Disease prevalences are given as a percentage in parentheses. For binary phenotypes, let h 2 obs denote the heritability on the observed 0-1 scale (this is h 2 in Figure 2). Let P denote the proportion of cases in the study: P = Ncases/(Ncases + N controls ). Then the heritability on the log-odds-ratio scale reported in (Zhang et al., 2018) is h 2 log = h 2 obs /(P (1 − P )). The transformation between the observed and liability scale is given by Eq. 39 in (Holland et al., 2020). † N is the total sample size for quantitative traits; for qualitative traits, N ef f = 4/(1/Ncases + 1/N controls ) -see main text.  Table B2: Comparison of the basic model estimates (Holland et al., 2020) of the number of causal SNPs (here denoted n β ) and heritability (here denoted h 2 β ) with the corresponding net contributions nc + n d and h 2 c + h 2 d from the large effects components ("c" and "d") of the extended model. (All disease heritabilities here are on the observed scale; see Table 3.) For BMI , it is the "b" component of the extended model that dominates, and a comparison with the basic model gives: n β = 7.5E3, n b = 1.7E4; h 2 β = 0.07, h 2 b = 0.08. * Model C for Height 2018.