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

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

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 . 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 . 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 ). 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. 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 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 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. 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 ). 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. 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). 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). We assumed that 1 SD = 4.5 kg/m 2 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 = 2 f (1- f ) β 2 . 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/m 2 ), normal weight (18.5 kg/m 2 ≤ BMI < 25 kg/m 2 ), overweight (25 kg/m 2 ≤ BMI < 30 kg/m 2 ) and obesity (BMI ≥ 30 kg/m 2 ) 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. 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 > 95 th percentile at any point in their childhood. Controls (1,412 boys, 1,143 girls) were defined as having < 50 th 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. 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. 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. 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). 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. 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.

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 max =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 max =72,613) and UK Biobank ( N max =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, N max =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. 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 2 =1; D’ =1) and two in ACHE ( r 2 =0.98; D’ =0.99) were in high LD, whereas the two SNVs in GIPR ( r 2 =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. 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 ). 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 ). 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/m 2 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/m 2 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 (effect mean =0.019 SD/allele, ~0.086 kg/m 2 or ~0.247 kg/allele) of which the largest effect is seen for the FTO locus (0.08 SD/allele, ~0.35 kg/m 2 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/m 2 , or 0.130 to 0.728 kg) 5 . 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/m 2 ) 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 (N max =119,781), we compared the prevalence of normal-weight (18.5 kg/m 2 ≤ BMI < 25 kg/m 2 ) and obesity (BMI ≥ 30 kg/m 2 ) 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 cases =4,395; N controls =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 binomial =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 . 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 WHR adjBMI 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. 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. 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. 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. 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 (r 2 =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.

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 WHR adjBMI 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 .