DNA methylation profiles in sibling pairs discordant for intrauterine exposure to maternal gestational diabetes

ABSTRACT

Intrauterine exposure to hyperglycemia is reported to confer increased metabolic risk in later life, supporting the ‘developmental origins of health and disease’ hypothesis. Epigenetic alterations are suggested as one of the possible underlying mechanisms. In this study, we compared pairwise DNA methylation differences between siblings whose intrauterine exposure to maternal gestational diabetes (GDM) were discordant. Methylation of peripheral blood DNA of 18 sibling pairs was measured using Infinium HumanMethylation450 BeadChip assays. Of the 465,447 CpG sites analyzed, 12 showed differential methylation (false discovery rate <0.15), including markers within genes associated with monogenic diabetes (HNF4A) or obesity (RREB1). The overall methylation at HNF4A showed inverse correlations with mRNA expression levels, though non significant. In a gene set enrichment analysis, metabolism and signal transduction pathways were enriched. In conclusion, we found DNA methylation markers associated with intrauterine exposure to maternal GDM, including those within genes previously implicated in diabetes or obesity.

KEYWORDS:

• DNA methylation
• epigenome-wide association study
• gestational diabetes mellitus
• siblings
• intrauterine environment

Abbreviations

CTNND2	=	Catenin delta 2
GDM	=	Gestational diabetes mellitus
HM450K	=	Illumina Infinium HumanMethylation450 BeadChip assays
HNF4A	=	Hepatocyte nuclear factor 4 α
RREB1	=	Ras responsive element binding protein 1

Introduction

Hyperglycemia during pregnancy is associated with increased risk of obesity and type 2 diabetes (T2D) in the offspring.^1,2 It has been reported that offspring born to mothers with gestational diabetes mellitus (GDM) have 8-fold increased risk of T2D compared with general population at early adulthood.³ This could be due to shared genetic risk factors between GDM mothers and their offspring. Recent advances in genetic association studies have improved our understanding on genetics of T2D, with identification of at least 80 susceptibility loci.⁴ However, these common genetic variants explain only a limited proportion of the T2D heritability.⁵ This implies that environmental factors including early-life exposures may also contribute to the pathogenesis of T2D in later life, as suggested by the hypothesis of ‘developmental origins of health and disease’.⁶⁻⁸

Intrauterine hyperglycemia is suggested as one of the environmental factors that confer life-long increased risk of metabolic diseases. It has been reported that in Pima Indians, offspring conceived after the mother developed T2D had increased risk of obesity and T2D compared with the sibling born before their mother was diagnosed with T2D.⁹ The detailed mechanism by which intrauterine hyperglycemia affects later life adverse metabolic outcomes is not fully understood. One of the possible mechanisms is that intrauterine hyperglycemia might induce changes in DNA methylation.¹⁰ During this period, rapid cell division and dynamic DNA methylation ensue and it is vulnerable to adverse environmental factors. Epigenetic alterations in response to early-life exposures could affect a wide range of tissues and maintained through cellular memory during subsequent cell divisions.^11,12 In addition, recent studies showed that DNA methylation in peripheral blood leukocytes can represent DNA methylation of metabolically relevant tissues, such as liver, adipose tissue, and pancreas.^13,14

There have been several reports on differential DNA methylation associated with maternal hyperglycemia in placental tissues or cord blood samples of offspring.¹⁵⁻¹⁸ However, the sample size was small and the findings were not confirmatory. Other studies investigated DNA methylation in tissues relevant for diabetes, such as skeletal muscle and adipose tissue in adult offspring.^19,20 In these studies, transcription level as well as DNA methylation of PPARGC1A, which plays an important role in adipogenesis and metabolic regulation, and adipokines, such as LEP, ADIPOQ, and RET were investigated. Recently, a DNA methylation study using peripheral blood was performed in offspring exposed to T2D during pregnancy.²¹ Several CpGs in/near LHX3, PRDM16, AK3, PCDHGA4, and STC1 were strongly associated with exposure to maternal T2D. Most of the previous studies used a case-control study design with unrelated offspring exposed to either hyperglycemic or normoglycemic environments in utero. One of the limitations of the previous studies was the potential confounding by different genetic backgrounds. Investigating sibling pairs discordant for exposure to maternal GDM might enable us to elucidate epigenetic change that is specific to intrauterine hyperglycemia without confounding the effect of genetic variation.

In this study, we conducted a discordant sibship study with siblings born to women who had experienced both GDM and non-GDM pregnancies. Using 18 discordant sibling pairs who were in their childhood, we aimed to identify peripheral leukocyte DNA methylation changes associated with exposure to maternal GDM.

Results

Clinical characteristics of study participants

A total of 18 sibling pairs participated in this study. In each sibling pair, one was conceived during GDM pregnancy and the other during non-GDM pregnancy. Clinical characteristics of the study participants are displayed in Table 1. The mean ages were 5.8 ± 1.4 and 9.4 ± 2.5 y for offspring from GDM pregnancy and non-GDM pregnancy, respectively (P < 0.0001). Siblings exposed to maternal GDM were younger than the other sibling, except for one sibling pair. There was no significant difference in sex distribution between the 2 groups. Among the 18 sibling pairs, 10 were sex concordant and 8 were sex discordant. There was no significant pairwise difference in BMI Z-scores between the 2 groups (P = 0.594). The mean maternal age at GDM pregnancy was 34.4 years, while the mean age at non-GDM pregnancy was 30.8. All the mothers were never smokers. The glucose values measured during 100 g oral glucose tolerance test (OGTT) are displayed in Table 1. Of the 18 mothers, 6 received insulin treatment during GDM pregnancy, and none of the mothers used oral anti-diabetic medications.

Table 1. Characteristics of the study population.

Characteristics of offspring	Siblings exposed to GDM in utero (n = 18)	Siblings not exposed to GDM in utero (n = 18)	P value
Current age (years)	5.8 ± 1.4	9.4 ± 2.5	<0.0001
Sex (male/female)	11/7	10/8	0.999
Maternal age at delivery (years)	34.4 ± 2.5	30.8 ± 2.4	<0.0001
Delivery (NFSD/PSD/C-sec)	6/0/12	5/1/11	0.999
Birth weight (kg)	3.3 ± 0.4	3.2 ± 0.6	0.671
Current BMI Z-score	0.19 ± 0.73	0.31 ± 0.86	0.594
Current HbA1C (%)	5.3 ± 0.3	5.3 ± 0.2	0.326
Characteristics of mothers^*	During GDM pregnancy (n = 18)
Age (years)	34.4 ± 2.5
BMI before pregnancy (kg/m^2)	22.7 ± 2.4
BMI at delivery (kg/m^2)	28.6 ± 8.2
Cigarette smokers	None
Gestational weeks at 50 g OGTT (weeks)	26 ± 2
Glucose levels at 50 g OGTT (mg/dL)	170 ± 26
Gestational weeks at 100 g OGTT (weeks)	28 ± 2
Glucose levels at 100 g OGTT (mg/dL)
0-hour	97 ± 18
0.5-hour	206 ± 40
2-hour	178 ± 39
3-hour	154 ± 34
Insulin treatments	6 treated (33%)

* Maternal data during GDM pregnancy only were available.

Continuous data are presented as mean ± standard deviation and P values were calculated using paired t-test. Categorical data were presented in counts and tested with Pearson's Chi-square test. BMI Z-score were calculated according 2007 Korean National Growth Charts. NFSD, normal full-term spontaneous delivery; PSD, preterm spontaneous delivery; C-sec, Caesarean section.

Differentially methylated CpG sites

For differential methylation analyses, a total of 465,447 CpG sites, which were at detection P<0.01, outside of single nucleotide variations, in CpG contexts, and located at autosomal chromosomes, were included. We investigated pairwise DNA methylation differences using M-values for its association with exposure to maternal GDM pregnancy. The volcano plot of the CpG markers is shown in Fig. 1A. After adjusting for age, sex, and experimental batch (Model 1), one CpG site, cg04988367, which was located 39 kb downstream of catenin delta 2 (CTNND2) showed significant difference at q < 0.10. Additional 11 CpG sites had pairwise methylation difference that was significant at q < 0.15 (Table 2). Among these, cg08407434 (coefficient for log₂ fold difference 1.003, P = 7.17 × 10⁻⁷, q = 0.109) was located at intron of hepatocyte nuclear factor 4 α (HNF4A), known as the causative gene for maturity onset diabetes of the young (MODY) type 1 and was the most significant differentially methylated marker after additional adjustment for cell type composition (Model 2). Another CpG marker cg18255813 was located at intron of ras responsive element binding protein 1 (RREB1), which was reported to be associated with visceral fat mass ²² and glycemic traits.²³

Figure 1. Differential DNA methylation analysis in sibling pairs discordant for intrauterine exposure to maternal GDM. (A) Volcano plot for the differential DNA methylation analysis of the 465,446 CpG sites. The effect size of GDM exposure on DNA methylation levels at each of the CpG sites is shown in x axis and the significance levels are shown in y axis. Coefficients > 0 refer to increased methylation in the GDM group compared with the non-GDM group. (B) Heatmap of adjusted β values for the top 100 differentially methylated CpG sites associated with intrauterine exposure to maternal GDM. The 36 samples (column-wise) and 100 CpG sites (row-wise) associated with intrauterine exposure to GDM were displayed using unsupervised hierarchical clustering. The upper panel and the lower panel represent sibling pairs and group information (black: GDM group; gray: non-GDM group), respectively. Blue to red gradient correspond to increasing methylation. Manhattan distance measures were used to compute inter-sample distances based on DNA methylation levels. S1-S18 denotes sibling IDs and G and N represent the GDM group and the non-GDM group, respectively.

Table 2. Differential DNA methylation in siblings discordant for their exposure to maternal GDM (q < 0.15).

CpG ID	Position (hg19)	Adjusted differences,*^,† (M-values)	Adjusted differences,*^,† (β values)	Raw β differences	P value^†	q value	P value^‡	Nearby gene	CGI information	Gene property
cg04988367	chr5:10,933,006	1.098	0.039	0.026	2.45 × 10^−7	0.09	3.12 × 10^−4	CTNND2	Open sea	Intergenic
cg08407434	chr20:43,001,691	1.003	0.039	0.013	7.17 × 10^−7	0.109	1.51 × 10^−5	HNF4A	Open sea	Body
cg01617280	chr8:8,745,938	1.094	0.077	0.031	1.08 × 10^−6	0.109	1.74 × 10^−5	MFHAS1	Shelf	Body
cg13911801	chr9:35,860,912	0.981	0.078	0.042	1.46 × 10^−6	0.109	1.36 × 10^−4	LOC92973	Open sea	Body
cg09304293	chr6:163,475,252	1.186	0.066	0.041	1.56 × 10^−6	0.109	2.16 × 10^−4	PACRG	Open sea	Body
cg23109897	chr8:48,871,883	0.927	0.158	0.065	1.90 × 10^−6	0.109	5.79 × 10^−4	PRKDC	Shore	Body
cg06319822	chr16:215,960	−0.652	−0.016	−0.005	2.08 × 10^−6	0.109	2.30 × 10^−4	HBM	Island	TSS200
cg06928695	chr17:6,384,119	0.713	0.086	0.059	2.39 × 10^−6	0.11	1.49 × 10^−4	PITPNM3	Open sea	Body
cg18255813	chr6:7,195,966	1.131	0.029	0.007	3.10 × 10^−6	0.127	3.95 × 10^−4	RREB1	Open sea	Body
cg26952618	chr16:10,912,545	−0.785	−0.124	−0.063	4.34 × 10^−6	0.142	7.63 × 10^−4	FAM18A	Island	5′UTR
cg19565171	chr6:73,935,200	1.17	0.053	0.032	4.58 × 10^−6	0.142	5.16 × 10^−4	KHDC1L	Open sea	TSS200
cg10431340	chr1:161,279,108	1.127	0.063	0.030	4.65 × 10^−6	0.142	2.76 × 10^−4	MPZ	Open sea	Body

*log₂ fold difference in methylation levels (sibling exposed to GDM – control sibling). Positive values indicate higher DNA methylation levels in the sibling exposed to maternal GDM. ^†Adjusted for age, sex and experimental batch (Model 1) ^‡Adjusted for age, sex, experimental batch, and cell type composition (Model 2). CHR, chromosome; CGI, CpG island; hg19, human genome version 19; UTR, untranslated region; TSS200, within 200 base pair upstream of transcription start site.

Clustering analysis

We performed unsupervised clustering analysis using top 100 differentially methylated CpG sites after adjustment for covariates (age, sex, and experimental batch). One large cluster mainly consisted of GDM offspring (GDM:non-GDM = 12:3), while another cluster consisted of a larger number of non-GDM offspring (non-GDM:GDM = 15:6) (Fig. 1B). Non-GDM offspring that failed to be clustered into the non-GDM-dominant cluster were in the same cluster with their sibling counterparts and the same was true for the GDM offspring that failed to belong to the GDM-dominant cluster. GDM offspring in the non-GDM-dominant cluster were closely clustered with each other.

Gene set enrichment analysis

Gene set enrichment analysis was performed using the ranked list of differentially methylated genes and pre-defined pathway database gene sets from ConsensusPathDB.²⁴ Genes with differentially methylated CpG sites preferentially participated in ‘metabolism’ and ‘signal transduction pathways’. Genes involved in ‘regulation of insulin like growth factor transport and uptake by insulin like growth factor binding proteins’ (enrichment score 0.770, P = 0.0024) were most significantly enriched in the ‘signal transduction pathways’. In addition, genes involved in several pathways including ‘superpathway of conversion of glucose to acetylCoA and entry into the TCA cycle’ (enrichment score 0.711, P = 0.0379) were significantly enriched in ‘metabolism’ pathways (Table 3).

Table 3. Gene set enrichment analysis of genes associated with differentially methylated CpG sites.

Gene sets	Number of genes	Enrichment score	Normalized enrichment score	P value	q value
Signal transduction pathways
Regulation of insulin like growth factor transport and uptake  by insulin like growth factor binding proteins (Reactome)	15	0.770	1.538	0.0024	1
TRK receptor signaling mediated by the MAPK pathway (PID)	18	0.740	1.528	0.0049	1
β3 integrin cell surface interactions (PID)	20	0.714	1.483	0.0154	1
Retinoic acid receptors mediated signaling (PID)	15	0.694	1.417	0.0358	0.7638
IGF signaling (INOH)	22	0.673	1.433	0.0219	0.6932
Regulation of RAS family activation (PID)	21	0.665	1.405	0.0379	0.7746
FGF signaling pathway (PID)	23	0.663	1.408	0.0351	0.8159
GPCRs other (WikiPathways)	35	0.634	1.445	0.0199	0.7079
Human thyroid stimulating hormone signaling pathway  (WikiPathways)	39	0.632	1.459	0.011	0.8781
Signaling by FGFR (Reactome)	35	0.597	1.366	0.0382	0.7321
Fibroblast growth factor 1 (NetPath)	40	0.593	1.383	0.0427	0.754
MAPK signaling pathway (WikiPathways)	73	0.547	1.33	0.0427	0.7644
Signaling by GPCR (Reactome)	457	0.446	1.164	0.034	0.8706
Metabolism
Superpathway of conversion of glucose to acetylCoA and  entry into the TCA cycle (HumanCyc)	15	0.711	1.418	0.0379	0.7746
C21 steroid hormone biosynthesis and metabolism (EHMN)	25	0.664	1.449	0.0186	0.8017
Androgen and estrogen biosynthesis and metabolism  (EHMN)	39	0.663	1.54	0.0033	1
Central carbon metabolism in cancer (KEGG)	33	0.652	1.46	0.009	0.9084
Metabolism of nucleotides (Reactome)	29	0.646	1.443	0.0174	0.7035
Selenocysteine synthesis (Reactome)	41	0.623	1.448	0.0154	0.7547
Glycerophospholipid metabolism (KEGG)	36	0.608	1.372	0.0446	0.75
Pi metabolism (Reactome)	31	0.605	1.378	0.0382	0.7201
Biological oxidations (Reactome)	82	0.576	1.399	0.0093	0.7397
Purine metabolism (KEGG)	62	0.562	1.355	0.0295	0.7328
Phospholipid metabolism (Reactome)	83	0.551	1.359	0.0115	0.7141
Metabolism of carbohydrates (Reactome)	119	0.532	1.338	0.0192	0.7502
Metabolism of amino acids and derivatives (Reactome)	143	0.507	1.294	0.0193	0.8099

Gene set enrichment analysis was performed with differentially methylated genes and gene sets from ConsensusPathDB. Gene set names and their source database (shown in parenthesis) with nominal P value lower than 0.05 are listed. FDR q values were presented. PID, pathway interaction database; INOH, integrating network objects with hierarchies; EHMN, Edinburgh

Overall differential methylation of HNF4A and gene expression

As this study had limited statistical power, we focused on a candidate site that had biologic plausibility. HNF4A, well known as the causative gene for MODY1, plays a crucial role in pancreatic β cell development and differentiation. Therefore, the DNA methylation pattern of HNF4A was further investigated. First, we performed pyrosequencing of cg08407434 and found a significant correlation with the methylation values measured using Illumina Infinium HumanMethylation 450 BeadChip (HM450K; Pearson's coefficient 0.469, P = 0.01) (Supplementary Figure S2.). Second, we compared DNA methylation levels of cg08407434 within each of the sibling pairs using β values that were adjusted for age, sex, and experimental batch. We found that all sibling pairs had increased methylation in the offspring exposed to maternal GDM pregnancy (Fig. 2A). Third, the average DNA methylation levels of all CpG markers within the HNF4A region showed a marginal trend of increased DNA methylation in siblings exposed to GDM except for one pair (P = 0.052) (Fig. 2B). Finally, there was a negative correlation between overall HNF4A methylation and mRNA expression, albeit statistically insignificant (Pearson's coefficient -0.364, P = 0.068) (Fig. 2C).

Figure 2. Within sibling pair comparison of DNA methylation levels at cg08407434 and HNF4A gene and correlation with mRNA expression. (A) Within sibling pair differences in covariate-adjusted β values (GDM – non-GDM) at cg08407434 show consistent increased methylation in siblings exposed to GDM. (B) Within sibling pair differences of covariate-adjusted overall methylation levels of all CpG sites across the HNF4A region also show a marginal trend of increased methylation in siblings exposed to GDM (P = 0.053). (C) Overall DNA methylation level of HNF4A showed a negative correlation with mRNA expression albeit non significant (Pearson's coefficient: -0.364, P = 0.068).

Discussion

In this study, we performed epigenome-wide association analysis in 18 sibling pairs discordant for their exposure to maternal GDM pregnancy. We identified 12 CpG sites from offspring peripheral leukocytes suggested to be associated with exposure to maternal GDM (q < 0.15). These CpG markers include those located at HNF4A and RREB1, which are known for its association with diabetes and obesity, respectively. We found that there were higher methylation levels overall in the HNF4A region in the affected siblings. It showed a negative correlation with its mRNA expression levels, though non significant. Several pathways relevant to metabolic regulation were enriched by genes harboring differentially methylated CpG markers. To the best of our knowledge, this is the first study to investigate discordant sibling pairs to identify DNA methylation markers that are associated with exposure to maternal GDM pregnancy. This unique study design enables us to assess effects of intrauterine hyperglycemia by reducing genetic and environmental confounding effects on DNA methylation.

It has been previously reported that siblings born after the mother developed T2D had significantly elevated risk of developing T2D with odds ratio of 3.7 compared with their sibling born before the mother developed T2D.⁹ In addition, BMI of siblings born after the mother developed T2D was 2.6 kg/m² higher.⁹ However, differences in epigenetic markers between siblings discordant for their exposure to GDM have not been reported previously. Unsupervised hierarchical clustering showed that one large cluster consisted mainly of offspring exposed to GDM. However, some sibling pairs were clustered into the same cluster, not separated according to their exposure status. These findings imply that both sibship and exposure to GDM should be considered together.

There were several characteristics of the study participants that are worthwhile to note. We investigated offspring of GDM, which is usually a milder form of hyperglycemia compared with pregestational diabetes. GDM usually develops in the third trimester and about 20–40% of GDM women require insulin therapy. On the other hand, women with pregestational diabetes have hyperglycemia at the beginning of conception and the degree of hyperglycemia is more severe. It would be of interest to investigate differential DNA methylation as well as phenotypic differences between offspring born to GDM and offspring born to maternal pregestational diabetes. Another point that should be noted is that children rather than adults were investigated for DNA methylation. We were interested in DNA methylation markers that could potentially mediate the effect of intrauterine hyperglycemia on later life metabolic dysregulation. Longitudinal investigation showed that in utero DNA methylation could be altered during childhood and adolescence.²⁵ Whether the difference in DNA methylation identified in our study persists in adulthood needs to be further studied.

There were several CpG markers of potential interest. The most significant differentially methylated CpG site was located within CTNND2, which encodes δ-catenin and is implicated in development of brain and eye. It is overexpressed in cancers and plays a role in disruption of adherence junction and tumorigenesis.²⁶ However, this gene has not been implicated in pathogenesis of diabetes or glucose homeostasis before.

The second most significant differentially methylated CpG site (cg08407434) in Model 1 was located at intron of HNF4A. This gene serves as a transcriptional factor that regulates pancreatic β cell development and differentiation.²⁷ It has been shown in rodents that maternal nutrition and aging can induce epigenetic change in HNF4A and regulate its expression level in offspring islets.¹⁰ When monozygotic twins discordant for T2D were investigated, there was increased methylation of HNF4A in adipose tissue of T2D individuals.²⁸ We showed that increased methylation of this marker in siblings exposed to maternal GDM was consistent in all pairs. However, the difference was small, and it is not known whether this small consistent change in DNA methylation can be biologically meaningful. Further replication studies are required to validate this finding. In addition, a follow-up study on how these small yet significant alterations persist through lifetime and affect metabolic health in later life would be of interest.

Another differentially methylated CpG site of potential interest is located within the gene body of RREB1. It encodes a zinc finger transcription factor that binds to RAS-responsive elements. Genetic variants in RREB1 were associated with several metabolic traits including waist-hip-ratio,²⁹ uric acid level,³⁰ and fasting glucose concentrations.²³ In addition, a previous report showed that RREB1 plays an important role in transcriptional activation of NeuroD1 target genes, such as insulin and glucokinase, in pancreatic endocrine cells.³¹ It would be worthwhile to investigate functional consequences of increased methylation of this gene on metabolic traits.

In the gene set enrichment analysis using genes that harbor differentially methylated CpG sites, we identified multiple important pathways relevant to GDM and T2D. The ‘regulation of insulin like growth factor transport and uptake by insulin like growth factor binding proteins’ was the most significantly enriched pathway. Insulin like growth factor axis is important in fetal growth and glucose homeostasis. There are lines of evidence that it also plays an important role in the pathogenesis of diabetes.³² It has been shown that levels of insulin like growth factor and their binding proteins are altered in T2D and obesity.³³ In addition, it has been reported that GDM is associated with lower insulin-like growth factor binding proteins in maternal and cord blood.³⁴ The ‘superpathway of conversion of glucose to acetylCoA and entry into the TCA cycle’ is crucial in glucose homeostasis and energy metabolism. Alterations in this pathway may result in T2D and related metabolic abnormalities. Enrichment of this pathway is in accord with the hypothesis that intrauterine exposure to GDM increases risk of metabolic disorders in offspring via epigenetic changes. However, none of the pathways reached multiple testing-corrected significance of FDR <0.25 probably due to small sample size. Pathways with a small gene set are prone to inflated scores and, thus, further analyses and cautious interpretation are warranted.

There are several limitations in this study. First, despite the uniqueness and strengths of the study design, the sample size was small and statistical power was limited. The threshold for statistical significance was relaxed in this study and it might have introduced false positive findings. Therefore, the results should be interpreted with consideration of possible type I errors. As this study had limited statistical power, we focused on candidate regions that had biologic plausibility, such as HNF4A and RREB1. Second, we used peripheral leukocyte samples to assess differential DNA methylation associated with intrauterine exposure status. Early-life exposure can affect DNA methylation in a wide range of tissues that can persist throughout the lifetime while later-life exposure can affect only one or few tissues.^11,35−37 This enables DNA methylation of peripheral leukocyte samples to be used as potential markers of intrauterine exposure. There are also reports on consistent DNA methylation changes in leukocytes as a predictor of development of T2D.¹³ Third, there was a significant difference in offspring ages and their maternal ages at gestation between cases and controls. Most of the participants’ mothers experienced GDM pregnancy after having had a normal pregnancy. This is in accord with the fact that advanced maternal age is a major risk factor for GDM. To take into account the age difference between groups, we included offspring age as a covariate in analyses. According to our further analysis on relationship between age and DNA methylation at top CpG markers, offspring age was significantly correlated with DNA methylation at several markers including cg04988367, cg13911801 and cg09304293 (Supplementary Figure S3). However, this trend tended toward null when the analysis was performed within each group. Meanwhile, for cg08407434 located within HNF4A, there was no significant association with offspring age or maternal age.

In summary, we found several suggestive epigenetic markers, including cg08407434 located at HNF4A, which are associated with intrauterine exposure to maternal GDM. This supports the hypothesis that intrauterine exposure to diabetic environment can induce epigenetic changes in offspring. Further large scale replication studies and functional investigations are required to validate our findings.

Materials and methods

Study design and population

This study examined 18 discordant sibling pairs to investigate DNA methylation markers associated with intrauterine exposure to GDM. In each family, an affected sibling was conceived during GDM pregnancy and the control sibling during non-GDM pregnancy. Each of the sibling pairs was full-siblings and the main difference was their exposure to intrauterine hyperglycemia. All the mothers had at least one GDM pregnancy and delivered healthy singleton full-term neonates at Seoul National University Bundang Hospital between December 2003 and August 2006. GDM was diagnosed by a 2-step approach using 50 g glucose load test (cut-off >140 mg/dL) followed by 100 g OGTT with glucose thresholds established by Carpenter and Coustan.³⁸ None of the mothers had overt diabetes before pregnancy. Sibling pairs of these mothers were recruited to be enrolled in this study when the offspring born during GDM pregnancy was at age 5–8 y. The Institutional Review Board of the Seoul National University Bundang Hospital (IRB No. B-1007–105–007) approved this study and written informed consent was obtained from the parents and children.

DNA methylation analysis

Genomic DNA was extracted from the peripheral blood leukocytes of the offspring immediately after sampling using Gentra Puregene Blood Kit (Qiagen, Hilden, Germany). The DNA samples were treated with sodium bisulfite using Zymo EZ-96 DNA methylation kit (Zymo Research, Irvine, CA, USA) and applied to the HM450 K (Illumina, San Diego, CA, USA) assay. The raw intensity data were processed using an R package RnBeads.³⁹ The data can be obtained from the GEO accession number GSE102177. Methylation status at each of the CpG sites was quantified as a β value, which is defined as the ratio of the methylated probe intensity and the total intensity, ranging from 0 (fully unmethylated) to 1 (fully methylated). CpG probes with a detection P > 0.01, associated with single nucleotide variations, outside of CpG context or located at sex chromosome were excluded from analysis. β-mixture quantile normalization (BMIQ)⁴⁰ was performed to reduce technical differences between 2 different bead chemistry of HM450K array (Supplementary Figure S1.). White blood cell composition was inferred using the method proposed by Houseman et al.⁴¹ to account for cell type heterogeneity in assessing differential DNA methylation.

Epigenome-wide association analysis was conducted to identify differentially methylated CpG sites associated with exposure to maternal GDM pregnancy. M-values ($\mathrm{M}\mathrm{ }={\log }_2\left(\frac{\mathrm{\beta }}{1-\mathrm{\beta }}\right)$) which are logit transformed β values were used to account for heteroscedasticity and to allow the assumption of Gaussian distribution.⁴² An empirical Bayes approach as implemented in the R package limma⁴³ was applied. M-values were fitted using a design matrix in which age, sex, experimental batch, GDM exposure status, and sibling IDs (Model 1), and additionally estimated cell type composition (Model 2), were defined. Pairwise contrasts between GDM and non-GDM samples were made at each of the CpG sites after adjusting for confounding factors of age, sex, experimental batch and additionally leukocyte composition. P values were estimated using moderated t-statistics⁴³ and false discovery rate (q value) was used to adjust for multiple comparison.⁴⁴ The same statistical model and estimation method were applied to compare overall DNA methylation of HNF4A between cases and controls. To make comparisons more interpretable between affected and unaffected siblings, we used DNA methylation levels that were adjusted for known confounders of age, sex, and experimental batches.

Unsupervised hierarchical clustering analysis was then performed using covariate-adjusted DNA methylation values of the 100 most significant probes. Distance between all pairs of the samples was estimated using Manhattan distance. Clustering was performed using average linkage algorithms.

Pyrosequencing analysis

Bisulfite pyrosequencing was performed to validate the finding of differential methylation of the CpG site cg08407434. Each primer was designed using Pyrosequencing Assay Design Software v2.0 (Qiagen). The primers are provided in Supplementary Table 1. PCR reaction was performed and confirmed by electrophoresis staining. ssDNA template was prepared from 16–18 µL biotinylated PCR product using streptavidin Sepharose HP beads (Amersham Biosciences, Sweden) following the PSQ 96 sample preparation guide. The respective sequencing primers (15 pM) were added for analysis. Sequencing was performed on a PyroMark ID system with the Pyro Gold reagents kit (Qiagen) according to the manufacturer's instruction without further optimization. The methylation level was quantified by the average of the degree of methylation at the CpG sites formulated in pyrosequencing.

Gene expression profiling

Buffy coat was snap-frozen and stored in liquid nitrogen until further processing. Using the TRIzol LS Reagent (Thermo Fisher Scientific, Waltham, MA, USA), RNA was isolated according to the manufacturer's protocol. Gene expression levels were investigated using Agilent SurePrint G3 Human GE 8 × 60K Microarray (Agilent Technologies, Santa Clara, CA, USA). A total of 27 samples were available for analysis with RNA integrity number 2.5 or above after excluding 9 samples (5 for sample unavailability and 4 for poor RNA quality). Raw data of the microarray was loaded and processed using the limma R package.⁴³ We applied background correction using a convolution of normal and exponential distributions and normalization using the Cyclic Loess method.⁴⁵ Irregular replicate probes were averaged and log₂-transfromed data was used for analysis.

Gene set enrichment analysis

A total of 8,545 genes with at least one nominally significant differentially methylated probe (P<0.05) were selected. These genes were ranked by their mean difference of β values using the most significant CpG site directly linked with the gene. Gene set enrichment analysis⁴⁶ was performed using the ranked list of genes and predefined 4,307 pathway gene sets downloaded from the ConsensusPathDB (http://cpdb.molgen.mpg.de/).²⁴

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

We are grateful to Ministry of Health & Welfare, Seoul National University Hospital and Seoul National University Bundang Hospital for giving us grants. We also thank the participants of this study.

Additional information

Funding

This study was supported by grants from the Korean Health Technology R&D Project, Ministry of Health & Welfare, Republic of Korea (Grant No. HI15C3131), Seoul National University Hospital (Grant No. 30–2015–0240) and Seoul National University Bundang Hospital (Grant No. 14–2014–030).