Protein-altering variants associated with body mass index implicate pathways that control energy intake and expenditure in obesity.

Genome-wide association studies (GWAS) have identified >250 loci for body mass index (BMI), implicating pathways related to neuronal biology. Most GWAS loci represent clusters of common, noncoding variants from which pinpointing causal genes remains challenging. Here we combined data from 718,734 individuals to discover rare and low-frequency (minor allele frequency (MAF) < 5%) coding variants associated with BMI. We identified 14 coding variants in 13 genes, of which 8 variants were in genes (ZBTB7B, ACHE, RAPGEF3, RAB21, ZFHX3, ENTPD6, ZFR2 and ZNF169) newly implicated in human obesity, 2 variants were in genes (MC4R and KSR2) previously observed to be mutated in extreme obesity and 2 variants were in GIPR. The effect sizes of rare variants are ~10 times larger than those of common variants, with the largest effect observed in carriers of an MC4R mutation introducing a stop codon (p.Tyr35Ter, MAF = 0.01%), who weighed ~7 kg more than non-carriers. Pathway analyses based on the variants associated with BMI confirm enrichment of neuronal genes and provide new evidence for adipocyte and energy expenditure biology, widening the potential of genetically supported therapeutic targets in obesity.

Obesity is a heritable disease and represents a major unmet public health problem with only a few safe and long-term effective therapies[1] and intervention strategies[2]. To understand the genetic basis of obesity and identify potential targets for new therapies, genome-wide association studies (GWAS) for body mass index (BMI) and obesity risk have identified >250 common variants over the past decade[3-7]. Consistent with single-gene disorders of obesity[8], tissue expression and gene-set enrichment analyses for genes in BMI-associated loci have shown that the central nervous system (CNS) plays a critical role in body weight regulation[5]. While the numerous GWAS loci have provided insight into broad biological mechanisms underlying body weight regulation, pinpointing the causal gene(s)/variant(s) remains a major challenge[9], as GWAS-identified variants are typically non-coding and may affect genes at long distance. The association of intronic FTO variants with BMI illustrates the challenges of identifying causal regulatory effects. The proposed causal variant was found to regulate the expression of nearby RPGRIP1L in some studies[10-12], whereas others found that it regulates distant IRX3/IRX5 genes in specific cell types[13,14].

To expedite mapping of obesity-related genes, we performed an exome-wide search for low-frequency (LF, minor allele frequency [MAF]=1–5%) and rare (R, MAF<1%) single nucleotide variants (SNVs) associated with BMI using exome-targeted genotyping arrays. A total of 125 studies (Nindividuals=718,734) performed single-variant association between up to 246,328 SNVs and BMI. In addition, we performed gene-based meta-analyses to aggregate rare and LF (R/LF) coding SNVs across 14,541 genes. Using genetic, functional and computational follow-up analyses, we gained insights into the function of BMI-implicated genes, and the biological pathways through which they may influence body weight.

RESULTS

Fourteen rare and low-frequency coding variants in 13 genes

Our study comprises a discovery and a follow-up stage (Supplementary Figure 1, Supplementary Tables 1–3, Online Methods). In our primary analysis, the discovery stage includes data from 123 studies (N=526,508) across five ancestry groups, predominantly European (~85%). Each study performed single-variant association analyses of coding variants present on the exome array, including up to 13,786 common (MAF>5%) and 215,917 R/LF coding SNVs (exons and splicing sites). Summary statistics were combined using fixed-effect meta-analyses. SNV-associations of R/LF variants that reached suggestive significance (P<2.0×10−6) were taken forward for follow-up in two European cohorts, deCODE (N=72,613) and UK Biobank (N=119,613 [interim release]). Overall significance was assessed after combining results of discovery and follow-up studies into a final meta-analysis (all-ancestries, sex-combined, additive model, Nmax=718,734); SNV-associations that reached P<2×10−7 were considered array-wide significant[15,16] (Table 1, Supplementary Table 4, Supplementary Figures 2–4). In secondary analyses, we performed sex-specific analyses, analyses limited to individuals of European ancestry, and analyses using a recessive model.

Table 1

Rare and low-frequency coding variants significantly associated with BMI

Chr:position	Variant	Codinglocus	Allele		Amino acidchange	EAF(%)	β(SD/allele)	SE	P-value	N	Explainedvariance(%)
			Effect	Other
All-ancestries additive
1:154987704	rs141845046	ZBTB7B*	T	C	p.Pro190Ser	2.44%	0.048	0.006	7.73E-18	718,628	0.011%
7:100490797	rs1799805	ACHE*	T	G	p.His353Asn	3.90%	0.029	0.005	2.82E-10	707,448	0.006%
12:48143315	rs145878042	RAPGEF3*	G	A	p.Leu300Pro	1.10%	0.066	0.008	1.56E-15	700,852	0.010%
12:49399132	rs1126930	PRKAG1	C	G	p.Thr38Ser	3.22%	0.034	0.005	3.98E-12	712,354	0.007%
12:72179446	rs61754230	RAB21*	T	C	p.Ser224Phe	1.74%	0.040	0.007	1.33E-09	693,373	0.005%
12:117977550	rs56214831	KSR2	T	C	p.Arg525Gln	0.82%	0.057	0.010	1.08E-08	655,049	0.005%
12:123345509	rs34149579	HIP1R	T	G	p.Cys938Phe	4.54%	−0.032	0.004	2.00E-14	716,253	0.009%
16:72830539	rs62051555	ZFHX3*	G	C	p.Gln1100His	4.34%	−0.024	0.004	4.01E-08	690,637	0.005%
18:58039478	rs13447324	MC4R	T	G	p.Tyr35Ter	0.01%	0.542	0.086	2.26E-10	631,683	0.006%
19:46178020	rs139215588	GIPR	A	G	p.Arg190Gln	0.11%	−0.148	0.028	1.25E-07	695,800	0.005%
19:46180976	rs143430880	GIPR	G	A	p.Glu288Gly	0.13%	−0.153	0.028	2.96E-08	599,574	0.006%
20:25195509	rs6050446	ENTPD6*	A	G	p.Lys185Glu	2.71%	−0.034	0.005	2.40E-10	717,084	0.006%
All-ancestries sex-specific additive (women only)
19:3813906	rs45465594	ZFR2*	C	A	p.Ile718Met	2.55%	−0.040	0.008	1.94E-07	373,848	0.008%
European Ancestry additive
9:97062981	rs12236219	ZNF169*	T	C	p.Arg381Cys	4.23%	−0.029	0.005	8.78E-10	612,396	0.007%

Array-wide significant is defined as P < 2×10−7.

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

Alleles (effect/other), effect allele frequency (EAF), beta (b), standard error (SE) and P values are based on the meta-analysis of Discovery Stage (GIANT) and Validations stage (deCODE, UKBiobank) studies. Effect allele is always the minor allele. Effects (b) are expressed in SD, assuming mean=0 and SD=1.

The amino acid change from the most abundant coding transcript is shown in this table (see Supplementary Table 25 for more details on protein annotation based on VEP tool and transcript abundance from GTEx database).

Novel gene; i.e. not previously implicated in human obesity

In our primary analysis of R/LF variants, we identified five rare SNVs in three genes (KSR2, 2 in MC4R, 2 in GIPR) and nine LF SNVs in eight genes (ZBTB7B, 2 in ACHE, RAPGEF3, PRKAG1, RAB21, HIP1R, ZFHX3, ENTPD6) (Table 1, Box 1, Supplementary Table 5, Supplementary Figure 3a). In secondary analyses, we identified two additional LF SNVs; one in all-ancestry women-only (ZFR2) and one in European ancestry only analyses (ZNF169) (Table 1, Supplementary Tables 6–8, Supplementary Figures 3b, 3c). Of these 16 SNVs, located in 13 genes, the two SNVs in MC4R (r=1; D’=1) and two in ACHE (r=0.98; D’=0.99) were in high LD, whereas the two SNVs in GIPR (r=0; D’=0.16) were independent of each other. Hence, the 16 SNVs represent 14 independent SNVs (4 rare, 10 LF), of which eight locate in genes not previously implicated in BMI (ZBTB7B, ACHE, RAPGEF3, RAB21, ZFHX3, ENTPD6, ZFR2, ZNF169), and six are located in five loci that were previously identified by GWAS (PRKAG1/BCDIN3D, HIP1R/CLIP1, MC4R, GIPR/QPCTL)[5] and/or through sequencing of severe early-onset obesity cases (MC4R, KSR2)[17-19] (Supplementary Figure 5). Conditional analyses established that coding SNVs in PRKAG1, MC4R and GIPR are independent of the common lead variants in GWAS loci (rs7138803, rs17782313, rs2287019, respectively), whereas the SNV in HIP1R and GWAS locus near CLIP1 (rs11057405) represent the same signal (Online Methods, Supplementary Tables 9, 10, Supplementary Figure 5).

Next, we performed gene-based association tests (SKAT, VT; broad, strict) in up to 14,541 genes[20] to examine whether these aggregated analyses would yield new evidence for multiple R/LF coding SNVs in the same gene affecting BMI (Online Methods). Using broad SNV inclusion criteria, associations for 13 genes reached array-wide significance (P<2.5×10−6)[15,16], four of which had not been highlighted in single-variant analyses (Table 2, Supplementary Table 11). Conditional analyses showed that only for GIPR was the gene-based association driven by multiple SNVs (Table 2, Supplementary Table 12). For all other genes, associations were driven by a single SNV only, but these SNVs had not reached array-wide significance in single-variant analyses.

Table 2

Genes significantly associated with BMI in a gene-based meta-analyses, aggregating R/LF coding SNVs

Gene	Location longest codingtranscript	Testd	Nvariants	P-value	Conditioned P-valuea	Single variant
						Top variant	P-value
All-ancestries sex-combined
SLC6A17	chr1:110693132–110744823	SKAT	13	2.73E-07	0.13	rs41313405	4.45E-07
RAPGEF3	chr12:48128453–48152889	SKAT	19	8.91E-15	0.20	rs145878042	5.16E-14
PRKAG1	chr12:49396055–49412629	SKAT	4	2.75E-12	0.53	rs1126930	2.63E-12
RAB21	chr12:72148643–72187256	SKAT	5	4.81E-08	0.27	rs61754230	4.96E-08
KSR2	chr12:117890817–118406028	SKAT	7	7.15E-09	0.19	rs56214831	4.59E-08
MAP1A	chr15:43809806–43823818	SKAT	25	9.42E-07	0.16	rs55707100	1.01E-06
MC4R	chr18:58038564–58040001	VT	4	3.72E-09	0.01	rs13447325	2.97E-11
GIPR	chr19:46171502–46186982	VT	10	8.24E-09	1.12E-04	rs143430880	5.76E-06
All-ancestries sex-specific
ALDH3A1 (men only)	chr17:19641298–19651746	SKAT	15	3.24E-07	0.003	rs142078447	8.62E-06
ZFR2 (women only)	chr19:3804022–3869027	SKAT	19	1.81E-07	0.82	rs45465594	3.64E-07
European sex-combined
ACHE	chr7:100487615–100493592	SKAT	6	3.30E-10	0.12	rs386545548	7.22E-10
European sex-specific
ANGPTL7 (men only)	chr1:11249346–11256038	VT	3	2.50E-06	0.008	rs202182115	2.56E-05
ZNF169 (women only)	chr9:97021548–97064111	SKAT	9	1.89E-07	0.24	rs12236219	1.06E-06

Array-wide significant gene-based association is defined as P<2.5×10−6. P-values are based on the meta-analysis of Discovery Stage studies.

Gene-based analyses were performed with SKAT and VT; results shown are from the test (SKAT or VT) for which the significance exceeded P<2.5×10−6. Only results using the "broad" SNV inclusion criteria reached array-wide significance.

Transcript positions are reported according to Build 37 for the longest coding transcript supported by RefSeq (as displayed in USCS Genome Browser).

P-value after conditioning on the most significant (top) single variant aggregated in the gene-based test.

Taken together, we identified 14 R/LF coding SNVs in 13 genes that are independently associated with BMI; four rare SNVs in three genes, and 10 LF SNVs in 10 genes. One SNV (ZFR2) showed a sex-specific effect, whereas no ancestry-specific effects were observed (Supplementary Note, Supplementary Tables 6–8, Supplementary Figure 6). Eight (ACHE, ENTPD6, RAB21, RAPGEF3, ZBTB7B, ZFHX3, ZFR2, ZNF169) of these 13 genes have not been previously implicated in body weight regulation (Table 1).

Novel common coding variants associated with BMI

Although the main focus of our study was on R/LF coding SNVs, we also identified 92 common coding variants (P<2.0×10−7; Supplementary Tables 4; Supplementary Figures 4, 7), of which 41 were novel (Supplementary Table 9, Supplementary Note). These novel common loci had not been identified in previous GWAS efforts, because our current sample size is more than twice as large as the most recent GWAS meta-analysis[5], and also because some SNVs were not tested before, as they were not present on the HapMap reference panel and/or were on the X-chromosome, which was not analyzed. Because of the increased samples size, effect sizes of the 41 novel common loci are smaller (on average 0.014 SD/allele, [range: 0.010–0.024]) than of previously established common loci (0.021 SD/allele, [0.010–0.050]) (Supplementary Figure 7).

Impact of R/LF SNVs on BMI and obesity risk

The minor allele for half of the 14 R/LF SNVs is associated with lower BMI (Table 1, Figure 1). The effects of LF SNVs range between 0.024 and 0.066 SD/allele, equivalent to ~0.11 to 0.30 kg/m2 in BMI or ~0.315 to 0.864 kg in body weight for a 1.7m tall person. Effects of rare SNVs range between 0.06 and 0.54 SD per allele, equivalent to 0.26 to 2.44 kg/m2 or 0.74 kg to 7.05 kg per allele (Table 1, Figure 1). By comparison, these rare SNV effect sizes are on average ten times larger than those for previously identified GWAS loci (effectmean=0.019 SD/allele, ~0.086 kg/m2 or ~0.247 kg/allele) of which the largest effect is seen for the FTO locus (0.08 SD/allele, ~0.35 kg/m2 or 1 kg/allele) and those for other GWAS loci range between 0.010 and 0.056 SD/allele (~0.045 to 0.25 kg/m2, or 0.130 to 0.728 kg)[5].

Figure 1

Effect sizes (y-axis) of the 14 BMI-associated R/LF coding variants by their minor allele frequency

Effect sizes are expressed in body weight (kg) per allele, assuming a SD of 4.5 kg and an average-sized person of 1.7m tall. Solid markers indicate that the minor allele is associated with higher BMI, and clear markers indicate that the minor allele is associated with lower BMI. Variants were identified in all-ancestry analyses (light blue diamonds), the European ancestry analyses (dark blue square) and women-only analyses (pink diamond). Effect sizes for previously identified GWAS loci are shown in navy blue diamonds. The dotted line represents 80% power, assuming α = 2×10−7 and N= 525,000 (discovery sample size).

Effect sizes increase as MAF decreases, in particular for SNVs with a MAF<0.5% (~1 heterozygote carrier in 100 people), consistent with the statistical power of our sample (Figure 1). For example, the nonsense p.Tyr35Ter MC4R SNV (rs13447324, MAF=0.01%) is present in ~1 in 5,000 individuals and results in a ~7 kg higher body weight for a 1.7m tall person. The two GIPR SNVs contribute independently to a lower body weight; carriers (1 in ~455 individuals) of p.Arg190Gln (rs139215588) weigh ~1.92 kg (0.148 SD BMI) less than non-carriers and carriers (1 in ~385 individuals) of p.Glu288Gly (rs143430880) weigh ~1.99 kg (0.153 SD BMI) less. Among 115,611 individuals of the UK Biobank, one apparently healthy 61-year-old woman, with no reported illnesses, carried both rare GIPR alleles and weighed ~11.2 kg less (equivalent to −0.86 SD BMI or 3.87 kg/m2) than the average non-carrier of the same height (Supplementary Figure 8). The possible synergistic effect of the two GIPR alleles needs confirmation by additional individuals that carry both variants.

Even though effect sizes of LF and, in particular, rare SNVs tend to be larger than those of common GWAS-identified loci[5], the 14 SNVs combined explain <0.1% of BMI variation, because of their low population frequency (Table 1, Online Methods). Also, although the effects of the four rare SNVs (KSR2, MC4R, 2 in GIPR) are large by GWAS standards, penetrance for obesity is still expected to be low. Indeed, using data from the UK Biobank (Nmax=119,781), we compared the prevalence of normal-weight (18.5 kg/m2 ≤ BMI < 25 kg/m2) and obesity (BMI ≥ 30 kg/m2) between carriers and non-carriers (Supplementary Table 13, Online Methods). For GIPR (p.Arg190Gln, p.Glu288Gly), both BMI-decreasing SNVs, carriers tended (P<0.05) to have a lower obesity prevalence (21.2%, 20.1%, respectively), compared to non-carriers (25.1%, 25%). For MC4R p.Tyr35Ter and KSR2 p.Arg525Gln, the prevalence of obesity between carriers (30%, 25.7%, resp.) and non-carriers (25.1%, 25.3%) was not significantly different.

We examined whether R/LF SNVs affect obesity risk early on in life by combining data from three case-controls studies of childhood obesity (N=4,395; N=13,072) (Online Methods, Supplementary Table 14). Associations for 10 of 13 SNVs were directionally consistent with those observed for BMI in adults (77%, P=0.046), three of which (ZBTB7B, PRKAG1, RAB21) reached nominal significance (P<0.05). While no carriers of the MC4R mutations were available for analyses, the role of MC4R in body weight regulation in childhood was established almost two decades ago[17,19,21].

Impact of R/LF SNVs on cardiometabolic and other traits

To examine whether identified SNVs affect other traits, we obtained results from multiple large-scale genetic consortia (GIANT[15], MAGIC, GoT2D/T2D-GENES[16], GLGC, ICBP[22], REPROGEN[23]) (Supplementary Table 15, Supplementary Figure 9), and performed phenome-wide association (PheWAS) analyses using electronic medical record (EMR) data from BioVu and UK Biobank (Online Methods, Supplementary Table 16). The BMI-increasing allele of ZBTB7B p.Pro190Ser is associated with greater height, and those of PRKAG1, ACHE, and RAPGEF3 SNVs are associated with shorter height, but association with other traits differ. Specifically, PRKAG1 p.Thr38Ser Ser-allele carriers appear heavier and shorter, have lower HDL-cholesterol levels, earlier age at menarche (reported before[23]) and higher systolic blood pressure, which is in agreement with PheWAS analyses showing an increased risk of “malignant essential hypertension” and “hypertension” (Supplementary Table 16). While carriers of the RAPGEF3 p.Leu300Pro Pro-allele are also heavier and shorter, they have a lower WHRadjBMI[24] and lower fasting insulin levels (Supplementary Table 15), consistent with PheWAS results that show lower odds of “secondary diabetes mellitus” (Supplementary Table 16). Thus, while all SNVs are associated with BMI, their patterns of association with other traits suggest they may affect different physiological pathways.

Gene set enrichment analyses

To test whether the R/LF variants implicate biological pathways, we performed gene set enrichment analyses. Similar to our previous analysis of GWAS for BMI[5], we analyzed coding variants that reached P<5×10−4, using a DEPICT version adapted for exome-array analysis[15] (Online Methods, Supplementary Note). We used 50 R/LF coding variants as input (all P<5×10−4; Online Methods) and observed significant enrichment (Figure 2, Supplementary Table 17, Supplementary Figure 10a). Many of these relate to neuronal processes, such as neurotransmitter release and synaptic function (e.g. glutamate receptor activity, regulation of neurotransmitter levels, synapse part), consistent with previous findings from GWAS[5]. When we excluded variants near (+/− 1Mb) previously identified GWAS loci, we still observed 29 significantly enriched gene sets (in 12 meta-gene sets) (Supplementary Table 18, Supplementary Figure 10b), thereby providing an independent confirmation of the GWAS gene set enrichment results. In addition to neuronal-related gene sets, the analyses with R/LF coding variants newly identified a cluster of metabolic pathways related to insulin action and adipocyte/lipid metabolism (e.g. enhanced lipolysis, abnormal lipid homeostasis, increased circulating insulin level; Figure 2). Finally, we observed that R/LF BMI-associated coding variants are more effective at identifying enriched gene sets compared to common coding variants. Specifically, adding 192 common coding SNVs (all P<5×10−4) to the analysis decreased the number of enriched gene sets from 471 (106 meta-gene sets) seen with R/LF coding SNVs to 62 (24 meta-gene sets) (Supplementary Table 19, Supplementary Figure 10c). We observed fewer significant genes sets with the combined common and R/LF analysis, despite including more total coding variants and a higher fraction of array-wide significant coding variants. One possible explanation is that R/LF coding variants may fall in the causal gene more often than do common coding variants, which suggests that the R/LF variants are more likely to be causal, rather than simply in LD with causal variants.

Figure 2

Heatmap showing DEPICT gene set enrichment results for suggestive and significant rare and low-frequency coding SNVs

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), based on the gene’s Z-score for gene set inclusion in DEPICT’s reconstituted gene sets (red indicates a higher, 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 Table 17) correspond to the most significantly enriched gene set within the meta-gene set. Annotations for genes indicate (1) whether it has an OMIM annotation as underlying a monogenic obesity disorder (black/grey), (2) the MAF of the significant ExomeChip (EC) variant (blue), (3) whether the variant’s P-value reached array-wide significance (<2×10−7) or suggestive significance (<5×10−4) (purple), (4) whether the variant was novel, overlapping “relaxed” GWAS signals from Locke et al.[5] (GWAS P<5×10−4), or overlapping “stringent” GWAS hits (GWAS P<5×10−8) (pink), and (5) whether the gene was included in the gene set enrichment analysis or excluded by filters (orange/brown) (Online Methods, Supplementary Note). Annotations for gene sets indicate if the meta-gene set was significant (green; FDR <0.01, <0.05, or not significant) in the DEPICT analysis of GWAS results[5]. Here, two regions of particularly strong gene set membership are shown (see full heat map in Supplementary Figure 10a).

We also used gene set enrichment analysis to prioritize candidate genes. Among the genes with R/LF coding variants associated with BMI at P<5×10−4, a subset is prominently represented in the CNS-related enriched gene sets (Figure 2) and is proposed to influence neurotransmission and/or synaptic organization, function and plasticity. These include genes in regions with suggestive evidence of association from GWAS (e.g. CARTPT, MAP1A, ERC2) and genes in regions not previously implicated by GWAS (e.g. CALY, ACHE, PTPRD, GRIN2A). The non-neuronal metabolic gene sets implicate two genes (CIDEA, ADH1B) that are markers of brown or “beige” adipose tissue[25,26], providing new supporting evidence for a causal role of this aspect of adipocyte biology.

Drosophila fly results

To test for potential adiposity-driving effects of gene regulation, we performed tissue-specific RNAi-knockdown experiments in Drosophila. We generated adipose-tissue (cg-Gal4) and neuronal (elav-Gal4) specific RNAi-knockdown crosses for nine of the 13 candidate genes for which fly orthologues exist (Supplementary Table 20) and performed whole body triglyceride analysis in young adult male flies. Triglycerides, the major lipid storage form in animals, were chosen as a direct measure of fly adiposity. Both neuronal and fat-body knockdown of zfh2, the orthologue of ZFHX3, resulted in significantly increased triglyceride levels. Adipose-tissue specific, but not neuronal, knockdown of epac (RAPGEF3) was lethal. Tissue-specific loss-of-function of the other seven genes tested did not affect triglyceride levels.

R/LF coding SNVs in monogenic and syndromic genes

We identified 39 genes in the literature that have been convincingly implicated in monogenic obesity or syndromes of which obesity is one of the main features (Supplementary Table 21, 22, Supplementary Figure 11). Of the 652 R/LF SNVs in these 39 monogenic and/or syndromic genes, five R/LF SNVs were significantly associated with BMI (Bonferroni-corrected P-value = 7.7×10−5 (=0.05/652)). Beside SNVs in MC4R (p.Tyr35Ter, Asp37Val) and KSR2 (Arg525Gln), already highlighted in the single-variant analyses, we identified an additional SNV in MC4R (p.Ile251Leu) and one in BDNF (p.Glu6Lys). MC4R p.Ile251Leu has been previously shown to protect against obesity[27], whereas BDNF p.Glu6Lys, independent of previously GWAS-identified SNVs (r2=0.01, D’=1.0)[5], has not been implicated in body weight regulation before. We examined whether the 652 R/LF SNVs showed enrichment for association with BMI compared to R/LF coding SNVs in all other genes, but found no evidence to support this.

DISCUSSION

In this meta-analysis of exome-targeted genotyping data, we identified 14 R/LF coding variants in 13 genes associated with BMI. Eight of these genes (ACHE, ENTPD6, RAB21, RAPGEF3, ZBTB7B, ZFHX3, ZFR2, ZNF169) have not been previously implicated in human obesity, but evidence from animal studies provides support for a role in energy metabolism for some of these, such as ACHE[28,29], RAPGEF3[30-33], and PRKAG1[34-39]. Others fall into established BMI GWAS loci (PRKAG1/BCDIN3D, HIP1R/CLIP1, MC4R, GIPR/QPCTL)[5] and/or were previously implicated in severe early-onset obesity (MC4R, KSR2)[17-19] and using this exome-targeted approach, we pinpoint R/LF variants in these loci that play a role in obesity in the general population. Pathway analyses confirm a key role for neuronal processes, and newly implicate adipocyte and energy expenditure biology.

Consistent with other polygenic traits[15,23,40-43], we show that large sample sizes are needed to identify R/LF variants. Observed effect sizes reflect the statistical power of our sample size, and are particularly large for SNVs with a MAF < 0.05%. The existence of rare alleles with larger effects on BMI than have been observed for common alleles might reflect negative or stabilizing selection on the extremes of BMI. However, rare variants with smaller effects almost certainly exist; larger samples will be needed to uncover these. Our study was limited to coding variants on the exome-array; large-scale sequencing studies will be needed to test for variants not covered by exome-arrays.

The strongest association was observed for a stop-codon (p.Tyr35Ter, rs13447324, MAF= 0.01%) in MC4R, with carriers weighing on average 7kg more than non-carriers. MC4R is widely expressed in the CNS and is an established key player in energy balance regulation[44,45]. Mouse and human studies showed already two decades ago that MC4R-deficiency results in extreme obesity, mainly through increased food intake[46-49]. p.Tyr35Ter, which results in MC4R-deficiency[51], was one of the first MC4R mutations discovered in monogenic cases of obesity[17,19], in whom the mutation is >20× more prevalent than in the general population[17,50,52,53]. Here, we show that p.Tyr35Ter plays a role outside the setting of early-onset and extreme obesity. Despite its large effect, penetrance is low, and does not fit the model of a fully penetrant Mendelian variant.

While significant R/LF coding variants are strong candidates for being causal, the strongest implication of causal genes is provided by association with multiple independent coding variants, as we demonstrate for GIPR. We identified two rare variants in GIPR (p.Arg190Gln, rs139215588, MAF=0.11%; p.Glu288Gly, rs143430880, MAF=0.13%) independently associated with lower BMI; carriers of either variant weigh ~2 kg less than non-carriers. Common variants in/near GIPR have been found to associate with lower BMI[55] and delayed glucose and insulin response to an oral glucose challenge[54]. However, the two rare variants influence BMI independently of these common ones and are not associated with type 2 diabetes or glycemic traits tested. Rodent models have provided strong evidence for a role of GIPR in body weight regulation. Gipr-deficient mice are protected from diet-induced obesity[56] and have an increased resting metabolic rate[57]. Blocking GIP-signaling using a vaccination approach in mice on a high-fat diet reduces weight gain, mainly through reduced fat accumulation, mediated through increased energy expenditure[58]. Manipulation of incretins (GIP, GLP1) and their receptors has complex effects on obesity and insulin secretion/action that may differ between human and mice[59]. The human genetic data suggest that inhibition of GIPR-signaling might present a therapeutic target for the treatment of obesity[60].

A fourth rare variant, in KSR2, (p.Arg525Gln, rs56214831, MAF=0.82%) increases body weight by ~740g/allele. KSR2 is another gene previously implicated in energy metabolism and obesity[18,61,62]. In a recent study, mutation carriers were hyperphagic, had a reduced basal metabolic rate and severe insulin resistance[18]. Consistent with human data, Ksr2 mice were obese, hyperphagic, and had a reduced energy expenditure[18,61-63]. KSR2 is almost exclusively expressed in the brain and interacts with multiple proteins[64], including AMP-activated protein kinase (AMPK), a key regulator of energy homeostasis[61,62]. Interestingly, KSR2 is one of the first genes implicated in severe, early-onset obesity in which mutations not only affect food intake but also basal metabolic rate, and is thought to act via neuronal effects[18] (Figure 2).

Despite convincing associations of these four rare variants in MC4R, GIPR and KSR2, their penetrance for obesity is low (Supplementary Table 13). This is consistent with the polygenic and multifactorial nature of obesity, where variants across a range of frequencies and effect sizes contribute to the phenotype in any one person. Despite low predictive power, it remains possible that the identities of particular variants in any one person may contribute to different balances of underlying physiologies and hence, different responses to treatments. This was illustrated in two patients with monogenic obesity due to POMC mutations; these patients lack the main activator of MC4R and were effectively treated with an MC4R-agonist[65].

Of the coding variants in newly identified genes, some have well-known connections to obesity. For example, PRKAG1 encodes the γ1-subunit of AMPK, a critical cellular energy sensor[34]. In the hypothalamus, AMPK integrates hormonal and nutritional signals with neuronal networks to regulate food intake and whole-body energy metabolism[35-37]. Furthermore, hypothalamic AMPK is a key regulator of brown adipose tissue in mice[36,38,39]. The BMI-decreasing allele at the associated PRKAG1 variant (p.Thr38Ser, rs1126930, MAF=3.22%) has additional beneficial effects on blood pressure, providing additional genetic support for modulation of AMPK as an ongoing therapeutic avenue for treatment.

ACHE, in which p.His353Asn (rs1799805, MAF=3.9%) is associated with increased BMI, is another candidate gene related to neuronal biology, involved in the signaling of acetylcholine at neuromuscular junction and brain cholinergic synapses[67,68]. Inhibitors of ACHE, used to treat moderate-to-severe Alzheimer’s Disease[69], results in weight loss in humans and Ache-deficient mice have delayed weight gain[28,29]. However, these may be indirect consequences of adverse gastrointestinal and neuromuscular effects, respectively[28,29,70,71].

Another LF coding variant (p.Leu300Pro, rs145878042, MAF=1.1%) is located in RAPGEF3, and has strong effects on multiple other phenotypes. The BMI-increasing 300Pro-allele is associated with shorter height, lower WHRadjBMI and lower insulin levels, suggesting that this variant has multiple physiologic consequences. Data from animal models also suggest complex effects of RAPGEF3 on adipocyte biology, energy balance and glucose metabolism[30-33]. For example, in one study, global deletion of Rapgef3 in mice on a high-fat diet are resistant to obesity due to reduced food intake and have an increased glucose tolerance[31]. However, in a similar study, Rapgef3 mice develop severe obesity, increased respiratory exchange ratio and impaired glucose tolerance[33]. Adipose tissue-specific Rapgef3 knockout mice on a high-fat diet are also more prone to obesity, show increased food intake, reduced energy expenditure, impaired glucose tolerance, and reduced circulating leptin levels[72]. More research is needed to understand the consequences of RAPGEF3 manipulation.

The remaining genes with significant associations, ENTPD6, HIP1R, RAB21, ZFR2, ZBTB7, and ZFHX3, have no clear prior evidence for a role in energy homeostasis, and in-depth functional follow up is needed to gain insight in how they affect body weight. Here, we performed gene set enrichment analyses to better understand the biology implicated by our genetic data, and confirm the importance of neuronal processes, in particular synaptic function and neurotransmitter release, providing an independent validation of previous GWAS findings[5]. The combination of gene set enrichment and association analyses of coding variants also enables us to highlight candidate genes that are both within these gene sets and show association with BMI at R/LF coding variants. These include genes reaching array-wide significance (e.g. ACHE, ZFR2), and others with clear prior evidence for a role in body weight regulation (e.g. CARTPT[73]), but that had not been highlighted in our single-variant or gene-based association analyses. Of note, the enrichment signals were stronger with R/LF coding variants only than with all coding variants, suggesting that R/LF variants are more likely to be causal and may more often point directly to relevant genes, whereas common coding variants may more often be proxies for common noncoding variants that affect nearby genes.

In addition, our gene set enrichment analyses now provide supporting evidence for a role of non-neuronal mechanisms as well. Specifically, CIDEA and ADH1B are both strongly predicted to be members of enriched gene sets related to insulin action and adipocyte biology, and both are markers that distinguish brown from white fat depots in mice[25] and humans[26]. CIDEA is predominantly expressed in adipose tissue and known as a key regulator of energy metabolism[25]. Cidea-deficient mice are resistance to diet-induced obesity with increased lipolysis and mitochondrial uncoupling[25]. The connection of ADH1B to obesity is less clear, but the gene is highly expressed in human adipocytes, has been implicated by gene expression analyses in obesity and insulin resistance, and functions early in a potentially relevant metabolic pathway (retinoid biosynthesis)[25,26,74,75]. Similar pathways were implicated by recent work dissecting the signal near FTO[13]. However, because SNV-association signals at ADH1B and CIDEA did not reache array-wide significance, additional genetic analysis of their role in obesity would be warranted.

In summary, we performed association analyses between R/LF variants and BMI in >700,000 individuals, and identified 14 variants in 13 genes, in 5 known and 8 novel genes. While each variant contributes little to BMI variation in the general population, they may have substantial impact on body weight at an individual level. Furthermore, prior literature for these genes and unbiased gene set enrichment analysis indicate a strong role for neuronal biology and also provide new support for a causal role of aspects of adipocyte biology. The identified genes provide potential targets that may lead to new and more precise approaches for the treatment of obesity, which has seen minimal innovation in the past 30 years[1].

ONLINE METHODS

Study design & participants

The discovery cohort consisted of 123 studies (163 datasets) comprising 526,508 adult (≥18yrs) individuals of the following ancestries (Supplementary Figure 1): 1) European (N = 449,889), 2) South Asian (N = 29,398), 3) African (N = 27,610), 4) East Asian (N = 8,839), and 5) Hispanic (N = 10,772). All participating institutions and coordinating centers approved this project and informed consent was obtained from all study participants. Discovery meta-analyses were carried out in each ancestry separately and in the All-ancestries combined group, for both sex-specific and sex-combined analyses. SNVs for which associations reach suggestive significance (P<2.0×10−6) in the discovery analyses, were taken forward for follow-up in 192,226 individuals of European ancestry from the UK BioBank and deCODE. Conditional analyses were conducted in the All-ancestries and European descent groups. Study-specific design, sample quality control and descriptive statistics are provided in Supplementary Tables 1–3.

Phenotype

Body mass index (BMI: weight [in kilograms] / height [in meters][2]) was corrected for age, age[2] and genomic principal components (PC, derived from GWAS data, the variants with MAF > 1% on ExomeChip, or 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 non-related individuals, residuals were calculated separately by sex, whereas for family-based studies sex was included as a covariate in the model. Additionally, residuals for case/control studies were calculated separately. Finally, residuals were subject to inverse normal transformation[96].

Genotype calling

The majority of studies followed a standardized protocol and performed genotype calling using the designated manufacturer software, which was then followed by zCall[97]. For 10 studies, participating in the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium, the raw intensity data for the samples from seven genotyping centers were assembled into a single project for joint calling[98]. Study-specific quality control (QC) measures of the genotyped variants were implemented before association analysis (Supplementary Table 2).

Statistical analyses

Study-level association analyses

Individual cohorts were analyzed separately for each ancestry, in sex-combined and sex-specific groups, with either RAREMETALWORKER (see URL links at the end of the Online Methods) or RVTEST[99] (Supplementary Table 2), to associate inverse normal transformed BMI with genotype accounting for potential cryptic relatedness (kinship matrix) in a linear mixed model. These software tools are designed to perform score-statistics based rare-variant association analyses, can accommodate both unrelated and related individuals, and provide single-variant results and variance-covariance matrices. The covariance matrix captures linkage disequilibrium (LD) relationships between markers within 1 Mb, which is used for gene-level meta-analyses and conditional analyses[100]. Single-variant analyses were performed for both additive and recessive models.

Centralized quality-control

A centralized quality-control procedure, implemented in EasyQC[101], was applied to individual cohort association summary statistics to identify cohort-specific problems: (1) assessment of possible problems in BMI transformation, (2) comparison of allele frequency alignment against 1000 Genomes Project phase 1 reference data to pinpoint any potential strand issues, and (3) examination of quantile-quantile (QQ) plots per study to identify any problems arising from population stratification, cryptic relatedness and genotype biases.

Meta-analyses

Meta-analyses were carried out by two different analysts at two sites in parallel. We excluded variants with a call rate < 95%, Hardy-Weinberg equilibrium P-value < 1×10−7, or large allele frequency deviations from reference populations (> 0.6 for all-ancestry analyses and > 0.3 for ancestry-specific population analyses). Significance for single-variant analyses was defined at the array-wide level (a Bonferroni-corrected threshold of P < 2×10−7 for ~250,000 SNVs). To test for sex-differences of the significant variants (P < 2×10−7), we calculated the P-diff for each SNP, which tests for differences between women-specific and men-specific beta estimates using EasyStrata[102]. For gene-based analyses, we applied the sequence kernel association test (SKAT)[103] and the Variable Threshold (VT)[104] gene-based methods using two different sets of criteria (broad and strict) to select predicted damaging R/LF variants with MAF < 5%, based on coding variant annotation from five prediction algorithms (PolyPhen2 HumDiv and HumVar, LRT, MutationTaster and SIFT)[20]. Our broad gene-based tests included nonsense, stop-loss, splice site, and missense variants that are 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. Statistical significance for gene-based tests was set at a Bonferroni-corrected threshold of P < 2.5×10−6 for about 20,000 genes[16,105]. Singe-variant and gene-based meta-analyses were both performed using RareMETALS R-package[106]. As our secondary analyses are nested and/or highly correlated with our primary analysis, we chose the same, already stringent, Bonferroni-corrected significance threshold for both analyses.

Genomic inflation

Although the overall λGC value is in the normal range for all coding variants (λGC = 1.1, Supplementary Table 23), we observed a marked genomic inflation of the test statistics even after adequate control for population stratification (linear mixed model) arising from common markers (λGC = 1.99, Supplementary Figure 2a and Supplementary Table 23). Such inflation is expected for a highly polygenic trait like BMI, as was previously confirmed for height[15], and is consistent with our very large sample size[5,107]. Furthermore, some of the inflation may be due to the design of the ExomeChip, which besides R/LF coding SNVs also contains (common and non-coding) SNVs that include previously identified GWAS loci for all traits, including for BMI and BMI-related traits, reported in the GWAS catalogue at the time of its design.

After removing established loci (+/− 1Mb), the excess of significant associations is markedly reduced and inflation reduced (Supplementary Figures 2c and 2d).

Furthermore, to exclude the possibility that some of the observed associations between BMI and R/LF SNVs could be due to allele calling problems in the smaller studies, we performed a sensitivity meta-analysis with primarily European ancestry studies totaling >5,000 participants. We found very concordant effect sizes, suggesting that smaller studies do not bias our results (Supplementary Figure 12).

Follow-up Analysis

We sought additional evidence for association of the top signals (P<2.0×10−6) identified in the discovery meta-analysis using two independent studies from the UK (UK Biobank, interim release, N = 119,613) and Iceland (deCODE, N = 72,613), respectively (Supplementary Tables 1–3). We used the same QC and analytical methodology as described above. We used the inverse-variance weighted fixed effects meta-analysis in METAL[108], to combine the discovery and follow-up association results. Significant associations were defined at P < 2×10−7 in the combined meta-analysis of discovery, UK Biobank and deCODE results.

Effect of study design

To investigate the potential effect of study design of the participating studies, we tested for heterogeneity between population-based, all case-control studies; T2D case-control studies (Supplementary Table 26). None of these comparisons showed significant evidence of heterogeneity (P<7.4×10−5, correcting for multiple testing).

Conditional analyses

The RareMETALS R-package[106] was used to identify independent BMI associated signals across the all-ancestry meta-analysis results in the discovery phase. RareMETALS performs conditional analyses by using covariance matrices from each individual cohort to distinguish true signals from the shadows of adjacent significant variants in LD. The conditional associations of all the variants within 1Mb of each R/LF coding variant were analyzed to identify [1] nearby secondary signals and [2] to determine independence from nearby non-coding variants or previously identified GWAS loci (previously defined as a window of 1Mb surrounding the lead SNP). Gene-based conditional analyses were also performed in RareMETALS.

Due to the selective coverage of variants on the ExomeChip, we also conducted the respective conditional analyses in the UK Biobank dataset that included 847,441 genome-wide genotyped markers, and 72,355,667 variants imputed against UK10k haplotype reference panel, merged with the 1000 Genomes Phase 3 reference panel. Where available, directly genotyped variants where used for conditional analyses. Otherwise, imputed variants with good imputation quality (IMPUTE2 info score > 0.6) were used. We used QCTOOL to extract variants of interest from the original imputed data set. Subsequently, GTOOL was used to convert to PLINK format (genotype calling threshold 0.99) and merged with the directly genotyped variants for conditional analyses in PLINK v1.90b3.35 64-bit (25 Mar 2016).

Conversions of effect size and explained variants

We assumed that 1 SD = 4.5 kg/m2 BMI-units, based on population based data, and 1.7m as the average height of a person to convert effects sizes in SD-units into body weight. The variance explained by each variant was calculated using the effect allele frequency (f) and beta (β) from the meta analyses using the formula[109] of explained variance = 2f(1-f)β2.

Penetrance analysis

We examined the penetrance for the four rare SNVs, p.Arg525Gln (rs56214831) in KSR2, p.Tyr35Ter (rs13447324) in MC4R, and p.Arg190Gln (rs139215588) and p.Glu288Gly (rs143430880) in GIPR in European ancestry data from the UKBiobank (N up to 120,000). For each variant, we compared the prevalence of underweight (BMI < 18.5 kg/m2), normal weight (18.5 kg/m2 ≤ BMI < 25 kg/m2), overweight (25 kg/m2 ≤ BMI < 30 kg/m2) and obesity (BMI ≥ 30 kg/m2) of non-carriers with non-carriers. We used a Pearson χ2 test to test for difference between distributions, and a χ2 for linear trend to test whether distributions of carriers were shifted compared to non-carriers. For p.Arg525Gln in KSR2 and p.Tyr35Ter in MC4R, we hypothesized that obesity prevalence was higher in carriers than in non-carriers, whereas for the two GIPR variants, we hypothesized that the prevalence of normal weight was higher in carriers than non-carriers.

Associations with obesity for the coding rare and low-frequency loci in children

For each of the 14 R/LF SNVs, we tested for association with childhood obesity in the CHOP cohort (Childhood Obesity: Early Programming by Infant Nutrition), the Severe Childhood Onset Obesity Project (SCOOP), the UK Household Longitudinal Study (UKHLS) and INTERVAL Study (INTERVAL). Summary statistics across the studies were combined using a fixed effects inverse-variance meta-analysis with METAL[108].

In the CHOP study, cases (1,358 boys, 1,060 girls) were defined as having a BMI > 95th percentile at any point in their childhood. Controls (1,412 boys, 1,143 girls) were defined as having < 50th percentile consistently through throughout childhood. The BMI percentiles are based on the CDC 2000 Growth Charts. All children were classified based on their BMI measurements between the ages of 2 and 18. All individuals are of European ancestry and were collected at the Children’s Hospital of Philadelphia. Informed consent was obtained from all study participants and study protocols were approved by the local ethics committees. Genotypes were obtained using the HumanHap550v1, HumanHap550v3, and Human610-Quad high-density SNP arrays from Illumina. The intersection of all SNPs on the arrays was used in all subsequent pre-imputation analyses. Before imputation, we excluded SNPs with a Hardy-Weinberg equilibrium P-value < 1.0×10−6, call rate of < 95% or MAF of < 1%. The genotypes were then pre-phased using Shapeit2 and imputed using the 1000 Genomes Phase 1 integrated variant set with Impute2. After imputation, SNPs were excluded if the INFO score was < 0.4. Boys and girls were analyzed separately using a logistic regression of case and control status, adjusting for three eigenvectors, and summary statistics were combined using a fixed effects inverse-variance meta-analysis with METAL[108].

SCOOP is a sub-cohort of the Genetics Of Obesity Study (GOOS) cohort. It includes >1,500 UK European ancestry individuals with severe, early onset obesity (BMI Standard Deviation Score > 3 and obesity onset before the age of 10 years), in whom known monogenic causes of obesity have been excluded (cases with MC4R mutations were excluded). Two case-control analyses with SCOOP cases were performed: 1) SCOOP vs. UKHLS for which array (Illumina HumanCoreExome) data was available, and 2) SCOOP vs. INTERVAL, for whom whole-exome sequencing data was available.

For the array based analyses, UKHLS controls were genotyped on the Illumina HumanCoreExome-12v1-0 Beadchip. SCOOP cases and 48 UKHLS controls were genotyped on the Illumina HumanCoreExome-12v1-1 Beadchip. The 48 overlapping UKHLS samples were used for quality control to ensure there were no systematic differences and bias between the two versions of the chip. SCOOP and UKHLS samples were phased with SHAPEITv2, and imputed with IMPUTE2 using the combined UK10K-1000G Phase III reference panel. For the WES analyses, SCOOP vs. INTERVAL controls were WES within the UK10K-EXOME project (Agilent v3) and the INTERVAL project (Agilent v5) respectively and were then jointly called and QC-ed on the union of the sequencing baits. Individuals overlapping or related between the array based and WES studies were removed.

After QC, 1,456 SCOOP and 6,460 UKHLS (BMI range 19–30), and 521 SCOOP and 4,057 INTERVAL individuals were available for the two analyses; all were unrelated, of high quality, and of European ancestry. For both analyses (i.e. SCOOP vs. UKHLS and SCOOP vs. INTERVAL), a maximum likelihood frequentist association test with the additive genetic model was implemented in SNPTEST v2.5. In the SCOOP vs. UKHLS analysis, sex and the first six PCs were included as covariates and variants with a SNPTEST INFO score <0.4 and HWE p<10−6 were removed. For the SCOOP vs INTERVAL analysis, we performed an unadjusted analysis (adjustment for PCs did not change sufficiently the results) and variants were limited to those covered at ≥7× in at least 80% of each sequencing cohort, meeting the VQSR threshold of –2.52, missingness <80%, HWE P-value<10−8, and GQ ≥30.

Cross-trait analyses

We evaluated each of the 14 R/LF SNVs for their association with other relevant obesity-related traits and conditions. We performed lookups in ExomeChip meta-analysis results from other consortia, including; our own GIANT consortium (height[15], WHR adjusted for BMI[24]), MAGIC (HbA1c, Fasting Insulin, Fasting Glucose, 2-hour glucose), GLGC (HDL-cholesterol (HDL-C), LDL-cholesterol (LDL-C), triglycerides and total cholesterol)), IBPC[40] (systolic and diastolic blood pressure), REPROGEN[23] (age at menarche and menopause) and GoT2D/T2D-GENES[16] (type 2 Diabetes). Associations were considered significant at P < 2.0×10−5, accounting for multiple testing.

Phenome-wide association analysis (PheWAS)

To evaluate the potential for pleiotropic effects for SNPs discovered from primary analyses, we performed phenome-wide association studies (PheWASs) using genotype and phenotype data from two independent sources of electronic health records (EHR): Vanderbilt University Medical Center Biorepository (BioVU) and the United Kingdom BioBank (UKBB). Phenotype selection and analysis strategy were synchronized across sites. A total of 1502 hierarchical phenotype codes from EHRs were curated by grouping International Classification of Disease, Ninth Revision (ICD-9) clinical/billing codes as previously described[110]. Phenotype codes with 20 or more cases and with minor allele count of 5 or greater in cases and controls were eligible for analysis. Series of logistic regression analyses were then performed in individuals of European ancestry for each eligible phenotype-genotype combination while adjusting for 5 genetic ancestry PCs. Odds ratios from genotype-phenotype combinations present in both BioVU and UKBB were then aggregated using inverse-variance weighted fixed-effects meta-analysis. Associations with p-values corresponding to false discovery rate (FDR) cut off of less than 10% were considered statistically significant.

Gene set enrichment analysis

We adapted DEPICT, a gene set enrichment analysis method for GWAS data, for use with the ExomeChip (‘EC-DEPICT’). DEPICT’s primary innovation is the use of “reconstituted” gene sets, where many different types of gene sets (e.g. canonical pathways, protein-protein interaction networks, and mouse phenotypes) were extended through the use of large-scale microarray data (see[111] 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 coding variants and only the genes directly containing those variants, rather than all genes within a specified amount of linkage disequilibrium (Supplementary Note).

Four analyses were performed for the BMI EC variants: [1] all coding variants with P<5×10−4, [2] all coding variants with P<5×10−4 independent of known GWAS variants[5], [3] all coding R/LF variants with P<5×10−4, and [4] all coding R/LF variants with P<5×10−4 independent of known GWAS variants. Affinity propagation clustering[3] 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 as representative for purposes of visualization (Supplementary Note). DEPICT for ExomeChip was written using the Python programming language (See URLs).

Drosophila RNAi knockdown experiments

For each of the 13 genes in which R/LF coding variants were associated with BMI, we searched for its corresponding orthologues in Drosophila in the ENSEMBL orthologue database. Orthologues were available for nine genes, but missing for ZBTB7B, MC4R, GIPR, and ZNF169. For each of the nine genes, we generated adipose-tissue (cg-Gal4) and neuronal (elav-Gal4) specific RNAi-knockdown crosses, leveraging upstream activation sequence (UAS)-inducible short-hairpin knockdown lines, available through the Vienna Drosophila Resource Center (VDRC). We crossed male UAS-RNAi flies and elav-GAL4 or CG-GAL4 virgin female flies. All fly experiments were carried out at 25 °C. Five-to-seven-day-old males were sorted into groups of 20, weighed and homogenated in PBS with 0,05% Tween with Lysing Matrix D in a beadshaker. The homogenate was heat-inactivated for 10 min in a thermocycler at 70 °C. 10µl of the homogenate was subsequently used in triglyceride assay (Sigma, Serum Triglyceride Determination Kit) which was carried out in duplicates according to protocol, with one alteration: the samples were cleared of residual particulate debris by centrifugation before absorbance reading. 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.

Enrichment analysis in monogenic genes of obesity

We identified 39 genes with strong evidence that disruption causes monogenic or syndromic forms of obesity (Supplementary Table 21). To test whether these genes are enriched for R/LF coding variant associations with BMI, we conducted simulations by matching each of the 39 genes with other genes based on gene length and number of variants tested, to create a matched set of genes. We generated 1,000 matched gene sets from our data and assessed how often the number of R/LF coding variants that exceeded given significance thresholds was greater in our monogenic/syndromic obesity gene set compared to the matched gene sets.

DATA AVAILABILITY

Summary statistics can be downloaded from http://portals.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium