Protein-coding variants implicate novel genes related to lipid homeostasis contributing to body-fat distribution.

Body-fat distribution is a risk factor for adverse cardiovascular health consequences. We analyzed the association of body-fat distribution, assessed by waist-to-hip ratio adjusted for body mass index, with 228,985 predicted coding and splice site variants available on exome arrays in up to 344,369 individuals from five major ancestries (discovery) and 132,177 European-ancestry individuals (validation). We identified 15 common (minor allele frequency, MAF ≥5%) and nine low-frequency or rare (MAF <5%) coding novel variants. Pathway/gene set enrichment analyses identified lipid particle, adiponectin, abnormal white adipose tissue physiology and bone development and morphology as important contributors to fat distribution, while cross-trait associations highlight cardiometabolic traits. In functional follow-up analyses, specifically in Drosophila RNAi-knockdowns, we observed a significant increase in the total body triglyceride levels for two genes (DNAH10 and PLXND1). We implicate novel genes in fat distribution, stressing the importance of interrogating low-frequency and protein-coding variants.

Central body fat distribution, as assessed by waist-to-hip ratio (WHR), is a heritable and a well-established risk factor for adverse metabolic outcomes[1-6]. Lower values of WHR are associated with lower risk of cardiometabolic diseases like type 2 diabetes (T2D)[7,8], or differences in bone structure and gluteal muscle mass[9]. These epidemiological associations are consistent with our previously reported genome-wide association study (GWAS) results of 49 loci associated with WHR (after adjusting for body mass index, WHRadjBMI)[10]. Notably, genetic predisposition to higher WHRadjBMI is associated with increased risk of T2D and coronary heart disease (CHD), which appears to be causal[9].

Recently, large-scale studies have identified ~125 common loci for multiple measures of central obesity, primarily non-coding variants of relatively modest effect[10-16]. Large scale interrogation of coding and splice site variants, including both common (minor allele frequency [MAF]≥5%) and low frequency or rare (LF/RV, MAF<5%), may lead to additional insights into the etiology of central obesity. Herein, we identify and characterize such variants associated with WHRadjBMI using ExomeChip array genotypes.

RESULTS

Protein-coding and splice site variation associations

We conducted a 2-stage fixed-effects meta-analysis testing additive and recessive models to detect protein-coding genetic variants that influence WHRadjBMI (Online Methods, Figure 1). Stage 1 included up to 228,985 variants (218,195 LF/RV) in up to 344,369 individuals from 74 studies of European, South and East Asian, African, and Hispanic/Latino descent individuals (Supplementary Data 1–3). Stage 2 assessed 70 suggestive (P < 2 × 10−6) Stage 1 variants in two cohorts, UK Biobank (UKBB) and deCODE for a total Stage 1+2 sample size of 476,546 (88% European). Of the 70 variants considered, two common and five LF/RV were not available in Stage 2 (Tables 1–2, Supplementary Data 4–6). Variants are considered novel and statistically significant if they were greater than one megabase (Mb) from a previously-identified WHRadjBMI SNP[10-16] and achieve array-wide significance (P < 2 × 10−7, Stage 1+2).

Figure 1.

Summary of meta-analysis study design and workflow.

Abbreviations:EUR- European, AFR- African, SAS- South Asian, EAS- East Asian, and HIS- Hispanic/Latino ancestry.* Novel variants include those that are >1MB from a previously published WHRadjBMIGWAS tag SNP.¥ Independent (INDEP) includes variants that are nearby known WHRadjBMI GWAS tag variants, but were determined independent after conditional analysis.

Table 1.

Association results for Combined Sexes.

Association results based on an additive or recessive model for coding variants that met array-wide significance (P < 2 × 10−7) in the sex-combined meta-analyses.

Locus (+/−1 Mb of a given variant)

Chr:Position (GRCh37)[b]

rsID

EA

OA

Gene[c]

Amino Acid Change[c]

If locus is known, nearby (< 1 MB) published variant(s)[d]

N

EAF

β[e]

SE

P-value

P-value for Sex-heterogeneity[f]

Other Criteria For Sig[h]

Abbreviations: GRCh37=human genome assembly build37;rsID=based on dbSNP; VEP=Ensembl Variant Effect Predictor toolset; GTEx=Genotype-Tissue Expression project;SD=standard deviation; SE=standard error;N=sample size; EAF=effect allele frequency; EA=effect allele; OA=other allele.

Coding variants refer to variants located in the exons and splicing junction regions.

Variant positions are reported according to Human assembly build 37 and their alleles are coded based on the positive strand.

The gene the variant falls in and amino acid change from the most abundant coding transcript is shown (protein annotation is based on VEP toolset and transcript abundance from GTEx database).

Previously published variants within +/−1Mb are from Shungin et al.[10], except for rs6976930 and rs10786152 from Graff et al.[14] and rs6499129 from Ng. et al [16].

Effect size is based on standard deviation (SD) per effect allele

P-value for sex heterogeneity, testing for difference between women-specific and men-specific beta estimates and standard errors, was calculated using EasyStrata: Winkler, T.W. et al. EasyStrata: evaluation and visualization of stratified genome-wide association meta-analysis data. Bioinformatics 2015: 31, 259–61.PMID: 25260699. Bolded P-values met significance threshold after bonferonni correction (P-value<7.14E-04; i.e. 0.05/70 variants).

rs1334576 in is a new signal in a known locus that is independent from the known signal, rs1294410; rs139745911 in is a new signal in a known locus that is independent from all known signals rs11961815, rs72959041, rs1936805, in a known locus (see Supplementary Table 4).

Each flag indicates a that a secondary criterion for significance may not be met, G- P-value > 5×10–8 (GWAS significant), C- Association Signal was not robust against collider bias; S- variant was not available in stage 2 studies for validation of Stage 1 association.

Table 2.

Association results for Sex-stratified analyses.

Association results based on an additive or recessive model for coding variants that met array-wide significance (P < 2 × 10−7) in the sex-specific meta-analyses and reach Bonferonni corrected P-value for sex heterogeneity (Psexhet< 7.14 × 10−4).

Locus (+/−1Mb of a given variant)	Chr:Position (GRCh37)[c]	rsID	EA	OA	Gene[d]	Amino Acid Change[d]	In sex-combined analyses[e]	If locus is known, nearby (< 1 MB) published variant(s)[f]	P-value for Sex-heterogeneity[g]	Men					Women					Other Criteria For Sig[j]
										N	EAF	β[h]	SE	P	N	EAF	β[h]	SE	P
Variants in Novel Loci
All Ancestry Additive model Men only analyses
1	13:96665697	rs148108950	A	G	UGGT2	P175L	No	-	1.5E-06	203,009	0.006	0.130	0.024	6.1E-08	221,390	0.004	−0.044	0.027	1.1E-01	G
2	14:23312594	rs1042704	A	G	MMP14	D273N	No	-	2.6E-04	226,646	0.202	0.021	0.004	2.6E-08	250,018	0.197	0.002	0.004	6.1E-01
All Ancestry Additive model Women only analyses
3	1:205130413	rs3851294	G	A	DSTYK	C641R	No	-	9.8E-08	225,803	0.914	−0.005	0.005	3.4E-01	249,471	0.912	0.034	0.005	4.5E-11
4	2:158412701	rs55920843	T	G	ACVR1C	N150H	Yes	-	1.7E-07	210,071	0.989	0.006	0.015	7.2E-01	245,808	0.989	0.113	0.014	1.7E-15
5	19:8429323	rs116843064	G	A	ANGPTL4	E40K	No	-	1.3E-07	203,098	0.981	−0.017	0.011	1.4E-01	243,351	0.981	0.064	0.011	1.2E-09
Variants in Previously Identified Loci
All Ancestry Additive model Women only analyses
1	1:154987704	rs141845046	C	T	ZBTB7B	P190S	Yes	rs905938	7.9E-07	226,709	0.975	0.004	0.010	6.9E-01	250,084	0.977	0.070	0.010	2.3E-13
2	2:165551201	rs7607980	T	C	COBLL1	N941D	Yes	rs1128249, rs10195252, rs12692737, rs12692738, rs17185198	3.0E-30	173,600	0.880	−0.018	0.005	5.8E-04	216,636	0.878	0.062	0.005	6.7E-39
3	3:129137188	rs62266958	C	T	EFCAB12	R197H	Yes	rs10804591	9.3E-05	226,690	0.937	0.018	0.006	3.1E-03	250,045	0.936	0.051	0.006	8.1E-18
	3:129284818	rs2625973	A	C	PLXND1	L1412V	Yes		1.6E-05	226,650	0.736	0.005	0.003	1.9E-01	250,023	0.730	0.025	0.003	8.2E-14
	3:129293256	rs2255703	T	C		M870V	Yes		5.0E-04	226,681	0.609	0.003	0.003	3.1E-01	250,069	0.602	0.018	0.003	1.9E-09
4	4:89625427	rs1804080	G	C	HERC3	E946Q	Yes	rs9991328	4.1E-06	222,556	0.839	0.008	0.004	6.6E-02	223,877	0.837	0.034	0.004	2.1E-16
	4:89668859	rs7657817	C	T	FAM13A	V443I	Yes		9.6E-05	226,680	0.816	0.006	0.004	1.5E-01	242,970	0.815	0.026	0.004	5.9E-12
5	6:127476516	rs1892172	A	G	RSPO3	synonymous	Yes	rs11961815, rs72959041, rs1936805	7.7E-09	226,677	0.541	0.018	0.003	5.6E-10	250,034	0.545	0.042	0.003	3.4E-48
	6:127767954	rs139745911[i]	A	G	KIAA0408	P504S	Yes		2.0E-04	188,079	0.010	0.057	0.017	6.8E-04	205,203	0.010	0.143	0.016	5.9E-19
6	11:64031241	rs35169799	T	C	PLCB3	S778L	Yes	rs11231693	1.3E-04	226,713	0.061	0.016	0.006	9.6E-03	250,097	0.061	0.049	0.006	6.7E-16
7	12:124265687	rs11057353	T	C	DNAH10	S228P	Yes	rs4765219, rs863750	2.7E-08	226,659	0.370	0.005	0.003	8.3E-02	250,054	0.376	0.029	0.003	3.1E-22
	12:124330311	rs34934281	C	T		T1785M	Yes		3.1E-08	226,682	0.891	0.006	0.005	1.9E-01	250,066	0.887	0.043	0.005	1.4E-20
	12:124427306	rs11057401	T	A	CCDC92	S53C	Yes		5.5E-11	223,324	0.701	0.013	0.003	4.3E-05	244,678	0.689	0.043	0.003	1.0E-41

Abbreviations: GRCh37=human genome assembly build 37;rsID=based on dbSNP; VEP=Ensembl Variant Effect Predictor toolset; GTEx=Genotype-Tissue Expression project; SD=standard deviation; SE=standard error;N=sample size; EA=effect allele; OA=other allele; EAF=effect allele frequency.

Coding variants refer to variants located in the exons and splicing junction regions.

Bonferonni corrected Pvalue for the number of SNPs tested for sex-heterogeneity is <7.14E-04 i.e. 0.05/70 variants.

Variant positions are reported according to Human assembly build 37 and their alleles are coded based on the positive strand.

The gene the variant falls in and amino acid change from the most abundant coding transcript is shown (protein annotation is based on VEP toolset and transcript abundance from GTEx database).

Variant was also identified as array-wide significant in the sex-combined analyses.

Previously published variants within +/−1Mb are from Shungin D et al. New genetic loci link adipose and insulin biology to body fat distribution. Nature 2015; 518, 187–196 doi:10.1038/nature14132 (PMID 25673412).

P-value for sex heterogeneity, testing for difference between women-specific and men-specific beta estimates and standard errors, was calculated using EasyStrata: Winkler, T.W. et al. EasyStrata: evaluation and visualization of stratified genome-wide association meta-analysis data. Bioinformatics 2015: 31, 259–61. PMID: 25260699.

Effect size is based on standard deviation (SD) per effect allele

rs139745911 in KIAA0408 is a new signal in a known locus that is independent from all known signals rs11961815, rs72959041, rs1936805, in a known locus (see Supplementary 8A/B).

Each flag indicates a that a secondary criterion for significance may not be met, G- P-value > 5×10–8 (GWAS significant), C- Association Signal was not robust against collider bias; S- variant was not availabel in Stage 2 studies for validation of Stage 1 association.

In our primary meta-analysis, including all Stage 1+2 samples, we identified 48 coding variants (16 novel) across 43 genes, 47 assuming an additive model, and one under a recessive model (Table 1, Supplementary Figures 1–4). Due to possible heterogeneity, we also performed European-only meta-analysis. Here, four additional coding variants were significant (three novel) assuming an additive model (Table 1, Supplementary Figures 5–8). Of these 52 significant variants, eleven were LF/RV and displayed larger effect estimates than many previously reported common variants[10], including seven novel variants in RAPGEF3, FGFR2, R3HDML, HIST1H1T, PCNXL3, ACVR1C, and DARS2. Variants with MAF ≤ 1% had effect sizes approximately three times greater than those of common variants (MAF ≥ 5%). Despite large sample size, we cannot rule out the possibility that additional LF/RV with smaller effects exist (See estimated 80% power in Figure 2). However, in the absence of common variants with similarly large effects, our results point to the importance of investigating LF/RV.

Figure 2.

Minor allele frequency compared to estimated effect. This scatter plot displays the relationship between minor allele frequency (MAF) and the estimated effect (β) for each significant coding variant in our meta-analyses. All novel WHRadjBMI variants are highlighted in orange, and variants identified only in models that assume recessive inheritance are denoted by diamonds and only in sex-specific analyses by triangles. Eighty percent power was calculated based on the total sample size in the Stage 1+2 meta-analysis and P = 2 × 10−7. Estimated effects are shown in original units (cm/cm) calculated by using effect sizes in standard deviation (SD) units times SD of WHR in the ARIC study (sexes combined = 0.067, men = 0.052, women = 0.080). WHR; waist-to-hip ratio

Given established sex differences in the genetic underpinnings of WHRadjBMI[10,11], we also performed sex-stratified analyses. We detected four additional novel variants that exhibit significant sex-specific effects (Psexhet < 7.14 × 10−4, Online Methods) in UGGT2 and MMP14 for men; and DSTYK and ANGPTL4 for women (Table 2, Supplementary Figures 9–15); including LF/RV in UGGT2 and ANGPTL4 (MAFmen = 0.6% and MAFwomen = 1.9%, respectively). Additionally, 14 variants from the sex-combined meta-analyses displayed significantly stronger effects in women, including the novel, LF/RV in ACVR1C (rs55920843, MAF=1.1%). Overall, 19 of the 56 variants (32%) identified across all meta-analyses (48 from all ancestry, 4 from European-only and 4 from sex-stratified analyses) showed significant sex-specific effects on WHRadjBMI: 16 variants with significantly stronger effects in women, and three in men (Figure 1).

In summary, we identified 56 array-wide significant coding variants (P < 2.0 × 10−7); 43 common (14 novel) and 13 LF/RV (9 novel). For the 55 significant variants from the additive model, we examined potential collider bias[17,18] (Online Methods, Supplementary Table 1, Supplementary Note). Overall, 51 of 55 variants were robust to collider bias[17,18]. Of these, 25 variants were nominally associated with BMI (PBMI < 0.05), yet effect sizes changed little after correction for potential biases (15% change in effect estimate on average). For four of the 55 SNPs (rs141845046, rs1034405, rs3617, rs9469913), attenuation following correction was noted (Pcorrected > 9 × 10−4, 0.05/55), including one novel variant, rs1034405 in C3orf18, demonstrating a possible overestimation of these effects in the current analysis.

Using Stage 1 results, we then aggregated LF/RV across genes and tested their joint effect with SKAT and burden tests[19] (Supplementary Table 2, Online Methods). None of the five genes that reached array-wide significance (P < 2.5 × 10−6, 0.05/16,222 genes tested: RAPGEF3, ACVR1C, ANGPTL4, DNAI1, and NOP2) remained significant after conditioning on the most significant single variant.

Conditional analyses

We next implemented conditional analyses to determine (1) the total number of independent signals identified, and (2) whether the 33 variants near known GWAS signals (< +/− 1 Mb) represent independent novel associations. We used approximate joint conditional analyses to test for independence in Stage 1 (Online Methods; Supplementary Table 3)[19]. Only the RSPO3-KIAA0408 locus contains two independent variants 291 Kb apart, rs1892172 in RSPO3 (MAF = 46.1%, Pconditional = 4.37 × 10−23 in the combined sexes, and Pconditional = 2.4 × 10−20 in women) and rs139745911 in KIAA0408 (MAF = 0.9%, Pconditional = 3.68 × 10−11 in combined sexes, and Pconditional = 1.46 × 10−11 in women; Figure 3a). For the 33 variants within one Mb of previously identified WHRadjBMI SNPs, sex-combined conditional analyses identified one coding variant representing a novel independent signal in a known locus [RREB1; Stage 1 meta-analysis, rs1334576, EAF=44%, Pconditional = 3.06 × 10−7, (Supplementary Data 7, Figure 3b); UKBB analysis, rs1334576, RREB1, Pconditional = 1.24 × 10−8, (Supplementary Table 4).

Figure 3.

Regional association plots for known loci with novel coding signals identified by conditional analyses. Point color reflects r2 calculated from the ARIC dataset. In a) there are two independent variants in RSPO3 and KIAA0408, based on results from the stage 1 All Ancestry women (N = 180,131 for RSPO3 and 139,056 for KIAA0408). In b) we have a variant in RREB1 that is independent of the GWAS variant rs1294421, based on results from the stage 1 All Ancestry sex-combined individuals (N = 319,090).

In summary, we identified 56 WHRadjBMI-associated coding variants in 41 independent association signals, 24 of which are new or independent of known GWAS-identified tag SNPs (either > 1 MB +/− or array-wide significant following conditional analyses) (Figure 1). Thus, we identified 15 common and 9 LF/RV novel and independent variants following conditional analyses.

Gene set and pathway enrichment analysis

To determine if significant coding variants highlight novel or previously identified biological pathways, we applied two complementary methods, EC-DEPICT (ExomeChip Data-driven Expression Prioritized Integration for Complex Traits)[20,21] and PASCAL[22] (Online Methods). For PASCAL all variants were used, for EC-DEPICT we examined only 361 variants with suggestive significance (P < 5 × 10−4)[10,23] from the all ancestries combined sexes analysis (which after clumping and filtering became 101 lead variants in 101 genes). We separately analyzed variants that exhibited significant sex-specific effects (Psexhet < 5 × 10−4).

The sex-combined analyses identified 49 significantly enriched gene sets (FDR < 0.05) that grouped into 25 meta-gene sets (Supplementary Note, Supplementary Data 8–9). We noted a cluster of meta-gene sets with direct relevance to metabolic aspects of obesity (“enhanced lipolysis,” “abnormal glucose homeostasis,” “increased circulating insulin level,” and “decreased susceptibility to diet-induced obesity”); we observed two significant adiponectin-related gene sets within these meta-gene sets. While these pathway groups had previously been identified in the GWAS DEPICT analysis (Figure 4), many of the individual gene sets within these meta-gene sets were not significant in the previous GWAS analysis, such as “insulin resistance,” “abnormal white adipose tissue physiology,” and “abnormal fat cell morphology” (Supplementary Data 8, Figure 4, Supplementary Figure 16a), but represent similar biological underpinnings implied by the shared meta-gene sets. Despite their overlap with the GWAS results, these analyses highlight novel genes that fall outside known GWAS loci and with strong contributions to the significantly enriched gene sets related to adipocyte and insulin biology (e.g. MLXIPL, ACVR1C, and ITIH5) (Figure 4).

Figure 4.

Heat maps showing DEPICT gene set enrichment results from the stage 1 All Ancestry sex-combined individuals (N = 344,369). For any given square, the color indicates how strongly the corresponding gene (x-axis) is predicted to belong to the reconstituted gene set (y-axis). This value is based on the gene’s z-score for gene set inclusion in DEPICT’s reconstituted gene sets, where red indicates a higher and blue a lower z-score. To visually reduce redundancy and increase clarity, we chose one representative “meta-gene set” for each group of highly correlated gene sets based on affinity propagation clustering (Online Methods, Supplementary Note). Heatmap intensity and DEPICT P-values (Supplementary Data 8–9) correspond to the most significantly enriched gene set within the meta-gene set. Annotations for the genes indicate (1) the minor allele frequency of the significant ExomeChip (EC) variant (blue; if multiple variants, the lowest-frequency variant was kept), (2) whether the variant’s P-value reached array-wide significance (< 2 × 10−7) or suggestive significance (< 5 × 10–4) (shades of purple), (3) whether the variant was novel, overlapping “relaxed” GWAS signals from Shungin et al.[10] (GWAS P < 5 × 10−4), or overlapping “stringent” GWAS signals (GWAS P < 5 × 10−8) (pink), and (4) whether the gene was included in the gene set enrichment analysis or excluded by filters (shades of brown/orange) (Online Methods, Supplementary Note). Annotations for the gene sets indicate if the meta-gene set was found significant (shades of green; FDR < 0.01, < 0.05, or not significant) in the DEPICT analysis of GWAS results from Shungin et al.[10]

Also, we conducted pathway analyses after excluding variants from previous WHRadjBMI analyses[10] (Supplemental Note). Seventy-five loci/genes were included in the EC-DEPICT analysis, and we identified 26 significantly enriched gene sets (13 meta-gene sets). Here, all but one gene set, “lipid particle size”, were related to skeletal biology, likely reflecting an effect on the pelvic skeleton (hip circumference), shared signaling pathways between bone and fat (such as TGF-beta) and shared developmental origin[24] (Supplementary Data 9, Supplementary Figure 16b). These previously identified GWAS DEPICT significant findings provide a fully independent replication of their biological relevance for WHRadjBMI.

We used PASCAL (Online Methods) to further distinguish between enrichment based on coding-only variant associations (this study) and regulatory-only variant associations (up to 20 kb upstream of the gene from a previous GIANT study[10]), finding 116 significantly enriched coding pathways (FDR < 0.05; Supplementary Data 10). We also compared the coding pathways to those identified in the total previous GWAS effort (using both coding and regulatory variants) identifying a total of 158 gene sets. Forty-two gene sets were enriched in both analyses, and we found high concordance in the -log10 (p-values) between ExomeChip and GWAS gene set enrichment [Pearson’s r (coding vs regulatory only) = 0.38, P < 10−300; Pearson’s r (coding vs coding+regulatory) = 0.51, P < 10−300)]. Nonetheless, some gene sets were enriched specifically for variants in coding regions (e.g., decreased susceptibility to diet-induced obesity, abnormal skeletal morphology) or unique to variants in regulatory regions (e.g. transcriptional regulation of white adipocytes) (Supplementary Figure 17).

The EC-DEPICT and PASCAL results showed a moderate but strongly significant correlation (for EC-DEPICT and the PASCAL max statistic, r = .28, P = 9.8 × 10−253; for EC-DEPICT and the PASCAL sum statistic, r = 0.287, P = 5.42 × 10−272). Common gene sets strongly implicate a role for skeletal biology, glucose homeostasis/insulin signaling, and adipocyte biology (Supplementary Figure 18).

Cross-trait associations

To assess the clinical relevance of our identified variants with cardiometabolic, anthropometric, and reproductive traits, we conducted association lookups from existing ExomeChip studies of 15 traits (Supplementary Data 11, Supplementary Figure 19).[21,25-29] Variants in STAB1 and PLCB3 displayed the greatest number of significant associations with seven different traits (P < 9.8 × 10−4, 0.05/51 variants tested). Also, these two genes cluster together with RSPO3, DNAH10, MNS1, COBLL1, CCDC92, and ITIH3. The WHR-increasing alleles in this cluster exhibit a previously described pattern of increased cardiometabolic risk (e.g. increased fasting insulin, two-hour glucose [TwoHGlu], and triglycerides; and decreased high-density lipoprotein cholesterol [HDL]), but also decreased BMI.[30-36] The impact of central obesity may be causal, as a 1 SD increase in genetic risk of central adiposity was previously associated with higher total cholesterol, triglycerides, fasting insulin and TwoHGlu, and lower HDL.[9]

We conducted a search in the NHGRI-EBI GWAS Catalog[37,38] to determine if our variants are in high LD (R2 > 0.7) with variants associated with traits or diseases not covered by our cross trait lookups (Supplementary Data 12). We identified several cardiometabolic traits (adiponectin, coronary heart disease, etc.), diet/behavioral traits potentially related to obesity (carbohydrate, fat intake, etc.), behavioral and neurological traits (schizophrenia, bipolar disorder, etc.), and inflammatory or autoimmune diseases (Crohn’s Disease, multiple sclerosis, etc.).

Given the established correlation between total body fat percentage and WHR of up to 0.483[39-41], we examined the association of our top exome variants with both total body fat percentage (BF%) and truncal fat percentage (TF%) available in a sub-sample of UKBB (N = 118,160) (Supplementary Tables 5–6). Seven of the common novel variants were significantly associated (P < 0.001, 0.05/48 variants examined) with both BF% and TF% in the sexes-combined analysis (COBLL1, UHRF1BP1, WSCD2, CCDC92, IFI30, MPV17L2, IZUMO1) and two with TF% in women only (EFCAB12, GDF5). Only rs7607980 in COBLL1 is near a known BF% GWAS locus (rs6738627; R = 0.1989, distance = 6,751 bp, with our tag SNP)[42]. Of the nine SNPs associated with at least one of these two traits, all variants displayed much greater magnitude of effect on TF% compared to BF% (Supplementary Figure 20).

Previous studies have demonstrated the importance of examining common and LF/RV within genes with mutations known to cause monogenic diseases.[43,44] Thus, we assessed enrichment of WHRadjBMI variants within monogenic lipodystrophy and/or insulin resistance genes.[43,44] (Supplementary Data 13). No significant enrichment was observed, possibly due in part to the small number of implicated genes and the relatively small number of variants in monogenic disease-causing genes (Supplementary Figure 21).

Genetic architecture of WHRadjBMI coding variants

We used summary statistics from our Stage 1 primary meta-analysis results to estimate the phenotypic variance explained by subsets of variants across various significance thresholds (P < 2 × 10−7 to 0.2) and conservatively using only independent SNPs (Supplementary Table 7, Online Methods, and Supplementary Figure 22). For only independent coding variants that reached suggestive significance in Stage 1 (P < 2 × 10−6), 33 SNPs explain 0.38% of the variation. The 1,786 independent SNPs with a liberal threshold of P<0.02 explain 13 times more variation (5.12%), however, these large effect estimates may be subject to winner’s curse. When considering all coding variants on the ExomeChip in combined sexes, 46 SNPs with a P < 2 × 10−6 and 5,917 SNPs with a P < 0.02 explain 0.51% and 13.75% of the variance in WHRadjBMI, respectively. As expected given the design of the ExomeChip, the majority of the variance explained is attributable to rare and low frequency coding variants. However, for LF/RVs, those that passed significance in Stage 1 explain only 0.10% of the variance in WHRadjBMI. We also estimated variance explained for the same SNPs in women and men separately and observed a greater variance explained in women compared to men (PRsqDiff < 0.002 = 0.05/21, Bonferroni-corrected threshold) at each significance threshold considered (differences ranged from 0.24% to 0.91%).

We conducted penetrance analysis using the UKBB (both sexes combined, and men- and women-only) to determine if there is a significant accumulation of the minor allele in either the centrally obese or non-obese groups (Online Methods). Three rare variants (MAF ≤ 1%) with larger effect sizes (effect size > 0.90) were included in the penetrance analysis using World Health Organization cut-offs for central obesity. Of these, one SNP (rs55920843-ACVR1C; Psex-combined = 9.25 × 10−5; Pwomen= 4.85 × 10−5) showed a statistically significant difference in the number of carriers and non-carriers of the minor allele in the combined and female-only analysis (sex-combined obese carriers = 2.2%; non-obese carriers = 2.6%; women obese carriers = 2.1%; non-obese women carriers = 2.6%, Supplementary Table 8, Supplementary Figure 23).

Drosophila Knockdown

Considering the genetic evidence of adipose and insulin biology in determining body fat distribution[10], and the lipid signature of the variants described herein, we examined whole-body triglyceride levels in adult Drosophila, a model organism in which the fat body is an organ functionally analogous to mammalian liver and adipose tissue as triglycerides are the major source of fat storage[45]. Of the 51 genes harboring our 56 significantly associated variants, we identified 27 Drosophila orthologues for functional follow-up analyses. We selected genes with large changes in triglyceride levels (> 20% increase or > 40% decrease, as chance alone is unlikely to cause changes of this magnitude) from an existing large-scale screen with ≤2 replicates per knockdown strain.[45] Two orthologues, for PLXND1 and DNAH10, met these criteria and were subjected to additional knockdown experiments with ≥5 replicates using tissue-specific drivers (fat body [cg-Gal4] and neuronal [elav-Gal4] specific RNAi-knockdowns) (Supplementary Table 9). A significant (P < 0.025, 0.05/2 orthologues) increase in the total body triglyceride levels was observed in DNAH10 orthologue knockdown strains for both the fat body and neuronal drivers. Only the neuronal driver knockdown for PLXND1 produced a significant change in triglyceride storage. DNAH10 and PLXND1 both lie within previous GWAS identified regions. Adjacent genes have been highlighted as likely candidates for the DNAH10 association region, including CCDC92 and ZNF664 based on expression quantitative trait locus (eQTL) evidence. Of note, rs11057353 in DNAH10 showed suggestive significance after conditioning on the known GWAS variants in nearby CCDC92 (sex-combined Pconditional = 7.56 × 10−7; women-only rs11057353 Pconditional = 5.86 × 10−7, Supplementary Table 4) thus providing some evidence of multiple causal variants/genes underlying this signal. Further analyses are needed to determine whether the implicated coding variants from the current analysis are the putatively functional variants.

eQTL Lookups

We examined the cis-association of variants with expression level of nearby genes in subcutaneous and visceral omental adipose, skeletal muscle, and pancreas tissue from the Genotype-Tissue Expression (GTEx)[46] project, and assessed whether exome and eQTL associations implicated the same signal (Online Methods, Supplementary Data 14–15). The lead exome variant was associated with expression level of the gene itself for DAGLB, MLXIPL, CCDC92, MAPKBP1, LRRC36 and UQCC1. However, for MLXIPL, MAPKBP1, and LRRC36, the lead variant is also associated with expression of additional nearby genes. At three additional loci, the lead exome variant is only associated with expression level of nearby genes (HEMK1 at C3orf18; NT5DC2, SMIM4 and TMEM110 at STAB1/ITIH3; and C6orf106 at UHRF1BP1). Thus, although detected with a missense variant, these results are also consistent with a regulatory mechanism of effect, and the association signal may well be due to linkage disequilibrium (LD) with nearby regulatory variants.

Some of the coding genes implicated by eQTL analyses are known to be involved in adipocyte differentiation or insulin sensitivity: e.g. for MLXIPL, the encoded carbohydrate responsive element binding protein is a transcription factor, regulating glucose-mediated induction of de novo lipogenesis in adipose tissue, and expression of its beta-isoform in adipose tissue is positively correlated with adipose insulin sensitivity[47,48]. For CCDC92, the reduced adipocyte lipid accumulation upon knockdown confirmed the involvement of its encoded protein in adipose differentiation[49].

Biological Curation

To investigate the possible functional role of the identified variants, we conducted thorough searches of the literature and publicly available bioinformatics databases (Supplementary Data 16–17, Box 1, Online Methods). Many of our novel LF/RV are in genes that are intolerant of nonsynonymous mutations (e.g. ACVR1C, DARS2, FGFR2; ExAC Constraint Scores >0.5). Other coding variants lie within genes that are involved in glucose homeostasis (e.g. ACVR1C, UGGT2, ANGPTL4), angiogenesis (RASIP1), adipogenesis (RAPGEF3), and lipid biology (ANGPTL4, DAGLB).

DISCUSSION

Our analysis of coding variants from ExomeChip data in up to 476,546 individuals identified a total of 56 array-wide significant WHRadjBMI associated variants in 41 independent association signals, including 24 newly identified (23 novel and one independent of known GWAS signals). Nine of these variants were LF/RV, indicating an important role for such variants in the polygenic architecture of fat distribution. While, due to their rarity, these coding variants explain a small proportion of the trait variance at a population level, they may be more functionally tractable than non-coding variants and have a critical impact at the individual level. For instance, the association between a LF/RV (rs11209026; R381Q; MAF < 5% in ExAC) located in the IL23R gene and multiple inflammatory diseases[50-53] led to development of new therapies targeting IL23 and IL12 in the same pathway.[54-56] Thus, we are encouraged that our LF/RV displayed large effect sizes; all but one of the nine novel LF/RV display larger effects than the 49 SNPs reported in Shungin et al. 2015[10], and some of these effects were up to 7-fold larger than those previously reported for GWAS. This finding mirrors results for other cardiometabolic traits[57], and suggests variants of possible clinical significance with even larger effect and rarer variants will likely be detected with greater sample sizes.

We continue to observe sexual dimorphism in the genetic architecture of WHRadjBMI[11]. We identified 19 coding variants with significant sex differences, of which 16 (84%) display larger effects in women compared to men. Of the variants outside of GWAS loci, we reported three (two LF/RV) that show a significantly stronger effect in women and two (one LF/RV) that show a stronger effect in men. Genetic variants continue to explain a higher proportion of the phenotypic variation in body fat distribution in women compared to men.[10,11] Of the novel female (DSTYK and ANGPTL4) and male (UGGT2 and MMP14) specific signals, only ANGPTL4 implicated fat distribution related biology associated with both lipid biology and cardiovascular traits (Box 1). Sexual dimorphism in fat distribution is apparent[58−60] and at sexually dimorphic loci, hormones with different levels in men and women may interact with genomic and epigenomic factors to regulate gene activity, though this remains to be tested. Dissecting the underlying molecular mechanisms of the sexual dimorphism in body fat distribution, and how it is correlated with – and causing – important comorbidities like cardiometabolic diseases will be crucial for improved understanding of disease pathogenesis.

Overall, we observe fewer significant associations, pathways, and cross-trait associations between WHRadjBMI and coding variants on the ExomeChip than Turcot et al for BMI[25]. One reason for this may be smaller sample size (NWHRadjBMI = 476,546, NBMI = 718,639), and thus, lower statistical power. Power is likely not the only contributing factor, as trait architecture, heritability (possibly overestimated in some phenotypes), and phenotype precision all likely contribute to our study’s capacity to identify LF/RV with large effects. Further, it is possible that the comparative lack of significant findings for WHRadjBMI may be a result of higher selective pressure against genetic predisposition to cardiometabolic phenotypes, thus rarer risk variants.[61] The ExomeChip is limited by the variants present on the chip, which was largely dictated by sequencing studies in European-ancestry populations and MAF detection criteria of ~0.012%. It is likely that through increased sample size, use of chips designed to detect variation across a range of continental ancestries, and high quality, deep imputation with large reference samples future studies will detect additional variation from the entire allele frequency spectrum that contributes to fat distribution.

The collected genetic and epidemiologic evidence has demonstrated that increased central adiposity is correlated with risk of T2D and CVD, and that this association is likely causal with potential mediation through blood pressure, triglyceride-rich lipoproteins, glucose, and insulin[9]. This observation yields an immediate follow-up question: Which mechanisms regulate depot-specific fat accumulation and are risks for disease driven by increased visceral and/or decreased subcutaneous adipose tissue mass. Pathway analysis identified several novel pathways and gene sets related to metabolism and adipose regulation, bone growth and development and adiponectin, a hormone which has been linked to “healthy” expansion of adipose tissue and insulin sensitivity.[62] Similarly, expression/eQTL results support the relevance of adipogenesis, adipocyte biology, and insulin signaling, supporting our previous findings for WHRadjBMI.[10] We also provide evidence suggesting known biological functions and pathways contributing to body fat distribution (e.g., diet-induced obesity, angiogenesis, bone growth/morphology, and lipolysis).

The ultimate aim of genetic investigations of obesity-related traits is to identify dysregulated genomic pathways leading to obesity pathogenesis that may result in a myriad of downstream illnesses. Thus, our findings may enhance the understanding of central obesity and identify new molecular targets to avert its negative health consequences. Significant cross-trait associations are consistent with expected direction of effect for several traits, i.e. the WHR-increasing allele is associated with higher values of triglycerides, DBP, fasting insulin, total cholesterol, LDL and T2D across many significant variants. However, it is worth noting that there are some exceptions. For example, rs9469913-A in UHRF1BP1 is associated with both increased WHRadjBMI and increased HDL. Also, we identified two variants in MLXIPL (rs3812316 and rs35332062), a well-known lipids-associated locus, in which the WHRadjBMI-increasing allele also increases all lipid levels, risk for hypertriglyceridemia, SBP and DBP. However, our findings show a significant and negative association with HbA1C, and nominally significant and negative associations with two-hour glucose, fasting glucose, and Type 2 diabetes, and potential negative associations with biomarkers for liver disease (e.g. gamma glutamyl transpeptidase). Other notable exceptions include ITIH3 (negatively associated with BMI, HbA1C, LDL and SBP), DAGLB (positively associated with HDL), and STAB1 (negatively associated with total cholesterol, LDL, and SBP). Therefore, caution in selecting pathways for therapeutic targets is warranted; we must look beyond the effects on central adiposity to the potential cascading effects of related diseases.

A major finding from this study is the importance of lipid metabolism for body fat distribution. In fact, pathway analyses that highlight enhanced lipolysis, cross-trait associations with circulating lipid levels, existing biological evidence from the literature, and knockdown experiments in Drosophila, point to novel candidate genes (ANGPTL4, ACVR1C, DAGLB, MGA, RASIP1, and IZUMO1) and new candidates in known regions (DNAH10[10] and MLXIPL[14]) related to lipid biology and their role in fat storage. ACVR1C, MLXIPL, and ANGPTL4, all of which are involved in lipid homeostasis, all are excellent candidate genes for central adiposity. Carriers of inactivating mutations in ANGPTL4 (Angiopoietin Like 4), for example, display low triglycerides and low risk of coronary artery disease[63]. ACVR1C encodes the activin receptor-like kinase 7 protein (ALK7), a receptor for the transcription factor TGFB-1, well known for its central role in general growth and development[64-68], and adipocyte development particularly[68]. ACVR1C exhibits the highest expression in adipose tissue, but is also highly expressed in the brain[69-71]. In mice, decreased activity of ACVR1C upregulates PPARγ and C/EBPα pathways and increases lipolysis in adipocytes, thus decreasing weight and diabetes.[69,72,73] Such activity suggests a role for ALK7 in adipose tissue signaling and a possible therapeutic target. MLXIPL, also important for lipid metabolism and postnatal cellular growth, encodes a transcription factor which activates triglyceride synthesis genes in a glucose-dependent manner.[74,75] The lead exome variant in MLXIPL is highly conserved, most likely damaging, and associated with reduced MLXIPL expression in adipose tissue. Furthermore, in a recent longitudinal, in vitro transcriptome analysis of adipogenesis in human adipose-derived stromal cells, gene expression of MLXIPL was up-regulated during the maturation of adipocytes, suggesting a critical role in the regulation of adipocyte size and accumulation.[76] However, given our cross-trait associations with variants in MLXIPL and diabetes-related traits, development of therapeutic targets must be approached cautiously.

Our 24 novel variants for WHRadjBMI highlight the importance of lipid metabolism in the genetic underpinnings of body fat distribution. We continue to demonstrate the critical role of adipocyte biology and insulin resistance for central obesity and offer support for potentially causal genes underlying previously identified fat distribution loci. Notably, our findings offer potential new therapeutic targets for intervention in the risks associated with abdominal fat accumulation and represents a major advance in our understanding of the underlying biology and genetic architecture of central adiposity.

ONLINE METHODS

Studies

Stage 1 included 74 studies (12 case/control, 59 population-based, and five family) comprising 344,369 adults of European (N=288,492), African (N=15,687), South Asian (N=29,315), East Asian (N=6,800), and Hispanic (N=4,075) descent. Stage 1 meta-analyses were conducted in each ancestry and in all ancestries together, for both sex-combined and sex-specific analyses. Follow-up analyses were performed in 132,177 individuals of European ancestry from deCODE and the UK Biobank, Release 1[112] (UKBB) (Supplementary Data 1–3). Informed consent was obtained by the parent study and protocols approved by each study’s institutional review boards.

Phenotypes

For each study, WHR (waist circumference divided by hip circumference) was corrected for age, BMI, and genomic principal components (derived from GWAS data, the variants with MAF >1% on the ExomeChip, and ancestry informative markers available on the ExomeChip), as well as any additional study-specific covariates (e.g. recruiting center), in a linear regression model. For studies with unrelated individuals, residuals were calculated separately by sex, whereas for family-based studies sex was included as a covariate in models with both men and women. Residuals for case/control studies were calculated separately. Finally, residuals were inverse normal transformed and used as the outcome in association analyses. Phenotype descriptives by study are shown in Supplementary Data 3.

Genotypes and QC

The majority of studies followed a standardized protocol and performed genotype calling using the algorithms indicated in Supplementary Data 2, which typically included zCall[3]. For 10 studies from the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium, raw intensity data for samples from seven genotyping centers were combined for joint calling[4]. Study-specific quality control (QC) of the genotyped variants was implemented before association analysis (Supplementary Data 1–2). To assess whether any significant associations with rare and low-frequency variants could be due to allele calling in smaller studies, we performed a sensitivity meta-analysis of all large studies (>5,000 participants) compared to all studies. We found very high concordance for effect sizes, suggesting that smaller studies do not bias our results (Supplementary Fig. 24).

Study-level statistical analyses

Each cohort performed single variant analyses for both additive and recessive models in each ancestry, for sexes combined and sex-specific groups, with either RAREMETALWORKER (see URLs) or RVTESTs (see URLs). to associate inverse normal transformed WHRadjBMI with genotype accounting for cryptic relatedness (kinship matrix) in a linear mixed model. Both programs perform score-statistic rare-variant association analysis, accommodate unrelated and related individuals, and provide single-variant results and variance-covariance matrices. The covariance matrix captures linkage disequilibrium (LD) between markers within 1 Mb, which is used for gene-level meta-analyses and conditional analyses[113,114].

Centralized quality-control

Individual cohorts identified ancestry outliers based on 1000 Genomes Phase 1 reference populations. A centralized QC procedure implemented in EasyQC[115] was applied to individual cohort summary statistics to identify cohort-specific problems: (1) possible errors in phenotype residual transformation; (2) strand issues, and (3) inflation due to population stratification, cryptic relatedness and genotype biases.

Meta-analyses

Meta-analyses were carried out in parallel by two analysts at two sites using RAREMETAL[113]. We excluded variants if they had call rate <95%, Hardy-Weinberg equilibrium P-value <1×10−7, or large allele frequency deviations from reference populations (>0.6 for all ancestries analyses and >0.3 for ancestry-specific population analyses). We also excluded markers not present on the Illumina ExomeChip array 1.0, Y-chromosome and mitochondrial variants, indels, multiallelic markers, and problematic variants based on Blat-based sequence alignment. Significance for single-variant analyses was defined at an array-wide level (P<2×10−7). For all suggestive significant variants (P<2×10−6) from Stage 1, we calculated Psexhet for each SNP, testing for differences between women-specific and men-specific beta estimates and standard errors using EasyStrata[116]. Each SNP that reached Psexhet<0.05/# of variants tested (70 variants brought forward from Stage 1, Psexhet<7.14×10−4) was considered significant. Additionally, while each individual study was asked to perform association analyses stratified by ancestry and adjusted for population stratification, all study-specific summary statistics were combined in our all ancestry meta-analyses. To investigate potential heterogeneity across ancestries, we examined ancestry-specific meta-analysis results for our top 70 variants from Stage 1 and found no evidence of significant across-ancestry heterogeneity for any of our top variants (I2 values noted in Supplementary Data 4–6).

For the gene-based analyses, we applied two sets of criteria to select variants with a MAF<5% within each ancestry based on coding variant annotation from five prediction algorithms (PolyPhen2, HumDiv and HumVar, LRT, MutationTaster, and SIFT)[117]. Our broad gene-based tests included nonsense, stop-loss, splice site, and missense variants annotated as damaging by at least one algorithm mentioned above. Our strict gene-based tests included only nonsense, stop-loss, splice site, and missense variants annotated as damaging by all five algorithms. These analyses were performed using the sequence kernel association test (SKAT) and variable threshold (VT) methods in RAREMETAL[113]. Statistical significance for gene-based tests was set at a Bonferroni-corrected threshold of P<2.5×10−6 (0.05/~20,000 genes).

Genomic inflation

We observed marked genomic inflation of the test statistics even after controlling for population stratification arising mainly from common markers; λGC in the primary meta-analysis (combined ancestries and combined sexes) was 1.06 for all variants and 1.37 for common coding and splice site markers, respectively (Supplementary Figures 3, 7 and 13, Supplementary Table 10). Such inflation is expected for a highly polygenic trait like WHRadjBMI, for studies using a non-random set of variants across the genome, and is consistent with our very large sample size[115,118,119].

The RAREMETAL R-package[113] was used to identify independent WHRadjBMI association signals across all ancestries and European meta-analysis results. RAREMETAL performs conditional analyses using covariance matrices to distinguish true signals from shadows of adjacent significant variants in LD. First, we identified lead variants (P<2×10−7) based on a 1Mb window centered on the most significant variant. We then conditioned on the lead variants in RAREMETAL and kept new lead signals at P<2×10−7 for conditioning in a second round of analysis. The process was repeated until no additional signal emerged below the pre-specified P-value threshold (P<2×10−7).

To test if the associations detected were independent of previously published WHRadjBMI variants [10,14,16], we used RAREMETAL to perform conditional analyses in the Stage 1 discovery set if the GWAS variant or its proxy (r2≥0.8) was on the ExomeChip. All variants identified in our meta-analysis and the previously published variants were available in the UKBB dataset[112], which was used as a replacement dataset if a good proxy was not on the ExomeChip. All conditional analyses in the UKBB were performed using SNPTEST[120-122].

The conditional analyses were carried out reciprocally, conditioning on the ExomeChip variant and then the previously published variant. An association was considered independent if it was significant prior to conditional analysis (P<2×10−7) with both the exome chip variant and the previously published variant, and the observed association with our variant remained significant upon conditional analysis. Conditional p-values between 9×10−6 and 0.05 was considered inconclusive, while those < 9×10−6 were considered suggestive.

Stage 2 meta-analyses

In Stage 2, we sought to validate 70 Stage 1 variants (P<2×10−6) in two independent studies, UKBB (N=119,572) and deCODE (N=12,605), using the same QC and analytical methodology. Genotyping, study descriptions and phenotype descriptives are provided in Supplementary Data 1–3. Stage 1+ 2 meta-analysis was performed using the inverse-variance weighted fixed effects method. Significant associations were defined as those nominally significant (P<0.05) in Stage 2 when available in Stage 2, and array-wide significance for Stage 1+2 at P<2×10−7 (0.05/~250,000 246,328 variants tested). Variants are considered novel and statistically significant if they were greater than one megabase (Mb) from a previously-identified WHRadjBMI lead SNP[10-16] and achieved a significance threshold of P<2×10−7.

Pathway enrichment analyses: EC-DEPICT

We adapted DEPICT, a gene set enrichment analysis method for GWAS data, for use with the ExomeChip (‘EC-DEPICT’) described further in a companion manuscript[21]. DEPICT uses “reconstituted” gene sets, where different types of gene sets (e.g. canonical pathways, protein-protein interaction networks, and mouse phenotypes) were extended through large-scale microarray data (see Pers et al.[20] for details). EC-DEPICT computes p-values based on Swedish ExomeChip data (Malmö Diet and Cancer (MDC), All New Diabetics in Scania (ANDIS), and Scania Diabetes Registry (SDR) cohorts, N=11,899) and, unlike DEPICT, takes as input only genes directly containing significant (coding) variants rather than all genes within a specified LD (Supplementary Note).

Two analyses were performed for WHRadjBMI ExomeChip: one with all variants p<5×10−4 (49 significant gene sets in 25 meta-gene sets, FDR <0.05) and one with all variants > 1 Mb from known GWAS loci[10] (26 significant gene sets in 13 meta-gene sets, FDR <0.05). Affinity propagation clustering[123] was used to group highly correlated gene sets into “meta-gene sets”; for each meta-gene set, the member gene set with the best p-value was used for visualization (Supplementary Note). EC-DEPICT was written in Python (see URLs).

Pathway enrichment analyses: PASCAL

We also applied PASCAL pathway analysis[22] to summary statistics from Stage 1 for all coding variants. PASCAL derives gene-based scores (SUM and MAX) and tests for over-representation of high gene scores in predefined biological pathways. We performed both MAX and SUM estimations for pathway enrichment. MAX is sensitive to genesets driven by a single signal, while SUM is better for multiple variant associations in the same gene. We used standard pathway libraries from KEGG, REACTOME and BIOCARTA, and also added dichotomized (Z-score>3) reconstituted gene sets from DEPICT[20]. To accurately estimate SNP-by-SNP correlations even for rare variants, we used the UK10K data (TwinsUK[124] and ALSPAC[125], N=3781). To distinguish contributions of regulatory and coding variants, we also applied PASCAL to summary statistics of only regulatory variants (20 kb upstream) and regulatory+coding variants from the Shungin et al[10] study. In this way, we could investigate what is gained by analyzing coding variants.

Monogenic obesity enrichment analyses

We compiled two lists consisting of 31 genes with strong evidence that disruption causes monogenic forms of insulin resistance or diabetes; and eight genes with evidence that disruption causes monogenic forms of lipodystrophy. To test for association enrichment, we conducted simulations by matching each gene with others based on gene length and number of variants tested to create 1,000 matched gene sets and assessed how often the number of variants exceeding set significance thresholds was greater than in our monogenic obesity gene set.

Variance explained

We estimated phenotypic variance explained by Stage 1 associations in all ancestries for men, women, and combined sexes[126]. For each associated region, we pruned subsets of SNPs within 500 kb of SNPs with the lowest P-value and used varying P-value thresholds (ranging from 2×10−7 to 0.02) from the combined sexes results. Additionally, we examined all variants and independent variants across a range of MAFs. The variance explained by each subset of SNPs in each stratum was estimated by summing the variance explained by individual top coding variants. To compare variance explained between men and women, we tested for significant differences assuming the weighted sum of χ2-distributed variables tend to a Gaussian distribution following Lyapunov’s central limit theorem.[126,127]

Cross-trait lookups

To evaluate relationships between WHRadjBMI and related cardiometabolic, anthropometric, and reproductive traits, association results for the 51 WHRadjBMI coding SNPs were requested from seven consortia, including ExomeChip data from GIANT (BMI, height), Global Lipids Genetics Consortium (GLGC) (total cholesterol, triglycerides, HDL-cholesterol, LDL-cholesterol), International Consortium for Blood Pressure (IBPC)[128] (systolic and diastolic blood pressure), Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC) (glycemic traits), and DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) consortium (type 2 diabetes). )[21,25-29]. For coronary artery disease, we accessed 1000 Genomes Project-imputed GWAS data released by CARDIoGRAMplusC4D[129] and for age at menarche and menopause, we used a combination of ExomeChip and 1000 Genomes Project-imputed GWAS data from ReproGen. Heatmaps were generated with gplots (R v3.3.2) using Euclidean distance based on p-value and direction of effect and complete linkage clustering (see URLs).

GWAS Catalog Lookups

To determine if significant coding variants were associated with any related cardiometabolic or anthropometric traits, we also searched the NHGRI-EBI GWAS Catalog for previous associations near our lead SNPs (+/− 500 kb). We used PLINK to calculate LD using ARIC European participants. All GWAS Catalog SNPs within the specified regions with an r2 > 0.7 were evaluated[37]. Consistent direction of effect was based on WHRadjBMI-increasing allele, LD, and allele frequency. We do not comment on direction of effect when a GWAS Catalog variant was not identical or in high LD (r2 > 0.9) with the WHR variant, and MAF >45%.

Body-fat percentage associations

We performed body fat percent and truncal fat percent look-ups of 48 of the 56 WHRadjBMI identified variants (Tables 1 and 2) available in UKBB. GWAS for body fat percent and truncal fat percent in UKBB excluded pregnant or possibly pregnant women, individuals with BMI < 15, or those without genetically confirmed European ancestry, resulting in a sample size of 120,286. Estimated body fat percent and truncal fat percent were obtained using the Tanita BC418MA body composition analyzer (Tanita, Tokyo, Japan). Participants were non-fasting and did not follow any specific instructions prior to bioimpedance measurements. SNPTEST was used to perform the analyses based on residuals adjusted for age, 15 principal components, assessment center and the genotyping chip[120].

Collider bias

To evaluate SNPs for possible collider bias[17], we used results from a recent GIANT BMI GWAS[25]. For each significant SNP from our additive models, WHRadjBMI associations were corrected for potential bias due to associations between each variant and BMI (Supplementary Note). Variants meeting Bonferroni-corrected significance (Pcorrected<9.09×10−4, 0.05/55 variants examined) were considered robust against collider bias.

Drosophila RNAi knockdown experiments

For each gene with WHRadjBMI-associated coding variants in the final combined meta-analysis (P < 2×10−7), its corresponding Drosophila orthologues were identified in the Ensembl ortholog database (see URLs), when available. Drosophila triglyceride content values were mined from a publicly available genome-wide fat screen data set[45] to identify potential genes for follow-up knockdowns. Estimated values represent fractional changes in triglyceride content in adult male flies. Data are from male progeny of crosses between male UAS-RNAi flies from the Vienna Drosophila Resource Center (VDRC) and Hsp70-GAL4; Tub-GAL8ts virgin females. (Supplementary Note). The screen comprised one to three biological replicates. We followed up each gene with a >0.2 increase or >0.4 decrease in triglyceride content.

Orthologues for two genes were brought forward for follow-up, DNAH10 and PLXND1. For both genes, we generated adipose tissue (cg-Gal4) and neuronal (elav-Gal4) specific RNAi-knockdown crosses to knockdown transcripts in a tissue specific manner, leveraging upstream activation sequence (UAS)-inducible short-hairpin knockdown lines, available through the VDRC (Vienna Drosophila Resource Center). Specifically, elav-Gal4, which​ drives expression of the RNAi construct in post mitotic neurons starting at embryonic stages all the way to adulthood, was used. Cg drives expression in the fat body and hemocytes starting at embryonic stage 12, all the way to adulthood. (Supplementary Note). Resulting triglyceride values were normalized to fly weight and larval/population density. We used the non-parametric Kruskall-Wallis test to compare wild type with knockdown lines.

Expression quantitative trait loci (eQTLs) analysis

We queried the significant variant (Exome coding SNPs)-gene pairs associated with eGenes across five metabolically relevant tissues (skeletal muscle, subcutaneous adipose, visceral adipose, liver and pancreas) with at least 70 samples in the GTEx database[46]. For each tissue, variants were selected based on the following thresholds: the minor allele was observed in at least 10 samples, and MAF ≥ 1%. eGenes, genes with a significant eQTL, are defined on a false discovery rate (FDR)[130] threshold of ≤0.05 of beta distribution-adjusted empirical p-value from FastQTL. Nominal p-values were generated for each variant-gene pair by testing the alternative hypothesis that the slope of a linear regression model between genotype and expression deviates from 0. To identify all significant variant-gene pairs associated with eGenes, a genome-wide empirical p-value threshold[64] (pt) was defined as the empirical p-value of the gene closest to the 0.05 FDR threshold. pt was then used to calculate a nominal p-value threshold for each gene based on the beta distribution model (from FastQTL) of the minimum p-value distribution f(pmin) obtained from the permutations for the gene. For each gene, variants with a nominal p-value below the gene-level threshold were considered significant and included in the final list of variant-gene pairs[64]. For each eGene, we also listed the most significantly associated variants (eSNP). Only these exome SNPs with r2 > 0.8 with eSNPs were considered for biological interpretation (Supplementary eQTL GTEx).

We also performed cis-eQTL analysis in 770 METSIM subcutaneous adipose tissue samples as described in Civelek, et al.[131] A false discovery rate (FDR) was calculated using all p-values from the cis-eQTL detection in the q-value package in R. Variants associated with nearby genes at an FDR less than 1% were considered to be significant (equivalent p-value < 2.46 × 10−4).

For loci with more than one microarray probeset of the same gene associated with the exome variant, we selected the probeset that provided the strongest LD r2 between the exome variant and the eSNP. In reciprocal conditional analysis, we conditioned on the lead exome variant by including it as a covariate in the cis-eQTL detection and reporting the p-value of the eSNP and vice versa. Signals were considered coincident if both the lead exome variant and the eSNP were no longer significant after conditioning on the other and the variants were in high LD (r2 > 0.80).

For loci that also harbored reported GWAS variants, we performed reciprocal conditional analysis between the GWAS lead variant and the lead eSNP. For loci with more than one reported GWAS variant, the GWAS variant with the strongest LD r2 with the lead eSNP was reported.

Penetrance analysis

Phenotype and genotype data from UKBB were used for penetrance analysis. Three of 16 rare and low frequency variants (MAF ≤ 1%) detected in the final Stage 1+2 meta-analysis were available in UKBB and had relatively larger effect sizes (>0.90). Phenotype data for these three variants were stratified by WHR using World Health Organization (WHO) guidelines, which consider women and men with WHR greater than 0.85 and 0.90 as obese, respectively. Genotype and allele counts were used to calculate the number of carriers of the minor allele. The number of obese vs. non-obese carriers for women, men and sexes combined was compared using a χ2 test. Significance was determined using a Bonferroni correction for the number of tests performed (0.05/9=5.5×10−3).