Epigenetic associations of type 2 diabetes and BMI in an Arab population.

BACKGROUND: The prevalence of type 2 diabetes (T2D) and obesity has dramatically increased within a few generations, reaching epidemic levels. In addition to genetic risk factors, epigenetic mechanisms triggered by changing environment are investigated for their role in the pathogenesis of these complex diseases. Epigenome-wide association studies (EWASs) have revealed significant associations of T2D, obesity, and BMI with DNA methylation. However, populations from the Middle East, where T2D and obesity rates are highest worldwide, have not been investigated so far.
METHODS: We performed the first EWAS in an Arab population with T2D and BMI and attempted to replicate 47 EWAS associations previously reported in Caucasians. We used the Illumina Infinium HumanMethylation450 BeadChip to quantify DNA methylation in whole blood DNA from 123 subjects of 15 multigenerational families from Qatar. To investigate the effect of differing genetic background and environment on the epigenetic associations, we further assessed the effect of replicated loci in 810 twins from UK.
RESULTS: Our EWAS suggested a novel association between T2D and cg06721411 (DQX1; p value = 1.18 × 10(-9)). We replicated in the Qatari population seven CpG associations with BMI (SOCS3, p value = 3.99 × 10(-6); SREBF1, p value = 4.33 × 10(-5); SBNO2, p value = 5.87 × 10(-5); CPT1A, p value = 7.99 × 10(-5); PRR5L, p value = 1.85 × 10(-4); cg03078551, intergenic region on chromosome 17; p value = 1.00 × 10(-3); LY6G6E, p value = 1.10 × 10(-3)) and one with T2D (TXNIP, p value = 2.46 × 10(-5)). All the associations were further confirmed in the UK cohort for both BMI and T2D. Meta-analysis increased the significance of the observed associations and revealed strong heterogeneity of the effect sizes (apart from CPT1A), although associations at these loci showed concordant direction in the two populations.
CONCLUSIONS: Our study replicated eight known CpG associations with T2D or BMI in an Arab population. Heterogeneity of the effects at all loci except CPT1A between the Qatari and UK studies suggests that the underlying mechanisms might depend on genetic background and environmental pressure. Our EWAS results provide a basis for comparison with other ethnicities.

Background

The Qatari population is one of the understudied Arab populations in type 2 diabetes (T2D) and obesity research, despite the high prevalence of these diseases among Qataris, with an estimated prevalence of ~23 % for T2D [44] and ~42 % for obesity (World Health Organization 2015) [58]. The increased prevalence of T2D and obesity has occurred during a short period of time (2–3 generations), suggesting an important contribution to the disease risk by changing environment and life style factors, whose effects are potentially mediated by the epigenome. Epigenetic modifications are changes that do not alter the primary DNA sequence itself and include DNA methylation, histone modifications, and other changes in chromatin structure that may affect the regulation of gene expression [24, 39, 40]. These epigenetic modifications are thought to provide a link between environmental exposures and clinical phenotypes and are also suspected to contribute to the unexplained heritability of complex diseases [30].

The recent development of genome-wide DNA methylation arrays, and of sequencing technologies, coupled with bisulfide treatment presents novel opportunities to investigate the role of DNA methylation in complex diseases through epigenome-wide association studies (EWASs). Only a small number of EWASs have been published so far on T2D [6, 7, 12, 14, 26, 52, 56] and on obesity or BMI ([1, 2, 15, 16, 18]; Sun.D. et al. [49]; [57]). These studies identified new potential T2D- and obesity-associated genes. Furthermore, efforts have been made to study the effect of DNA methylation on gene expression and on metabolic profiles in order to provide better understanding of disease mechanisms [38, 40]. Most EWASs with T2D and obesity conducted so far were focused on Caucasian populations, and little is known about whether their findings translate to other ethnicities and genetic backgrounds. We performed here the first EWAS for T2D and BMI in an Arab population using 123 individuals from 15 Qatari families.

Previous genetic studies have shown that translatability between populations does not always hold. For example, genetic variants at the PPARϒ locus that associate with T2D in individuals of European descent seem not to exert any effect on T2D risk in the Qatar population [5]. More discrepancies might be expected for epigenetic risk factors, which are additionally under strong environmental influence. To address this question, we investigated whether the effect of the replicated loci was homogeneous between Caucasians and the Qatari population by using data from 810 twins from UK, therefore having different genetic backgrounds and under different life style/environmental pressures.

Methods

The studies were conducted in concordance with the Helsinki Declaration of ethical principles for medical research involving human subjects. The studies were approved by the relevant institutional review boards in Qatar (Institutional Review Board of Weill Cornell Medical College in Qatar, ethical approval numbers 2012–003 and 2012–0025) and in the UK (Guy’s and St. Thomas’ Hospital Ethics Committee). Written informed consent was obtained from every participant in each study.

Qatari family study

The methylation data used in this study was obtained from whole blood (the only easily accessible type of sample) and has been previously described in [59]. Details on whole-genome sequencing can be found in [27]. Briefly, the study group consisted of 123 subjects of Qatari descent from 15 families of various sizes and structures. The dataset included 72 females with mean age 39 (±16.9) years and 51 males with mean age 36.3 (±17.2) years. The average BMI of the females was 28.3 (±6.2) kg/m2 and of the males was 29.2 (±7.2) kg/m2. A total of 30 individuals consisting of 19 females and 11 males were previously diagnosed with T2D, ascertained by the diabetes clinic at Hamad Medical Corporation. T2D subjects were receiving treatment for diabetes, and no other major diseases were reported. DNA methylation profiling was performed by Illumina using the Infinium HumanMethylation450K BeadChip platform and reported in the form of beta values. After quality control and exclusion of individual probes containing SNPs within a region of ±110 base pairs of the CpG site, based on 40× whole-genome sequencing data [27], a total of 468,472 probes were selected for this study. Normalization was carried out using the Lumi:QN+BMIQ pipeline, using the smoothQuantileNormalization method (Additional file 1: Figure S1). Blood cell type coefficients of monocytes, granulocytes, NK-cells, B cells, CD8+-T-cells, and CD4+-T-cells were estimated from the methylation data using the method described by Houseman et al. [22].

TwinsUK cohort

The TwinsUK cohort was established in 1992 to recruit monozygotic and dizygotic twins [34]. More than 80 % of participants are healthy female Caucasians (age range from 16 to 98 years old). The cohort includes more than 13,000 twin participants from all regions across the United Kingdom, and many have had multiple visits over the years. The TwinsUK cohort has been used in many epidemiological studies and is representative of the general UK population for a wide range of diseases and traits [3]. DNA methylation was measured for 877 individuals randomly selected from the TwinsUK cohort, 810 of who had both BMI and T2D information. All 810 subjects were female Caucasians. The average BMI was 27.8 (±5.2) kg/m2 and 32 individuals were previously diagnosed with T2D. The Infinium HumanMethylation450 BeadChip (Illumina Inc, San Diego, CA, USA) was used to measure DNA methylation. Details of experimental approaches have been previously described [54]. Normalization was carried out using the “minfi” R package [4], a procedure equivalent to the Lumi:QN+BMIQ pipeline. DNA methylation probes that mapped incorrectly or to multiple locations in the reference sequence and probes with detection p value of >0.05 or missing values were removed, resulting in 452,874 probes. Blood cell type coefficients were estimated from the methylation data using the method described by Houseman et al. [22].

Statistical analysis

Associations of T2D and BMI with DNA methylation levels in the Qatari family study and in the TwinsUK cohort were carried out within a variance component framework to model the resemblance among family members. Specifically, the association between the phenotypic traits and DNA methylation levels was evaluated by using a linear mixed model where the total phenotypic variance was partitioned into polygenic and environmental variance, the latter including also measurement errors. DNA methylation levels were modeled as fixed effects, whilst the polygenic and environmental effects were modeled as random components. The phenotypic covariance matrix between subjects was modeled using the matrix of the expected proportion of alleles shared IBD over the genome between each pair of individuals. The significance of the associations was evaluated by comparing the likelihood of a full model including the methylation status in the fixed effect, and the likelihood of a null model where the effect of DNA methylation values was constrained to zero. Age, sex (only for the Qatari dataset), smoking status, and the six Houseman blood cell type coefficients (for B cells, granulocytes, monocytes, natural killer cells, CD8+-T-cells, and CD4+-T-cells) were included in the association model. Additionally, BMI association analysis included T2D status as confounder and vice versa. BMI values were standardized to have zero mean and one standard deviation. Given the limited sample size, and to avoid potential inflation of the association statistics by directly carrying out the study on selected probes, we preliminarily performed an EWAS in the Qatari family sample using the whole set of probes. False discovery rate was evaluated using Storey’s method [47].

Selection of CpG sites for replication

At the time this study was conducted, two large EWASs for T2D [12, 26] and two for obesity and BMI [15, 16] were available. We attempted to replicate the most significant CpG probe for each reported locus that reached genome-wide significance in at least one of these four studies, resulting in eight CpG probes for T2D (Table 1) and 39 for obesity and BMI (Table 2). Conservative Bonferroni method was used to correct for multiple testing, considering an association replicated with T2D if its p value was lower than 6.25 × 10−3 (0.05/8) and with BMI if its p value was lower than 1.28 × 10−3 (0.05/39).

Table 1

Replication of T2D-DNA methylation associations in the Qatari family study. The betas represent the slope of the regression model indicating the rate of change in the dependent variable (trait) as independent variable (methylation beta value) changes. Coordinates are in hg19

Probe	Chr	Position	Gene symbol	Beta	SE	p value	Reference
cg19693031	1	145441552	TXNIP	−2.41	0.57	2.46 × 10−5	[12, 26]
cg00574958	11	68607622	CPT1A	−3.77	2.07	0.068	[12, 26]
cg11024682	17	17730094	SREBF1	1.99	1.30	0.124	[12]
cg09152259	2	128156114	PROC	−1.02	0.71	0.148	[12]
cg06500161	21	43656587	ABCG1	2.17	1.52	0.153	[12]
cg04999691	7	150027050	C7orf29	1.87	1.40	0.182	[12]
cg02650017	17	47301614	PHOSPHO1	−3.22	2.58	0.212	[12]
cg18181703	17	76354621	SOCS3	−0.66	0.91	0.465	[12]

Table 2

Replication of BMI-DNA methylation associations in the Qatari family study. The betas represent the slope of the regression model indicating the rate of change in the dependent variable (trait) as independent variable (methylation beta value) changes. Coordinates are in hg19

Probe	Chr	Position	Gene symbol	Beta	SE	p value	Reference
cg18181703	17	76354621	SOCS3	−10.78	2.34	3.99 × 10−6	[12]
cg11024682	17	17730094	SREBF1	14.56	3.56	4.33 × 10−5	[15]
cg07573872	19	1126342	SBNO2	−9.96	2.48	5.87 × 10−5	[15]
cg00574958	11	68607622	CPT1A	−21.04	5.33	7.99 × 10−5	[15]
cg07136133	11	36422377	PRR5L	−10.43	2.79	1.85 × 10−4	[15]
cg03078551	17	41656298	NA	−19.23	5.85	1.00 × 10−3	[15]
cg13123009	6	31681882	LY6G6E	12.41	3.80	1.10 × 10−3	[15]
cg09349128	22	50327986	NA	−9.60	3.04	1.60 × 10−3	[15]
cg08972190	7	2138995	MAD1L1	11.11	3.54	1.70 × 10−3	[15]
cg06192883	15	52554171	MYO5C	6.68	2.27	3.20 × 10−3	[15]
cg06500161	21	43656587	ABCG1	11.90	4.17	4.30 × 10−3	[12]
cg27243685	21	43642366	ABCG1	11.01	4.11	7.30 × 10−3	[15]
cg06946797	16	11422409	NA	−6.24	2.36	8.10 × 10−3	[15]
cg12992827	3	101901234	NA	−5.60	2.14	8.90 × 10−3	[15]
cg23998749	1	154968781	NA	7.47	2.87	9.30 × 10−3	[15]
cg26354221	22	24822802	ADORA2A	12.11	4.77	0.011	[15]
cg11592786	15	89533581	NA	−28.45	12.46	0.022	[15]
cg26403843	5	158634085	RNF145	3.63	1.61	0.024	[12]
cg26033520	10	74004071	NA	4.92	2.20	0.025	[15]
cg01844514	7	149557121	ZNF862	6.55	3.48	0.060	[15]
cg14017402	2	86225602	NA	3.98	2.35	0.090	[15]
cg22891070	19	46801642	HIF3A	1.64	1.05	0.120	[16]
cg08857797	17	40927699	VPS25	3.59	2.33	0.125	[15]
cg06876354	2	121020189	RALB	7.28	4.76	0.126	[15]
cg15871086	18	56526595	NA	4.08	2.86	0.155	[15]
cg07814318	15	31624584	KLF13	3.64	2.60	0.161	[16]
cg04816311	7	1066650	C7orf50	2.07	1.98	0.295	[15]
cg04927537	17	76976091	LGALS3BP	1.95	1.87	0.297	[15]
cg04869770	1	164561550	PBX1	2.46	2.50	0.325	[15]
cg25178683	17	76976267	LGALS3BP	2.23	2.38	0.350	[15]
cg20954977	2	232260116	B3GNT7	1.65	1.93	0.391	[15]
cg18568872	15	90606494	ZNF710	2.41	3.55	0.497	[15]
cg00863378	16	56549757	BBS2	1.37	2.51	0.584	[15]
cg17560136	8	21915510	EPB49	1.18	2.64	0.654	[15]
cg13708645	12	121974305	KDM2B	0.67	1.71	0.697	[15]
cg15695155	12	121973871	KDM2B	0.81	3.44	0.813	[15]
cg27614723	15	92399897	SLCO3A1	0.81	4.02	0.840	[15]
cg09664445	17	2612406	CLUH	−0.86	5.28	0.871	[16]
cg18307303	5	158757456	IL12B	−0.24	3.39	0.943	[15]

Meta-analysis

Meta-analyses between the Qatari and UK samples were carried out using the Genome-Wide Association Meta-Analysis (GWAMA) software [29]. Specifically, we used a fixed effects model with inverse variance to combine the regression coefficients of each study and their standard errors. Inter-study heterogeneity was estimated by using the Cochran’s Q test and by measuring the proportion of variability that is explained by between-trial heterogeneity (I2 estimates, [21]), both implemented in GWAMA.

Results

Our EWAS (Additional file 2: Figure S2 and Additional file 3: Figure S3) identified one CpG association with T2D that reached genome-wide Bonferroni significance (p value <1.07 × 10−7) (cg06721411 at DQX1; p value = 1.18 × 10−9). No methylation probes were significantly associated with BMI after Bonferroni correction for multiple testing, the strongest association being at cg17501210 (RPS6KA2; p value = 4.90 × 10−7). The full EWAS association data is available as Additional files 4 and 5. Q-Q plots (Additional file 3: Figure S3) of the EWASs for BMI and T2D showed mild inflation of the p value statistics (the genomic inflation factor was 1.10 for T2D and 1.09 for BMI). We also replicated the association of our top T2D CpG, cg06721411 (DQX1), in the TwinsUK cohort (p value = 9.00 × 10−3).

We calculated the heritability of DNA methylation at these probes in the Qatari families. At 1 % FDR, 41,374 (about 10 %) methylation levels showed segregation between family members (median heritability = 0.70; first to third quartile = 0.31–1.00). The replicated loci showed heritability between 0.46 and 0.96, apart from cg11024682 (SREBF1) and cg07573872 (SBNO2), which were not significant at 1 % FDR.

We replicated in the Qatari family study the associations between T2D and cg19693031 (TXNIP; p value = 2.46 × 10−5) (Table 1) and between BMI and cg18181703 (SOCS3; p value = 3.99 × 10−6), cg11024682 (SREBF1; p value = 4.33 × 10−5), cg07573872 (SBNO2; p value = 5.87 × 10−5), cg00574958 (CPT1A; p value = 7.99 × 10−5), cg07136133 (PRR5L; p value = 1.85 × 10−4), cg03078551 (intergenic region on chromosome 17; p value = 1.00 × 10−3), and cg13123009 (LY6G6E; p value = 1.10 × 10−3) (Table 2). Boxplots and scatter plots of these associations are in Figs. 1 and 2. The distributions of the methylation values for these eight CpG sites are in Additional file 6: Figure S4. Although we decided to adopt Bonferroni correction for the replication study, 12 additional associations with BMI showed nominal level of significance and the same direction of the associations as the original EWASs (Table 2).

Fig. 1

Boxplot of methylation beta values at cg19693031 (TXNIP) against diabetes state. The middle lines show the medians of the data while the boxes show the 25th to 75th percentiles (Q1 and Q3). The whiskers extend to include 99 % of the data above while circles represent outliers. Beta value distributions at this probe in the diabetics and healthy individuals showed a difference in the levels of background methylation (Wilcoxon test p value = 2.10 × 10−3)

Fig. 2

Scatter plots of BMI against methylation beta values. Red lines represent the slopes of the regression model. BMI values were corrected for age and sex

The eight replicated associations were analyzed in the TwinsUK cohort, and effects were combined in meta-analyses. The meta-analysis of BMI with the TwinsUK results indicated a moderate presence of heterogeneity between the two studies for cg00574958 (CPT1A; I2 = 56.8 %; Cochran’s heterogeneity statistic’s p value >0.05) and increased the significance of this association to p value = 7.32 × 10−14. On the other hand, a considerable presence of heterogeneity between the two studies was identified for all the other associations (Table 3; Cochran’s heterogeneity statistic’s p value <0.05), despite association at these loci was significant in both populations and with concordant direction (Table 3).

Table 3

Meta-analyses of the replicated loci in the Qatari study with the TwinsUK results

Trait	Probe (gene)	Qatari cohort			TwinsUK cohort			Meta-analysis
		Beta	SE	p value	Beta	SE	p value	Beta (95U/95 L)	SE	p value	I 2 (%)
T2D	cg19693031 (TXNIP)	−2.41	0.57	2.46 × 10−5	−0.34	0.12	6.74 × 10−3	−0.43 (−0.66/−0.20)	0.12	2.71 × 10-4	92.1
BMI	cg18181703 (SOCS3)	−10.78	2.34	3.99 × 10−6	−2.90	0.71	5.45 × 10−5	−3.56 (−4.90/−2.23)	0.68	1.59 × 10−7	90.4
	cg11024682 (SREBF1)	14.56	3.56	4.33 × 10−5	5.86	0.91	2.12 × 10−10	6.39 (4.67/8.12)	0.88	4.28 × 10−13	82.2
	cg07573872 (SBNO2)	−9.96	2.48	5.87 × 10−5	−2.64	0.81	1.16 × 10−3	−3.35 (−4.85/−1.84)	0.77	1.41 × 10−5	87.3
	cg00574958 (CPT1A)	−21.04	5.33	7.99 × 10−5	−12.43	1.90	1.08 × 10−10	−13.40 (−16.91/−9.89)	1.79	7.32 × 10−14	56.8
	cg07136133 (PRR5L)	−10.43	2.79	1.85 × 10−4	−3.82	0.78	8.02 × 10−7	−4.36 (−5.84/−2.89)	0.75	6.43 × 10−9	80.4
	cg03078551 (NA)	−19.23	5.85	1.00 × 10−3	−7.35	1.43	3.36 × 10−7	−8.02 (−10.74/−5.30)	1.39	7.97 × 10−9	74.3
	cg13123009 (LY6G6E)	12.41	3.80	1.10 × 10−3	3.06	1.06	3.81 × 10−3	3.74 (1.73/5.74)	1.02	2.56 × 10−4	82.2

The table reports the results we obtained using a fixed effects model with inverse variance to combine the regression coefficients of each study and their standard errors. P values, effect sizes (beta), and their standard errors (SE) are reported for both studies and for the meta-analysis results. For the meta-analysis, we also reported: upper/lower 95 % CI for beta (beta 95U/95L), and heterogeneity estimates (I 2)

Discussion

The high prevalence of T2D and obesity in Qatar motivated the initiation of genetic and epigenetic research in this country. To the best of our knowledge, this is the first association study of CpG methylation with T2D and BMI in an Arab population. We conducted a full EWAS and attempted to replicate in the Qatari population eight CpG sites associated with T2D (Table 1) and 39 CpG sites associated with BMI (Table 2) in Caucasians. Our EWAS with T2D (Additional file 2: Figure S2 and Additional file 3: Figure S3) identified one significantly associated CpG site at cg06721411 (DQX1; p value = 1.18 × 10−9), while the strongest association with BMI at cg17501210 (RPS6KA2; p value = 4.90 × 10−7) did not reach genome-wide significance. The inflation shown in the Q-Q plots is possibly due to hidden confounders, including potentially reduced folate levels in diabetic subjects. However, the observed inflation is only moderate and does not substantially affect our conclusions. As only cg06721411 (DQX1) in the T2D EWAS achieved Bonferroni significance criteria in our discovery cohort, we also replicated this locus in the TwinsUK cohort (p value = 9.00 × 10−3). The effect was in the same direction. DQX1 (DEAQ box RNA-dependent ATPase 1) is a protein coding gene located on chromosome 2 and is classified as a member of the DEAD/H family. The highest expression of the DQX1 is found in the muscle and liver [25].

Using a conservative Bonferroni correction, we replicated eight of the 47 associations: SOCS3, SREBF1, SBNO2, CPT1A, PRR5L, an intergenic region on chromosome 17, and LY6G6E with BMI; and TNXIP1 with T2D, while nominal significance was reached for a further 12 loci associated with BMI (Table 2). Despite association with methylation at SOCS3 being previously reported for both T2D and BMI [12, 26], only the association with BMI was replicated in our study. However, the association between SOCS3 and T2D was not significant in the study of Chambers and colleagues after adjustment for BMI, suggesting that the observed association with T2D in their study may be driven by differences in adiposity between their T2D cases and controls.

Although the mechanisms linking DNA methylation of SOCS3, SREBF1, SBNO2, CPT1A, PRR5L, and LY6G6E to BMI, and of TXNIP to T2D are not fully established yet, some of these genes have been already functionally linked to metabolic phenotypes. The replicated methylation sites are within proximity of known genes suggesting a regulatory role of the methylation. However, because expression data is not available for this population, we used data from prior studies to confirm the functional relevance of methylation to phenotypes and the expression of these genes. For instance, TXNIP is a pro-apoptotic beta-cell factor and encodes for a protein that acts as a regulator of metabolism and an inhibitor of the antioxidant thioredoxins. A recent study showed that TXNIP is involved in glucose regulation by controlling insulin sensitivity in the periphery of the human body, and its expression is elevated in the skeletal muscles in T2D patients [37] indicating a linkage to phenotype. We observed concordant results in this study: individuals diagnosed with T2D show lower levels of TXNIP methylation, thus suggesting higher TXNIP expression. Also, SOCS3 belongs to the SOCS protein family, which is rapidly induced by cytokines, and acts as an inhibitor of various cytokine signaling pathways. Previous studies have shown that SOCS3 is linked to phenotype by being a negative regulator of leptin [9, 10] and insulin signaling [17, 42, 45]. In addition, there is evidence for association between variants located near SOCS3 with glucose homeostasis, BMI, and other obesity traits [50, 51]. It was also shown that expression of SREBF1 was reduced in adipose and skeletal muscle of diabetic subjects [12, 26]. SREBF1 was shown to regulate carbohydrate metabolism and synthesis in an animal model of obesity and T2D [43]. In another study, qPCR experiments showed that CPT1A expression is correlated with the methylation status of CPT1A gene with p value = 4.1 × 10−14 and replicated in the Framingham Heart Study with p value = 3.1 × 10−13 [23]. Differential methylation at CPT1A was also found to influence gene expression in Dick et al. ([15, 16]; Sun.D. et al. [49]). Also, [55] showed an increase of the suppressor of the cytokine signaling proteins including SOCS3 in the liver, muscle, and fat, in obesity. SOCS3 overexpression in the fat cells was accompanied by glycogen synthesis and activation of glucose uptake. Finally, we also used a recently available database [11] to understand whether our T2D EWAS hit and the eight replicated were associated with any gene expression level. We found that all methylation sites had more than one association with expression (cis-meQTL, cis-eQTM, and/or trans-meQTL) at 5 % FDR. Therefore, differential methylation may suggest a regulatory role.

Interestingly, two CpG sites replicated by our study (CPT1A and TXNIP) and a third CpG we attempted to replicate (ABCG1) were also the only probes significantly associated with alpha-hydroxybutyrate (AHB) in our recent EWAS with blood serum metabolomics traits (cg00574958 in CPT1A, p value = 1.3 × 10−10; cg06500161 in ABCG1, p value = 7.8 × 10−6; cg19693031 in TXNIP, p value = 7.2 × 10−8) [38]. AHB is a sub-product of the ketones metabolism; elevated AHB levels indicate potential insulin resistance [20], and this biomarker is part of the metabolomics-based pre-diabetes test Quantose™[53].

Other biomarkers of T2D and pre-diabetes have been previously associated with methylation at ABCG1, TXNIP, and CPT1A, including 3-methyl-2-oxovalerate [33], glycine [19], several lipid traits, including phosphatidylcholines (PCs) [48], chylomicrons, and their remnants, very-low-density lipoprotein (VLDL) and IDL cholesterol particles [38]. They all showed diabetes-related effect directions that are in agreement with the effect directions observed in this study. Most interestingly, the list of metabolic traits associated with CpGs [38] also includes the product of CPT1A itself, palmitoylcarnitine. Furthermore, higher levels of cg00574958 (CPT1A) methylation were also associated with higher levels of related long-chain fatty acids in the EWAS reported by Petersen et al., including palmitate (16:0), stearate (18:0), and oleate (18:1n9) (see Supplementary Table 5 of [38]). Higher levels of cg06500161 (ABCG1) methylation were also associated in a recent EWAS with higher levels of chylomicrons and VLDL-cholesterol (see Supplementary Table 5 of [38] and Table 4).

Table 4

CpG–metabotype associations at the three replicated loci. Only associations with metabotypes that were significant at P value <1.3 × 10−5 (Bonferroni correction of testing multiple metabolic traits) are shown; effect size (beta′), P value of the linear model, and number of samples (N) (data from Supplementary Table 5 of [38]; for details on this dataset, see there)

Metabolic trait	cg06500161 (ABCG1)			cg00574958 (CPT1A)			cg19693031 (TXNIP)
	Beta′	P value	N	Beta′	P value	N	Beta′	P value	N
1-Oleoylglycerol (1-monoolein)	9.435	2.84 × 10−11	1676	−0.983	4.37 × 10−11	1676	−0.841	1.02 × 10−14	1674
Alpha-hydroxybutyrate (AHB)	1.834	7.80 × 10−6	1749	−0.926	1.30 × 10−10	1749	−0.574	7.24 × 10−8	1747
3-methyl-2-oxovalerate	0.771	1.31 × 10−5	1749	−0.605	4.50 × 10−5	1749	−0.472	6.81 × 10−13	1747
Glycine	−0.531	2.91 × 10−6	1744	4.064	7.58 × 10−9	1744	0.656	4.91 × 10−6	1742
Palmitoylcarnitine	1.729	1.35 × 10−5	1737	−0.849	2.28 × 10−6	1737	−0.579	3.11 × 10−8	1735
PC aa C36:4	2.549	2.82 × 10−5	1781	−0.947	2.07 × 10−8	1781	−0.628	1.66 × 10−6	1779
PC aa C42:0	−0.883	3.91 × 10−8	1781	22.276	3.36 × 10−6	1781	6.181	1.11 × 10−13	1779
PC aa C42:1	−0.844	4.63 × 10−7	1781	12.688	4.18 × 10−5	1781	4.631	4.93 × 10−12	1779
PC ae C44:6	−0.918	8.52 × 10−12	1781	17.438	4.88 × 10−6	1781	5.075	5.48 × 10−13	1779
Chylo-A (nM)	594.7	8.44 × 10−14	1773	−1.000	1.35 × 10−13	1773	−0.996	1.11 × 10−21	1771
Chylo-B (nM)	29.53	7.09 × 10−6	1766	−1.000	3.24 × 10−12	1766	−0.988	8.30 × 10−18	1764
Chylo-Rem (nM)	416.1	1.34 × 10−11	1772	−1.000	2.38 × 10−12	1772	−0.991	5.40 × 10−15	1770
IDL (nM)	30.86	1.32 × 10−9	1773	−0.993	6.00 × 10−7	1773	−0.868	2.22 × 10−7	1771
VLDL-A (nM)	150.7	2.72 × 10−13	1773	−1.000	9.23 × 10−14	1773	−0.985	5.73 × 10−19	1771
VLDL-B (nM)	166.0	4.22 × 10−13	1773	−1.000	9.83 × 10−12	1773	−0.979	9.24 × 10−16	1771

The direction of the associations for all metabolites at these three loci is coherent with the association of CPT1A and TXNIP being in one direction (lower methylation values associated with T2D or obesity) and that of ABCG1 in the opposite one (higher methylation values being associated with obesity). Taken together, these observations support the claim that lower methylation of the CPT1A and TXNIP loci and increased methylation of the ABCG1 locus associate with a well-defined diabetes-specific metabolic phenotype, which is mirrored by the association of the loci with the respective clinical phenotypes, obesity, and diabetes.

Replicated associations identified in this study were also confirmed in the TwinsUK cohort (Table 3). Meta-analysis increased the significance of the associations, but highlighted heterogeneity of the effect sizes for all loci but CPT1A. Some heterogeneity of effects between our results and what was reported in the original papers might be expected, as they could be driven by potential differences in the normalization pipeline of the array data or by the correction of the methylation values using different confounders. However, despite that there were no differences in the normalization pipeline or in the use of confounders between the Qatari family sample and the TwinsUK cohort, we still observed significant effect heterogeneity. This heterogeneity may partly be explained by the different environmental pressures. While the standardized BMI distribution was not different between the two samples (Kolmogorov–Smirnov p value >0.05), the distribution of six out of eight methylation values at the tested probes was different in either location or shape (Kolmogorov–Smirnov p value <0.05; Additional file 6: Figure S4).

There are some limitations to our present study that we are aware of. First, the use of Illumina Infinium HumanMethylation450K arrays targets only a subset of methylation sites across the human genome. Array-based technologies are sensitive to artifacts induced by genetic variants (SNPs) within probe binding sites. This problem is commonly addressed by excluding probes that contain known SNPs, based on the annotations given in the Illumina manifest. However, these annotations are based on SNP tagging technologies and might provide incomplete information, in particular, in less studied non-Caucasian populations. One of the strengths of our present study is the ability to fully remove such potentially confounding genetic effects, based on the availability of whole-genome sequencing data with deep coverage.

Second, DNA methylation was measured using DNA extracted from whole blood that was the only accessible type of sample and may not be representative of more disease-relevant tissues for the diseases under study, such as pancreatic cells and adipose tissues. Studies of methylation in obesity or T2D based on disease-relevant tissues such as the skeletal muscle, adipose tissue, or pancreatic islets are interesting but only exist for relatively small studies [8, 13, 28, 31, 32, 35, 36, 41, 46]. Since DNA methylation can be strongly tissue dependent and as our data was obtained from blood samples, for consistency, we only selected methylation probes for replication from EWAS studies that were also performed in the blood. In addition, the blood consists of various cell types (including B cells, granulocytes, monocytes, natural killer cells, and T cells subset) that may bias methylation estimates. Estimation of cell type coefficients from the methylation data was performed using a method described by [22], and correction for these coefficients in the association model is common practice and it is also applied here.

Conclusions

Given the early state of the epigenome-wide technologies, the number of published EWASs on T2D and obesity so far is small. However, the availability of the technology along with the availability of novel computational tools is expected to accelerate the increase in the number of studies conducted in this field. To the best of our knowledge, this study is the first EWAS of T2D and obesity in an Arab population. Our EWAS identified one new CpG association with T2D at DQX1 that reached genome-wide Bonferroni significance. We also replicated eight previously reported T2D and BMI associations, although they were not genome-wide significant, confirming the relevance of these CpG sites to these phenotypes.