Evaluation of the MC4R gene across eMERGE network identifies many unreported obesity-associated variants

Background/Objectives Melanocortin-4 receptor (MC4R) plays an essential role in food intake and energy homeostasis. More than 170 MC4R variants have been described over the past two decades, with conflicting reports regarding the prevalence and phenotypic effects of these variants in diverse cohorts. To determine the frequency of MC4R variants in large cohort of different ancestries, we evaluated the MC4R coding region for 20,537 eMERGE participants with sequencing data plus additional 77,454 independent individuals with genome-wide genotyping data at this locus. Subjects/Methods The sequencing data were obtained from the eMERGE phase III study, in which multisample variant call format calls have been generated, curated, and annotated. In addition to penetrance estimation using body mass index (BMI) as a binary outcome, GWAS and PheWAS were performed using median BMI in linear regression analyses. All results were adjusted for principal components, age, sex, and sites of genotyping. Results Targeted sequencing data of MC4R revealed 125 coding variants in 1839 eMERGE participants including 30 unreported coding variants that were predicted to be functionally damaging. Highly penetrant unreported variants included (L325I, E308K, D298N, S270F, F261L, T248A, D111V, and Y80F) in which seven participants had obesity class III defined as BMI ≥ 40 kg/m2. In GWAS analysis, in addition to known risk haplotype upstream of MC4R (best variant rs6567160 (P = 5.36 × 10−25, Beta = 0.37), a novel rare haplotype was detected which was protective against obesity and encompassed the V103I variant with known gain-of-function properties (P = 6.23 × 10−08, Beta = −0.62). PheWAS analyses extended this protective effect of V103I to type 2 diabetes, diabetic nephropathy, and chronic renal failure independent of BMI. Conclusions MC4R screening in a large eMERGE cohort confirmed many previous findings, extend the MC4R pleotropic effects, and discovered additional MC4R rare alleles that probably contribute to obesity.


Introduction
The melanocortin-4 receptor (MC4R, OMIM:155541) is a G protein-coupled receptor that is critical in leptin-melanocortin pathway. The protein MC4R is predominantly expressed in the hypothalamus and is involved in regulation of satiety, feeding behavior, and energy homeostasis [1,2]. MC4Rknockout mice are hyperphagic with a reduced metabolic rate and elevated body weight [3] with body mass increased by 7-45% in heterozygous and 50-100% in homozygous genotypes in comparison to wild type [4,5]. In humans, more than 170 distinct rare variants are associated with early onset obesity and hyperphagia [6]. MC4R-associated obesity is the most common monogenic form of obesity with a reported prevalence of up to 6% [7]. However other reports have shown lower prevalence [8,9]. In the last two decades, in vitro assays have been developed in order to determine the biosynthesis, cell-surface expression, ligand binding, and Gs activation (measurement of cAMP) of MC4R and its naturally occurring variants [10,11]. However, the observed functional defects are sometimes discordant with obesity due to complex gene environmental interactions and variable expressivity of dominant inheritance [12]. Moreover, gain-of-function variants are also known. From a meta-analysis, V103I appears to be protective against obesity (OR 0.69) [13,14]. The V103I variant may reduce receptor internalization [15]. Another missense variant I251L originally classified as functionally neutral shows increased MC4R basal activity through alteration of cAMP signal transduction that protects against obesity (OR = 0.5) [16]. Indeed, in a large United Kingdom Biobank study of 450 K individuals, 11 out of 61 studied variants had potential gain-of-function properties that await confirmation [15]. Because of this heterogeneity and that most original studies were limited to severely obese families that may overestimate the penetrance estimates, larger studies from unselected populations are needed to better elucidate and confirm prior results. Moreover, apart from rare coding variants, common noncoding regulatory variants upstream of MC4R gene are known determinants of BMI variation at the population level that add complexity [17][18][19].
The Electronic Medical Record (EMR) is a rich source of clinical information. In 2007, The electronic MEdical Records and GEnomics (eMERGE) network was initiated by the National Human Genome Research Institute (NHGRI) to explore the utility of DNA biobanks linked to EMRs for research [20]. Recently, the network developed protocols to perform genetic sequencing of 109 of the most clinically relevant genes including MC4R [21].
In this study, we evaluate and catalog the MC4R sequencing and genotyping data from eMERGE III participants to further study this genomic region and its association with obesity and related phenotypes.

Study cohort and sequencing
Protocols for this study were approved by the Institutional Review Boards at each institution; all included participants provided written informed consent at study enrollment.

eMERGE-seq
The sequencing data obtained from the eMERGE III Network consist of 24,956 participants from 11 US sites. The network developed protocols to sequence 109 genes including MC4R [21]. Two CAP/CLIA (College of American Pathologists (CAP), Clinical Laboratory Improvement Amendments (CLIA)) certified DNA sequencing laboratories, Baylor College of Medicine Human Genome Sequencing Center and Broad Institute and Partners Laboratory for Molecular Medicine were responsible for sequencing, data harmonization and quality control (QC) analyses. Details of these procedures have been reported previously [21].

eMERGE-GWAS
Apart from sequencing data, postimputation whole-genome genotyping data for additional 77454 independent participants from eMERGE network were available to us (dbGAP (phs000888.v1.p1)). The imputation process and genotype QC in eMERGE followed guidelines that have been published previously [22]. The QC process included sample call rates, sample relatedness, population stratification, and sex inconsistency as well as marker quality (i.e., marker call rate, minor allele frequency (MAF), and Hardy-Weinberg equilibrium (HWE)). All common variant analyses were limited to participants with call rates ≥ 98%, variants with call rates ≥ 99%, as well as variants with MAF ≥ 1% and HWE P ≥ 0.00001.

Phenotype data
Demographic and anthropomorphic measures from all participants were accessed from eMERGE coordinating center and included in Table 1(a), (b). First, available BMI (kg/m 2 ) data across all sites with a total of 3,368,260 entries were investigated for missing data, lab entry errors, and data inconsistency. The algorithm and initial screening indicated that BMI values above 100 kg/m 2 most likely were due to lab entry errors in these cohorts and therefore excluded. On average, there were 30 BMI records per participant with life span duration of 10 years. Next, the mean and median BMI per participant was calculated. The median BMI per individual was then used for common variant genetic analyses. For participants with rare MC4R variants, all sites were recontacted to retrieve any additional missing BMI reports, and check for data inconsistency. For this group, the highest post-QC BMI was used for penetrance estimation and we made sure that temporary physiologic conditions that influence BMI such as pregnancy, ascites, or edema were not present at the same time. BMI-for-age percentile was separately calculated for pediatric participants using CDCbased online resource (https://www.cdc.gov). After removing missing BMI data, 20,537 individuals from eMERGEseq cohorts and 77,454 from eMERGE-GWAS were used for the analyses. As shown in Table 1(a), the mean age of adults and pediatrics were 58.47 and 11.73, respectively, (Table 1(a)).

Penetrance
Since, the eMERGE participants were not preselected for obesity, population allele frequency information was used  to estimate penetrance as described previously (URL: http:// cardiodb.org/allelefrequencyapp/) [23]. The penetrance was defined as case allele frequency divided by control allele frequency multiplied by disease prevalence with a penetrance range of 0-1. We used traditional binary classifications to define obesity (BMI ≥ 30, pediatrics BMI ≥ 95%) and overweight (BMI ≥ 25, pediatrics BMI percentile ≥ 85%) groups, respectively. The overall prevalence of obesity in eMERGE participants was 30% consistent with US data used for penetrance estimation [24].

Annotation and genetic analysis
All detected coding MC4R variants were annotated using latest version of variant annotation program (Annovar) [25]. Annotations per variant include latest ClinVar reports, evidence of previous publications in connection with obesity, exonic function, conservative prediction algorithm scores of SIFT, PolyPhen-2, MutationTaster, and PRO-VEAN [25]. The corresponding allele frequency per ancestry and all together for the Exome Aggregation Consortium (ExAC) and The Genome Aggregation Database (gnomAD) were also cataloged for comparison.

GWAS analyses
For evaluation of common variants surrounding MC4R and genome-wide association of BMI, quantitative linear regression analyses were performed adjusting for site of genotyping (11 sites), principal components derived from genomic data (ten PCs), age, and sex using second generation of PLINK [26]. Variants with MAF above 1% that overlapped with eMERGE-GWAS successfully merged and re-imputed and only variants with high imputation info score (r 2 quality score > 0.7) were selected for GWAS analyses. The HAPMAP (haplotype map of the human genome) reference population frequency and linkage disequilibrium (LD) in different ancestries were obtained using LocusZoom (http://csg.sph.umich.edu/locuszoom) and LD statistics (r 2 ) between a pair of SNPs was calculated using PLINK [26].
At the next step, selected GWAS hits were functionally evaluated using FUMA (Functional Mapping and Annotation of Genome-Wide Association Studies) and Haplo-R [27,28]. For the unreported variants in this study we also evaluate the effect of mutation on protein stability using STRUM, a machine learning-based mutation stability predictor algorithm [29]. This new algorithm measures the unfolding free energy difference between the wild type and mutant protein (ΔΔG) (https://zhanglab.ccmb.med.umich. edu/STRUM) [29]. In general, ΔΔG values below zero means that the mutation contributes to protein destabilization.

Phenome-Wide Association Study (PheWAS) analyses
A PheWAS was also performed in order to evaluate pleotropic effects of variants in MC4R region with any other trait. The detail of methodology is described in previous publications [30]. We used the PheWAS package in R version 3.5.1 [30]. The trait definition in PheWAS approach is based on billing International Classification of Diseases (ICD) codes. ICD9 codes are collapsed into PheWAS codes (phecodes) according to the PheWAS map; then cases and controls are determined according to the code under study. We used a threshold of at least 20 cases for the code to be included in the model to have sufficient power for 1789 phecodes. Next, for each PheWAS code, a logistic regression model was created and adjusted for age, sex, and PCs. A false discovery rate (FDR) of 0.05 using the Benjamini-Hochberg method was then used to correct the threshold for multiple hypotheses testing.

Burden test approach
In the eMERGE-seq population, PheWAS analysis was performed on extremely rare variants (MAF < 0.1%) using burden test (SNP-set (sequence) kernel association test (SKAT-O)) procedure in R, which aggregates individual score test statistics of variants while adjusting for covariates [31]. To avoid extreme case-control ratios, our criteria include subgroup of ICD9 code of ≥20 sample size in which ≥2 MC4R carriers, regardless of variant, were present for each trait code. We excluded the four frequent MC4R coding variants with reported protective or no effect (V103I, I251L, F202L, I198I). ICD9 codes related to injuries or poisoning were removed leaving a total of 1967 ICD9 codes. In this exploratory analysis, statistical significance was determined using the Bonferroni correction in which 1967 ICD9 codes remained for analyses with a P value threshold of (0.05/1967 = 2.55 × 10 −5 ) for a single gene set.

Power analyses
QUANTO software was used for power estimate of individual variants [32]. In GWAS analyses using BMI in linear regression model, given this small genomic region with only two independent haplotypes and large sample size (N = 97,991) we had >90% power to detect associations for rare variants (MAF ≥ 0.01 at Beta ≥ 0.5) and 99% power for common variants (MAF ≥ 0.2 at Beta ≥ 0.3) in an additive model. In PheWAS analyses for the binary outcomes in logistic regression, we calculated power for each variantphecode pair given type 1 error rate of (α = 0.05/1789 × (2 haplotypes) = 1.39 × 10 −5 ). We had more than 80% power for all reported phenotypes with MAF > 1%. Phe-WAS of extreme rare variants, using SKAT, was considered exploratory due to limited statistical power (12-30% power) and we reported all findings with the Bonferroni correction mentioned above in Supplementary Table S11.

Results
After removing individuals with missing BMI measurements, this study included 20,537 eMERGE-seq participants with sequencing data for MC4R and additional 77,454 independent participants with postimputed genotyping data from eMERGE-GWAS. The demographic and detail of BMI distribution per race ethnicity and sex after QC are shown in Table 1(a), (b). The median BMI histogram plot for all participants is shown in Supplementary Fig. S1, which closely matched statistics reports in the United States (overall mean = 27.50, SD = 7.49) [24].

Sequencing results
Sequencing data of MC4R revealed 125 variants among 1839 eMERGE-seq participants (out of the possible 24,956 sequenced, 7.3%). From these, 20,537 had post-QC BMI data. The demographic distribution of these participants as well as the frequency of BMI ≥ 30 (obesity) and BMI ≥ 25 (overweight) per ancestry are included in Supplementary  Table S1. The detected variants in MC4R include 85 nonsynonymous, 28 synonymous, two frameshifts, two start-loss, five stop-gains, one at 3′ UTR, and two in 5′ UTR region. Sixty two of these have been previously reported or studied in connection with obesity. A comprehensive annotated overview of all variants is included in Supplementary Table S2. Overall, the allele frequencies of all variants of study participants were comparable to the reported public resources per ancestry such as ExAC and gnomAD as shown Supplementary Table S2. For clarity, we classified these variants into the three major groups:

Group 1: pathogenic or likely pathogenic variants according to ClinVar classification
We identified 17 variants in this category for 74 MC4R carrier participants; 39 with BMI ≥ 30 (60%), 61 with BMI ≥ 25 (94%), and 9 with missing BMI (Supplementary Table S3). As shown, highly penetrant variant for obesity (BMI ≥ 30) include I269N (2 out of 2: meaning 2 participants had BMI ≥ 30 out of possible 2 carrier individuals), P299H (3 out of 3), I170V (7 out of 12) (pediatrics (5 out of 5)), R156Q (2 out of 2), and Q156X (3 out of 4). Consistent with prior reports, the pathogenic Y35X stop-gain variant was in complete LD with the D37V substitution for which we detected seven compound heterozygotes [33]. As shown, four had obesity and all seven were overweight (BMI ≥ 25). Apart from obesity, three of the seven had lipid disorders, the most frequent shared diagnosis (ICD9 272.4). Of note, one participant was a triple heterozygote (D37V-Y35*-I69K) with obesity (BMI = 33.1) and hyperlipidemia. The list of all compound heterozygotes and estimated penetrance per ancestry are in Supplementary Table S6. Group 2: Undetermined, conflicting according to ClinVar, but previously linked to obesity in at least one published study This group consists of 45 variants (2 synonymous, 42 nonsynonymous, and 1 frameshift deletion) in which four variants had allele frequencies above 0.1% (Supplementary Table S4). The detected four frequently observed variants in MC4R coding region include three nonsynonymous variants of V103I (MAF = 1.6%), I251L (MAF = 0.8%), F202L (MAF = 0.1%), and one synonymous variant I198I (MAF = 0.6%). Of note, the frequency of synonymous coding variant I198I in AA was 3.7% comparable to public estimates (ExAC African = 3.4%, gnomeAD_African = 4.1%, (Supplementary Table S4)). In addition, consistent with prior reports, the F202L variant was exclusively present in combination with I198I disregarding race in which we identified 90 compound heterozygote participants (Supplementary Table S6). No significant difference in frequency of obesity between I198I-only group and I198I-F202L compound heterozygote subgroup were observed. However, we found two exceptional cases: one African American individual who was a double homozygote for I198I-F202L, an 18-year-old male with BMI of 40.2 kg/m 2 or 99.6 percentile who also carried the diagnosis of chronic asthma. Another 2.5-year-old African American male with congenital nystagmus was a triple heterozygote for I198I-F202L-I251L with BMI in the 95th percentile.

Group 3: previously unreported variants
We identified 30 variants (26 nonsynonymous, 1 frameshift, 2 start-loss, and 1 stop-gain) that to our knowledge have not been previously linked to obesity and are predicted to result in loss of function by at least one functional algorithm or result in a frameshift with premature truncation (Table 2(a)). In addition, protein stability scores (ΔΔG) were negative for most of these variants indicating destabilization (Table 2(a)). Two additional nonsynonymous variants were identified, predicted as tolerant and listed in Table 2(a) with obesity    (L325I, A114V) and a negative ΔΔG. As shown in Supplementary Table S2, some of these unreported variants previously were detected in ExAC or gnomeAD databases, but to our knowledge were not previously linked to obesity in any publication. Overall, these unreported variants were observed in 36 participants.

GWAS analysis
We used post-QC median BMI as a quantitative trait and performed linear regression GWAS analyses in the MC4R region for all eMERGE participants (N = 97,991) adjusted for principal components (ten PCs), age, sex, and site of genotyping. This analysis identified one common risk haplotype and another novel and relatively rare protective haplotype ( Fig. 1, Supplementary Table S7). The common and known risk haplotype upstream of MC4R gene spans 170 KB in Europeans with the best variant rs6567160 in our study (P = 5.36 × 10 −25 , Beta = 0.37) (Fig. 1a, Supplementary Table  S7). Other previously published variants are viable proxies (r 2 ≥ 0.9) with this variant including rs571312, rs523288 rs12967135, and rs17782313. We also analyzed separately pediatric only population and report the results in Supplementary Table S7. Of note, the effect was consistent in the pediatric only population (best variant rs1942860, P = 2.07 × 10 −12 , Beta = 0.48, Supplementary Table S7). Ancestry specific GWAS analyses also are shown in Fig. 1a-d. Similar effect was detected in African American but at lower magnitude where the best marker was rs11664369 (P = 4.89 × 10 −6 , Beta = 0.51). An effect upstream of MC4R was not detectably significant in Asians population and was weak in Hispanic ethnicity participants (Fig. 1c, d).
In addition to this common risk haplotype, a novel rare protective haplotype spanning 63 Kb near MC4R was detected with no LD with the upstream common risk  haplotype (Fig. 1e). This haplotype had overall frequency of 1.8% across all ancestries and encompasses the coding variant rs2229616 (V103I) (P = 6.23 × 10 −8 , Beta = −0.62) ( Supplementary Table S7). Importantly, because of lack of LD between rs2229616 (V103I) and the main common risk variant rs6567160 (r 2 = 0.0001), the protective effect remained significant after conditioning on rs6567160 (P = 8.77 × 10 −7 ). The list of variants in this rare haplotype are shown in Supplementary Table S7; 38 intergenic markers are proxies for rs2229616 (V103I). Next, we evaluated the potential functional effects of these variants and those included in Supplementary Table  S8. Allele specific expression according to resources for brain tissue, psychENCODE, and the CommonMind Consortium indicate negative eQTL effects of the risk alleles of top markers (rs523288-T, rs6567160-C, rs11664369-T, rs17782313-C) for MC4R in brain tissue. These variants are also listed as eQTL in GTEx/v8, but no results for brain tissue are reported. In addition, markers with chromatin or histone mark regulatory effects in brain tissue according to HaploReg V4.1 includes rs2229616 (V103I) and listed in Supplementary Table S8 [28].

PheWAS analyses
Pleotropic effects of the common and rare variants in MC4R region against available EMR disease traits in all eMERGE participants were evaluated by PheWAS.
For common variants, PheWAS analyses independently replicated the GWAS findings with strong association of multiple MC4R variants to ICD9 codes related to overweight, obesity, and morbid obesity (best P = 6.74 × 10 −13 ) ( Table 4,  Supplementary Table S9). Apart from obesity, PheWAS association with common upstream variants includes dysmetabolic syndrome X, chronic venous insufficiency, and malignant neoplasm of small intestine (Table 4).
Furthermore, the protective effect of V103I (rs2229616), against obesity was replicated using this approach (P = 1.39 × 10 −6 , Beta = −0.38). This variant also shows protective effect against type 2 diabetes, diabetic nephropathy, hypertensive nephropathy, and chronic kidney disease (Table 4). Another coding variant, the marker rs121913563 (A175T) was associated with the presence of casts and cells in urine (Supplementary Table S9). In order to avoid intercorrelated phenotypic associations that are often linked to obesity, we then reanalyzed the PheWAS controlling for BMI as another covariate. The protective effect of V103I against type 2 diabetes (P = 6.88 × 10 −6 ) and kidney disease (P = 8.22 × 10 −6 ) remained significant after controlling for BMI indicating multiple independent effects for this variant (Supplementary Table S10).
For extreme rare variants (MAF < 0.1%) that were limited to the eMERGE-seq population, PheWAS was performed using (SKAT-O) a burden test procedure which aggregates individual score test statistics to improve power. Abnormal glucose test (N = 33), colonic polyps (N = 29), dysmetabolic  Supplementary Table S5.

Discussion
In the present study, we performed detailed annotation and penetrance estimation of 125 rare variants in MC4R detected in 1839 of 24,956 eMERGE participants with sequencing data. In addition, we studied 77,454 independent eMERGE participants with whole-genome genotyping data for GWAS and PheWAS of common variants in this region. MC4R variants represent the most frequent cause of monogenic obesity. In our collection, we found that~7.3% of the total population and 11.3% of obese participants (BMI ≥ 30), carried at least one coding variant in MC4R coding region, higher than previous reports of 5 to 6% in obese patients [7] (Table 5). However, representation of coding variants varies between studies influencing the  overall prevalence. Furthermore, not all rare MC4R variants are deleterious, hence the difference in prevalence of rare MC4R variants in obese and nonobese can be driven by both pathogenic and benign variants. As shown in Table 5, excluding the known V103I common coding variant with protective effect, gives the overall rate of 4.2% in our population, while increasing the fold difference of prevalence estimates especially when comparing normal vs extreme obese subgroup (OR = 7.17, 95% CI = 5.81-8.87, Table 5). Further restricting to only 17 known pathogenic variants detected in this study, yields a more than 95% overweight rate (Supplementary Table S3). Because of this heterogeneity, and other confounding factors such as the effect of race or age, we estimate penetrance per each variant in each race and age group and include in Supplementary Table S2 to contribute to future studies for collective penetrance estimations.
In this study, we detected 30 naturally occurring MC4R variants that to our knowledge have not been previously reported to be associated with obesity and are predicted to be damaging by at least one functional algorithm. These unreported variants are observed in 36 carrier individuals with 54% obesity rate. Importantly, seven adult participants carrying novel rare variants of L325I, E308K, D298N, F261L, T248A, D111V, and Y80F had obesity class III (BMI ≥ 40) and two unreported variants (L325I and S270F) were detected in two obese participants. As shown in Table 2(a), apart from predicted to be deleterious based on conservation score algorithms (SIFT and others), most of these variants also induce protein destabilization effects with negative ΔΔG scores [29].
In African Americans, the synonymous I198I variant can be seen up to 4%. In our study population, it has a frequency of 3.7% in African American in which 141 out of 242 individuals were obese (58%). Consistent with previous reports, all individuals harboring the F202L variant were compound heterozygotes with I198I. In vitro functional studies for F202L suggest both decrease or no influence on MC4R receptor activity and were observed in both obese and nonobese individuals [15,34]. We found an exceptional double homozygote African American male (I198I-F202L) with BMI of 40.2 kg/m 2 consistent with an allelic dosage effect. Other detected variants in African Americans include R305W, N240S, N123S, and R7H, all previously linked to obesity in the same ancestry [35].
In Hispanics, G323E previously reported in Iberians, was also detected in three Hispanics including an obese pediatric individual (Supplementary Table S2), [36]. I269N is another known pathogenic variant exclusively seen in Hispanic ethnicity with penetrance of 100% for obesity in our collection (Supplementary Table S2) [15].
Variants exclusively present in Asians include C277X, I195T, S171G, I169V, Q115X, Y35C, and R7R. C277X has been previously linked to severe obesity in the Chinese Han [37]. Similarly, we detected this variant in one adult Asian overweight male (BMI = 26.4) with hyperlipidemia (Supplementary Table S2). Y35C has been previously detected in both obese and nonobese Chinese, two of our five Asian carriers were obese (Supplementary Table S2) [38]. However, in our collection, most of Asian carriers were either lean or overweight rather than obese indicating lower penetrance of these variants in Asian or influence of other genetic or environmental factors.
Of note, in Asian populations, the WHO has recognized lower BMI cutoffs as a trigger for increased health risks [39].
In GWAS analyses, the known common risk haplotype upstream of MC4R was confirmed with the best intergenic marker rs6567160 (P = 5.36 × 10 −25 , Beta = 0.37). This variant as well as its proxy markers (e.g., rs17782313) near the MC4R gene have been associated with obesity in previous publications [18,19]. To test the hypothesis that the common variants upstream of MC4R also are important determinants of BMI in children, we separately analyzed the pediatric population. The results were consistent with the adult findings with no evidence of heterogeneity (for marker rs6567160, P = 3.02 × 10 −12 , Beta = 0.48 (Supplementary  Table S7)). Indeed, previous GWAS studies suggest a strong overlap between the genetic architecture of childhood and adult BMI [40]. This justifies our approach to combine all data together while adjusting for age and other covariates. However, we caution that BMI measures in children that are under active growth and development, even when adjusted using best methods may not be precisely comparable to adult values for the purposes of defining groups for analysis.
Interestingly we noted negative eQTL of the BMIincreasing risk alleles for MC4R expression in brain tissue which could functionally explain a milder form of MC4R deficiency for these common variants, consistent with monogenic mutation effects. To our knowledge, the only functional study for MC4R upstream variants was performed in nonbrain tissue (intestinal tissue) showing increased eQTL expression of MC4R for BMI-risk allele (rs17782313-C) [41].
PheWAS of extreme rare variants (MAF < 0.1%), using burden test methodology, provided an initial exploratory association of the MC4R coding variants with abnormal glucose test, presence of casts in urine, dysmetabolic syndrome X, pituitary hypofunction, and neurofibromatosis. In at least one study, signs and symptoms associated with neurofibromatosis have been reported in one case with MC4R deficiency [42]. Additional studies are needed to explore these findings in larger cohorts in order to detect true signals with sufficient statistical power.
We studied a large, unbiased, geographically diverse population, confirmed previous findings, extended associations to other related phenotypes, and added 30 new findings that warrant further evaluation in independent cohorts. The penetrance estimation for BMI and obesity from this collection will assist variant reclassification for MC4R and improve clinical decision making in the context of personalized medicine. These results are also relevant for targeted drug development of improved next generation melanocortin agonist therapies. Author contributions All authors agree to be accountable for all aspects of the work, and read and approved the final manuscript. BN, IBS, and DRC: study concept and design, data review and comments, interpretation of data, data analysis, study supervision, and manuscript preparation. BN, TL, TM, FDM, BB, OD, XN, NS, AHS, DRC, RS, JGN, IBS, and DJC: data acquisition, data preparation, and data analyses. BN, MSW, JCD, DRC, IJK, HH, GPJ, and JBH: organization, data preparation, manuscript review, and critique.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.