Genome-wide genetic homogeneity between sexes and populations for human height and body mass index

Abstract

Introduction

There has been a long-standing question about the contribution of autosomal sex-specific genetic effects to complex trait variation in humans (1,2). For many traits, including height, body mass index (BMI), waist circumference, hip circumference and waist-to-hip ratio, twin studies find systematically lower correlations of opposite-sex twin pairs compared with same-sex dizygotic twin pairs (3). For example, Schousboe et al. (4) concluded from a large study of European twin pairs aged 20–40 years that ‘the sets of genes contributing to variation in BMI are not identical among men and women’, with an estimate of genetic correlation between men and women <0.5 (averaged across cohorts and age groups). However, genome-wide association studies (GWAS) using very large sample sizes have not identified any autosomal genetic variant that showed a significant sex-specific effect on height or BMI (5). In addition, pedigree studies in other species consistently report very large genetic correlations in body size between males and females (6–10), suggesting that little phenotypic variation is contributed by genotype–sex interactions. In this study, we address the question about genotype-by-sex and genotype-by-country interaction variation using a powerful population design in which unrelated individuals have genome-wide single nucleotide polymorphism (SNP) data and a phenotypic measure on height and BMI. We estimate the correlation in genome-wide additive effects that are captured by SNP data with high precision.

Results

Using the bivariate genomic restricted maximum likelihood (GREML) method (11) as implemented in the genetic complex trait analysis (GCTA) software (12), we estimated the genetic correlation between men and women for height and BMI in a combined GWAS data set comprising 44 126 unrelated individuals of European descent and ∼1.2 M SNPs (see Material and Methods section). Since the bivariate GCTA-GREML method does not require traits measured on the same set of individuals, we regarded height (or BMI) in men as a distinct trait from that in women and used genome-wide SNP data to estimate genetic correlation between the ‘two traits’ (e.g. men's height versus women's height). The expected value of the genetic correlation (r_g) should, per definition, range from −1 or 1. In the analysis, we allowed the estimate to go beyond this range so that the estimate was unbiased (an estimate constrained between −1 and 1 is a biased estimated). The estimate of r_g was 1.02 for height with a standard error (SE) of 0.023, and 1.01 (SE = 0.064) for BMI (Table 1), both of which were not significantly different from 1 (Wald's test P = 0.48 for height and P = 0.86 for BMI).

A recently proposed alternative method, called bivariate linkage disequilibrium score (LDSC) regression analysis (13), is also able to estimate genetic correlations between traits measured on different samples. This method is built on the univariate LDSC regression analysis (14), which was developed to distinguish the effects due to population stratification from polygenic effects in GWAS. The LDSC regression method only requires GWAS summary data so that it can be applied to data with very large sample sizes. We used summary data from the sex-stratified Genetic Investigation of ANthropometric Traits (GIANT) meta-analyses (5) (up to 60 586 men and 73 137 women) for height and BMI (see Material and Methods section). The estimate of r_g from the bivariate LDSC regression analysis using summary data was 0.957 (SE = 0.023) for height and 0.879 (SE = 0.035) for BMI (Table 2). The estimate for BMI was significantly different from 1 (Wald's test P = 5.9e-4), which was inconsistent with the result from bivariate GREML analysis on individual-level data. We therefore investigated what could cause this apparent discrepancy.

Unlike the GREML analysis where there were roughly the same number of men and women in each cohort (Table 3), the bivariate LDSC regression analysis used summary data from a meta-analysis where there were cohorts that only contributed samples of one sex, e.g. all > 20 000 samples from the Women's Genome Health Study (WGHS) were women (see the Supplementary Material, Table S1 of Randall et al. (5)). Therefore, if there is genetic heterogeneity between cohorts (populations), the estimate of r_g between sexes from the bivariate LDSC regression analysis will be biased downwards, and this may explain the discrepancy between the estimates of r_g from bivariate GREML and LDSC regression. We therefore used the bivariate GREML analysis to test whether there is genetic heterogeneity between populations in the combined GWAS sample with individual-level data (n = 44 126). We stratified the data (from both sexes) into two groups, i.e. samples from the USA (ARIC, GENEVA-T2D, HRS) and Europe (TwinGene, Lifelines and EGCUT), and then used bivariate GCTA-GREML approach to estimate the genetic correlation between USA and Europe for height and BMI. The estimates of between-population r_g were 0.961 (SE = 0.030) for height and 0.894 (SE = 0.058) for BMI (Table 4). These estimates are remarkably consistent with the between-sex r_g estimates of 0.957 for height and 0.879 for BMI (Table 2) from bivariate LDSC regression using summary data. The results are consistent with the hypothesis that there is genetic heterogeneity between populations for BMI that confounds the estimate of r_g between sexes in the bivariate LDSC regression analysis. The European samples (TwinGene, Lifelines and EGCUT) used in our bivariate GREML analyses are from Sweden, the Netherland and Estonia, which are not a perfect representative of the whole Europe, whereas the US samples (ARIC, GENEVA-T2D, HRS) are of descent of almost all Europeans (see the PCA plot in the Supplementary Material, Fig. S19 of Yang et al. (15)). Therefore, the between-population genetic heterogeneity for BMI could be due to genetic difference (e.g. difference in causal variants or difference in linkage disequilibrium between SNPs and causal variants) across European populations. Of course, we could not rule out the possibility that there is a genotype-by-environment interaction for BMI, where population (Europe or USA) is a proxy for environmental difference.

We noticed that the estimates of proportion of variance explained by all SNPs (h²_g) from the LDSC regression analysis (Table 2) were much smaller than those from bivariate GREML analysis for both height and BMI (Table 1). The discrepancy could be partly attributed to the error of approximating the LDSC in the sample used for analysis by the estimate from the 1000 Genome Project (1KGP) data because the regression slope (from which h²_g is estimated in the LDSC regression analysis) is inversely proportional to the variance of the estimated LDSCs, which consists of the true variation of LDSCs and the variation of estimation errors. This was demonstrated by the analysis of applying the LDSC regression analysis in the same data that were used in the GREML analysis (see Material and Methods section), where the estimates of h²_g from LDSC regression (Supplementary Material, Table S1) were consistently smaller than those from GREML (Table 1). However, the estimates of r_g are generally consistent between the two methods with the standard errors from LDSC regression being about twice larger than those from GREML (Table 1 and Supplementary Material, Table S1). Another explanation of the discrepancy between estimate of h²_g from GREML using individual-level data (Table 1) and that from LDSC regression using summary data (Table 2) is that the summary statistics available in the public domain had been corrected by the genomic control (GC) approach (16), i.e. reported SE equals to original SE multiplied by $λGC$ in each cohort and further in the whole meta-analysis sample, where $λGC$ is the genomic inflation factor (median of χ² statistics divided by 0.455). The LDSC regression analysis was based on z-statistics that are shrunk by GC correction, resulting in an underestimated h²_g. This is demonstrated by the decrease in the estimate of h²_g from LDSC regression as a result of GC correction in the combined GWAS sample with individual-level data (Supplementary Material, Table S1). Fortunately, GC correction did not affect the estimate of r_g because the genetic variance and the genetic covariance were scaled by the same factor, which cancels out when calculating r_g (Supplementary Material, Table S1). We cannot rule out the possibility that the deflated estimates of h²_g from LDSC regression using meta-analysis summary data were partly due to technical factors [e.g. difference in design, genotyping platform and/or quality control (QC) criteria] that led to artifactual discrepancies in SNP calls. In addition, summary data from meta-analysis based on a large number of cohorts may suffer from more experimental noise, such that the implicit h²_g is lower than that from an equivalent sample with individual-level data. In addition, we used an IMPUTE-INFO threshold of 0.3 for SNP inclusion in the analysis of the combined GWAS data. We demonstrated by additional analysis that the GREML estimates using an IMPUTE-INFO threshold of 0.6 (Supplementary Material, Table S2) were highly consistent with those using 0.3 (Table 1).

We further estimated the between-sex genetic correlation for BMI-adjusted waist circumference (WCadjBMI), hip circumference (HIPadjBMI) and waist–hip ratio (WHRadjBMI) using individual-level data in the combined GWAS sample (bivariate GREML analysis, n = up to 29 091) and using summary-level data from the sex-specific GIANT meta-analysis (bivariate LDSC regression analysis, n = up to 82 057) (see Material and Methods section). The waist circumference, hip circumference and waist–hip ratio (WHR) phenotypes were adjusted for BMI to capture body fat distribution independent of overall adiposity (5). Although the standard errors of the estimates were relatively large, the results from both analyses suggested that there was a significant between-sex genetic heterogeneity for WHRadjBMI, r_g = 0.607 (SE = 0.112) from GREML analysis and r_g = 0.770 (SE = 0.108) from LDSC regression, consistent with WHR being a sexually dimorphic trait under differential selection for men versus women (17,18). Our results also seem to suggest that there is between-sex genetic heterogeneity for WCadjBMI but not for HIPadjBMI (Tables 1 and 2). These estimates are consistent with the GWAS results reported in the Randall et al. study that out of 9 traits analysed, only WHRadjBMI and WCadjBMI had significant sex-specific effects (detected at 7 loci) (5).

Discussion

The search for genetic variants with sex-specific effects and thereby sex-genotype interaction variation is often motivated by large differences in average phenotypes (as in height) or biological processes that are known to differ between the sexes (e.g. hormonal effects). However, mean differences do not imply that sex-by-genotype interactions exist. For example, a recent study (19) reported a mean difference in genetic profile for height across European nations, consistent with the observed difference in mean height phenotype. This finding does not conflict with the result from a previous study (20) that there is no evidence for a difference in heritability between European populations or result from our study (Table 4) that the between-population genetic correlation for height is not significantly different from unity. One caveat of our study is that we estimate the genetic correlation from effects captured by common SNPs. It is possible that rare variants of large effects individually have a different effect in men and women. Nevertheless, it is difficult to see how the missing part of genetic (co)variance would cumulatively lead to a much lower overall genome-wide correlation, since common SNPs imputed to the 1KGP reference capture the majority of genetic variation for height and BMI (15).

In conclusion, we provide precise estimates of the genetic correlation between men and women for height and BMI using very large samples of unrelated individuals. The results are consistent with the hypothesis that there is no between-sex genetic heterogeneity for both height and BMI and that there is some between-population genetic heterogeneity for BMI, suggesting that obesity in men and women is influenced by the same set of autosomal genes in the same population, and the effect of these genes may vary across populations. The analyses performed in this study in principle can be applied to detect genome-wide heterogeneity between sexes, populations or environments for other complex traits and diseases in humans and even in other species. Results from these kinds of analyses provide important information and a baseline for investigators seeking to identify individual genetic variants or genes with sex-, population- and/or environment-specific effects for complex traits.

Material and Methods

Bivariate GCTA-GREML analysis

We used individual-level data from seven GWAS cohorts, i.e. ARIC, NHS, HPFS, TwinGene, HRS, EGCUT and Lifelines. Informed consent was obtained from all subjects. QCs of the genotype data and imputation have been detailed elsewhere (15). In brief, we performed QCs of the genotype data in each cohort following the criteria as listed in the Supplementary Material, Table S4 of Yang et al. (15), and imputed the genotype data to the 1KGP reference panels (21) using IMPUTE2 (22). After imputation, we further excluded SNPs with Hardy–Weinberg equilibrium (HWE) test P < 1e-6, minor allele count < 3 or IMPUTE-INFO (the metric outputted from IMPUTE2 software as a measure of imputation accuracy) < 0.3. We used an IMPUTE-INFO threshold of 0.3 because a recent study (15) suggests that the proportion of variance explained by imputed variants decreases substantially when filtering variants with an IMPUTE-INFO threshold of >0.3. We included only in the analysis the variants with MAF ≥ 0.01 and those in common with HapMap Project Phase 3 reference panels (HM3) because the HM3 SNPs were optimized to capture common genetic variation (23). We used the GCTA (12,24) to estimate the genetic relatedness between all possible pairs of samples from SNP data (∼1.2 M HapMap3 SNPs passed QCs), removed one of each pair of individuals with estimated genetic relatedness > 0.05 and retained 44 126 unrelated individuals (19 323 men and 24 803 women). All the individuals are of European descent as demonstrated by the principal component analysis in the Supplementary Material, Figure S19 of Yang et al. (15). We adjusted the phenotypes for age using linear regression and standardized the residuals to z-score in each gender group of each cohort to remove differences in the mean or variance between gender groups or cohorts. We then used the bivariate GCTA-GREML approach (11,12) to estimate the genetic correlation between men and women for height and BMI. We also performed the bivariate GREML analysis for WCadjBMI, HIPadjBMI and WHRadjBMI, where WC, HIP and WHR phenotypes were adjusted by BMI in each sex group of each cohort by linear regression to capture body fat distribution independent of overall adiposity (5). The sample sizes for the analyses of WC, HIP and WHR were smaller than those for height and BMI (Table 1) because the WC and HIP phenotypes in the HRS and GENEVA-T2D cohorts were not available to us.

Bivariate LDSC regression analysis

We accessed the summary data from GWAS meta-analyses for height, BMI, WCadjBMI, HIPadjBMI and WHRadjBMI (n = up to 133 723). These analyses were stratified by sex, and the summary data are available in the public domain (see the URLs section). We used the bivariate LDSC regression method (13), which is able to estimate genetic correlation between two traits using GWAS summary data. The method is implemented in LDSC software (see the URLs section). The LDSC data (provided in accompany with the LDSC regression software tool) were estimated from the 1KGP data with a window size of ±1 Mb. For consistency, we included in the bivariate LDSC analysis the SNPs in common with those used in the bivariate GCTA-GREML analysis. We further removed SNPs with reported sample size < 10 000 and those that did not exist in the LDSC file provided by LDSC software. We retained ∼0.95 M SNPs for LDSC regression analysis. In addition, for method comparison, we also performed the bivariate LDSC regression analysis in the combined GWAS data (n = 44 126), with or without GC correction of the summary statistics, for height, BMI, WCadjBMI, HIPadjBMI and WHRadjBMI.

URLs

GIANT summary data: https://www.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium_data_files;

GCTA software: http://cnsgenomics.com/software/gcta/index.html;

LDSC software: https://github.com/bulik/ldsc

Funding

This research was supported by the Australian National Health and Medical Research Council (APP1078037), the Australian Research Council (DP130102666) and the Sylvia & Charles Viertel Charitable Foundation.

Acknowledgements

We thank William G. Hill for helpful discussion. This study makes use of data from the database of Genotypes and Phenotypes (dbGaP) under accessions phs000090.v1.p1, phs000091.v2.p1 and phs000428.v1.p1, and the EGCUT, Lifelines and TwinGene (see the Supplementary Material for the full set of acknowledgements for these data).

Conflict of Interest statement. None declared.

References

Author notes