Genome-partitioning of genetic variation for complex traits using common SNPs

Variance explained by all autosomal SNPs for height, BMI, vWF and QTi

We selected 14,347 individuals from three population-based GWAS, i.e. the Health Professionals Follow-up Study (HPFS), the Nurses’ Health Study (NHS) and the Atherosclerosis Risk in Communities (ARIC) study^6–8, and estimated the genetic relationship matrix (GRM) of all the individuals using 565,040 autosomal SNPs which passed quality control (Online Methods). We excluded one of each pair of individuals with an estimated genetic relationship > 0.025 (i.e. more related than third- to fourth-cousins) and retained a subset of 11,586 unrelated individuals. The reason for excluding related pairs is to avoid the possibility that the phenotypic resemblance between close relatives could be due to non-genetic effects (e.g. shared environment) and causal variants not tagged by SNPs but captured by pedigree^10,11. We then fitted the GRM in a mixed linear model (MLM) to estimate the proportion of variance explained by all the autosomal SNPs ( $h_{\mathrm{G}}^2$) for height, BMI, vWF and QTi in each cohort and the combined data where applicable (Online Methods, Table 1 and Supplementary Table 1). The traits vWF and QTi were available from the ARIC sample only. We show that 44.8% (s.e = 2.9%) of phenotypic variance for height can be explained by all the autosomal SNPs, in line with an estimate of 44.5% (s.e. = 8.3%) in a similar analysis of an Australian cohort (3,925 unrelated individuals genotyped by 294,831 SNPs on Illumina arrays, in contrast to the Affymetrix arrays in the present study)⁴. We show for the first time that 16.5% (s.e. = 2.9%), 25.2% (s.e. = 5.1%) and 20.9% (s.e. = 5.0%) of variances for BMI, vWF and QTi can be explained by all the autosomal SNPs, approximately 10-fold, 2-fold and 3-fold larger than the variance explained by all known validated loci found by GWAS for BMI^12–15, vWF¹⁶ and QTi¹⁷, respectively. We note that the ABO blood group locus on chromosome 9 is known to explain approximately 10% of phenotypic variation for vWF¹⁶, through modification of the amount of H antigen expression on the circulating vWF glycoprotein^18,19. The estimate of $h_{\mathrm{G}}^2$ for weight is 18.6% (s.e. = 2.8%). Because of the high phenotypic correlation between BMI and weight (r = 0.92), results for these two traits are very similar. We therefore report results for BMI in the following sections and for completion give all results for weight in the supplementary online material.

Genome-partitioning of genetic variation

Next, we estimated the GRM from the SNPs on each autosome and partitioned the total genetic variance onto individual chromosomes by fitting the GRMs of all the chromosomes simultaneously in a joint analysis (Online Methods). We observed a strong linear relationship between the estimate of variance explained by each chromosome ( $h_{\mathrm{C}}^2$) and chromosome length (L_C, in Mb units) for height (P = 1.4×10⁻⁶ and R² = 0.695) and QTi (P = 1.1×10⁻³ and R² = 0.422) (Fig. 1 and Supplementary Tables 2 and 3). We mapped SNPs to 17,787 genes according to positions on the UCSC Genome Browser hg18 assembly (http://genome.ucsc.edu)²⁰, 17,652 of which had at least one SNP within ± 50 Kb of the 5′ and 3′ untranslated regions (UTRs). There was also a significant correlation between the estimate of $h_{\mathrm{C}}^2$ and the number of genes on each chromosome (N_g(C)) for height and QTi (Supplementary Table 3). Since L_C and N_g(C) are correlated (r = 0.628), we performed a multiple regression analysis of the estimate of $h_{\mathrm{C}}^2$ on L_C and N_g(C), and fitted models in which chromosome length was fitted after the number of genes and vice versa. When including both L_C and N_g(C) in the regression model, N_g(C) was not significant and L_C was still significant for height and QTi (Supplementary Table 3). The regression of the estimate of $h_{\mathrm{C}}^2$ on either L_C or N_g(C) was not significant for BMI and vWF. These results are consistent with variance explained by each chromosome for height and QTi (but less so for BMI and vWF) being proportional to the proportion of the genome being considered. Although longer chromosomes harbour more genes that are implicated in abnormal growth or skeletal development, the relationship between variance explained for height and chromosome length remains significant after fitting the number of such genes (Supplementary Fig. 1). We provide evidence that the linear relationship between the estimate of $h_{\mathrm{C}}^2$ and L_C cannot be attributed to the fact that longer chromosomes have more SNPs and thereby smaller sampling errors when estimating genetic relationships between individuals (Supplementary Note and Supplementary Figs. 2 and 3).

However, genes vary greatly in size and when we considered the length of genes we observed that the estimate of $h_{\mathrm{C}}^2$ for height and QTi was also proportional to the total length of genes on each chromosome (L_g(C)) where gene length is defined as physical distance between the beginning and end of UTRs (Supplementary Fig. 4). Since the correlation between L_C and L_g(C) is extremely high (r = 0.97), we were unable to discriminate whether L_C or L_g(C) is causative by multiple regression, i.e. the regression of $h_{\mathrm{C}}^2$ on L_C is not significant fitted after L_g(C) and vice versa (Supplementary Table 3). Therefore, a different analysis was required. We asked whether we could still observe a significant regression of $h_{\mathrm{C}}^2$ on L_g(C) when chromosome length was held constant. We implemented this by dividing the genome into segments with the same length of 50 or 30 Mb and then estimated the variance explained by each segment ( $h_{\mathrm{S}}^2$) in a joint analysis (Online Methods). We found that the regression of $h_{\mathrm{S}}^2$ on the total gene length per segment (L_g(S)) remained significant for height, P = 1.7×10⁻³ for 50 Mb segments and P = 1.2×10⁻⁴ for 30 Mb segments (Supplementary Fig. 5). The regressions of $h_{\mathrm{S}}^2$ on the number of genes, the total length of exons and the number of exons on each segment were also significant in some cases, but none of them was significant when fitted after L_g(S) while L_g(S) was always significant fitted after any of them (Supplementary Table 4). These results suggest that, at least for height, genomic regions explain variation in proportion to their genic content.

To quantify these effects genome-wide, we partitioned the variance explained by all the SNPs onto genic ( $h_{\mathrm{Gg}}^2$) and intergenic ( $h_{\mathrm{Gi}}^2$) regions of the whole genome (Online Methods). We defined the gene boundaries as ± 0 Kb, 20 Kb and 50 Kb of the 3′ and 5′ UTRs. A total of 213,509, 282,058, and 336,127 SNPs were located within the boundaries of 13,406, 17,277 and 17,652 protein-coding genes for the three definitions (± 0 Kb, 20 Kb and 50 Kb) respectively, which covered 35.8%, 49.4% and 58.7% of the genome. Some genes did not have any SNPs within them, especially if the most stringent definition of gene boundary (± 0 Kb) was used. We tested the estimates of $h_{\mathrm{Gg}}^2$ and $h_{\mathrm{Gi}}^2$ against the expected values from the genic and intergenic coverages by a goodness-of-fit test. We found strong evidence for height and vWF, and less so for BMI and QTi, that genic regions proportionally explain more variation than intergenic regions (legends of Fig. 2 and Supplementary Fig. 6). As an example we consider the case of genes ± 20 Kb where genic and intergenic coverages are roughly equal (49.4% vs. 50.6%). The estimates of $h_{\mathrm{Gg}}^2$ vs. $h_{\mathrm{Gi}}^2$ are 32.8% vs. 12.6% (P = 2.1×10⁻¹⁰) for height, 22.7% vs. 4.0% (P = 5.1×10⁻⁴) for vWF, 11.7% vs. 4.7% (P = 0.022) for BMI and 13.5% vs. 7.5% (P = 0.251) for QTi. We further partitioned the genetic variance onto the genic and intergenic regions of each chromosome (Online Methods). In general, the results agree with that of the whole-genome partitioning analysis, in that the genic regions proportionally explained more variation (Fig. 2 and Supplementary Fig. 6). The variance attributable to chromosome 9 for vWF is dominated by the genic regions, which is expected because the ABO gene on this chromosome explains ~10% of its variance¹⁶. However, there appear to be exceptions, for example, the intergenic regions of chromosome 2 and chromosome 5 seemed to be more important for BMI and QTi, respectively. These results are not conclusive because the standard errors of the estimates are large. Despite these special cases, overall the results are consistent with causal variants being more likely to be located in the vicinity of functional genes.

Variance attributed to population structure

Mixed linear model methods are useful to control for population structure in GWAS^21,22. Population structure in the data causes correlations of SNPs on different chromosomes. Consequently, fitting only one chromosome in the model (separate analysis) also captures some of the variance due to other chromosomes, so that the estimate of variance explained by each chromosome from the separate analysis ( $h_{\mathrm{C}}^2(\mathrm{sep})$) is biased upwards (Online Methods). The joint analysis has the advantage of protecting against such inter-chromosomal correlations because the estimate of each $h_{\mathrm{C}}^2$ is conditional on the other chromosomes in the model so that the estimates of variance explained by different chromosomes are independent of each other. We therefore can calculate the variance attributable to population structure by comparing the estimates between $h_{\mathrm{C}}^2(\mathrm{sep})$ and $h_{\mathrm{C}}^2$. The inter-chromosomal SNP correlations occur for two reasons: 1) cryptic relatedness (e.g. unexpected cousins) because closely related individuals will share SNPs identical by descent on more than one chromosome; 2) systematic difference in allele frequencies between subpopulations (population stratification). We modelled the variance attributed to these two forms of population structure as $h_{\mathrm{C}}^2(\mathrm{sep})-h_{\mathrm{C}}^2=b_0+b_1{L}_{\mathrm{C}}+\epsilon$, where the slope b₁ allows the possibility that longer chromosomes track population structure better than smaller chromosomes. To illustrate the effect of population structure, we also estimated $h_{\mathrm{C}}^2(\mathrm{sep})$ and $h_{\mathrm{C}}^2$ in the entire sample of 14,347 individuals (i.e. without removing close relatives). The intercept b₀ appears to be due to cryptic relatedness because when we eliminate relatives with a relationship > 0.025, b₀ declines to zero (Fig. 3). We therefore predicted that cryptic relatedness accounted for 1.5%, 0.084%, 0.22% and 0.065% (not significant) of the phenotypic variance for height, BMI, vWF and QTi, respectively, in the entire sample. The variance attributed to cryptic relatedness is irrespective of chromosome length because it does not require many SNPs per chromosome to detect close relatives. Conversely, the regression slope b₁ appears to be due to population stratification because longer chromosomes are likely to have more ancestry informative markers (AIMs), assuming that the AIMs are randomly distributed across the genome. We then predicted that population stratification accounted for 6.9×10⁻⁵L_C, and 7.2×10⁻⁶L_C, −1.92×10⁻⁶L_C (not significantly different from zero) and 2.3×10⁻⁵L_C of variance for height, BMI, vWF and QTi, respectively, in the entire sample and a similar amount in the data set of unrelated individuals (Fig. 3). The difference between $h_{\mathrm{C}}^2(\mathrm{sep})$ and $h_{\mathrm{C}}^2$ represents the overall effect of all the other 21 chromosomes on one chromosome. Therefore, the proportion of variance attributed to population structure (cryptic relatedness and population stratification) across the whole genome is approximately equal to $b_{0}22/21+b_1{\displaystyle {\sum }_{\mathrm{C}=1}^{22}}{L}_{\mathrm{C}}/21$, which is (1.6% + 0.91%), (0.088% + 0.095%), (0.23% + 0.0%) and (0.068% + 0.30%) for height, BMI, vWF and QTi, respectively, in the entire sample. Hence, we provide a simple approach to estimate and partition the variance attributed to population structure for complex traits. The variances due to cryptic relatedness and population stratification depend on the data structure in the sample. Therefore, the estimates we present above are specific for the data in this study.

It is common to fit eigenvectors (PCs) from principal component analysis (PCA) in single SNP association studies to correct for possible population structure^23,24. We show that fitting the first 10 PCs and one chromosome at a time or fitting all chromosomes simultaneously without fitting PCs, led to similar estimates of the variance explained by each chromosome (Supplementary Fig. 7), which suggests that the majority of variance attributed to population structure is well captured by the first 10 PCs in these data.

Estimation of variance explained by the X-chromosome

We estimated the GRM for the X-chromosome and parameterized it under three assumptions of dosage compensation¹⁰: 1) equal X-linked variance for males and females; 2) no dosage compensation i.e. both X-chromosomes are active for females; 3) full dosage compensation i.e. one of the X-chromosomes is completely inactive for females. We fitted the parameterized GRMs for the X-chromosome in an MLM whilst simultaneously estimating $h_{\mathrm{G}}^2$ in the model to capture the genetic variation on the autosomes and variation due to possible population structure. For all the traits, the full-dosage compensation model fits the data best and the no-dosage compensation model is the worst with the equal-variance model in between (Supplementary Table 5). However, the differences in estimates are relatively small and none of them is statistically significant. Larger datasets will be required to distinguish such small differences. Under the assumption of full dosage compensation, the variance attributable to the X-chromosome for females was 0.61% (s.e. = 0.32%), 0.82% (s.e. = 0.35%), 0.57% (s.e. = 0.52%) and 0.0% (s.e. = 0.48%) for height, BMI, vWF and QTi, respectively. To verify that we detect heterogeneous variances on the X-chromosome rather than autosomal variance differences between males and females, we fitted the same dosage compensation models for the autosomes. As expected, the equal variance model fitted the data best and the full dosage compensation model was the worst fit for all the traits (Supplementary Table 6). Therefore, the data are consistent with twice as much additive genetic variation for height, BMI, and vWF on the X-chromosome in males as in females, which is predicted from theory under the assumption of random X-inactivation²⁵. While there are syndromic examples illustrating the phenotypic effect of the Lyon hypothesis (e.g., Turner’s syndrome and Kleinfelter syndrome), to our knowledge, this is the first empirical evidence from genotype-phenotype associations on complex traits that the amount of genetic variation on the X-chromosome appears consistent with X-chromosome inactivation. However, the evidence is indirect and not overwhelming. Larger samples sizes and the detection of multiple associated loci on the X-chromosome will be necessary to investigate the expression of genes on the X-chromosome that affect the traits studied.

Comparison with variants known to be associated with BMI, vWF and height

To quantify the effect of known associated variants on the results, we included the FTO SNP rs9939609 on chromosome 16 for BMI and the ABO SNP rs612169 on chromosome 9 for vWF as a covariate when estimating $h_{\mathrm{C}}^2$ by the joint analysis of all autosomes. FTO was the first locus to be detected through GWAS that is associated with BMI¹⁴ and ABO is a major determinant of vWF¹⁹. When compared to the result without adjustment, the estimate of variance due to chromosome 16 ( $h_{16}^2$) for BMI decreased from 1.19% to 0.61%, in line with an estimate of ~0.34% to ~1% of variance explained by the FTO locus for BMI in previous GWAS^12,14,15 and an estimate of ~0.46% from association analysis in the present study; the estimate of $h_9^2$ for vWF decreased by 11.8%, consistent with an estimate of ~10% of variance for vWF explained by the ABO locus in GWAS¹⁶; and the estimates for the other chromosomes remained the same (Supplementary Fig. 8).

The meta-analysis of ~133,000 individuals by the GIANT consortium has identified 180 independent loci associated with genetic variation of height²⁶. The estimate of $h_{\mathrm{C}}^2$ by joint analysis in our study shows a high correlation (r = 0.715 and P = 1.8×10⁻⁴) with the sum of the variance explained at the associated loci on each chromosome from the GIANT meta-analysis (Fig. 4).

Additional models

We fitted a number of other models to quantify the effect of having multiple phenotypic observations per individual, to test for genotype-sex interaction effects and for the effect of sample ascertainment. We also estimated the genetic correlation between height and weight.

As both height and BMI have repeated measures in the ARIC cohort and BMI also has repeated measures in the HPFS and NHS cohorts, we fitted a repeatability model for these repeated records, assuming that the genetic correlation between repeated observations is unity. We estimated the repeatability of height as ~0.99 in the ARIC study and the repeatabilities of weight and BMI of > 0.93 in all cohorts (Supplementary Table 7). The estimates of $h_{\mathrm{G}}^2$ by repeatability model analyses are similar to those using the mean of the repeated measures as in all other analyses (Table 1 and Supplementary Table 7). The ratio of an estimate of $h_{\mathrm{G}}^2$ based upon m repeated records to that of a single observation is 1/[(1 − ρ)/m + ρ] with ρ the repeatability²⁷. In this study, ρ > 0.93 and m is small, so that our inference based upon the mean of 3 or 4 observations is very similar to that from a single observation.

We fitted a genotype-sex interaction effect in the model, and did not find any significant interaction effect for all the traits (Supplementary Table 8). We also estimated the variance explained by all the autosomal SNPs in each gender group in each cohort separately, and did not observe clear evidence of genetic heterogeneity between males and females (Supplementary Table 9). However, there are exceptions, for example, the p-value for a test of heterogeneity between males and females for vWF was 0.007 (0.024 after PC-adjustment), which may indicate real genetic difference between males and females for this trait. When taking multiple testing into account, the test was not significant after a Bonferroni correction. For weight and BMI, the p-values only verge on significance even when $h_{\mathrm{G}}^2$ are estimated to be zero in males, suggesting a lack of power to detect an interaction.

The HPFS and NHS cohorts were ascertained by type-II diabetes (T2D), which might affect the estimate of genetic variance because BMI is a risk factor for T2D²⁸. We performed analyses of the variance explained by all the autosomal SNPs separately in T2D cases and controls of the HPFS and NHS cohorts, and did not observe significant differences in estimated genetic variance between cases and controls for all the three traits (Supplementary Table 10). Moreover, we estimated $h_{\mathrm{G}}^2$ in the combined dataset excluding the T2D cases, and the estimate was not different from that using all samples (Supplementary Table 11). These two additional analyses show little, if any, impact of the ascertainment of T2D cases on the estimate of variance explained by the SNPs.

All analyses we have performed are univariate, i.e. a single phenotype at a time. However, multivariate models fit easily in the same analysis framework and the only limitation is computational. As an example, we approximated a full bivariate analysis of height and weight by using logarithms of the phenotypes, exploiting the relationship log(BMI) = log(Weight) − 2log(Height), so that from three univariate analyses we can estimate the genetic correlation between log(Weight) and log(Height). We estimated a genetic correlation of 0.45 (s.e. = 0.17) between log(Weight) and log(Height) (Supplementary Table 12). Although this is on the logarithm scale, the estimate on the observed scale is unlikely to be very different. This estimate indicates that for genetic variation tagged by common SNPs there is a substantial overlap in genome-wide additive factors for height and weight.

GWAS samples and quality control

Details of the HPFS, NHS and ARIC cohorts have been described previously^6–8. The GWAS data in terms of study design, sample selection and genotyping have been detailed by Qi et al.³⁷ for the HPFS and NHS cohorts and by Psaty et al.⁸ for the ARIC cohort. All three cohorts have been studied as part of the GENEVA (The Gene, Environment Association Studies) project³⁸ and this study has benefitted by using data from this consortium that have been generated and cleaned using a common protocol. We selected 6,293 individuals (2,745 T2D cases and 3,148 controls) from the NHS and HPFS cohorts and 15,792 individuals from the ARIC cohort. All of these selected individuals were genotyped using the Affymetrix Genome-Wide Human 6.0 array.

Of the 909,622 SNP probes, 874,517 (HPFS), 879,071 (NHS) and 841,820 (ARIC) passed quality control analysis by the Broad Institute and the GENEVA Coordinating Center (excluding SNPs with missing call rate ≥ 5% or plate association P < 1×10⁻¹⁰)⁹. We further excluded SNPs with missing rate ≥ 2%, > 1 discordance in the duplicated samples, Hardy-Weinberg equilibrium P < 1×10⁻³ or minor allele frequency < 0.01. A total of 687,398 (27,578), 665,163 (24,108) and 593,521 (23,664) autosomal (X-chromosome) SNPs were retained for the HPFS, NHS and ARIC cohorts, respectively, 565,040 (21,858) of which were in common across the three cohorts.

We included only one of each set of duplicated samples and one of each pair of samples which were identified as full siblings by an initial scan of relatedness in PLINK³⁹. We investigated population structure by PCA of all the autosomal SNPs that passed QC, and included only samples of European ancestry (Supplementary Fig. 9). We excluded samples with gender misidentification by examining the mean of the intensities of SNP probes on the X and Y chromosomes. We also excluded samples with missing call rate ≥ 2% and samples on two plates which showed extremely high level of mean inbreeding coefficients. A total of 2,400 (HPFS), 3,265 (NHS) and 8,682 (ARIC) samples were retained for analysis respectively, with a combined set of 14,347 samples.

Phenotypes

Summary statistics of the phenotypes of height, weight, BMI, vWF and QTi are shown in Supplementary Table 13. There are three measures of weight and a single measure of height in both HPFS and NHS cohorts, four measures of weight and three measures of height in the ARIC cohort, and single measures of vWF and QTi in the ARIC cohort. For height, weight and BMI, we used the mean of repeated measures in all the analyses except for the analysis of the repeatability model. We adjusted the phenotypes (or the mean phenotype) for age and standardized it to a z-score in each gender group in each of the three cohorts separately.

Statistical analysis

We estimated the GRM of all individuals in the combined data from all the autosomal SNPs using the method we recently developed^4,10, and excluded one of each pair of individuals with an estimated genetic relationship > 0.025. We then estimated the variance explained by all autosomal SNPs by restricted maximum likelihood analysis of an MLM y = Xβ + g_G + ε, where y is a vector of phenotypes, β is a vector of fixed effects (e.g. the first 10 PCs) with its incidence matrix of X, g_G is a vector of aggregate effects of all autosomal SNPs with $\mathrm{var}({\mathbf{g}}_{\mathrm{G}})={\mathbf{A}}_{\mathrm{G}}{\sigma }_{\mathrm{G}}^2$ and A_G is the GRM estimated from all autosomal SNPs. The proportion of variance explained by all autosomal SNPs is defined as $h_{\mathrm{G}}^2={\sigma }_{\mathrm{G}}^2/{\sigma }_{\mathrm{P}}^2$ with ${\sigma }_{\mathrm{P}}^2$ being the phenotypic variance.

Furthermore, we estimated the GRM from the SNPs on each chromosome (A_C) and estimated the variance attributable to each chromosome by fitting the GRMs of all the chromosomes simultaneously in the model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+{\displaystyle \sum _{\mathrm{C}=1}^{22}}{\mathbf{g}}_{\mathrm{C}}+\mathbf{\epsilon }$ where ${\mathbf{g}}_{\mathrm{C}}$ is a vector of genetic effects attributable to each chromosome and $\mathrm{var}({\mathbf{g}}_{\mathrm{C}})={\mathbf{A}}_{\mathrm{C}}{\sigma }_{\mathrm{C}}^2$ (joint analysis). The proportion of variance explained by each chromosome is defined as $h_{\mathrm{C}}^2={\sigma }_{\mathrm{C}}^2/{\sigma }_{\mathrm{P}}^2$. We also fitted one chromosome at a time in the model y = Xβ + g_C + ε (separate analysis). If there is an effect of population structure, SNPs on one chromosome will be correlated with the SNPs on the other chromosomes such that $h_{\mathrm{C}}^2(\mathrm{sep})$ will be overestimated in the separate analysis.

We extended the joint analysis of chromosomes to that of genomic segments. We divided the genome evenly into N_S segments with each of d_S Mb length, and then estimated the GRM using the SNPs on each segment. We estimated the variance explained by each segment ( $h_{\mathrm{C}}^2$) by fitting the GRMs of all the segments in an MLM $\mathbf{y}=\mathbf{X}\mathbf{\beta }+{\displaystyle \sum _{\mathrm{S}=1}^{N_{\mathrm{S}}}{\mathbf{g}}_{\mathrm{S}}}+\mathbf{\epsilon }$ where g_s is a vector of genetic effects attributable to each segment.

We further partitioned the variance explained by all the SNPs onto genic and intergenic regions of the whole genome ( $h_{\mathrm{Gg}}^2$ and $h_{\mathrm{Gi}}^2$) as well as that of each chromosome ( $h_{\mathrm{Cg}}^2$ and $h_{\mathrm{Ci}}^2$). The gene boundaries were defined as ± d_g Kb away from the 3′ and 5′ UTRs. We estimated $h_{\mathrm{Gg}}^2$ and $h_{\mathrm{Gi}}^2$ by fitting all the genic and intergenic SNPs in an MLM y = Xβ + g_Gg + g_Gi + ε, and estimated $h_{\mathrm{Cg}}^2$ and $h_{\mathrm{Ci}}^2$ by fitting the genic and nongenic SNPs on individual chromosomes in the model $\mathbf{y}=\mathbf{X}\mathbf{\beta }+{\displaystyle \sum _{\mathrm{C}=1}^{22}{\mathbf{g}}_{\mathrm{Cg}}}+{\displaystyle \sum _{\mathrm{C}=1}^{22}{\mathbf{g}}_{\mathrm{Ci}}}+\mathbf{\epsilon }$.

We estimated the variance attributable to the X chromosome using the method we recently developed¹⁰. In brief, we estimated the GRM for the X-chromosome (A_X) using the following equations

• ${\widehat{A}}_{\mathit{jk}}^{\mathrm{M}}={\displaystyle \sum _i^N\frac{(x_{\mathit{ij}}^{\mathrm{M}}-p_i)(x_{\mathit{ik}}^{\mathrm{M}}-p_i)}{p_i(1-p_i)}}$ for a male-male pair

• ${\widehat{A}}_{\mathit{jk}}^{\mathrm{F}}={\displaystyle \sum _i^N\frac{(x_{\mathit{ij}}^{\mathrm{F}}-2p_i)(x_{\mathit{ik}}^{\mathrm{F}}-2p_i)}{2p_i(1-p_i)}}$ for a female-female pair, and

• ${\widehat{A}}_{\mathit{jk}}^{\mathrm{MF}}={\displaystyle \sum _i^N\frac{(x_{\mathit{ij}}^{\mathrm{M}}-p_i)(x_{\mathit{ik}}^{\mathrm{F}}-2p_i)}{\sqrt{2}{p}_i(1-p_i)}}$ for a male-female pair,

where $x_{\mathit{ij}}^{\mathrm{M}}$ and $x_{\mathit{ij}}^{\mathrm{F}}$ are the number of copies of the reference allele for an X-chromosome SNP for a male and a female, respectively, p_i is the frequency of the reference allele and N is the number of SNPs. Assuming the male-female genetic correlation to be 1, the X-linked phenotypic covariance is ${\mathrm{cov}}_{\mathrm{X}}(y_j^{\mathrm{M}},y_k^{\mathrm{M}})=A_{\mathit{jk}}^{\mathrm{M}}{\sigma }_{\mathrm{X}(\mathrm{M})}^2$ for a male-male pair, ${\mathrm{cov}}_{\mathrm{X}}(y_j^{\mathrm{F}},y_k^{\mathrm{F}})=A_{\mathit{jk}}^{\mathrm{F}}{\sigma }_{\mathrm{X}(\mathrm{F})}^2$ for a female-female pair, or ${\mathrm{cov}}_{\mathrm{X}}(y_j^{\mathrm{M}},y_k^{\mathrm{F}})=A_{\mathit{jk}}^{\mathrm{MF}}{\sigma }_{\mathrm{X}(\mathrm{M})}{\sigma }_{\mathrm{X}(\mathrm{F})}$ for a male-female pair^25,40, where ${\sigma }_{\mathrm{X}(\mathrm{M})}^2$ and ${\sigma }_{\mathrm{X}(\mathrm{F})}^2$ are X-linked genetic variances for males and females, respectively. Assumptions about inactivity of the X-chromosome (dosage compensation) imposed a relationship between ${\sigma }_{\mathrm{X}(\mathrm{M})}^2$ and ${\sigma }_{\mathrm{X}(\mathrm{F})}^2$, which allow a single variance component ${\sigma }_{\mathrm{X}(\mathrm{F})}^2$ to account for the X-linked genetic variance for both sexes. Therefore, we can express the X-linked phenotypic covariances as ${\mathrm{cov}}_{\mathrm{X}}(y_j^{\mathrm{M}},y_k^{\mathrm{M}})=d^2{A}_{\mathit{jk}}^{\mathrm{M}}{\sigma }_{\mathrm{X}(\mathrm{F})}^2$, ${\mathrm{cov}}_{\mathrm{X}}(y_j^{\mathrm{F}},y_k^{\mathrm{F}})=A_{\mathit{jk}}^{\mathrm{F}}{\sigma }_{\mathrm{X}(\mathrm{F})}^2$, and ${\mathrm{cov}}_{\mathrm{X}}(y_j^{\mathrm{M}},y_k^{\mathrm{F}})=d{A}_{\mathit{jk}}^{\mathrm{MF}}{\sigma }_{\mathrm{X}(\mathrm{F})}^2$, where d is the lyonization coefficient, ${\sigma }_{\mathrm{X}(\mathrm{M})}=d{\sigma }_{\mathrm{X}(\mathrm{F})}$, which takes 1 under the hypothesis of equal X-linked genetic variance for both sexes, takes $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{2}$}\right.$ under the hypothesis of no dosage compensation (both X-chromosomes are active for females) and takes $\sqrt{2}$ under the hypothesis of full dosage compensation (complete inactivity of one X-chromosome for females) (Supplementary Note). In the analysis of MLM, we took the lyonization coefficient into the account by parameterizing the raw A_x matrix, i.e. ${\mathbf{A}}_{\mathrm{X}}^{\mathrm{P}}=d^2{\mathbf{A}}_{\mathrm{X}}$ for male pairs, ${\mathbf{A}}_{\mathrm{X}}^{\mathrm{P}}={\mathbf{A}}_{\mathrm{X}}$ for female pairs and ${\mathbf{A}}_{\mathrm{X}}^{\mathrm{P}}=d{\mathbf{A}}_{\mathrm{X}}$ for male-female pairs. We estimated ${\sigma }_{X(F)}^2$ under the three hypotheses by fitting the parameterized GRM for the X-chromosome ( ${\mathbf{A}}_{\mathrm{X}}^{\mathrm{P}}$) conditional on the GRM estimated from all autosomal SNPs in an MLM $\mathbf{y}=\mathbf{X}\mathbf{\beta }+{\mathbf{g}}_{\mathrm{X}}+{\mathbf{g}}_{\mathrm{G}}+\mathbf{\epsilon }$, where g_X is a vector of X-linked genetic effects with $\mathrm{var}({\mathbf{g}}_{\mathrm{X}})={\mathbf{A}}_{\mathrm{X}}^{\mathrm{P}}{\sigma }_{\mathrm{X}(\mathrm{F})}^2$.