Glucose-Induced Changes in Gene Expression in Human Pancreatic Islets: Causes or Consequences of Chronic Hyperglycemia

The function of pancreatic islets is critical for maintaining glucose homeostasis, and dynamic changes of gene expression is part of the islets’ response to blood glucose changes. In patients with type 2 diabetes (T2D), islet function declines progressively. Although the initial pathogenic trigger of impaired β-cell function is still unknown, elevated glucose levels are known to further aggravate β-cell function, a condition referred to as glucotoxicity, which can stimulate apoptosis and lead to reduced β-cell mass (1–5). Prolonged exposure to hyperglycemia also can induce endoplasmic reticulum (ER) stress and production of reactive oxygen species (6), which can further impair islet function and thereby the ability of islets to secrete the insulin needed to meet the increased demands imposed by insulin resistance and obesity (7). Although these changes are likely to contribute to deterioration of islet function in patients with manifest disease, they are less likely to explain the development of T2D in individuals with normoglycemia.

In previous studies, we analyzed the gene expression profile in individuals with chronically elevated glucose as measured by elevated HbA_1c levels (8,9). However, these studies could not demonstrate whether the changes in gene expression are the cause or the consequence of hyperglycemia. One way to address this question is to compare gene expression in islets chronically exposed to hyperglycemia (prediabetes or diabetes) with gene expression changes after short-term exposure to hyperglycemia, with the assumption that gene expression changes seen in islets from patients with T2D but not after short-term hyperglycemia are the cause rather than the consequence of hyperglycemia (i.e., contributing to the pathogenesis of T2D). Thus, we performed RNA sequencing of human islets incubated at physiological (5.5 mmol/L) and high (18.9 mmol/L) glucose concentrations and compared the glucose-regulated genes with the gene expression profile seen in islets from hyperglycemic donors (8).

Human islets of Langerhans were obtained from the Human Tissue Laboratory (Lund University), which is funded by the Excellence Of Diabetes Research in Sweden (EXODIAB) network (www.exodiab.se/home) in collaboration with The Nordic Network for Clinical Islet Transplantation Program (www.nordicislets.org). All islet donors had given consent for donation of organs for medical research, and the procedures were approved by the ethics committee at Lund University (Malmö, Sweden; permit number 2011263). The islets were prepared from cadaver donors by using enzymatic digestion and density gradient separation. Islet preparation purity and count were determined as described previously (10). Between 3,000 and 5,000 islets each from 45 donors were divided into two pools that were incubated for 24 h in CMRL 1066 medium (ICN Biomedicals) containing 5.5 mmol/L or 18.9 mmol/L glucose at 37°C. The study design is described in Fig. 1A, and the clinical characteristics of the donors can be found in Supplementary Table 1. RNA from the islets was extracted by using the miRNeasy Mini Kit (QIAGEN). RNA quantity and integrity was assessed with an ND-1000 spectrophotometer (NanoDrop) and on a 2100 Bioanalyzer or a 2200 TapeStation instrument (Agilent Technologies).

With RNA sequencing, we performed a truly global analysis without a priori assumptions (11). A total of 1-μg high-quality RNA (RNA integrity number ≥8) was used as input for the TruSeq RNA Library Preparation Kit (Illumina). The resulting libraries were quality controlled on a 2200 TapeStation instrument before sequencing on a HiSeq 2000 system (Illumina) for an average depth of 32.4 mol/L paired-end reads (2 × 100 base pairs).

Insulin secretion capacity of the islet preparations was evaluated by stimulatory index (SI). Twenty handpicked islets were perifused with low glucose (1.67 mmol/L) for 42 min, high glucose (20 mmol/L) for 48 min, and then low glucose again. Fractions were collected at 6-min intervals, and the secreted insulin was measured by ELISA. SI was defined as the ratio between the areas under the curve calculated for the low and high glucose concentrations (12). SI was used for all estimates of insulin secretion in the in vitro experiments.

To test the effect of high glucose incubation for 24 h on islet viability and apoptosis, we incubated islets from three donors in normal and high glucose for 24 h as described above. More details can be found in the Supplementary Data.

The data were aligned to hg19 with STAR (Spliced Transcripts Alignment to Reference) (13) by using the GENCODE (Encyclopedia of Genes and Gene Variants) (14) v20 gene annotation to guide the alignment, and featureCounts (15) was used to count the number of reads aligned to the genes. Samples with <10 million counts and genes with <2 counts per million (CPM) in ≥10 samples were excluded from further analysis. The raw counts were transformed to log2 CPM, and the mean-variance trend was identified by using mean-variance modeling at the observational level (voom) (16) to allow for linear modeling after batch correction with ComBat (17). Differential expression was analyzed by using the paired data for 31 pairs of normoglycemic donor islets (HbA_1c ≤6% [42 mmol/mol]) and 14 pairs of hyperglycemic donor islets (HbA_1c >6% [42 mmol/mol]) separately by fitting a linear model with linear models for microarray and RNA sequencing data (limma) (18). To exclude that the number of genes responding to glucose was not due to differences in group size, we permuted 14 normoglycemic samples 1,000 times and calculated the mean number of changed genes by using the R function sample (19).

Linear modeling also was used to analyze differential gene expression (adjusted for age, sex, and days in culture) in untreated islets from 81 donors with normoglycemia (HbA_1c <6% [42 mmol/mol]), 19 donors with hyperglycemia (HbA_1c <6.5% [48 mmol/mol]), and 16 donors with T2D (HbA_1c ≥6.5% [48 mmol/mol]), including the samples in the glucose-treated data set and overlapping a previous study by us (8). The HbA_1c cutoffs for the groups were chosen to reflect the increased risk seen in individuals with HbA_1c 6.0–6.5% and to follow the recommendations by the International Expert Committee Report on the Role of the A1C Assay in the Diagnosis of Diabetes (20).

To test the effect of gene expression on insulin secretion in vivo in humans, we used exome genotype array data to find potential loss-of-function and gain-of-function variants in filtered genes. Genotypes with a frequency >5% from the human Infinium Exome-24 v1.1 array (Illumina) were obtained from the Prevalence, Prediction and Prevention-Botnia study (21), and association with corrected insulin response (CIR) was performed, adjusted for age, sex, and BMI, in 3,720 samples; unadjusted P <0.05 was considered nominally significant. The single nucleotide polymorphisms (SNPs) in coding parts of the gene were considered as potential loss-of-function/gain-of-function variants.

Islet expression quantitative trait loci (eQTLs) were used to investigate genetic findings from the RNA sequencing experiments. More details can be found in the Supplementary Data and Supplementary Table 2.

The aim of the project was to identify genes whose expression was altered in islets from donors with hyperglycemia but not changed by exposure to acute hyperglycemia or whose expression was altered in the opposite direction, assuming that such genes are more likely to be a cause than a consequence of hyperglycemia and thus involved in the pathogenesis of impaired insulin secretion leading to T2D. Altogether, expression of 717 genes (Fig. 1B and Supplementary Table 3) differed among donors with normoglycemia (HbA_1c <6%), hyperglycemia (HbA_1c ≥6% or <6.5%), and T2D (HbA_1c ≥6.5%) islets at a false discovery rate (FDR) of 5%. Of them, 392 were not affected by acute glucose exposure (see below and Supplementary Tables 4 and 5).

Islets from 45 human cadaver donors (31 normoglycemic and 14 hyperglycemic) were incubated at normal (5.5 mmol/L) or high glucose (18.9 mmol/L) concentrations for 24 h (Fig. 1A). RNA sequencing yielded an average 32.4 million paired-end reads per sample. A total of 15,958 genes were identified as expressed (defined as >2 CPM in >10 samples), corresponding to 27% of the genes in the GENCODE v20 gene annotation and 70% of the protein coding genes (n = 13,999). Incubation of islets from the 31 donors with normoglycemia (HbA_1c ≤6%) in high glucose resulted in up- or downregulation (FDR-corrected P value [q] < 0.05) of 4,658 genes (Supplementary Table 6 and Fig. 1B). Information on the top 10 genes whose expression increased or decreased after short-term hyperglycemia is provided in Fig. 2. In islets from donors with hyperglycemia (HbA_1c >6%), expression of 107 genes changed in response to exposure to glucose (Fig. 1B and Supplementary Table 7), which is far less than the 4,658 genes whose expression changed in normoglycemic islets. Part of this finding could be explained by differences in the number of donors (n = 31 vs. 14). However, after permuting 14 normoglycemic islets 1,000 times by using the sample package in R (19), the median number of differentially expressed genes was 501 (interquartile range 499) and significantly higher than that of donors with hyperglycemia (501 vs. 107; P = 0.04), demonstrating the unresponsiveness to glucose of islets from donors with hyperglycemia (Supplementary Fig. 1). Of these 107 genes, only 18 were affected by acute glucose in islets from donors with hyperglycemia (Supplementary Table 8).

To ensure that the gene expression changes we observed in the islets after glucose incubation is not due to differences in cell viability and/or cell numbers, we tested the effect of high glucose incubation on islet function by measuring apoptosis and cell viability. Islets from four normoglycemic islets were incubated for 24 h at a high (18.9 mmol/L) or normal (5.5 mmol/L) glucose concentration. Neither apoptosis (n = 3; P = 0.20) (Supplementary Fig. 2B) nor the number of viable cells (n = 3; P = 0.46) (Supplementary Fig. 2C) was significantly affected by high glucose incubation for 24 h. The glucose-stimulated insulin secretion is affected already after 24 h of preincubation in high glucose, with increased basal insulin secretion (n = 4; P = 0.01) (Supplementary Fig. 2A).

To identify genes possibly involved in T2D pathogenesis, we screened for genes whose expression was unaffected by acute hyperglycemia, correlated negatively with insulin secretion, and showed increased expression in donors with hyperglycemia compared with donors with normoglycemia (or vice versa if protective). Expression of 392 genes differed between the donor groups (FDR <0.05) (Supplementary Tables 4 and 5) and was not affected by acute glucose exposure in the same direction. Of them, 92 genes correlated with insulin secretion (Table 1 and Fig. 3).

Noncoding genetic variants influencing gene expression are referred to as eQTLs and can be used to explore the effects of gene expression on a phenotype. If the variant was associated with the expression of a gene whose expression correlates with insulin secretion, we used that variant as a proxy for insulin secretion. For the eQTLs analysis, we used genome-wide association study (GWAS) genotype data in combination with RNA sequencing data from 191 donors. We identified one or more eQTLs in 16 of the 392 genes that differed between donors with normoglycemia and donors with hyperglycemia and was not affected by acute glucose exposure in the same direction, including the SID1 transmembrane family member 1 (SIDT1) and forkhead box E1 (FOXE1) genes; these variants also were nominally associated with SI in the islets. The same eQTLs in the SIDT1 gene were nominally associated (P < 0.05) with T2D in the DIAGRAM (Diabetes Genetics Replication and Meta-analysis) study (22) (Supplementary Table 9). An eQTL SNP in the rhomboid 5 homolog 1 gene (RHBDF1) was nominally associated with CIR.

To increase the likelihood that a gene would be involved in the pathogenesis of or protection from T2D, we searched for potential loss- or gain-of-function variants in the 92 genes whose expression correlated with SI and explored whether such variants would influence in vivo insulin secretion (CIR). Five genes (transmembrane protein 132C [TMEM132C], ERO1LB, dedicator of cytokinesis protein 10 [DOCK10], proline-rich 14-like protein [PRR14L], and IGSF11) (Fig. 3D and Supplementary Table 10) harbored variants, which were associated with CIR corrected for age, sex, and BMI. Expression of ERO1LB, DOCK10, PRR14L, and IGSF11 was positively correlated with insulin secretion and downregulated in donors with hyperglycemia and, thus, were potentially protective against T2D. For example, the expression of ERO1LB was lower in donors with hyperglycemia (β = −0.37, q = 0.02) and correlated positively with insulin secretion (Pearson r = 0.29; P = 0.002). A coding variant (with a possible effect on expression) in this gene was nominally associated with decreased CIR (rs2477599, allele T, β = −0.03; P = 0.03), suggesting that expression of ERO1LB is required for normal insulin secretion. In contrast, expression of TMEM132C was higher in donors with hyperglycemia (β = 0.40, q = 0.02) and correlated negatively with insulin secretion (Pearson r = −0.23; P = 0.01). A coding variant in TMEM132C was associated with impaired CIR (rs11059681, allele G, β = −0.03; P = 0.03), making it a potential candidate to contribute to the pathogenesis of T2D.

To obtain some insight into the potential role of genes/loci previously reported to be associated with T2D or glycemic traits, we analyzed whether expression of 134 genes associated with T2D or glycemic traits was influenced by hyperglycemia and whether they harbored variants influencing expression (eQTLs) (22–29). Expression of 21 genes harboring such variants (Supplementary Table 11) changed in response to short-term hyperglycemia, whereas expression of 7 genes differed between individuals with normal glucose tolerance (NGT) and abnormal glucose tolerance (AGT) (i.e., hyperglycemia [SLC2A2 (β = −0.66, q < 0.01), RASGRP1 (β = −0.38, q < 0.01), PDX1 (β = −0.33, q < 0.05), PCSK1 (β = −0.32, q < 0.05), ZBED3 (β = −0.20, q < 0.05), PROX1 (β = −0.18, q < 0.05), FTO (β = −0.12, q < 0.05)]). Expression of PROX1 and ZBED3 also correlated positively with SI (r = 0.19 [P < 0.05] and 0.20 [P < 0.05], respectively). The BCAR1 T2D-associated SNP rs7202877 was also an eQTL in the islets. The SNP was associated with SI (effect allele G, β = 1.87; P < 0.05) but was not an eQTL for the BCAR1 but for the pseudogene RP11-331F4.4. Another 13 T2D GWAS SNPs influenced expression of nearby genes (cis-eQTLs) whereof 3 influenced more than one gene (rs3132524, rs1167800, and rs3829109). Seven of these 13 eQTLs were not eQTLs for the gene suggested in the GWAS, such as rs1046896, which tags fructosamine 3 kinase–related protein (FN3KRP) (Table 2). The expression of FN3KRP was negatively correlated with SI (r = −0.19; P < 0.05). None of these eQTL genes were differentially expressed in islets from donors with hyperglycemia compared with islets from donors with normoglycemia, although expression of STARD10 was upregulated by acute glucose exposure (log2 fold change [log2FC] = 0.20, q < 0.05). The SNP rs10830963, which is one of the strongest eQTLs for the melatonin receptor 1B gene (MTNR1B) in human islets (30) was nominally associated with CIR (G allele, β = −0.11; P < 0.05). The same was seen for rs11257655 in the calcium/calmodulin-dependent protein kinase ID (CAMK1D) gene (T allele, β = −0.14; P < 0.01).

We assumed that similar changes in gene expression seen after acute and long-term (prediabetes and diabetes) glucose exposure could possibly be attributed to glucotoxicity. A total of 325 genes (46% of the genes changed in donors with hyperglycemia and T2D) had similar gene expression changes both after 24-h glucose exposure and in islets from donors with hyperglycemia (Supplementary Table 12). Of them, the expression of 73 genes correlated with SI (48 positively and 25 negatively), which is more than expected by chance (22% correlate with SI vs. 9% of genes in general; Fisher exact test P < 0.0001) (Fig. 4 and Table 3). The genes A-kinase anchoring protein 6 (AKAP6), GALNT2, and FERMT1 also harbored coding variants associated with CIR (Supplementary Table 13). Finally, expression of the SIPA1L2, HRK, TMED132D, MBP, and CPEB1 genes was also affected by acute glucose exposure in islets from donors with hyperglycemia; gene expression of TMEM132D and MBP correlated with SI (Supplementary Table 14). Of note, MBP had an eQTL SNP (rs1667903, effect allele C, β = −0.50, q < 0.05) that was associated with both SI and CIR (β_SI = −1.52 [P < 0.05] and β_CIR = 0.14 [P < 0.05]).

Differences in gene expression between human pancreatic islets from donors with hyperglycemia and normoglycemia may be a physiological consequence or the cause of elevated glucose. These effects have been difficult to separate in cross-sectional studies, which have not explored changes in gene expression in response to short-term (acute) hyperglycemia. We provide such a map in islets from both donors with normoglycemia and donors with hyperglycemia exposed to high glucose for 24 h. Not surprisingly, because glucose is a transcriptional modulator (31), ∼30% of the genes expressed in the human islets are acutely regulated by glucose (Supplementary Table 6). Only 20% of genes changed in response to short-term hyperglycemia in islets from donors with hyperglycemia (Supplementary Table 7), suggesting that they already are turned on or off by chronic hyperglycemia. Of note, these islets were subjected to a washout period where they were cultured at 5.5 mmol/L glucose for several days after isolation, so the lack of acute effect of glucose might also be the consequence of a metabolic memory.

To identify genes that might be causally involved in the pathogenesis of T2D, we postulated that they would not only show differences in expression between islets from donors with hyperglycemia or normoglycemia but also show that expression would not change in the same direction after short-term hyperglycemia. Furthermore, we expected that these changes in gene expression would correlate with insulin secretion. One of these genes was ERO1LB that encodes an oxidoreductase involved in the folding of proinsulin in the ER. Knockdown of ERO1LB in the insulin misfolding–prone Akita mouse resulted in islet destruction and development of diabetes (32). In support of this, expression of Ero1b was reduced in islets from two diabetic mouse models. However, overexpression of Ero1b also seems to induce ER stress in β-cells, leading to impaired glucose-stimulated insulin release and suggesting a possible U-shaped effect of ERO1LB expression on ER stress (33). Four other putatively causative genes (TMEM132C, DOCK10, PRR14L, and IGSF11) harbored potential loss-of-function variants associated with in vivo insulin secretion. TMEM132C is an interesting candidate, but little is known about its role in islet function. Expression of DOCK10 and PRR14L were both possibly protective in the islets. Dock10 is a potential guanine nucleotide exchange factor required for the activation of the rho family of GTPases, and Dock10 interacts with Cdc42 to promote the formation of the guanosine triphosphate–bound activated form (34). Cdc42 activated by glucose (35) has previously been associated with insulin secretion in rodent and human islets (36). Small interfering RNA–mediated silencing of Cdc42 in isolated islets results in impaired second phase insulin secretion (37). DOCK10 deficiency could, therefore, potentially impair insulin secretion through its effect on GTPases, such as Cdc42.

In addition, two genes, SIDT1 and RHBDF1, that showed different expression in islets from donors with hyperglycemia and normoglycemia harbored eQTLs that were associated with insulin secretion. The A allele of eQTL SNP rs11929640 was associated with increased expression of SIDT1 in the islets and with increased in vitro insulin secretion. Expression of SIDT1 was decreased in islets from the donors with hyperglycemia and correlated positively with SI, suggesting a link between decreased expression and impaired insulin secretion. SIDT1 encodes an RNA transporter that can transfer small RNA molecules, such as micro RNA, inside the cell (38). The RHBDF1 gene showed increased expression in islets from donors with hyperglycemia and correlated negatively with insulin secretion. In contrast, the T allele of the eQTL SNP rs9930775 was associated with decreased expression and in vivo insulin secretion, suggesting that the first finding was a consequence of hyperglycemia, whereas the eQTL could reflect a causal genetic mechanism. These findings will require validation by genetic engineering, predominantly in human islets.

Thirteen of the 134 established T2D/glycemic trait loci had eQTLs in human islets, with most of them showing cis effects (only 7 trans). Expression FN3KRP correlated negatively with insulin secretion in islets. Of note, the T allele of the eQTL SNP rs1046896 in FN3KRP has been previously associated with increased HbA_1c levels (23). Because expression of FN3KRP was decreased in T allele carriers, the data suggest that normal expression of FN3KRP is required for maintaining normal insulin secretion. Of note, the GWAS SNP is located in the 3′ untranslated region of FN3KRP, but the gene suggested in the literature to be affected by the SNP is the paralog FN3K. Both genes encode enzymes involved in deglycation of proteins, but they have different substrates. FN3K, but not FN3KRP, can phosphorylate fructosamine, which is glycated by glucose (39).

In line with other published results from our group (40), the rs10830963 T2D risk allele was a strong eQTL for increased expression of MTNR1B and was associated with decreased in vivo insulin secretion. Paradoxically, rare loss-of-function variants in MTNR1B have also been associated with elevated glucose levels (41). Although the common rs10830963 is the strongest eQTL observed in human islets, rare loss-of-function variants would most likely exert their effects in most organs where MTNR1B is expressed and, thereby, possibly increase the risk of diabetes by other mechanisms than the common variant.

The T2D risk locus rs11257655 located between CDC123 and CAMK1D was an eQTL for the CAMK1D gene. The risk allele was associated with increased expression of CAMK1D and associated with in vivo insulin secretion. This risk allele has been shown to increase transcriptional activity possibly through FOXA1 and FOXA2 (42). The current results agree that as with eQTLs observed in other tissues, the locus regulates the expression of CAMK1D, not CDC123. Although the role of CAMK1D in the pathogenesis of T2D remains speculative, it may involve regulation of gene transcription through phosphorylation and activation of CREBP (43). CREBP is important for many aspects of β-cell function, including insulin exocytosis (44).

This study also provides information on genes sensitive to glucotoxicity, assuming that genes changed similarly by both acute and chronic hyperglycemia are susceptible to glucotoxicity (Fig. 4 and Table 3). Three of these genes harbored coding variants associated with in vivo insulin secretion: AKAP6, GALNT2, and FERMT1. AKAP6 could possibly influence insulin secretion through modulation of protein kinase A because mice with disruption of another A-kinase anchoring protein (AKAP150) show impaired insulin secretion (45). Expression of AKAP6 was positively correlated with insulin secretion, suggesting that glucose-induced downregulation of AKAP6 could result in impaired insulin secretion. In support of this, a coding variant in the gene was negatively associated with insulin secretion in vivo. Moreover, expression of TMEM132D and MBP was decreased after both acute and chronic hyperglycemia and correlated positively with insulin secretion, supporting the view that their downregulation could be involved in impaired insulin secretion. This is further supported by previous work showing that MBP stimulates insulin secretion (46). Both genes also were affected in islets from donors with hyperglycemia exposed to short-term hyperglycemia (Supplementary Table 14), suggesting that they are sensitive to glucotoxicity.

Studies of human islets are subject to a number of difficulties. A considerable number of factors will affect the islets’ phenotype connected to the nature of death of the donor, harvest of the organs, preparations of the islets, and so forth. In this study, we attempted to address these issues by using paired data where the changes we see by treating the islets with glucose are observed within the islet preparation and the effect is steady across the data set. This approach will decrease the contribution of potential effects of external and individual factors in each islet preparation. We also used a considerable number of samples in the study to present stable changes and differences. Although most of the results at this stage are descriptive, they provide a catalog of changes in gene expression in human pancreatic islets after acute and chronic exposure to glucose that can serve as a resource for the dissection of the molecular mechanisms leading to or protecting from T2D.

Acknowledgments. Human pancreatic islets were obtained from The Nordic Network for Clinical Islet Transplantation, supported by the Swedish national strategic research initiative EXODIAB and JDRF.

Funding. Support for this research was provided by the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement no. 269045 awarded to L.G. The work was also supported by project grants from the Vetenskapsrådet (Swedish Research Council) to L.G. (Dnr 2010-3490) and Pfizer.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. E.O.-L. and U.K. performed experiments. E.O.-L., U.K., P.S., R.B.P., N.O., E.A., J.F., and P.V. performed analysis of the data. E.O.-L., U.K., P.S., R.B.P., N.O., E.A., J.F., O.H., L.G., and P.V. critically revised and contributed to the manuscript. E.O.-L., L.G., and P.V. wrote the draft of the manuscript. O.H., L.G., and P.V. designed the study. P.V. 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.