Interpreting SNP heritability in admixed populations

SNP heritability hsnp2 is defined as the proportion of phenotypic variance explained by genotyped SNPs and is believed to be a lower bound of heritability (h2), being equal to it if all causal variants are known. Despite the simple intuition behind hsnp2, its interpretation and equivalence to h2 is unclear, particularly in the presence of population structure and assortative mating. It is well known that population structure can lead to inflation in hˆsnp2 estimates because of confounding due to linkage disequilibrium (LD) or shared environment. Here we use analytical theory and simulations to demonstrate that hsnp2 estimates can be biased in admixed populations, even in the absence of confounding and even if all causal variants are known. This is because admixture generates LD, which contributes to the genetic variance, and therefore to heritability. Genome-wide restricted maximum likelihood (GREML) does not capture this contribution leading to under- or over-estimates of hsnp2 relative to h2, depending on the genetic architecture. In contrast, Haseman-Elston (HE) regression exaggerates the LD contribution leading to biases in the opposite direction. For the same reason, GREML and HE estimates of local ancestry heritability hγ2 are also biased. We describe this bias in hˆsnp2 and hˆγ2 as a function of admixture history and the genetic architecture of the trait and show that it can be recovered under some conditions. We clarify the interpretation of hˆsnp2 in admixed populations and discuss its implication for genome-wide association studies and polygenic prediction.


Introduction
The ability to estimate (narrow-sense) heritability (h 2 ) from unrelated individuals was a major advance in genetics.Traditionally, h 2 was estimated from family-based studies in which the phenotypic resemblance between relatives could be modeled as a function of their expected genetic relatedness [1].However, this approach was limited to analysis of closely related individuals where pedigree information is available and the realized genetic relatedness is not too different from expectation [2].With the advent of genome-wide association studies (GWAS), we hoped that many of the variants underlying this heritability would be uncovered.However, when genome-wide significant SNPs explained a much smaller fraction of the phenotypic variance, it became important to explain the missing heritability -were family-based estimates inflated or were GWAS just underpowered, limited by variant discovery?Yang et al. (2010) [3] made the key insight that one could estimate the portion of h 2 tagged by genotyped SNPs, regardless of whether or not they were genome-wide significant, by exploiting the subtle variation in the realized genetic relatedness among apparently unrelated individuals [3][4][5].This quantity came to be known colloquially as 'SNP heritability' (h 2 snp ) and it is believed to be equal to h 2 if all causal variants are included among genotyped SNPs [3].Indeed, estimates of h 2 snp explain a much larger fraction of trait heritability than GWAS SNPs [3], approaching family-based estimates of h 2 when whole genome sequence data, which captures rare variants, are used [6].This has made it clear that GWAS have yet to uncover more variants with increasing sample size.Now, h 2 snp has become an important aspect of the design of genetic studies and is often used to define the power of variant discovery in GWAS and the upper limit of polygenic prediction accuracy.
Despite the utility and simple intuition of h 2 snp , there is much confusion about its interpretation and equivalence to h 2 , particularly in the presence of population structure and assortative mating [7][8][9][10][11][12].
But much of the discussion of heritability in structured populations has focused on biases in ĥ snp -the estimator -due to confounding effects of shared environment and linkage disequilibrium (LD) with other variants [7,[9][10][11]13].There is comparatively little discussion, at least in human genetics, on the fact that LD due to population structure also contributes to genetic variance, and therefore, is a component of heritability [1] (but see [14][15][16] for a rigorous discussion).We think this is at least partly due to the fact that most studies are carried out in cohorts with primarily European ancestry, where the degree of population structure is minimal and large effects of LD can be ignored.However, that is not the case for diverse, multi-ethnic cohorts, which have historically been underrepresented in genetic studies, but thanks to a concerted effort in the field, are now becoming increasingly common [17][18][19][20][21][22][23].The complex structure in these cohorts also brings unique methodological challenges and it is imperative that we understand whether existing methods, which have largely been evaluated in more homogeneous groups, generalize to more diverse cohorts.
Our goal in this paper is to study the behavior of ĥ snp in admixed populations.What is its interpretation in the ideal situation where causal variants are known?Is it an unbiased estimate of h 2 ?To answer these questions, we derived a general expression for the genetic variance in admixed populations, decomposing it in terms of the contribution of population structure, which influences both the genotypic variance at individual loci and the LD across loci.We used theory and simulations to show that ĥ snp of copies of the trait-increasing allele.Importantly, the effect sizes are fixed quantities and differences in genetic values among individuals are due to random variation in genotypes.Note, that this is different from the model assumed by GREML where genotypes are fixed and effect sizes are random [14].
We denote the mean, variance, and covariance with E(.), V(.), and C(.), respectively, where the expectation is measured over random draws from the population rather than random realizations of the evolutionary process.We can express the additive genetic variance of a quantitative trait as follows: Here the first term represents the contribution of individual loci (genic variance) and the second term is the contribution of linkage disequilibrium (LD contribution).We make the assumption that loci are unlinked and therefore, the LD contribution is entirely due to population structure.We describe the behavior of V g in a population that is a mixture of two previously isolated populations A and B that diverged from a common ancestor.To do this, we denote θ as the fraction of the genome of an individual with ancestry from population A. Thus, θ = 1 if the individual is from population A, 0 if they are from population B, and θ ∈ (0, 1) if they are admixed.Then, V g can be expressed in terms of ancestry as (Appendix): where f A i and f B i are the allele frequencies in populations A and B, and E(θ) and V(θ) are the mean and variance of individual ancestry.The sum of the first three terms represents the genic variance and the last term represents the LD contribution.

Demographic history
From Eq. 1, it is clear that, conditional on the genetic architecture in the source populations (β, f A , f B ), V g is a function of the mean, E(θ), and variance, V(θ), of individual ancestry in the admixed population.
We consider two demographic models that affect E(θ) and V(θ) in qualitatively different ways.In the first model, the source populations meet once t generations ago (we refer to this as t = 0) in proportions p and 1 − p, after which there is no subsequent admixture (Fig. 1A).In the second model, there is continued gene flow in every generation from one of the source populations such that the mean overall amount of ancestry from population A is the same as in the first model (Fig. 1A).For brevity, we refer to these as the hybrid-isolation (HI) and continuous gene flow (CGF) models, respectively, following Pfaff et al. (2001) [25].V(θ) is also affected by ancestry-based assortative mating, where individuals are more likely to partner with others of similar ancestry.We refer to this simply as assortative mating for brevity and model this following Zaitlen et al. (2017) using a parameter P ∈ (0, 1), which represents the correlation of the ancestry of individuals across mating pairs in the population [26].
Under these conditions, the behavior of E(θ) and V(θ) has been described previously [26,27] (Fig. 1B and C).Briefly, in the HI model, E(θ) remains constant at p in the generations after admixture as there is no subsequent gene flow.V(θ) is at its maximum at t = 0 when each individuals carries chromosomes either from population A or B, but not both.This genome-wide correlation in ancestry breaks down in subsequent generations as a function of mating, independent assortment, and recombination, leading to a decay in V(θ), the rate depending on the strength of assortative mating (Fig. 1C).In the CGF model, both E(θ) and V(θ) increase with time as new chromosomes are introduced from the source populations.But while E(θ) continues to increase monotonically, V(θ) will plateau and decrease due to the countervailing effects of independent assortment and recombination which redistribute ancestry in the population, reaching equilibrium at zero if there is no more gene flow and the population is mating randomly.V(θ) provides an intuitive and quantitative measure of the degree of population structure (along the axis of ancestry) in admixed populations.

Genetic variance in admixed populations
To understand the expectation of genetic variance in admixed populations, it is first worth discussing its behavior in the source populations.In Eq. 1, the first term represents the within-population component (V gw ) and the last three terms altogether represent the component of genetic variance between populations A and B (V gb ).Note that V gb = ( ḡA − ḡB ) 2 2 is positive only if there is a difference in the mean genotypic values (Fig. 2).This variance increases with genetic divergence since the expected values increase monotonically with increasing divergence, ) is expected to be zero under neutrality because the direction of frequency change will be uncorrelated across loci.In this case, the LD contribution, i.e., (1.4), is expected to be zero and V gb = (1.1)+ (1.2) + (1.3).However, this is true only in expectation over the evolutionary process and the realized LD contribution may be non-zero even for neutral traits.| ḡA − ḡB | is expected to be greater than that expected under genetic drift.For traits under stabilizing selection, | ḡA − ḡB | will be less than that expected under genetic drift, and zero in the extreme case.
For traits under selection, the LD contribution is expected to be greater or less than zero, depending on the type of selection.Under divergent selection, trait-increasing alleles will be systematically more frequent in one population over the other, inducing positive LD across loci [28,29], increasing the LD contribution, i.e., term (1.4).Stabilizing selection, on the other hand, induces negative LD [30,31].In the extreme case, the mean genetic values of the two populations are exactly equal and V gb = (1.2) + (1.3) + (1.4) = 0.For this to be true, (1.4) has to be negative and equal to (1.2) + (1.3), which are both positive, and the total genetic variance is reduced to the within-population variance, i.e., term (1.1) (Fig. 2).This is relevant because, as we show in the following sections, the behavior of the genetic variance in admixed populations depends on the magnitude of V gb between the source populations.
We illustrate this by tracking the genetic variance in admixed populations for two traits, both with the same mean F ST at causal loci but with different LD contributions (term 1.4): one where the LD contribution is positive (Trait 1) and the other where it is negative (Trait 2).Thus, traits 1 and 2 can be thought of as examples of phenotypes under divergent and stabilizing selection, respectively, and we refer to them as such from hereon.To simulate the genetic variance of such traits, we drew the allele frequencies (f A and f B ) in populations A and B for 1,000 causal loci with F ST ∼ 0.2 using the Balding-Nichols model [32].We drew their effects (β) from N (0, where f is the mean allele frequency between the two populations, m is the number of loci.To simulate positive and negative LD, we permuted the effect signs across variants 100 times and selected the combinations that gave the most positive and negative LD contribution to represent the genetic architecture of traits that might be under directional (Trait 1) and stabilizing (Trait 2) selection, respectively (Methods).We simulated the genotypes of 10,000 individuals under the HI and CGF models for t ∈ {10, 20, 50, 100} generations post-admixture and calculated genetic values for both traits using g = m i=1 β i x i , where m = 1, 000 (Method).The observed genetic variance at any time can then be calculated simply as the variance in genetic values, i.e.V g = V(g).
In the HI model, E(θ) does not change (Fig. 1B) so terms (1.1) and (1.2) are constant through time.
Terms (1.3) and (1.4) decay towards zero as the variance in ancestry goes to zero and V g ultimately converges to (1.1) + (1.2) (Fig. 3).This equilibrium value is equal to the E(V g |θ) (Appendix) and the rate of convergence depends on the strength of assortative mating, which slows the rate at which V(θ) decays.V g approaches equilibrium from a higher value for traits under divergent selection and lower value for traits under stabilizing selection because of positive and negative LD contributions, respectively, at t = 0 (Fig. 3).In the CGF model, V g increases initially for both traits with increasing gene flow (Fig. 3).This might seem counter-intuitive at first because gene flow increases admixture LD, which leads to more negative values of the LD contribution for traits under stabilizing selection (Fig. S1).But this is outweighed by positive contributions from the genic variance -terms (1.1) + (1.2) + (1.3) -all of which initially increase with gene flow (Fig. S1).After a certain point, the increase in V g slows down as any increase in V(θ) due to gene flow is counterbalanced by recombination and independent assortment.
Ultimately, V g will decrease if there is no more gene flow, reaching the same equilibrium value as in the HI model, i.e., E(V g |θ) = (1.1)+ (1.2).Because the loci are unlinked, we refer to the sum (1.3) + (1.4) as the contribution of population structure.

GREML estimation
In their original paper, Yang et al. (2010) defined h 2 snp as the variance explained by genotyped SNPs and not as heritability [3].This is because h 2 is the genetic variance explained by causal variants, which are unknown.Genotyped SNPs may not overlap with or tag all causal variants and thus, h 2 snp is understood to be a lower bound of h 2 , both being equal if causal variants are known [3].Our goal is to demonstrate that this may not be true in structured populations and quantify the bias in ĥ snp , even in the ideal situation when causal variants are known.
We used GREML, implemented in GCTA [3,5], to estimate the genetic variance for our simulated traits.GCTA assumes the following model: y = Zu + where Z is an n × m standardized genotype matrix such that the genotype of the k th individual at the i th locus is , f i being the ) and ∼ N (0, Iσ 2 ) is random environmental error.Then, the phenotypic variance can be decomposed as: where ZZ m is the genetic relationship matrix (GRM), the variance components σ 2 u and σ 2 are estimated using restricted maximum likelihood, and ĥ snp is calculated as u +σ 2 .We are interested in asking whether σ2 u is an unbiased estimate of V g .To answer this, we constructed the GRM with causal variants and estimated σ2 u using GCTA [3,4].
GCTA under-and over-estimates the genetic variance in admixed populations for traits under divergent (Trait 1) and stabilizing selection (Trait 2), respectively, when there is population structure, i.e., when V(θ) > 0 (Fig. 4A).One reason for this bias is that the GREML model assumes that the effects are independent, and therefore the LD contribution is zero.This, as discussed in the previous section, is not true for traits under divergent or stabilizing selection between the source populations, and only true for neutral traits in expectation.Because of this, σ2 u does not capture the LD contribution, i.e. term (1.4) (Fig. 4A).But σ2 u can be biased even if the LD contribution is zero if the genotypes are scaled with 2f i (1 − f i ) -the standard practice -where f i is the frequency of the allele in the population.This scaling assumes that V(x i ) = 2f i (1 − f i ), which is true only if the population were mating randomly.
In an admixed population V( , where f i , f A i , and f B i correspond to frequency in the admixed population, and source populations, A and B, respectively (Appendix).
Alternatively, if the genotypes are scaled, V( st where st is the F st at the i th locus.
We show that this assumption biases σ2 u downwards by a factor of 2 genotypes are scaled) -term (1.3) (Fig. 4B, Appendix).Thus, with the standard scaling, σ2 u gives a biased estimate in the presence of population structure, even of the genic variance.
The overall bias in σ2 u is determined by the relative magnitude and direction of terms (1.3) and (1.4), both of which are functions of V(θ), and therefore, of the degree of structure in the population.
The contribution of term (1.3) will be modest, even in highly structured populations (Fig. S1) and therefore, the overall bias is largely driven by the LD contribution.If there is no more gene flow, V(θ) will ultimately go to zero and V g will converge towards σ2 u .Thus, σ2 u is more accurately interpreted as the genetic variance expected if the LD contribution were zero and if the population were mating randomly.In other words, E(σ 2 u ) = (1.1)+ (1.2) = V g (Fig. 4B).
In principle, we can recover the missing components of V g by scaling the genotypes appropriately.
For example, we can recover term (1.3) by scaling the genotype at each variant i by its sample variance, (Fig. 4C) (Appendix).We can also recover term (1.4) by scaling the genotypes with the covariance between SNPs, i.e., the LD matrix, as previously proposed [33,34] (Methods).In matrix form, the 'LD-scaled' genotypes can be written as Z = (X − 2P )U −1 where P is an n × m matrix such that all elements of the i th column contain the frequency of the i th SNP and U is the (upper triangular) square root matrix of the LD matrix, i.e., Σ = U U [33].GREML recovers the LD contribution under

F-H). (A-B & E-F) show the behavior of σ2
u for the default scaling, (C, G) shows σ2 u when the genotype at a locus is scaled by its sample variance (V(x) scaled), and (D, H) when it is scaled by the sample covariance (LD scaled).
this scaling, resulting in unbiased estimates of V g for both traits (Fig. 4D, Appendix).
In practice, however, the LD contribution may not be fully recoverable for two reasons.One, the LD-scaled GRM requires computing the inverse of Σ or U which may not exist, especially if the number of markers is greater than the sample size -the case for most human genetic studies.Second, it is common to include individual ancestry or principal components of the GRM as fixed effects in the model to account for inflation in heritability estimates due to shared environment.This should also have the effect of removing the components of genetic variance along the ancestry axes, the residual variance being equal to E{V(g|θ)} = (1.1)+ (1.2) − (1.3) (Appendix).Indeed, this is what we observe in Fig. 4H.Thus, if ancestry is included as a fixed effect, we expect V g to be underestimated in the presence of population structure, regardless of genetic architecture.

HE estimation
Haseman-Elston (HE) regression also assumes a random-effects model but uses a method-of-moments approach, as opposed to GREML, which maximizes the likelihood to estimate V g .Previous work has shown that as long as all causal variants are included in the GRM calculation, the HE estimator will not be biased, even if they are in LD with each other [35].We show that in the presence of positive and negative LD between causal loci, as exemplified by traits under divergent and stabilizing selection, respectively, the HE estimates of V g are biased upwards and downwards, respectively (Fig. 5A-B).To understand this discrepancy and the source of bias in our simulations, recall that HE estimates V g from the regression of the (pairwise) phenotypic covariance between individuals on their genotypic covariance [24].More specifically, if we denote Y kl = y k y l as the product of the (centered) phenotypes of k th and l th individuals, and ψ kl as the k th and l th entry of the GRM, then the HE estimator can be written as: Where the first and second terms represent the genic and LD components, respectively, of the estimate.Population structure induces correlations between the alleles at a given locus as well as across loci (i.e., LD).But the LD may not be directional, i.e., trait-increasing alleles may be as likely to be co-inherited with each other as they are to trait-decreasing alleles, and vice versa -implicit under the standard random-effects model.Thus, in the absence of directional LD, the second term is zero and the first term is unaffected as long as all causal variants are included in the GRM, because the increase in the numerator due to population structure is proportional to the denominator [35].Directional LD does not affect the first term but exaggerates the contribution from the second term, i.e., the LD component (see Appendix section A3.2).Consequently, HE regression over-and under-estimates V g for traits with positive and negative LD, respectively.Note that this bias is in the opposite direction of the bias observed with GREML, which fails to capture the LD contribution.Scaling the genotype at a locus by its LD with other loci, as discussed in the previous section, corrects for the bias in HE regression regardless of genetic architecture, yielding estimates consistent with GREML (Fig. 5C).Thus, GREML and HE regression are guaranteed to yield the same estimates only if the underlying model specifying the distribution of effects is consistent with the true architecture of the trait.
The practice of including individual ancestry as a covariate in HE regression to account for shared environment [11] reduces the bias from exaggerated LD contributions (Fig. 5D-F).But, as with GREML, this also removes any genetic variance that may exist along the ancestry axis, yielding underestimates of V g , regardless of genetic architecture.

Local ancestry heritability
A related quantity of interest in admixed populations is local ancestry heritability (h 2 γ ), which is defined as the proportion of phenotypic variance that can be explained by local ancestry.Zaitlen et al. (2014) [36] showed that this quantity is related to, and can be used to estimate, h 2 in admixed populations.
The advantage of this approach is that local ancestry segments shared between individuals are identical by descent and are therefore, more likely to tag causal variants compared to array markers, allowing one to potentially capture the contributions of rare variants [36].Here, we show that in the presence of population structure, (i) the relationship between h 2 γ and h 2 is not straightforward and (ii) ĥ2 γ may be a biased estimate of local ancestry heritability under the random effects model for the same reasons that ĥ2 snp is biased.
We define local ancestry γ i ∈ {0, 1, 2} as the number of alleles at locus i that trace their ancestry to population A. Thus, ancestry at the i th locus in individual k is a binomial random variable with E(γ ik ) = 2θ k , θ k being the ancestry of the k th individual.Similar to genetic value, the 'ancestry value' of an individual can be defined as is the effect size of local ancestry (Appendix).Then, the genetic variance due to local ancestry can be expressed as (Appendix): and heritability explained by local ancestry is simply the ratio of V γ and the phenotypic variance.Note ) and therefore its behavior is similar to V g in that the terms (1.3) and (1.4) decay towards zero as V(θ) → 0, and V γ converges to (1.2) (Fig. S2).Additionally, the dependence of V γ on both E(θ) and V(θ) precludes a straightforward derivation between local ancestry heritability and GREML estimation of ĥ2 γ is similar to that of ĥ2 snp , the key difference being that the former involves constructing the GRM using local ancestry instead of genotypes [36].The following model is assumed: ancestry effects, and ξ ∼ N (0, Iσ 2 ξ ).Note that σ 2 ξ captures both environmental noise as well as any genetic variance independent of local ancestry.The phenotypic variance is decomposed as is the local ancestry GRM and σ 2 v is the parameter of interest, which is believed to be equal to V γ -the genetic variance due to local ancestry.
We show that, in the presence of population structure, i.e., when V(θ) > 0, GREML σ2 v is biased downwards relative to V γ for traits under divergent selection and upwards for traits under stabilizing selection because it does not capture the contribution of LD (Fig. 6A).But there is another source of bias in σ2 v , which tends to be inflated in the presence of population structure if individual ancestry is not included as a covariate, even with respect to the expectation of V γ under equilibrium (seen more clearly in Fig. 6B-C).We suspect this inflation is because of strong correlations between local ancestry -local ancestry disequilibrium -across loci that inflates σ2 v in a way that is not adequately corrected even when all causal variants are included in the model [4,10].Scaling local ancestry by its covariance removes this bias and recovers the contribution of LD (Fig. 6D) presumably because this accounts for the correlation in genotypes across loci.Including individual ancestry as a fixed effect also corrects for the inflation in σ2 v (Fig. 6E-H).But as with σ2 u , this practice will underestimate the genetic variance due to local ancestry in the presence of population structure because it removes the variance along the ancestry axis (Fig. 6E-H).
Based on the above, GREML ĥ2 γ and corresponding estimates of h 2 are more accurately interpreted as the heritability due to local ancestry and heritability, respectively, expected in the absence of population structure.We believe ĥ2 γ is still useful in that, because it should capture the effects of rare variants, it can be used to estimate the upper bound of ĥ snp .
In a previous paper, we suggested that local ancestry heritability could potentially be used to estimate the genetic variance between populations [37].Our results suggest this is not possible for two reasons.
First, the GREML estimator of local ancestry heritability, as we show in this section is biased and does not capture the LD contribution.But even if we were able to recover the LD component, our decomposition shows that local ancestry is equal to the genetic variance between populations (V gb ) only when E(θ) = 0.5 and V(θ) = E(θ){1−E(θ)} = 0.25, which is only possible at t = 0 in the HI model.After admixture, V(θ) decays and the equivalence between V γ and V gb is lost, making it impossible to estimate the latter from admixed populations, especially for traits under divergent or stabilizing selection, even if the environment is randomly distributed with respect to ancestry.We note that this conclusion was recently reached independently by Schraiber and Edge (2023) [38].How much does LD contribute to V g in practice?
In the previous sections, we showed theoretically that ĥ snp may be biased in admixed populations even if the causal variants are known and in the absence of confounding by shared environment.GREML fails to capture the LD contribution whereas HE regression overestimates it.The extent to which ĥ snp is biased because of this reason in practice is ultimately an empirical question, which is difficult to answer because the true genetic architecture -the LD contribution in particular -is unknown.In this section, we develop some intuition for this contribution among individuals with mixed African and European ancestry using a combination of simulations and empirical data analysis.
First, we simulated a neutral trait using genotype data from the African Americans (ASW) from the 1,000 Genomes Project (1KGP) [39].To do this, we sampled m ∈ {10, 100, 1, 000} causal loci from a set of common (MAF > 0.01), LD pruned variants and assigned them effects such that , i.e., the expected genic variance is E{ m i=1 β 2 i V ar(x i )} = 1 (Methods).We computed the genic and LD contributions and repeated this process 1,000 times where each replicate can be thought of as an independent realization of the genetic architecture of a neutrally evolving trait.We show that the LD contribution may be zero in expectation but can be substantial for a given trait (up to 50% of the genic variance, Fig. S4), even in the absence of selection.
Second, we estimated the LD contribution of genome-wide significant SNPs for 26 quantitative traits from the GWAS catalog [40].To do this, we decomposed the variance explained in ASW into the four components in Equation 1 using allele frequencies (f A and f B ) from the YRI and CEU and the mean (E(θ) ≈ 0.77) and variance (V(θ) ≈ 0.02) of individual ancestry from ASW (Methods).We show that for skin pigmentation -a trait under strong divergent selection -the LD contribution, i.e. term (1.4), is positive and accounts for ≈ 40 -50% of the total variance explained.This is because of large allele frequency differences between Africans and Europeans that are correlated across skin pigmentation loci, consistent with strong polygenic selection favoring alleles for darker pigmentation in regions with high UV exposure and vice versa [37,[41][42][43][44].But for most other traits, LD contributes relatively little, explaining a modest, but non-negligible proportion of the genetic variance in height, LDL and HDL cholestrol, mean corpuscular hemoglobin (MCH), neutrophil count (NEU), and white blood cell count (WBC) (Fig. 7).
Because we selected independent associations for this exercise (Methods), the LD contribution is driven entirely due to population structure in ASW.The contribution of population structure to the genic variance, i.e., term (1.3) is also small even for traits like skin pigmentation and neutrophil count with large effect alleles that are highly diverged in frequency between Africans and Europeans [42,43,[45][46][47].Overall, this suggests that population structure contributes relatively little, as least to the variance explained by GWAS SNPs.

Discussion
Despite the growing size of GWAS and discovery of thousands of variants for hundreds of traits [40], the heritability explained by GWAS SNPs remains a fraction of twin-based heritability estimates.causal variants but assumes that they are numerous and are more or less uniformly distributed across the genome (the infinitesimal model), their contributions to the genetic variance 'tagged' by genotyped SNPs [3].h 2 snp is now routinely estimated in most genomic studies and at least for some traits (e.g.height and BMI), these estimates now approach twin-based heritability [6].But despite the widespread use of ĥ snp , its interpretation remains unclear, particularly in the presence of admixture and population structure.
It is generally accepted that ĥ snp can be biased in structured populations because of confounding effects of unobserved environmental factors and LD between causal variants [4,7,[9][10][11]48].But ĥ snp may be biased even in the absence of confounding because of misspecification of the underlying random-effects model, i.e., if the model does not represent the genetic architecture from which the trait is sampled [14-16, 49, 50].
Under the standard GREML model, SNP effects are assumed to be uncorrelated and the total genetic variance can be represented as the sum of the variance explained by individual loci, i.e. the genic variance [14][15][16].In admixed populations, there is substantial LD, which can contribute to the genetic variance, and can persist for a number of generations, despite recombination, due to continued gene flow and/or ancestry-based assortative mating.GREML does not capture this LD contribution [12,15], and therefore, may lead to biased estimates of h 2 snp .The LD contribution can be negative for traits under stabilizing selection, and positive for traits under divergent selection between the source populations, leading to over-or under-estimates, respectively.Thus, GREML estimates of h 2 snp , assuming genotypes are scaled properly (see below), is better interpreted as the proportion of phenotypic variance explained by the genic variance.Estimates of local ancestry heritability ( ĥ2 γ ) [36,51] should be interpreted similarly.
We show that with GREML, ĥ snp can be biased even when the LD contribution is zero if the genotypes are scaled by 2f (1 − f ) -the standard approach, which implicitly assumes a randomly mating population.In the presence of population structure, the variance in genotypes can be higher and ĥ snp does not capture this additional variance, which we show can be recovered by scaling genotypes by the SNP variance ( V ar(x)).In principle, the LD contribution can also be recovered by scaling genotypes by the SNP covariance, i.e., the LD matrix, as previously suggested [33,34].But this approach is limited to situations where the sample size is much larger than the number of markers.
We also investigated the behavior of another widely used approach to estimate h 2 snp -Haseman-Elston regression.We show that ĥ snp estimated with HE regression is also biased, but for different reasons and in the opposite direction of the bias observed with GREML.HE regression exaggerates the LD contribution, leading to over-and under-estimates of h 2 snp for traits where the causal loci are in positive and negative LD, respectively.Approaches that correct for population structure [35] should remove this source of bias but would also remove any genetic variance in the trait along the ancestry axis, including the LD contribution.This results in underestimates of h 2 snp , regardless of trait architecture.
One limitation of this paper is that we have focused on random-effects estimators of h 2 snp because of their widespread use.Estimators of h 2 snp can be broadly grouped into random-and fixed effect estimators based on how they treat SNP effects [35].Fixed effect estimators make fewer distributional assumptions but they are not as widely used because they require conditional estimates of all variants -a highdimensional problem where the number of markers is often far larger than the sample size [52].This is one reason why random effect estimators, such as GREML, are popular -because they reduce the number of parameters that need to be estimated by assuming that the effects are drawn from some distribution where the variance is the only parameter of interest.Fixed effects estimators, in principle, should be able to capture the LD contribution but this is not obvious in practice since the simulations used to evaluate the accuracy of such estimators still assume uncorrelated effects [35,52,53].Further research is needed to clarify the interpretation of the different estimators of h 2 snp in structured populations under a range of genetic architectures.

Does the LD contribution to the genetic variance have practical implications?
The answer to this depends on the context in which SNP heritability is used.ĥ snp can be useful in quantifying the power to detect variants in GWAS where the quantity of interest is the genic variance.But ĥ snp can lead to misleading conclusions if used to measure the extent to which genetic variation contributes to phenotypic variation, in predicting the response to selection, or in defining the upper limit of polygenic prediction accuracy [2] -applications where the LD contribution is important.
Ultimately, the discrepancy between ĥ snp and h 2 in practice is an empirical question, the answer to which depends on the degree of population structure (which we can measure) and the genetic architecture of the trait (which we do not know a priori ).We show that for most traits, the contribution of population structure to the variance explained by GWAS SNPs is modest among African Americans.Thus, if we assume that the genetic architecture of GWAS SNPs represents that of all causal variants, then despite incorrect assumptions, the discrepancy between ĥ snp and h 2 should be fairly modest.But this assumption is unrealistic given that GWAS SNPs are common variants that in most cases cumulatively explain a fraction of trait heritability.What is the LD contribution of the rest of the genome, particularly rare variants?This is not obvious and will become clearer in the near future through large sequence-based studies [54].While these are underway, theoretical studies are needed to understand how different selection regimes influence the directional LD between causal variants -clearly an important aspect of the genetic architecture of complex traits.

Simulating genetic architecture
We first drew the allele frequency (f 0 ) of 1,000 biallelic causal loci in the ancestor of populations A and B from a uniform distribution, U (0.001, 0.999).Then, we simulated their frequency in populations A and B (f A and f B ) under the Balding-Nichols model [32], such that where F = 0.2 is the inbreeding coefficient.We implemented this using code adapted from Lin et al.
(2021) [55].To avoid drawing extremely rare alleles, we continued to draw f A and f B until we had 1,000 loci with f A , f B ∈ (0.01, 0.99).
We generated the effect size (β) of each locus by sampling from N (0, , where m is the number of loci and f is the mean allele frequency across populations A and B. Thus, rare variants have larger effects than common variants and the total genetic variance sums to 1.Given these effects, we simulated two different traits, one with a large difference in means between populations A and B (Trait 1) and the other with roughly no difference (Trait 2).This was achieved by permuting the signs of the effects 100 times to get a distribution of V gb -the genetic variance between populations.This has the effect of varying the LD contribution without changing the F ST at causal loci.We selected the maximum and minimum of V gb to represent Traits 1 and 2.

Simulating admixture
We simulated the genotypes, local ancestry, and phenotype for 10,000 admixed individuals per generation under the hybrid isolation (HI) and continuous gene flow (CGF) models by adapting the code from Zaitlen et al. (2017) [26].We denote the ancestry of a randomly selected individual k with θ, the fraction of their genome from population A. At t = 0 under the HI model, we set θ to 1 for individuals from population A and 0 if they were from population B such that E(θ) = p ∈ {0.1, 0.2, 0.5} with no further gene flow from either source population.In the CGF model, population B receives a constant amount q from population A in every generation starting at t = 0.The mean overall proportion of ancestry in the population is kept the same as the HI model by setting q = 1 − (1 − p) 1 t where t is the number of generations of gene flow from A. In every generation, we simulated ancestry-based assortative mating by selecting mates such that the correlation between their ancestries is P ∈ {0, 0.3, 0.6, 0.9} in every generation.We do this by repeatedly permuting individuals with respect to each other until P falls within ±0.01 of the desired value.It becomes difficult to meet this criterion when V(θ) is small (Fig. 1C).To overcome this, we relaxed the threshold up to 0.04 for some conditions, i.e., when θ ∈ {0.1, 0.2} and t ≥ 50.We generated expected variance in individual ancestry using the expression in Zaitlen et al. (2017) [26].At time t under the HI model where P measures the strength of assortative mating, i.e, the correlation between the ancestry between individuals in a mating pair.Under the CGF model, We sampled the local ancestry at each i th locus as Bin(1, θ f ) and θ m and θ f represent the ancestry of the maternal and paternal chromosome, respectively.
The global ancestry of the individual is then calculated as , where m is the number of loci.We sample the genotypes x im and x if from a binomial distribution conditioning on local ancestry.
For example, the genotype on the maternal chromosome is where f A i and f B i represent the allele frequency in populations A and B, respectively.
Then, the genotype can be obtained as the sum of the maternal and paternal genotypes: We calculate the genetic value of each individual as g = m i=1 β i x i and the genetic variance as V(g).

Heritability estimation with GREML
We used the --reml and --reml-no-constrain flags in GCTA [5] to estimate σ 2 u and σ 2 v , the genetic variance due to genotypes and local ancestry, respectively.We could not run GCTA without noise in the genetic values so we simulated individual phenotypes with a heritability of h 2 = 0.8 by adding random noise . We computed three different GRMs, which correspond to different transformations of the genotypes: (i) standard, (i) Variance or V (x) scaled, and (ii) LD-scaled.
For the standard GRM, the genotypes at the i th SNP are standardized such that .
For the variance scaled GRM, we computed where V(x i ) is the sample variance of the genotypes at the i th SNP.The LD-scaled GRM conceptually corresponds to standardizing the genotypes by the SNP covariance, rather than its variance.Let X represent the n × m unstandardized matrix of genotypes and P represent an n × m matrix where the i th column contains the allele frequency of that SNP.Let U be the upper triangular 'square root' matrix of the m × m SNP covariance matrix Σ such that Σ = U U .Then, the standardized genotypes are computed as Z = (X − 2P )U −1 and the GRM becomes (X − 2P )Σ −1 (X − 2P ) [33].Similarly, the three GRMs for local ancestry were computed by scaling local ancestry with (i) 2γ i (1 − γi ) where we denote γi as the mean local ancestry at the i th SNP, or with the (ii) variance, or (iii) covariance of local ancestry, respectively.We estimated σ 2 u and σ 2 v with and without individual ancestry as a fixed effect to correct for any confounding due to genetic stratification.This was done by using the --qcovar flag.

Heritability estimation with HE regression
Haseman-Elston regression with and without ancestry correction was implemented using custom scripts in R [56].To estimate V g without ancestry correction, we first computed the cross-product of the centered phenotypes (y), resulting in an n × n matrix yy .We stacked the upper-triangular matrix of yy into a vector and regressed it on the corresponding elements of the GRM (ψ), taking the slope as an estimate of V g : To correct for individual ancestry, we followed the approach of Min et al. ( 2022) [35].To do this, we first regressed out the effect of individual ancestry (θ) on phenotype.The regression coefficient can be expressed as θ(θ θ) −1 θ and the residuals as y * = (I − θ(θ θ) −1 θ )y.Then, we fit the following model: where θθ represents the cross-product of individual ancestry, δ represents its corresponding regression coefficient, and V g represents the parameter of interest, i.e., the genetic variance and V e , the residual variance.
To demonstrate that the bias in HE estimates arises because of a bias in the estimate of LD contribution, not the genic variance, we carried out a simple simulation where half of the individuals in the population derive their ancestry from population A and the rest from population B. This is equivalent to the meta-population under the HI model at t = 0 where E(θ) = 0.5.We simulated genotypes for 1, 000 individuals for m = 100 loci where the allele frequencies in populations A and B were set to f A = 0.1 and f B = 0.8, respectively.We standardized the genotypes at each locus i using the square-root of the sample variance and assigned effect sizes such that the total genetic variance explained by all loci is equal to 1, i.e., the effect of the scaled genotype at the i th locus is u i = 1 √ m .This is equivalent to the effect size of the unscaled genotypes being where V(x i ) is the sample variance at the i th locus.
We introduced randomness in the direction of the effect by assigning a negative or positive sign to each locus uniformly at random 100 times to generate 100 traits with the same genic variance but varying LD contributions.Then, for each trait we computed the two terms in Eq. 2, which should converge to the genic variance and LD contributions, which represent the genic and LD components to the HE regression estimate.Fig. S5 shows that in the presence of directional LD, the overall bias is in the HE regression estimate is due to an exaggerated estimate of the LD contribution.

Estimating variance explained by GWAS SNPs
To decompose the variance explained by GWAS SNPs in African Americans, we needed four quantities: (i) effect sizes of GWAS SNPs, (ii) their allele frequencies in Africans and Europeans, and (iii) the mean and variance of global ancestry in African Americans (Equation 1).
We retrieved the summary statistics of 26 traits from GWAS catalog [40].Full list of traits and the source papers [44,[57][58][59][60][61][62][63][64] are listed in Table S1.To maximize the number of variants discovered, we chose summary statistics from studies that were conducted in both European and multi-ancestry samples and that reported the following information: effect allele, effect size, p-value, and genomic position.For birth weight, we downloaded the data from the Early Growth Genetics (EGG) consortium website [61] since the version reported on the GWAS catalog is incomplete.For skin pigmentation, we chose summary statistics from the UKB [65] released by the Neale Lab (http://www.nealelab.is/uk-biobank)and processed by Ju and Mathieson [44] to represent effect sizes estimated among individuals of European ancestry.We also selected summary statistics from Lona-Durazo et al. (2019) where effect sizes were meta-analyzed across four admixed cohorts [57].Lona-Durazo et al. provide summary statistics separately with and without conditioning on rs1426654 and rs35397 -two large effect variants in SLC24A5 and SLC45A2.
We used the 'conditioned' effect sizes and added in the effects of rs1426654 and rs35397 to estimate genetic variance.
We selected independent hits for each trait by pruning and thresholding with PLINK v1.90b6.21[66] in two steps as in Ju et al. (2020) [44].We used the genotype data of GBR from the 1000 genome project [39] as the LD reference panel.We kept only SNPs (indels were removed) that passed the genome-wide significant threshold (--clump-p1 5e-8 ) with a pairwise LD cutoff of 0.05 (--clump-r2 0.05 ) and a physical distance threshold of 250Kb (--clump-kb 250 ) for clumping.Second, we applied a second round of clumping (--clump-kb 100 ) to remove SNPs within 100kb.
When GWAS was carried out separately in different ancestry cohorts in the same study, we used inverse-variance weighting to meta-analyze effect sizes for variants that were genome-wide significant (p-value < 5 × 10 −8 ) in at least one cohort.This allowed us to maximize the discovery of variants such as the Duffy null allele that are absent among individuals of European ancestry but polymorphic in other populations [47].
We used allele frequencies from the 1000 Genomes CEU and YRI to represent the allele frequencies of GWAS SNPs in Europeans and Africans, respectively, making sure that the alleles reported in the summary statistics matched the alleles reported in the 1000 Genomes.We estimated the global ancestry of ASW individuals (N = 74) with CEU and YRI individuals from 1000 genome (phase 3) using ADMIXTURE 1.3.0[67] with k=2 and used it to calculate the mean (proportion of African ancestry = 0.767) and variance (0.018) of global ancestry in ASW.With the effect sizes, allele frequencies, and the mean and variance in ancestry, we calculated the four components of genetic variance using Equation 1and expressed them as a fraction of the total genetic variance.
Initially, the multi-ancestry summary statistics for a few traits (NEU, WBC, MON, MCH, BAS) yielded values > 1 for the proportion of variance explained.This is likely because, despite LD pruning, some of the variants in the model are not independent and tag large effect variants under divergent selection such as the Duffy null allele, leading to an inflated contribution of LD.We checked this by calculating the pairwise contribution , i.e., , of all SNPs in the model and show long-range positive LD between variants on chromosome 1 for NEU, WBC, and MON, especially with the Duffy null allele (Fig. S6A-C).A similar pattern was observed on chromosome 16 for MCH, confirming our suspicion.This also suggests that for certain traits, pruning and thresholding approaches are not guaranteed to yield independent hits.To get around this problem, we retained only one association with the lowest p-value, each from chromosome 1 (rs2814778 for NEU, WBC, and MON) and chromosome 16 (rs13331259 for MCH) (Fig. S6D).For BAS, we observed that the variance explained was driven by a rare variant (rs188411703, MAF = 0.0024) of large effect (β = −2.27).We believe this effect estimate to be inflated and therefore, we removed it from our calculation.
As a sanity check, we independently estimated the genetic variance as the variance in polygenic scores, calculated using --score flag in PLINK, [66] in ASW individuals.We compared the first estimate of the genetic variance to the second (Fig. S7) to confirm two things: (i) the allele frequencies, and mean and variance in ancestry are estimated correctly, and (ii) the variants are more or less independent in that they do not absorb the effects of other variants in the model.We show that the two estimates of the genetic variance are strongly correlated (r ∼ 0.85, Fig. S7).The 95% confidence intervals were calculated by sampling individuals with replacement 10,000 times.
where E θ represents the expectation taken over θ.
We derive V(x i |θ) by further conditioning on the local ancestry at each locus.
where E γ represents expectation taken over local ancestry.Since we are interested in the variance at a single locus, we will ignore the subscript i and denote the frequency of the trait-increasing allele in populations A and B with f A and f B , respectively.
To derive V{E γ (x|γ, θ)}, note that Note, that we can also express V (x i ) as: where the second term is the contribution of population structure to the genetic variance at locus i.
We can derive C(x i , x j ) using the law of total covariance: E θ {C(x i , x j |θ)} = 0 because we assume that the loci are unlinked and therefore, x i and x j are conditionally independent.Putting this all together, we get the genetic variance in admixed populations as presented in the main text: The only difference being that in the main text we use E instead of E θ for simplicity.With two 'unadmixed' source populations with equal number of individuals, E(θ) = 0.5 and V(θ) = E(θ){1−E(θ)} = 0.25 and V g reduces to: A3 The effect of genotype scale on Vg In the main text, we showed that both GREML and Haseman-Elston regression estimates of V g depend on how the genotypes are scaled.We provide an explanation of this behavior using the Haseman-Elston (HE) regression estimator, which is asymptotically equivalent to the GREML estimator if effects are uncorrelated [69] but which, unlike GREML, has a closed-form solution.First, let's assume a genetic architecture where all loci contribute equally to the genetic variance and there is no LD contribution.With the standard scaling, the genotype at a given locus i is z i = xi−2fi 2fi(1−fi) where f i is the frequency of the allele in the population.Under the random-effects model, this is equivalent to saying that the unscaled effects are: u being the parameter of interest.In a panmictic population, In an admixed population, to the genic variance The HE estimator of V g is based on the regression of products of (centered) phenotypes y k y l for all pairs of individuals k = l on the corresponding entries of the GRM (ψ) where ψ kl = m i=1 z ik z il m and z ik is the centered and scaled genotype of individual k for locus i: Where E kl represents the expectation over all k × l pairwise comparisons between individuals.It is simpler to express the HE estimator in terms of the scaled effects u i ∼ N (0, Where E(u i ) and E(u i u j ) represent expectations over random realizations of effect sizes.Thus, the last line follows from our assumption that the effect sizes are independent in expectation, i.e., E(u i u j ) = 0.
Note, that the estimate is still biased since it does not capture the contribution of population structure.

A3.1.2 Scaling by V(x i )
Next, we consider the case where the genotypes are standardized instead by the sample variance, i.e., such that V(z i ) = 1.We can derive E(u 2 i ) corresponding to this scaling by noting that the genetic variance is invariant under linear transformations of the genotype [14]: Then, the HE estimator becomes: Which provides an unbiased estimate of the genic variance.It's important to note that even though we assumed effect sizes under a random-effect model, the above result holds under a fixed-effect model as long as there is no directional LD.We discuss the implications of directional LD in the following section.

A3.2 Directional LD
Under the standard random-effect model, the effect sizes are assumed to be independent in expectation.
We discussed in the main text how certain processes (e.g.selection and assortative mating) can induce directional LD across causal loci.But directional LD might arise even for neutral traits and under the random-effects model for any given realization of effects.This can lead to biases both HE and GREML estimates of V g , though the direction and reason for bias is different for the two methods.GREML does not have a closed-form solution so the exact estimand is difficult to derive.Instead, we develop some intuition based on HE regression.

A3.2.1 Scaling by V(x i )
To do this, let u = [u 1 , u 2 , ..., u m ] represent the vector of a given realization of (fixed) effects corresponding to the standardized genotypes such that each locus contributes equally to σ 2 u , the genic variance, i.e., Let there be positive LD across loci such that all cross-product terms are u i u j = σ 2 u m .Then, the genetic variance explained by all loci is: where C(z i , z j ) is the LD between the i th and j th loci that ranges from 0 (no LD) to 1 (perfect LD).
Thus, the LD contribution to V g ranges from 0 to (m − 1)σ 2 m .In comparison, the HE estimator is: This shows that the bias due to directional LD in the HE estimate of V g does not come from the genic, but the LD component.When there is no LD, e.g. if the population has reached equilibrium after generations of random mating, this component goes to zero and both the estimate and V g converge to the same value -the genic variance.The LD component is maximum when the i th and j th loci are in perfect LD.In this case, i and j are exchangeable and the second term of the estimator reduces to (m − 1)σ 2 m .Thus, HE regression should give an unbiased estimate of V g , even in the presence of directional LD, but only when LD is perfect.For any other value 0 < C(z i , z j ) < 1, the estimate is biased (Fig. A1).An interpretable, analytical derivation of the second term in Eq. 4 is complicated but we illustrate the bias with simulations below.
For unlinked markers, C(z i , z j ) is a function of 4 V(θ)(f

A3.2.2 Scaling by LD
In the main text, we showed that standardizing the genotypes at a locus by its covariance with other loci accounts for the bias for GREML and HE estimators.More specifically, the 'LD-scaled' genotypes can be written as Z = (X − 2P )U −1 where P is an n × m matrix such that all elements of the i th column contain the frequency of the i th SNP and U is the (upper triangular) square root of the LD matrix, i.e, Σ = U U .Under this scheme, the standardized genotypes are uncorrelated and therefore, the second term in Eqs.Which captures both the genic and LD contributions and therefore, provides an unbiased estimate of V g .

A4 Genetic variance after correction for individual ancestry
In the main text, we stated that including individual ancestry as a fixed effect in GREML can lead to an underestimate of V g in the presence of population structure.Mixed effect models deal with fixed effects (ancestry in our case) by projecting them out of the phenotypes, and estimating the residual variance.This is conceptually equivalent to measuring the residual variance of the regression between phenotype and ancestry.As a result, any variance in the phenotype that is explained by ancestry is removed.To

Figure 1 :
Figure 1: The behavior of mean and variance of individual ancestry as a function of admixture history.(A) Shows the demographic models under which simulations were carried out.Admixture might occur once (Hybrid Isolation, HI, left column) or continuously (Continuous Gene Flow, CGF, right column).(B) The mean individual ancestry, E(θ) remains constant over time in the HI model and increases in the CGF model with continued gene flow.(C) The variance in individual ancestry, V(θ) is maximum at t = 0 in the HI model, decaying subsequently.V(θ) increases with gene flow in the CGF model and will subsequently decrease with time.P measures the strength of assortative mating, which slows the decay of V(θ).P=0.6 is missing for simulations run for 50 and 100 generations and θ ∈ {0.1, 0.2} due to the difficulty in finding mate pairs (Methods).

Figure 2 :
Figure 2: Decomposing genetic variance in a two-population system.The plot illustrates the expected distribution of genetic values in two populations under different selective pressures and the terms on the right list the total (V g ) and between-population genetic variance (V gb ) expected over the evolutionary process.For neutrally evolving traits (top row), we expect there to be an absolute difference in the mean genetic values (| ḡA − ḡB |) that is proportional to F ST .For traits under divergent selection (middle),

Figure 3 :
Figure 3: Genetic variance in admixed populations under the (A) HI and (B) CGF models.Dotted lines represent the expected genetic variance based on Eq. (1) and solid lines represent results of simulations averaged over ten replicates.Red and blue lines represent traits under divergent and stabilizing selection, respectively.P = 0.6 is missing for simulations run for 50 and 100 generations and θ ∈ {0.1, 0.2} due to the difficulty in finding mate pairs (Methods)

Figure 4 :
Figure 4: The behavior of GREML estimates of the genetic variance (σ 2 u ) in admixed populations under the HI (left column) and CGF (right column) models either without (A-D) or with (E-H) individual ancestry as a fixed effect.The solid lines represent estimates from simulated data averaged across ten replicates with red and blue colors representing estimates for traits under divergent and stabilizing selection, respectively.P indicates the strength of assortative mating.The shaded area represents the 95% confidence bands generated by bootstrapping (sampling with replacement 100 times) the point estimate reported by GCTA.The dotted lines either represent the expected variance in the population based on Eq. 1 (A & B) or the expected estimate for three different ways of scaling genotypes (B-D &

Figure 5 :
Figure 5: Genetic variance ( Vg ) estimated with HE regression in admixed populations under the HI (left column) and CGF (right column) models either without (A-C) or with (D-F) adjustment for individual ancestry.The solid lines represent estimates from simulated data averaged across ten replicates with red and blue colors representing estimates for traits under divergent and stabilizing selection, respectively.P indicates the strength of assortative mating.(A & D) show the behavior of Vg for the default scaling,(B, E) shows Vg when the genotype at a locus is scaled by its sample variance (V(x) scaled), and (C, F) when it is scaled by the sample covariance (LD scaled).The dotted lines in A-E represent the expected V g in the population based on Eq. 1 and in F, represent the expected V g after removing any genetic variance along the ancestry axis.The shaded areas represent the 95% bootstrapped confidence bands of the estimate.

Figure 6 :
Figure 6: The behavior of GREML estimates of the variance due to local ancestry (σ 2 v ) in admixed populations under the HI (left column) and CGF (right column) models either without (A-D) or with (E-H) individual ancestry included as a fixed effect.The solid lines represent estimates from simulated data averaged across ten replicates with red and blue colors representing estimates for traits under divergent and stabilizing selection, respectively.P indicates the strength of assortative mating.The dotted lines either represent the expected variance in the population (A & B) or the expected estimate for three different ways of scaling local ancestry (B-D & F-H).(A-B & E-F) show the behavior of σ2 v for the default scaling, (C, G) shows σ2 v when local ancestry is scaled by the sample variance, and (D, H) when it is scaled by the sample covariance.Shaded regions represent the 95% confidence bands.Some runs in (D & H) failed to converge as seen by the missing segments of the solid lines because the expected variance in such cases was too small.

Figure 7 :
Figure 7: Decomposing the genetic variance explained by GWAS SNPs in the 1000 Genomes ASW (African Americans from Southwest).We calculated the four variance components listed in Eq. 1, their values shown on the y-axis as a fraction of the total variance explained (shown as percentage at the bottom).The LD contribution, which can be positive or negative, is shown in yellow.The number of variants used to calculate variance components for each trait is also shown at the bottom.

σ 2 u 2 uj
m z ik z il m w=1 z wk z wl /m E kl ( m i=1 z ik z il /m m w=1 z wk z wl /m) m z ik z jl m w=1 z wk z wl /m E kl ( m i=1 z ik z il /m m w=1 z wk z wl /m) =i z ik z jl z wk z wl E kl ( m i=1 z ik z il m w=1 z wk z wl )

A i +f B i 2 = 0 . 5 .
1).Perfect LD arises when (i) both f A i − f B i = 1 and f A j − f B j = 1 and (ii) V(θ) is maximum, which occurs at the time of admixture when source populations mix equally, i.e, E(θ) = 0.5.To generate a range of LD, we simulated an admixed population (N = 1, 000) with equal number of individuals from populations A and B. Thus,4 V(θ) = 4 E(θ){1 − E(θ)} = 1.We simulated genotypes for each individual at 50 'causal' loci where the difference between the frequencies in the source populations,f A i − f B i ∈ [0,1] with the condition that f We assigned each locus the same effect size (on the variance-standardized scale) of +1/ √ m summing up to a genic variance of 1.The positive sign ensures positive LD across loci, i.e, all off-diagonal elements of uu are set to 1/m.For each simulation, we computed the expected and estimated LD component using the second terms in Eqs. 3 and 4, respectively, and averaged the results over 100 replications.

Figure A1 :
Figure A1: The behavior of the LD contribution (y-axis) to the genetic variance (red) and the Haseman-Elston regression estimate (blue) as a function of LD, i.e.C(z i , z j ) (x-axis).Each point represents the contribution calculated from a random draw of genotypes, given C(z i , z j ) ∝ (f A i − f B i )(f A j − f B j ).The red line represents the expected LD contribution and the black dashed line represents the contribution expected in the case of perfect LD.

β
3 and 4 are zero.This reduces the estimator to the first term, representing the sum of squares of effect sizes, i.e. u u = m i=1 u 2 i .The effect sizes corresponding to the LD scaled genotypes are u = U β and the sum of squares is:36u u = (U β) (U β) = β U U β = βΣβ = i β j C(x i , x j )

Figure S2 :
Figure S2: The behavior of the genetic variance due to local ancestry in admixed populations under the (A) HI and (B) CGF models.The solid lines represent values observed in simulations averaged across ten replicates and the dotted lines represent the expected values based on Eq. 1 of the main text.The red and blue lines represent values for traits 1 and 2, respectively.P indicates the strength of assortative mating.P=0.6 is missing for simulations run for 50 and 100 generations and θ ∈ {0.1, 0.2} due to the difficulty in finding mate pairs (Methods).

Figure S3 :
Figure S3: The behavior of GREML estimates of SNP heritability ( ĥ snp ) in admixed populations under the HI (left column) and CGF (right column) models either without (A-C) or with (D-F) individual ancestry as a fixed effect.The solid lines represent ĥ snp averaged across ten replicates, with red and blue colors representing estimates for traits under divergent and stabilizing selection, respectively.(A, D) show the behavior of ĥ snp for the default scaling, (B, E) shows ĥ snp when the genotype at a locus is scaled by its sample variance (V(x) scaled), and (C, F) when it is scaled by the sample covariance (LD scaled).The shaded area represents the 95% confidence bands generated by bootstrapping (sampling with replacement 100 times) the point estimate reported by GCTA.The black dotted lines represent the expected heritability value given the simulation settings (h 2 = 0.8).P indicates the strength of assortative mating

Figure S5 :
Figure S5: The effect of directional LD on Haseman-Elston estimate of genetic variance (V g ).Each individual point is an independent simulation where the effects were drawn from a normal distribution and applied to genotypes from an admixed population (Methods).The solid red line shows the y = x line and the color of each point represents the contribution of LD to V g .

Figure S6 :
Figure S6: The LD contribution to the variance explained by variant pairs for (A) neutrophil counts (NEU), (B) white blood count (WBC), (C) monocyte count (MON), and (D) mean corpuscular hemoglobin (MCH).Only chromosomes where we suspected there was a disproportionate contribution to the variance explained are shown.

Figure S7 :
Figure S7: Expected variance explained estimated using Equation 1 (x-axis) vs the observed variance in polygenic scores (y-axis) in ASW are strongly correlated.Confidence intervals were generated by non-parametric bootstrap (Methods).