In addition to genetic predisposition, environmental and lifestyle factors contribute to the pathogenesis of type 2 diabetes (T2D). Epigenetic changes may provide the link for translating environmental exposures into pathological mechanisms. In this study, we performed the first comprehensive DNA methylation profiling in pancreatic islets from T2D and non‐diabetic donors. We uncovered 276 CpG loci affiliated to promoters of 254 genes displaying significant differential DNA methylation in diabetic islets. These methylation changes were not present in blood cells from T2D individuals nor were they experimentally induced in non‐diabetic islets by exposure to high glucose. For a subgroup of the differentially methylated genes, concordant transcriptional changes were present. Functional annotation of the aberrantly methylated genes and RNAi experiments highlighted pathways implicated in β‐cell survival and function; some are implicated in cellular dysfunction while others facilitate adaptation to stressors. Together, our findings offer new insights into the intricate mechanisms of T2D pathogenesis, underscore the important involvement of epigenetic dysregulation in diabetic islets and may advance our understanding of T2D aetiology.
Type 2 diabetes (T2D) has developed into a major public health concern. While previously considered as a problem primarily for western populations, the disease is rapidly gaining global importance, as today around 285 million people are affected worldwide (IDF, 2009). Lifestyle and behavioural factors play an important role in determining T2D risk. For example, experimentally induced intrauterine growth retardation as well as nutrient restriction during pregnancy in rats have been shown to result in development of T2D in offspring (Inoue et al, 2009) while chronic high‐fat diet in fathers programs β‐cell dysfunction in female rat offspring (Ng et al, 2010). In humans, a reduced birth weight together with an accelerated growth in infancy has been associated with impaired glucose tolerance (IGT) in adulthood (Bhargava et al, 2004). The pancreatic islets of Langerhans are of central importance in the development of T2D. Under normal conditions, increasing blood glucose levels after a meal trigger insulin secretion from the pancreatic islet β‐cells to regulate glucose homeostasis. β‐Cell failure marks the irreversible deterioration of glucose tolerance (Cnop et al, 2007b; Tabak et al, 2009) and results in T2D (UKPDSG, 1995). The unbiased genome‐wide search for T2D risk genes (Saxena et al, 2007; Scott et al, 2007; Sladek et al, 2007; Zeggini et al, 2007, 2008) has placed the insulin‐producing β‐cells at centre stage. These approaches have also inadvertently highlighted the complexity of the biological mechanisms critical to T2D development. Most T2D risk genes identified in these genome‐wide association studies (GWAS) affect β‐cell mass and/or function (Florez, 2008). While the majority of studies in the field have characterised diabetes aetiology on the basis of genetics, new findings suggest the potential involvement of epigenetic mechanisms in T2D as a crucial interface between the effects of genetic predisposition and environmental influences (Villeneuve and Natarajan, 2010). Epigenetic changes are heritable yet reversible modifications that occur without alterations in the primary DNA sequence. DNA methylation and histone modifications are the main molecular events that initiate and sustain epigenetic modifications. These modifications may therefore provide a link between the environment, that is, nutrition and lifestyle, and T2D but only few studies so far have documented aberrant DNA methylation events in T2D (Ling et al, 2008; Park et al, 2008).
DNA methylation occurs as 5‐methyl cytosine mostly in the context of CpG dinucleotides, so‐called CpG sites. It is the best‐studied epigenetic modification and governs transcriptional regulation and silencing (for review, see Suzuki and Bird, 2008). Unlike the relatively study genome, the methylome changes in a dynamic way during development, tissue differentiation and aging. Pathologically altered DNA methylation is well described in various cancers (reviewed in Jones and Baylin, 2007) and its role is starting to be revealed in several other diseases such as multiple sclerosis (Casaccia‐Bonnefil et al, 2008), Alzheimer's disease (Mastroeni et al, 2009) and systemic lupus erythematosus (Javierre et al, 2010). About 75% of human gene promoters are associated with CpG islands (CGIs) (Jones and Baylin, 2007; Suzuki and Bird, 2008), which are clusters of 500 bp to 2 kb length with a comparatively high frequency of CpG dinucleotides. They usually harbour low levels of DNA methylation but can become hypermethylated; this CGI hypermethylation was demonstrated to abrogate transcription of tumour suppressor genes during tumourigenesis (Jones and Baylin, 2007). Lately, DNA methylation changes in CpG sites adjoining yet outside of CGIs, so‐called CGI shores (Irizarry et al, 2009), are gaining increased attention. Intriguingly, CpG sites in these shore sequences, in addition to those within CGIs, are proposed to display differential DNA methylation between cancer and normal cells as well as between cells of different tissues (Irizarry et al, 2009).
The goal of the present work was to clarify the hitherto poorly understood connection between DNA methylation and T2D pathogenesis and to determine whether identified epigenetic changes translate into functional effects that impinge on pancreatic β‐cell function. For this, we have explored DNA methylation landscapes in islets isolated from T2D patients and non‐diabetic individuals.
Identification of the T2D‐related differential DNA methylation profile
We performed DNA methylation profiling to analyse the methylomes of freshly isolated islets from 16 human cadaveric donors of similar age, BMI and ethnicity (5 diabetic and 11 non‐diabetic Caucasian donors; Table I). Using electron microscopy (EM), we examined the purity of the islet preparations (see Supplementary data). In diabetic islets, decreases in the percentage of β‐cells have been reported (Sakuraba et al, 2002; Rahier et al, 2008). As shifts in the composition of islet cell types (especially β‐cells that constitute about two thirds of the islet cell mass) with different epigenomes might overlay T2D‐related changes in the cells’ DNA methylation patterns, we used EM to estimate the percentage of β‐cells. The amount of β‐cells in three randomly chosen control islet preparations was 66.3±0.9%. Diabetic islets (n=3, randomly chosen) contained only marginally less β‐cells, accounting for 59.7±1.7% of total islet cell number (cf. Supplementary Table S1, Materials and methods, Supplementary data).
To perform DNA methylation profiling, we used the recently developed Infinium Methylation Assay (Illumina® Infinium® HumanMethylation27 BeadChip; Supplementary Figure S1; Bibikova et al, 2009). This assay interrogates the methylation status of 27 578 CpG sites corresponding to 14 475 consensus coding sequences and well‐known cancer genes (Bibikova et al, 2009). Although the Infinium methylation assay is not a genome‐wide DNA methylation technology, it is a useful screening tool that is sensitive, specific and highly reproducible (Bibikova et al, 2009) allowing for analysis of a defined set of CpG sites in a large number of samples. The analysed CpGs are primarily located in proximal promoter regions (and preferentially inside CGIs) and the array encompasses probes with an average coverage of 1.9 CpG sites per promoter.
As an initial step, we evaluated whether DNA methylation changes in T2D were sufficient to distinguish the diabetic from the control samples when comparing complete methylation profiles. For this, an unsupervised hierarchical clustering was performed which placed the diabetic islet samples as one self‐contained group distinct from the control samples in the resulting dendrogram (Supplementary Figure S2). This outcome highlights two facts: first, diabetic DNA methylation profiles are more similar to each other than to the methylation profile of any control sample, indicating the possibility of a T2D‐specific DNA methylation profile; second, the existence of a single branch containing the five diabetic samples shows that the DNA methylation changes are sufficiently pronounced (even in the unfiltered data sets) to distinguish diabetic from control samples. To further substantiate the correlation of variations of DNA methylation with T2D, we performed a principal component analysis on the data set (Materials and methods; Supplementary Figure S3) and compared the resulting principal components with the T2D state in a point‐biserial correlation (Lowry, 1998–2011). We found that PC #1 and PC #3, together accounting for 37% of the variance in the data set, correlate well with T2D state (Supplementary Figure S3). This is indicative of distinct, T2D‐specific changes in the epigenome of pancreatic islets.
Following this initial data analysis, we identified T2D‐related methylation changes by filtering the data sets for CpG sites showing significant differences in DNA methylation levels between control and T2D groups (cf. Materials and methods and Supplementary Table S2). The results of the filtering are shown as a heatmap (Figure 1A). The depicted methylation profiles discriminate between control samples (left side, indicated by yellow bar above heatmap) and T2D samples (right side, indicated by blue bar). It is already apparent from the above data that there are marked DNA methylation changes in T2D islets. The number of differentially methylated CpG loci in T2D islets is in the same range as reported for other non‐malignant conditions analysed with the same technology platform (19 in T1D‐related nephropathy; Bell et al, 2010); 84 and 360, respectively, in analyses of ageing in cells and tissue specimens (Bork et al, 2009; Rakyan et al, 2010).
We then set out to evaluate the descriptive power of the CpG sites in the filtered data set to differentiate diabetic from non‐diabetic specimens in sample‐wise comparisons. We therefore extracted the methylation values for each sample and performed a supervised clustering (Figure 1B, cf. Materials and methods). As expected, the resulting dendrogram shows that samples group together in two clusters containing exclusively control (CTL, yellow bar) or diabetic (T2D, blue bar) samples, indicating that class identity (CTL, T2D) is the most important separation criterion (Figure 1B, left‐most branch). To assess clustering confidence in an unbiased way and to overcome inherently subjective visual interpretation of the results depicted in the heatmap (Figure 1A), a bootstrapping analysis was carried out after dendrogram computation (cf. Materials and methods). The obtained bootstrap of 0.85 indicates significant statistical support for the bipartite distribution between diabetic and non‐diabetic samples based on the analysis of the CpGs contained in the filtered data set. The occasional high bootstrap values adjoined to sample pairs illustrate similarities in the DNA methylation profiles of these samples. These data demonstrate that human pancreatic islets undergo DNA methylation alterations in T2D that are discernible by means of DNA methylation profiles.
T2D‐related aberrations encompass mostly promoter‐specific DNA hypomethylation
The above experiments enabled us to collect the first comprehensive DNA methylation data set for T2D human islets. We identified 276 CpG sites, affiliated to 254 gene promoters, showing differential methylation between normal and diseased samples (Figure 1C; Supplementary Table S2). Strikingly, 266 of these 276 CpGs (96%) showed decreased methylation levels, while only 10 were hypermethylated (Figure 1C). This unexpected finding contrasts with the well‐known DNA methylation changes observed in cancers, where gene‐specific hypomethylation and hypermethylation are more or less evenly distributed (Jones and Baylin, 2007). With respect to global DNA methylation, cancers generally display hypomethylation (Jones and Baylin, 2007; Tost, 2010), primarily in repetitive DNA.
To test whether the observed T2D‐related changes are gene specific or whether they reflect global hypomethylation in the genome of islet cells, we measured DNA methylation levels of the repetitive LINE‐1 element in control and diabetic samples with bisulphite pyrosequencing (BPS). Analysing DNA methylation of LINE‐1, which makes up ∼20% of human genome, provides an accurate estimate of global DNA methylation changes (Yang et al, 2004). Figure 1D shows that repetitive elements are not differentially methylated in T2D, as substantiated by the strong overlap between CTL and T2D samples, indicating the absence of global hypomethylation in T2D islets.
As an additional quality control, we examined the set of 276 differentially methylated CpG sites for overlap with known single‐nucleotide polymorphism (SNP) positions to be excluded from further data analysis. We found no overlap with the reported 180 potentially problematic CpG sites contained in the Humanmethylation27 array (Bell et al, 2011) and therefore continued our analyses with the full set of 276 CpGs.
BS validation of T2D‐related differential DNA methylation
To corroborate the observed Infinium measurements (cf. Figure 1 and Supplementary Table S2), we applied BPS and in some cases conventional BS to randomly selected, differentially methylated CpG sites. In all 19 cases tested, differential DNA methylation at the respective CpG sites was confirmed by BPS (Figure 2; Supplementary Figure S4). Where implemented, BS also confirmed the DNA methylation profiling data (Figure 2A; Supplementary Figure S4A). Figure 2A depicts an exemplary analysed gene, ALDH3B1, for which the Infinium data were confirmed by BS and BPS. Additional validated genes are shown in Figure 2B (CASP10) and Figure 2C (PPP2R4 alias PP2A). We discovered two differentially methylated CpG sites inside a CGI in the IGF2/IGF2AS locus. The differential methylation of one of the CpGs in this region was tested and confirmed by BPS (Figure 2D). Further examples are shown in Supplementary Figure S4. A direct comparison of methylation percentages obtained by the Infinium Methylation assay and BPS (Figure 2E) yielded a highly positive correlation (Pearson's correlation R=0.927) confirming the validity of the data. BPS analysis of three negative controls constituting high (>90%), intermediate (∼40%) and low (<10%) levels of DNA methylation showed expectedly no differential DNA methylation between sample groups (Supplementary Figure S4M–O). Of note, CpG sites neighbouring the Infinium‐assayed site are often encompassed in the BPS experiments (cf. Figure 2 and Supplementary Figure S4); these adjacent CpG sites also displayed similar DNA methylation levels. Hence, our validation experiments indicate that individual CpGs from the Infinium Methylation array can be used as informative markers for the methylation status of the respective surrounding regions. Taken together, BPS and BS results confirm the methylation values obtained from the Infinium‐based assay, thus validating the T2D‐specific DNA methylation profiles.
Genomic features associated with differential DNA methylation in T2D
Having shown that the observed differential DNA methylation in T2D is gene‐specific rather than global (cf. Figure 1D), we set out to determine whether this aberrant DNA methylation is (a) located within or outside CGIs, (b) prevalent in distinct promoter classes (that are based on CpG density) or (c) correlated with specific regulatory elements.
Previous work that investigated the localisation of differential DNA methylation in cancer and different tissues has suggested that substantial DNA methylation differences occur in CGI shores (Doi et al, 2009; Irizarry et al, 2009). We used bioinformatic tools to determine the distance of differentially methylated CpGs to the nearest CGI. Utilising CGI prediction (Bock et al, 2007), we found that ∼50% of the differentially methylated CpG sites are located >2 kb from the nearest CGI (‘other CpGs’ in Figure 3A and B and Supplementary Table S3). Additionally, about one quarter of CpGs resides in CGI shores (1–2000 bp from a CGI border) while another quarter of CpG sites was located inside CGIs (Figure 3A; Supplementary Tables S2 and S3). A more detailed representation of the CpG distribution reveals the sites of differential methylation inside the CGI shores to be located preferentially close to the CGIs (bar chart in Figure 3B). This distribution of differentially methylated sites inside CGI shores is similar to regions showing differential methylation in cancer (c‐DMRs; Irizarry et al, 2009), between tissues (t‐DMRs; Irizarry et al, 2009) and during differentiation/cell reprogramming (r‐DMRs; Doi et al, 2009). However, an overall comparison of differential methylation locations between these DMRs (Doi et al, 2009; Irizarry et al, 2009) and our data (Figure 3A and B) shows significant disparity, with a predominance of DNA methylation changes >2 kb away from CGIs (‘other CpGs’ in Figure 3A and B and Supplementary Table S3), thus distinguishing the localisation of the DNA methylation profile in T2D islets from those found in tumours, between tissues and in stem cells.
We next analysed the occurrence of these differentially methylated CpG sites in relation to CpG density of the affiliated gene promoters. Saxonov et al (2006) discovered a bipartite distribution of gene promoters with minor overlap between both classes when categorising promoter sequences by means of their CpG content. They discovered that promoters are either relatively depleted of CpG sites (low CpG promoters, LCPs) or enriched in CpG sites (high CpG promoters, HCPs), preferentially around the transcription start site (TSS). Weber et al (2007) introduced a third class of promoters called intermediate CpG promoters (ICPs), to account for the overlap between the classes mentioned above. They also developed precise classification criteria for the three classes, which we utilised in an adapted form: in this study, we considered positions −700 to +500 relative to the TSS for the promoter classification (cf. Materials and methods), since both CpG density and differential DNA methylation are distributed symmetrically around the TSS (Saxonov et al, 2006). Figure 3C shows that most of the differentially methylated CpG sites from T2D islets are located in LCP and ICP class promoters. ICP class promoters have been described as regions of dynamic DNA methylation changes (Weber et al, 2007), while LCP class promoters have seldomly been investigated. Their role as sites of hypomethylation in T2D therefore remains to be explored.
Looking directly at the CpG ratio of the promoter sequences (cf. Materials and methods), the Infinium array resembles the distribution found in a comprehensive set of human promoters (cf. Weber et al, 2007 and Figure 3D, red bars). The CpG ratio of the promoters displaying differential methylation in T2D islets (Figure 3D, blue bars) is clearly distinct from the probes represented on the Infinium array and also contrasts with the CpG distribution found in promoters throughout the genome. Following the distribution of the blue bars in Figure 3D, it becomes apparent that most T2D‐linked differentially methylated promoters are relatively depleted in CpG sites (CpG ratio <0.5) whereas only few CpG‐rich promoters gain or lose methylation in T2D islets.
In terms of genomic features, we third evaluated whether the observed aberrant DNA methylation in T2D islets correlated with specific regulatory elements. The use of computational approaches to extract functional meaning from the annotated genome has gained importance in recent years. Roider et al (2009) underscored the need to first separate gene promoters on the basis of their CpG content before analysing the presence of enriched transcription factor‐binding motifs. These authors found that promoters depleted of CpG sites often contain tissue‐specific transcription factor‐binding sites. Therefore, we first extracted CpG‐depleted promoters by selecting all differentially methylated promoters with a CpG ratio lower than 0.5 (Saxonov et al, 2006). We then used the extracted set of 172 CpG‐poor promoters as a starting point to detect putative transcriptional regulatory signals using the Pscan (Zambelli et al, 2009) software and the TRANSFAC transcription factor motifs database (Matys et al, 2006). As shown in Figure 3E and Supplementary Figure S5A, we identified significant enrichment of GATA transcription factor family binding sites, namely GATA1 and GATA2 binding sites, as well as a common binding site assigned to all GATA proteins (P‐values for enrichment ranged from 2.38 × 10−15 to 3.28 × 10−13; Supplementary Figure S5). We then repeated this analysis using all differentially methylated LCP and ICP promoters (cf. Supplementary Tables S2 and S3) to enlarge the above mentioned set of 172 ‘CpG ratio <0.5’ promoters, yielding similar results. To increase precision and to assess validity of our predictions, we performed a negative control by conducting the detection routine described above with 100 sets of 172 promoters randomly sampled from CpG‐depleted genes (CpG ratio <0.5) that are not differentially methylated but present on the Infinium array. Importantly, the highest significance of enrichment (best P‐value) encountered in the searches with these random gene promoter sets never exceeded the best P‐value obtained with the set of differentially methylated gene promoters (Supplementary Figure S5B). Combined, these data reveal a statistically significant enrichment of GATA transcription factor‐binding motifs in our set of promoters (CpG ratio <0.5) differentially methylated in T2D islets.
Differential DNA methylation observed in T2D islets is not inducible by high glucose
The differential DNA methylation observed in T2D islets might be either causative or instead secondary to the hyperglycaemia inherent to diabetes. Recent reports have shown that transient exposure of aortic endothelial cells to high glucose induced persistent epigenetic changes (DNA methylation, H3K4 mono‐methylation, H3K9/H3K16 acetylation; El‐Osta et al, 2008; Pirola et al, 2011). To determine whether high‐glucose stress modifies DNA methylation in islets, we exposed non‐diabetic human islets to 28 mM glucose for 72 h and subsequently analysed DNA methylation. Similar culture conditions (30 mM glucose, 48 h) have been demonstrated to result in determinate DNA methylation changes in vascular cells (Pirola et al, 2011). We chose CpG sites that displayed alterations in DNA methylation in T2D islets and which we previously confirmed by BPS (cf. Figure 2 and Supplementary Figure S3). Islets from the same donor incubated for the same duration at 6 mM glucose served as control. As shown in Figure 4A, the high‐glucose treatment did not significantly affect DNA methylation of 16 CpGs in 13 gene promoters tested, particularly in relation to the magnitude of DNA methylation change observed in T2D islets (Figure 4A, right). Figure 4B illustrates representative genes, namely SIRT6, GRB10 and CHAC1, for which methylation changes measured between islets exposed to high glucose and control islets were minor compared with the absolute levels of methylation at these positions. To assess global DNA methylation after high‐glucose stress, we analysed the repetitive LINE‐1 element by BPS. As shown above for the gene‐specific methylation changes, we did not find significant alterations in LINE‐1 methylation between islets incubated under normal glucose conditions and those exposed to high‐glucose concentrations (Figure 4C). Together, these findings suggest that the T2D‐related differential DNA methylation is probably not secondary to hyperglycaemia.
The T2D‐related differential methylation pattern is a feature of diseased islets
An important question regarding the differential DNA methylation detected in T2D islets is whether these epigenetic changes are unique for pancreatic islets or whether they are part of a general pattern occurring in different tissues.
To examine this, we analysed peripheral blood leukocytes from 12 T2D patients and compared their DNA methylation profiles with 12 non‐diabetic controls matched for age and BMI (Supplementary Table S4). The two groups had significantly different glycaemia and HbA1c (Supplementary Table S4). We quantified DNA methylation levels of 11 CpG loci in 9 gene promoters that exhibited differential methylation in T2D islets and were validated by BPS (cf. Figure 2 and Supplementary Figure S4). We detected no significant differences in DNA methylation levels between control and T2D blood samples in any of the 11 CpG loci analysed (Supplementary Figure S6). This impelled us to examine DNA methylation broadly, using again the Infinium Humanmethylation27 array already employed for the islets. Differently from the islets, however, both groups displayed very similar DNA methylation profiles (linear regression R=0.9989; Supplementary Figure S7A). As a matter of fact, we detected almost no T2D‐related differential DNA methylation in blood surpassing the cutoff (±15%, P<0.01). Only one CpG locus in the promoter of the CIDEB gene showed significant hypermethylation (+15.9%, P=0.00066). CIDEB influences obesity and liver steatosis and is a negative regulator of insulin sensitivity (Li et al, 2007). Regarding the 276 CpG loci differentially methylated in T2D islets, these showed very limited DNA methylation changes between non‐diabetic and T2D blood cells (Supplementary Figure S7B).
In conclusion, the T2D‐related DNA methylation changes detected in pancreatic islets are essentially absent from whole blood DNA. Indeed, we detected no T2D‐related differential methylation satisfying our significance criterion except for a single CpG in the promoter of CIDEB that, in turn, displays no significant differential DNA methylation in T2D pancreatic islets. These data suggest that (1) the methylation pattern observed in islets is apparently not a general phenomenon; (2) blood is not a suitable surrogate tissue for studying T2D‐related epigenetic changes in pancreatic islets.
Differential DNA methylation can be correlated with changes in gene expression in a subset of genes
Only few studies to date have reported gene expression profiling in human pancreatic islets. An example of such a study is the work by Bhandare et al (2010) that described gene expression in islets from 31 non‐diabetic donors. It was of interest to examine whether the differential DNA methylation observed in our study occurred in promoters of expressed genes identified by Bhandare et al (2010) or whether it was associated with inactive genes that may become transcriptionally activated when hypomethylated. By comparison of Entrez gene IDs, 196 genes (266 probes) from the reported expression array could be matched to our set of differentially methylated genes (for expression data cf. Supplementary Table S2 columns AP ff.). Expression of all matched genes was above background levels, asserting that differential methylation occurs at promoters of genes that are active in islets. As expected, their absolute expression levels covered several orders of magnitude with no significant correlation between expression and promoter methylation level (cf. Supplementary Table S2), that is, highly active genes in islets are not necessarily devoid of DNA methylation in their promoters. Our comparison of methylation and expression data therefore strongly suggests that the observed changes in promoter DNA methylation levels are not restricted to silent or lowly expressed genes but are also occurring in promoters of expressed genes. In T2D islets, we observed hypomethylation in the promoters of these active genes.
A recent study assessed gene expression in different islet cell types including the insulin‐producing β‐cells (Dorrell et al, 2011). A comparison showed that 240 of our 254 genes are covered by the microarray used by these authors. In all, 170 of these genes have a positive presence call in β‐cells. This indicates that the majority of the genes we detected as differentially methylated in T2D islets are expressed in non‐diabetic β‐cells to a sufficient amount to be reliably detected by microarrays, that is, these are genes actively transcribed in β‐cells.
We next examined whether the altered promoter DNA methylation observed in T2D as compared with control islets could be correlated with changes in gene expression. For this purpose, we analysed gene expression in six islet preparations from male islet donors (three non‐diabetic and three T2D) using an Affymetrix HG‐U133A microarray (cf. Materials and methods). Of the 254 genes showing differential DNA methylation in T2D (cf. Figure 2, Supplementary Figure S4 and Supplementary Table S2), 181 had associated probes on the microarray, with a total of 320 matches between Infinium and expression assays (due to multiple methylation or expression probes affiliated to the same genes; Supplementary Table S4). Forty‐one of these genes showed significant differential gene expression between CTL and T2D islets (BH‐adjusted P<0.05), of which about 80% (34 genes) displayed an anticorrelation between changes in DNA methylation and expression, for example, a decrease of DNA methylation coinciding with an increase of expression (Figure 5A and B; Supplementary Table S6). This suggests that, at least for a subset of T2D‐modified genes, changes in the promoter DNA methylation correlate inversely with the transcriptional activity of the gene. Regarding genomic features, the distribution of the differentially methylated CpGs in the promoters of these genes relative to the TSS is similar to the distribution of all 276 differentially methylated CpGs in our study (Supplementary Figure S8A). Also, the distribution of promoter classes is comparable between the 254 differentially methylated genes and the subgroup displaying an inverse correlation of methylation and expression (Supplementary Figure S8B; cf. also Figure 3B and Supplementary Table S3).
The observation that only a subset of genes displays such inverse correlation between DNA methylation and expression changes could be partially explained by the fact that comparisons between promoter DNA methylation and expression of the affiliated genes are often hampered by positioning of array probes. For instance, the Infinium array achieves an average coverage of 1.9 CpG sites per promoter. However, due to the occasionally high number of alternative promoters in one gene, in most cases the Humanmethylation27 array does not interrogate all promoters of a given gene (see Figure 5C for examples). On the other hand, microarray probes are mainly positioned in common exons that are present in most or all alternative transcripts of a gene, thereby potentially masking expression changes in one or few isoforms (Dai et al, 2005).
Consequently, our data are consistent with other studies that were either unable to establish a connection between differential DNA methylation and gene expression (Illingworth et al, 2008) or found a small but significant set of genes with anticorrelation between changes in promoter methylation and gene expression (Rakyan et al, 2008). Our work is thus in agreement with the emerging notion that, for most genes, the methylation–expression relationship cannot be broken down to the simple formula that decreased methylation equals elevated expression (Eckhardt et al, 2006; Illingworth et al, 2008; Suzuki and Bird, 2008).
Biological pathways associated with differential DNA methylation in diabetic islets: at least some of the uncovered differentially methylated genes possess pathophysiological significance for β‐cell function
To determine the biological relevance of the differential DNA methylation in pancreatic islets, we followed a two‐winged strategy consisting of a combined computational and literature approach on one hand, and on the other hand in‐vitro experiments to uncover the functions of genes that were previously unknown in the context of pancreatic islets and/or β‐cells.
For the first approach, we assessed whether the differentially methylated genes have any overlap or other association with known T2D risk genes. Then, we carried out an Ingenuity Pathway Analysis (IPA; Figure 6A) to identify pathways that are epigenetically affected in T2D islets according to our methylation profiling data. This was augmented by a manual search for the differentially methylated genes in scientific literature reporting on the general biology as well as T2D‐related functions of these genes or the pathways they are part of (Figures 6 and 7). For the second approach, we knocked down expression of several genes by RNA interference and tested the functional consequence of their depletion in β‐cells (Figure 8). For two selected genes, we explored their functional role more extensively in isolated β‐cells and human islets (Figure 9).
At the outset, we analysed the differentially methylated genes with regard to their reported functions and the biological pathways they are part of. To date, large‐scale GWAS have uncovered a number of T2D susceptibility loci (Saxena et al, 2007; Scott et al, 2007; Sladek et al, 2007; Zeggini et al, 2007, 2008). Several of these genes have a role in β‐cells. Therefore, we compared our data with the known T2D risk genes and asked whether these loci are also targets of epigenetic dysregulation. None of the established T2D susceptibility loci (Voight et al, 2010) revealed significant differential DNA methylation in our analyses. The only gene previously detected in GWAS as a potential T2D‐associated locus is GRB10 (Di Paola et al, 2006; Rampersaud et al, 2007), an imprinted gene (cf. Supplementary data) that displays hypomethylation in diabetic islets (Supplementary Figure S4L). We did find aberrantly methylated genes with similar biological functions to the known T2D susceptibility genes (Supplementary Table S2). For example, KCNQ1 and KCNJ11 (SNP variants of which are associated with higher T2D incidence; Scott et al, 2007) were not significantly altered in their methylation levels but three other potassium channel genes, KCNE2, KCNJ1 and KCNK16, had changed promoter methylation states (Supplementary Table S2). Other examples of T2D susceptibility loci for which we identified genes with related functions are SLC30A8 and CDKAL1. In our data sets, SLC30A8 was not differentially methylated but two other zinc transporter genes, SLC39A5 and ZIM2, were hypomethylated. For CDKAL1 we found its methylation state unchanged in T2D islets, while its target gene CDK5R1 exhibited pronounced hypomethylation (Supplementary Figure S4C). In summary, although the promoter methylation of established T2D susceptibility loci was unchanged in our profiling approach, other genes with similar biological function (e.g., potassium and zinc transporters) or part of the same regulatory networks (e.g., CDK5R1 in the CDK5 pathway and GRB10 in the insulin signalling pathway) displayed aberrant DNA methylation.
The analyses described above found only few common T2D candidate genes among the differentially methylated genes uncovered in this study. This could imply that T2D pathogenesis in islets is partially mediated by previously unappreciated genes. To decipher their roles in the context of T2D islets, as a first step we performed an IPA to determine which canonical pathways were overrepresented in our set of genes (Figure 6A). Inflammation‐related processes were highly enriched, in particular the acute phase response and IL‐8 signalling. Other enriched pathways, such as apoptosis and death receptor signalling, emphasise the role of β‐cell loss in T2D. Enrichment for pathways involved in metabolism and internal and external cell structure (e.g., actin cytoskeleton and integrin signalling) may be indicative of altered islet function and architecture.
Second, we performed an extensive manual curation according to a previously described β‐cell‐targeted annotation (Kutlu et al, 2003; Ortis et al, 2010). In partial agreement with the IPA, we found these genes to fall into three broad categories: (1) genes related to β‐cell dysfunction and death, (2) genes potentially facilitating the adaptation of the pancreatic islets to the altered metabolic situation in T2D and (3) genes whose role in disease pathogenesis remains to be unearthed (Figure 6B). The adaptation‐related gene category contains few metabolism‐associated genes (e.g., HK1, FBP2; Figure 6B, right part, Figure 7) and many more genes involved in signal transduction or encoding hormones, growth factors (e.g., EGF, FGF1, IGF2/IGF2AS; Figure 7), or transcription factors involved in important regulatory networks (for instance, FOXA2/HNF3B, PAX4 and SOX6) (Figure 6B, right part, Figure 7). In the β‐cell dysfunction and death category, there were hypomethylated genes related to DNA damage and oxidative stress (e.g., GSTP1, ALDH3B1; Figure 7), the endoplasmic reticulum (ER) stress response (NIBAN, PPP2R4, CHAC1), and apoptosis (CASP10, NR4A1, MADD; Figure 6B, left part, Figure 7).
Some genes of interest from the highlighted categories are depicted in Figure 7. Their annotated functions provide possible explanations of how the epigenetic dysregulation of these genes in diabetic islets is connected to T2D pathogenesis. Numerous genes that were identified by our methylation profiling approach have been functionally implicated in insulin secretion. Examination of the available literature on the function of these genes revealed three aspects of insulin secretion with which they interfere: some of these genes influence the expression of the insulin gene, like MAPK1 and SOX6, or its post‐translational maturation, like PPP2R4 (cf. Figure 7 and references therein). Others can deregulate the process of insulin secretion itself (SLC25A5, Ahuja et al, 2007; RALGDS, Ljubicic et al, 2009) or influence synthesis as well as secretion (vitronectin, Kaido et al, 2006). A third group of differentially methylated genes affects (i) signalling processes in the β‐cell leading to insulin secretion or (ii) glucose homeostasis in β‐cells, thereby modulating insulin response upon stimulation. GRB10 (Yamamoto et al, 2008), FBP2 and HK1 (Figure 7) are examples for these genes. Additional genes found in our study have been implicated in the β‐cells' capability to secrete insulin, though the mechanisms have not yet been fully established. The putative functions of these genes indicate a potential epigenetic impact on insulin secretion at multiple levels, namely signalling, expression/synthesis and secretion.
Functional analysis of differentially methylated genes reveals their involvement in specific biological processes in islets and β‐cells
We next sought to provide experimental proof that altering the activity of some of the differentially methylated genes translates into functional effects that impinge on β‐cell viability. We first performed a screening for cellular outcome after knockdown of eight selected differentially methylated genes. β‐Cell survival was examined under basal conditions and after exposure to the saturated fatty acid palmitate or the ER stress‐inducer cyclopiazonic acid (CPA, a specific inhibitor of the sarcoendoplasmic reticulum Ca2+ ATPase pump (SERCA; Pahl et al, 1996)) following silencing of MKNK2 (MAP kinase‐interacting serine/threonine kinase 2 isoform 2), GUCA2B (guanylate cyclase activator 2B), PER2 (period isoform 2), SFRS2IP (splicing factor arginine/serine‐rich 2, interacting protein), NR4A1, BCL2, CHAC1 and NIBAN by RNA interference (Figure 8). Knockdown of BCL2, an antiapoptotic protein (Gurzov and Eizirik, 2011), NR4A1, an immediate‐early response gene, and MKNK2, a kinase that functions downstream of p38 and ERK, controlling cell survival and negatively regulating global protein synthesis, resulted in increased basal cell death in the rat β‐cell line INS‐1E (Figure 8). NR4A1, BCL2 and NIBAN also were antiapoptotic under lipotoxic conditions, as their knockdown sensitised β‐cells to palmitate. Four of the selected genes exerted proapoptotic functions in β‐cells, illustrated by the improved survival following their knockdown and palmitate or CPA exposure. These were the clock gene PER2, GUCA2B, a ligand for guanylate cyclase receptor C that transduces signals for insulin secretion and growth and differentiation, SFRS2IP, a protein required for pre‐mRNA splicing, and CHAC1.
The biological role of two of these differentially methylated genes with previously unknown function in β‐cells, namely NIBAN and CHAC1, was further studied. NIBAN and CHAC1 have been suggested to be part of the ER stress response in other cell types (Sun et al, 2007; Mungrue et al, 2009). Such ER stress response can initiate apoptosis and has been implicated in β‐cell demise in diabetes (Eizirik et al, 2008).
To explore the role of NIBAN and CHAC1 in ER stress in greater depth, human islets were treated with the physiological ER stressors oleate and palmitate (Cunha et al, 2008) or the synthetic ER stressors thapsigargin (THA), tunicamycin (TUN) or brefeldin (BRE). Expression of NIBAN was induced about two‐fold by the saturated fatty acid palmitate (Figure 9A, column 3) but not with oleate (Figure 9A, column 2), which is a less potent inducer of ER stress (Cunha et al, 2008). A similar induction of expression was observed for CHAC1 when islets were treated with palmitate (Figure 9D, column 3) but not with oleate (Figure 9D, column 2). The synthetic ER stressors THA, TUN and BRE also increased NIBAN and CHAC1 mRNA expression in human islets (Figure 9A and D, right panels). THA, which causes ER calcium depletion by blocking the SERCA pump (Pahl et al, 1996), TUN, which blocks glycosylation of nascent proteins (Hickman et al, 1977), and especially BRE, which inhibits ER‐to‐Golgi transport (Misumi et al, 1986), induced NIBAN and CHAC1 gene expression. The magnitude of ER stress induced by these three chemicals and palmitate, as measured by ATF3, CHOP and BiP mRNA expression (Supplementary Figure S9), was closely correlated with the NIBAN and CHAC1 induction.
We next determined the functional role of NIBAN and CHAC1 expression on the outcome of ER stress in β‐cells. For this purpose, we exposed INS‐1E cells to palmitate or CPA. ER stress induced by these agents increased Niban mRNA expression (columns 1, 3, 5 in Figure 9B). As shown in Figure 9B, expression of Niban is efficiently diminished by a specific small interfering RNA (siRNA) (siNiban; compare columns 1 and 2, 3 and 4, 5 and 6). RNAi‐mediated knockdown of Niban increased apoptosis induced by palmitate (Figure 9C; columns 3 and 4) as well as by CPA (columns 5 and 6 in Figure 9C). This is mirrored by an elevated caspase 3 activation (Figure 9C, lower panel). Taken together, the increase of Niban expression during ER stress and augmented apoptosis after its knockdown indicates an antiapoptotic role for Niban during the β‐cell ER stress response.
Chac1 mRNA expression was also induced by palmitate and CPA in INS‐1E cells (Figure 9E). Chac1 knockdown by siRNA (siChac) significantly reduced its expression upon incubation with palmitate (Figure 9E, columns 3 and 4) or with CPA (columns 5 and 6 in Figure 9E). We next analysed whether the knockdown of Chac1 affected ER stress‐induced apoptosis stimulated by palmitate and CPA (Figure 9F). We found that Chac1 knockdown protected β‐cells from apoptosis (cf. columns 3, 4 and 5, 6 in Figure 9F), which was confirmed by a lessened activation of caspase 3 (Figure 9F, lower panel). The induction of Chac1 expression by ER stressors and the fact that its knockdown protected β‐cells from ER stress‐induced apoptosis suggest that Chac1 is involved in apoptosis triggered after the different branches of ER stress response fail to re‐establish ER homeostasis in β‐cells.
These experiments investigated two novel genes with opposite effects on ER stress‐induced apoptosis in human islets, CHAC1 and NIBAN. Expression of both genes is strongly induced by ER stress, elicited by different agents. While NIBAN protects β‐cells against apoptosis, the activation of CHAC1 indicates failure of the processes restoring β‐cell functionality and heralds ER stress‐induced apoptosis (Eizirik et al, 2008; Eizirik and Cnop, 2010).
In conclusion, by looking at genes known to be involved in T2D/islet pathology as well as by loss of function experiments, our analyses indicate that we did not pick up some irrelevant variations in DNA methylation levels but that the uncovered epigenetic changes indeed affect genes vital for β‐cell function and survival. These genes might therefore be important, and provide previously unappreciated players in T2D pathogenesis.
In this study, we performed the first comprehensive DNA methylation profiling of human T2D pancreatic islets identifying 276 differentially methylated CpG sites that are affiliated to 254 genes. These generated data illustrate, for the first time, the epigenetic dimension of T2D in islets and how this is associated with transcriptional alterations. The pathophysiological significance of eight genes found to be differentially methylated in their promoters was validated by RNAi experiments.
It has been suggested that progressively occurring DNA methylation errors lead to diminished gene responsiveness to external stimuli and might thus contribute to the development of T2D (Gallou‐Kabani and Junien, 2005). Our findings of prevalent promoter hypomethylation in T2D islets are indicative of active biological processes involved in adaptation to the diabetic environment as well as biological pathways associated with β‐cell dysfunction and apoptosis (Figures 6B and 7). The functional relevance of some of the differentially methylated genes in β‐cells was documented by screening for β‐cell survival/death following RNAi and subsequent exposure to stresses relevant to T2D (Figure 8). Given the increased evidence that ER stress‐induced apoptosis is one of the mechanisms of β‐cell loss in T2D (Eizirik et al, 2008), it was of interest to further assess the biological functions of two putative ER stress‐related genes that we found to be hypomethylated in T2D islets, namely NIBAN and CHAC1. We observed that these two genes are upregulated by synthetic ER stressors and by the more physiologically relevant saturated fatty acid palmitate in human islets, while knockdown of their expression by specific RNAi demonstrated their modulatory role in apoptosis (cf. Figure 9). While NIBAN protects against ER stress‐induced apoptosis, CHAC1 seems to contribute to cell death. The hypomethylation observed at both genes could be explained by competing proapoptotic and antiapoptotic processes during ER stress response in diabetic islets. NIBAN is a negative regulator of translation initiation factor eIF2α (Sun et al, 2007). Therefore, its hypomethylation may indicate an attempt to re‐establish ER homeostasis by reduction of protein synthesis (Eizirik et al, 2008). Pending the outcome of these attempts, ER stress‐induced apoptosis may be triggered by CHAC1 and other proapoptotic genes.
An important question with regard to epigenetic changes is: are the observed DNA methylation changes reflected in gene activity? By comparing the obtained DNA methylation profiles with microarray gene expression data, we were able to determine that a high proportion of genes in whose promoter T2D‐related differential DNA methylation occurs are actively transcribed in pancreatic islets. A comparison with expression data of islet cell types (Dorrell et al, 2011) showed that most of the differentially methylated genes are expressed in β‐cells. This allowed us to conclude that T2D‐related aberrant DNA methylation partially happens in the promoters of active genes. One has to keep in mind though that the expression studies in islets as well as in the β‐cells analysed non‐diabetic material. We observed mostly DNA hypomethylation in diabetic islets, not infrequently accompanied by elevated gene expression. Therefore, it can be assumed that the T2D‐related hypomethylation leads, in part, to the induction of formerly silent genes.
Regarding differential gene expression in T2D islets, we observed an inverse correlation between differential promoter methylation and differential gene expression for a subset of genes. It is worth mentioning that for a significant proportion of differentially methylated genes, no statistically significant differential expression was observed (compare Supplementary Tables S4 and S5). This may partially be due to the incompatibilities between methylation assays and the design of expression microarrays mentioned above that hamper in‐depth comparisons between methylation and expression. However, for many genes the link between differential methylation and gene activity may be quite complex. Methylation of the cytosine base changes the topology of the major DNA groove, which may affect the binding of many transcription factors and DNA‐binding proteins. Supporting this possibility, binding of two transcription factors whose target genes are differentially methylated in T2D islets has been described as methylation sensitive, namely CTCF (which binds to IGF2/IGF2AS, Bell and Felsenfeld, 2000; Filippova et al, 2005) and YY1 (which binds to ZIM2, Kim et al, 2003). Additionally, we found several genes encoding chromatin‐associated proteins and transcription factor genes differentially methylated in T2D islets (cf. Figure 6B and Supplementary Table S2). Should the expression of these genes and/or their binding motifs be influenced by differential DNA methylation in T2D islets, it might further add to the complexity of the methylation–expression relationship. This could potentially explain the observations made in this study as well as others (Eckhardt et al, 2006; Illingworth et al, 2008; Suzuki and Bird, 2008) that the relationship between DNA methylation and gene expression is rather complex.
In terms of genomic features, we detected T2D‐associated differential DNA methylation mainly in LCP and ICP, while HCPs are underrepresented. Analysing LCP and a subset of ICP genes (CpG ratio <0.5), we discovered GATA family transcription factors that are predicted to regulate a significant subset of these genes. Interestingly, the GATA transcription factor family members are critical regulators in endocrine development, function and pathologies (Viger et al, 2008). The physiological roles of many differentially methylated loci in T2D can be described as genes responding to (external) stimuli and to stress. Of note, Saxonov et al (2006) found that a disproportionately high percentage of genes affiliated to these biological functions possess promoters with a low CpG density. This might indicate a general principle with regard to the promoter class of the differentially methylated gene loci: while in chronic diseases such as T2D and lupus (Javierre et al, 2010) LCP genes are overrepresented, in diseases associated with cellular overgrowth (such as cancer) there is increased prevalence of HCP and relatively few LCP genes (Richter et al, 2009; Martin‐Subero et al, 2009a, 2009b). Further studies are required to test this intriguing possibility.
A key issue is whether the methylation changes we report play a causal role in T2D or are secondary to the diabetic condition. Indeed, the hypomethylation observed in oxidative stress, ER stress and apoptotic pathways may result from chronic exposure to the stressful metabolic environment of T2D, for example, high‐glucose concentrations (Cnop et al, 2005). An interesting example in this respect is CASP10: we found significant hypomethylation in its promoter (Figure 2B) and since caspase 10 is inducible by advanced glycation end products (Lecomte et al, 2004; Obrenovich and Monnier, 2005), this hypomethylation could be indicative of gene activation caused by chronically elevated blood glucose levels and consequently heightened non‐enzymatic glycosylation events. Interestingly, experimental exposure of islets from non‐diabetic donors to high‐glucose concentrations (28 mM) for 72 h did not induce differential DNA methylation in any of the genes that display methylation changes in T2D islets. Even though these findings do not exclude an impact of chronic exposure to stressors like hyperglycaemia on the islet epigenome, they do make it unlikely that the observed alterations in DNA methylation are merely a consequence of relatively short metabolic insults. By inferring from the functions of the differentially methylated genes, it is possible that some of the identified epigenetic changes play a role in the progressive islet dysfunction in T2D, that is, they have potentially been acquired at different time points during pathological decline. Thus, the hypomethylation observed at some genes, like CASP10, may be a consequence of T2D and severe and long‐lasting hyperglycaemia. On the other hand, some genes, for example, those related to insulin secretion, may have obtained alterations in promoter methylation much earlier. For instance, defects in acute insulin response to glucose (AIRg) are among the earliest impairments and even precede the onset of pre‐diabetic IGT (Bogardus and Tataranni, 2002; Fukushima et al, 2004; Leahy, 2005; Bunt et al, 2007). If aberrant AIRg arises from epigenetic aberrations at genes involved in insulin secretion (which is an established function for many genes in our study), these defects should manifest ahead of clinical T2D development. Whether such changes can be designated ultimately causal for the decline into T2D will remain to be proven. Overall, from the data at hand the changed methylation in the promoters of some genes identified in our study might thus be consequential and represent reactions to the diabetic environment. At other genes, the methylation aberrations could be interpreted to play a causal role, driving the islet dysfunction and T2D pathogenesis. Future, large‐scale studies involving multiple stages of T2D development will be needed to elucidate the role of the epigenetic changes in the various stages of T2D pathogenesis. Due to medical ethics, it is impossible to obtain repeated pancreatic biopsies. Therefore, these studies will need to rely on surrogate tissues that remain to be validated. The availability of the presently described human islet methylation profiling will allow future search and validation of surrogate tissues.
However, identification and validation of tissues whose T2D‐related DNA methylation profiles can serve as a proxy for pancreatic islets might prove difficult. The apparent absence of significant T2D‐related differential DNA methylation in blood raises the possibility that T2D‐related epigenetic aberrations are tissue‐specific although more tissues will have to be screened to substantiate this. The finding of almost no differential DNA methylation in blood cells of T2D patients versus the significant changes in pancreatic islets implies the question whether the observed blood–islet difference is attributable to the different lifespan of the cells, for blood cells being days to months while β‐cells have a lifespan of many decades (Cnop et al, 2010). The validity of blood for epigenetic analysis has, however, been established by previous studies that uncovered differential methylation in DNA isolated from whole blood of individuals that were prenatally exposed to famine (Heijmans et al, 2008; Tobi et al, 2009). Further investigations into T2D‐related epigenetic changes in surrogate tissues for pancreatic islets might elucidate their causative role or expose them as consequences of the disease.
A possible confounding factor for the identification of T2D‐related epigenetic profiles is the medication that T2D patients receive and that may influence gene regulation. Histone deacetylase inhibitors (HDACi), for example, have been demonstrated to increase insulin sensitivity in muscle and liver and partially thwart diabetic nephropathy and retinopathy (Christensen et al, 2011). It is possible that diabetes drugs like rosiglitazone, a PPAR‐γ agonist, or metformin will alter gene activity patterns and confound profiling approaches. Adequately powered epigenetic profiling studies of surrogate tissues that consider the patients' medication may yield new insight of relevance for drug‐based T2D therapy.
As acknowledged by McCarthy and Zeggini (2009), the >40 gene variants of T2D susceptibility genes known to date cannot fully explain T2D predisposition. Our study points to the involvement of epigenetic alterations in T2D thus underscoring the previously established contribution of lifestyle habits to its development. Combining the advantages of genome‐scanning techniques and epigenome analyses might pave the way to better comprehend the pathogenesis of T2D. It will be of great interest to examine SNPs in the differentially methylated genes in T2D described in this study since the interplay between SNPs and differential (allele‐specific) DNA methylation has recently been described (Shoemaker et al, 2010). Linked to the topic of allele‐specific DNA methylation, it is noteworthy that a number of the genes found to display differential methylation are also reported to be imprinted (Supplementary data). It could hence be speculated that at least a partial loss of imprinting occurs in T2D islets.
In conclusion, we report the first comprehensive and detailed analysis of epigenetic changes in T2D, specifically an altered DNA methylation profile in the pancreatic islets of T2D patients with a major preponderance of hypomethylation in sequences outside CGIs. These aberrant methylation events affect over 250 genes, a subset of which is also differentially expressed. The dysregulation of these genes in T2D may notably be linked to β‐cell functionality, cell death and adaptation to metabolic stress. Examination of two genes identified by methylation profiling, NIBAN and CHAC1, revealed their biological functions in distinct processes of the ER stress response. Furthermore, our data highlight genes belonging to biological processes whose involvement in T2D is not yet fully understood, such as inflammation and ion transporters/channels/sensors. Importantly, it can be envisaged that the uncovered DNA methylation changes might be, on one part, indicative of reactions of the islet cells to the diabetic condition and on another part, might be causal of T2D. A challenge in the future is to provide further evidence for the primary effects of methylation changes in the diabetic condition. Taken together, our DNA methylation study on human islets thus lays the ground to further unravel the biological complexity of T2D and outlines an unexpected level of epigenetic regulation in islets, which must be taken into account in future studies aiming to understand the pathogenesis of T2D.
Materials and methods
Isolation of pancreatic islets
From September 2004 to November 2009, pancreatic islets of Langerhans were isolated from pancreata of 5 T2D and 11 non‐diabetic male cadaveric donors in Pisa, Italy, with the approval of the local Ethical Committee, and as described previously (Del Guerra et al, 2005). Glucose‐induced insulin secretion was measured as described. The diagnosis of T2D was based on the previously described clinical criteria (ADA, 1997; Genuth et al, 2003).
Islet purity and β‐cell content
Reliable purity assessment for diabetic islets is challenging. In T2D, the degranulated β‐cells contain less insulin and zinc (Ostenson et al, 2006) and the qualitative dithizone assessment (which targets zinc) therefore underestimates T2D islet purity. Hence, we used EM to analyse islet purity in some of the samples used in the methylation profiling after >2 days in culture (Supplementary data; n=3 for T2D and non‐diabetic samples, respectively) as described (Welsh et al, 2005).
After obtaining written informed consent, blood was sampled from 12 male T2D patients and 12 male age‐ and BMI‐matched controls in K2EDTA tubes.
Methylation profiling using the Infinium assay
Genomic DNAs were isolated from pancreatic islets using the Wizard® SV Genomic DNA kit (Promega Corp.) and from 200 μl of blood using the QIAamp DNA Mini kit (Qiagen, Hilden, Germany). In all, 1 μg of genomic DNA was treated with sodium bisulphite using the EZ DNA Methylation™ kit (Zymo Research) according to manufacturer's procedure, respecting the recommended alterations in protocol for consecutive Infinium methylation analysis. Methylation status of 27 578 distinct CpG sites was analysed using the Infinium HumanMethylation27 BeadChip array (Illumina, Inc., San Diego) according to manufacturer's protocol. Data acquisition was done with the Illumina BeadArrayReader and quality control was performed using the Methylation module (version 1.0.5) of the GenomeStudio™ software (version 126.96.36.19906). Data normalisation was omitted due to resulting in essentially unchanged data sets (background normalisation) or due to inapplicability (Loess, quantile, Bayes) to methylation datasets because of their heteroscedasticity (Cancer Genome Atlas Research Network, 2011). For quality control, a standard quality analysis was performed for each array assessing Bisulphite conversion efficiency, hybridisation efficiency and specificity, single base extension rate, target removal as well as staining for negative and non‐polymorphic probes (GenomeStudio Controls Dashboard). Data handling, comparisons and so on were performed with the Methylation module of the GenomeStudio software package (Illumina), MS Excel, R 2.8.0‐2.11.1, Openstat and MicroarrayAnalyse v1.0 (Graessler, 2008).
CGI and promoter class annotation
Annotations for the Infinium HumanMethylation27 provided by Illumina were augmented with respect to (i) the position of the analysed CpG relative to the nearest CGI (inside a CGI, in CGI shore or >2 kb away from an island) and (ii) the promoter class of the gene affiliated to the evaluated CpG (high/intermediate/low CpG content promoter). For all annotations, the human genome build 36.1 (hg18, March 2006) provided the basis.
For classification of the CpG position relative to CGIs, the CGI map provided by Bock et al (2007) (combined epigenetic score >0.5; genome assembly hg18/NCBI36) was used as reference; the CpGs were classified into three categories according to Irizarry et al (2009). Designation of the CpGs is as follows: ‘inside CGI’ if the CpG was inside a CGI, ‘CGI shore’ if the CpG was located within a 2‐kb region around a CGI and ‘other CpG’ otherwise (distance to closest CGI >2 kb).
Promoters of the gene loci affiliated to the analysed CpG sites were classified according to their CpG content. First, we extracted sequences ranging from positions −700 to +500 relative to the TSS from UCSC genome browser database, then calculated the CpG ratio and the GC content of these sequences in sliding windows of 500 with 5 bp offsets. For classification criteria, we followed the definition by Weber et al (2007). In short, promoters were defined as HCPs if at least one 500 bp window contained a CpG ratio >0.75 together with a GC content >0.55 whereas in LCPs no 500 bp window reached a CpG ratio of at least 0.48. All promoters not fitting in either of the above promoter classes were termed ICPs. Five differentially methylated gene promoters (and a total of 54 gene promoters on the Infinium array) could not be classified due to great distance to TSS or lack of annotation.
For unsupervised hierarchical clustering, the data sets were filtered for probes/CpG sites with a detection P‐value of <0.05. CpGs not fitting into this criterion in at least one sample were excluded from further analysis. The clustering was performed with the average β‐values (equals methylation percentage/100) of 24 422 CpGs for each sample using the cluster analysis function in GenomeStudio applying Euclidean distance metrics.
For UPGMA (Unweighted Pair Group Method using Arithmetic averages) clustering, the data sets for computation were assembled as follows: a group‐wise ∣Δβ∣>0.15 and P<0.01 (Mann–Whitney) were set as filtering criteria. In all, 276 probes fitted into these criteria. The methylation percentage of the CpG site corresponding to each probe was extracted for each sample. Then, the methylation values were categorised into 10 equal classes and imported into MEGA4 (Tamura et al, 2007) in which the phylogenetic analysis was conducted. The dendrogram was computed using the UPGMA method applying the ‘number of differences’ model (Sneath and Sokal, 1973). To determine the validity of the sample clustering based on the methylation data, a bootstrap test (10 000 sampling steps) was used to calculate the percentage of replicate trees in which the associated samples clustered together (Felsenstein, 1985). Bootstrap values of 0.7 or higher were considered significant and are shown next to the branches in Figure 1B.
Conventional BS and BPS
In all, 750 ng genomic DNA was subjected to bisulphite conversion using the Epitect® Bisulfite Kit (Qiagen) or the EZ DNA Methylation kit (Zymo Research) according to manufacturer's protocol. Elution of the converted DNA was generally performed with 26 μl elution buffer and 8 μl of the eluted DNA were used as template in subsequent PCRs. To ensure sufficient amount of product, amplifications were generally performed as nested PCRs.
PCR and sequencing primers for BPS were deduced using the PyroMark® Assay Design 2.0 software (Qiagen). Primers for pre‐amplification and conventional BS were designed manually or with the help of BiSearch primer design tool (http://bisearch.enzim.hu) and evaluated using the GeneRunner software (v3.05 Hastings Software, Inc.). Primers were obtained from Eurogentec S.A. or Sigma‐Aldrich Corp. Biotinylated primers were ordered HPLC purified, all other primers desalted. PCR and sequencing primers are listed in Supplementary Table S7.
The pre‐amplification PCR was conducted with primers (see EF, ER primers in Supplementary Table S7) amplifying 400–720 bp spanning the CpG of interest and additionally as many as possible neighbouring CpG sites. CpG sites in the annealing positions of the PCR primers were avoided where possible; otherwise primers were ordered with ambiguities at the respective positions. PCR was conducted with 3 mM MgCl2, 1 mM of each dNTP, 12% (v/v) DMSO, 500 nM of each primer and optionally 500 mM Betaine in heated‐lid thermocyclers under the following conditions: 95°C 3:00; 25 × (94°C 0:30; 51°C 0:40; 72°C 1:30); 72°C 5:00 and cooled afterwards to 10°C.
For conventional BS, nested PCRs were conducted as described above using 35 PCR cycles and a decreased elongation time of 1 min (for PCR primers cf. Supplementary Table S7). PCR products were separated on a 1% agarose gel and single bands were cut and eluted from the gel. Cloning of the nested PCR products was performed with TOPO TA Cloning® kit (Invitrogen Corp.) and the plasmids were sequenced by Genoscreen.
For pyrosequencing, a nested PCR was performed with primers designed by the PyroMark Assay Design software (Qiagen) using the HotStarTaq PCR kit (Qiagen) according to manufacturer's recommendations. Reactions were performed in heated‐lid thermocyclers under the following conditions: 95°C 15:00; 45 × (94°C 0:30; 55°C 0:30; 72°C 0:30); 72°C 10:00 and finally cooled to 8°C. Sample preparation and pyrosequencing reactions were performed with the Pyromark Q24 system (Qiagen). For validation of Infinium assay‐derived DNA methylation by BPS, usually three to five randomly chosen samples from each group (CTL and T2D) were analysed and DNA methylation degrees were averaged.
Microarray gene expression studies
In all, 100 ng total RNA was prepared and hybridised onto Affymetrix Human HG‐U133A chips, according to the protocols described in the Affymetrix GeneChip® Expression Analysis Manual (Affymetrix, Santa Clara, CA). Chips were scanned in an Affymetrix GeneChip Scanner 3000 and their quality verified with the Microarray Analysis Suite 5.0 (Affymetrix) software and functions from the R/Bioconductor affy and affyPLM packages (Gautier et al, 2004; Bolstad et al, 2005) (R version 2.8.0; Bioconductor version 2.3). Raw gene expression data were normalised using the Robust Multiarray Average (RMA) method described by Irizarry et al (2003) implemented in the affy package. Gene expression ratios were calculated using the limma package by fitting a linear model on each gene (Smyth, 2005).
For comparison of differentially methylated loci with expression profiles of non‐diabetic islets (Bhandare et al, 2010), we utilised the extracted expression data (downloaded from ArrayExpress; accession number E‐MTAB‐191; http://www.ebi.ac.uk/microarray-as/ae/).
INS‐1E cell and human islet culture
The rat insulin‐producing INS‐1E cell line (a kind gift from Professor C Wollheim, Centre Medical Universitaire, Geneva, Switzerland) was cultured in RPMI‐1640 (with 2 mM GlutaMAX‐I) containing 5% FBS (Asfari et al, 1992; Cnop et al, 2007a) and used at passages 59–73. Human islets were isolated from 11 organ donors (age 69±6 years; body mass index 26±1 kg/m2) in Pisa, Italy, as described above. The islets were cultured in Ham's F‐10 medium containing 6.1 or 28 mM glucose as previously described (Cunha et al, 2008; Igoillo‐Esteve et al, 2010; Ladriere et al, 2010). The percentage of β‐cells, assessed in dispersed islet preparations following staining with mouse monoclonal anti‐insulin antibody (1:1000, Sigma) and donkey anti‐mouse IgG Rhodamine (1:200, Jackson Immuno Research Europe, Soham, Cambridgeshire, UK), was 53±3%. Palmitate and oleate (Sigma‐Aldrich, Schnelldorf, Germany) were dissolved in 90% ethanol, and used at a final concentration of 0.5 mM in the presence of 1% BSA (Cunha et al, 2008). The chemical ER stressors THA (diluted in DMSO and used at a final concentration of 1 μM), CPA (diluted in DMSO and used at final concentration of 25 μM), TUN (diluted in PBS and used at a final concentration of 5 μg/ml) and BRE (diluted in ethanol and used at a final concentration of 0.1 μg/ml) were obtained from Sigma‐Aldrich. The control condition contained similar dilutions of vehicle.
Assessment of β‐cell death
Cell death was measured using the neutral red kit (Sigma, TOX4) following manufacturer's instructions. Briefly, cells were incubated with 17 mM neutral red for 3 h at 37°C, washed and the dye extracted for absorbance measurement in a spectrophotometer. Quantitative evaluation of INS‐1E cell apoptosis was done by fluorescence microscopy following staining with the DNA‐binding dyes propidium iodide (5 μg/ml) and Hoechst 33342 (5 μg/ml) (Cnop et al, 2007a). Caspase 3 activation was assessed by western blot, as previously described (Gurzov et al, 2009), using anti‐cleaved caspase 3 antibody (1:1000; from Cell Signaling, Beverly, MA, USA).
Genes were knocked down using siRNA. The Niban siRNA was SMARTpool (L‐080179‐01 from Dharmacon, Chicago, IL, USA). Stealth RNAi (Invitrogen, Carlsbad, CA) was used for CHAC1 (RSS324745), MKNK2 (RSS312292), PER2 (RSS329646), BCL2 (RSS340652) and NR4A1 (RSS330510). The siRNAs for SFRS2IP (s161636) and GUCA2B (s133784) were Silencer Select from Ambion (Austin, TX). A negative control (siCTL) of 21 nucleotide duplex RNA with no known sequence homology was obtained from Qiagen. Lipid–RNA complexes were formed in Optimem1 with 1.5 μl Lipofectamine 2000 (Invitrogen) to 150 nM siRNA and added at a final concentration of 30 nM siRNA for transfection as described (Cunha et al, 2008). Transfected cells were cultured for 2 days and subsequently treated. The achieved knockdown was 74.8±10.4% for MKNK2 (P=0.03), 86.4±3.5% for GUCA2B (P=0.003), 34.2±12.9% for PER2 (P=0.098), 58.3±5.9% for SFRS2IP (P=0.04), 72.4±3.2% for CHAC1 (P=0.02), 72.5±3.6% for NR4A1 (P=0.001), 47.2±3.8% for Bcl‐2 (P=0.01) and 50.7±4.8% for NIBAN (P=0.02), as measured at the mRNA level, except for Bcl‐2 that was measured by western blotting.
Poly(A)+ RNA was isolated and reverse transcribed as previously described (Chen et al, 2001). The PCR was done in 3 mM MgCl2, 0.5 μM forward and reverse primers, 2 μl SYBR Green PCR master mix (Qiagen) and 2 μl cDNA. Standards for each gene were prepared using appropriate primers in a conventional PCR. The samples were assayed on a LightCycler instrument (Roche Diagnostics, Mannheim, Germany) and their concentration was calculated as copies per μl using the standard curve (Overbergh et al, 1999). The expression level of the gene of interest was corrected for the expression of the housekeeping gene glyceraldehyde‐3‐phosphate dehydrogenase (Gapdh, for INS‐1E cells) or β‐actin (for human islets). PCR primers are listed in Supplementary Table SXXX. The different treatments utilised in the study did not change expression of the housekeeping gene (Igoillo‐Esteve et al, 2010; data not shown).
Significance of group‐wise differences in DNA methylation profiles was measured by Mann–Whitney rank‐sum test, P<0.01 was considered significant. Taking into account inter‐individual differences in methylation levels and following Illumina Inc. recommendations, a 15% group‐wise difference of methylation levels was set as a cutoff additional to PMann–Whitney<0.01. Correlation between methylation values by Infinium and BPS was computed using Pearson's correlation test. Differences between distributions (CpG localisation, promoter class) were calculated with χ2 goodness‐of‐fit test (R 2.11.1); P‐values were estimated from the resulting χ2 value. Significance of gene expression differences was tested by Bayes moderated t‐test and P‐values were FDR adjusted using Benjamini‐Hochberg method (R package limma; Smyth, 2005); adj. P<0.05 was considered significant. Differences in glucose‐stimulated insulin secretion, methylation as analysed by BPS and gene expression as analysed by RT–qPCR were assessed by Student's t‐test, P<0.05 was considered significant. Data are represented as mean±s.d. unless indicated otherwise.
DNA methylation data sets for pancreatic islets and whole blood have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo) under accession numbers GSE21232 and GSE34008, respectively.
Supplementary data are available at The EMBO Journal Online (http://www.embojournal.org).
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary Table S2
Supplementary Table S5
Supplementary Table S6
Supplementary Table S7
We thank Drs Francoise Fery and Bernard Corvilain for their precious help in the recruitment of T2D patients and healthy controls. This work was supported by grants from the Action de Recherche Concertée (ARC04/09‐309 and ARC 2010‐2015–Eizirik), the Interuniversity Attraction Poles (IUAP P6/28 and P6/40), the F.N.R.S. (1.1.371.10.F, 1.2.136.06.F, 3.4559.05, 3.4.542.09.F, 7.4559.05 and 3.4583.12), the Belgian Télévie, the European Union (FP6 and FP7) projects Eurodia (LSMH‐CT‐2006‐518153), CEED3 (FP7‐HEALTH‐F2‐2008‐223211) and Naimit (HS‐F2‐2009‐241447) and the EFSD/Novo Nordisk Programme 2011.
Author contributions: MV designed experiments, performed research and interpreted data. Infinium methylation assays were done by EC and SD; MV and MD interpreted the Infinium methylation assay data. Bisulphite genomic sequencing and bisulphite pyrosequencing were done by MV, RD and EC; MV, UV, MD and NN conducted bioinformatic and statistical analyses. MB and PM performed the array analysis. DAC conducted gene functional analyses. For clinical samples: selection was done by SDG, PM and MC; preparation was done by SDG, MB, PM; characterisation was done by DC, MIE, SDG, MM, PM and MC. MV, MC and DLE designed experiments and interpreted data. MC, DLE and FF directed the study. MV, MNN, MC, DLE and FF wrote the manuscript.
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2012 European Molecular Biology Organization