DNA methylation of blood cells is associated with prevalent type 2 diabetes in a meta-analysis of four European cohorts

Background Type 2 diabetes (T2D) is a heterogeneous disease with well-known genetic and environmental risk factors contributing to its prevalence. Epigenetic mechanisms related to changes in DNA methylation (DNAm), may also contribute to T2D risk, but larger studies are required to discover novel markers, and to confirm existing ones. Results We performed a large meta-analysis of individual epigenome-wide association studies (EWAS) of prevalent T2D conducted in four European studies using peripheral blood DNAm. Analysis of differentially methylated regions (DMR) was also undertaken, based on the meta-analysis results. We found three novel CpGs associated with prevalent T2D in Europeans at cg00144180 (HDAC4), cg16765088 (near SYNM) and cg24704287 (near MIR23A) and confirmed three CpGs previously identified (mapping to TXNIP, ABCG1 and CPT1A). We also identified 77 T2D associated DMRs, most of them hypomethylated in T2D cases versus controls. In adjusted regressions among diabetic-free participants in ALSPAC, we found that all six CpGs identified in the meta-EWAS were associated with white cell-types. We estimated that these six CpGs captured 11% of the variation in T2D, which was similar to the variation explained by the model including only the common risk factors of BMI, sex, age and smoking (R2 = 10.6%). Conclusions This study identifies novel loci associated with T2D in Europeans. We also demonstrate associations of the same loci with other traits. Future studies should investigate if our findings are generalizable in non-European populations, and potential roles of these epigenetic markers in T2D etiology or in determining long term consequences of T2D.


Background
Type 2 diabetes (T2D) is a heterogeneous disease with both environmental and genetic factors implicated in disease onset and progression [1,2]. The contribution of environmental factors (e.g. diet, physical inactivity, smoking and ethnicity) in T2D etiology is well-known [3][4][5][6][7][8]. Recent large Genome Wide Association Studies (GWAS) have identified > 400 genetic variants associated with T2D that together explain around 15-18% of T2D estimated heritability [9,10]. However, the pathophysiology linking many of these environmental and genetic factors with disease onset and progression is less well understood. Consequently, there is growing interest in understanding the role of epigenetic mechanisms surrounding T2D [7].
In epidemiological studies, DNA methylation (DNAm) is the most widely studied epigenetic mechanism, partly due to the fact that it can be measured at scale [11][12][13]. DNAm is a heritable mark associated with regulation of gene expression and high-order DNA structure [14]. Because DNAm can be modified in response to lifestyle and environmental factors [15] and is associated with genetic variants [16][17][18], the study of disease-related dysregulation in DNAm could ultimately reveal novel mechanisms in the pathophysiology of T2D, provide new drug targets, and facilitate the discovery of prognostic biomarkers in non-invasive tissues such as peripheral blood [15,19]. Existing evidence demonstrates the association between T2D and DNAm in metabolically relevant tissues [20], and in blood [1,2,14,15,21]. So far, prospective longitudinal studies have discovered up to 18 blood-based CpG sites associated with future liability for T2D [19,22,23]. Many of these CpG sites are also associated with prevalent T2D [1,2,14,21,23,24]. However, most of the latter studies have been conducted in relatively modest sample sizes, rarely using meta-analysis to increase sample size and power to detect differential methylation. To reveal new loci (which could be informative of disease onset or progression) and to confirm previously reported associations, we conducted a large meta-analysis of prevalent T2D using epigenome-wide DNAm from blood samples in four European cohorts. We then implemented functional analyses to investigate possible mechanisms explaining the association between identified markers and prevalent T2D.

Study characteristics
Four cohorts carried out independent EWAS and provided summary statistics for meta-analysis. These were: The Avon Longitudinal Study of Parents and Children (ALSPAC), The Lothian Birth Cohort of 1936 (LBC1936), and two sub-cohorts of the Rotterdam Study (RSIII-1 and RS-Bios) ( Table 1). Of 3,428 total participants, 340 (10%) had diabetes. Individuals with diabetes had on average higher BMI, fasting glucose (or the hemoglobin A1c (HbA1c) in LBC1936), systolic blood-pressure, LDL and triglyceride levels compared with controls (Additional file 1: Table S2). In most cohorts there was little evidence of difference in age, sex and smoking between cases and controls. We identified < 5 (< 0.9%), 37 (11.2%), 41 (6.4%) and 24 (3.9%) new cases of T2D among 558, 329, 643 and 612 controls with available follow-up data in ALSPAC, LBC1936, RSIII-1, and RS-Bios, respectively. New cases were detected after a follow-up period ranging from 5-years (ALSPAC) to 12-years (LBC1936).

T2D is strongly associated with peripheral blood DNA methylation at six CpGs
In the meta-EWAS (model 1, lambda = 1.33), we identified 58 CpGs associated with T2D at p < 1.0 × 10 -5 . Six associations were at p < 1.3 × 10 -7 (Table 2). Associations between T2D and CpGs cg11024682 (SREBF1) and cg18181703 (SOCS3) have been reported previously [22,23,25]. The median absolute difference in effect size amongst the top 58 CpGs from the meta-analysis was 0.8% (range 0.1% to 1.9%) (Fig. 1b), or an absolute difference in 0.2 SDs of DNAm (range 0.0002 SDs to 0.14 SDs) between T2D cases and controls (only in ALSPAC). Adjustment for BMI (model 2) attenuated associations detected in the top six CpG sites identified by p value (i.e., 67% less significant signals). Only the associations at cg19693031 (TXNIP) and cg16765088 (near SYNM) remained significant (p < 1.3 × 10 -7 ) in the BMI model. We observed that among the CpGs attenuated, one was associated with BMI (cg06500161 in ABCG1) and another with waist-circumference (cg00574958 in CPT1A) in minimally adjusted regressions conducted in ALSPAC (see below). Distribution of DNAm by T2D  CpGs discovered in association with incident T2D, and later validated in studies of prevalent T2D by Cardona et al. [23] c CpG identified in association with incident T2D by Chambers et al. [22] d In RS-Bios, estimates of the EWAS for the CpG cg06114363 (ZNF683) did not surpass QC prior to meta-analysis  status for top six CpGs identified in the meta-analysis is presented in the Additional file 1: Figure S1.
We assessed robustness of our findings to interstudy heterogeneity using different strategies. First, we looked at the I 2 statistic, which showed weak heterogeneity with I 2 < 40% in 4/6 T2D-associated CpG sites in the metaanalysis. Heterogeneity was high for CpGs in CPT1A and ABCG1 (I 2 > 70%, p < 0.05) ( Table 2). Second, we used a random-effect meta-analysis adjusted for same covariates as in model 1, finding small changes in the effect estimate compared to the fixed-effect analysis at the T2Dassociated CpGs (absolute change in effect range: 0% to 21%) (Additional file 1: Table S4). In total, 4/6 differentially methylated CpGs remained associated with T2D at p < 1.3 × 10 -7 in the random-effect analysis: cg19693031 (TXNIP), cg00144180 (HDAC4), cg16765088 (near SYNM) and cg24704287 (near MIR23A). Finally, forestplots showed consistency in the direction of association between studies at our six differentially methylated CpGs (see Fig. 2), and the leave-one-out analysis revealed that no single study was consistently having a disproportionately large influence in the combined effect at these sites (Additional file 1: Figure S2).

DMR analyses
We identified 77 regions associated with T2D based on the main meta-analysis (Additional file 1: Table S5). Among these regions, we found an overrepresentation for hypomethylated DMRs (n = 55/77 DMRs) in association with T2D. The DMR with the smallest Sidakcorrected p value overlapping a meta-EWAS signal was identified in an intron of CPT1A (estimate = − 0.01, Sidak p = 1.11 × 10 -9 ). This DMR showed lower methylation values in T2D cases compared with controls. In addition, several DMRs mapped to loci or included CpGs that have been previously associated with prevalent or incident T2D: cg21766592 in SLC1A5 [2], cg14476101 in PHGDH [23] and PFKB3 [23]. Our findings are directionally consistent with these studies.

Association between DNAm and phenotypic traits in ALSPAC at key CpGs identified in EWAS meta-analysis
T2D-associated CpG sites were also associated (at α = 0.05/23 traits analyzed or p < 2.0 × 10 -  Table S6). All six CpGs from the metaanalysis were also associated with white cell-types, and the CpG in MIR23A was exclusively associated with  Table S7-Table S12). Overall, we observed that the six T2D-associated CpGs from the meta-analysis showed directionally concordant associations with fasting glucose, fasting insulin, 2-h glucose and HOMA scores among non-diabetic control participants in ALSPAC ( Table 2, Additional file 1: Table S6).
Odds of T2D in the ALSPAC cohort were calculated for the six CpG sites and one DMR identified in the metaanalysis (Table 3). We observed consistency in the direction and strength of associations in ALSPAC compared to results of the meta-analysis. Using the Nagelkerke's R 2 statistic, we identified small variation in T2D captured by the individual CpG sites (R 2 range 1.3% to 5.7%) or by CpG sites within the DMR. Combining the six differentially methylated CpG sites, we explained 11% of the total Fig. 2 Forest-plots showing results of the meta-EWAS of T2D across four European cohorts for six differentially methylated CpGs identified at p < 1.33 × 10 -7 . For each CpG site we show association estimates (effect-size and 95% confidence interval) of the EWAS conducted by each cohort, and the combined estimate using a fixed-effect inverse-variance weighted meta-analysis (diamond at the bottom). Results were adjusted for age, sex, SVs, cellular heterogeneity and smoking

Table 3 Odds of T2D for six CpG sites and one DMR identified in the meta-EWAS of T2D
Highlighted in bold are associations with p < 0.05. Associations were conducted in a subsample of ALSPAC (n = 1050, N = 48 T2D cases). Main model was adjusted for age, sex, SVs, 6-Houseman cells and smoking, and a second model was additionally adjusted for BMI a Variation calculated using the Nagelkerke's R 2 statistic derived from an unadjusted logistic regression variation in T2D using ALSPAC data. This was similar in LBC1936, where combining the six differentially methylated CpG sites explained 13.6% of the total variation in T2D. For comparison, variation attributed to the combined common risk factors of age, sex, BMI and smoking was R 2 = 10.6% in ALSPAC and 12.2% in LBC1936. Variation attributed to these common risk factors and methylation at six differentially methylated CpG sites combined was R 2 = 22% in ALSPAC and 21.3% in LBC1936. Adding fasting glucose to the predictive model with clinical risk factors and CpG sites captured 68% of the variation in T2D in ALSPAC. We could not perform the same estimation in LBC1936 due to lack of fasting glucose measures and high collinearity between HbA1c, the available glucose trait in this cohort, and T2D.

Pathway enrichment and cross-tissue comparison of DNAm at CpG sites and DMRs associated with T2D
We ran pathway analyses using a total of 80 unique CpG sites, including the six T2D-associated CpGs from the meta-analysis, and 77 index CpGs with the smallest p value identified within each DMR identified. After applying correction for multiple testing, none of the GO terms or KEGG pathways [26] were enriched (Additional file 1: Table S13). For the in-silico comparison of DNAm across tissues, we observed positive correlations between DNAm levels in blood cells and five tissues (r range 0.81 to 0.97, p < 0.05) at the six identified CpGs in the metaanalysis (Additional file 1: Table S14), and at the 80 CpG sites identified across analyses (Additional file 1: Figure  S3). In the enrichment analysis using LOLA [27], we did not identify regulatory elements overlapping with the genomic position of hyper-(n = 23 out of 80) or hypomethylated sets (n = 57 out of 80) associated with T2D. Information from publicly available resources revealed different molecular markers associated with some of our six differentially methylated CpGs from the meta-analysis. Using the BIOS QTL browser [28], we identified expression quantitative trait methylation sites (eQTMs) (at FDR < 0.05) for CpG sites in cg19693031 (TXNIP), cg06500161 (ABCG1) and cg00574958 (CPT1A), where peripheral blood DNAm was inversely associated with gene expression in blood of specific transcripts within the same genes (Additional file 1: Table S15). Using data from the Genetics of DNAm consortium (GoDMC) [29], we found six methylation quantitative trait loci (meQTL) in blood (five in cis and one in trans) associated with five of our six T2D-associated CpGs in cg19693031 (TXNIP), cg06500161 (ABCG1), cg00144180 (HDAC4) (n = 2 meQTL), cg16765088 (SYNM) and cg24704287 (MIR23A) (Additional file 1: Table S16). Using GWAS data, we found weak evidence that the trans meQTL rs6657798 for TXNIP (p = 3.5 × 10 -202 ) was associated with 2-h glucose [30], and that the cis meQTL rs220182 for ABCG1 (p = 3.5 × 10 -202 ) was associated with the homeostasis model assessments for β-cell function (HOMA-B) [31] (Fig. 3, Additional file 1: Table S17). For the SNP identified in common, the same effect allele had opposite effects on DNAm and on the glycemic trait. In both cases, associations between SNPs and traits were observed with unadjusted association p < 0.05. Additionally, no formal causal analysis or colocalization method was applied to establish if an association was present between DNAm at the CpGs and the glucose trait due to a shared causal variant. While of interest, these findings therefore need to be interpreted with caution.

Discussion
Our meta-EWAS of prevalent T2D achieves a sample size of more than double that of previous studies [1,14]. Even though the proportion of cases in our study was relatively low (10%), it was consistent with the percent prevalence of T2D reported in other population-based European studies (i.e., 4.6-16% prevalence) [1,14]. In this study, we found that prevalent T2D was associated with blood DNAm at six individual CpGs in the EWAS and at 77 genomic regions in the DMR analysis, some of them mapping to loci or including CpGs identified in Fig. 3 Overlap of a cis-meQTL for cg06500161 (ABCG1, SNP rs220182) and a trans-meQTL for cg19693031 (TXNIP, SNP rs6657798), with a a GWAS SNP for HOMA-B and b a GWAS SNP for 2-h glucose, respectively. meQTL were retrieved from the Genetics of DNAm consortium (GoDMC, www.godmc .org.uk/) [29] at p < 10 -8 for cis-meQTL (SNPs within 1 Mb from CpG position) and at p < 10 -14 for trans-meQTL (SNPs > 1 Mb or in different chromosomes from CpG position). GWAS SNPs for the glycemic traits were retrieved from the MAGIC consortium (https ://www.magic inves tigat ors.org/) [30,31]. Associations of peripheral blood DNAm with HOMA-B and 2-h glucose were estimated using linear regressions adjusted for age and sex, when appropriate, in two subsamples of diabetes-free individuals in ALSPAC (n = 622 for HOMA-B (only females) and n = 1002 for 2-h glucose) previous studies of prevalent or incident T2D. For CpGs that have been reported in other studies that were not replicated in our meta-analysis, we found that hypermethylation in cg11024682 (SREBF1) and hypomethylation in cg11376147 (SLC43A1), cg18181703 (SOCS3) and cg14476101 (PHGDH), were nominally associated with T2D cases versus controls, in agreement with the original studies [22,23,25]. Additional analyses in the ALSPAC cohort revealed that most of our top six sites from the meta-analysis were also associated with categories of glucose tolerance, age, sex, white cell-types and other clinical phenotypes. Finally, we interrogated publicly available databases to investigate the potential functional role of the CpG sites identified. We observed strong correlation between DNAm levels in blood and in other tissues of relevance for T2D. In addition, we found meQTL, eQTMs and GWAS SNPs for HOMA-B and 2-h glucose related to some of the CpGs detected in the meta-analysis.
For the strongest meta-EWAs association at TXNIP, T2D cases had on average 1.6% lower DNAm levels compared with controls, which represents approximately 3.2% of the overall variance in methylation observed at this site (TXNIP DNAm β-values range: 50-100%). Direction of the association at TXNIP was consistent with previous findings [2,14,[22][23][24][25], but the effect size we identified was slightly smaller than previously reported (effect range between − 3 and − 5%) [14,[23][24][25]. In follow-up analyses in ALSPAC, we demonstrated that per 1% increase in methylation at TXNIP was associated with 6% lower risk of prevalent T2D. This effect was similar to that reported by Chambers et al. [22] for incident T2D. Furthermore, we observed that with respect to controls, participants with diabetes in ALSPAC had 1.2% lower DNAm at TXNIP compared to those with prediabetes. Overall, our results at TXNIP are consistent with previous studies suggesting that hypomethylation of this CpG could be a good prognostic marker for T2D as well as indicating prevalent disease [1,14].
DNAm at cg06500161 annotated at ABCG1 was associated with T2D but not independently of BMI. DNAm at ABCG1 was higher in T2D cases compared with controls, which is consistent with previous evidence at this locus in blood [2,19,22], and in adipose tissue of monozygotic twins discordant for T2D [19,23]. Among diabetes-free participants in ALSPAC, ABCG1 was positively associated with BMI, waist-circumference, fasting insulin, HOMA scores, CRP and triglyceride levels, but negatively associated with HDL levels. Direction of association of ABCG1 methylation with anthropometric [32], glycemic traits [19,[33][34][35][36], HDL [34] and triglyceride levels [19,34,35] was consistent with previous findings in blood. Thus, differences in ABCG1 methylation may appear before T2D, but due to our cross-sectional study design we can't confirm true direction of this association. The lookup in the BIOS QTL browser revealed that DNAm at ABCG1 was inversely correlated with gene expression of the same gene, and this finding agrees with previous observations in blood [22,[34][35][36] and in primary target tissues for T2D [22], suggesting that expression of ABCG1 may be under epigenetic control. No genetic variants predisposing for T2D have been previously reported at ABCG1 in GWAS studies [37]. However, Hidalgo et al. [33] previously identified an association between a cis meQTL for ABCG1 and fasting insulin and the homeostasis model assessment for insulin resistance (HOMA-IR). In our study, we reported a nominal association between a cis meQTL for ABCG1 retrieved from GoDMC, and HOMA-B at GWAS p value < 0.05 (unadjusted p value).
Lower DNAm at cg00574958 in CPT1A was associated with T2D, and this association also appeared to be dependent on BMI. Analyses in ALSPAC showed that per 1% increase in methylation at CPT1A was weakly associated with 21% lower risk of T2D. Direction of the association at CPT1A methylation was consistent with previous findings at this locus [23,24]. We also showed that CPT1A methylation was inversely associated with waist-circumference and triglycerides levels in the continuous analysis, and with BMI, fasting insulin, HOMA scores and CRP in the stratified analysis by quartiles of DNAm using control samples in ALSPAC. Associations of CPT1A with anthropometric traits [32,38] and triglyceride levels [39] were directionally consistent with previous studies in blood. Other associations with cg00574958 in CPT1A have been identified in EWAS of fasting blood lipids [39][40][41], adiponectin [42], the metabolic syndrome [43], and cardiovascular disease (CVD) risk [44]. Using the BIOS QTL browser, we identified an eQTM in blood associated with cg00574958 in CPT1A, indicating that DNAm at this site was inversely associated with expression of a transcript in the same gene, in line with findings by Irvin et al. [39]. CPT1A encodes for the carnitine palmitoyl-transferase 1A enzyme, a protein essential in the oxidative metabolism of fatty acids in the mitochondrion, and in the secretion of glucagon in pancreatic islets [45]. Some evidence of a causal effect of DNAm at CPT1A on T2D was recently demonstrated by Cardona et al. [23] using a Mendelian randomization analysis. However, further evidence is required to validate this result.
We reported three novel CpG sites associated with T2D: cg00144180 (HDAC4), and two intergenic CpGs in cg16765088 (near SYNM) and cg24704287 (near the micro RNA MIR23A). However, the role of these CpG sites in T2D is yet unknown. In disease-free individuals in ALSPAC, these CpGs were also associated with CRP (HDAC4), and with cell-types (HDAC4, MIR23A and SYNM), suggesting a potential role of DNAm at these sites in T2D through inflammatory pathways. The CpG cg24704287 near MIR23A was exclusively associated with cell-types, indicating that this site may be specifically tagging differences in cell composition between T2D cases and controls. In addition, other EWAS have identified an association between HDAC4 and BMI [38], and between cg24704287 (MIR23A) and levels of the soluble Tumor Necrosis Factor Receptor-2 [46] and smoking [47]. Histone deacetylases (HDACs) are enzymes that catalyze the removal of acetyl groups from lysine residues of non-histone and histone proteins, facilitating gene transcription [48]. Animal models have shown that overexpression of HDAC4 causes reduction of β-cell mass, and ongoing clinical trials are evaluating the utility of inhibitors and activators of HDACs in T2D therapy [48]. We also observed weaker associations between T2D and cg11024682 (SREBF1), cg18181703 (SOCS3), cg14476101 (PHGDH) and cg11376147 (SLC43A1) at p values in the order of 10 -5 , helping to confirm previous reports of these CpG sites in studies of incident and prevalent T2D [22,23].
The enrichment analysis did not uncover any potential novel pathways or gene regulatory elements related with T2D-associated hypo-or hypermethylated sets. An in-silico analysis demonstrated high correlation between methylation in blood and five internal target tissues based on CpGs identified in the meta-EWAS (n = 6 sites) and the DMR analysis (n = 77 index CpGs). Therefore, it is unclear whether differences in DNA methylation identified from blood cells are biologically relevant to T2D disease processes, but patterns in blood do appear to represent those seen in disease relevant tissues. Overall, CpG sites identified offer utility in improving our understanding of underlying disease mechanisms. Further validation of these sites in prospective studies will examine them as predictive markers of incident disease, or of adverse complications of T2D.
Results of the model with additional adjustment for BMI revealed a decrease in the effect size and increase in model p value for most of the associations identified in the meta-analysis, but not for associations at cg19693031 (TXNIP) and cg16765088 (near SYNM). Thus, for the remaining CpGs, the association with T2D may have been influenced by underlying differences in BMI between T2D cases and controls. By using a reverse analysis, we were able to estimate the proportion of variation in T2D explained by CpGs in the meta-analysis, which was similar in magnitude to the variation attributed to the common risk factor of age, sex, BMI & smoking, but much lower than that considering fasting glucose as a predictor. In combination, CpGs and clinical risk factors explained up to 68% of the variation in T2D. However, these results should be cautiously interpreted as variation in T2D attributed to the discovered CpG sites was estimated in two cohorts included within the meta-analysis.

Strengths and limitations
Our study benefits from the large sample-size included in the meta-analysis, which allowed us to confirm previous associations with T2D and to report novel findings. By establishing an analysis plan shared across cohorts, we were able to minimize potential sources of technical noise in the data. In addition, we conducted different exploratory methods to test for robustness of our findings to inter-study heterogeneity and reported the association of our top signals with metabolic and anthropometric traits of relevance in T2D. Furthermore, the use of blood as the source of DNAm markers allowed us to identify novel non-invasive signatures of prevalent disease. Future research should explore the application of these signals in T2D incidence, prognosis and monitoring of complications, irrespective of their role in T2D etiology [49]. One of the limitations of this study was the use of a cross-sectional study design, meaning that we cannot establish if observed variation in methylation occurred as a cause or a consequence of T2D. Thus, further studies using a longitudinal approach to assess incident T2D, or implementing causal inference methods such as Mendelian randomization [50][51][52], are required to establish true direction of causality in the association between DNAm and T2D. Another limitation was the use of samples included in the main analysis to estimate proportion of variation in T2D explained by the identified CpGs, which could give biased results. Ideally, T2D variation should be estimated in an independent sample using a weighted methylation score to assess in combination, the predictive ability of the identified markers. Finally, because our sample was restricted to participants of European origin, the generalizability of our novel markers to other populations is still unknown.

Conclusions
We detected cross-sectional associations between blood DNAm and T2D that were consistent with results in previous studies of incident and prevalent T2D in participants of European and non-European ancestry. Novel markers were also identified. Assessment of the etiological role of markers reported in this study may benefit from the use of causal inference methods such as Mendelian randomization, and from the incorporation of tissuespecific analyses. CpG sites and regions identified in this study could potentially be used as prognostic biomarkers for disease onset or complications.

Phenotypic measurement
Prevalent T2D was defined based on medical diagnosis, use of medication to lower blood glucose, fasting blood glucose levels > 7.0 mmol/l or hemoglobin A1c (HbA1c) levels > 6.5%. Further detail of case ascertainment in each study is described in the Additional file 1: section I. Controls were participants without medical diagnosis of T2D, fasting blood glucose < 7.0 mmol/L, HbA1c < 6.5% (when available), and no reported use of glucose-lowering drugs. Age, sex and smoking were extracted from questionnaire data. Missing values for smoking in one of the cohorts were predicted using a methylation score (Additional file 1: section II). Body mass index (BMI) was calculated from study clinic assessed weight and height (kg/ m 2 ).

DNA methylation measurement
As far as possible, cohorts followed a uniform procedure to generate data. Briefly, purified DNA from whole blood samples was bisulfite converted using the Zymo EZ DNA Methylation ™ kit (Zymo Research Corporation, Irvine, USA) and hybridized to the Illumina Infinium HumanMethylation 450 BeadChip (HM450) (Illumina, CA, USA). Samples were excluded if signal detection rate across probes was < 95%. Further detail of quality control steps applied to the methylation data at the probe and sample level, are described in the Additional file 1: section III, Table S1. To control for differences in DNAm arising from cellular heterogeneity, cell proportions for six leucocyte subtypes were calculated from DNAm data [61]. In RS-Bios, direct counts for lymphocytes, monocytes and granulocytes were initially available for all participants. Using directly measured cell counts instead of Houseman cell count estimates in RS-Bios did not impact on the results observed. Cohorts corrected for batch effects using surrogate variable (SV) analysis [62]. Ten SVs were included as covariates in the EWAS provided they were individually independent of T2D status (p > 0.05). Lastly, samples with extreme measures of DNAm were identified using the Tukey method [63] and set to missing values before conducting the EWAS. This latter step involved trimming DNA methylation beta values for each CpG site by removing observations that were more than three times the interquartile range below or above the 25 th and 75 th percentiles, respectively [64]. This method allows to exclude outliers of DNA methylation caused by technical artefacts or rare genetic variants that can bias the analysis [65]. This method has been widely implemented in other studies [65][66][67].

Epigenome-wide association study of prevalent T2D
Each cohort conducted independent EWAS using the meffil R package [68], where DNAm was modeled as the outcome and T2D as the exposure in multivariable linear regression models. The basic model (model 1) was adjusted for age, sex, SVs, cell counts and smoking status, in model 2 we additionally included BMI. The genomicinflation factor (λ) was used to report systematic bias in EWAS results. Before meta-analysis, EWAS estimates were inspected for potential sources of bias using the QCEWAS R package [69] (Additional file 1: section IV). High quality probes obtained in the standardized QC process were included in a fixed-effect inverse variance weighted meta-analysis in METAL (version 2011-03-25) [70]. On average, 368,208 autosomal probes (probe range 349,413 to 376,820) were meta-analyzed. We applied correction for multiple testing using the Bonferroni method, considering associations at p value < 1.33 × 10 -7 . Heterogeneity was evaluated using the I 2 statistic (substantial heterogeneity was defined as I 2 > 40% and heterogeneity p < 0.05), and by visualizing results in forest plots and conducting "leave-one-out" analyses using the metafor R package [71]. To further assess robustness of our findings to inter-study heterogeneity, we conducted a randomeffect meta-analysis using a modified version of METAL [70].
CpGs with p values below an arbitrary p-threshold < 1.0 × 10 -5 were reported. For these CpGs, effect estimates were described using the median absolute difference in methylation β-values between T2D cases and controls. Furthermore, we re-ran analyses for these CpGs using z-values of methylation (mean = 0 and SD = 1) based on data from the ALSPAC study to provide effect estimates standardized by baseline differences in DNAm. Standardized effects were described as the median absolute difference in standard deviations (SDs) of DNAm between groups. We investigated DMR associations with prevalent T2D using comb-p [72] based on summary statistics from the meta-analysis, and using default analysis parameters. All analyses were conducted in R (version 3.3.3) [73].

Association of DNAm with phenotypic traits in ALSPAC
Additional analyses in ALSPAC investigated whether the individual CpG sites identified in the meta-analysis (at p < 1.33 × 10 -7 ) were associated with various clinical phenotypes (Additional file 1: section V). We tested associations using multivariable linear and logistic regressions adjusted for age and sex, where DNAm at the identified CpGs was modelled as the exposure (continuous) against the traits. Analyses were restricted to diabetes-free participants to avoid bias by T2D treatment. We repeated analyses stratifying methylation at each CpG into quartiles to assess robustness of associations. We applied Bonferroni adjustment for multiple testing assuming 23 independent tests (i.e. number of different traits) and α = 0.05. Furthermore, we investigated the proportion of variation in T2D explained by CpG sites discovered in the meta-analysis (p < 1.33 × 10 -7 ), and by CpGs within the single DMR with the smallest corrected p value, using the Nagelkerke's R 2 statistic. This statistic was derived from adjusted logistic regressions between DNAm at the individual CpGs in the meta-EWAS or the average methylation at CpGs within the top DMR as the exposure, and T2D as the outcome (reverse model from the EWAS, Additional file 1: section V).

Pathway enrichment and cross-tissue comparison of DNAm at CpGs identified in the meta-analysis
We conducted pathway analyses using CpGs identified in the meta-analysis (p < 1.33 × 10 -7 ), and the CpG with the smallest p value identified in each DMR (Sidak p value < 0.05). We performed an enrichment analysis for biological pathways using gene ontology (GO) terms and KEGG pathways implemented in the missMethyl R package [26]. p values of enrichment were adjusted for multiple testing using the Bonferroni method (i.e. p < 0.05/# genes in a pathway). For the same list of CpGs, we compared the level of DNAm between blood and five internal tissues relevant to diabetes (liver, skeletal muscle, pancreas, visceral fat or omentum and subcutaneous fat) using a publicly available dataset [74]. We performed an enrichment analysis for regulatory elements using LOLA [27] and identified expression quantitative trait methylation sites (eQTMs) and methylation quantitative trait loci (meQTL) associated with our CpGs from public databases (BIOS QTL browser [28], GoDMC [29]). In addition, we looked up blood-based meQTL associated with our CpGs in expression quantitative trait loci (eQTL) data [75], and in GWAS for T2D [76] and glycemic traits [30,31,77,78] (Additional file 1: section VI) to understand the possible function of the identified CpGs.