Gestational diabetes and maternal obesity are associated with epigenome-wide methylation changes in children

Gestational diabetes mellitus (GDM) is associated with increased risk of type 2 diabetes (T2D) and associated cardiometabolic diseases in the offspring (1–9). In a recent substudy (n = 1,234) of the Danish National Birth Cohort (DNBC), we found that offspring of women with GDM, already at the preadolescent age, displayed increased adiposity and associated dysmetabolic traits (10). However, the extent to which dysmetabolic traits and increased risk of T2D among offspring of GDM women are due to genetic variation, shared environment, and/or epigenetic variation remains unclear. Considerable data support the latter possibility, revealing the time in utero as a sensitive period of susceptibility linking adverse environmental exposure to adverse offspring health outcomes, even much later in life (11).

Epigenetic mechanisms, including DNA methylation, histone variants/modifications, chromatin-modifying proteins, and noncoding RNAs, work in combination to regulate gene expression. Accumulating evidence suggests epigenetic dysregulation of gene functions in early life as a link between prenatal environment and later metabolic disease (12–16). Indeed, previous smaller studies of offspring of women with GDM reported differential DNA methylation changes in cord blood or placenta tissue samples, using either a candidate gene approach (17–23) or in epigenome-wide approaches (24–26). Interestingly, the epigenome-wide association studies (EWAS) consistently suggested methylation changes in genes involved in metabolic pathways among GDM versus control offspring. Despite this, results have been generally inconclusive due to lack of (a) statistical power due to insufficient sample size, (b) consideration of potential confounders, such as prepregnancy maternal BMI (mBMI), and (c) replication studies.

For environmentally mediated epigenetic variation arising in utero to be considered as potential biomarkers for later-onset disease, it is a great strength if they are persisting and detectable postnatally. Although showing some tissue specificity, DNA methylation profiles in blood often reflect epigenetic fingerprints in other tissues (27, 28). To this end, there is a great potential to use DNA methylation in blood as biomarkers to predict later disease development. In support of this idea, studies of blood cell DNA have showed that individuals who were prenatally exposed to famine during the Dutch hunger winter presented altered methylation levels 6 decades later, at several genes including IGF2, when compared with their unexposed, same-sex siblings (13, 29). Recently, a study investigating DNA methylation in peripheral blood from adolescent Pima Native Americans reported genome-wide methylation differences in offspring of mothers with T2D prior to pregnancy, which were not associated to prepregnancy mBMI (30). Whether pregnancies affected by GDM have similar long-lasting effects on the offspring epigenome remains unclear. Using a genome-wide discovery approach in a nested case-control subcohort (n = 188) and subsequent replication of most promising changes in the remaining cohort (n = 905), we have identified DNA methylation in peripheral blood of 9- to 16-year-old offspring of GDM women relative to non-GDM controls (Figure 1).

Figure 1

Overview of study design.

Maternal, birth, and offspring characteristics of the DNBC GDM subcohort. Characteristics of the study population are outlined in Table 1. GDM mothers had a higher prepregnancy mBMI but a smaller pregnancy weight gain (2 kg less) than controls (P < 0.001). They were also, on average, 1 year older, more likely to be smokers, and among lower socio-occupational status classes compared with controls (P ≤ 0.015). Furthermore, the outcome of pregnancies in women with GDM showed a higher proportion of cesarean section delivery, shorter gestational age (GA), higher birthweight, and a shorter period of breastfeeding when comparing to the controls (P < 0.001).

Table 1

Maternal, paternal, birth, and offspring characteristics

As previously described (10), offspring of GDM pregnancies showed a range of adverse metabolic characteristics, including increased fasting plasma glucose, insulin, and C-peptide levels; HOMA-IR; total body fat percentage (total BF%); BMI; and resting heart rate (HR) — as well as a decreased total muscle mass (Table 1).

The clinical characteristics of the subgroup of 188 offspring included in the epigenome-wide analysis (n = 93 GDM and 95 control offspring) are shown in Supplemental Table 1 (supplemental material available online with this article; https://doi.org/10.1172/jci.insight.122572DS1). Generally, most of the clinical and metabolic characteristics in the discovery subgroup were representative of the entire cohort of 1,234 offspring (n = 608 GDM and 626 control offspring), as well as the maternal characteristics including prepregnancy mBMI (Supplemental Table 1 and Table 1).

Epigenome-wide DNA methylation analysis of blood from GDM versus control offspring. To investigate whether DNA methylation in blood differed in 9- to-16-year-old GDM versus control offspring, epigenome-wide DNA methylation was analyzed in a subgroup of 93 GDM and 95 control offspring using the Infinium HumanMethylation450 BeadChip (HM450K array). In total, 425,344 probes (referred to as CpGs) on the HM450K array were retained for subsequent analysis.

Analysis using the removal of unwanted variation for methylation data (RUVm) method from the missMethyl Bioconductor package, with a FDR less than 5% (P < 0.05), identified 76 significantly differentially methylated CpGs between GDM offspring and controls not including any covariates (Table 2). The methylation differences ranged from ~1%–5.2%, and the majority of the 76 CpGs (92%) showed lower DNA methylation in the GDM offspring (Table 2). The genomic location of these sites varied, with 58 annotated to 56 unique genes, 2 of which (cg09452568 and cg00992687) mapped to ESM1 (located in close proximity in exon 3) and another 2 (cg19739596 and cg14328641) mapped to MS4A3 (located in the distal promoter and exon 1, respectively). A total of 5 CpGs were located within a CpG island, 9 at a CpG island shore, 11 within a CpG island shelf, and the remaining 51 in “open sea” regions. Of the 76 CpGs, 42% were localized to gene body regions (Figure 2). The full ranked list of differentially methylated CpGs, including their genome location, is provided in Supplemental Table 2.

Figure 2

Genomic features and distribution of the 76 differentially methylated CpGs.

Table 2

List of differentially methylated CpGs between GDM offspring and controls

GDM-associated methylation differences remained after the adjustment of offspring BMI. Multivariate linear regression analyses were applied to address potential confounding effects of offspring sex, age, and adiposity on the association between GDM and offspring methylation level of the 76 differentially methylated CpGs from the discovery study. Overall, when only adjusting for offspring age and sex, the same group-dependent differences in methylation as found in the RUVm analysis were observed for all CpGs, although P values for significance were higher (Supplemental Table 3, Model 1). Importantly, when offspring BMI was included as a covariate together with offspring age and sex, the significant differences between GDM offspring and controls remained for all 76 CpGs (Supplemental Table 3, Model 2). In contrast, differential methylation at cg13246235 (PHACTR1), cg25348105 (BAHCC1), and cg00042683 (SHF) were found to be associated with offspring age (P ≤ 0.04) (data not shown). None of the CpGs were significantly associated to sex (data not shown).

Prepregnancy mBMI confounded methylation differences between GDM and control offspring. We next investigated whether prepregnancy mBMI confounded the differences observed in methylation between cases and controls. When prepregnancy mBMI was included as a covariate, the significant differences between cases and controls became statistically insignificant for 63 of 76 CpGs. Instead, the majority of these CpGs showed significant associations to their mothers’ prepregnancy BMI (Supplemental Table 3, Model 3; and Figure 3). The remaining 13 CpGs remained significantly associated to GDM status. The genes associated with these CpGs are listed in Figure 3. We further tested whether the following potential maternal and birth-related confounders were associated with methylation degree, also under adjustment for GDM status and offspring age and sex: maternal age, pregnancy weight gain, maternal smoking during pregnancy, GA, cesarean section, birth weight, and breastfeeding. None of these factors were associated with offspring methylation degree at any of the 76 CpGs (data not shown).

Figure 3

Maternal GDM and/or obesity -associated CpGs/genes. Genes and CpGs associated to GDM and/or maternal prepregnancy BMI after adjustment for GDM status, prepregnancy mBMI, and offspring sex and age (P ≤ 0.05), using multivariate, linear modeling. ^aAssociation to GDM status at P ≤ 0.10; ^bassociation to prepregnancy mBMI at P ≤ 0.10. Underline represents CpG sites that were used in subsequent validation/replication analyses. Multivariate linear regression analysis, adjusting for GDM status, offspring age, sex, and prepregnancy maternal BMI; n = 175.

Biological relevance of differential methylation in GDM offspring. To further understand the functional relevance of the identified differentially methylated CpGs, a literature search was performed to investigate potential involvement of the 76 CpGs in metabolism, environmental exposures, and related diseases. The search was performed in the PubMed database (https://www.ncbi.nlm.nih.gov/pubmed/), using each CpG number or gene name and the following terms: DNA methylation or metabolism. This showed that either the CpG DNA methylation or expression of the associated gene were for 9 of 76 (11.8%) previously associated with T2D, obesity, diabetic nephropathy, or coronary heart disease. These CpG-associated genes included PDE6A, PRKCZ, PVT1, GALNT2, MS4A3, IL1RN, and BTD (Supplemental Table 2).

Pathway analysis. In order to gauge some further insights into the cellular pathways showing epigenetic variation in association with exposure to maternal GDM, the list of CpGs showing differential methylation was expanded to include those with an unadjusted threshold of P < 0.0001 in the RUVm analysis. This produced a list of 292 CpG sites, annotated to 176 genes. Using ingenuity pathway analysis (IPA; http://www.ingenuity.com), these were enriched for 5 important functional networks, with lipid metabolism being the highest ranked (Table 3).

Table 3

Ingenuity Pathway Analysis

Technical validation of 3 differentially methylated CpGs from the discovery study. To technically validate the results of the discovery study, we performed pyrosequencing of 4 CpGs including cg00992687 and cg09452568 (ESM1), cg14328641 (MS4A3), and cg09109411 (PDE6A) within the 186 offspring from the discovery study. Selection criteria of these CpG sites are explained in detail in Methods.

Direct comparison of the HM450K array and pyrosequencing methylation data of the 4 CpGs (Supplemental Figure 2) revealed highly significant positive correlations (P < 0.0001, Pearson’s r coefficient ≥ 0.91), highlighting the sensitivity and reproducibility of the HM450K approach (Supplemental Figure 3, A–D). Three additional CpG sites in the vicinity of cg00992687 (ESM1), 1 CpG site in the vicinity of MS4A3, and 3 additional CpG sites in the vicinity of PDE6A were also analyzed in the pyrosequencing assays, revealing trends toward lower methylation in association with the GDM, in accord with the HM450K array discovery data (Supplemental Figure 4, A–C).

Replication of 3 differentially methylated CpGs identified in the discovery study. For biological replication of the results of the discovery study, we performed pyrosequencing of cg00992687 and cg09452568 (ESM1), cg14328641 (MS4A3), and cg09109411 (PDE6A) in the remaining 905 offspring from the original cohort not included in the discovery approach. We performed multivariate analysis with adjustment for GDM status, mBMI, and offspring sex and age (n = 445–788). These analyses were performed independently on pyrosequencing data for 10 CpG sites: 4 identified from the discovery study and 6 flanking sites.

Methylation at cg09452568 (ESM1) and the adjacent 2 CpG sites was significantly ~0.3%lower per each kg/m² increase in prepregnancy mBMI (P ≤ 0.03) (Table 4). Cg00992687 (ESM1) was also borderline significantly 0.2% lower per each kg/m² increase in prepregnancy mBMI (P = 0.08). A similar pattern was present for cg14328641 (MS4A3) and the adjacent position 2 CpG, which both decreased significantly by 0.2% per each kg/m² increase in prepregnancy mBMI (P ≤ 0.03). No independent association between GDM status and methylation was found for either ESM1 or MS4A3 CpGs sites (P ≥ 0.51). Interestingly, we observed that methylation at cg09109411 (PDE6A) and all 3 adjacent CpGs were significantly associated with GDM status, increasing by ~2% in GDM offspring. Only 1 site (position 1) of PDE6A showed weak evidence of lower methylation with increasing mBMI (P = 0.06) (Table 4).

Table 4

Association between maternal GDM, prepregnancy BMI, and offspring methylation degree at ESM1, MS4A3, and PDE6A CpG sites

In addition, we observed that methylation at both ESM1 and MS4A3, but not the PDE6A, were associated with sex, and all were lower in females (P ≤ 0.01). Methylation at all 4 ESM1 sites were negatively associated with age (P ≤ 0.02) (Table 4).

To our knowledge, this is the largest genome-wide DNA methylation study in GDM offspring and controls yet performed. We showed that exposure to GDM is associated with DNA methylation variation in 9- to 16-year-old offspring. The list of identified differentially methylated genes includes several genes previously associated to metabolic disease. Interestingly, the majority of differentially methylated CpGs among offspring of women with GDM appear to be confounded by and associated with prepregnancy mBMI rather than GDM in isolation, suggesting a role of maternal obesity in programming of epigenetic traits among offspring of women with GDM in pregnancy. Three of these mBMI-associated CpGs identified in the discovery cohort (cg00992687 and cg09452568 of ESM1 and cg14328641 of MS4A3) were replicated in our larger cohort (n = 445–788), while cg09109411 (PDE6A) was found to be specifically associated with GDM status. Furthermore, our results show that the maternal exposures of GDM and obesity occurring 9–16 years ago are more associated with offspring DNA methylation profile than the offspring’s own BMI at the adolescent age. Lipid metabolism and endocrine disorders were among the most significant biological functions identified in the pathway analysis of the genes harboring the differentially methylated CpGs, supporting the notion that GDM-associated methylation may play an important role in transmission of metabolic risk from GDM mother to offspring. Indeed, this is consistent with the mounting data showing increased risk of T2D, obesity, or associated cardiometabolic disease in offspring affected by diabetes in pregnancy (1, 31, 32).

Previous genome-wide methylation studies in GDM offspring performed in birth samples (cord blood and placenta) and in 8- to 12-year-old offspring (24–26) did not examine the potential confounding effects of maternal obesity, nor the infant’s or child’s own adiposity. Our results for ESM1 and MS4A3 from the discovery study, which were confirmed in the replication study, suggest that prepregnancy mBMI, rather than GDM in isolation, is associated with the offspring epigenetic profile. Indeed, in our previous study of adult offspring of women with and without GDM, we observed that prepregnancy mBMI was associated with PPARGC1A DNA methylation in s.c. adipose tissue of the offspring (33). A recent PACE consortium meta-analysis found small effects of prepregnancy mBMI on cord blood methylation, which were mainly explained by cell heterogeneity and/or confounding by shared maternal-child environment. Still, several associations between prepregnancy mBMI and methylation at birth were replicated in an independent large sample of adolescents (34). None of the CpGs associated with GDM or mBMI in our study were associated with mBMI in the PACE study. However, our study participants were recruited based on GDM status; therefore, different methylation changes may be in play, compared with the PACE study, where recruitment were based on maternal obesity. Furthermore, the stressful perinatal environment may cause acute and more transient methylation changes that make the true programming changes undetectable. However, interestingly, we found that 2 of the 76 CpGs (PRKCZ, a serine-threonine kinase involved in proliferation, differentiation and secretion, and PRKAR1B, the regulatory subunit of cyclic AMP-dependent protein kinase A that is involved in signaling pathways of many cellular events, including ion transport, metabolism, and transcription; refs. 35 and 36) were both also reported to be differentially methylated in cord blood by Quilter et al. (24).

Previous studies of GDM-associated DNA methylation in offspring have been inconclusive, likely due to lack of power, but they all indicate methylation changes associated with genes involved in metabolic pathways (24–26). Indeed, our pathway analysis results suggest that exposure to GDM or maternal obesity impacts multiple gene-associated networks, including lipid metabolism, developmental disorders, endocrine disorders, and organ development. These incorporate key genes to a number of intracellular signaling pathways including the TREM1 signaling pathway involved in inflammatory responses through monocyte and neutrophil activation, as well as genes involved in the signaling pathways of IL-1, a proinflammatory cytokine produced, for example, by activated macrophages and lymphocytes. Furthermore, key genes in GPCR-mediated nutrient sensing in enteroendocrine cells and genes involved in sphingosine-1-phosphate signaling important in regulation of embryogenesis, cardiogenesis, and vascular development were also affected, further suggesting that our findings of differentially methylated CpGs offer a window into pathophysiological processes. In addition, Hajj et al. reported endocrine system disorder and developmental processes as top disorders and biological functions affected (20), and Ruchat et al. reported the top pathways gastrointestinal disease, metabolic disease, and endocrine system disorder affected (25) — all pathways that were replicated in our study. The pathway networks — tissue morphology, cellular maintenance, and organismal injury — were also among top hits in our analysis and may be speculated to be a reflection of DNA damage and repair mechanisms, which could be associated with diabesity. Functional studies are required to further examine this association.

By exploring the chromatin states of 76 CpGs using ENCODE data, we found high presentation of promoters and regulatory elements. This further suggests that the differentially methylated CpGs may have functional importance in regulating transcription; unfortunately, RNA was not available in the present study to explore this hypothesis further.

Several of the 76 CpGs from the discovery study have previously been reported to show differential methylation in association with risk for T2D, obesity, associated cardiometabolic diseases, or aging. These include the top-ranked cg09109411, which is localized to PDE6A, encoding the cyclic-GMP–specific phosphodiesterase 6A α subunit, involved in transmission and amplification of the visual signal (37). This CpG was previously reported to show decreased methylation in association with obesity shown in healthy young adult female and male monozygotic twin pairs discordant for BMI (n = 30 twin pairs) (38). Interestingly, in our discovery study, we also find decreased methylation among the GDM offspring. The same PDE6A CpG site (cg09109411) has also been shown to increase in methylation degree with GA, with preterm infants showing lower methylation levels (39). Furthermore, the second-ranked CpG (cg11638399) was also found in the above-mentioned twin study to exhibit decreased methylation in the obese group (38), with the same direction of methylation change as we report in the GDM offspring. Cg19739596 (MS4A3, encoding a member of the membrane-spanning 4A gene family, involved in immune signal pathways) showed decreased methylation in adipose tissue of T2D subjects (40), again in the same direction as we found in the GDM offspring. Moreover, cg26930596 (PRKCZ) methylation degree in CD4⁺ T cell DNA has been associated to exposure to sunlight (41), enforcing that methylation degree of this site is influenced by environmental factors. Additionally several CpGs in the PRKCZ promoter have been shown to be hypermethylated in leucocyte DNA from adult T2D patients (35). Of the majority of the remaining CpGs, the functionality related to DNA methylation of the associated gene has not yet been recognized.

There are limitations to this study. DNA from peripheral blood samples consists of more than 1 cell type; nonetheless, we did not observe any differences in cell composition. In this study, we have only measured the average DNA methylation degree across all the different blood cell types of the samples collected. Therefore, the individual methylation changes of the specific cell types may most likely be of a higher extent than the average differences we can report. A small methylation difference in blood might reflect a larger difference in more metabolic relevant tissue, which is known from previous studies (27). Also, it is unknown to which extent an effect of 1%–5% change in methylation may have across a whole lifetime or how important the accumulative effects are. Finally, it is important to recognize that methylation of the circulating heterogeneous mixed blood cells at a single point in time may not reflect a fixed methylation change at any other potentially more important tissue level. Investigating methylation changes at such prime tissues in large populations such as the current is, however, not feasible, especially in childhood studies.

The potential exists that these changes reflect target organ disease mechanisms in diabetes and/or are suitable as biomarkers to predict diabetes. However, validation of all 76 differentially methylated CpGs from the discovery study would be needed to state that the associations with the offspring dysmetabolic phenotype in the larger cohort are associated with methylation changes of the enriched pathways found in the discovery cohort.

Unfortunately, we have in this cohort no information on the severity of GDM and how it was treated, and we acknowledge that treatment-related differences cannot be excluded as a potential confounder. Finally, in this study design, it is not possible to conclude on the stability of these CpGs; longitudinal studies are required for this. Further validating in cord blood samples or across tissue types would be beneficial in elucidating possible roles of these CpGs in the disease risk of GDM offspring.

Conclusions. Our results suggest that GDM and prepregnancy mBMI have long-term effects on offspring epigenetic profile, which supports that DNA methylation is involved in fetal metabolic programming. We provide evidence of potential postnatal stability of previously published differentially methylated CpGs shown at birth among GDM offspring. Further longitudinal studies, commencing prior to birth, with detailed clinical information and biospecimens are required to unravel the potentially complex interplay between prepregnancy mBMI and GDM in both offspring metabolic health and epigenetic profile. The potential exists that methylation changes in circulating blood cells reflect target organ disease mechanisms in diabetes and are suitable as biomarkers to predict metabolic disease in offspring of women with GDM.

Study cohort. The DNBC enrolled 91,827 women (101,045 total pregnancies) between January 1996 and October 2002, of which 1,350 had GDM. This GDM subgroup, together with a random control group consisting of 2,629 pregnancies within DNBC, was established to study the long-term health implications of glucose intolerance in pregnancy (42). GDM cases were defined by a GDM-related diagnosis recorded in the Danish National Patient Registry (International Classification of Diseases–10 [ICD-10] codes O24.4 and O24.9) and/or a self-reported GDM event from study interviews at 30 weeks of gestation or 6 months postpartum. The diagnosis of suspected GDM in the DNBC has previously been described (43).

Between March 2012 and May 2014, 608 offspring of mothers with GDM and 626 offspring of mothers from the control group not affected by GDM were included in the study. Clinical examinations and biospecimen collection was conducted at local temporarily established clinics across all major regions of Denmark.

Clinical examinations. The offspring of GDM women and controls underwent a clinical examination including anthropometry measurements, puberty Tanner status examination (44), blood pressure, HR, and fasting blood samples for later measurement of biomarkers including glucose, insulin, C-peptide, and lipid profiles. Standard assays were used for biomarker analysis as previously described (10). HOMA-IR was calculated as: ([(fasting plasma insulin × fasting plasma glucose)/22.5] × 0.144). Study participants examined at the Copenhagen University Hospital (n = 666) were also studied for body composition using Dual-Energy X-ray Absorptiometry (DXA) scanning (Lunar, Scanex Medical Systems).

DNA sample collection and extraction. A strictly enforced standardized sample collection procedure was used across all study examination sites to secure comparable results. EDTA tubes (10 ml, BD Biosciences) were placed directly on ice after collection of venous blood and then, within 10 minutes, centrifuged for 10 minutes at 4^oC, at 1,100 g (Eppendorf Centrifuge, 5810R, MilliporeSigma). Buffy coat was prepared from EDTA tubes and stored immediately at –20^oC with subsequent transfer (within 8 hours) to –80^oC until later DNA isolation. Buffy coats were collected from 536 GDM offspring and 555 control offspring: in total 1,091 offspring. DNA samples were extracted from all 1,091 buffy coat samples using the QIAamp 96 DNA blood kit (Qiagen).

DNA methylation profiling. We designed our DNA methylation study as a genome-wide discovery study performed in a subgroup of GDM and control offspring, followed by a technical validation in the same samples, and finally biological validation in the remaining samples from the entire cohort (Figure 1). Epigenome-wide DNA methylation measurement was performed using the Illumina’s Infinium HumanMethylation450 BeadChip (HM450K array) on a subset of 95 GDM offspring (46 female, 49 male) and 95 control offspring (44 female, 51 male), matching age, sex, and time of the year of the clinical examinations. The discovery cohort was selected carefully to ensure the GDM diagnosis in the GDM offspring group and, furthermore, to minimize the risk of unrecognized glucose intolerance/diabetes in the mothers of the control offspring group. Thus, GDM mothers were defined using strict criteria of (a) fasting glucose level ≥7 mmol/l, (b) a 2-hour oral glucose tolerance test (OGTT) glucose level ≥7.8 mmol/l, and/or (c) with a clinicians verified diagnosis of GDM. GDM offspring of mothers who had either prediabetes or diabetes based on the OGTT at time of the follow-up study were prioritized. Additionally, we only included offspring of control mothers with no prediabetes or diabetes at the time of the follow-up study in order to reduce risk of including false negatives.

DNA from the 190 buffy coat samples was processed at Service XS (Leiden) for methylation analysis. To minimize batch effects, samples were randomized per study group, age, and sex, as well as GA, previously shown to influence epigenetic variation (39). Bisulfite (BS) conversion of 500 ng DNA was performed (EZ-DNA Methylation Gold kit, Zymo), and conversion quality control (QC) was conducted by quantitative PCR (qPCR) and melting curve analysis.

BS-converted DNA was hybridized to the Illumina BeadChips according to the Infinium Protocol (Illumina). The BeadChip images were scanned using the iScan system, and raw β values were exported as iDAT files with Methylation Module (v.1.9.0, Illumina), using default analysis settings.

Preprocessing and analysis of HumanMethylation450 BeadChip Arrays. Raw methylation data was preprocessed using the Minfi (45) package from the Bioconductor project (http://www.bioconductor.org) in the R statistical environment (http://cran.r-project.org, version 3.0.2). QC was performed by assessment of sample-dependent and sample-independent control probes for BS conversion, staining, hybridization, extension, and specificity. Probes with detection P > 0.01 in 1 or more samples and samples with a mean detection P > 0.05 were removed prior to statistical analysis. In total, 2 samples were removed, as they did not pass 1 or more QC measures. Data filtering was conducted to remove: (a) X/Y chromosome probes that may lead to sex-specific bias, (b) cross-reactive probes (>47 bases matched to more than 1 genomic target), (c) probes containing a single nucleotide polymorphism (SNP) at the interrogated CpG or the single-base extension site, and (d) non-CpG targeting probes (46). Following QC, the size of the final data set was 188 samples (93 GDM offspring, 95 control offspring) with methylation data from 425,344 probes (referred to as CpGs) across all samples. Cell count estimation was performed empirically using surrogate variable analysis (47, 48), by the “estimate cell composition” function in the Minfi package. No differences in cell composition were identified, and cell composition was deemed unlikely to be a confounder; therefore, we did not correct for cell composition in following analyses (Supplemental Figure 5).

Following data preprocessing, β values were converted to M-values (M = log₂ [unmethylated signal intensity/methylated signal intensity]) (49). We performed multidimensional scaling analysis (MDS) and principal component analysis (PCA) to examine the sources of variation in the data and to determine whether these correlate with participant factors of interest or known covariates, which included offspring sex and age, GA, birthweight, and maternal age.

To identify differentially methylated CpGs between GDM and control offspring, the 2-step RUVm method from the MissMethyl Bioconductor package was used (50, 51). Briefly, in the first step, the 613 Illumina negative controls were used in a linear regression to derive the components of unwanted variation that should be adjusted, such as technical variation including batch, slide, and array effects. By including these in the first pass analysis, we reduced that type of unwanted variation, which improves the ability to rank the probes with regard to their association with the factor of interest: GDM status. However, this initial analysis does not capture unwanted biological variation. Therefore, in the second step of the analysis, the probes least likely to be associated with the factor of interest are used as adjustment (the probes in the bottom half of the first list) to capture further components of unwanted variation, including biological factors such as cellular heterogeneity. All P values were adjusted for multiple testing using the Benjamini-Hochberg method, with a false discovery rate threshold of 5% (52). To facilitate biological interpretation, percentage methylation was expressed as β values (β = Meth/[Unmeth + Meth + 100]).

IPA. The identified differentially methylated CpGs between GDM and control offspring were assigned to linked genes using the HM450K array annotation file. These were then subjected to pathway analysis using the IPA tools based on curated databases of previously ascertained biological interactions. This approach, however, does not take into account the bias toward cancer genes and promoter regions on the HM450K array.

ENCODE data analysis. To determine the likely functional relevance of CpG sites of significant differential methylation, we interrogated the ENCODE data (ChromHMM, cell line GM12878, Broad Institute) (53).

Validation and replication. We chose to validate 3 different genes from the discovery study. These 3 were selected by, firstly, prioritizing genes with more than 1 differentially methylated CpG site in the same gene (there were, in total, 2 genes with 2 differentially methylated CpGs) and, secondly, prioritizing P value significance (the top 1 CpG/gene on the list. These genes and CpG sites were as follows: ESM1 (cg00992687 and cg09452568), MS4A3 (cg14328641), and PDE6A (cg09109411). For ESM1, both CpG sites that were found to be differentially methylated in the discovery study were located in close proximity; therefore, it was possible to include both sites in the same pyrosequencing assay. An overview of ESM1, MS4A3, and PDE6A CpG sites investigated are shown in Supplemental Figure 1. We conducted both a technical validation of the array samples (n = 188) and a biological validation/replication, using the remaining samples from the original cohort as a perfectly matched replication cohort (n = 905) (Figure 1).

DNA BS conversion was performed with 500 ng using the EZ-96 Methylation Gold Kit (Zymo). Measurement of DNA methylation was performed using the Pyrosequencing technique (PyroMark Q96ID, Qiagen) with PyroMark Gold Q96 reagents. PCR and sequencing primers were designed using PyroMark Assay Design 2.0 (see Supplemental Table 4 for primer sequences). Data was analyzed with PyroMark Q96 (version 2.5.8) software. After QC, 788 samples passed for the ESM1 CpG assay, 787 samples passed for the MS4A3, and 445 samples passed for the PDE6A assay and were included in the final analyses.

Statistics. Continuous variables were tested for normality. Normally distributed data are presented in mean ± SD and compared by 2-tailed Student’s t test. Irregularly distributed data are presented as median and interquartile range (IQR) and compared by the Mann-Whitney U test. The χ² test was used for analysis of frequencies.

The statistical analysis of the genome-wide methylation data by the RUVm method is described above.

Multivariate linear regression models were applied to address potential confounding effects on the association between GDM status and offspring methylation level of the differentially methylated CpG sites from the discovery study. In order to test whether the association of GDM with offspring methylation level of the differentially methylated CpGs was influenced by the offspring phenotype, we conducted regression analyses by adjusting for offspring sex, age, and BMI in the discovery cohort. In a secondary analysis, we investigated if prepregnancy mBMI had confounding effects on the results of the differentially methylated CpGs, also adjusting for offspring sex and age.

In the replication study of the selected ESM1, MS4A3, and PDE6A CpG sites, multivariate regression modeling was applied to address the association between GDM status, prepregnancy mBMI, and DNA methylation levels, with adjustment for offspring sex and age.

In all regression models, we estimated mean differences as β coefficients and 95% CI. Assumptions of equal variance and normally distributed residuals were visualized in quantile-quantile plots (QQ plots) and histograms. All statistical analyses were performed using the SAS 9.3 statistical software (SAS Institute Inc.) and a P ≤ 0.05 was considered statistically significant.

Study approval. Consent from both parents was essential for the participation of the child in the study. The study design and protocol were approved by the Regional Scientific Ethical Committee of the municipalities of Copenhagen and Frederiksberg (H-4-2011-045 and H-4-2013-129) and conformed to the Declaration of Helsinki.