Integrated epigenome, whole genome sequence and metabolome analyses identify novel multi-omics pathways in type 2 diabetes: a Middle Eastern study

Background T2D is of high prevalence in the middle east and thus studying its mechanisms is of a significant importance. Using 1026 Qatar BioBank samples, epigenetics, whole genome sequencing and metabolomics were combined to further elucidate the biological mechanisms of T2D in a population with a high prevalence of T2D. Methods An epigenome-wide association study (EWAS) with T2D was performed using the Infinium 850K EPIC array, followed by whole genome-wide sequencing SNP-CpG association analysis (> 5.5 million SNPs) and a methylome-metabolome (CpG-metabolite) analysis of the identified T2D sites. Results A total of 66 T2D-CpG associations were identified, including 63 novel sites in pathways of fructose and mannose metabolism, insulin signaling, galactose, starch and sucrose metabolism, and carbohydrate absorption and digestion. Whole genome SNP associations with the 66 CpGs resulted in 688 significant CpG-SNP associations comprising 22 unique CpGs (33% of the 66 CPGs) and included 181 novel pairs or pairs in novel loci. Fourteen of the loci overlapped published GWAS loci for diabetes related traits and were used to identify causal associations of HK1 and PFKFB2 with HbA1c. Methylome-metabolome analysis identified 66 significant CpG-metabolite pairs among which 61 pairs were novel. Using the identified methylome-metabolome associations, methylation QTLs, and metabolic networks, a multi-omics network was constructed which suggested a number of metabolic mechanisms underlying T2D methylated genes. 1-palmitoyl-2-oleoyl-GPE (16:0/18:1) – a triglyceride-associated metabolite, shared a common network with 13 methylated CpGs, including TXNIP, PFKFB2, OCIAD1, and BLCAP. Mannonate – a food component/plant shared a common network with 6 methylated genes, including TXNIP, BLCAP, THBS4 and PEF1, pointing to a common possible cause of methylation in those genes. A subnetwork with alanine, glutamine, urea cycle (citrulline, arginine), and 1-carboxyethylvaline linked to PFKFB2 and TXNIP revealed associations with kidney function, hypertension and triglyceride metabolism. The pathway containing STYXL1-POR was associated with a sphingosine-ceramides subnetwork associated with HDL-C and LDL-C and point to steroid perturbations in T2D. Conclusions This study revealed several novel methylated genes in T2D, with their genomic variants and associated metabolic pathways with several implications for future clinical use of multi-omics associations in disease and for studying therapeutic targets. Supplementary Information The online version contains supplementary material available at 10.1186/s12916-023-03027-x.


Background
The worldwide prevalence of Type 2 diabetes (T2D) in 2021 was estimated to be 9.8% [1].However, the prevalence in countries from the Middle East and North Africa was 18.1%.The prevalence of T2D in Qataris has been estimated to be 14-17% [2], similar to other Middle Eastern populations.Such high prevalence of T2D in Qatar compared to western populations increases the risk of diabetes complications, including cardiovascular disease, retinopathy, nephropathy and early mortality.More than 400 independent genetic loci have been found to be associated with T2D across different ethnicities and yet explain only a modest portion of the prevalence [3][4][5].T2D is caused by complex interactions of multiple factors including genetic, epigenetic and environmental influences.Current evidence indicates that lifestyle and environmental factors can influence gene expression and clinical phenotype through epigenetic mechanisms such as changes in DNA methylation.Identifying the epigenetic factors driving T2D are thus important to our understanding of T2D etiology.
DNA methylation in blood or pancreatic beta cells has been associated with T2D in multiple populations [6][7][8][9][10][11][12][13].Most used a previous version of the Infinium EPIC array utilizing 450K CpGs (rather than the newer 850K array) and involved western populations.A previous study done on Qataris (n = 123) demonstrated heterogeneity in methylation between Qatari and UK populations, indicating the importance of studying multiple world populations with varying prevalence of T2D [14].
Genetic variants have also been associated with methylation sites, referred to as methylation quantitative trait loci (meQTLs), and have further elucidated factors affecting methylation.For example, 4.7 million cis-and 630,000 trans-meQTL variants were identified in one study using the 450K array [15].Using the 850K EPIC arrays, another study identified a large number of significant SNP-CpG pairs related to various features of diabetic kidney disease [16].Those, and other previous studies, used imputed genotypes from available SNP arrays.
Circulating metabolites have also been associated with T2D and its complications [17,18].Since metabolite levels are influenced by both genetics [19] and epigenetics, and methylation affects gene expression which in turn affects metabolic pathways, T2D methylation would be expected to be associated with metabolites and metabolic pathways related to T2D.This study focuses on the epigenetic associations with T2D in a large number of subjects from the Qatari population to highlight the similarities and differences in T2D methylation compared to other populations.The integration of genomics and metabolomics with epigenetics was used to better understand the biological mechanisms underlying methylation changes associated with diabetes.First, the larger Infinium 850K EPIC array was used to identify novel methylation associations with T2D using Qatari samples.Second, whole genome-wide sequencing (> 5.5M SNPs) was used to identify methylation quantitative trait loci (meQTLs) for the T2D-associated CpGs to highlight genetically driven CpGs.Third, an untargeted panel of metabolites (Metabolon) was used to identify the associations of the T2D CpGs with the metabolome.Finally, combining these results with metabolic correlations, multiomics networks were constructed to find functional metabolic pathways connecting the methylated genes and to understand their biological relation to T2D in this high-risk Middle Eastern population.

Subjects
The Qatar Biobank (QBB) is Qatar's National Repository Centre for biological samples and health information [2].T2D patients were selected based on the following criteria: a) The patient replied "Yes" to the question "Have you ever been diagnosed by a doctor as diabetic;" or b) HbA1c was ≥ 6.5); and c) The selfreported age of diabetes onset was above 30 years or missing.Controls were selected to exclude any diabetics (T2D or T1D) according to the above criteria.
Samples were delivered to Weill Cornell Medicine's genomics core for methylation in three batches as kits became available.Subjects in cohorts 1 (454 samples after QC, 39% with T2D) and 2 (381 samples after QC, 48% with T2D) were selected randomly within disease category with a target to select approximately 40% of the total sample (actual final percent was 43%) with T2D to increase statistical power.A third cohort, Cohort 3 (191 samples after QC, 33% with T2D), was selected after the T2D EWAS analyses of cohorts 1 and 2 were completed (batch 3) and was combined with these cohorts to increase the power of the methylation-SNP (meQTL) association analyses.

DNA methylation profiling, quality control and statistical methods
A total of 1056 blood samples were collected and profiled using an Infinium Methylation EPIC 850K beadchip (Illumina, Inc).The R package and library "minfi" were used for quality control.Samples that failed the "qc" function and sex assignment were removed.Samples or sites with a median detection p value < 0.05 were kept.~ 30,000 CpG sites were removed after dropping CpGs in a SNP locus ("droploci" function).Normalization was done using "Funnorm" [20] and beta values were used for the analysis.A total of 813,660 CpGs in autosomal chromosomes, and 1026 samples remained after QC (details in Additional file 1).
For the T2D EWAS, CpGs associated with T2D in cohort 1 were tested for replication in cohort 2 and those discovered in cohort 2 were tested for replication in cohort 1.Each time, the discovery p-value of 5.8 × 10 -8 and the replication p-value of 0.05/(#EWAS significant CpGs identified in the discovery cohort) were used.This experiment was repeated two times, with and without adjusting for BMI, and all CpGs with significant replication in both experiments were combined.
Linear regression models were used for T2D EWAS associations, with CpG as the dependent variable and T2D the independent variable.Covariates included were sex, two principal components of the cell counts measured in the lab (neutrophils, basophils, eosinophils, monocytes, and lymphocytes), plate number, batch effects, sample well position, smoking surrogate (AHRR CpG site cg05575921) and three whole genome principal components to correct for population stratification [21] and relatedness.Batch effects were determined according to groups of samples that were profiled together, and included differences in recording the gender code in some batches, which was included as a batch-gender interaction effect.Because the random selection of subjects for cohorts 1 and 2 resulted in different distributions of BMI among cases and controls, we used both BMIadjusted and BMI-unadjusted association models.As age was found to be collinear with T2D in our cohorts, an age residual model was used for age adjustment in the T2D EWAS, where the CpGs were first regressed with age in controls and then CpGs in all samples were adjusted using the β0 (intercept) and β1 coefficients resulting from regression.Moreover, the CpGs that were identified as associated with T2D in the EWAS were further tested for association with age at a later stage where a larger set of controls was available.Only two CpGs were found to be associated with age as indicated in Additional file 2: Table S1.
Finally, pathway analysis was done using the Enrichr online tool [22,23], where all genes of the 66 identified CpGs were submitted and KEGG pathways were selected.All KEGG results were used to assign the genes to pathways in order to discuss their biological relevance to diabetes.

Whole genome sequence, quality control and statistical methods (methylation quantitative trait loci)
A total of 5 million variants remained for analysis after quality control (Additional file 1) for 1026 samples.The subjects were divided into discovery (n = 703) and replication cohorts (n = 323) according to the genomic profiling batches.For the meQTL analysis (SNP-CpG), mixed models were used and included kinship, age, sex, T2D status, BMI, two principal components of the actual cell counts (neutrophils, basophils, eosinophils, monocytes, and lymphocytes), plate number, batch effects, sample well position, smoking surrogate (AHRR CpG site cg05575921) and three whole genome principal components as covariates.This analysis was done using the "Genabel" package [24] in the R statistical package, specifically the "polygenic" and "mmscore" functions [25].

Metabolomics, quality control and statistical methods (methylation-metabolite association analysis)
A total of 1160 serum metabolites were measured using an untargeted metabolomics platform by Metabolon Inc. in 2985 samples.Data was quality controlled by log transforming the values, removing outliers (above or below 3 standard deviations from the mean over all samples for each metabolite) and z-scoring the metabolite values.A total of 936 metabolites remained for association with T2D after quality control.For the T2D MWAS analyses, the samples were randomly divided into a discovery (n = 1791, 70%) and a replication cohort (n = 1194, 30%).Additional file 2: Table S12 shows the distribution of gender, age and BMI in the discovery and replication cohorts.
For CpG-metabolite association analyses, 708 samples from combined cohorts 1, 2 and 3, that had both metabolomics and methylation were divided by batch into two sets, and METAL software [26] was used for the metaanalysis.Associations with heterogeneity at p < 0.05 were excluded and corresponded to Isq > 74.
First, a T2D metabolome wide association study (MWAS) was performed to identify associations between metabolites and T2D, using linear regression models, with the metabolite as the dependent variable and T2D as the independent variable, including age, sex, BMI, and three whole genome principal components as covariates.Second, using the metabolites identified as significant from the MWAS and the CpGs significant from the T2D EWAS, CpG-metabolite associations were computed using the same model and covariates (excluding BMI) from the T2D EWAS analyses, replacing T2D status with T2D metabolites.

Mendelian randomization
A two sample mendelian randomization (2SMR) analysis was used where CpG-HbA1c associations were computed in our cohort for CpGs that had meQTL SNP-HbA1c association statistics reported in the GWAS catalogue or Stanford Bio Bank.Linear regression models for CpG-HbA1c associations included all covariates considered in the EWAS analysis except for T2D.CpGs in HK1 and PFKFB2 were found significantly associated with HbA1c using cohort 1 and cohort 2 as discovery and replication cohorts, respectively.The causal associations were then tested for those genes with HbA1c.2SMR statistics were computed using "mr_maxlik" (maximum likelihood) and "ivw" (inverse variance weighted) functions in "Mendeli-anRandomization" package in R.

Network analysis
Metabolites were corrected for all covariates mentioned previously and residuals were used for the partial correlation analysis, using "GeneNet" package in R, and utilizing a Bonferroni p-value threshold.Pairs of all identified CpG associations and significantly correlated metabolites were combined and visualized using Cytoscape software.

Results
Figure 1 summarizes the overall study design.Three cohorts of samples were obtained from the Qatar Biobank over three years with 43.5% being T2D patients (Table 1 shows the cohort characteristics).Cohorts 1 and 2 were used for the T2D Epigenome Wide Association Study (EWAS) (n = 835) and cohort 3 was added to improve the statistical power of SNP-CpG association analysis (n = 1026).Samples from the 3 cohorts that had both metabolomics and methylation were used for the CpG-metabolite association analyses (n = 708).

CpG sites in 48 genes were associated with T2D
A T2D EWAS analysis was performed using cohorts 1 and 2, in a two-way discovery and replication analysis.A total of 47 CpGs were significant at an EWAS significance threshold of p-value p < 5.8 × 10 -8 in cohort 1 and replicated in cohort 2.An additional 10 CpGs were significant when using cohort 2 as the discovery cohort and cohort 1 as the replication cohort.When including BMI as a covariate, 33 CpG sites were identified in cohort 1 and replicated in cohort 2 and 14 CpG sites were identified as significant in cohort 2 and replicated in cohort 1.In total, 74 CpGs from models with or without adjustment for BMI were identified and replicated, out of which 66 CpGs in 48 genes had the same direction of association in both cohorts (Table 2, Additional file 2: Table S1, Fig. 2).Twenty-seven of the 66 CpGs were significant both with and without BMI adjustment.Of the 66 CpGs, only two CpGs in TXNIP were lower in T2D patients compared to nondiabetics; all other CpGs were higher in T2D compared to nondiabetics.
Three of the 66 CpGs (TXNIP, POR, and DQX1) were previously reported at an EWAS significance, and 63 CpGs are thus novel.Thirty one of the 66 identified CpGs have been previously reported to be methylated in T2D or related traits at a nominal p-value (p < 0.05), including fasting glucose, HOMA IR, HbA1c [6,16], or BMI (Additional file 2: Table S2).An additional 15 T2D CpGs were reported to be associated with kidney function, albuminuria and kidney function decline in [16] (Additional file 2: Table S2).We also replicated T2D methylation associations [6-8, 27, 28] at a replication p-value in the larger cohort using the "no-BMI" model (Additional file 2: Table S3).
KEGG pathways containing the genes associated with the significant CpGs included fructose and mannose metabolism (PFKFB2, HK1) and insulin signaling (PPP1R3E, HK1, PRKAR2A).HK1 is linked to galactose, starch and sucrose metabolism, carbohydrate absorption and digestion, glycolysis and gluconeogenesis, and amino sugar and nucleotide sugar metabolism (Additional file 2: Table S5).

whole genome meQTLs in 27 loci were identified for T2D CpGs
The 66 significant CpG sites were tested for association with whole genome sequence SNPs (> 5.5 million SNPs) in 703 samples and tested for replication in 308 samples (see Methods).A total of 688 associations were significant in the discovery cohort at p < 8.7 × 10 -10 (Fig. 3, Additional file 2: Table S6) and constituted 22 unique CpG sites that were associated with 27 unique genes, in 27 unique genic and intergenic loci.Out of 688 SNP-CpG pairs (27 loci), 222 pairs (7 loci) were either replicated at the exact SNP position or in the same locus (within 500 kb of the SNP position) in the Qatari replication cohort, and 369 pairs confirmed previous pairs or loci in Sheng et al. [16] or in BIOS QTL [30].Collectively, 591 pairs (86% of pairs) in 19 loci (75% of loci) were either      S6), and 51 pairs were in novel loci, replicated in Qataris (Table 3, Additional file 2: Table S6).These loci were in TXNIP-SLC2A1, OCIAD1, and SERPINF1 genes.A total of 542 (i.e., 92% of 591) pairs were cis pairs (defined as less than 1 Mb distance between the SNP and the CpG; however, all identified cis pairs were within < 134 kb distance, with an average of 32 kb) and 49 pairs were trans associations (24 pairs were on the same chromosome at > 100 Mb distance, and 25 pairs were on another chromosome) in TXNIP-SLC2A1 and OCIAD1.SNPs in the 591 pairs are annotated as follows: 416 in introns, 84 intergenic, 48 in a downstream/upstream gene, 14 in 3'UTR or 5'UTR, 14 non-coding exons, 6 synonymous, 3 intron-near-splice, and 6 missense and stop-gained SNPs.Three meQTLs harbored nonsynonymous and 5'UTR variants, specifically in SERPINF1, DOCK10 and SLC2A1-TXNIP (Fig. 4).Four missense SNPs in the SER-PINF1-SMYD4 locus were associated with cg11692409 (chromosome 17:1665181), among which rs1136287 (17:1673276) and rs1804145 (17:1674434) associated with this CpG at p = 4.0 × 10 -22 and p = 5.4 × 10 -17 respectively.The former (cg11692409-rs1136287) was also reported in [15] at p = 3.9 × 10 -165 and in [16] at p = 3.8 × 10 -12 .Moreover, rs1136287 had been associated with pigment epithelium-derived factor (PEDF, p = 3.9 × 10 −35 ), where PEDF was associated with diabetic nephropathy (p < 0.001) and sight threatening diabetic retinopathy (p < 0.001) [31].Additionally, [32] showed an inverse correlation between cg11692409 and the mRNA expression of its gene (PEDF) (correlation coefficient = − 0.38, p < 0.001).Another gene, DOCK10, harbored a stop-gained (rs12328236), a missense (rs4674940) and an intron-near-splice (rs77952666) variant associated with cg12331557, which were considered novel associations compared to previously known meQTLs [16].It is worth noting that rs80176144 (~ 20 kb away from this locus) was associated with coronary artery calcified atherosclerotic plaque in T2D at p = 7 × 10 -6 [33] and is ~ 50 kb away from the identified meQTL SNP.The third gene, SLC2A1, harbored a 5'UTR SNP (rs11537640) that associated with cg19693031 in TXNIP, two genes that are biologically linked [34,35].Located * Epigenome-wide (EWAS) significant.⊥ indicates a previously known T2D methylated gene.Beta column indicates the effect size from regression analysis.Cohort 1 and cohort 2 were used in a 2-way discovery and replication procedure.The 66 identified CpGs were a combination of CpGs identified in cohort 1 and replicated in cohort 2 (asterisk in the cohort 1 column), and CpGs identified in cohort 2 and replicated in cohort 1 (asterisk in the cohort 2 column), using models either adjusted or unadjusted for BMI.The presented statistics were obtained from the statistical model indicated (including or excluding BMI as a covariate) when both models were significant.More details are in Additional file 2: Table S1 6738 bp upstream of rs11537640, rs12407920 had been associated with diabetic nephropathy (DN) [36], whereas cg19693031, besides its association with T2D, has been associated with kidney function [37].

meQTLs with GWAS associations and causal relationships
To investigate the biological functions of meQTL genes and identify whether meQTL SNPs are in loci that harbor GWAS associations with T2D relevant traits, we searched the GWAS catalogue [29] and the Stanford Biobank [38] for phenotypic associations with the meQTL genes.From the GWAS catalogue, 64 unique SNPs were identified at < 1kb away from meQTL SNPs, and were associated with several traits, some of which were relevant to T2D.Among those, 29 SNPs were exact matches of meQTL SNPs associated with 10 T2D CpGs (Fig. 5, Table 4, Additional file 2: Table S7).Among the most significant GWAS associations were meQTL SNPs in HK1 with HbA1c.Using the Stanford portal, we identified 20 unique SNPs located at < 1kb away from meQTL SNPs, that were associated with several traits (Table 4, Additional file 2: Table S8).Among those, 7 were exact matches of meQTL SNPs associated with 3 T2D CpGs.Both HK1 and PFKFB2 GWAS associations with HbA1c were significant in this database.For meQTL SNPs in HK1 and PFKFB2 that had associations with HbA1c, causal relationships of methylated genes with HbA1c as a biomarker of T2D were studied.The associations of methylation sites in those genes with HbA1c were first confirmed to be significant in our cohort.cg08992189 in HK1 was associated with HbA1c at a discovery p-value of p = 2.39 × 10 -5 and replicated at a p-value of p = 1.6 × 10 -4 , whereas cg22699725 in PFKFB2 was associated with HbA1c at a discovery p-value of p = 2.68 × 10 -6 and replicated at a p-value of p = 2.53 × 10 -5 .Afterwards, a two-Sample Mendelian Randomization (2SMR) analysis was conducted using CpG-SNP statistics from our cohort (meQTL results) and SNP-HbA1c statistics from the GWAS catalogue and Stanford Biobank.Table 5 shows the 2SMR results for PFKFB2 and HK1 with HbA1c (See Methods).S9).All 66 T2D CpGs identified in this study were tested for associations with the 75 T2D metabolites in 708 individuals.Those were divided into two cohorts of 364 and 344 individuals which were used for metaanalysis of the CpG-metabolite associations obtained from each.We identified 77 significant CpG-metabolite associations (p < 0.05/66*75), out of which 11 associations were removed due to a high heterogeneity, resulting in 66 significant associations in 24 unique CpGs (19 genes) and 25 unique metabolites (Table 6, Fig. 6, Additional file 2: Table S10).Of these, 61 were considered novel associations.

Multi-omics networks in T2D
A multi-omics network was constructed from all associations obtained and metabolic networks to identify pathways combining methylation and metabolism.Integrating the 66 CpG-metabolite pairs with 9 meQTL SNPs of the CpGs in those pairs and the T2D metabolites that had significant partial correlations with the metabolites involved in those 66 CpG-metabolite pairs, we obtained a network with 105 nodes (SNP/CpG/metabolite) and 165 edges (associations or correlations).The largest subnetwork had 75 nodes and 120 edges and harbored all CpGs that had metabolic associations (Fig. 7, Additional file 3: Figure S1).
To enrich the metabolic pathway information, the associations of metabolites with BMI, lipoprotein cholesterols and triglycerides were identified.Of the 43 T2D metabolites in this network, 16 showed higher significance in association with BMI, LDL-C, HDL-C or triglycerides compared to their association with T2D (Additional file 2: Table S11), among which 1-palmitoyl-2-oleoyl-GPE (16:0/18:1) was associated with triglycerides and linked to 13 CpGs.

Discussion
This study used a large number of subjects from the Qatari population that has a high prevalence of diabetes and used the largest methylation EPIC array.Whole genome sequencing and metabolomics data were incorporated to enrich the identification of biological pathways associated with methylation and diabetes.
A total of 66 CpG sites were significantly associated with T2D, of which 63 CpGs were novel in the insulin signaling pathway, fructose and mannose pathway and other T2D relevant pathways.Twenty-two of the  The most significantly-associated T2D methylated gene, thioredoxin-interacting protein (TXNIP), plays a role in insulin sensitivity, is a tumor suppressor that is upregulated in diabetes, is induced when glucose levels are elevated, and its deficiency improves glucose tolerance and increases insulin sensitivity in high fat diet-induced obesity [34,39].TXNIP suppresses glucose uptake by directly binding and suppressing the glucose transporter, GLUT1 (SLC2A1), by facilitating its clathrin-mediated endocytosis [34,40].When elevated, TXNIP induces β-cell apoptosis, while its deficiency protects against type I and type II diabetes by promoting β-cell survival [41,42].It is also known to be involved in kidney injury and with IL-1 β in the pathogenesis of T2D [43,44].Thus, findings from the meQTL analysis, linking SLC2A1 to TXNIP and its several metabolic associations support and add to those known functions of the gene.Besides TXNIP, pathway analysis and the known biology of the novel methylated genes support their relevance to T2D.For example, PPP1R3E, PRKAR2A and HK1 are involved in insulin signaling and resistance (Additional file 2: Table S5).Together with the several T2D-related pathways reported for HK1, further evidence of its relevance to T2D has been provided by GWAS associations of HK1 genetic variants with HbA1c [45], with the strongest signals reported in [46].HK1 is in the family of hexokinases, which are known to phosphorylate glucose to produce glucose-6-phosphate, the first step in most glucose metabolism pathways [47].Its association with the insulin signaling pathway and hyperinsulinemia [48], along with being part of the carbohydrate absorption, glycolysis, fructose and mannose metabolism pathways further emphasizes the statistically inferred causality relation.Other methylated genes have been associated with T2D incidence [49] and cardiovascular disease [50].For example, CACNA2D2 is involved in pathways related to T2D-associated coronary heart disease, cardiomyopathy, cardiac muscle contraction [51].Variants of DOCK6 have been associated with GWAS studies on coronary heart disease and total and HDL cholesterol levels [52,53].THBS4 has been associated with more advanced peripheral artery disease and T2D [54].It has also been reported to mediate breast cancer inflammation and growth in mouse models in response to hyperglycemia and TGF-beta [55].IDO2, one of the isoforms of IDO (Indoleamine 2, 3-dioxygenase) which forms part of the tryptophan metabolism and catalyzes the conversion of tryptophan to kynurenine, is upregulated in T2D patients with baseline tryptophan being associated with higher risk of incident T2D [56].Variants of PBRM1 have been reported to be associated with BMI-related traits and diabetes pathways [57,58].While the pathways and biological function of the identified genes could partially be explained from previous findings, our study reveals other biological functions through linking methylation with genomic variants and metabolic pathways.
The meQTL analysis indicated the effect of genomic variation on methylation levels of 22 T2D CpGs, whereas the remaining 44 T2D CpGs are more likely affected by Table 5 Mendelian randomization results for associations of HK1 and PFKB2 with HbA1c * x is considered the exposure (CpG) obtained from meQTL statistics in our cohort and y is considered the outcome obtained from GWAS databases (GWAS catalogue or Stanford portal).Bx and Bxse stand for effect size and standard error for x, and similar symbols are used for y † Results for the 3 SNPs are from the inverse variance weighted method rather than the maximum likelihood method used for single SNPs.Heterogeneity score (I2) is 90.34007, p = 2.4E-20  environment/lifestyle.The association of the missense SNP rs1136287 in the SERPINF1/SMYD4 locus with cg11692409 aligns with the previously reported GWAS association of SMYD4 with T2D that is harbored in the SRR locus [4].The DOCK10 association with diabetic atherosclerosis [33] and its gene expression association with glucose in pancreatic islets [59] suggest that the nonsynonymous meQTL SNPs in DOCK10 may link the gene's methylation to T2D genetic factors, diabetes complications or mechanisms in islet cells.The presence of a 5'UTR variant in SLC2A1 that codes the GLUT1 protein and that controls methylation in TXNIP is another novel finding that confirms previous biological links found between TXNIP and SLC2A1 [34,35,40] as discussed above.GWAS associations of meQTL SNPs with T2D-related traits such as HbA1c (HK1 and PFKFB2), BMI (EFEMP2, CFL1), eGFR and creatinine (PFKFB2), diabetic nephropathy, urate levels (EFEMP2), bilirubin (HK1), and white blood cell counts (DOCK10, EFEMP2, SLC2A1, CFL1, PFKFB2) [51,60,61], further emphasize the involvement of the meQTL genes in the pathogenesis of T2D.In other words, SNP variation may alter methylation levels in T2D CpGs inducing various biological perturbations and complications in T2D patients.Furthermore, the mendelian randomization analysis suggested a causal relationship of DNA methylation in HK1 with HbA1c levels, supporting previous findings of HK1's GWAS associations with HbA1c and HK1 pathways involving glucose [47], insulin signaling and hyperinsulinemia [48].

Gene
This study reports novel T2D methylation-metabolite associations.Compared to a few previously reported metabolic associations with TXNIP CpG cg19693031, this study reveals a much larger role of TXNIP in T2D metabolic pathways that span several pathways of interest, including branched chain amino acids, alanine, lipid metabolism, and sugars among others.DQX1, previously identified to be methylated with T2D in an Arab population, has no meQTL associated with it, suggesting that the metabolic associations of DQX1 with lipid metabolism and ribitol could be largely affected by lifestyle factors.Metabolic associations with genes known to associate with bladder cancer (BLCAP), ovarian cancer (OCIAD1), and infertility caused by azoospermia (DAZAP1) could a indicates associations previously reported in other studies.*indicates compounds that have not been officially confirmed based on a standard, but Metabolon is confident in its identity highlight the importance of integrating both methylation and metabolomics for assessing T2D complications.
Integrating methylation with both genomics and metabolomics has an important role in understanding the cause-effect direction of methylation-metabolism in T2D, i.e., whether methylation drives metabolic changes through gene expression or changes in metabolites cause changes in methylation.For example, methylation of TXNIP is controlled by SNPs in SLC2A1, known to be involved in glucose metabolism, and thus TXNIP association with the glycolysis pathway could be through this meQTL.Similarly, we showed that the meQTL of PFKFB2 had SNPs associated with HbA1c, eGFR, and creatinine, and at the same time the gene had associations with carbohydrates, lipids, phenylalanine and alanine that are involved in both T2D [17, 18] and kidney function [62,63], suggesting that the metabolic perturbations could be an end effect of the meQTL.Linking methylated genes to metabolism may help future investigation of the population-specificity of some methylated genes to Qataris.Perturbation of the methylation-metabolic pathway by changing lifestyle factors or nutrition in Qataris may result in different methylation patterns in Qataris compared to the other populations and thus may help explain T2D mechanisms in Qataris.
Multi-omics networks have facilitated the identification of metabolic pathways linked to the methylated genes, thus revealing their functionality through metabolites that were associated with T2D risks/complications.One example is the alanine subnetwork that links PFKFB2 and TXNIP to the urea cycle and to amino acids, where amino acids were associated with triglycerides (alanine association with triglycerides: discovery p-value of p = 6.4 × 10 -12 , replication p-value of p = 9.3 × 10 -13 ) and where both urea cycle and alanine have been associated with kidney failure [64].Moreover, 1-carboxyethylphenylalanine in the alanine subnetwork was associated with triglycerides (p-value discovery: 1.42 × 10 -34 , p-value replication: 4.6 × 10 -21 ) and is known to be associated with hypertension [65].Both links extend the GWAS findings for PFKFB2 described above and reveal the involvement of triglycerides in T2D methylation.The association of 1-palmitoyl-2-oleoyl-GPE(16:0/18:1) with triglycerides (discovery p-value p = 6.47 × 10 -100 , replication p-value p = 1.39 × 10 -80 ) and its link to 11 methylated genes suggests a possible cause of methylation linked to a diabetes risk/complication that is common between those genes.Mannonate, a xenobiotic (food component/plant) is another example that shares a network with six genes, namely TXNIP, BLCAP, THBS4, PEF1, DAZAP1, and KIAA1211L and may suggest a possible xenobiotic effect on methylation.The networks also highlighted a STYXL1-POR-sphingosines/ceramides-LDL,HDL-steroids connection, where the unique association of POR with a lipid subnetwork of sphingosines, GPCs, and ceramides is supported by the fact that gene's production of the enzyme cytochrome P450 oxidoreductase is required for the synthesis of cholesterol and steroid hormones (https:// medli neplus.gov/ genet ics/ gene/ por/, and references therein).Moreover, since reproductive and fertility complications may arise from diabetes [66], it is interesting to note that STYXL1, the meQTL gene of the cg01676795 (POR) is known to be associated with seminal vesicle tumor and male reproductive organ benign neoplasm [67], that may help in understanding the link of POR to a Fig. 7 Multi-omics network of T2D methylated genes associations with metabolites, and their meQTLs integrated with metabolic networks.Oval nodes are CpGs, diamonds are SNPs, and hexagons are metabolites network of sphingosines and ceramides which play a role in forming steroids [68].
Limitations of the study include measuring methylation from whole blood rather than pancreatic islets as that was the only available tissue from the Qatar Bio Bank.Another limitation was the inability to correct for medication as that information was not yet annotated by the provider at the time of the study.

Conclusions
A total of 66 CpG sites were identified to be associated with T2D, of which 63 were novel and linked to biological pathways of T2D.Genomic drivers of the CpG sites were identified in 688 significant CpG-SNP pairs, among which 181 CpG-SNP pairs were novel associations or in novel loci.Several GWAS associations in meQTLs' genes revealed various underlying factors in pathways linked to T2D.The study identified novel nonsynonymous and 5'UTR variants associated with T2D methylation in SERPINF1, DOCK10, and TXNIP, as well as the statistically-inferred causal relation between HbA1c and each of the HK1 and PFKFB2 methylation sites.A total of 61 novel methylation-metabolite associations revealed association of several methylated genes with T2D metabolic pathways including population specific genes such as DQX1, among others.The multi-omics network suggested a number of T2D biological mechanisms associated with methylation, including the association of triglyceride metabolism and xenobiotics to 11 and 6 methylated genes respectively, indicating a common factor among those methylated genes and may have a role in their methylation.The pathways identified linked to steroid metabolism, and to metabolites associated with BMI, lipoprotein cholesterols and kidney function, namely STYXL1-POR, SLC2A1-TXNIP, and PFKFB2, may reveal the association of T2D risk factors or complications with methylation mechanisms.Finally, this study revealed several novel methylated genes related to T2D, with their associated genomic variants and metabolic pathways.

Fig. 2
Fig. 2 Manhattan plot of 66 CpGs associated with T2D.The red line indicates the Bonferroni p-value threshold for significance

Fig. 5
Fig. 5 Bonferroni significant associations for associated with SNPs within 1kb of the meQTL SNPs.Traits are clustered according to the phenotype in > 9 clusters, indicated on the horizontal axis, that could be relevant to T2D and the blue bars are irrelevant to T2D.The marked clusters (*) have at least one exact matching meQTL and GWAS SNPs.47 associations are exact matches, with further details shown in Additional file 2: TableS7

Table 1
Sample characteristics for 3 cohorts selected from the Qatar Biobank

Table 2 66
CpG sites associated with T2D from two regression models (with and without BMI adjustment) replicated in the Qatari replication cohort or confirmed in previous studies with an exact pair match or in the same locus.A total of 130 pairs were considered novel CpG-SNP associations located in previously published loci and replicated in Qataris (Additional file 2: Table

Table 3
Fifty-one novel SNP-CpG pairs in novel loci a : chr and position of the CpG.*Beta column indicates the effect size from regression analysis 1-carboxyethylphenylalanine, alanine, erythronate and phosphatidylethanolamines, and POR (P450 Cytochrome Oxidoreductase) with 1,5 anhydroglucitol and a hexosylceramide among others.

Table 4
Significant GWAS associations with HbA1c and exact meQTL SNP matches a Repeated associations with the same SNP and traits with less significance are not shown

Table 6
66 significant CpG-metabolite associations sorted by gene and significance