Epigenome wide association study in peripheral blood of pregnant women identifies potential metabolic pathways related to gestational diabetes

ABSTRACT

Gestational diabetes mellitus (GDM) increases the risk of developing metabolic disorders in both pregnant women and their offspring. Factors such as nutrition or the intrauterine environment may play an important role, through epigenetic mechanisms, in the development of GDM. The aim of this work is to identify epigenetic marks involved in the mechanisms or pathways related to gestational diabetes. A total of 32 pregnant women were selected, 16 of them with GDM and 16 non-GDM. DNA methylation pattern was obtained from Illumina Methylation Epic BeadChip, from peripheral blood samples at the diagnostic visit (26–28 weeks). Differential methylated positions (DMPs) were extracted using ChAMP and limma package in R 2.9.10, with a threshold of FDR <0.05, deltabeta >|5|% and B >0. A total of 1.141 DMPs were found, and 714 were annotated in genes. A functional analysis was performed, and we found 23 genes significantly related to carbohydrate metabolism. Finally, a total of 27 DMPs were correlated with biochemical variables such as glucose levels at different points of oral glucose tolerance test, fasting glucose, cholesterol, HOMAIR and HbA1c, at different visits during pregnancy and postpartum. Our results show that there is a differentiated methylation pattern between GDM and non-GDM. Furthermore, the genes annotated to the DMPs could be implicated in the development of GDM as well as in alterations in related metabolic variables.

KEYWORDS:

• Gestational diabetes mellitus
• Epigenome-wide association
• epigenetics
• diabetes
• metabolism
• glucose homoeostasis

Background

GDM is defined as a carbohydrate intolerance with onset during pregnancy [1]. GDM can increase the risk of prematurity, C-section delivery, preeclampsia, the long-term risk of type 2 diabetes mellitus (T2DM), metabolic syndrome, renal disease and cardiovascular disease (CVD). Furthermore, offspring from women with GDM have an increased risk of foetal macrosomia, higher infant fat mass and long-term obesity [2]. Globally, the prevalence of GDM varies depending on the country and the diagnostic criteria, ranging from 1% to 31% [3]. In Europe, the prevalence ranges from 1.5% in Sweden to 10% in Italy [4,5]. Obesity, increasing reproductive age of women, sedentary lifestyle, and environmental factors are some of the principal risk factors for developing GDM [6].

During normal pregnancy, beta-cells secrete more insulin to compensate and maintain normal glucose levels, due to the insulin resistance, especially at the end of the gestation. However, in GDM, certain alterations exist that produce reduced insulin sensitivity and impaired insulin secretion [7], developing carbohydrate intolerance. The molecular mechanisms that produce these alterations are still unknown. GDM development is influenced by genetic, epigenetic, and environmental factors [8–11]. Epigenetics is the link between genetics and the environment, and it could explain how environmental factors affect the development of diseases. There are different epigenetic mechanisms involved such as histone modification, DNA methylation and some authors have also included non-coding RNAs. DNA methylation is the most extensively studied, and it consists of the addition of a methyl group to the five positions of cytosine residues [12]. There are different techniques to study DNA methylation, but one of the best ways to have a global view of methylation is with an epigenome-wide association study (EWAS), performing methylation arrays.

A few EWAS studies on GDM have been reported with different approaches such as the epigenetic changes in cord blood and the effect on the offspring and epigenetic patterns on the foetus-placental endothelium [13–15]. However, there are not many epigenetic studies on blood samples of mothers with GDM. According to an epigenetic study in peripheral blood in GDM vs non-GDM pregnancies, it has been revealed that five of the most differentiated methylated CpG in the following genes, Solute Carrier Family 9 Member A3 (SLC9A3), Male-Enhanced Antigen 1 (MEA1), Kelch Domain Containing 3 (KLHDC3) and Calmodulin Binding Transcription Activator 1 (CAMTA1) were correlated to glucose levels and fasting insulin in the first visit (at 23 weeks of pregnancy) [10]. Another study on peripheral blood of GDM vs. non-GDM presented five CpGs within COP9 signalosome Subunit 8 (COPS8), Phosphoinositide-3-Kinase Regulatory (PIK3R), 3-Hydroxyanthranilate 3,4- Dioxygenase (HAAO), Coiled-Coil Domain Containing 124 (CCDC124) and Chromosome 5 Open Reading Frame 34 (C5orf34) with biomarker potential for early therapeutic intervention [16]. Also, a recent study in peripheral blood identified three CpGs in Long Intergenic Non-Protein Coding RNA 917 (LINC00917), Trafficking Protein Particle Complex Subunit 9 (TRAPPC9) and Lymphoid Enhancer Binding Factor 1 (LEF1) genes, respectively, that may be implicated in the development of GDM and postpartum abnormal glucose tolerance [17].

In this study, we have compared by EWAS the epigenome of GDM vs non-GDM women in the moment of the diagnosis and studied the association with clinical variables related to diabetes in an attempt to elucidate the molecular mechanisms of the development of gestational diabetes.

Methodology

Subjects

The participants of this study are part of the EPI_DG study, which started at the beginning of 2019. The recruitment was done in the Diabetes and Pregnancy Unit at the University Hospital Virgen de la Victoria, Málaga, Spain. The inclusion criteria were age 18–45 years, singleton pregnancy and being in the second or third trimester of pregnancy. The existence of a GDM diagnosis at <14 weeks of pregnancy and multiple pregnancy were exclusion criteria. All the participants were of Caucasian Mediterranean ethnicity.

Sixteen GDM pregnant women and 16 non-GDM pregnant women (controls) were studied. The GDM diagnosis was made according to the National Diabetes Data Group NDDG criteria in two-step screening [18]: firstly, an O’Sullivan test (50 g glucose overload) was performed in primary care centres. If post-load glucose was ≥7.8 mmol/L, women were referred to our Diabetes and Pregnancy Unit for performing the 100 gr oral glucose tolerance test (OGTT-100 g). GDM was diagnosed if glucose values were higher than the threshold, at least in two points: fasting >5.8 mmol/l (OGTT0.0), after 1 h >10.6 mmol/l (OGTT60.0), after 2 h >9.2 mmol/L (OGTT120.0), and after 3 h >8.0 mmol/l (OGTT180.0). Pregnant women with normal OGTT-100 were considered as controls.

The study protocol includes a baseline visit (at the moment of performance of 100 gr OGTT, at 24–29 weeks of pregnancy; T0), a pre-labour visit (at 36–37 weeks of pregnancy; T1), and a pos-labour visit (8 weeks after birth; T2), all of them at our Diabetes and Pregnancy Unit. Height was measured at T0, and weight and blood pressure (BP) were measured at every visit, according to standardized procedures. At T0, we collected data regarding pre-pregnancy weight to calculate ‘pre-pregnant body mass index (BMI)’ as weight (kg)/height² (m²). Weight gain was calculated as the difference between the weight at T1 and at T0.

At T2, we collected anthropometric data—weight, length, and head circumference— from the digital medical history of newborns (these parameters were measured by the neonatologist during the first evaluation of newborns according to local protocols and compiled at the medical history).

Samples

Peripheral blood samples were collected at three different times: T0 during the diagnostic visit (24–29 weeks); T1 pre-labour visit (37–38 weeks) and T2 post-birth visit (8 weeks after birth). The 75 gr OGTT at T2 (OGTT.2), measuring fasting glucose (OGTT0.2) and 2 h postload glucose (OGTT120.2), was only performed in the GDM group to test these women for prediabetes or diabetes according to ADA recommendation [19]. Blood and serum were stored at −80ºC until the analysis.

Biochemical analysis

Glucose was measured by the oxidase method (Bayer, Leverkusen, Germany). Insulin levels were measured by a radioimmunoassay method using BioSource International Inc. (Camarillo, CA, USA). We calculated the homoeostasis model assessment of insulin resistance index (HOMA-IR) as described by Matthews et al. [20]. Serum 25(OH)D levels were determined by the Enzyme-Linked ImmunoAssay (ELISA) kit (Immundiagnostik, Bensheim, Germany). Cholesterol, triglycerides, and HDL were measured using commercial enzymatic methods. The blood pressure was measured twice with a sphygmomanometer with an interval of 5 min between measurements and the average of the two measurements was used in the analyses.

DNA isolation and bisulphite conversion

DNA was isolated using QIAamp DNA Blood Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Quality and concentration of DNA were measured using Qubit 3.0 Fluorometer with Qubit dsDNA HS Assay Kit Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). A total of 500 ng of genomic DNA from T0 was bisulphite treated with Epitect Bisulphite Kit (Qiagen, Germany) for posterior DNA methylation analysis.

Epigenome-wide DNA methylation

DNA methylation was tested using the Infinium Mehylation EPIC BeadChip Kit. R statistical software (https://www.rstudio.com/) was used to extract the data, using ChAMP package version 2.9.10 [21]. The probes were deleted using the following threshold: a detection p-value above 0.01 in one or more samples, with a bead count less than 3 in at least 5% of the samples, probes non-cpg and probes with SNPs [22,23]. Intra-cell type normalization was done using the beta-mixture quantile normalization (BMIQ) method. Houseman correction was used to correct the differences in methylation resulting from differences in cellular heterogeneity [24]. Finally, inflation levels were calculated with RaMWas package in R studio.

Methylation data analysis

β-values and M-values were calculated to obtain the methylation levels. M-values are a logical transformation of β-value, and they were used to obtain the differentially methylated positions (DMPs) by performing an eBayes moderated t-statistic with limma package in R statistical software. β-value is the estimate of the methylation level using the ratio of the methylation probe intensity and the overall intensity. The DMPs were filtered using the False Discovery Ratio (FDR) <0.05, the difference of β level between groups (deltaβ) higher than 5% in absolute terms, and the log odds that the CpG site is deferentially methylated (B) higher than 0. The models were adjusted by variables previously associated with GDM such as age, pre-pregnant BMI, foetal sex, weight gain, gestational age, and HOMAIR [25].

Functional analysis

A functional enrichment analysis was performed looking for the potential mechanism that may be altered in GDM compared with non-GDM. Once we obtained the DMPs with annotated genes, a gene ontology (GO) and Kyoto encyclopaedia of genes and genomes (KEGG) were performed. GO terms and KEGG pathways were obtained using Missmethyl package in R software. In addition, a String net (https://string-db.org/) was performed to identify protein–protein interactions (PPI). The significant threshold was an interaction score >0.7 and FDR <0.05. Clustering was performed in Cytoscape (https://cytoscape.org/) with MCODE app. Finally, each cluster was analysed by DAVID (https://david.ncifcrf.gov/summary.jsp) to identify the potential function.

Statistical analysis

Due to the small sample size and the fact that some variables did not show a normal distribution, non-parametric tests were used to compare quantitative variables and for correlation analysis. R software (4.0.4) was used to study the differences in clinical and biochemical variables using the U-Mann Whitney test. A correlation analysis by spearman test was performed between methylation levels and biochemical variables. Genes with more than 1 DMP, and those genes related with carbohydrate metabolism function in the functional analysis were chosen for this analysis. The correlation was made with clinical variables at different times (T0, T1 and T2). The correlation matrix was obtained and represented with the package GGally in R studio (4.0.4), and the threshold was r > 0.5 and p < 0.05.

Results

Clinical characteristic of the subjects

The characteristics of the subjects are represented in Table 1. Both groups were matched by age, gestational age, and pre-pregnancy BMI. The control group showed higher weight gain during pregnancy compared with GDM. The rest of the variables were not statistically different between groups. In the GDM group, after diet and exercise recommendation, six pregnant women required additional pharmacological treatment (insulin).

Table 1. Clinical characteristic of the subject included in the study.

Comparison	Basal (24/28 weeks)			Pre-Labor (36–38 weeks)			Pos-labour (8 weeks)
	Non-GDM	GDM	P-value	Non-GDM	GDM	P-value	Non-GDM	GDM	P-value
Age	34.2 ± 4.5	33.8 ± 4.1	0.808
Gestational age (weeks)^†	29 ± 3	28 ± 2	0.579	36 ± 2	36 ± 1	0.282
O’sullivan (mg/dl)	160.7 ± 16.3	172.9 ± 22.2	0.085
Weight (Kg)	75.4 ± 11.3	76.0 ± 13.5	0.892	78.7 ± 11.9	76.5 ± 13.1	0.622	71.9 ± 12	67.3 ± 11.5	0.335
Pre-pregnant BMI^†	24.6 ± 6.9	23.8 ± 6.9	0.862
SBP (mm Hg)	104.2 ± 9.9	111.4 ± 15.4	0.131	106.9 ± 12.8	109.1 ± 12.8	0.643	104.1 ± 7.32	108.3 ± 8.9	0.179
DBP (mm Hg)	69.1 ± 8.2	70.1 ± 7.8	0.710	71.2 ± 6.5	73.6 ± 8.8	0.393	78.4 ± 6	75.9 ± 5.4	0.242
Cholesterol (mg/dl)	272 ± 48.1	258.5 ± 51.1	0.446	273.6 ± 49.1	254.7 ± 54	0.308	214.4 ± 15.3	193.5 ± 35.3	0.09
HDL (mg/dl)	82.7 ± 18.2	77.8 ± 15.7	0.424	74.6 ± 12.3	73.2 ± 16.8	0.785	66.3 ± 15.3	58.1 ± 15	0.156
Tg (mg/dl)^†	194 ± 68	183 ± 38	0.829	249 ± 112	242 ± 79	0.749	89 ± 38	61 ± 57	0.890
HbA1C (%)	5.1 ± 0.3	5.3 ± 0.4	0.151	5.3 ± 0.3	5.8 ± 1.6	0.045	5.4 ± 0.3	5.4 ± 0.41	0.756
HOMA-IR^†	1.5 ± 1.1	1.7 ± 1.4	0.147	1.5 ± 1.4	2.1 ± 3.2	0.883	1.0 ± 0.7	0.8 ± 1.3	0.964
Weight_NB (g)							3455 ± 430.7	3309 ± 354.2	0.312
Length_NB (cm)							50.4 ± 1.4	50.3 ± 1.6	0.853
HC_NB (cm)							34.2 ± 1.2	34.1 ± 1.6	0.916
Weight gain (kg)				3.3 ± 1.4	0.49 ± 2.6	0.001

Data are expressed as the mean ± standard deviation or median ± interquartile range (IQR) for variables with no normal distribution. P-values for continuous data were calculated using U-Mann Whitney test. O’sullivan: glucose level in O’sullivan test. BMI: body mass index. HDL: high density lipoprotein cholesterol. SBP: diastolic blood pressure. SBP: systolic blood pressure. Tg: triglycerides. NB: newborn. HC: head circumference. Weight gain: calculated as the difference between the weight at the visit and at basal visit. ^†Median ± IQR.

Differential methylation analysis

A total of 1.141 differentially methylated positions (DMPs) were found comparing GDM vs. non-GDM, with a threshold of FDR <0.05, deltaβ >|0.5|, and B ≥0. Slight level of inflation was found in our EWAS (λ = 1.136) (Supplementary fig. S1). The PCA analysis shown in Figure 1a represents the two populations. As shown in Figure 1b, cg13824270 in pre-mRNA processing factor 4B (PRPF4B) gene, cg20409752 in CAMTA1 and cg050373962 are the most hypomethylated, and cg04248279 in rabphilin 3A like (RPH3AL) gene, cg20758759 in non-annotated gene, cg02518222 in hydroxysteroid 17-beta dehydrogenase 11 (HSD17B11) gene are the most hypermethylated (Supplementary table S1).

Figure 1. a) Principal component analysis between GDM (blue, cluster 2) and non-GDM (Red, cluster 1) b) Volcano plot, green: CpGs most hypometylated, red: most hypermethylated, blue: non-significant. c) DMPs according with feature position in the gene d) DMPs according with island position (Cgi position).

According to the position of the CpGs respecting the gene (feature position), 36% are located in the body, 37% in the intergenic region (IGR), around 20% in the promoters area (transcription start site (TSS) 1500 1500, 5’UTR, and TSS200) and the rest are in the beginning of the exon and in 3’UTR (Figure 1c). Related to the cgi position, most of the DMPs were in the Open sea position (67%). The rest of the DMPs were in shore (19%), island (7%) and shelf (7%) positions, respectively (Figure 1d).

Functional enrichment analysis

In order to understand what potential mechanisms are involved in GDM, a functional enrichment analysis was performed. First, a gene ontology (GO) and a Kyoto encyclopaedia of genes and genomes (KEGG) analyses were performed in all the CpGs with genes annotated (714). The most significant GO terms are represented in Figure 2a. Among the significant GO, we identified the following mechanisms: the regulation of insulin secretion, cellular response to amino acid stimulus and response to nutrient levels. KEGG analysis showed that the most significant pathways were the insulin signalling pathway, growth hormone synthesis, ErbB signalling pathway, and T2DM (Figure 2b).

Figure 2. a) GO terms most enrichment pval<0.01 b) KEGG pathways most enriched. c)String net PPI (p-value<1.0e-16) with 23 genes related to glucose metabolism.

In addition, a String net was performed in genes annotated (714) to identify protein–protein interactions (PPI) and their related mechanisms. The resulting net was clustered in Cytoscape using MCODE app and analysed in DAVID. As shown in Figure 2c, a total of 15 clusters were defined using the default values of MCODE app. The three most significant (enrichment p-value <1e-16) clusters were related to carbohydrate metabolism pathways: positive regulation of glycolytic process, MAPK signalling pathway, insulin pathway, cellular response to nutrient levels, type 1 diabetes, T2DM, and insulin resistance. A total of 23 of our annotated genes were included in these clusters, and there was also interaction between them.

Correlation analysis

Within the top 100 DMPs, genes with more than one DMP and those genes related with carbohydrate metabolism according to our enrichment analysis were selected to perform a correlation analysis with the biochemical variables at different times T0, T1, and T2. A total of 27 CpGs were significantly correlated to biochemical variables (p < 0.05, r>|0.5|) (Table 2).

Table 2. CpGs most correlated to biochemical variables.

CpG	Genes	Variables	r	P value
cg07894331	MAPK1	BMI_p	0.5187	0.002
cg12080079	PAX7	Chl.1	0.5809	0.0005
cg09418268	MTOR	Glucose.0	0.5391	0.001
cg12854232	ERBB4	Glucose.0	−0.5263	0.002
cg20344956	LOC646214	HbA1c.2	0.53	0.003
cg20409752	CAMTA1	HOMAIR.0	−0.5521	0.001
cg07420274	GAD1	HOMAIR.0	−0.5227	0.002
cg14688342	EGFR	OGTT0.2	−0.6216	0.023
cg04217177	ESRRB	OGTT0.2	0.6725	0.01
cg17434154	ESRRG	OGTT0.2	0.5535	0.04
cg21873971	PDIA5	OGTT0.2	0.7919	0.001
cg05894391	PRKAA2	OGTT0.2	0.6177	0.02
cg11055991	PTPRN2	OGTT0.2	0.5803	0.037
cg14479198	RASGRF2	OGTT0.2	−0.5752	0.039
cg19891476	SYK	OGTT0.2	0.5576	0.047
cg13487126	AKT2	OGTT120.2	0.5939	0.032
cg10211745	MINOS1-NBL1	OGTT120.2	0.6742	0.011
cg19234140	MINOS1-NBL1	OGTT120.2	0.6898	0.009
cg14579430	NBL1	OGTT120.2	0.6161	0.024
cg15589641	NBL1	OGTT120.2	0.6261	0.022
cg12474394	NBL1	OGTT120.2	0.6912	0.008
cg14715327	OR4A16	OGTT120.2	0.6704	0.012
cg20230784	GNAQ	OGTT180.0	0.5602	0.0008
cg17816038	CAMTA1	Tg.1	−0.5287	0.002
cg00856432	CACNA1D	VitaminD.2	0.5914	0.0007

r: Spearman correlation. OGTT0: Oral Glucose Tolerance Test first point. OGTT120: Oral Glucose Tolerance Test, point 120 min. OGTT180: Oral Glucose Tolerance Test, point 180 min. Tg: Triglycerides. BMI_p: pre-pregnancy body mass index. Chl: cholesterol levels. .0: basal visit. .1: pre labour visit. .2: post labour visit..

The highest correlations were for OGTT0.2 in epidermal growth factor receptor (EGFR) (cg14688342), protein disulphite isomerase family A member 5 (PDIA5) (cg21873971), protein kinase AMP-activated catalytic subunit alpha 2 (PRKAA2) (cg05894391) genes, and in non-annotated region cg04217177; and for OGTT120.2 in DAN Family BMP antagonist (NBL1) gene (cg14579430, cg15589641, cg12474394, and cg14715327) and in the mitochondrial inner membrane organizing system 1-NBL1 (MINOS1-NBL1) gene (cg10211745, cg19234140) (Figure 3). There is a strong correlation between CpGs in the same gene (r > 0.9), and in the same direction (Figure 3b).

Figure 3. Correlation-plots a)the most correlated CpGs with OGTT0.2 (r>|0.6|). b) the most correlated CpGs with OGTT120.2(r>|0.6|). *: p-value<0.05, **: p-value<0.01, ***: p-value<0.005, ****: p-value<0.001. OGTT0: Fasting glucose; OGTT120: Oral Glucose Tolerance Test, point 120 min. .2: post labour visit.

Discussion

In this work, we found a total of 1041 DMPs in the peripheral blood of GDM vs. non-GDM women. Based on functional analysis of the 714 genes annotated, we found 23 genes related to positive regulation of glycolytic process, MAPK signalling pathway, insulin pathway, cellular response to nutrient levels, T2DM, and insulin resistance. Furthermore, the correlation analysis showed that some of the most enriched genes are correlated to glucose and insulin variables (fasting glucose, different points of OGTT, HbA1C, HOMAIR) at baseline and post-labour clinical analytics.

Currently, there is a growing interest in the study of genetic and epigenetic patterns in GDM pregnant women and their offspring [14,26]. Our results agree with recent studies. Ballesteros et al. [17] identified 3 DMPs within LINC00917, TRAPPC9, and LEF1 genes in GDM women with abnormal glucose tolerance, and also LINC00917 and TRAPPC9 were associated with abnormal glucose levels after 4 years of postpartum. In another study, carried out in maternal blood before the GDM diagnosis, they found five DMPs in COPS8, PIK3R5, HAAO, CCDC124, and C5orf34 genes with the potential to be clinical biomarkers of GDM development [16]. In another study, in a group of South African women with GDM and non-GDM, the five CpG more differentiated were in SLC9A3, MEA1, KLHDC3, CAMTA1, and RASA3 genes. According to their functional analysis, these genes were related to signal transduction, cell growth, proliferation, differentiation and apoptosis, insulin resistance, glucose metabolism, inflammation, neurological signalling, and oncogenesis [10].

CAMTA1 was one of the most hypomethylated genes and is correlated with relevant metabolic variables. CAMTA1 seems to be implicated in the beta cell insulin regulation and secretion [27]. Guerra-Mollet suggests that in response to a high level of sugar there is a decreased level of CAMTA1 downregulating the expression of INS1 and INS2, and decreasing insulin secretion. Recent studies highlight the importance of CAMTA1 at multiple levels, from the detection of glucose to insulin production and secretion, and beta cell function. As it is well known, insulin regulates the glucose and lipid metabolism [28]. The hypomethylation observed in the GDM group could imply an altered expression of CAMTA1 affecting the insulin secretion, and in consequence in the regulation of glucose and lipid metabolism [28], as we observed with a negative correlation with HOMA-IR and triglyceride levels. However, this effect should be evaluated with additional experiments. Dias et al. found, as we did, a lower methylation of CAMTA1 in GDM vs non-GDM [29]. Although our DMPs are different from their work, the location of both our CPGs and theirs are in the body region, meaning that the biological response may be similar. Finally, different studies have shown an association between polymorphisms in CAMTA1 and carbohydrate metabolism such as blood lipid traits and blood pressure [30], fasting glucose [31] and T2DM [32].

It is known that insulin secretory defect, high insulin resistance, and beta cell dysfunction may be implicated in the development of GDM [7,33,34]. Accordingly, in our functional enrichment analysis of the DMPs, we found that the most significant GO and KEGG were related to the regulation secretion of insulin, insulin signalling pathway, and T2DM. All of these pathways are implicated in the development of T2DM [35–37]. Therefore, the epigenetic pattern of these genes may be implicated in the development/risk of GDM.

We found that some CpGs were correlated with some clinical variables at different times. Most of these associations are related to glucose levels and the OGTT result after labour. Several genes seem to be associated with the response to the glucose overload. For example, AKT serine/threonine kinase 2 (AKT2) and mechanistic target of rapamycin kinase (MTOR) are known for their implication in the insulin signalling pathways. Some studies in mice suggest that downregulation of AKT2 is related to high levels of glucose and insulin resistance [38]. In humans, certain polymorphisms associated with decreased levels of AKT2 were related with impaired glucose, insulin resistance, and T2DM risk [39,40]. In general, hypermethylation in promoters regions is associated with downregulation, whereas gene body methylation is associated with expression in several studies, although this latter is still unclear [41]. We have found that AKT2 is hypermethylated in the promoter region (5’UTR), suggesting that this gene could be downregulated and therefore contributing to the increase in glucose levels. Regarding MTOR, overexpression is associated with anabolic process and insulin secretion [42,43]. In our study, CpGs of MTOR were hypermethylated in the body region, which could involve a higher expression of this gene. PRKAA2 was associated with glucose levels after OGTT, and this gene seems to be related to insulin sensitivity [44] and some polymorphisms have been associated with T2DM [45]. Furthermore, PRKAA2 is closely related to MTOR regulation, i.e., downregulation of PRKAA2 induces MTOR expression and an increase in protein synthesis and insulin resistance in rats [46]. Perhaps higher methylation in the TSS1500 region of PRKAA2, as we have found, could mean a downregulation of this gene resulting in a decrease in insulin secretion and higher glucose levels after OGTT.

In our result NBL1 and MINOs-NBL1 had more than 1 CpG correlated with the response to OGTT at 120 min, after labour. NBL1 could be implicated in the regulation of β cell function and in diabetic kidney disease acting as a bone morphogenetic protein (BMP) antagonist, which are part of the TGF-β pathway [47,48]. Therefore, lower level of NBL1 may induce a dysregulation in beta cell function. As a result, an altered function of β cell, maintained in the time, could result in an increased risk of high glucose levels in response to an OGTT. Also, there are recent studies that suggest the implication of the BMP family in obesity and glucose metabolism (potential functions of the BMP family in bone, obesity, and glucose metabolism). Therefore, the antagonist role of NBL1 through BMP could be implicated in the glucose level after an OGTT in the post-labour visit. However, more studies would be necessary to understand the role of NBL1 gene in glucose metabolism.

Our study presents some limitations that we must comment. Although the differences found are relevant and significant beyond the FDR threshold, the sample size is small, and the results should be interpreted with caution. Furthermore, validation in a larger cohort or with other technologies would reinforce our findings. Nevertheless, we consider that the probability of false positives is low due to the small genomic inflation lambda factor obtained, but more studies are necessary to confirm our results. Regarding the characteristics of our population, the control group are pregnant women with a previous positive O’Sullivan test, so this issue should be considered when extrapolating these results to other populations. Finally, with our study design, we cannot demonstrate that the associations found imply causality.

Regarding the strengths of our study, we highlight the homogeneous selection of the sample, emphasizing that the differences found are mostly due to the different epigenetic marks and not to other possible confounding variables, which were also considered in the study. The main difference between the groups was that the control group showed higher weight gain during pregnancy compared with GDM. This finding is not surprising and indicates good adherence from GDM pregnant women to therapeutic measures implemented at diagnosis. Furthermore, our results are in accordance with previous studies that found a different methylation pattern between GDM and non-GDM. In fact, our findings were strongly significant at a threshold of FDR <0.05.

Conclusion

In addition, these epigenetic differences correspond to genes involved in the glucose metabolism and are associated with some biochemical variables after labour. These results help to understand which mechanism could be modified in GDM and how the methylation profile at the diagnostic visit can be associated with glucose metabolism-related variables in the future.

Data availability of statement

The data sets used during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

This study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of the Hospital Universitario Virgen de la Victoria, Málaga.

Authors’ contributions

TML-P: methodology, analysis, interpretation of the data and wrote original draft, NP-M: processed the samples and performed the experiments prior to the DNA methylation, CG-R: participated in the design of the work, FL-M: sample collection, AS-P: revision of the manuscript, FJ-T: supervision and revision of the paper, MM-V: conceptualization and acquired data from the patients and collected samples, MJ-P: conceptualization, funding acquisition and contributed to the discussion, SM: conception/design, funding acquisition, interpretation of the data and has substantively revised the work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2023.2211369.

Additional information

Funding

This study was supported by the Juan Rodés program from “Instituto de Salud Carlos III” (JR20-00040 to MM-V), Miguel Servet Type I program from the ISCIII-Madrid, Spain (CP20/00066 to CG-R), PFIS program (FI19/00178 to TML-P), Garantía Juvenil program (POEJ-0039-18) from Ministerio de trabajo y economía social de España, and the Nicolas Monardes Program from the “Servicio Andaluz de Salud, Junta de Andalucía”, Spain (RC-0008-2021 to SM). In addition, this study was supported by the “Centros de Investigación Biomédica en Red” (CIBER) of the Institute of Health Carlos III (ISCIII) (CB06/03/0018) and research grants from the ISCIII (PI18/01175, PI21/01864) and from “Servicio Andaluz de Salud”, Junta de Andalucía (PI-0283-2018, PI-0419-2019). This study has been co-funded by FEDER funds (“A way to make Europe”).