Systematic Functional Characterization of Candidate Causal Genes for Type 2 Diabetes Risk Variants

Type 2 diabetes risk is determined by a complex interplay between environmental and genetic factors, with heritability estimates ranging from 20% to 80% (1). Over the past decade, genome-wide association studies (GWAS) of ever-increasing size have discovered more than 100 regions of the genome (loci) associated with type 2 diabetes risk (2). Studies in individuals with diabetes have demonstrated that a large number of these association signals exert their effects on disease susceptibility through pancreatic islet dysfunction (3).

Despite these advances, progress in translating genetic findings into disease biology has been relatively slow. The majority of risk variants are located in noncoding regions of the genome and pinpointing the underlying causal genes or “effector transcripts” has proved challenging (4). Recent efforts have focused on identifying structural or functional links between association signals and regional genes (5,6). A complementary strategy uses candidate gene biology to prioritize genes located near association signals. High-throughput screening could facilitate the identification of genes implicated in β-cell function and thereby highlight potential effector transcripts at type 2 diabetes GWAS loci. To date, such approaches have been limited by the inadequacies of available human cellular models and the high cost of insulin immunoassays (∼$2 per data point), the gold standard for measuring insulin. To circumvent these issues, previous studies have relied on rodent β-cell models and either used reporter assays as a proxy for insulin measurements or focused on cellular proliferation (7–11).

Recently, the first glucose-responsive human β-cell line, EndoC-βH1, was generated (12,13). The line is derived from fetal pancreatic buds matured in vivo and displays modest but robust induction of insulin secretion in response to glucose and secretagogues. Detailed characterizations have shown the cell line to be an authentic model system for studying stimulus-coupled secretion (14–16).

To accelerate the discovery of causal genes for type 2 diabetes, the current study performed and validated a genetic screen in the EndoC-βH1 cell line. We identified genes at half of the type 2 diabetes–associated loci studied (37/75) where small interfering RNA (siRNA)-mediated silencing resulted in β-cell dysfunction. This demonstrates the feasibility of performing systematic screening for insulin secretion in a human β-cell model, with implications for both high-throughput genetic and chemical compound screening. Our results can be integrated with existing lines of evidence to prioritize effector transcripts at GWAS loci and highlight potential roles for ARL15, ZMIZ1, and THADA in the regulation of insulin secretion.

The EndoC-βH1 cell line was cultured as previously described and grown to near confluency (12). RNA was then TRIzol extracted and sequenced at the Oxford Genomics Centre (Wellcome Trust Centre for Human Genetics, University of Oxford) (see Supplementary Fig. 4 for details). The raw sequencing data have been deposited at the European Nucleotide Archive (http://www.ebi.ac.uk/ena) under accession number PRJEB15283.

Cellular phenotypes were adapted for automated screening on a PerkinElmer Janus liquid handling workstation based on previously described assays (Supplementary Fig. 1A) (17). Briefly, 20,000 cells/well were reverse transfected in 96-well format at final siRNA concentrations of 25 nmol/L preincubated with 0.2 μL RNAiMAX in Opti-MEM. Custom libraries of siRNAs (ON-TARGETplus SMARTpools [Dharmacon] for the primary screen and Silencer Select [Thermo Fisher Scientific] for follow-up validation) were designed based on criteria described in Supplementary Table 2. In each case, nontargeting (NT) sequences based on the same chemistries were used as negative controls. Three days after transfection, cells were starved overnight in complete media containing 2.8 mmol/L glucose followed by 1-h starvation in 0 mmol/L media. Static insulin secretion assays were then performed for 1 h in complete media under the indicated conditions, after which cells were counted as described below.

Following secretion assays, supernatants were analyzed for insulin content using AlphaLISA Human Insulin Immunoassays (PerkinElmer) on a PHERAstar FS plate reader (BMG Labtech). Supernatants (50–250 nL) and beads prediluted in water (500 nL) were dispensed into 384–shallow well microplates with an Echo 550 (Labcyte) acoustic liquid handler before manual addition of immunoassay buffer to a final volume of 5 μL. Cell counts were measured using the CyQUANT Direct Cell Proliferation kit (Thermo Fisher Scientific) on an EnVision plate reader (PerkinElmer). All responses were normalized as indicated (see relevant figure legends) and expressed as a percentage of NT control for each phenotype. Effect sizes are given as the percentage difference from NT (Response_Gene − Response_NT) and the absolute values hereof (|Response_Gene − Response_NT|).

Data analysis was performed using R 3.0.2. To identify significant responses, cell counts and normalized insulin secretion measurements for each gene were compared with NT control using Student two-sample t test. The false discovery rate (FDR) was controlled at 5% by applying the Benjamini-Hochberg procedure to produce adjusted P values (q values) for each phenotype. The Z-factor measuring the control response for each phenotype was calculated as

We first developed an automated assay for disease-relevant phenotypes in the human β-cell line EndoC-βH1 (Supplementary Fig. 1A). Selected targets were silenced in a parallel format using RNA interference. Cells were then assessed for effects on cell number and insulin secretion under four different conditions: low glucose (1 mmol/L), high glucose (20 mmol/L), and high glucose with the sulfonylurea tolbutamide (100 μmol/L) or with the phosphodiesterase inhibitor IBMX (100 μmol/L). Low- and high-glucose conditions were included to provide information on the effect of gene silencing under conditions representing the fasted and fed states in vivo. Tolbutamide and IBMX act on the depolarizing and the potentiating pathways of insulin secretion, respectively, and were included to provide additional mechanistic insights through modulation (e.g., synergy or pharmacological rescue) of any primary defects observed in low or high glucose.

To reduce the cost of sample analysis, we made use of acoustic liquid handling to miniaturize insulin immunoassays. This generalizable method enabled us to maintain high sensitivity for insulin measurements (coefficient of variation <3%) (Supplementary Fig. 2A and B), while obtaining a 10-fold reduction in the cost of sample analysis ($0.20 per data point). Using the insulin gene (INS) as a positive control, we confirmed that we were able to robustly detect effects of gene silencing on the phenotypes of our assay (mean Z′ = 0.6 across conditions) (Supplementary Fig. 3).

On the basis of this combined analysis and assay pipeline, we designed a primary screen to assess the role of positional candidate genes for type 2 diabetes GWAS loci in β-cell function (Supplementary Fig. 1B). For target selection, we considered all protein-coding genes located within 1 Mb of a type 2 diabetes association signal. To exclude genes not expressed in our cellular model, we performed whole-genome RNA sequencing of the EndoC-βH1. Our expression data strongly correlated with published sequencing data for enriched primary β-cells (ρ = 0.78) (Supplementary Fig. 4) and showed robust expression of key β-cell genes (12,18) (Supplementary Table 1). We included only genes expressed in both EndoC-βH1 and primary β-cells (Supplementary Table 2), resulting in inclusion of 300 positional candidates from 75 type 2 diabetes GWAS loci.

We next performed our primary screen in triplicate and derived standardized scores for each phenotype. Knockdown was visibly confirmed using PLK1, an essential gene, which caused extensive cell death across all conditions. In a representative subset of 16 genes, we assessed knockdown efficiency at the transcript level and found the median residual expression to be 43% (Supplementary Fig. 5), roughly equivalent to monoallelic loss of function. To account for differences in plating efficiency and proliferation, cell counts were used to normalize insulin secretion data on a per-well basis. Two criteria were then applied to identify robust effects (“hits”): 1) an FDR-adjusted q value <0.05 and 2) an absolute effect size among the top 5% (Supplementary Fig. 6). This identified a total of 67 hits (15 for cell count and 52 for insulin secretion phenotypes) between 45 genes at 37 loci (Table 1).

For cell numbers, effect sizes for each gene were estimated based on 12 independently plated replicates (four conditions in triplicate) and therefore likely represent true differences in cellular proliferation and/or viability rather than random plating effects. Aside from KIF11, a gene with a known role in cell division, the largest effect sizes compared with NT control (coefficient of variation = 4% for cell numbers) were observed for ZMIZ1 (−15.2%; q value = 6.5 × 10⁻⁵) and PRDX3 (+16.6%; q value = 9.2 × 10⁻⁵).

For the insulin secretion data, we first performed an enrichment analysis for genes implicated in maturity-onset diabetes of the young (MODY). MODY describes a collection of monogenic subtypes of diabetes characterized by insufficient release or production of insulin. As would be expected for a set of bona fide regulators of β-cell function, we observed a strong enrichment of MODY genes among the significant hits (Fisher exact test, P = 5.5 × 10⁻⁹). Aggregating absolute effect sizes for MODY and non-MODY genes revealed this enrichment to be driven by altered insulin secretion and not by effects on cell numbers (Fig. 1).

Further validating our secretion data, we observed strong positive correlations between the normalized responses across conditions (P < 2.2 × 10⁻¹⁶) (Fig. 2) and found 10 of 35 genes to cause significant effects under two or more conditions. This included four known MODY genes and ZMIZ1, which, independently of the effect on cell numbers, was one of the strongest hits for reduced insulin secretion (q value <0.01 for low and high glucose). Knockdown of the ABCC8 gene, which encodes a subunit of the ATP-sensitive potassium channel, was found to significantly increase insulin secretion under low-glucose and IBMX stimulation. As expected, the depolarization caused by this was masked under high glucose (as cells are already partially depolarized) and fully rescued by tolbutamide (due to pharmacological depolarization of the cells). The pattern of modulation by secretion conditions can thus be used to pinpoint specific biological pathways affected by gene silencing. To explore the relationship between conditions in greater detail, we performed clustering analysis on Z-scores derived from the normalized secretion values. This revealed high glucose and tolbutamide to be most similar in terms of modulating knockdown effects, with low glucose and tolbutamide being most dissimilar (Supplementary Fig. 7).

Finally, we assessed the contribution of off-target effects by performing a small-scale validation experiment using siRNAs designed with an alternate algorithm. The sequences were confirmed to be different from those of the primary screen and could thus be used to establish the biological relevance of positive hits. We selected eight target genes, representing hits for both positive and negative defects across the four conditions, and confirmed that knockdown efficiency was satisfactory (median residual expression = 19.3%) (Supplementary Fig. 8). Compared with insulin secretion results from the primary screen, we observed an excellent linear correlation (ρ = 0.85, P = 6.7 × 10⁻¹⁰) (Supplementary Fig. 9) and 88% directional consistency in normalized responses. The validated hits included several genes with limited prior evidence of a role in the regulation of β-cell function, including ARL15 and ZMIZ1, which were found to significantly reduce insulin secretion across conditions (q values <0.05) (Fig. 3A and B), and THADA, which modestly elevated insulin secretion across three conditions, though the effect under low glucose was not observed in the primary screen (q value = 5.6 × 10⁻³) (Fig. 3C). Interestingly, gene silencing of the known MODY gene HNF4A was confirmed to cause a paradoxical increase in insulin secretion across all four conditions tested (q values <0.001) (Fig. 3D).

High-throughput screens for β-cell dysfunction offer the opportunity to systematically characterize the role of genes in a disease-relevant tissue for type 2 diabetes. Previous efforts have focused on nonhuman model systems (7–10), reporter-based proxy measurements for insulin (7,8), and/or phenotypes not directly related to insulin production and secretion (10,11). Here, we report a genetic screening strategy for the interrogation of multiple disease–relevant phenotypes in the human β-cell line EndoC-βH1. In a primary screen of 300 positional candidates, we successfully identified 15 genes regulating cell number (proliferation and/or viability) and 35 genes regulating insulin secretion. This is, to our knowledge, the first systematic, large-scale effort to identify genes involved in insulin secretion. Importantly, the identified hits can be used to prioritize novel effector transcripts for type 2 diabetes GWAS loci and may shed further light on the mechanisms underlying genes previously implicated in β-cell dysfunction.

The known MODY gene HNF4A was unexpectedly observed to cause a consistent increase (>40%) in insulin secretion across all conditions. HNF4A encodes the transcription factor HNF4α and is mutated in about 10% of all MODY cases (19). HNF4A loss-of-function mutations that cause monogenic diabetes later in life have also been associated with increased birth weight (indicative of increased fetal insulin secretion) and congenital hyperinsulinism in early infancy (20). The underlying reason for this switch from elevated to reduced insulin secretion is unknown, but it has been speculated that gradual β-cell exhaustion or, alternatively, a shift in the modulating cofactors of HNF4α may underlie this phenomenon (21,22).

Among the hits with limited prior evidence of a role in β-cell function, we independently validated ZMIZ1, ARL15, and THADA. Overexpression and knockdown of ZMIZ1, encoding ZMIZ1, has recently been shown to negatively impact on insulin secretion in primary human islets (6). Moreover, a nearby type 2 diabetes association signal overlaps a cis–expression quantitative trait loci for the gene, supporting its candidacy as the regional effector transcript (6). ARL15 encodes ARL15, a relatively uncharacterized member of the ARF family of proteins involved in regulation of vesicle trafficking and biogenesis. The gene is highly expressed in β-cells and located downstream of an islet-active enhancer bound by key β-cell transcription factors (18,23) (Supplementary Fig. 10A). THADA encodes the protein THADA and contains a coding disease-association signal that has also been associated with reduced β-cell function (24) (Supplementary Fig. 10C). Consistent with the directionality of our findings, expression profiling has shown the gene to be more highly expressed in patients with type 2 diabetes compared with control subjects (25). All three genes thus emerge as strong candidates for future studies.

While successfully enabling unbiased functional characterization, our current screening strategy has a number of limitations. False negatives (i.e., true causal genes not identified as hits) could arise as a result of primary effects of the causal gene on non–β-cell tissues or through effects on genes expressed at different developmental stages. Likewise, overexpression or greater knockdown efficiency may in some cases be required to expose a disease-relevant phenotype. Among the targets analyzed for silencing efficiency, a variable range of knockdown was observed (34%–88%), and some genes might remain undetected due to insufficient silencing. Conversely, false-positive effects (i.e., non–β-cell regulators identified as hits) also cannot be excluded, and unexpected findings should be further functionally validated (e.g., SLC2A4 effect on IBMX-stimulated insulin secretion). Though the EndoC-βH1 cell line has been found to recapitulate many aspects of β-cell function, it remains a possibility that some findings would not translate directly into human physiology. Finally, a subset of the identified hits may represent true β-cell regulators that are independent of any disease risk variants and, though still of biological importance, not genuine effector transcripts for type 2 diabetes. In addition to the possibility of more than a single effector transcript per locus, this phenomenon likely also contributes to the relatively high proportion of multihit loci observed in the primary screen (8/37).

Despite these limitations, our screening strategy successfully replicated well-established biological mechanisms and identified genes involved in β-cell function at half of the loci investigated. This demonstrates, for the first time, the feasibility of performing scalable screens for insulin secretory defects in human pancreatic β-cells and opens up the possibility for not only large-scale genetic manipulations but also compound high-throughput screening to therapeutically manipulate human β-cells. Insights from this and subsequent functional screens can be integrated with complementary lines of evidence from exome-wide association studies, chromatin conformation capture, and cis–expression quantitative trait loci studies to prioritize genes for follow-up studies. Ultimately, this could accelerate the translation of genetic association signals into molecular mechanisms for β-cell dysfunction, insulin insufficiency, and type 2 diabetes.

Funding. This study was funded in Oxford by the Wellcome Trust (095101/Z/10/Z and 098381). S.K.T. is a Radcliffe Department of Medicine Scholar. M.I.M. is a Wellcome Trust Senior Investigator. A.L.G. is a Wellcome Trust Senior Fellow in basic biomedical science.

Duality of Interest. M.v.d.B. is supported by a Novo Nordisk postdoctoral fellowship run in partnership with the University of Oxford. No other potential conflicts of interest relevant to this article were reported.

Author Contributions. S.K.T., A.C., D.E., M.I.M., and A.L.G. conceived and designed the study. S.K.T., M.v.d.B., M.I.M., and A.L.G. analyzed and interpreted the data. S.K.T., A.C., C.B., and A.B. performed the experiments. R.S. provided protocols. S.K.T., M.I.M., and A.L.G. wrote the manuscript. S.K.T., A.C., M.v.d.B., C.B., A.B., R.S., D.E., M.I.M., and A.L.G. edited and approved the manuscript. A.L.G. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Prior Presentation. Parts of this study were presented at the 52nd European Association for the Study of Diabetes Annual Meeting, Munich, Germany, 12–16 September 2016.