Assortative mating biases marker-based heritability estimators

Many traits are subject to assortative mating, with recent molecular genetic findings confirming longstanding theoretical predictions that assortative mating induces long range dependence across causal variants. However, all marker-based heritability estimators implicitly assume mating is random. We provide mathematical and simulation-based evidence demonstrating that both method-of-moments and likelihood-based estimators are biased in the presence of assortative mating and derive corrected heritability estimators for traits subject to assortment. Finally, we demonstrate that the empirical patterns of estimates across methods and sample sizes for real traits subject to assortative mating are congruent with expected assortative mating-induced biases. For example, marker-based heritability estimates for height are 14% – 23% higher than corrected estimates using UK Biobank data.

P ositive primary phenotypic assortative mating (hereafter simply "AM"), the phenomenon whereby mate-choice is based on phenotypic similarity, has been observed for a variety of heritable traits in human and non-human animals [1][2][3][4] . A century ago, Fisher demonstrated that AM induces long-range positive correlations between trait-increasing allele counts at causal loci across the genome, thereby increasing genetic variance across successive generations until it approaches a stable equilibrium 5 . Since Fisher's time, it has been established that many human traits are subject to AM, and that estimates of genetic and environmental variance from twin and family designs, which assume random mating, can be biased in the presence of AM 6 . Recently, there has been increased focus on how AM can influence results based on measured genetic data, in part due to the recognition that AM may contribute to observed discrepancies between population-based versus family-based estimates in genome-wide association studies (GWAS) [7][8][9][10] and Mendelian randomization studies 11,12 . However, despite hundreds of publications reporting estimates of SNP-heritability, no study has yet investigated how AM influences these estimates. Moreover, many benchmark traits central to the scientific discourse regarding the "still-missing heritability"-for example, height and educational attainment-are precisely those for which phenotypic and genetic data is consistent with AM 4,13-15 , raising the possibility that conclusions about SNP-heritability could be biased in systematic but as yet poorly understood ways.
Here, we characterize the impact of AM on the two families of marker-based heritability estimators, which we generically denoteĥ 2 SNP . The first uses the method of moments (MoM) and is typified by univariate Haseman-Elston (HE) regression (ĥ 2 HE ) 16 but includes PCGC regression 17 and linkage disequilibrium (LD) score regression 18 . The second uses residual maximum likelihood (REML) 19 to estimate heritability (ĥ 2 REML ) and is implemented in software such as by GCTA 20 and BOLT-REML 21 . We assume Fisher's classical model of AM, which describes the equilibrium properties of a heritable trait for which mates' genotypes are conditionally independent given the heritable components of their phenotypes and for which offspring environments are independent of parental phenotypes. This model has formed the theoretical foundation for recent investigations of AM using measured genetic data 4,14,22 . We provide mathematical and simulation-based arguments demonstrating that AM induces a modest yet non-negligible bias in both classes of estimators that is not addressed by conventional methods of accounting for population structure. In the process, we extend results in random matrix theory and classical quantitative genetics by characterizing the higher-order moments of causal variants and the limiting spectral distribution of the genomic relatedness matrix (GRM) under AM, which clarifies why genomic principal components fail to capture the effects of AM. Additionally, we provide empirical results using data from the UK Biobank that are congruent with our theoretical predictions. Finally, we provide guidelines for estimating and interpretingĥ 2 SNP when traits are subject to AM.

Results
Our theoretical results depend on several key parameters: r denotes the phenotypic correlation between mates on an additive phenotype y with a heritable component Zu; h 2 0 denotes the panmictic heritability, what the narrow-sense heritability of the phenotype would be in the absence of AM; h 2 1 denotes the equilibrium narrow-sense heritability approached under multiple generations of AM; and Z denotes the n m matrix of n unrelated individuals' standardized genotypes at m causal loci with effects vector u. We initially assume that all causal variants are present in Z (i.e., that the heritability is fully explained by measured variants); the reason for this is that the dynamics of AM across generations depend on the true heritability of the trait, irrespective of the fraction of genetic variation tagged by SNPs. We later employ simulations to demonstrate that our major qualitative conclusions remain valid when this assumption is relaxed.
The rows of Z (individuals' genotypes) are independent random vectors with m m covariance matrix ϒ, which quantifies the correlation between loci. Under random mating, ϒ:¼ϒ 0 is approximately block diagonal such that causal variants are largely (aside from LD between nearby variants) stochastically independent. However, under equilibrium AM, ϒ:¼ϒ 1 is dense due to the presence of positive long-range correlations among traitincreasing allele counts within and across chromosomes. As the elements of ϒ 1 agree in sign with the corresponding elements of uu T (i.e., trait increasing alleles are positively correlated), the equilibrium additive genetic variance under AM, σ 2 g;1 , is considerably greater than the panmictic additive genetic variance, Fig. 1 and Supplementary Notes 1, 2).
Haseman-Elston regression estimates under AM. We first characterize the influence of AM onĥ 2 HE , perhaps the simplest marker-based MoM heritability estimator. Let e y denote the standardized phenotype.ĥ 2 HE is the slope of the subdiagonal elements of the phenotypic outer product, e ye y T , regressed on the subdiagonal elements of the GRM, m À1 ZZ T . Assuming that all genetic variance is explained by measured variants, we establish the following general result (see Supplementary Note 3): Under AM, the bracketed quantity is greater than one, and thusĥ 2 HE is upwardly biased relative to both h 2 0 and h 2 1 . Intuitively, this bias occurs because AM influences the outcome (e ye y T ) differently than the predictor (m À1 ZZ T ). The AM-induced positive correlations between all causal variants inflate genetic (co-)variance, and this inflation is accurately reflected in e ye y T . However, the elements of m À1 ZZ T represent the average correlation between pairs of individuals at homologous loci taken one at a time. Thus, the elements of m À1 ZZ T are governed almost completely by the correlations within loci; correlations between different loci have almost no influence. Because the vast majority of the inflation in genetic (co-)variation is caused by the correlations between and not within loci, AM on phenotype has a negligible impact on the GRM for polygenic traits in a genetically homogenous population. This latter point is itself noteworthy because some studies have claimed that inflated genomic similarity between friends or mates 22 is evidence for phenotypic assortment, but the actual increase is trivial (of order O m À1 À Á relative to the increase in genetic variance) and such observations have more probable explanations, such as imperfectly controlled stratification 23 .
Under the stricter assumption of exchangeable loci (i.e., each causal variant explains equal phenotypic variance), we derive the following approximate expression, dependent only on r and h 2 0 (see Supplementary Note S3): Large-scale simulations using realistic genotype data (see Online Methods) with p ¼ 10 6 SNPs of which m ¼ 10 4 were causal with effects u $ i:i:d: Nð0; σ 2 g;0 =mÞ demonstrate that these approximations are accurate even when the exchangeable loci assumption is violated (Figs. 1,2). Results were nearly identical when setting m ¼ 10 5 causal variants. In addition, our simulations confirmed that LDSC, which is mathematically equivalent to HE regression when LD scores are exact 24 , is similarly biased upwards under AM (Fig. 3a). However, the impact of this bias in real world applications depends not only on the extent of AM, but also on the degree to which estimated LD scores reflect the true LD structure in a given population.
Residual maximum likelihood estimates under AM. In contrast toĥ 2 HE , the value ofĥ 2 REML changes as a function of sample size under AM. Using the same simulation approach described above, we found thatĥ 2 REML exhibits upward biases similar in magnitude to those ofĥ 2 HE in small samples (e.g., n % 10;000). In large samples (e.g., n > 200;000), however,ĥ 2 REML drops below h 2 1 (Fig. 2a). More formally, under the assumption of exchangeable loci, we prove that  Fig. 1). In essence, the parameter values that maximize the residual likelihood function depend only on the eigenvalues of the GRM, and the long-range dependence among causal variants induced by AM is "weak" in the sense that the distributions of the eigenvalues (i.e., spectral distributions) of the GRM under random mating and under AM are asymptotically equivalent. Thus, in large samples,ĥ 2 REML converges to what the heritability would be if the causal variants were independent, i.e., to h 2 0 . However, in finite samples, we show via simulation that this convergence can be extremely gradual, requiring samples approaching millions of individuals beforeĥ 2 REML approaches h 2 0 (Fig. 2b).

Conventional methods do not mitigate AM-induced bias.
Inclusion of ancestral PCs as covariates fails to mitigate the AMinduced bias in both the MoM and the REML estimates (Fig. 3a). Indeed, we demonstrate that AM has a negligible effect on the spectral distribution of the GRM in large samples with many variants (see Supplementary Note 2). Similarly, these biases are not mitigated by modeling multiple genetic variance components that partition SNPs according to LD score and minor allele frequency (Fig. 3a).
AM-induced bias persists when not all causal variants are measured. In real-world applications, measured genotypes will often include only a fraction of the total number of causal variants (or at least good proxies in high LD). Previous work has demonstrated that the extent to which marker data imperfectly tag causal variants will biasĥ 2 REML estimates in random mating populations 25,26 ; these circumstances are further complicated under AM as imperfect tagging will not only corrupt the relationship between a given marker and proximal causal variants, but between a given marker and all other causal variants. To assess the impact of AM when not all variants are measured, we comparedĥ 2 HE andĥ 2 REML in simulated data that randomly dropped π = 0, 50, or 75% of all p ¼ 10 6 variants, including both causal (m ¼ 10 4 ) and non-causal (p À m ¼ 990;000) SNPs (Fig. 3b). As expected, this resulted in attenuated estimates commensurate with, but not proportional to the fraction of missing data: on average across methods,ĥ 2 SNP was 8.4% (s.e. 0.009%) and 17.6% (s.e. 0.009%) lower after randomly dropping 50% and 75% of SNPs, respectively (the degree of attenuation was smaller than the proportion of SNPs dropped due to LD between retained SNPs and discarded causal variants). Though future work is required to fully characterize the relationship between heritability estimates and the degree of missing heritability under AM, the pattern of bias due to AM when some of the heritability is missing is qualitatively similar to that when all causal variants are present.
The differential influence of sample size onĥ 2 HE andĥ 2 REML in real traits. To investigate whether our theoretical predictions are observable in real data, we examined the relationship between sample size andĥ 2 SNP in a sample of 335,551 unrelated Europeanancestry individuals in the UK Biobank 27 . We a priori selected four phenotypes based on evidence (height, years of education) or lack of evidence (body mass index, bone mineral density) for primary phenotypic AM in a previous study 14 . We then computedĥ 2 HE andĥ 2 REML in pairs of small (n small = 16,000) versus large (n large = n-16,000) disjoint subsamples, where n denotes the total available sample size for each phenotype. Consistent with theory and simulation results, the effect of sample size onĥ 2 REML was significantly different from the effect of sample size onĥ 2 HE for height (p = 5.24e−4) and for years of education (p = 3.94e−4); this was not observed for body mass index (p = 0.094) or for bone mineral density (p = 0.302; see "Methods" section, Fig. 4). These We computed multiple estimates per sample size for each estimator and parameter combination by applying estimators to independent sub-samples. Whereas HE regression estimates are upwardly biased independent of sample size, REML estimates slowly converge to the panmictic heritability as sample sizes increase. b Extended simulations demonstrating high-dimensional behavior of the REML estimator as a function of sample size. Forward time simulations required a larger population size (N sim ¼ 3 10 6 ) to obtain samples of up to n ¼ 648;000 unrelated individuals. Obtaining REML estimates for samples larger than this was not computationally feasible, but the dashed red line shows predicted values for larger sample sizes extrapolated from a regression model including first and second order log-linear components. Results are consistent with theoretical predictions that the REML estimator converges to the panmictic heritability in very large samples.
height results are congruent with those of Mandal and colleagues, who independently observed thatĥ 2 REML decreased with increasing sample size 28 . Unlike previous approaches which quantified the degree of AM for traits by correlating polygenic scores between mates 4 or within-individuals but between chromosomes 14 , our approach does not require the calculation of polygenic scores and is agnostic with respect to variant direction and effect size. Thus, calculatingĥ 2 REML across subsamples of varying sizes provides an alternative and independent way to detect genomic signatures of AM.
Interpretation ofĥ 2 SNP in the presence of AM. As we have shown, AM complicates the interpretation ofĥ 2 SNP . This complication is due not only to the differential effect AM has onĥ 2 HE andĥ 2 REML as a function of sample size, but also because AM changes the true heritability in the population. This means that h 2 SNP can be compared to two different population parameters: the equilibrium heritability h 2 1 , and the panmictic heritability h 2 0 . We inputĥ 2 HE and previously reported spousal correlations for height 2 and educational attainment 4  Although these estimates of h 2 1 and h 2 0 rely on the assumption that AM has reached equilibrium, theoretically-sound bounds for h 2 can be derived for traits where AM is at disequilibrium. Specifically, the current heritability for a population having undergone t generations of AM, h 2 t , is bounded in expectation from above byĥ 2 HE , with equality at t ¼ 0 (panmixis), and from below byĥ Additionally, simulations confirm that LDSC is subject to equivalent biases. Single component: standard single genomic variance component models. Single comp. + 10 PCs: included the first ten within-sample PCs as covariates. Partitioned: included four annotation-based variance components generated by median splits of within-sample minor allele frequencies and LD scores. b Simulations demonstrate that conclusions regarding estimator bias do not change when some of the influence of causal variants is not captured by measured SNPs. Simulations employed the same parameters as above except that 0, 50, or 75% of randomly selected SNPs (both causal and non-causal) were dropped. As expected, estimates were attenuated when SNPs were dropped but overall patterns remained consistent.  (Fig. 1).

Discussion
Summary of findings. In the last decade, there has been much interest in the interpretation ofĥ 2 SNP and the factors that can influence it 29 . Such factors include the number of causal variants and their typical levels of LD compared to the LD of SNPs used in the analysis 30,31 , the relationship between causal variant effect sizes and minor allele frequencies 32 , and passive gene-environment correlation arising from population stratification 33 or genetic nurture 7 . Despite this activity and the recent evidence corroborating long-standing theoretical expectations that AM alters the genetic architecture of heritable traits 14 , the effects of AM onĥ 2 SNP have remained a matter of speculation 9,30,34 . In the present investigation, we clarify how AM influencesĥ 2 SNP and show that these influences are different for MoM versus REML estimators as a function of sample size. By characterizing the full joint distribution of causal variants at equilibrium, we prove that REML produces a consistent estimator of h 2 0 (the heritability in the ancestral random-mating population), not h 2 1 (the present generation heritability for a trait at equilibrium), in large samples (see Supplementary Notes 2 and 3). However, under AM,ĥ 2 REML behaves oddly in finite samples. It is higher than h 2 1 in sample sizes typical of those published in the literature but drops below h 2 1 in sample sizes that are large by current standards (e.g., n > 200;000). On the other hand, MoM estimators yield estimates that are biased upwards with respect to both h 2 1 and h 2 0 and remain stable across sample sizes. Using UK  Spousal correlations (r) as previously reported in British cohorts 2,4 and heritability estimates for selected UK Biobank phenotypes (n denotes sample size). Height and years of education were selected a priori on the basis of previous evidence for primary phenotypic AM whereas body mass index and bone mineral density were selected a priori as negative controls. Assuming equilibrium, the equilibrium (ĥ 2 ) and panmictic (ĥ Biobank data, we observed this differential behavior ofĥ 2 HE and h 2 REML for two traits with previous evidence of primary phenotypic AM but not for two negative control traits. We further showed how the two population parameters of interest, h 2 1 and h 2 0 , can be estimated under the assumption of equilibrium given spousal correlations andĥ 2 HE , and that, at disequilibrium, the likely ranges of h 2 in the current generation and h 2 0 can also be estimated. Implications. The impact that AM has onĥ 2 SNP may often be negligible simply because AM is very low or nonexistent for many traits. For other traits, however, this does not appear to be the case, including "benchmark" traits like height and educational attainment as well as many behavioral and psychiatric traits. For example, the tetrachoric correlations between Swedish spouses for attention deficit hyperactive disorder (r = 0.31), autism spectrum disorder (r = 0.28), schizophrenia (r = 0.26), and substance use disorder (r = 0.30) are all higher than those typically observed for height 35 . Furthermore, many social attitudes such as conservatism-liberalism (r = 0.61) and religiosity (r = 0.71) exhibit extremely high spousal correlations, and available evidence is more consistent with a primary phenotypic model than social homogamy or convergence 36  for traits subject to AM when characteristics of the populations, samples, or designs differ. First, the effect that AM has on causal variant correlations is rapid-about half of its effect occurs after a single generation and it approaches equilibrium within 5-15 generations; as such, both the true heritability and the extent to whichĥ 2 SNP is biased will differ across populations that vary in strength or duration of AM. Second, because the expectation of h 2 REML decreases monotonically in n, comparingĥ 2 REML across samples of different sizes can create artifactual discrepancies. Third, though bothĥ 2 REML andĥ 2 HE are expected to increase as the fraction of causal variants included in the model increases, the relative bias ofĥ 2 REML will further increase if the inclusion of additional markers comes at the cost of decreased sample size. Therefore, the difference in REML estimates from a smaller sample of whole-genome sequence data (which ostensibly includes previously missing causal variants) versus those derived from a larger sample of SNP array data will appear artificially greater, though we note that additional factors, such as differences in the extents to which array versus sequence data are corrupted by subtle forms of population structure, might also contribute to such discrepancies. Fourth, our results demonstrate that the controlling for genomic principal components does not mitigate AM-induced biases; indeed, the asymptotic distributions of the eigenvalues of the GRM under random versus assortative mating are equivalent. As AM is in essence a form of population structure, our results reveal that our current understanding of the extent to which residual population structure might impact estimators of quantitative genetic parameters is incomplete. Finally, our results suggest that AM complicates comparisons ofĥ 2 SNP to estimates of heritability derived from family-based methods (mostly from twins, henceĥ depending on the design, biasing it downwards in "ACE" models that estimate additive genetic and shared environmental variance and upwards in "ADE" models that estimate additive genetic and dominance genetic variance. Together with the previously discussed dependence ofĥ 2 SNP on method and sample size, this further problematizes efforts to quantify the "still-missing" heritability (ĥ 2 twin Àĥ 2 SNP ).
Limitations and future directions. There are several limitations of the current approach. First among these are assumptions inherent to the primary phenotypic AM model. Some of these assumptions, including equilibrium and constancy of the phenotypic mating correlations across generations, aid with mathematical tractability but are to some extent inessential to the resulting phenomena. For example, while the problem of characterizing the joint distribution of causal variants is more difficult in a population subject to disequilibrium AM, we observed that estimators behave in a similar, albeit less extreme, fashion relative to their behavior in an equilibrium population. Other assumptions, such as the absence of gene-environment correlation (which may occur, for example, in the presence of vertical transmission) and the conditional independence of mates' genotypes given the heritable components of their phenotypes (which may be violated in structured populations), are more difficult to evaluate and deserve consideration in future investigations. For instance, when environmental influences are transmitted from parents to offspring, the dependence between causal variants will extend to non-heritable factors such that, in subsequent generations, all trait-increasing allele-counts and trait-increasing environmental factors will become positively correlated, and the linear mixed model assumed by both the REML and MoM estimators will be further misspecified.
We have also not investigated the impact that AM has on polygenic scores or within-family classes of estimators. For example, LD-score regression can be estimated based upon regression estimates from within-sibling analyses. The estimated genetic variance from such an analysis would be that in the base population (because it is based on the segregation variance), but the estimated phenotypic variance would be that in in the current population. This would also lead to a bias in heritability estimation (and similarly for heritability estimates using relatedness disequilibrium regression 10 ), as was pointed out in ref. 38 . This bias is different than those described here in that it is lower than either the base population or equilibrium population heritability. Additionally, we have not found closed-form solutions toĥ 2 REML that would allow us to correct REML estimates in the same way we have corrected HE estimates under the assumption of equilibrium AM. Further, we have not investigated the impact that natural selection has on the long-range correlations between causal variants and SNPs and what impact this might have onĥ 2 SNP . Finally, we have not considered the effects of multivariate AM or the impact of AM on marker-based genetic correlation estimators, which will also be misspecified under cross-trait AM. As such, these results should provide motivation and a starting place for the development of new methods that can provide unbiased estimates of genomic variance in the presence of AM and other factors that can influence the long-range correlations between SNPs.

Theoretical framework
The primary phenotypic assortment model. Here we introduce the model of AM as proposed by Fisher 5 and further developed by Nagylaki and others 39,40 (see Supplementary Note 1 for a detailed exposition). Briefly, we consider a phenotype as a random vector composed of independent heritable and non-heritable components: where the rows of Z, representing individuals' standardized genotypes, are independent m-dimensional random vectors following a multivariate discrete distribution with finite moments and finite support and which we assume are independent under panmixis. The vector of allele substitution effects u, which we treat as fixed, is such that u T u ¼ σ 2 g;0 . Further, we assume that 1. parent-parent-offspring trios' phenotypes are jointly Gaussian; 2. the phenotypic correlation between mates, r is constant across generations; and 3. there exists c 0 2 0; 1 ð Þsuch that max k¼1; ;m u k ≤ c 0 Á m À1=2 ; that is, as traits become increasingly polygenic, the maximal variance attributable to individual variants decreases commensurately.
The equilibrium distribution of causal variants. Over successive generations, the correlation between mates' phenotypes induces positive correlations across trait increasing allele counts independent of physical position on the genome and thereby increases the total genetic variance of the trait. The genetic variance rapidly approaches a stable equilibrium after several generations (typically within ten generations), at which point the within-individual and cross-mate correlations among causal variants are equal to one another. Using the results of Nagylaki 39 , we can express the equilibrium covariance matrix between causal variants as a low rank perturbation of a diagonal matrix of the form: where ϕ is a known vector-valued function of the substitution effects and mating correlation (see Supplementary Note 1.2) with elements ϕ k ¼ Oðm À1=2 Þ uniformly.
Higher order moments and the limiting spectral distribution of GRM. Employing tools from the study of thermodynamic equilibria, we extend these classical results to bound moments of higher orders (see Supplementary Note 2). Using these results, we extend the widely-known Marčenko-Pastur theorem, which describes the limiting distribution of the spectrum of sample covariance matrices corresponding to random matrices with independent sub-Gaussian elements 41 , to the case of random matrices with independent rows meeting particular moment conditions (see Supplementary Note 2). Together, these results establish the limiting spectral distribution of the sample GRM (i.e., the distribution of the eigenvalues of the sample GRM as both sample size and polygenicity increase) under AM, providing the necessary theoretical foundation to characterize the asymptotic behavior of the REML estimator. Further, these results explain why controlling for principal components fails to remove AM-induced biases: the impact of AM on the spectrum of the GRM is asymptotically negligible.
Haseman-Elston regression under AM. The HE regression estimator 42 is obtained by regressing the subdiagonal elements of the standardized phenotypic outer product e ye y T on the subdiagonal elements of the GRM m À1 ZZ T . Whereas elements of the outcome (the phenotypic outer product) reflect the dependences among all pairs of causal loci: elements of the GRM only capture the dependences within loci: As a result, the variance of the outcome increases whereas the variance of predictor remains largely unaffected, leading to overestimation of h 2 SNP ( Fig. 1; see Supplementary Note 3 for further details and proof).
REML and the spectrum of the GRM under AM. The REML estimator 19 models the phenotype as a random vector with marginal distribution, where X is an n c matrix of covariates with fixed effects β and the covariance structure is comprised of a heritable component (σ 2 g times the GRM) and a nonheritable component (σ 2 e times the identity). The heritability estimatorĥ 2 REML 1 σ 2 g =ðσ 2 g þσ 2 e Þ is derived by finding the values of the variance components that satisfy the equation, where l denotes the marginal log likelihood of the transformed random variable A T y for A T : R n ! ðcol XÞ ? R nÀc , A T A ¼ I. The conditional expectation of ∇l given the genotypes is a function of the eigenvalues of the GRM and, as a result, the asymptotic behavior ofĥ 2 REML is governed by the asymptotic distribution of the eigenvalues of m À1 ZZ T . A foundational result in random matrix theory states that for zero-mean unit-variance sub-Gaussian random matrices W 2 C n m with independent elements, the empirical spectral distribution function of m À1 WW T converges almost surely to the Marčenko-Pastur distribution 41 . Employing this result, Jiang and colleagues 43 demonstrated that, in the case of independent causal variants, REML consistently estimates the true heritability in high dimensional settings and is robust to certain forms of model misspecification. In Supplementary Note 2, we demonstrate that even though AM induces dependence among causal variants, this dependence is "weak" in the sense that it doesn't change the limiting spectral distribution of the GRM, thereby allowing us to apply arguments in line with those of Jiang and colleagues' (see Supplementary Note 3). Intuitively our result can be summarized as follows: as the sample size and the number of causal variants become large, the eigenvalues of the GRM under AM behave as if the causal variants were independent (as is largely the case under random mating). The behavior of the REML estimator is determined by the behavior of the eigenvalues of the GRM, and thusĥ 2 REML converges to what the heritability would be if the causal variants were independent, i.e., the panmictic heritability.
Simulation studies. We employed a realistic forward-time simulation framework to generate genotypic and phenotypic data. We then used these data to motivate and verify theoretical results. Below, we describe the general framework and specific simulations we performed.
Simulation framework. Given a recombination map and n input individuals' phased biallelic genotypes at p diploid loci as input, we divided the genome into k ( p contiguous, non-overlapping 50kB intervals to obtain a collection of blocks, the disjoint union of which comprises the genome. Recombination events, which occurred with probabilities dictated by the recombination map, were restricted to interval boundaries, thus dramatically reducing the number of haplotypes that to keep track of while maintaining high genomic resolution and realistic LD. To achieve a target population size N sim > n input , N sim pairs of the n input individuals were non-monogamously 'mated' (i.e., matched and subject to meiosis), resulting in a new generation of N sim individuals whose genomes were could be represented in terms of the which of the n input haplotype blocks they inherited at each of the k intervals. We then repeated this random mating procedure for an additional five generations, resulting in N sim chimeric combinations of the original n input genotypes while maintaining the LD structure of the original data. These chimeric genotypes comprised the input for the principal AM simulations.
At the beginning of each particular AM simulation with prespecified panmictic heritability h 2 0 ¼ σ 2 g;0 =ðσ 2 g;0 þ σ 2 e Þ; phenotypic correlation between mates r, p SNPs, and m diploid causal loci z 1 ; ; z m , m ( p, the standardized allele substitution effects u 1 ; ; u m were independently drawn from a Gaussian distribution with expectation zero and variance σ 2 g;0 =m. Unless otherwise stated, all simulations used p = 10 6 SNPs. At each generation, phenotypes were constructed via y ¼ Zu þ e where e was i.i.d. Gaussian with zero expectation and variance σ 2 e . Next, mates were matched according to their respective phenotypes y i , y j such that corr y i ; y j % r 44 . This was achieved by drawing N sim independent doubles were constructed such that i; j À Á k were the positions of w Ã k and w ÃÃ k after concatenating and sorting each element of f w Ã ; w ÃÃ ð Þ T g N sim k¼1 . Similarly, N sim indices l ¼ l 1 ; ; l N sim were constructed such that l k indexed the kth largest of the N sim simulated phenotypes. Finally, each kth mating pair was determined by taking the l bi k =2c th and l bi j =2c th replicates. Having chosen mates, meiosis occurred as previously detailed to construct the next generations' genotypes. All simulations were conducted in Python v3.6.8 45  Simulations using UK Biobank data. For each simulation, the input data were derived from phased, imputed genotypes at p = 10 6 randomly selected imputed SNP loci in a sub-sample of n input = 435,301 European UK Biobank participants 27 . Retained SNPs met the following criteria: minor allele frequency >0.01, Hardy-Weinberg p-value >10 −6 , INFO score >0.95, and presence on the 1000 Genomes Phase 3 (1KG3) reference panel 49 . Genotype data were then phased to the 1KG3 reference panel in batches of 40,000 individuals using Eagle v2.4 50 . Data were then grown to a population of N sim = 10 6 chimeric genotypes and subjected to an additional five generations of random mating as described above before being subject to AM. We conducted AM simulations for varying mating correlations, r 2 f0; :25; :5; :75g and numbers of causal variants, m 2 10 4 ; 10 5 È É , with panmictic heritability fixed at h 2 0 ¼ :5. Each simulation consisted of fifteen generations of AM and produced results congruent with classical theory. Prior to heritability estimation, close relativesπ ≥ :05, were removed using GCTA v1.93.1 20 , resulting in an average sample size of 141,667 across simulated datasets. Additionally, we ran a limited number of larger, more computationally intensive simulations (N sim = 3 × 10 6 ) with mating correlations fixed at r ¼ :5 to investigate the large sample behavior of the REML estimator, resulting in at least 648,000 unrelated individuals across simulated datasets. There were no apparent differences across simulations as a function of m, p, or N sim (Supplementary Fig. 1). We note that these simulations assume a homogenous population without any population structure apart from that induced by primary phenotypic AM; thus, the irrelevance of p here does not necessarily speak to the consequences of increasing the number of variants under consideration by decreasing allele frequency or imputation quality score thresholds, as is common practice, if the additional variants capture additional fine-scale population structure.
Heritability estimation in simulated data. We split each simulated genotype-phenotype dataset into collections of random subsamples mutually exclusive within collection but not across collections, yielding 16 samples of 16,000 individuals, eight samples of 16,000, four samples of 32,000 individuals, two samples of 64,000 individuals, and one sample of 128,000 individuals. We then performed HE regression and single-component REML for each subsample (Fig. 2a). We used GCTA v1.91.3b 20 to construct genomic related matrices and perform HE regression. We obtainedĥ 2 REML using BOLT-LMM v2.3.4 51 for computational efficiency; though BOLT-LMM uses a randomized algorithm, its numerical accuracy is comparable to that of the exact algorithm implemented GCTA 52 .
We also performed a variety of supplementary analyses for a limited set of simulation parameters (r ¼ :5, h 2 0 ¼ :5, and m 2 10 4 ; 10 5 È É ; N ¼ 10 6 ). To demonstrate that including genomic PCs as covariates does not mitigate the impact of AM, we included ten PCs as covariates in the HE regression and REML analyses. For the former, HE regression was conducted in LDAK v5.0 53 , as the HE regression implementation in GCTA cannot accommodate covariates. To demonstrate that the behavior of LDSC under AM is equivalent to that of HE regression (assuming that the LD scores accurately reflect the LD structure of the sample), we used PLINK v1.9 54 to obtain GWAS summary statistics and LDSC v1.0.1 18 to estimate within-sample LD scores using a one centiMorgan sliding window and to perform LDSC regression (Fig. 4a). To demonstrate that multiple variance component (also known as partitioned approaches 30,55 ) do not mitigate the impact of AM, we fit multicomponent HE regression and REML after partitioning SNPs by minor allele frequency and LD score (Fig. 4a). Finally, to assess the scenario wherein a nontrivial fraction of causal variants are missing from Z, we estimated HE regression and REML models after removing 50 or 75% of simulated SNPs at random (Fig. 4b).
Statistics Sampling procedures. We analyzed 1,211,273 biallelic 1KG3 SNPs with in-sample minor allele frequency >0.01, Hardy-Weinberg p-value >10 −6 , and INFO scores >0.95, in a sample of 335,551 unrelated European UK Biobank participants 27 . We selected phenotypes a priori on the basis of previous evidence for AM; we chose height (n ¼ 334,798) and years of education (n ¼ 332,198) as traits with previous evidence of primary phenotypic AM, whereas we chose body mass index (n ¼ 334,429) and bone mineral density (n ¼ 191,330) as negative control traits 14 . We measured years of education following the procedures detailed in ref. 9 . For the results reported in Table 1, we used previously reported estimates of phenotypic spousal correlations 2,4 . For height specifically, we combined the estimates from multiple British cohorts included in a meta-analysis 2 via inverse variance weighting.
Analysis/resampling. We tested for evidence of AM by comparingĥ 2 HE andĥ 2 REML in small and large samples. We randomly selected ten mutually exclusive subsamples of n small = 16,000 individuals for each trait and compared HE and REML estimates in each subsample to the non-overlapping complementary subsample comprised of the remaining n large ¼ n−16,000 individuals, controlling for sex, age, genotyping batch, testing center, and the first ten genomic ancestry PCs. To eliminate variance in heritability estimates due to chance differences in covariate effect estimates across subsamples, we adjusted genotypes and phenotypes in the full sample prior to all following analyses. To our knowledge, existing software is incapable of efficient REML analysis using adjusted genotypes (analogous to dosages) in large samples; e.g., BOLT-REML requires hard-calls as input, whereas GCTA and LDAK have cubic complexity in the number of individuals and markers and would require weeks to run on a high thread-count server. We therefore utilized a modified Python implementation of the REML algorithm presented in ref. 52 (available at https://github.com/rborder/SL_REML 56 ). We used LDAK v5.0 to obtain adjusted HE regression estimates 53 . In order to quantify the divergence ofĥ 2 HE andĥ 2 REML as a function of n, we performed the following test, analogous to a t-test of the interaction effect in a 2 × 2 within-subjects experimental design: where S denotes a given large subsample and S c its complementary small subsample. Though this procedure accounts for the dependence among estimates derived in the same subsamples, the individual observations were derived from various partitionings of the same data and do not constitute truly independent observations. Therefore, the p-values for this test reported in the main text should be interpreted as descriptive despite our application of inferential procedures.
Ethical approval. Ethics approval for the UK Biobank study was obtained by the UK Biobank team from the North West Center for Research Ethics Committee (11/ NW/0382). Access to the UK Biobank data was granted to principal investigator Dr. Matthew C. Keller (researcher ID 16651).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Data are available through the UK Biobank Access Management System (https:// www.ukbiobank.ac.uk/enable-your-research/apply-for-access). A full catalogue of the data examined, including raw materials and descriptive statistics, is available via the online showcase (https://biobank.ndph.ox.ac.uk/ukb). A full description of the study protocols has been provided by the UK Biobank team 60 .