Unsupervised Clustering of Missense Variants in HNF1A Using Multidimensional Functional Data Aids Clinical Interpretation.

Exome sequencing in diabetes presents a diagnostic challenge because depending on frequency, functional impact, and genomic and environmental contexts, HNF1A variants can cause maturity-onset diabetes of the young (MODY), increase type 2 diabetes risk, or be benign. A correct diagnosis matters as it informs on treatment, progression, and family risk. We describe a multi-dimensional functional dataset of 73 HNF1A missense variants identified in exomes of 12,940 individuals. Our aim was to develop an analytical framework for stratifying variants along the HNF1A phenotypic continuum to facilitate diagnostic interpretation. HNF1A variant function was determined by four different molecular assays. Structure of the multi-dimensional dataset was explored using principal component analysis, k-means, and hierarchical clustering. Weights for tissue-specific isoform expression and functional domain were integrated. Functionally annotated variant subgroups were used to re-evaluate genetic diagnoses in national MODY diagnostic registries. HNF1A variants demonstrated a range of behaviors across the assays. The structure of the multi-parametric data was shaped primarily by transactivation. Using unsupervised learning methods, we obtained high-resolution functional clusters of the variants that separated known causal MODY variants from benign and type 2 diabetes risk variants and led to reclassification of 4% and 9% of HNF1A variants identified in the UK and Norway MODY diagnostic registries, respectively. Our proof-of-principle analyses facilitated informative stratification of HNF1A variants along the continuum, allowing improved evaluation of clinical significance, management, and precision medicine in diabetes clinics. Transcriptional activity appears a superior readout supporting pursuit of transactivation-centric experimental designs for high-throughput functional screens.

Introduction

Precision medicine increasingly relies on an accurate interpretation of the consequence of genetic variation. Large-scale multi-ethnic genetic sequencing studies have challenged our understanding of the relationship between coding variants in Mendelian disease genes, including those involved in monogenic forms of diabetes such as HNF1A (MIM: 142410). Until relatively recently, the consensus has been that heterozygous highly penetrant loss-of-function alleles in HNF1A give rise to a clinically distinct diabetes subtype, characterized by an early age of onset (typically < 25 years), dominant inheritance, sensitivity to sulphonylureas, and non-obesity, and termed HNF1A maturity-onset diabetes of the young (HNF1A-MODY [MIM: 600496]). While this genotype-phenotype correlation is true for a subset of HNF1A variant carriers, it represents one end of a broad spectrum of HNF1A variant effects.2, 3, 4

Genome-wide association and next-generation sequencing studies of randomly ascertained individuals have challenged binary assumptions and overinflated pathogenicity estimates regarding variants in HNF1A (and other Mendelian disease genes) and identified common coding variants of low effect associated with increased risk of type 2 diabetes (MIM: 125853).5, 6, 7, 8 Whole-exome sequencing studies in populations of Mexican American ancestry have revealed a low-frequency missense variant (c.1522G>Α [p.Glu508Lys]) in HNF1A associated with a 5-fold increase in type 2 diabetes prevalence. These complex genomic insights warrant a more nuanced understanding of the phenotypic manifestation of HNF1A gene variants: some alleles are sufficient for early-onset sulfonylurea-responsive diabetes (HNF1A-MODY), although it should be noted that not all carriers of these alleles get early-onset diabetes; not even diabetes at all. Moreover, some alleles modify susceptibility for developing complex multifactorial hyperglycemia later in life (type 2 diabetes), and most alleles will likely manifest as benign and neutral.

A correct diabetes diagnosis is important because a mutation in HNF1A leads to clinical actions involving diagnosis, treatment, and genetic counselling. Individuals with rare, deleterious HNF1A alleles and young-onset diabetes (typically < 25 years) are sensitive to treatment with oral sulfonylureas and can often avoid insulin injections until late in life.,

The ubiquity of genetic sequencing means that more novel and incidentally detected variants of uncertain clinical significance (VUS) will be identified in individuals with less extreme phenotypes. The challenge today is in the ability to map HNF1A sequence-function relationships at high fidelity, using clinical and molecular characterization and analytical pipelines with sensitivity to capture the subtleties along the pathophysiological continuum. Rigorous functional follow-up of rare sequence-identified alleles in HNF1A is crucial to making correct assignments of pathogenicity. Indeed, functional data are considered a strong line of evidence for accurate clinical diagnostic classification of variants. Furthermore, it has been shown that diabetes severity in HNF1A variant carriers is influenced by allele position in the gene: the transactivation domain is more tolerant to genetic variation and variants in the latter exons8, 9, 10 are only present in hepatocyte-dominant isoforms and would thus not likely translate to a strong beta-cell phenotype.,

To understand the relationship between HNF1A sequence variation, molecular dysfunction, and clinical phenotype, we characterized the functional impact of a total of 73 HNF1A missense variants detected in the exomes of 12,940 multi-ethnic type 2 diabetes case subjects and control subjects using standard functional assays. Our primary objective was to develop an analytical approach that would enable (1) an unbiased and comprehensive evaluation of HNF1A variant behavior based on multiple molecular mechanisms and (2) sensitive mapping of multi-dimensional in vitro function to HNF1A glycemic phenotypes in vivo. We hypothesized that severity of molecular dysfunction in vitro (wild-type/wild-type-like, moderate/intermediate impact, loss-of-function/deleterious) would correlate positively with the severity of clinical phenotype (benign, increased type 2 diabetes risk, young-onset sulfonylurea-responsive hyperglycemia).

Material and Methods

The study was approved by the regional ethical committee in Bergen (#2009/2079). We investigated the function of all rare (MAF < 0.5%) and low-frequency (0.5% < MAF <. 5%) as well as three common (MAF > 5%) HNF1A nonsynonymous missense variants (n = 73) identified in an exome sequencing study of 12,940 type 2 diabetes case subjects and control subjects from five different ancestry groups (Figure S1, Table S1). Collectively, the variants did not enrich for a type 2 diabetes phenotype under any of the several variant filters used (MAF < 0.1%, conserved and predicted damaging [PolyPhen: SKAT p value = 0.30, BURDEN p value = 0.37]).

Bioinformatic Prediction

The following four in silico tools were used to evaluate the pathogenicity of the alleles: SIFT, PolyPhen-2, MutationTaster, and Combined Annotation Dependent Depletion (CADD). A CADD cut-off score of 15 was used (>15, pathogenic).

Functional Characterization

The individual effects of the 73 HNF1A missense variants were functionally investigated by two research teams at the Universities of Oxford (UK) and Bergen (Norway) using four different molecular assays (see detailed description of assays below). Using two laboratories allowed us to evaluate the robustness of the functional studies. Each laboratory assessed a unique set of exome-detected variants (n > 30), a shared subset of exome-detected variants (n = 5), shared type 2 diabetes risk variants (n = 2), as well as shared HNF1A-MODY reference variants (positive controls, n = 6) (Figures S2 and S3). The positive controls were selected on the basis of previously reported functional data supporting pathogenicity, clinical evidence for causality (sulfonyurea sensitivity in multiple carriers), and/or genetic (co-segregation) evidence to support their role pathogenesis (Table S2). Plasmid and HNF-1A variant constructs, transactivation assays, HNF-1A protein abundance, subcellular localization, and DNA binding are detailed below.

Plasmid and HNF-1A Variant Constructs

A construct encoding the human HNF1A cDNA (GenBank: NM_000545.6) in cDNA3.1 His/C plasmid was used as wild-type and template for introducing HNF1A variants using the QuikChange XL Site-directed Mutagenesis Kit (Stratagene). The wild-type sequence used in this study also harbors the common coding variant c.79A>C (p.Ile27Leu) (MAF = 34.8%) and a common synonymous variant c.51C>G (p.Leu17=) (MAF = 46.5%). All constructs were verified by DNA Sanger sequencing.

For transactivation experiments, two reporter constructs were used: (1) pGL3-RA, containing the promoter of the rat albumin gene (nucleotide −170 to +5) next to the Firefly luciferase gene in vector pGL3-Basic (Promega) (kindly provided by Prof. Graeme I. Bell, University of Chicago, Chicago, IL, USA) and (2) pGL3-HNF4AP2, which contains the human HNF4A (MIM: 600281) P2 promoter (nucleotide −418 to +13) next to the Firefly luciferase gene (kindly provided by Prof. Maria Angeles Navas, Madrid University, Madrid, Spain). The pRL-SV40 reporter vector containing the Renilla luciferase gene was used as internal control in the transactivation assay (Promega).

Transactivation Assays

Assessment of variant effects on transcriptional activity (TA) were performed in the HeLa cell line, representing cells negative for endogenous HNF-1A expression, and in the INS-1 (rat insulinoma cell line), representing cells positive for HNF-1A expression. HeLa and INS-1 cells were grown as previously described. Transient transfection of variant plasmids (wild-type or variant HNF1A), reporter, and control plasmids was performed using Lipofectamine 2000 (Life Technology), as reported. Luciferase activity was measured 24 h post-transfection with the Dual-Luciferase Assay System (Promega) in a Chameleon luminometer (Hidex) or using the Enspire platform (Perkin Elmer). Luciferase activity was normalized for transfection efficiency by the Renilla luciferase activity. The Bergen dataset included some variants previously reported (p.Ile27Leu, c.92G>A [p.Gly31Asp], c.142G>A [p.Glu48Lys], c.290C>T [p.Ala97Val], c.293C>T [p.Ala98Val], c.298C>A [p.Gln100Lys], c.341G>A [p.Arg114His], c.392G>A [p.Arg131Gln], c.965A>G [p.Tyr322Cys], c.1165T>G [p.Leu389Val], c.1405C>T [p.His469Tyr], c.1460G>A [p.Ser487Asn], c.1469T>C [p.Met490Thr], c.1541A>G [p.His514Arg], c.1544C>A [p.Thr515Lys], and c.1729C>G [p.His577Asp]).

HNF-1A Protein Abundance

The level of wild-type and individual HNF-1A variant protein expressions in total HeLa cell lysates was determined. For this purpose, the team at Bergen used 20 μL of HeLa cell lysates generated for transactivation assays as previously described. HNF-1A and actin protein levels were quantified by densitometric analysis using Quantity One 1-D software (Bio-Rad). The team at Oxford evaluated protein expression by transfecting HeLa cells with 5 μg of wild-type or variants plasmids after culturing for 24 h. Total protein quantification was carried out using Bradford reagent (Bio-Rad) and 10 μg of total protein was electrophoresed then immunoblotted with antibodies for HNF-1A (Santa Cruz Biotechnology) and beta-tubulin (Santa Cruz Biotechnology) and visualized using the ChemiDoc Imaging System (Bio-Rad). Densitometry (for western blots and EMSA) was carried out using Image Lab Software (Bio-Rad).

Subcellular Localization

For nuclear translocation assessments, the teams at Bergen and Oxford examined HNF-1A presence in nuclear versus cytosolic HeLa cell fractions. Cultured and plated cells were transiently transfected with wild-type or HNF1A variant plasmids. Sequential cell fractionation from each transfected sample was performed 24 h post-transfection as described. 20 μg total protein from each isolated compartment (nucleus and cytosol) was analyzed by SDS-PAGE and immunoblotting using an HNF-1A-specific antibody (Cell Signaling or Santa Cruz Biotechnology, respective to the two centers). GAPDH antibody (Santa Cruz Biotechnology) and Topoisomerase II-alpha antibody (Cell) were used as loading control for cytosol and nuclear compartments, respectively. The HNF1A variant c.589_615del (p.Leu197_Leu205del), denoted p.delB, was included as a positive control for impaired nuclear localization (cytosolic retention).

DNA Binding

DNA binding ability test was conducted for HNF1A variants that were located in DNA binding domain (1–287 aa), and those that demonstrated transactivation activity < 50%. At Bergen the HNF-1A proteins were expressed using an in vitro transcription and translation system (TNT Coupled Reticulocyte Lysate System, Promega) and equal amounts of synthesized protein were bound to a [γ-32P]-radiolabeled rat albumin oligonucleotide as described. DNA-protein bound complexes were separated by 6% DNA retardation gel electrophoresis (EMSA) (Life Technologies) followed by autoradiography (LAS-1000 Plus, Fujifilm Medical System). Level of DNA binding was assessed by quantification of the intensity of HNF-1A protein-oligonucleotide complexes by the program Image Gauge 3.12 (Fujifilm Medical Systems). The two HNF1A-MODY control variants c.335C>T (p.Pro112Leu) and c.608G>A (p.Arg203His) were included as positive controls for reduced DNA binding ability. DNA binding data for two variants p.Gln100Lys and p.Arg131Gln are the same as published previously. The team at Oxford used 10 μg of the cell lysate used in the protein abundance western and a CY5 labeled probe of the same promoter sequence as used by the Bergen team. However, they used the Odessey Infra-Red EMSA kit (Li-Cor Inc.) to conduct binding affinity. The percentage of HNF-1A-oligo complexes for each variant fraction was then calculated compared to wild-type.

Stratification of HNF1A Variants with Unsupervised Learning Methods

We designed an analytical pipeline to stratify functionally characterized HNF1A variants along the spectrum of glycemic phenotypes (Figure S3). Briefly, the pipeline begins with preparation of the dataset for analysis using unsupervised learning tools. This was followed by the addition of two scores to each variant, one to account for functional domain and the other for spatial variation in HNF1A isoform expression (exon location), as there are well-established correlations between variant position in HNF1A and clinical phenotype., The “polished” dataset was then processed using principal component analysis, k-means, and hierarchical clustering methods (Figure S3).The prcomp (principal components analysis) function in R (stats package v.3.5.0) was used to perform principal component analysis. Polished data matrices were zero centered and scaled to account for unit variance. We used the NbClust Package in R for determining the best number of clusters (distance measure set as “Euclidean”) to estimate the optimal number of k-means clusters in PC space (Figures S4A and S4B). Hierarchical clustering was performed on a Euclidean distance matrix comprised of PC scores for each allele from total principal components that explained >85% of data variance (Figures S4C and S4D) using the hclust (hierarchical clustering) function in R. The WARD minimum variance hierarchical clustering method (ward.d2) was selected as it yielded highest resolution clusters in a comparative analysis against complete, single, and average linkage methods (data not shown) based on (1) predicted grouping patterns of wild-type and MODY/type 2 diabetes risk reference variants and (2) known molecular function of variants which co-occupied clusters defined by wild-type, type 2 diabetes risk, and MODY reference variants. To compare the cluster dendrograms of variants shared between Oxford and Bergen, we used the untangle, tanglegram, and entanglement functions in R (part of the dendextend package) to untangle dendrogram lists and find the best alignment layout, plot the two dendrograms side by side, and compute the quality of alignment (entanglement coefficient), respectively.

Dataset Preparation for Clustering Analysis

The functional datasets were prepared for PCA and clustering analysis by harmonizing the number of variables across tested variants. DNA binding ability was interrogated for only a small subset of variants, so EMSA data were excluded. Functional data were available in three different formats: raw instrument data, data normalized to internal assay controls (renilla luciferase for TA assay, beta-tubulin/actin antibody for protein abundance assay, and nuclear:cytosolic ratio of raw protein abundance reads for nuclear localization) expressed as biological replicates, and fully processed summary data normalized to wild-type values. The most statistically suitable input format for PCA and unsupervised clustering is functional data normalized to internal assay controls (semi-processed) as intra-assay measurements are harmonized (versus raw instrument data) and the organic structure of the data is retained and uninfluenced by assumptions (versus wild-type normalized data). Further, this format yielded the most robust clustering trends based on distribution quality in multivariate space and known and expected sequence-function relationships. Scores for tissue-specific expression of HNF1A isoforms (implications for clinical phenotypic manifestation) and functional domain (varied levels of mutation tolerance) were assigned to each variant (Table S3). For stratification of HNF1A variants with unsupervised learning methods, see Figure S4 and Material and Methods.

Variant-Phenotype Mapping

We surveyed the UK MODY Diagnostic Registry (Royal Devon and Exeter NHS Foundation Trust, Exeter, UK) and the Norwegian MODY Registry (Haukeland University Hospital, Bergen, Norway) for functionally annotated HNF1A missense variants. A total of 162 and 53 HNF1A missense variants were documented in the UK and Norwegian diagnostic registries, respectively. Tables S4 and S5 show the list of clinical features that were available from the database for alleles which overlapped with the Oxford-Bergen dataset (not all features available for each variant). Sequence variants in the Norwegian MODY Registry were classified prior to the incorporation of the ACMG/AMP guidelines, as described using a 5-tier score system. Sequence variants in the Exeter MODY Registry had been classified using the ACGS guidelines from 2013 (see Web Resources), a 5-tier system used in the UK prior to the advent of the ExAC database and publication of the ACMG/AMP guidelines. Original clinical reports of carriers were accessed for additional details, particularly where clinical features were sparse—such as extra-pancreatic features, vascular complications, additional family history data, and whether, for example, other MODY genes were next-generation sequenced as part of a MODY gene panel. The classification system adopted by each center was used for reclassification of variants from its database.

Role of the Funding Source

The funders of the study had no role in study design, data collection, data analysis, data interpretation, or writing of the report. The corresponding author had full access to all the data in the study and had final responsibility for the decision to submit for publication.

Results

In Silico and In Vitro Functional Characterization of Variants

To resolve HNF1A genotype-phenotype complexity, we sought to evaluate the function of 73 HNF1A missense alleles which were observed almost 26K times in the exomes of ∼13K multi-ethnic type 2 diabetes case subjects and control subjects. The majority of HNF1A variants were identified in both type 2 diabetes case subjects and control subjects, and for the few observed exclusively in type 2 diabetes case subjects, they were identified with a frequency of either one or two case subjects per variant (Table S1). Although there was no evidence for HNF1A association with type 2 diabetes susceptibility in the study either at the single variant or at the gene level, there was a marginal aggregate association with type 2 diabetes risk in a multigene test (p = 0.023) which included rare coding alleles in a set of genes implicated in monogenic/syndromic diabetes or related glycemic traits. Consensus across multiple in silico tools for predicting pathogenicity was observed in only 38 of the 73 variants (53%) (Table S1). Based on CADD scores (suggested pathogenicity cut-off > 15), ∼70% of the missense variants would be bioinformatically classified as disease causing (i.e., sufficient to cause MODY).

The 73 HNF1A variants were divided between the two centers (Oxford and Bergen) and individually evaluated in terms of functional effect using a common pipeline including assays measuring variant effect on HNF-1A transcriptional activity, subcellular localization, protein expression level, and DNA binding ability (Figures S1–S3). The Oxford laboratory investigated variants predominantly of South-East Asian etiology, while the Bergen laboratory studied variants mainly of European etiology.

The variants demonstrated a wide range of functional effects from benign to damaging across assays and laboratories with highest variability in transactivation assessments, particularly through regulation of the rat albumin promoter in HeLa cells (activity range 30%–110% in Bergen data [Figure S5A], 52%–114% in Oxford data [Figure S6A]). Transcriptional activity was consistently higher for variants using HNF4A P2 promoter in INS-1 cells (versus rat albumin promoter in HeLa cells) (activity range 40%–90% Bergen data, 77%–158% Oxford data), and most likely due to interference of endogenous HNF-1A in INS-1 cells (2- to 4-fold higher basal promoter activity, Figures S5B and S6B). In assessments of protein abundance, >85% of all variants displayed adequate HNF-1A protein levels (>60%, Figures S5C and S6C). Similarly, in nuclear translocation assays, most variants were predominantly detected in the nucleus (level > 60%), the exception being five variants from the Bergen dataset (c.827C>A [p.Ala276Asp], p.Ser487Asn, c.1812C>G [p.Ser604Arg], c.1322C>A [p.Thr441Lys], p.Arg131Gln) (Figures S5D and S6D) from the Oxford dataset (c.185A>G [p.Asn62Ser], c.1552C>T [p.Leu518Phe], c.1605C>A [p.Ser535Arg], c.1610C>T [p.Thr537Met], c.1748C>A [p.Arg583Gln]) displaying < 50% level. The subsets of variants investigated by EMSA demonstrated overall normal DNA binding ability (∼95% of variants > 60%), with the exception of one variant from the Bergen dataset (p.Arg131Gln < 50%) and five from Oxford dataset (p.Asn62Ser, c.340C>T [p.Arg114Cys], c.467C>T [p.Thr156Met], c.481G>A [p.Ala161Thr], p.Ser535Arg < 50%) (Figures S5E and S6E).

Multi-Dimensional Data Analysis

We performed unsupervised stratification of the HNF1A variants using the multi-dimensional in vitro functional data (semi-processed data normalized to internal technical controls in each assay) supplemented with scores for isoform expression and functional domain, with the aim of mapping molecular dysfunction to clinical phenotype. Each of the two datasets were analyzed independently to minimize the interference of inter-laboratory variability with true biological signal. We used principal component analysis to facilitate informative dissection and visualization of multi-parametric functional data. Eigendecomposition of the data matrices revealed transcriptional activity as the greatest contributor to data variance and structure (Figure 1). To enable variant subgroup discovery for function-phenotype mapping, data were partitioned using (1) k-means clustering in PC space (Figure 2) and (2) hierarchical clustering using data coordinates from the total number of informative principal components for each dataset (Figure 3). The analysis yielded variant clusters neatly organized along the spectrum of HNF1A molecular dysfunction ranging from neutral/benign to intermediate to damaging. As such, we broadly annotated the known HNF1A in vivo spectrum, from benign to type 2 diabetes risk-modifying to HNF1A-MODY (inherited early-onset hyperglycemia, likely to be sulfonylurea-responsive based on MODY registry data), onto the principal components plots and dendrograms based on the spatial distribution of the variants along the in vitro data-derived functional spectrum, from wild-type/wild-type-like to intermediate to damaging (Figures 2 and 3). To understand (mis)alignment of HNF1A variants shared between the two centers, we visually compared the cluster dendrograms of shared variants only (Figure S7). We computed an entanglement coefficient (the quality of the alignment of the two dendrograms expressed as a value from 0 to 1 where lower values correspond to higher quality alignment) of 0.055 which indicated good high-quality alignment of shared variants (Figure S7).

Figure 1

Eigendecomposition of Principal Components Explaining > 85% of Variance

Shown are (A) Oxford and (B) Bergen datasets. TA INS1_P2 and TA HeLa_ALB are transcriptional activity data from INS-1 cells using HNF4A P2 promoter and from HeLa cells using rat albumin promoter, respectively. PE, protein expression; nuc loc, nuclear localization data.

Figure 2

K-Means Clustering

HNF1A missense alleles characterized at Oxford (A) and Bergen (B) in principal component (PC) space. Blue and green k clusters represent alleles with benign and benign-to-intermediate effects on function, respectively; purple k clusters represent alleles with intermediate functional impact; red k clusters indicate intermediate-to-damaging or functionally damaging alleles.

Figure 3

Hierarchical Clustering Analysis

HNF1A missense alleles characterized at Oxford (A) and Bergen (B). WARD minimum variance method was used and analysis performed using orthogonally transformed functional data from PC1-PC4 (>85% explained variance) from Oxford dataset and PC1-PC5 (>85% explained variance) form Bergen dataset. To optimize visualization of the function phenotype gradient, some branches were rotated. The numbers of the y axes of (A) and (B) refer to clustering height calculated as by Ward's criterion (total within-cluster variance).

Clinical Interpretation of HNF1A Variants

To assess the medical diagnostic utility of multi-tiered HNF1A sequence-function annotations, we examined their mapability to HNF1A clinical phenotype using clinical data from overlapping HNF1A missense variants in the UK and Norway MODY diagnostic registries (Tables S4 and S5).

Of the 31 total overlapping variants between our functional effort and the UK registry, 19 were originally classified as pathogenic/likely pathogenic and 15 as VUS/likely benign in the diagnostic registry. Three of the 31 overlapping variants (c.1816G>A [p.Gly606Ser], p.His469Tyr, c.871C>T [p.Pro291Ser]) were present under both pathogenic/likely pathogenic (where they were considered the MODY-causal variant in the case subjects) and VUS/likely benign (cases of co-occurrence with a pathogenic variant in HNF1A or another MODY gene) original clinical classifications. All 15 missense variants categorized as VUS/likely benign in the UK database demonstrated benign clustering patterns in our analysis (i.e., did not form subgroups with variants which exhibited impaired function). However, for 10 of the 19 variants clinically categorized as pathogenic/likely pathogenic in the UK diagnostic database (p.Ala161Thr, c.521C>T [p.Ala174Val], c.139G>C [p.Gly47Arg], p.Gly606Ser, p.His469Tyr, c.1235T>C [p.Met412Thr], p.Asn62Ser, p.Pro291Ser, p.Arg131Gln, c.29C>T [p.Thr10Met]), patterns of in vitro functional clustering patterns did not match clinical diagnostic variant interpretation. The variants co-occupied clusters either with known type 2 diabetes risk modifiers (some moderately impacted in functional assays) or with wild-type/neutral variants. Discordance between functional genotype and clinical variant interpretation prompted a thorough reassessment of variant pathogenicity.

The missense variant p.Asn62Ser (gnomAD allele count n = 33) was consistently dissimilar to dysfunctional variants in dendrograms and k-means derived clusters from both Bergen and Oxford datasets (Figures 2 and 3). In the UK MODY registry, it was identified in an obese individual who was diagnosed with diabetes at age 36 years (Table S6). The patient suffered from microvascular complications (MIM: 603933) (nephropathy and retinopathy) (Table S6). These features are inconsistent with neither HNF1A-MODY nor a type 2 diabetes phenotype (which might be familial considering the number of affected individuals in the carrier’s pedigree). HNF1A was the only MODY gene sequenced in this individual as genetic testing was performed before the advent of the targeted MODY gene exome sequencing panel which is the current diagnostic procedure. Alone, p.Asn62Ser allele frequency values are sufficient to confidently re-categorize the variant as VUS/likely benign in the context of MODY (Figure 4).

Figure 4

Distribution of Functionally Annotated HNF1A Missense Alleles

(A and B) As a function of frequency in the (A) UK MODY diagnostic registry and (B) Norway MODY diagnostic registry on the x axis and reported frequency in the genome aggregation database (gnomAD) on the y axis. Alleles are colored on the basis of the (re)classification scheme on the top right.

(C) Frequency of functionally characterized exome-detected HNF1A missense alleles in gnomAD. The red and orange dashed lines mark known ultra-rare, MODY pathogenic (allele count ≤ 2, AF < 0.0008%) and low frequency type 2 diabetes predisposing allele frequencies (allele count ≤ 121, AF < 0.04%) respectively.

The variants p.Ala174Val and p.Pro291Ser, characterized by both laboratories, were more difficult to interpret. In the Oxford dataset, these variants formed a separate outlying k-means-derived cluster (Figure 2). They also occupied an independent subgroup in hierarchical clustering which branches high on the height scale away from the larger cluster defined by wild-type and other neutral variants (Figure 3). Atypically high luciferase renilla values (internal luciferase reporter gene assay control used for normalization) were reported for these variants as well as for p.Ser535Arg, which have resulted in a potentially exaggerated reduction in transactivation values for these variants upon normalization to the internal assay reference in the Oxford dataset (Figure S6). In the Bergen dataset, these variants also consistently lie in the type 2 diabetes risk modifier zone (not pathogenic for MODY) (Figures S2 and S4). Not only are the activity profiles of p.Ala174Val and p.Pro291Ser dissimilar to those of pathogenic MODY variants, they also occur at a much higher frequency in the general population (Figure 4). In terms of clinical profiles, p.Ala174Val was detected in an individual with diet-controlled diabetes, which was diagnosed at age 25 years (Table S6). The p.Pro291Ser variant was detected in an overweight individual who was diagnosed with diabetes at age 42 years when it was classified as likely pathogenic/pathogenic (Table S6). In another unrelated individual, p.Pro291Ser was co-expressed with p.Gly31Asp and both were annotated as VUS/likely benign in the diagnostic database.

The clinical diagnostic classifications of p.Gly606Ser and p.Ala161Thr did not match clustering patterns in multivariate space (Figures 2 and 3). The variants did not impact HNF-1A function in the in vitro assays tested. The frequency associated with these variants (n = 12 alleles in gnomAD and n = 2 in the Exeter diagnostic clinic) are inconsistent with those of rare MODY-causing variants. The p.Gly606Ser variant has also been found in a single case of hyperinsulinemic hypoglycemia (on diazoxide treatment; MIM: 256450) in the UK registry and in this case was classified as VUS/likely benign.

In clustering analysis, p.His469Tyr occupied either the same or highly similar (adjacent) subgroups as wild-type (Leu27 and Ile27) (Figure 3). The clinical features of the variant carrier described in Table S6 are consistent with severe young-onset familial diabetes; however, the allele high frequency in gnomAD is 32. Another variant in HNF1A, c.620G>A (p.Gly207Asp) (not present in gnomAD), was, however, detected in the same individual. It was identified in three other case subjects (including co-occurrence with p.His469Tyr) and was classified as pathogenic/likely pathogenic each time it was identified in the UK MODY Registry. Thus, it is likely that p.Gly207Asp is the MODY-causal variant and that p.His469Tyr is either benign or potentially type 2 diabetes risk-modifying.

Despite alignment between clustering pattern and clinical diagnostic interpretation (Figures 2 and 3), c.1576G>A (p.Asp526Asn) was reclassified from pathogenic/likely pathogenic to VUS/likely benign. In transactivation assays, HNF1A-p.Asp526Asn was the most impaired of all tested exome-identified variants in the Oxford dataset (∼50% in HeLa and ∼80% in INS-1 cells; MODY reference variants exhibited transactivation range of 20%–40% in HeLa cells and 30%–50% in INS-1 cells in Oxford, with the exception of p.Arg203His and c.347C>T [p.Ala116Val], which yielded transactivation values of 40% and ∼60% in HeLa and ∼50 and ∼100% in INS-1 assays, respectively). The variant is observed only five times in gnomAD, which suggests it might not be causal for MODY (Figure 4). The clinical profile of the p.Asp526Asn carrier did not appear to be consistent with HNF1A-MODY, besides presence of diabetes in three generations of the carrier’s family. The variant carrier had BMI 32.4 kg/m2 and was diagnosed with diabetes at age 33 years. Other clinical features included dyslipidemia (MIM: 144250), polycystic ovary syndrome (MIM: 184700), insulin resistance (MIM: 610549), and hypertension (MIM: 145500). The variant was also found in a patient in the Norwegian MODY registry. This patient, diagnosed at 19 years of age, had normal BMI and C-peptide levels. Type 1 diabetes (MIM: 222100) autoantibody status and type 1 diabetes risk score were not known. The carrier was treated with metformin. His mother and the mother’s brother also have diabetes (treated with diet and insulin, respectively). Moreover, the patient was diagnosed with Crohn disease (MIM: 266600). Altogether, this suggests that the carriers might have a combination of type 2 diabetes and HNF1A-MODY, which is not uncommon, or a phenotype representing a possible continuum of diabetes sub-phenotypes from MODY to type 2 diabetes. Further, the variant is expressed in the hepatocyte-dominant isoform and is thus unlikely to manifest in a strong beta-cell phenotype despite its poor functionality.

Of the 19 HNF1A missense variants that overlapped with the Norwegian MODY Registry, 18 were originally classified as benign (class 1), likely benign (class 2), or VUS (class 3), and 1 (p.Ala276Asp) as likely pathogenic (class 4) in the diagnostic registry. The variant p.Ala276Asp consistently demonstrated impaired HNF-1A function in in vitro assays and clustered with the MODY reference variants in the unsupervised clustering analyses, supporting the clinical interpretation of this variant as pathogenic (Figures 2 and 3). It was also clinically classified as likely pathogenic/pathogenic in the UK MODY registry (Figure 4). All variants originally classified as benign/likely benign/VUS (class 1–3) in the Norwegian registry clustered in the benign or intermediate type 2 diabetes risk modifier zones, with the exception of four variants (c.1016C>T [p.Thr354Met], p.Thr441Lys, c.1745A>G [p.His582Arg], c.1756G>A [p.Ala586Thr]), which demonstrated variable trends across clustering methods (Figures 2 and 3).

In k-means clustering along principal component 1 and principal component 2, these variants co-occupy a hard cluster with MODY reference variants and variants which exhibited damaging in vitro function. This is not entirely unexpected for HNF1A-p.Thr441Lys which displayed reduced activity (∼50% on both promoters in INS-1 and HeLa cells) and with reduced (<40%) nuclear localization. Although, in hierarchical clustering, where (dis)similarity between variants was determined using principal component scores from all principal components contributing to >85% of overall variance, the trends were more consistent with clinical features and classification; p.Thr441Lys and p.Thr354Met are in the type 2 diabetes risk modifier space of the in vivo continuum, hierarchically distanced from the sub-cluster defined by the majority of MODY reference variants and pathogenic damaging variants p.Ala276Asp and c.1135C>A (p.Pro379Thr). The clinical phenotypic data of the p.Thr354Met variant carriers seems more consistent with type 1 diabetes. The p.Thr354Met variant was identified in two unrelated individuals. Upon revisiting clinical data on these two allele carriers, it was found that one of the carriers with insulin-treated diabetes from age 14 years was positive for GAD and IA2 autoantibodies. A sister of the proband had diabetes, but the parents were apparently unaffected. The other p.Thr354Met variant carrier had autoantibody-negative diabetes from age 12 years without strong family history of diabetes (grandmother only). Moreover, the population frequency (n = 18 in gnomAD) associated with this variant allele is inconsistent with rare, causal MODY alleles. Thus, the clinical phenotypic data of the p.Thr354Met variant carriers seem more consistent with type 1 diabetes. In the p.Thr441Lys variant carrier, another variant in HNF1A c.872dup (p.Gly292Argfs∗25) was considered the pathogenic MODY variant (Table S7). Moreover, the population frequency values of p.Thr441Lys (gnomAD allele count n = 18) and p.Thr354Met (gnomAD allele count n = 7) are slightly higher than expected for rare disease-causing variants (Figure 4). As for p.His582Arg (gnomAD allele count n = 14) and p.Ala586Thr (gnomAD allele count n = 20), in hierarchical clustering, these variants form a subgroup defined by liver isoform variants which demonstrated suboptimal function in one or more in vitro assays. Much like p.Asp526Asn, these variants are likely to be strong type 2 diabetes risk modifiers. The p.His582Arg variant carrier was diagnosed with diabetes age 11 years, she had a BMI of 29 at referral one year later, and C-peptide was measured to 1,000 pmol/L (Table S7). The p.Ala586Thr variant carrier was diagnosed at age 11, C-peptide positive (78 pmol/L), negative GADA, IA2A, ZnT8A, with no known family history of diabetes, and treated with insulin (Table S7).

Based on this comprehensive variant re-assessment effort, we changed the classification of 7 out of 31 variants shared with the UK MODY diagnostic database (p.Ala161Thr, p.Ala174Val, p.Gly606Ser, p.His469Tyr, p.Asn62Ser, Pro291Ser, p.Asp526Asn) from likely pathogenic to VUS/likely benign (Figure 4) and all five variants categorized as VUS (class 3) or VUS/likely benign (class 3−) in the Norway MODY registry (p.Tyr322Cys, p.Thr354Met, p.Thr441Lys, p.Asp526Asn, and p.His582Arg) to likely benign (class 2) (Figure 4). This represents ∼23% and ∼26% of total HNF1A missense variants in the UK and Norway MODY registries, respectively, that overlap with the functionally interrogated HNF1A missense variants detected in the exomes of ∼13K multi-ethnic type 2 diabetes cases and controls.

Discussion

In this study, we investigated the functional impact of 73 missense variants in HNF1A, detected by exome sequencing of a multi-ethnic type 2 diabetes case-control cohort from four different mechanistic angles (Figures S1 and S3). We developed an approach for the analysis of multi-parametric functional data which, in the context of HNF1A, has enabled (1) a holistic assessment of variant behavior by combining as many mechanistic dimensions as possible, (2) unbiased stratification along the spectrum of glycemic phenotypes ranging from neutral/benign effects, to modification of multifactorial polygenic diabetes risk, to deleterious and causal for early-onset sulfonylurea-responsive diabetes, (3) an assessment of the relative contributions of each functional parameter to molecular variability, and (4) rigorous phenotype mapping and a thorough re-evaluation of the clinical classifications of overlapping variants in two national MODY diagnostic registries.

Revisiting clinical variant classifications using HNF1A functional clusters led to the reclassification of ∼4% (7/162) and ∼9% (5/53) of all HNF1A missense variants in the UK and Norwegian MODY diagnostic registries, respectively. Decisions on variant reclassification were primarily motivated by the juxtaposition of allele frequency values in the general population (based on gnomAD allele counts) against their frequency in the MODY diagnostic registries (highest frequency values belonged to bona fide loss-of-function alleles used as MODY reference controls in this study) (Figure 4). This is based on the rules given by the ACMG guidelines for variant interpretation: “an allele frequency in a control population that is greater than expected for the disorder is considered strong support for a benign interpretation for a rare Mendelian disorder” (BS1). Information from other layers of variant annotation such as in vitro function (in the tested assays), clinical features, family history, ethnicity, and in silico prediction all helped to support reclassification decisions.

Dissection of the individual principal components revealed transactivation to be the primary contributor to the spatial distribution of multi-parametric data. This suggests that it might be a superior functional readout and potentially more informative than other molecular assays for assessing HNF1A variant pathogenicity, and in line with our previous experience on various functional assays of HNF1A variants. Since transactivation is a relatively all-encompassing measure of transcription factor protein function, we assumed defects in transcriptional activity would capture defects in its biochemical prerequisites (protein expression, nuclear transport, DNA binding). While this may have been the case for the majority of functionally interrogated HNF1A variants, p.Arg203His and p.Ala116Val highlight the limitations of this assumption. For these variants, severely impaired DNA binding ability was not adequately captured by transactivation. This might explain clustering of p.Ala116Val and p.Arg203His in the intermediate zones among type 2 diabetes risk modifiers, and not directly among MODY reference alleles, in both k-means and hierarchical clustering for both centers (Figure 3).

We were able to mitigate error associated with handling data from two centers with methodological differences by benchmarking several HNF-1A variants (benign, type 2 diabetes risk, and MODY) in both laboratories. Discrepancies between the two centers with respect to shared variants can be explained in part by variability in technical protocols between laboratories and the handling of samples by various individuals over the course of the study. The relative clustering position of the shared variants is impacted by the function trends observed in each dataset; while the majority of Oxford variants behaved wild-type-like, the Bergen dataset was more complex as variants demonstrated a wider range of effects. For instance, for an intermediate variant and a known type 2 diabetes risk modifier such as p.Glu508Lys, the dissimilarity to MODY reference alleles is more pronounced in the Bergen dataset where there are more data points between moderately impaired and damaging function (Figure S7).

While important, informative, and powerful first lines of evidence, functional annotations should not be treated as superior or stand-alone determinants of variant pathogenicity, which they have not in this study. Indeed, the same variant in a MODY gene can give rise to a spectrum of clinical phenotypes and exhibit variable penetrance depending on genomic (regulatory variants in cis or trans, or haplotype epistasis) and environmental (epigenomic) context which are difficult to capture in functional assessments.29, 30, 31, 32, 33 It is also entirely possible that some of the noise in functional-clinical mapping is a reflection of the heterogeneity in the clinical phenotypic manifestation of HNF1A variants. The same variant which has a mild effect on HNF1A function, and thus beta-cell function and ability to respond appropriately to a given level of glycemia, could play out differently in individuals who are already struggling to meet the insulin demand through insulin resistance and/or other genetically driven defects in their beta-cells.Another aspect to consider is the expected variation in clinical practice between the two centers in the UK and Norway to which diabetes patients have been referred. It would thus be naive to attempt to draw conclusions regarding variant effects in vivo from, for example, a single registry observation. The reality of phenotypic variant manifestation is often complex, context dependent, non-linear, and spectrum based. Developing a contextual and thorough understanding of variant behavior from diverse functional, clinical, biochemical, and demographic datasets is necessary to facilitate highly accurate interpretation.

The p.Asp526Asn variant in HNF1A is a perfect example that illustrates the need for nuanced evaluations of variant effects despite the availability of multiple layers of functional annotation from various cell systems. The variant was clinically classified as pathogenic/likely pathogenic and exhibited impaired in vitro functional activity (shared a cluster with known MODY-causal variants). Its impact on molecular function would be consistent with biomarker profiles (hsCRP and glycans) suggestive of HNF1A-MODY. Yet, re-evaluation of the clinical features of variant carriers in UK and Norway diabetes registries and the fact that it is present only in the longest HNF1A transcript isoform expressed predominantly in liver suggest that it is more likely to be a contributing factor to common multifactorial diabetes rather than a primary driver of early-onset sulfonylurea-responsive familial hyperglycemia. Integration of isoform weights into the unsupervised clustering model helped separate bona fide loss-of-function MODY variants from functional variants expressed in exons 8–10 at the lowest level of hierarchical clustering.

The recent advent of multiplexed assays of variant effects (MAVEs) has made it possible to interrogate the function of every possible sequence perturbation in a single experimental system. Successful and productive implementation of these technologies requires overcoming the technical and analytical complexities associated with scaling up, which represent the most significant barrier in the face of closing the chasm between variant resolution and variant interpretation. A meticulously designed MAVE for HNF1A would enable functional annotation of all possible missense variants (>12,000) in a single assay. A high-performance variant classifier built using MAVE-based data can then be used to generate an exhaustive catalog of variant effects which researchers and clinicians can consult upon sequence-identification of an HNF1A variant.

The performance and predictive utility of any model built using HNF1A MAVE-derived function scores would be enhanced immeasurably upon calibration against these multi-layered data and comprehensively annotated function-clusters. These data can also improve existing prediction algorithms which operate on the basis of multi-factorial probability and multi-data integration such as CADD (Combined Annotation-Dependent Depletion), MutationTaster, FitCons (fitness consequence), and VAAST (Variant Annotation, Analysis and Search Tool). Further, they can be incorporated into rigorous and collaborative multi-level annotation efforts led by the Clinical Genome Resource (ClinGen) program and evidence-based disease-specific variant classification databases. Lastly, our approach can assist in filtering HNF1A missense variants for gene burden testing of rare variants. An immediate example is its application to the ∼50K exomes in the UK Biobank not ascertained on the basis of diabetes. At present, in the UK Biobank dataset, the number of carriers of HNF1A variants that overlap with variants functionally investigated in our effort are insufficient to conduct a robust gene burden analysis. However, we have observed that seven variants present in the UK biobank and predicted by our functional data to be likely damaging (p.Asp526Asn, p.Arg171Gln, p.Thr354Met, p.Gly288Arg, p.His582Arg, p.Pro379Thr, p.Ser592Pro) were not identified in patients with diabetes; it is possible that these alleles represent very rare pathogenic variants causing MODY with reduced penetrance or are type 2 diabetes risk variants.

In conclusion, we have developed an analytical framework for robust and unbiased variant stratification using multi-dimensional functional follow-up data from the largest number of exome-identified missense variants in HNF1A ever studied. This allowed us to annotate functional clusters with clinical knowledge and identify discordant classifications between functional genotype and clinical phenotype. We believe our pipeline is an important proof-of-principle technical contribution on the path toward more reliable, scalable, and comprehensive mapping of sequence-function relationships: a significant factor in making well-informed initial judgements of allele pathogenicity in the context of individual phenotypic presentations.

Declaration of Interests

The views expressed in this article are those of the authors and not necessarily those of the NHS, the NIHR, or the Department of Health. M.I.M. has served on advisory panels for Pfizer, NovoNordisk, and Zoe Global, received honoraria from Merck, Pfizer, NovoNordisk, and Eli Lilly, and received research funding from Abbvie, Astra Zeneca, Boehringer Ingelheim, Eli Lilly, Janssen, Merck, NovoNordisk, Pfizer, Roche, Sanofi Aventis, Servier, and Takeda. As of June 2019, M.I.M. is an employee of Genentech and a holder of Roche stock. A.L.G. reports grants from Wellcome Trust, grants from NIHR Oxford Biomedical Research Centre, grants from Horizon 2020, grants from NIDDK, and grants from MRC during the conduct of the study and personal fees from NovoNordisk and Merck outside the submitted work.