Epigenome-Wide Meta-analysis Reveals Associations Between Dietary Glycemic Index and Glycemic Load and DNA Methylation in Children and Adolescents of Different Body Sizes

Over the past decades, the prevalence of obesity has increased worldwide, leading to increased morbidity and mortality (1), mainly due to emergence of obesity-related, noncommunicable diseases (NCDs) (2). There is strong evidence that obesity and obesity-related adverse health consequences have their origins in early life. Early childhood is a vulnerable phase when sustained obesity manifests (3). Aside from the prenatal and early childhood periods, late childhood and adolescence have been proposed as a second window of opportunity for establishing lasting lifestyle behaviors that affect the development of obesity and related NCDs (4). This is supported by findings from a large population-based study indicating that an increase in BMI between 7 years of age and early adulthood predicts an increased risk of type 2 diabetes (T2D) in adulthood (5). Growth during childhood and adolescence is affected by a variety of genetic and environmental factors, especially diet.

Diets high in glycemic index (GI) and glycemic load (GL) are associated with cardiovascular risk factors (6), hyperinsulinemia, and markers of insulin resistance (7) and obesity (8) in children, and with T2D, coronary heart disease, and metabolic syndrome in adults (9). The underlying (molecular) mechanisms remain largely unexplored, although they may help explain the known health benefits of a low GI and/or GL diet. DNA methylation (DNAm) is a key epigenetic mechanism involved in the regulation of gene activity. Accumulating evidence suggests that dietary factors, such as folate and resveratrol, can influence DNAm and may thus contribute to the effects of diet on long-term health (10). Recently, maternal dietary GI and GL in pregnancy have been associated with DNAm in cord blood (11). In that study, the findings differed between mothers with normal weight and those with overweight or obesity (hereafter, overweight/obesity), and most associations were observed in the offspring of mothers with overweight/obesity, potentially indicating modification by body size (11). Consistent with these findings, some of the observed associations between GI and, particularly, GL and cardiometabolic health were more pronounced in participants with an increased BMI (9).

So far, the effects of dietary GI and GL on the epigenomes of children and adolescents remain unknown. We hypothesized that dietary GI and GL are associated with blood DNAm in children and adolescents and, therefore, could serve as markers of dietary exposure. We further hypothesized that BMI subgroup analysis reveals a different epigenetic pattern in participants with normal weight than those with overweight/obesity, and we explored its potential role in the development of NCDs by investigating metabolism-related gene expression in blood and subcutaneous adipose tissue.

Six cohorts that contribute to the Pregnancy and Childhood Epigenetics (PACE) consortium (12) participated in this cross-sectional meta-analysis (N = 1,187) (for cohort-specific information, see Supplementary Material). With the exception of the Raine Study Gen2 (Raine) in Australia, all cohorts were from Europe (Spain: Infancia y Medio Ambiente [INMA]; Germany: The LIFE Child Study [LIFE Child]; Leipzig Atherobesity Childhood Cohort [Atherobesity]; the Study in Teens of the Natural Course of Type 1 Diabetes [TEENDIAB]; and Finland: the Northern Finland Birth Cohort 1986 [NFBC1986]). The analysis only included children and adolescents (age range: 4.5–17.0 years). There were no participants with diabetes in the analysis. To avoid the potential underlying genetic effects, only one participant per family was included, based on the completeness of the data, and by equal completeness, the selection was done randomly. Participants with a GI or GL above or below 5 times the SD of the individual cohort mean were also excluded (NFBC1986: GI, n = 2 and GL, n = 11; LIFE Child: GL, n = 1). Further descriptions of the cohorts are provided in the Supplementary Material. Ethical approval and informed consent were obtained from each cohort according to national and international standards (Supplementary Material).

A food’s GI is defined as the 2-hour incremental area under the curve of blood glucose after consumption of a defined portion, typically 50 g of available carbohydrate, expressed as percentage to the glycemic response of a reference food, usually glucose or white bread (13). Thus, a food’s GI ranges between 0 and 100. The GL is defined as the food-specific GI multiplied by the amount of carbohydrates in that food item (13). We used a harmonized method across cohorts to calculate the dietary GI and GL from food frequency questionnaires (FFQs) or dietary records (for cohort-specific information, including validation of assessment methods, see Supplementary Material), as described previously (11). Briefly, the carbohydrate content and GI (using glucose as the reference) of food items were retrieved from country-specific reference data (European countries: Diabetes, Obesity, and Genes [DIOGENES] [14]; Australia: Atkinson et al. [15]). The food-specific information was linked to individual food items consumed, based on the FFQ or dietary record data. If food items were missing from the country-specific reference data, the information was taken from those of other countries. All food items with carbohydrate contents <1 g/100 g were excluded. The individual dietary GI and GL values were calculated as the weighted average of all foods consumed according to common formulas (11).

We used blood DNAm as the outcome measure, determined with either Infinium HumanMethylation450 (450 K) or MethylationEPIC (EPIC) Bead-Chip arrays (Illumina Inc., San Diego, CA). Cohorts with the 450 K array were INMA, NFBC1986, and Raine; cohorts with the EPIC array were LIFE Child, Atherobesity, and TEENDIAB. Each cohort performed their own quality control and normalization procedures (Supplementary Table 1). For all analyses, we used normalized, untransformed β values, where 0 indicated no methylation and 1 indicated complete (100%) methylation. For better interpretation, effect estimates were thus multiplied by 100 to express the change in percentage of methylation. A common prespecified analysis plan with R script for performing the epigenome-wide association studies (EWAS) was provided to all cohorts. We excluded extreme DNAm values, using the Tukey method, as described previously (11). The cohort-specific analyses between dietary GI or GL and epigenome-wide DNAm were performed using robust linear regression (R package Limma) and were adjusted for age, sex, parental education (cohort-specific definition, see Supplementary Material), smoking, total energy intake, blood cell type estimates, and technical variables. The analyses were additionally adjusted for cohort-specific covariables, if needed (Supplementary Material). Total energy intake was unavailable in NFBC1986 and was not included in these analyses. Blood cell types were predicted using the Houseman method, implemented in the R package minfi (16). Because we considered BMI as a potential moderator, we analyzed the same models in the BMI stratum normal weight, as defined by the BMI SD score (SDS) (−2 ≤ BMI SDS ≤1) and overweight/obese (1 < BMI SDS), based on World Health Organization age- and sex-specific reference charts (17,18). All EWAS results were successfully checked for quality, using the R package QCEWAS (19).

Before running the meta-analyses, we excluded probes on sex chromosomes, known cross-reactive probes, and probes with sequence polymorphism (20,21). All remaining probes were meta-analyzed if they were contained in at least two cohorts. Given that the 450 K and EPIC arrays differ in terms of probe numbers, the sample size for array-specific probes varied. The EPIC array probes covered 91.9% of probes on the 450 K array. To correct for potential deflation or inflation and bias in the individual-cohort EWAS results, we applied the R package bacon with the following modifications: niter = 100,000, nburnin = 25,000 (22) (Supplementary Table 2). We then performed fixed-effects invariance-weighted meta-analyses with the bacon-corrected data using the R package metafor (23).

After each meta-analysis, we used Cochrane Q tests to evaluate the heterogeneity of the model. The cohorts Raine and Atherobesity introduced notable heterogeneity in the GL model (Supplementary Fig. 1), and these cohorts, therefore, were excluded from the GL analysis (remaining n = 693) to reduce the possibility of spurious findings. Because the models of the NFBC1986 cohort did not include the total energy intake, we checked the heterogeneity with and without this cohort. Because the heterogeneity did not change after inclusion of NFBC1986 (Supplementary Fig. 2), we performed all meta-analyses including NFBC1986. We report the I² value to indicate potential heterogeneity. To avoid possible human error, the meta-analyses were independently repeated by another analyst using the software METAL (24), and all results were confirmed. The final meta-analyses models showed no sign of deflation/inflation (λ) and bias (µ) (Supplementary Table 3). We used the false discovery rate (FDR) to correct for multiple testing, and P_FDR values < 0.05 were considered statistically significant. Annotation to genomic locations (reference genome hg19) was done with Illumina’s annotation files (450K: IlluminaHumanMethylation450kanno.ilmn12.hg19; EPIC: IlluminaHumanMethylationEPICanno.ilm10b4.hg19). The results of the meta-analyses have been uploaded to the public EWAS catalog (Zenodo; https://doi.org/10.5281/zenodo.7220349).

We checked whether our findings overlapped with significant DNAm sites in meta-EWAS of related traits: maternal dietary GI and GL (11) or glycemia (25) in pregnancy and cord-blood DNAm, and childhood BMI and blood DNAm (26). We performed look-ups within the EWAS catalog (27) and the Accessible Resource for Integrative Epigenomic Studies (ARIES) methylation quantitative trait loci (mQTL) database of blood from children aged 7 years (450 K array only) (28). For the overlapping mQTLs, we searched whether the genetic variants are linked to the following traits: childhood obesity, obesity, T2D, and cardiovascular disease. We performed mendelian randomization (MR) analysis (two-sample), using the R package TwoSampleMR, for the identified risk variant using the Wald ratio. We used the ARIES mQTL data as the exposure set and genome-wide association studies catalog data of these traits as the outcome set to explore potential causal effects of the mQTL–CpGs on obesity and disease risk (Supplementary Material).

To determine whether our findings were significantly related to tissue-specific gene expression, we queried the Human Early Life Exposome (HELIX) cis-expression quantitative trait methylation (cis-eQTM) database (adjusted for blood cell type) of blood from 832 children (mean age: 8.1 years; 450 K array only) (29), and checked for associations with cis-eQTMs in subcutaneous adipose tissue samples from 223 participants (mean age:10.5 years) in the Leipzig Childhood Adipose Tissue Cohort (ClinicalTrials.gov identifier NCT02208141; EPIC array; Supplementary Material).

We queried the Locus Overlap Analysis (LOLA) Core database to search for enrichments of chromosomal markers, enhancer or repressor regions, and transcription-factor binding sites (30).

This analysis included a total of 1,187 children and adolescents (Table 1). Across the cohorts, the mean dietary GI ranged from 50.8 to 61.9, and GL ranged from 99.6 to 161.6. Similar frequencies of sex and smoking and similar ranges of GI and GL were observed after stratification into the normal weight (n = 801) and overweight/obese (n = 386) categories (Table 1).

In the analysis of all available samples, dietary GI showed no association with blood DNAm after FDR correction. The top five CpGs based on P_nominal are presented in Supplementary Table 4. Among these, the association with the smallest SE was observed between dietary GI and cg01578632 (β = −2.5 × 10⁻²; SE = 5.1 × 10⁻³; I² = 0; P_nominal = 8.7 × 10⁻⁷), annotated to WW domain containing E3 ubiquitin protein ligase 2 (WWP2).

Dietary GL was positively associated with blood DNAm at cg20274553 (β = 2.3 × 10⁻³; SE = 4.2 × 10⁻⁴; I² = 0; P_FDR = 2.2 × 10⁻²) located in intron 13 of the WD repeat domain 27 (WDR27) gene (see Table 2 and Supplementary Table 5 for the top five CpGs based on P_nominal).

In the normal-weight stratum, the DNAm of 85 CpGs was significantly associated with GI and the DNAm of 17 CpGs with GL (Fig. 1A and B, Supplementary Tables 6 and 7). There was considerable heterogeneity (I² > 75) in 36 CpGs of the GI analysis and six CpGs of the GL analysis (Supplementary Figs. 3 and 4).

In the overweight/obese stratum, the DNAm of 136 CpGs was significantly related to GI and the DNAm of 298 CpGs to GL (Fig. 1C and D, Supplementary Tables 8 and 9). Considerable heterogeneity (I² > 75) was observed in 40 CpGs of the GI analysis and 47 CpGs of the GL analysis (Supplementary Figs. 5 and 6).

Among the CpGs with highest heterogeneity, mainly one cohort but not always the same one, had a strong effect. Leave-one-out analysis did not show one distinct cohort driving all the associations (Supplementary Figs. 7–10).

In total, the BMI-stratified analyses identified 536 significant DNAm sites (Supplementary Table 10). Of these, 222 were only available in the cohorts with the EPIC array. We found no overlap of significant CpGs across any model (full and BMI-stratified models). However, approximately 93% of the significant CpGs in the BMI-stratified analyses (n = 500 of 536) had the same effect direction as in the whole-sample meta-analysis, and effect estimates showed a strong correlation between the BMI-stratified and whole-sample results (Spearman r = 0.90; P < 2.2 × 10⁻¹⁶; Supplementary Table 10). When comparing the effect estimates of CpGs significant in one stratum and nonsignificant in the other, 245 of 536 CpGs had different directions of effects (Supplementary Table 10).

Table 2 shows the overlapping annotated genes of the significant CpGs across all models. Like the association with GL in the whole-sample analysis, DNAm at a CpG (cg25687360) located in intron 8 of the WDR27 gene was positively associated with GI in children and adolescents with overweight/obesity (Table 2). However, this association was largely driven by the NFBC1986 cohort, and did not remain significant without NFBC1986 (P_nominal = 7.4 × 10⁻¹). In participants with overweight/obesity, CpGs associated with both GI and GL were annotated to the genes karyopherin β1 (KPNB1) and hyaluronan synthase 3 (HAS3) (Table 2).

Across all significant findings (n = 537 CpGs with FDR < 0.05), we found no overlap with CpGs identified in meta-EWAS of related traits (see Research Design and Methods).

By querying the EWAS catalog using all 537 CpGs, we identified overlaps with 1,121 previously published associations (corresponding to 339 unique CpGs; Supplementary Table 11). Enrichment was found for the traits age, sex, tissue, and Alzheimer disease Braak stage (hypergeometric test, FDR < 0.05). Because of the large overlap with age and sex, we evaluated whether these variables were associated with the identified CpGs in the TEENDIAB cohort. Only one CpG was significantly associated with age, and eight were associated with sex (Supplementary Table 12). Of 1,121 overlapping associations with the EWAS catalog, 77 were attributed to environmental exposures, including alcohol consumption, breastfeeding, and sedentary behavior; 95 were attributed to pregnancy-related exposures (e.g., gestational age, maternal overweight/obesity), and seven associations were attributed to glycemic traits (e.g., maternal gestational diabetes) (Supplementary Table 11). Of the 339 unique CpGs, 77 were related to multiple traits (Supplementary Table 13), for example, cg10807894, which has been associated with maternal overweight/obesity, nitrogen dioxide exposure, maternal BMI, eosinophilia, gestational age, and smoking.

Of all 537 CpGs, 99 showed an overlap with mQTLs (Supplementary Table 14). None of the mQTLs overlapped with genetic variants associated with childhood obesity, obesity, and cardiovascular disease, but one, rs9560114, has been linked to T2D. MR analysis showed that the association between this single nucleotide polymorphism and T2D is putatively caused by DNAm of CpG cg24892433. Higher DNAm at cg24892433 was associated with a decreased risk for T2D (β = −6.8 × 10⁻²; SE = 1.2 × 10⁻²; P = 3.6 × 10⁻⁸; Supplementary Table 15). Interestingly, a higher GL was associated with lower DNAm at this CpG in participants with overweight/obesity (Supplementary Table 9).

For the look-up in the blood eQTM data of the HELIX project (450 K array), we used 315 CpGs from our findings located on the 450 K array. We observed overlaps with 76 eQTMs (corresponding to 39 unique CpGs; Supplementary Table 16). The overlapping eQTMs were annotated to 35 genes, and the top five CpG–RNA transcript pairs corresponding to the most significant identified CpGs in the meta-analyses are shown in Table 3.

For the look-up in subcutaneous adipose tissue eQTMs from the Leipzig Childhood Adipose Tissue Cohort (EPIC array), we used 502 CpGs from our findings available in the EPIC array. We found an overlap with 89 eQTMs (corresponding to 48 unique CpGs) associated with 84 unique RNA transcripts (Supplementary Table 17). Among these, the top five most significant CpGs are shown in Table 3. When we compared the overlapping eQTMs from both tissues, eight CpGs were found in both sources; however, these were related to different genes, apart from PLEK (cg13060970).

By querying the LOLA database using the significant CpGs of each BMI model separately, we found enrichments for various regulatory elements, including transcription factors (e.g., forkhead-box protein A2) and active (e.g., H3K9K14ac) or repressive (e.g., H3K9me3) histone modifications (Supplementary Table 18). Contrary to participants with normal weight, CpG sites identified in participants with overweight/obesity showed enrichments for the metabolism-related transcription factors cAMP-responsive element-binding protein 1 and GA-binding protein (Supplementary Table 18).

This study is the largest EWAS of dietary GI and GL and blood DNAm in children and adolescents to date, to our knowledge. Overall, we found one CpG associated with GL in the total group and 536 hits in the BMI-stratified models, particularly in the participants with overweight/obesity. Among these, many CpGs were significantly related to eQTMs in the blood and subcutaneous adipose tissue of same-age children and were located within genetic regions of various regulatory factors.

The CpG positively associated with GL in the total group (cg20274553) is annotated to WDR27, which is involved in cellular scaffolds for protein–protein interactions. We also found that a nearby CpG (cg25687360) in WDR27 was positively associated with GI in participants with overweight/obesity. Greater WDR27 gene expression was detected in a model of hepatic steatosis and, importantly, was linked to higher DNAm in the intronic region of WDR27 (31). Both of the identified CpGs are also located in introns of WDR27, and higher GI or GL exposure was associated with higher DNAm of these probes. Because there was no overlap with eQTM data, confirmation is needed that these CpGs affect WDR27 expression. In a study with patients with syndromic obesity, duplicated copy number variants of WDR27 were observed (32), potentially indicating a pathogenic role here. Therefore, WDR27 is an interesting candidate in the development of obesity and metabolic disease and may be regulated by intronic DNAm, which, in turn, may be susceptible to dietary exposure rich in simple carbohydrates.

Dietary GI and GL were associated with 536 CpGs in participants with normal weight or overweight/obesity, a pattern also observed previously (11). In participants with overweight/obesity, GI and GL were related to CpGs located at the genes KBNP1 and HAS3. Both genes potentially promote metabolic disorders (33,34). Greater hepatic KBNP1 expression has been found in mice fed a high-fat diet, contributing to enhanced pro-inflammatory cytokine levels and insulin resistance (33). HAS3 has been proposed as a potential therapeutic target for the activation of brown adipose tissue (BAT) thermogenesis. Inhibition of HAS3 in mice improved BAT’s thermogenic capacity and thereby prevented diet-induced obesity (34). However, because the respective CpGs were unavailable in the eQTM data, it remains open whether these CpGs, indeed, affect KBNP1 and HAS3 expression.

To shed more light on CpG–gene relations, we interrogated the overlap with eQTMs from blood and adipose tissues from participants of similar ages—an important aspect because many eQTMs differ between young people and adults (29). We observed 165 CpG–transcript pairs possibly affected by GI or GL exposure, such as FRAT regulator of WNT signaling pathway 1 (FRAT1), which prevents glycogen–synthase kinase 3 (GSK3) from inhibiting mammalian target of rapamycin complex 1 (mTORC1) (35). We found that higher GI intake was associated with higher CpG methylation, which was related to higher gene expression of FRAT1. This may lead to reduction in GSK3 activity, allowing increased mTORC1 action (35). The nutrient-sensor mTORC1 is involved in hepatic lipid metabolism promoting, for example, de novo lipogenesis, and inhibiting lipophagy (36). Moreover, mTORC1 stimulation may lead to hyperinsulinemia by impeding autophagy in pancreatic β-cells (37). Therefore, mTORC1 overactivation may contribute to metabolic disorders. Another gene identified was colony stimulating factor 3 (CSF3), which appears to play a role in obesity and insulin resistance. Patients with morbid obesity had a greater abundance of CSF3 in their adipose tissues (38), and exposure to higher levels of CSF3 induced marked insulin desensitization in human adipocytes and myotubes (39). Pleckstrin (PLEK) appeared in both eQTM look-ups and is involved in platelet activation. Although the knowledge about PLEK in the context of cardiovascular disease is limited, a potential role has been suggested in atherosclerosis (40). FRAT1, CSF3, PLEK, and other eQTMs overlapped with CpGs identified in the overweight/obese models, which may indicate that these effects are present and/or more pronounced in people with higher BMI (9).

Although the role of GI and GL on the pathways discussed above requires further investigation, the results of the MR analysis suggest that in children and adolescents with overweight/obesity, GL affects CpG cg24892433 methylation, which, in turn, is causally related to T2D risk.

We can only speculate about the nature of the observation that BMI stratification revealed a higher number of associations compared with the whole-sample analysis. A potential explanation is that certain selection factors (e.g., genetic and/or environmental factors) are more enriched in the stratified groups, leading to (stronger) associations between dietary GI or GL and DNAm. Indeed, the vast majority of identified CpGs had the same direction of effects in the stratified and the whole-sample analysis, whereas effect sizes in the latter were attenuated.

The major limitations of this study include its modest to low sample size, especially in the BMI-stratified analysis. We followed a consortium-wide effort to maximize the sample sizes, but the studies were partly limited by the availability of data on consumed foods, which are required to fully calculate the dietary indices. Heterogeneity arising from the various dietary resources across the cohorts (e.g., FFQ vs. dietary records, type of FFQ [non- or semi-quantitative], and/or the number of collected food items with carbohydrate content) may have influenced the GI and GL scores; we attempted to minimize this by using a harmonized calculation approach. Misreporting of food consumption is a known concern and may have affected the associations. A larger fraction of our findings showed considerable heterogeneity, which can be expected given the inherent heterogeneity of methods of the dietary assessments, the availability of certain probes across cohorts, the modest to low sample sizes per cohort, and the age range from 4.5 to 17 years, because age can influence DNAm. Last, the possibility remains that other factors not accounted for affected the associations.

The major strengths of the study include the harmonized approach to derive the GI and GL values across cohorts and the use of eQTM data from similar age groups to evaluate the functional properties of DNAm sites in blood and, especially, adipose tissue. Unfortunately, eQTM data for childhood and adolescence are rare and omit key metabolic tissues (e.g., the liver). Because key findings of our analysis involved genes with major roles in the liver, a look-up in hepatic eQTM data would have given more insights into tissue-specific CpG–transcript relations.

In conclusion, this meta-analysis revealed 537 blood DNAm sites associated with dietary GI and GL in young participants—the vast majority after BMI stratification. These may be considered potential markers of dietary response. Although the functional importance of the identified CpGs needs to be further defined, multiple CpGs appear to play regulatory roles in the expression of genes involved in the impairment of metabolism and obesity development (e.g., FRAT1, CSF3, WDR27). More analyses with larger sample sizes and complementary designs are required to support our observations and to explore (reverse) causality of these findings and (tissue-specific) functionality of identified CpG–gene relations.

Acknowledgments. Cohort-specific acknowledgments are stated in the Supplementary Material.

Funding. This work was funded by the Joint Programming Initiative—A Healthy Diet for a Healthy Life under proposal number 655 (PREcisE Project), in Germany; by the German Federal Ministry of Education and Research (grant FKZ 01EA1905); by the Medical Research Council and the Biotechnology and Biological Sciences Research Council (MR/S03658X/1) in the U.K. in Spain, by Instituto de Salud Carlos III (grant PCI2018-093147); in the Netherlands, by ZonMw (grant 529051023); and in France, by the French National Research Agency (grant ANR18-HDHL-0003-05). The funding sources of the contributing cohorts or cohort-independent researchers can be found in the Supplementary Material.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. R.O. and S.H. contributed to the research design and drafted the manuscript, and are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. R.O., R.S., H.H.H.-a., J.R., S.F.-B., H.K., and P.E.M. contributed to the cohort-specific analyses; R.O. and H.A. contributed to the meta-analysis; R.-C.H., M.V., W.K., A.K., S.S., M.-R.J., A.G.Z., and S.H. contributed to data acquisition; all authors contributed to the data interpretation and critical revision of the manuscript for important intellectual content, and gave final approval of the version to be published.