Whole-Genome Bisulfite Sequencing of Human Pancreatic Islets Reveals Novel Differentially Methylated Regions in Type 2 Diabetes Pathogenesis
Abstract
Current knowledge about the role of epigenetics in type 2 diabetes (T2D) remains limited. Only a few studies have investigated DNA methylation of selected candidate genes or a very small fraction of genomic CpG sites in human pancreatic islets, the tissue of primary pathogenic importance for diabetes. Our aim was to characterize the whole-genome DNA methylation landscape in human pancreatic islets, to identify differentially methylated regions (DMRs) in diabetic islets, and to investigate the function of DMRs in islet biology. Here, we performed whole-genome bisulfite sequencing, which is a comprehensive and unbiased method to study DNA methylation throughout the genome at a single nucleotide resolution, in pancreatic islets from donors with T2D and control subjects without diabetes. We identified 25,820 DMRs in islets from individuals with T2D. These DMRs cover loci with known islet function, e.g., PDX1, TCF7L2, and ADCY5. Importantly, binding sites previously identified by ChIP-seq for islet-specific transcription factors, enhancer regions, and different histone marks were enriched in the T2D-associated DMRs. We also identified 457 genes, including NR4A3, PARK2, PID1, SLC2A2, and SOCS2, that had both DMRs and significant expression changes in T2D islets. To mimic the situation in T2D islets, candidate genes were overexpressed or silenced in cultured β-cells. This resulted in impaired insulin secretion, thereby connecting differential methylation to islet dysfunction. We further explored the islet methylome and found a strong link between methylation levels and histone marks. Additionally, DNA methylation in different genomic regions and of different transcript types (i.e., protein coding, noncoding, and pseudogenes) was associated with islet expression levels. Our study provides a comprehensive picture of the islet DNA methylome in individuals with and without diabetes and highlights the importance of epigenetic dysregulation in pancreatic islets and T2D pathogenesis.
Introduction
Impaired insulin secretion is a key feature of type 2 diabetes (T2D). However, the molecular mechanisms underlying pancreatic islet dysfunction in T2D are largely unknown. Although genetic risk factors are known to contribute to T2D, <20% of the estimated T2D heritability can be explained by single nucleotide polymorphisms (SNPs) identified by genome-wide association studies (1). Hence, one must look elsewhere to find disease-causing mechanisms. Given the important role of environmental factors in T2D pathogenesis, mechanisms mediating gene-environment interactions, such as epigenetics, may be of particular significance. Studying epigenetic processes in the tissue of primary pathogenetic importance, the pancreatic islets, may reveal central mechanisms for T2D. Indeed, we and others have identified altered DNA methylation in islets from human T2D donors (2–4). However, the methods used in previous studies only covered up to ∼1.5% of genomic CpG sites (5). To obtain a more complete picture of the human islet methylome and to further dissect the impact of epigenetics in T2D, analyses covering the majority of CpG sites are needed.
Whole-genome bisulfite sequencing (WGBS) is the most comprehensive method to study DNA methylation on a single nucleotide resolution. This method is costly and most studies have so far been based on few samples in selected sets of human tissues (6–8). The largest effort to describe the human epigenome was performed by the National Institutes of Health Roadmap Epigenomics Consortium (9). However, previous WGBS studies did not include human islets or T2D case-control cohorts. To address this knowledge gap and to dissect epigenetic alterations in T2D, we performed WGBS of human islets from control subjects without diabetes and T2D donors.
Research Design and Methods
Human Islets
Human islets were obtained from the Nordic Network for Islet Transplantation. The islet cohort is described in the Supplementary Data.
WGBS
A total of 300 ng of islet DNA was bisulfite treated with the EZ DNA Methylation-Gold D5005 kit (Zymo Research Corp., Irvine, CA) according to the manufacturer’s protocol (Supplementary Data).
WGBS Data Analysis
WGBS data in FASTQ format generated using the Illumina HiSeq platform (Illumina, San Diego, CA) were used for further analyses (Supplementary Data).
Infinium 450K Array
DNA methylation of human islets (Table 1) was analyzed with the Infinium HumanMethylation450K BeadChip (5) (Supplementary Data).
Characteristics for donors of pancreatic islets included in the WGBS analysis
RNA-seq Data
The RNA library was prepared with the TruSeq kit (Illumina) and then sequenced with HiSeq2000 as previously described (10) (Supplementary Data).
Genomic Annotation
Genomic elements, such as transcription start sites (TSS), transcription end sites (TES), and exons and introns for 198,442 transcripts corresponding to 60,483 genes, were extracted from GENCODE version 22 (GRCh38). Each differentially methylated region (DMR) was annotated based on its position in relation to all transcripts above; hence, one DMR can have multiple annotations.
Pyrosequencing
DNA methylation for biological replication of a genomic region in the most significant PDX1 DMR (chr13:27921804:27925104) and regions in DMRs annotated to ARX, CADM1, PARK2, PID1, SLC2A2, and SOCS2 was analyzed by pyrosequencing and the PyroMark Q96 ID (Qiagen, Hilden, Germany). Primers were designed using PyroMark Assay Design 2.0 (Supplementary Table 1). Procedures were performed according to recommended protocols.
Luciferase Assays
Luciferase assays were performed as previously described (4) (Supplementary Data).
β-Cell Lines
INS-1 832/13 rat β-cells were used in functional experiments as previously described (2) (Supplementary Data).
Statistical Analysis
Donor characteristics, pyrosequencing data, and luciferase experiments were analyzed using unpaired t tests (Supplementary Data).
Results
WGBS in Human Islets
To characterize the methylome in human islets, we generated WGBS data from eight control subjects and six T2D donors (Table 1). An average of 74% of the resulting reads per sample was uniquely mapped to the human reference genome (hg38). Sequencing information and alignment statistics are reported in Supplementary Table 2. The samples were sequenced with an average coverage of 21× per base and methylation levels of ∼2.4 × 107 CpG sites (∼83% of all genomic CpG sites) on the forward strand were obtained for all samples.
We compared the methylation data obtained by WGBS with data generated with the Infinium 450K array for the same 14 samples. The replicates of each islet sample analyzed by both WGBS and microarray showed high reproducibility (Supplementary Table 2).
One islet sample was analyzed by WGBS using both Illumina type 3 and 4 chemistry and by two library preparation kits; the TruSeq (EpiGnome) DNA Methylation Kit and the NEXTflex Bisulfite Library Prep Kit. The high correlation between the WGBS data generated with the two different sequencing chemistries (r = 0.995) and with different library preparations (r = 0.987) further confirmed the data quality (Supplementary Fig. 1A).
Human Islet Methylome
We first sought to characterize the overall variability of the human islet methylome. The degree of DNA methylation throughout the genome was highly correlated among the analyzed islets samples (r = 0.973–0.989) and the average methylation level was 75.9%. The distribution of the methylation level in human islets is bimodal, with the highest peak at 90.2%, showing that most CpG sites are highly methylated, and the second peak at 1.4%, representing CpG sites with very low levels of methylation (Fig. 1A).
The DNA methylome of human pancreatic islets. A: The distribution of DNA methylation in human islets. B–D: Average DNA methylation in human pancreatic islets, separated by different Gencode transcript types (protein-coding, noncoding [long and small combined], and pseudogenes) and gene regions. Here, TSS 50 kb represents 1,501–50,000 bp upstream from the TSS, TSS 1500 represents 201–1,500 bp upstream from the TSS, TSS 200 represents 1–200 bp upstream from the TSS, and TES 10 kb represents 1–10,000 bp downstream of the TES. Donors with T2D and normoglycemic control subjects display no differences in genome-wide average methylation for any genomic region or transcript type (P = 0.4–1.0). E and F: Average DNA methylation for regions overlapping with different histone marks or transcription factor binding sites. G–I: Average methylation levels are significantly different between transcripts of different expression levels. Here, we identified the not expressed transcripts based on transcript per million <0.1 and then divided the expressed transcripts into three equally sized groups that we categorize into low-, medium-, and high-expressed transcripts. Additionally, different genomic regions display specific methylation patterns that are also dependent on transcript type. Data presented as mean ± SEM (***P < 0.0001, as analyzed by ANOVA).
Next, we computed the average genome-wide methylation level in relation to different genomic regions in human islets. While introns and exons had the overall highest degree of methylation (78.5% and 77.4%), regions close to the TSS, such as the 1st exon and promoter regions (TSS 200 and TSS 1500), had the lowest degree of methylation (34.7%, 25.4%, and 44.4%, respectively). Regions more distant to the TSS had a rather high methylation level (TSS 50 kb 73.0% and TES 10 kb 71.2%).
We also studied methylation based on annotation to different types of transcripts; protein-coding genes, noncoding RNAs, and pseudogenes (Fig. 1B–D). There was a striking difference in average DNA methylation between the different transcript types, where the typical drop in methylation seen close to the TSS for protein-coding genes (TSS 1500, TSS 200, and 1st exon) was less pronounced in noncoding RNAs and almost completely absent in pseudogenes. Additionally, there were no differences in average methylation for either genomic regions or transcript types between islets from T2D and control donors (P = 0.4–1.0).
To gain further insight into how the human islet epigenome is coordinated, we studied the relationship between DNA methylation and histone modifications across the genome. Here, we integrated our WGBS data with histone modifications generated in human islets by the Roadmap Epigenomics Consortium (9). Regions occupied by histone modifications associated with active chromatin had the lowest degree of methylation (11.8% for H3K9ac and 7.8% for H3K4me3), whereas regions occupied by modifications associated with repressive chromatin had higher methylation levels (49.6% for H3K27me3 and 80.2% for H3K9me3) (Fig. 1E). Additionally, regions occupied by histone modifications considered to be enriched at enhancer regions had either quite low (21.4% for H3K27ac) or high (69.8% for H3K4me1) degree of methylation (Fig. 1E). We also examined the relationship between DNA methylation levels and genomic binding of islet-specific transcription factors (PDX1, FOXA2, MAFB, NKX6.1, and NKX2.2) by combining our WGBS data with published ChIP-seq data (11). Regions occupied by transcription factors had a lower degree of methylation (range 20.9%–43.9%) compared with the whole genome (Fig. 1F). There were no differences in the methylation levels in any of the regions occupied by the studied histone marks or transcription factors in control versus T2D islets (Fig. 1E and F) (P = 0.4–1.0).
We proceeded to explore the relationship between DNA methylation and gene expression levels using WGBS and RNA-seq data from the same 14 human islet donors. Here, we categorized 60,483 transcripts into not expressed (TPM <0.1, 38,261 transcripts) and equally sized groups of low- (7,407), medium- (7,407), and high-expressed (7,408) transcripts. We also studied protein-coding genes, noncoding RNAs, and pseudogenes separately. We found an association between DNA methylation and the expression level in all genomic regions and for all types of transcripts (P < 0.0001) (Fig. 1G–I). Importantly, protein-coding genes not expressed were hypermethylated in regions close to the TSS (TSS 1500, TSS 200, and 1st exon) (Fig. 1G). However, protein-coding nontranscribed genes had significantly lower methylation levels in exon and intron regions compared with the transcribed genes (Fig. 1G), supporting the hypothesis that increased methylation in the gene body is associated with a higher gene transcription (12). More modest differences in methylation were seen when comparing the different expression levels of transcribed genes, i.e., dividing transcribed genes in low-, medium-, and high-expressed genes (Fig. 1G–I). Some similarities in the methylation pattern were seen for both noncoding RNAs and pseudogenes compared with protein-coding genes (Fig. 1H and I). However, regions close to the TSS were hypermethylated in nontranscribed genes and the methylation level in these regions was much higher in transcribed noncoding RNAs and pseudogenes (Fig. 1H and I).
Principal Component Analysis of the Islet WGBS Data
In an unsupervised principal component analysis of the islet WGBS data, the methylation data was segregated according to sex (Supplementary Fig. 1B), which is in agreement with our published 450K array data (13). Next, we correlated the top five principal components of the WGBS data with T2D, age, sex, and BMI. Here, T2D and sex correlated with one of the top five principal components (P < 0.002) (Table 2).
P values for correlations of the top five principal components for the WGBS data in human pancreatic islets with T2D, age, sex, and BMI
DMRs in Human T2D Islets
To address the epigenetic basis of T2D, we used BSmooth (14) to analyze the islet WGBS data. BSmooth determines methylation in WGBS data and identifies DMRs that account for biological variability. DMRs were defined as regions of three or more consecutive differentially methylated sites with an average absolute methylation difference ≥5% between groups. This analysis identified 25,820 DMRs in T2D islets (Supplementary Table 3). These were included in an unsupervised hierarchical clustering analysis presented as a heatmap in Fig. 2A, which shows distinction in methylation between donors with and without diabetes. A total of 13,696 DMRs showed average increased and 12,124 decreased levels of methylation in T2D islets. The mean DMR size was 414 bp (range 6–3,411 bp), and the mean CpG site count in the DMRs was 8.7 (range 3–164). The maximum methylation difference for a DMR was 27.5% when comparing diabetic versus control islets, and 692 DMRs had a methylation difference ≥10% (Supplementary Table 3). Among the DMRs found to have the largest absolute differences in methylation were regions annotated to ARX and TFAM (Fig. 2B and C), both important in islet function (15,16).
DMRs in human pancreatic islets from donors with T2D. A: Heatmap of the 25,820 T2D-associated DMRs. Among the DMRs with the largest absolute difference were regions upstream of ARX (B) and TFAM (C). Among the most significant DMRs were two large intergenic regions of PDX1, chr13:27921804-27925104 (D) and chr13:27926170-27928845 (E). F: Biological replication of the PDX1 region located in intronic regions and exon two and presented in panel D in human pancreatic islets from an independent cohort of 56 normoglycemic control and 19 T2D donors. G and H: Validation of DNA methylation in the PDX1 distal promoter in human pancreatic islets from T2D and control donors. G: Data from Yang et al. (4), produced using Sequenom's EpiTYPER technology. H: Data from current study, based on WGBS DMR chr13:27918705-27919232. Data presented as mean ± SEM (*P < 0.05).
Interestingly, two of the most significant DMRs covered 164 and 105 CpG sites that span 3,301 and 2,676 bp regions of PDX1 (Fig. 2D and E), a key islet transcription factor (17). In total, seven DMRs were annotated to PDX1 (Supplementary Table 3). Using pyrosequencing, we biologically replicated differential methylation of sites in the most significant DMR, located in PDX1 (chr13:27921804-27925104) (Fig. 2D and F), in an independent cohort (Supplementary Table 4). We previously reported increased methylation of 10 sites and decreased expression of PDX1 in human T2D islets (4). Importantly, a DMR from the current study (chr13:27918705-27919232) confirmed our previous finding, and comparisons of seven PDX1 sites covered by both studies validated the significant association with T2D (Fig. 2G).
We next examined the overlap between T2D-associated islet DMRs and 65 T2D candidate genes identified by the DIAGRAM (DIAbetes Genetics Replication And Meta-analysis) consortium (18). Among our islet DMRs, 159 were annotated to 43 known T2D candidate genes (Supplementary Table 5). Out of these, the DMR with the largest difference in methylation between donors with diabetes and donors without diabetes was located in ADCY5. T2D candidate genes with the highest number of DMRs were GLIS3, THADA, KCNQ1, and TCF7L2 (16, 14, 10, and 9, respectively).
We further investigated if our islet DMRs covered any of the SNPs reported to be associated with T2D (www.genome.gov/gwastudies, accessed 2 February 2016). Among these SNPs, rs163184 was covered by a DMR annotated to KCNQ1 (chr11:2825382:2826548) and rs11257655 by a DMR annotated to RP11 (chr10:12265391:12266540).
Distribution of T2D-Associated Islet DMRs
We then examined the genomic distribution of T2D-associated islet DMRs. We found ∼55% of the DMRs located within TSS 50 kb (1,501–50,000 bp upstream from the TSS), ∼1.5% located within TSS 1500 (201–1,500 bp upstream from the TSS), and ∼1.0% located within TSS 200 (1–200 bp upstream from the TSS) (Fig. 3A). Furthermore, ∼1.4% of all DMRs were located in the 1st exon, ∼3.0% in subsequent exons, and ∼20% in introns, whereas ∼12.5% of the DMRs were located within TES 10 kb (1–10,000 bp downstream of the TES). Finally, ∼5%–6% of all DMRs were located >50 kb from the nearest transcript and considered intergenic.
Genomic distribution of T2D-associated DMRs. A: Genomic distribution of T2D-associated islet DMRs, separated based on increased (n = 13,696) or decreased (n = 12,124) average DNA methylation. Each DMR can be annotated to several transcripts and will then be counted in all gene regions, except from intergenic DMRs that only include DMRs located more than 50 kb from a transcript and thereby are not annotated to any of the TSS, intragenic, or TES regions. B: Overlap between our T2D-associated DMRs and chromatin state, transcription factor binding, and active enhancer regions. Arrows represent significant over- or underrepresentation (P < 1 × 10−6). C: Enrichment of transcription factor recognition sequences in the T2D-associated islet DMRs based on HOMER (19). CTCF, CCCTC binding factor.
To explore the relationship between chromatin state and DNA methylation, we integrated our DMRs with published maps of histone modifications and enhancer regions in human islets (9,11). When comparing our T2D-associated islet DMRs with epigenetic marks generated in human islets by the Roadmap Epigenomics Consortium (9), 12.4% of our islet DMRs were occupied by modifications associated with active chromatin (3,203 by H3K4me3 and 3,194 by H3K9ac) (Supplementary Table 6 and Fig. 3B), which are significant enrichments for both (P < 1 × 10−6). Here, DMRs with higher methylation in donors with diabetes occupied 1,872 H3K4me3 and 1,851 H3K9ac marks. Moreover, 14.5% and 19.1% of islet DMRs were occupied by histone modifications enriched at enhancer regions (3,744 by H3K27ac and 4,935 by H3K4me1, respectively) (P < 1 × 10−6). Here, DMRs with higher methylation in donors with diabetes occupied 2,358 H3K27ac and 2,478 H3K4me1 marks. A smaller fraction of the DMRs were occupied by modifications associated with repressed chromatin (1,364 DMRs [5.3%] by H3K27me3 and 120 DMRs [0.5%] by H3K9me3). H3K27me3 was overrepresented, whereas H3K9me3 was underrepresented (P < 1 × 10−6) in this overlap. Moreover, DMRs with higher methylation in donors with diabetes occupied 583 H3K27me3 and 51 H3K9me3 marks. There was also an underrepresentation of DMRs (581, or 2.3%) overlapping with H3K36me3 (P < 1 × 10−6), and out of those, 381 showed higher methylation in donors with diabetes. Additionally, 618 islet DMRs overlapped with active enhancer regions identified in islets (Supplementary Table 6 and Fig. 3B) (11), which is more than expected by chance (P < 1 × 10−6), and out of those, 381 showed higher methylation in donors with diabetes.
Next, Hypergeometric Optimization of Motif EnRichment (HOMER) analysis (19) showed that motifs specific to key transcription factors in islets, including FOXA2, NeuroD1, MAFA, RFX, PDX1, and HNF1, and binding sites for the insulator CCCTC binding factor were significantly enriched in T2D-associated islet DMRs (Fig. 3C and Supplementary Table 7). We also identified novel motifs with statistical enrichment (Supplementary Table 7).
To gain insight into the relationship between the genomic binding of islet-specific transcription factors and differential methylation in islets from T2D donors, we used ChIP-sequencing data for five key transcription factors (FOXA2, MAFB, NKX2.2, NKX6.1, and PDX1) in human islets (11). We found an overrepresentation of DMRs overlapping with the binding sites of each of these transcription factors (P < 1 × 10−6) (Supplementary Table 8 and Fig. 3B). Interestingly, we found an overlap between T2D-associated DMRs annotated to SLC2A2, KCNJ11, and PDX1 and sites bound by PDX1 (Supplementary Table 8). These data suggest that differential methylation in transcription factor binding sites may be of importance in T2D.
Tissue-Specific DMRs
To identify DMRs in human islets that are also altered between other tissues or cell types, we analyzed the overlap between our T2D-associated islet DMRs and a set of 716,087 cross-tissue dynamic DMRs (20). Of our 25,820 DMRs, 12,911 (49.8%) overlapped with a dynamic DMR identified in other tissues (P < 1 × 10−6) (Supplementary Table 9).
Altered Expression of Genes Annotated to T2D-Associated DMRs
To determine whether genes annotated to T2D-associated DMRs also show altered expression in T2D versus control islets, we combined the 25,820 DMRs presented in Supplementary Table 3 with islet RNA-seq data of a previous publication (10). We identified 457 genes that had both significantly altered expression (false discovery rate <5%; q <0.05) and DMRs in T2D islets (Supplementary Table 10). These include genes important for islet function and metabolism, such as CACNA1D, CHL1, GLP1R, IGF1R, IL6, NR4A3, PARK2, PDX1, PID1, SEPT9, SIK2, SLC2A2 (also known as GLUT2), SOCS2, and SOX6 (Fig. 4A) (2,13,21–32). Moreover, 26 genes exhibited differential expression and DMRs with a methylation difference ≥10%.
Altered expression, T2D-associated DMRs, and functional consequences in β-cells. In total, 457 genes had both significantly altered expression in pancreatic islets from donors with T2D (q <0.05) (10) and significant DMRs (Supplementary Table 10). These include genes of importance in islet function and metabolism as depicted in panel A. B: Selected KEGG pathways based on enrichment of genes annotated to T2D-associated DMRs that also show altered RNA expression in human pancreatic islets from donors with T2D (10). A total of 457 genes were included in the pathway analysis and the boxes include genes contributing to the enrichment score for each pathway. A full list of significantly enriched KEGG pathways is found in Supplementary Table 11. C: DMRs from the four genes selected for functional studies; NR4A3, PID1, SOCS2, and PARK2. Overexpression of Nr4a3, Pid1, and Socs2 in INS-1 832/13 cells was verified by qPCR (D–F) and Western blot (G). ***P < 0.001 and **P < 0.01, as analyzed by one-tailed paired t tests (n = 7). H: Overexpression of Nr4a3, Pid1, and Socs2 resulted in perturbed insulin secretion (*P < 0.05 compared with pcDNA3.1 at 16.7 mmol/L glucose; #P < 0.01 compared with pcDNA3.1 at 2.8 mmol/L glucose, as analyzed by two-tailed paired t tests; n = 7). I: Altered fold change of insulin secretion (calculated as insulin secretion at 16.7 mmol/L glucose divided by insulin secretion at 2.8 mmol/L glucose) for each of the three genes (*P < 0.05). J: Small interfering RNA–mediated knockdown of Park2 was verified by qPCR (**P < 0.01, as analyzed by a one-tailed paired t test; n = 6). K: Park2 deficiency resulted in reduced insulin secretion at stimulatory glucose levels (*P < 0.05, as analyzed by a two-tailed paired t test; n = 6) but did not affect fold change of insulin secretion (L). siNC, negative control siRNA; siPark2, siRNA targeting Park2.
Next, we performed a KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis using WebGestalt (33) to identify biological pathways with enrichment of genes that had both significantly altered expression and DMRs (n = 457) in islets from T2D donors. These genes are significantly enriched for categories including Ribosome, Jak-STAT signaling pathway, Pathways in cancer, Type II diabetes mellitus, and Metabolic pathways (adjusted P < 0.006) (Fig. 4B and Supplementary Table 11), further supporting the hypothesis that epigenetic and transcriptional changes in islets may contribute to altered metabolism and T2D.
We then functionally studied the impact of altered methylation on the transcriptional activity using luciferase assays. We selected TMED6 and KIF3A, two genes where the T2D-associated DMRs cover promoter regions and with an inverse relation between methylation and expression (Supplementary Fig. 1C and Supplementary Table 10). TMED6 has also been linked to insulin secretion (10). Luciferase constructs containing respective promoter sequence were either methylated or mock-methylated and transfected into clonal β-cells. In line with our islet data, complete methylation resulted in reduced transcriptional activity of both promoters, whereas partial methylation only reduced KIF3A transcriptional activity (Supplementary Fig. 1D and E).
We then asked if genes with both DMRs and altered expression in human T2D islets have a functional role in β-cells. Genes were selected for functional follow-up based on having multiple and/or large DMRs and exhibiting differential expression (q <0.05) in diabetic islets (Supplementary Table 10 and Fig. 4A and C) together with a potential role in β-cell function (22,25,27,31). To model the situation in T2D islets, we overexpressed Nr4a3, Pid1, and Socs2 and silenced Park2 in rat clonal β-cells (Fig. 4D–G and J). We then measured insulin secretion at basal (2.8 mmol/L) and stimulatory (16.7 mmol/L) glucose levels. Overexpression of Pid1 and Socs2 resulted in reduced glucose-stimulated insulin secretion (GSIS), and all three overexpressed genes caused a slight increase in basal insulin secretion (Fig. 4H). These changes resulted in decreased fold change of insulin secretion (secretion at stimulatory divided by the secretion at basal glucose levels) in β-cells overexpressing either gene (Fig. 4I), which is in line with what is seen in donors with diabetes (2). Park2 deficiency resulted in reduced GSIS (Fig. 4J and K) but did not affect fold change of insulin secretion (Fig. 4L). These secretory defects were not due to altered insulin content (data not shown).
Replication of Differential DNA Methylation in DMRs
We finally used pyrosequencing to replicate differential DNA methylation of some CpG sites in DMRs annotated to ARX, CADM1, PARK2, PID1, SLC2A2, and SOCS2 (Supplementary Table 3) in an independent cohort (Supplementary Table 4). Here, we could replicate differential methylation in T2D versus control islets (Supplementary Table 12).
Additionally, as a technical replication, we correlated methylation data for CpG sites in islet DMRs (Supplementary Table 3) overlapping with sites analyzed by the Infinium 450K array for the same 14 islet samples. Here, the WGBS and 450K array data showed high reproducibility (Supplementary Table 13).
Discussion
Alternative approaches are needed to end the rapid increase in T2D incidence. This study provides valuable insights into T2D pathology and human biology, including the first comprehensive DNA methylation analysis of 24 million sites in a case-control cohort of human islets. We identified 25,820 DMRs in islets from T2D versus control donors. These cover loci annotated to genes important for islets and T2D pathogenesis as well as novel candidates. The identified DMRs were enriched in both experimentally identified and putative binding motifs for islet-specific transcription factors. Integrating our WGBS data with RNA-seq data further identified novel candidate genes that contribute to impaired insulin secretion.
Intriguingly, seven of the identified DMRs were annotated to PDX1, which encodes a transcription factor of key importance during pancreatic development and in mature β-cells where it regulates insulin expression (34). Additionally, mutations in PDX1 cause maturity-onset diabetes of the young 4 (35), knockout of Pdx1 in β-cells causes diabetes (21), and epigenetic changes of Pdx1 in islets of rats exposed to an impaired intrauterine environment predispose to future diabetes (36). Using a candidate gene approach, we previously found altered methylation and expression of PDX1 in T2D islets (4). One PDX1 DMR in the current study validated these previous data. We also found that glucose directly increases Pdx1 methylation in cultured β-cells (4), and a SNP in PDX1 associated with hyperglycemia alters PDX1 methylation in human islets (37). These data support the hypothesis that epigenetic modifications of PDX1 may contribute to diabetes.
The array-based methods used in previous islet studies only covered ∼450,000 and ∼27,000 CpG sites, respectively (2,3), whereas the current study covered ∼24 million sites. Additionally, we studied DMRs, whereas previous studies analyzed individual CpG sites, making it difficult to compare the results among the different studies. However, we found a strong correlation between DNA methylation analyzed with the 450K array and WGBS from the same islet samples. Additionally, regions/sites annotated to 47 genes were found to have differential methylation by all three studies.
The majority of SNPs associated with T2D impact insulin secretion, supporting their importance in islet function (38). However, these SNPs only explain a modest proportion of the estimated heritability of T2D, and it is possible that combinations of genetic and epigenetic variation contribute to disease susceptibility (37,39). This hypothesis is supported by the fact that a large number of our T2D-associated DMRs are located in the same regions as T2D candidate genes identified by genome-wide association studies, such as TCF7L2, ADCY5, KCNQ1, and GLIS3. Functional studies further show the importance of some of these candidate genes, e.g., TCF7L2 and ADCY5, in β-cells (40,41).
DMRs overlapping with enhancer regions and transcription factor binding sites may have important roles in gene regulation and disease development (42). Indeed, these regions were overrepresented in the DMRs we detected. A recent study proposed competition between DNA methylation and binding of transcription factors to DNA (43). In relation to their data, it is worth mentioning that we found 20.9%–43.9% methylation in regions occupied by transcription factors such as PDX1, FOXA2, MAFB, NKX6.1, and NKX2.2. It is possible that these transcription factors bind to DNA in a subset of islet cells with hypomethylated DNA, whereas they may not bind to DNA in the cells with methylated DNA. Methylation may thereby regulate gene expression differentially in different islet cells. The fact that motifs for several islet-specific transcription factors were enriched in T2D-associated islet DMRs further supports an important role of methylation in these regulatory regions.
When combining the identified DMRs with RNA-seq data from islets of T2D and control donors, we identified 457 genes with both altered methylation and expression. For example, SLC2A2 had increased methylation and decreased expression in T2D islets. SLC2A2 encodes GLUT2, which is a major glucose transporter in rodent islets, whereas other glucose transporters have been suggested to be more important in human β-cells (44). Nevertheless, Sansbury et al. (45) reported a loss of function of SLC2A2 as a cause of neonatal diabetes, suggesting a role for GLUT2 in human β-cells. Additionally, we could replicate increased methylation of SLC2A2 in T2D islets in an independent cohort.
To further model human T2D, we performed functional follow-up experiments of additional genes exhibiting DMRs and altered expression in diabetic islets. We overexpressed Nr4a3, which reduces insulin expression by modulating Pdx1 and NeuroD1 expression (46); Socs2, which regulates proinsulin processing and insulin secretion in mice (27); and Pid1, which has been implicated in mitochondrial dysfunction (22). Overexpression of all genes decreased the fold change of insulin secretion, which is in line with what is seen in human diabetic islets. We also silenced Park2 expression, which regulates the mitochondrial control system in β-cells (25). Again, this impaired GSIS. These experiments support the hypothesis that genes identified by WGBS and RNA-seq may contribute to islet dysfunction in T2D.
DNA methylation was initially thought to be a silencing mark; however, recent data show that its function may vary with genomic context and is more complex than initially thought (12). Here, we integrated the islet methylome with histone modifications analyzed by the Roadmap Epigenomics Consortium (9). Notably, Kundaje et al. (9) generated WGBS data in several tissues but only reduced representation bisulfite sequencing methylation data for human islets. Our study is hence the first to provide WGBS coverage in human islets. As may be expected, histone marks enriched around the TSS of actively transcribed genes (H3K9ac, H3K27ac, and H3K4me3) and active enhancers (H3K27ac) had a relatively low methylation level, whereas histone marks enriched around inactive TSS (H3K27me3) or inactive regions (H3K9me3) had a medium or high methylation level. We also found a relatively high methylation level in regions enriched with H3K4me1, which is found at enhancer regions and gene bodies of actively transcribed genes. Thus, the high methylation level in regions enriched for H3K4me1 is in line with the high methylation level that is often found in gene bodies (47). We also integrated the islet methylome and transcriptome and examined the degree of methylation in different genomic regions of protein-coding genes, noncoding RNAs, and pseudogenes. Interestingly, methylation of different transcript types shows different patterns. The relation to gene expression was also distinct when comparing nonexpressed and expressed transcripts, with reduced impact with increasing distance from TSS. Our data clearly show that the relationship between DNA methylation and gene expression depends on transcript type, expression level, and distance from the TSS.
As DMRs may be cell-type specific (20), it is important to consider the impact of cellular heterogeneity. Here, all samples had a high purity, i.e., a high endocrine content, without differences between the groups. Furthermore, we have previously shown that there is no difference in β-cell content in the pancreatic islets from donors with T2D versus control subjects without diabetes (2). These data suggest that the identified T2D-associated DMRs are not due to altered cell composition.
It has been suggested that minimal sequencing requirements in WGBS experiments starts from 5× (48). However, it was also stated that higher coverage is required to detect shorter DMRs with smaller methylation differences. On the basis of our previous findings, we expected some CpG sites to show absolute differences in methylation less than 10% between islets from T2D and control donors (2); hence, we implemented a higher coverage (21×) in our study design and used an average difference in methylation of 5% as a threshold for the DMRs. Moreover, T2D is a polygenic disease where modest effect sizes of a large number of genes are expected to contribute to the disease. Indeed, we discovered numerous DMRs that show absolute methylation differences between 5% and 10%. It should also be noted that some individual CpG sites in a DMR with an average methylation difference of 5% will show bigger differences in methylation.
This comprehensive study identified novel diabetes-related changes in DNA methylation that support an important role for epigenetics in T2D. We need to combine multiple layers of biological information to understand the pathogenesis and progression of T2D. Here, we combined WGBS and RNA-seq data from human islets with known regulatory elements such as histone marks, transcription factor binding sites, and enhancer regions. These integrated data advance our understanding of the etiology of T2D and highlight the importance of epigenetic dysregulation in pancreatic islets and T2D pathogenesis.
Article Information
Acknowledgments. Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala, Sweden. The facility is part of the National Genomics Infrastructure Sweden and Science for Life Laboratory. The SNP&SEQ Platform is also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. Human pancreatic islets were obtained through collaboration with Olle Korsgren at the Nordic Network for Islet Transplantation (Uppsala University, Sweden) and the Human Tissue Laboratory at Lund University Diabetes Centre, coordinated by Ulrika Krus. The authors thank Anna-Maria Veljanovska-Ramsay, of Lund University, for technical assistance and Manolis Kellis, of MIT, Boston, for support with the computational analysis.
Funding. This work was supported by grants from the Swedish Research Council, Region Skåne (ALF), Knut and Alice Wallenberg Foundation, Novo Nordisk Foundation, European Foundation for the Study of Diabetes/Lilly, Söderberg Foundation, Royal Physiographic Society in Lund, Swedish Diabetes Foundation, Påhlsson Foundation, EXODIAB, and Linné grant (B31 5631/2006). Part of the work described in this article was undertaken as part of the 2013-2014 BLUE ScY educational exchange program, which was supported by the Faculty of Medicine at Umeå University and EpiHealth (Lund University).
Duality of Interest. No potential conflicts of interest relevant to this article were reported.
Author Contributions. P.V. performed computational analyses. P.V. and C.L. designed the study. P.V., K.B., T.R., and C.L. wrote the manuscript. K.B., J.K.O., J.L.S.E., and T.R. performed experiments. P.V., K.B., J.K.O., J.L.S.E., and T.R. performed statistical analyses. P.V., K.B., J.K.O., J.L.S.E., L.E., T.R., and C.L. designed experiments and edited the manuscript. P.V., T.R., and C.L. are guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.
Footnotes
This article contains Supplementary Data online at http://diabetes.diabetesjournals.org/lookup/suppl/doi:10.2337/db16-0996/-/DC1.
- Received August 15, 2016.
- Accepted December 28, 2016.
- © 2017 by the American Diabetes Association.
Readers may use this article as long as the work is properly cited, the use is educational and not for profit, and the work is not altered. More information is available at http://www.diabetesjournals.org/content/license.
References
In this Issue
