Cord blood epigenome-wide meta-analysis in six European-based child cohorts identifies signatures linked to rapid weight growth

Background Rapid postnatal growth may result from exposure in utero or early life to adverse conditions and has been associated with diseases later in life and, in particular, with childhood obesity. DNA methylation, interfacing early-life exposures and subsequent diseases, is a possible mechanism underlying early-life programming. Methods Here, a meta-analysis of Illumina HumanMethylation 450K/EPIC-array associations of cord blood DNA methylation at single CpG sites and CpG genomic regions with rapid weight growth at 1 year of age (defined with reference to WHO growth charts) was conducted in six European-based child cohorts (ALSPAC, ENVIRONAGE, Generation XXI, INMA, Piccolipiù, and RHEA, N = 2003). The association of gestational age acceleration (calculated using the Bohlin epigenetic clock) with rapid weight growth was also explored via meta-analysis. Follow-up analyses of identified DNA methylation signals included prediction of rapid weight growth, mediation of the effect of conventional risk factors on rapid weight growth, integration with transcriptomics and metabolomics, association with overweight in childhood (between 4 and 8 years), and comparison with previous findings. Results Forty-seven CpGs were associated with rapid weight growth at suggestive p-value <1e−05 and, among them, three CpGs (cg14459032, cg25953130 annotated to ARID5B, and cg00049440 annotated to KLF9) passed the genome-wide significance level (p-value <1.25e−07). Sixteen differentially methylated regions (DMRs) were identified as associated with rapid weight growth at false discovery rate (FDR)-adjusted/Siddak p-values < 0.01. Gestational age acceleration was associated with decreasing risk of rapid weight growth (p-value = 9.75e−04). Identified DNA methylation signals slightly increased the prediction of rapid weight growth in addition to conventional risk factors. Among the identified signals, three CpGs partially mediated the effect of gestational age on rapid weight growth. Both CpGs (N=3) and DMRs (N=3) were associated with differential expression of transcripts (N=10 and 7, respectively), including long non-coding RNAs. An AURKC DMR was associated with childhood overweight. We observed enrichment of CpGs previously reported associated with birthweight. Conclusions Our findings provide evidence of the association between cord blood DNA methylation and rapid weight growth and suggest links with prenatal exposures and association with childhood obesity providing opportunities for early prevention. Supplementary Information The online version contains supplementary material available at 10.1186/s12916-022-02685-7.


Background
Childhood obesity was declared an epidemic by the World Health Organization (WHO) more than two decades ago [1]. Forty million children below 5 years of age were affected by overweight or obesity in 2016 [2]. The changes in eating behaviors, available food choices, and physical activity observed during the COVID-19 pandemic have further amplified this global health issue [3]. Being obese during childhood is associated with both short-term health consequences (psychosocial, including social stigma, depression, and anxiety [4], and physical, including asthma [5]) and long-term health consequences (including obesity, cardiovascular diseases, diabetes, and cancers in adulthood), which can lead to disability and premature death [6,7]. There is increasing evidence that the path to childhood obesity is established early in life, and rapid weight growth (RWG) in infancy has emerged as a major early risk factor [8][9][10], consisting of an upward centile crossing in weight growth charts and defined as a change greater than 0.67 in weight standard deviation (SD) scores. In utero and early life are critical time windows due to developmental plasticity. During these periods, exposure to detrimental factors may induce long-lasting alterations increasing susceptibility to diseases later in life, as postulated by the developmental origin of health and disease (DOHaD) theory [11]. RWG represents an early phenotype occurring as a thrifty adaptive response compensating the effects of adverse exposures, which in the long term becomes detrimental [12].
Biological mechanisms underlying RWG are poorly understood. Epigenetics, the mitotically inheritable changes in gene function not explained by changes in the DNA sequence, is a possible mechanism through which in utero exposures influence health and disease later in life [13]. Epigenome-wide association studies (EWAS) in large population-based child cohorts coordinated by the Pregnancy And Childhood Epigenetics (PACE) consortium found that several DNA methylation marks (N=914) at birth are associated with birthweight [14], in contrast with only one association with weight in childhood [15]. Previous studies investigating cord blood DNA methylation and early infancy weight and weight growth were limited to candidate genes (IGF2 [16], TAC-STD2 [17], MEG3 [18]) and gestational age acceleration (representing the difference between epigenetic gestational age predicted by CpG targets and actual gestational age determined using the last menstruation and/or ultrasounds) [19], and one previous EWAS was conducted in a small case-control study (N=40 children with RWG versus 40 without) [20]. Furthermore, none investigated concomitant changes spanning entire genomic regions, known as differentially methylated regions (DMRs), nor investigated which DNA methylation marks predict RWG. A better understanding of biological mechanisms acting at birth and underlying RWG is critical to creating effective policy and developing workable early-life prevention programs.
Therefore, in this study, we conducted a meta-analysis of six European-based child cohort (N=2003) EWAS to test the association of cord blood DNA methylation with RWG at 1 year at single CpG sites and CpG genomic regions. Furthermore, we investigated the association of gestational age acceleration with rapid weight growth via a fixed-effect meta-analysis similar to the main analysis. Then, we tested if the DNA methylation marks identified as related to RWG (i) improved the prediction of RWG by assessment of the predictive performance of models incorporating identified DNA methylation levels and conventional risk factors, including maternal education level [21], age at delivery [22], smoking during the index pregnancy [23], pre-pregnancy body mass index (BMI) [24], parity [25], and child gestational age [26]; (ii) were mediators of the effect of conventional risk factors on RWG; (iii) were associated with transcriptome and the metabolome, to guide functional interpretation; (iv) were associated with childhood overweight phenotype (between 4 and 8 years); and (v) overlapped with the previous literature findings.

Study population
The study population includes six European ancestry children cohorts: (i) the Avon Longitudinal Study of Parents And Children (ALSPAC) [27], (ii) the ENVironmental Influences ON early AGEing (ENVIRONAGE) study [28], (iii) the Generation XXI (GXXI) study [29], (iv) the INfancia y Medio Ambiente (INMA) cohort [30], (v) the Piccolipiù cohort [31], and (vi) the Rhea cohort [32,33]. However, the EWAS analyses are performed within four studies: (i) the ALSPAC, (ii) the ENVIRONAGE study, (iii) the GXXI study, and (iv) the EXPOsOMICS study [34], because the latter study is a combination of samples from four cohorts (ENVIRONAGE, INMA, Piccolipiù, and RHEA). EXPOsOMICS samples have been analyzed at the same moment, in the same laboratory, randomized across the different arrays. The ENVIRONAGE cohort is part of EXPOsOMICS and has additionally performed separate DNA methylation arrays; nevertheless, there are no samples in common between the two studies. Multiple births (non-singleton) were excluded from the analyses. Participants of the cohorts included in the main analyses had a combined sample size of 2003 children with completed information (Table 1).

DNA methylation
Cord blood DNA methylation was measured using the Infinium HumanMethylation450 BeadChip in ARIES and EXPOsOMICS [39,40] and using Infinium Meth-ylationEPIC BeadChip in ENVIRONAGE and GXXI studies. Each cohort independently performed quality control and normalization of DNA methylation data, as reported in Additional file 1: Supplementary Methods [27][28][29][30][31][32][33][34][35][36][37][38]. DNA methylation levels were trimmed using the Tukey method if the removal of outliers had not been performed by cohort-specific preprocessing, and were expressed as beta values. CpGs were annotated to the nearest gene by the annotation provided by Illumina.

Cell type estimation
Cell types were estimated through established de-convolution approaches using Gervin's in for ARIES and GENXXI [41] studies and Bakulski's method in ENVI-RONAGE and EXPOsOMICS studies [42].

DNA methylation gestational age estimation and gestational age acceleration
DNA methylation gestational age was estimated using the Bohlin [43] and Knigh's [44] epigenetic clocks via the methylclock R package (version 1.0.0) [45]. We used Spearman's correlation to select the clock that predicted best the chronological gestational age. Gestational age acceleration was calculated as the residuals of the regression of chronological gestational age on DNA methylation gestational age.

Rapid weight growth
In all the cohorts, weight measurements were obtained from obstetric records at birth, while measurements at later times were measured by trained staff (in ALSPAC and GXXI) or measured by trained staff or self-reported from parents (in ENVIRONAGE and EXPOsOMICS) as detailed in Additional file 1: Supplementary Methods [27][28][29][30][31][32][33][34][35][36][37][38]. A two-step prediction approach was used for calculating sex-and age-specific weight at exactly 1 year, using fractional polynomials of age by gender in each cohort, as previously described [46]. Weight gain (WG) at 1 year was calculated as the difference between sex-and age-adjusted WHO-SD scores of birthweight and predicted weight at 1 year. Children were classified as having RWG if WG was > 0.67 SD scores according to Ong et al. [47].

Covariates
In all the models, the following covariates, identified as a priori confounders or potential predictors depending on the performed analysis, were considered: maternal tobacco smoke during pregnancy [23] (categorized as smoker or non-smoker), education level at delivery [21] (categorized as low, medium, and high education based on cohort-specific preferential classification), prepregnancy BMI [24] (in kg/m 2 ), age at delivery [22] (in years), parity [25] (categorized as nulli-and multiparous), and child gestational age [26] (in weeks) based on last menstrual period or ultrasound, sex, and cohort membership (for EXPOsOMICS only). Bead array row and bisulfite conversion batch were considered as technical confounders. All the analyses were adjusted for cell types estimated as described before. Infants with low birthweight are more likely to experience RWG, but birthweight could be either a mediator or a confounder in the association between cord blood DNA methylation and RWG, as it is measured at the same time as DNA methylation. Hence, birthweight was not included as a confounder in the main analyses but only in sensitivity analyses, along with delivery mode [48] (categorized as natural delivery or cesarean section) which could also be a potential confounder. Additional covariates, which may also confound the association under study, were used to restrict analyses and included maternal gestational diabetes [49] (categorized as having or not gestational diabetes) and ethnicity [50] (categorized as white European or non-white European children). A detailed description of the covariates in each cohort is reported in Additional file 1: Supplementary Methods [27][28][29][30][31][32][33][34][35][36][37][38].

Gene expression
Gene expression levels were measured in the 200 cord blood samples of the ENVIRONAGE cohort participating in the EXPOsOMICS project [34]. Total RNA was extracted using the total RNA miRNeasy mini kit (Qiagen, Venlo, Netherlands) according to the manufacturer's protocol, as detailed previously [39,51]. In brief, samples were quality checked and further hybridized onto Agilent Whole Human Genome 8×60 K microarrays, and microarray signals were detected by an Agilent DNA G2505C Microarray Scanner. After preprocessing and quality control were performed using an in-house developed R pipeline, gene expression was log 2 transformed and normalized by quantile normalization using arrayQC. Residuals from linear regression between transcripts and hybridization date and leucocyte were used instead of the actual measures of gene expression to account for technical noise for a final sample size left available for further analysis of 29,164 transcripts for 165 children.

Metabolomics
Untargeted metabolomics was measured in 500 cord blood samples of EXPOsOMICS, as previously described [52,53]. Briefly, a reversed-phase liquid chromatography-quadrupole time-of-flight mass spectrometry system was used in positive ion mode with 499 of the 500 EXPOsOMICS samples successfully analyzed. After data preprocessing, 4712 features for 499 samples were left available for the subsequent analysis. Data were logtransformed and missing values were imputed using the impute. QRILC function within the imputeLCMD R package.

Childhood overweight
Weight and height measurements were available during childhood, between 4 and 8 years of age, in ALSPAC, GXXI, and EXPOsOMICS studies. BMI was determined as the ratio of weight (in kilograms) over squared height (in meters) self-reported by parents or measured by trained staff. When multiple BMI measurements were available, the closest measurement to 6 years of age was considered to assess childhood overweight. Childhood overweight (including obesity) was defined if the child's sex-and age-adjusted WHO-SD BMI score was >2 in children below 5 years of age and >1 in children older than 5 years of age, according to the WHO cut-offs [54].

Statistical analysis
The study workflow is depicted in Fig. 1.

Epigenome-wide association studies of rapid weight growth
First, EWAS were conducted for the single cohorts separately, or in the case of the EXPOsOMICS cohorts together as the separate sample sizes were too small. Logistic mixed models with random effects on bead array row and bisulfite conversion batch were used to test the association between the cord blood CpGs (as explanatory variables) and RWG (as dependent variable) via the lme4 R package (version 1.1.26). Models were adjusted for maternal tobacco smoke during pregnancy, maternal education level at delivery, pre-pregnancy BMI, age at delivery, parity, child gestational age, sex, cohort membership (for EXPOsOMICS only), and cell proportions. Models were adjusted for child sex, although outcomes were based on sex-and age-adjusted SD weight scores, because sex is known to be a major source of variation in DNA methylation. Quality control of the results was performed by visual inspection of scatter plots of coefficients and standard errors, quantile-quantile plots (QQ plots) of p-values, and calculating inflation using the bacon method (λ bacon ) via the R package bacon (version 1.14.0) [55].

Meta-analysis of EWAS of rapid weight growth
Then, single EWAS results were meta-analyzed via fixed-effects meta-analysis weighted by the inverse of the variance using the R package metafor (version 2.4.0). Meta-analysis was performed only on DNA methylation probes available in at least three EWAS (for a total of 398,036 CpG probes). Associations were deemed to be significant at the genome-wide level if p-values were below the Bonferroni-adjustment threshold (p Bonferroni ) of 1.25e−07 (0.05/398,036). We also investigated CpGs at a less stringent p-value threshold (p Suggestive ) of 1e−05. Results were presented as mean (and standard error) differences in log odds of RWG per 10% increase in methylation at each CpG. Quality control of the results was performed by visual inspection of coefficients, standard errors and p-values, and calculation of λ bacon . To assess the genome-wide robustness of the findings to interstudy heterogeneity, we used QQ plots of the p-values for heterogeneity (I 2 ). Probes with I 2 > 50% were excluded from further analyses. Results were inspected by forest plots, and meta-analyses leaving out one cohort at a time were performed. DMRs were identified by using ENmix-comb-p (version 1.22.6) [56] and DMRcate (version 2.0.7) [57] R packages. We used the results from the meta-analysis of EWAS (estimated coefficients, p-values and z-values) as inputs for DMR analyses. The default setting values were used for DMRcate (e.g., minimum number of CpGs in a region = 2 and minimum length of nucleotides = 1000). To correct for multiple comparisons, we used FDR correction in DMRcate and 1-step Siddak correction in ENmix-comb-p. DMRs were considered to be statistically significant if both FDR-adjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p were < 0.01. Sensitivity analyses were also performed by repeating the meta-EWAS, adding birthweight and delivery mode to the confounders, and removing estimated cell counts from the confounders. Furthermore, we restricted meta-EWASs by excluding children born from mothers affected by gestational diabetes and of non-white European ethnicity.

Analysis of DNA methylation gestational age acceleration and rapid weight growth
In each study, we analyzed the association between gestational age acceleration and RWG using logistic models adjusted for maternal tobacco smoke during pregnancy, maternal education level at delivery, pre-pregnancy BMI, age at delivery, parity, child sex, cohort membership (for EXPOsOMICS only), and blood cell estimations. Then, similarly to the analyses described above, results were meta-analyzed via fixed-effects meta-analysis weighted by the inverse of the variance. Results were presented as odds ratios (ORs) (and 95% confidence intervals (95% CI)) of RWG per 1-week increase in gestational age acceleration. We performed sensitivity analyses by adding the mode of delivery and removing estimated cell counts from the confounders. Analyses of DNA methylation gestational age acceleration were further restricted by excluding children born from mothers affected by gestational diabetes and children of non-white European ethnicity.

Prediction of rapid weight growth
We estimated how well RWG was predicted using identified signatures compared to conventional risk factors using the RandomForest R package (version 4.7.1.1). We used three sets of variables: (1) DNA methylation levels of CpGs significantly associated in the meta-analysis of EWAS with RWG at p Suggestive <1e−05, (2) conventional risk factors (including maternal education level [21], age at delivery [22], smoking during the index pregnancy [23], pre-pregnancy BMI [24], parity [25], and child gestational age [26] and sex), (3) DNA methylation levels of CpGs significantly associated in the meta-analysis of EWAS with RWG at p Suggestive <1e−05 in combination with conventional risk factors. Residuals from linear regression between DNA methylation levels and bead array row, bisulfite conversion batch, and cell proportions were used instead of the actual measures of DNA methylation levels to account for technical noise. Missing values were imputed using the missForest package. Data were split into training (including ALSPAC, ENVIRONAGE, and GXXI, being approximately 80% of the total study population) and test (including EXPOsOMICS cohorts for a total of approximately 20% of the total study population) sets. The model was trained on the training set with 10,000 trees. In the test set, the model performance was evaluated using the area under the curve (AUROC) to assess the classifier's goodness of fit. The model calibration was performed by visualizing the agreement between the observed and predicted values [58].
Using the same methodology, we also tested the prediction of RWG by all the CpG sites belonging to each DMR identified as associated with RWG with FDR-adjusted p-values in DMRcate and Siddak p-values in ENmixcomb-p < 0.01.

Mediation analysis
For each CpG associated in the meta-analysis of EWAS with RWG at p Suggestive < 1e−05, we tested in each cohort separately, or in the case of the EXPOsOMICS cohorts together as the separate sample sizes are too small, if DNA methylation was mediator (M) of the effect of conventional risk factors (including maternal education level [21], age at delivery [22], smoking during the index pregnancy [23], pre-pregnancy BMI [24], parity [25], and child gestational age [26]) here named exposures (E) on RWG (Y) via model-based single mediation analysis using the imputation approach [59] via the medflex R package (version 0.6.7). We accounted for technically induced variation by using the residuals from a preliminary linear model of each CpG (as outcome variable) adjusted for bead array row and bisulfite conversion batch, instead of the levels of DNA methylation. Mediation models were adjusted for conventional risk factors other than the exposure of interest, child sex and cohort membership (for EXPOsOMICS only) and blood cell estimations. Child gestational age was excluded from confounders as it could be a mediator of the effect of the other prenatal exposures. Mediation analysis of maternal education was adjusted only for child sex and cohort membership (for EXPOsOMICS only) as the other prenatal exposures may act as mediators of its effect. We estimated the total effect (TE), which was further decomposed into the natural indirect effect (NIE) operating via each mediator M and the natural direct effect (NDE). These effects were meta-analyzed via fixed-effects metaanalysis weighted by the inverse of the variance using the R package metafor and considered significant if p Bonferroni was < 1.14e−03 (0.05/44). Results were reported as ORs (and 95% CI). The TE ORs express the effects of Y on E. The NIE and NDE ORs express the effects of E on Y, mediated and unmediated by M.
We applied the same methodology to explore mediation via the DMRs associated with RWG (with FDRadjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p < 0.01) using all the CpGs belonging to each DMR as joint multiple mediators. Results were considered significant if p Bonferroni was< 3.12e−03 (0.05/16). Finally, we explored mediation via gestational age acceleration and considered significant results with a p-value <0.05.

Transcriptomics and metabolomics functional analysis
To better characterize the functional role of the CpG associated in the meta-analysis of EWAS with RWG at p Suggestive < 1e−05, we integrated CpG measurements with the levels of the transcriptome available for the children of ENVIRONAGE participating in EXPOsOMICS (N = 152) and of the metabolome available in the entire set of EXPOsOMICS participants (N = 444).
CpG methylation levels were regressed against transcript (n= 29,164 transcripts) and metabolite (n = 4712 metabolic features) signals (as outcomes) in linear mixed models with random effects on bead array row and bisulfite conversion batch and adjusted for maternal tobacco smoke during pregnancy, maternal education level at delivery and pre-pregnancy BMI, age at delivery, and parity, gestational age, child sex, cohort membership (for metabolomics analyses), and blood cell estimations. Results were presented as mean (and standard error) differences in transcript and metabolite levels per 10% increase in methylation level at each CpG. Results were considered significant if p Bonferroni was< 3.90e−08 (0.05/(44×29,164)) for transcriptomics analyses and < 2.41e−07 (0.05/(44×4712)) for metabolomics analyses.
Transcriptomics analyses were further restricted to cis transcripts (on the same chromosome up to 10Kb in both directions from CpGs' position).
Overrepresentation analyses (ORA) of the transcripts associated with CpGs in the functional analyses below the p Suggestive of 1e−05 were performed using Consensus-pathDB online tool (http:// conse nsusp athdb. org/), with a pathway considered significantly enriched if the p-value was smaller than 0.05 and included at least three genes.

Analyses of childhood overweight
To assess if the identified CpGs were associated with childhood overweight, we used information on overweight in childhood available in ALSPAC, GXXI, and EXPOsOMICS children (N=1916).
Association between cord blood DNA methylation and childhood overweight was assessed via single studyspecific logistic regression models that were integrated via fixed-effect meta-analysis using the same methodology adopted in the main analysis and adding age at the measurement of BMI in childhood among the confounders. Results were presented as mean (and standard error) differences in log odds of childhood overweight per 10% increase in methylation for each CpG. We performed a look-up of the differentially methylated CpGs in the meta-analysis of EWAS at p Suggestive < 1e−05, which were considered to be statistically significant if p Bonferroni < 1.14e−03 (0.05/44), while DMRs were considered if both FDR-adjusted p-values in DMRcate and Siddak p-values in comb-p were < 0.01.

Comparison with previous findings
We investigated whether CpG sites associated with childhood anthropometrics, as identified by our recent systematic review [60], were associated with RWG by a look-up of these hits (N=1526 CpGs) in our study population. Results were considered significant if p Bonferroni < 3.28e−05 (0.05/1526).
We compared our results at single CpG sites with those of a previous meta-EWAS of birthweight (N=8825 children) and we tested the enriched CpGs significantly overlapped with our analysis using the chi-square test.

Population
Of the total 2003 participants, 34% (677) were classified as having RWG (Table 1). In most of the studies, children born at a shorter gestational age, with lower birthweight, that were first born, and with mothers who smoked during the pregnancy were more likely to show RWG (Additional file 2: Table S1).
Forest plots showing coefficients per each single study did not show evidence of heterogeneity among CpGs with p Bonferroni < 1.25e−07 (Fig. 3). Visualization of I 2 p-values via QQ plot did not reveal heterogeneity (Additional file 3: Fig. S1). Evidence of high betweenstudy heterogeneity (I 2 > 50%) was detected for five (cg04677123, cg22341513, cg22807187 annotated to SNORD115-15, cg24335751 annotated to PRDM16, and cg07780199 annotated to CRCT1) of the 49 CpGs with p Suggestive < 1e−05 (Additional file 2: Table S2 and Additional file 3: Fig. S2), which were removed from subsequent analyses. Analyses leaving out one cohort at a time indicated that no single study had an influential effect on   meta-analysis results ( Fig. 3 and Additional file 3: Fig. S2). Quality control of a single EWAS of RWG revealed the inflation measured by the λ bacon ranged from 0.92 to 1.05 indicating little deflation or inflation from the expected p-values, as detailed in Additional file 3: Fig. S3. In the meta-analysis of EWAS, λ bacon was 1.05 suggesting little evidence of inflation (Fig. 2c).
In sensitivity analyses, we found that the direction of the associations at the 49 CpGs with p Suggestive < 1e−05 was stable, while p-values were attenuated. Adding birthweight as a confounder, the significance of all the associations faded above p Suggestive of 1e−05, except for six CpGs (cg13710814, cg05455036, cg00747477, cg26682904, cg23404711, cg06846833; Additional file 2: Table S3). Adding delivery mode as a covariate, p-values for the association of all the CpGs were still below p Suggestive of 1e−05, with the exception of one CpG (cg16072126), and p-values of all the three genome-wide significant CpGs in the main analyses were still below the p Bonferroni significance level of 1.25e−07 (Additional file 2: Table S3). Removing cell types from confounders, 21 CpGs had still p Suggestive < 1e−05, and the p-values of two of the genome-wide significant CpGs (cg25953130 and cg14459032) were still below the p Bonferroni significance level of 1.25e−07 (Additional file 2: Table S3). Excluding subjects with mothers having had gestational diabetes or unknown information (N = 113) and children that were non-white European or unknown information (N=53) (Additional file 2: Table S4), p-values of the association of 28 and 35 CpGs were still below p Suggestive of 1e−05, respectively, in each analysis, and p-value of one out of the three genome-wide significant CpGs faded above the p Bonferroni threshold of 1.25e−07 in the analyses excluding gestational diabetes cases (cg00049440) and nonwhite European children (cg14459032) (Additional file 2: Table S3). In sensitivity analyses of the DMR analyses, Plots show log odds of showing rapid weight growth per 10% increase of methylation levels and 95% confidence intervals from leave out one cohort at time analyses. 95% CI 95% confidence interval when birthweight was added as a confounder, the significance of all the DMRs faded above the FDR-and Siddak-adjusted p-value threshold of 0.01 in DMRcate and ENmix-comb-p, respectively, except for the THEM5 and AURKC DMRs (Additional file 2: Table S5). When adding delivery mode, all the regions, except for the AURKC DMR had FDR-and Siddak-adjusted p-values below of 0.01 in DMRcate and ENmix-comb-p, respectively (Additional file 2: Table S5). Excluding mothers with gestational diabetes or unknown information and children that were non-white European or unknown information (Additional file 2: Table S4), three (LOXL1-AS1, ACTG1, and ALOX12-AS1) and one (FETUB) DMRs, respectively, faded above the FDR-and Siddak-adjusted p-values of 0.01 in DMRcate and ENmix-comb-p, respectively (Additional file 2: Table S5). The p-value of the associations with rapid growth remained below FDR-and Siddak-adjusted p-value threshold of 0.01 in DMRcate and ENmix-comb-p, respectively, for all the regions, apart from seven (located at FETUB, LOXL1-AS1, PRDM16, C17orf64, STK10 and ALOX12-AS1, and GNMT), in models not adjusted for cell types (Additional file 2: Table S5).

Association between gestational age acceleration and rapid weight growth
In all the studies, few CpGs (ranging from 0 to 8) were missing for the calculation of gestational age acceleration using both clocks (Additional file 2: Table S6). The correlation of Bohlin's clock with chronological gestational age was stronger (correlation coefficient r range= 0.63-0.73) than Knight's (correlation coefficient r range = 0.33-0.55), which was excluded in the subsequent analyses (Additional file 2: Table S6 and Additional file 3: Fig. S4). The meta-analysis found that gestational age acceleration was associated with decreasing risk of showing RWG (OR per one gestational age acceleration week=0.71, 95% CI= 0.60-0.85, p-value=9.75e-04) (Fig. 4). Adding delivery mode, removing estimated cells from confounders, excluding mothers with gestational diabetes or unknown information and children that were non-white European or with unknown information, did not mitigate the results (Additional file 3: Fig. S5).

Rapid weight growth prediction
Using Random Forest classification, the rapid growth prediction model showed a predictive ability of an AUROC value of (  Receiver operating characteristic (ROC) mean value of Random Forest prediction models of rapid weight growth. a ROC curves of models including conventional risk factors, 44 CpGs associated with rapid weight growth at p Suggestive < 1e−05 in the meta-analysis of EWAS, and both as identified by the color legend. b ROC curves of models including conventional risk factors, 96 CpGs belonging to the 16 differentially methylated regions associated with rapid weight growth at FDR-adjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p were < 0.01, and both as identified by the color legend. DMRs differentially methylated regions; conventional risk factors include maternal tobacco smoke during pregnancy, education level at delivery, pre-pregnancy BMI, age, parity, and child sex and gestational age (Additional file 3: Fig. S6) or CpGs belonging to DMRs related to rapid growth in DMRcate and Siddak p-values in ENmix-comb-p (Additional file 3: Fig. S7) in addition to the conventional risk factors showed moderately low classification of RWG.

Mediation analysis
Among the 44 CpGs associated with RWG at p Suggestive < 1e−05, three CpGs (cg20068209, located on TMEM30A; cg25953130, located on ARID5B; and cg26433582, located on TPCN2) partly mediated the effect of gestational age on RWG with p Bonferroni < 1.14e−03 (0.05/44), yet the total effect was mainly explained through other pathways (Fig. 6). While some of the 44 CpGs were highly correlated between themselves, these three CpGs were not (Additional file 3: Fig. S8
No association was found between the 44 CpGs associated with RWG at p Suggestive < 1e−05 and the metabolome in the EXPOsOMICS population Fig. 6 Effects of gestational age on rapid weight growth mediated by three CpGs. The plot represents on the x-axis the point estimates odds ratio (dots) and 95% confidence intervals (bars) of the mediation analysis of the 1-week increase of gestation on rapid weight gain via DNA methylation. Only CpGs with p-values of the natural indirect effects below p Bonferroni threshold of 1.14e−03 (0.05/44 CpGs associated with rapid weight growth at p Suggestive < 1e−05 in the meta-analysis of EWAs) are shown. NDE natural direct effect, NIE natural indirect effect, TE total effect (minimum p-value= 8.36e−06 versus p-value Bonferroni threshold=2.41e−07, Additional file 3: Fig. S10A). The 96 CpGs belonging to the 16 DMRs associated with RWG at FDR-adjusted p-value in DMRcate and Siddak p-values in ENmix-comb-p <0.01 were not associated with the metabolome (minimum p-value= 1.01e−06 versus p-value Bonferroni threshold = 1.10e−07, Additional file 3: Fig. S10B).

Association between DNA methylation and childhood overweight
The study population for the analyses of childhood overweight included 1916 children, with a prevalence of overweight of 23.4% (N=435) at a mean age of 6.5 years (Additional file 2: Table S10). Among the 44 CpGs that were associated with RWG at p Suggestive < 1e−05, none was associated with childhood overweight at p Bonferroni < 1.14e−03 (0.05/44) (cg21844291 had the smallest p-value being equal to 0.05, Additional file 2: Table S11). Of the 16 DMRs associated with RWG at FDR-adjusted p-value in DMRcate and Siddak p-values in ENmixcomb-p <0.01, the AURKC DMR was also associated with childhood overweight (FDR-adjusted p-value in DMRcate=2.64e−13 and Siddak p-values in ENmixcomb-p=1.86e−15) (Additional file 2: Table S12). Gestational age acceleration was not associated with childhood overweight (OR per one gestational age acceleration week = 0.85, 95% CI= 0.70-1.04, p-value=0.12) (Additional file 3: Fig. S11).

Comparison with previous findings
In the look-up of the 1526 CpG sites associated with anthropometrics in children by a previous systematic review [60], no signal was significantly associated with RWG in our meta-analysis (Additional file 3: Fig. S12 [60], minimum p-value for cg00753924 located on RXRA = 8.72e−05 versus p Bonferroni < 3.28e−05 (0.05/1526)). Out of the 44 CpGs we identified as associated with RWG, 19 were among the 2423 CpGs (chi-squared p-value<1e−16) associated with birthweight with p-value<1e−05 (out of 473,864 CpGs) in the previous PACE meta-EWAS [14].

Discussion
In this meta-EWAS of six population-based cohorts, including a total of more than 2000 children, we found evidence of an association of cord blood DNA methylation with RWG during the first year of life at the DMR (N=16) level (using FDR-adjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p thresholds of 0.01, respectively), although only three CpGs reached genome-wide significance at the single CpG level (p Bonferroni threshold of 1.25e−07). Considering a less stringent significance level (p Suggestive threshold of 1e−05), our findings of 44 CpGs were enriched in CpGs previously reported in relation to birthweight. However, six were associated with RWG independently of birthweight, indicating that some epigenetic signatures might only be related to RWG while others are a shared mechanism with birthweight (or a correlate of this measure). Our findings also showed that higher gestational age acceleration was associated with a lower risk of experiencing RWG. A small improvement in predicting RWG was obtained by coupling the top CpGs associated with RWG with conventional risk factors. Translating our findings to the metabolomic and gene expression level indicated that the identified CpGs were associated with differential expression of different genes and long non-coding RNAs, but not with metabolic changes. Exploring a set of prenatal exposures previously associated with RWG [21][22][23][24][25][26], we found that gestational age was inversely associated with RWG and part of this association was mediated by the methylation level at three CpGs. AURKC DMR was associated with RWG and childhood obesity between 4 and 8 years old.
The identified genome-wide significant CpGs (cg14459032, cg25953130, and cg00049440) are here related to RWG for the first time; nevertheless, they have been previously associated with birthweight [14,61,62], which is not surprising given the fact that RWG is calculated based on the difference between weight at 1 year and birthweight. cg14459032 was not mapped to any gene, but enhancing annotation as previously described [63], we found the nearest known gene within 10 MB is PCSK5, encoding for a proprotein convertase subtilisin/kexin type 5, which genome-wide association studies have associated to height [64] and high-density lipoprotein cholesterol [65]. cg25953130 is located on the AT-rich interactive domain-containing protein 5B (ARID5B), which encodes for a transcriptional repressor of beige adipocyte biogenesis leading to a shift towards increasing energy-storing white adipocytes [66]. cg00049440 is located on the Kruppel-like factor 9 (KLF9), which encodes a zinc-finger transcription factor that regulates the adipocyte differentiation by binding to peroxisome proliferator-activated receptor γ2 (PPARγ2) promoter [67] and was associated with BMI in adults in a genome-wide association study [68]. We found that increasing cord blood DNA methylation at these sites was associated with a higher occurrence of RWG, which is consistent with the previous findings of an inverse relationship with birthweight, which in turn is inversely associated with RWG [61]. Since the temporal co-occurrence of birthweight and DNA methylation impedes assessing whether birthweight is a mediator or a confounder in the association between DNA methylation and RWG, our main analyses were not adjusted for birthweight. The involvement of birthweight in these associations was further demonstrated in sensitivity analyses, where adjustment for birthweight revealed that the strength of association decreased for the three CpG sites, and the significance faded at a crude p-value of 0.05 at cg25953130. Similar to single CpGs, only two DMRs were associated with RWG independently of birthweight. Our findings suggest that birthweight mostly has an influence on DNA methylation, rather than being a mediator of its effect on RWG. Studies in adults and children supported this hypothesis [69,70]. Along this line, DNA methylation between 10 and 18 years at the one gene (RPH3AL) was found to be predicted by BMI of 10-year-old children [71]. THEM5 and AURKC DMRs (among the 16 DMRs associated with RWG at FDR-adjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p< 0.01) and six CpGs (among the 44 CpGs associated with RWG at p Suggestive < 1e−05) were associated with RWG independently from birthweight. This finding may indicate that some epigenetic signatures are shared by birthweight and RGW, while others are only related to RWG. Thioesterase superfamily member 5 (THEM5) encodes a thioesterase involved in the cardiolipin remodeling process, which in animal models has a critical role in the development of fatty liver [72]. Aurora kinase C (AURKC) encodes a protein kinase that regulates multiple key steps during the mitotic cell division process and whose expression levels in placentas are reduced in early-onset fetal growth-restricted pregnancies [73].
Prediction of RWG was slightly improved by adding the CpGs associated with RWG with p Suggestive <1e−05 or the CpGs belonging to DMRs associated with RWG (using FDR-adjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p < 0.01) to conventional risk factors (c-index increased from 0.67 in models including only conventional risk factors to 0.70 and 0.69 in models including also the CpGs and the CpGs belonging to DMRs, respectively), which could be due to the CpGs laying on the casual paths linking the prenatal exposures to RWG. Our mediation analysis supports this hypothesis by the finding that three CpGs (g25953130, previously described, cg20068209 located on TMEM30A, and cg26433582 located on TPCN2) mediate the effect of gestational age on RWG, although most of the effect of gestational age on RWG is direct. Methylation levels of these CpGs at birth have been previously associated with gestational age [43] and birthweight [14]. cg25953130 in cord blood was also associated with maternal hypertensive disorders in pregnancy [74]. Our analyses could not exclude the possibility that maternal hypertensive disorders or other unmeasured confounders affect DNA methylation at birth, gestational age, and RWG later in infancy.
We identified 16 DMRs associated with RWG (at FDR-adjusted p-values in DMRcate and Siddak p-values in ENmix-comb-p < 0.01), by applying two independent methods (ENmix-comb-p [56] and DMRcate [57]), compared to three only single CpG sites. The difference in the number of identified signals may be due to the greater statistical power of DMR analyses. However, by identifying DMRs and CpGs, the mechanisms underlying the propensity to RWG involve multiple CpGs. In accordance with this finding, a previous study found that at birth DMRs but not single CpG were predictive of the waist-to-hip ratio of children at 5 years, including the PRMD16 DMR [75], which we also found in our study associated with RWG. PRMD16 encodes for a transcriptional regulatory factor that controls the differentiation of brown adipose tissue, and promotes the transition from white to beige adipose tissue and is closely related to obesity [76]. Methylation at single CpG sites of this gene was previously reported to be associated cross-sectionally to childhood obesity and severe obesity [15,77], and more recently with overweight and obesity in adulthood, both at single CpG sites and methylation haplotypes [78]. However, in our study, we were not able to identify an association of this DMR with childhood overweight. AURKC DMR was the only one we identified significantly associated with both RWG and childhood overweight.
Among the identified DMRs, two, LOXL1-AS1 and ALOX12-AS1 DMR, were mapped to genes putatively encoding long non-coding RNAs. DMRs are generally regarded as regions with functional roles. While most studies assume that DNA methylation influences proximal genes, recent multi-omics studies demonstrated that changes in DNA methylation could be associated with distal expression changes [79,80]. Here, we explored the relationships between identified DNA methylation signals and gene expression at the genome-wide level, finding that three DMRs were associated with seven transcripts at the transcriptome-wide level, three of which were cis transcripts (located in the proximity of the DMRs). In addition, investigating the functional relevance of single CpG sites associated with RWG, we found that two CpGs, including cg20068209 (identified as a possible mediator in mediation analyses), were negatively associated with ten transcripts at the transcriptome-wide level, including several long non-coding RNAs. All these transcripts were located far from the identified CpGs (only one transcript shared the same chromosomal location of the associated CpG), both of which are enhancers. This finding suggests that the regulation of transcription goes beyond classical repression of promoter methylation and suggests a possible and more complex interplay between multiple epigenetic layers (methylation, chromatin, and non-coding RNAs). This hypothesis is consistent with findings from a previous study that miRNA expression levels in the placenta were correlated with methylation levels, which in turn were associated with the offspring's anthropometry at 6 years [81].
Pathway analyses confirmed a functional role of the identified DNA methylation signals by identification of gene expression being involved in several pathways, including type II diabetes mellitus, insulin signaling, growth hormone synthesis, immune system, and cell signaling pathways.
Finally, our results indicate that children with greater gestational age acceleration have a lower risk of RWG (p-value<0.05), which is consistent with previous studies finding that greater gestational age acceleration was associated with higher birthweight [82] and lower weight growth from birth to 10 years (although this last association was not significant) [19]. These findings suggest that gestational age acceleration reflects factors related to gestation, such as growth and development, rather than general processes of aging (as adult age acceleration), and is advantageous in childhood. Maternal smoking and socioeconomic position have been associated with infant's DNA methylation age at birth [23,83] and RWG [21]. Our analysis found that maternal tobacco smoke during pregnancy, maternal age, and primiparity were associated with a higher risk of children showing RWG, but gestational age acceleration did not mediate these associations.
The main strengths of this study are the large sample size, the investigation of regional patterns of DNA methylation in addition to the single CpG sites, the incorporation of multiple omic layers in the analyses to explore downstream functionality of DNA methylation targets, and the mediation analyses that we have undertaken to disentangle the role of DNA methylation in the association between conventional risk factors and RWG. We first studied RWG independently of birthweight, as birthweight is part of the definition of RWG. Nevertheless, we tested all associations also taking birthweight into account and hereby further clarifying the role of birthweight in the epigenetic programming of RWG.
We acknowledge some limitations. We used DNA methylation in cord blood as an accessible collectible tissue for large populations, but we acknowledge that cord blood is a mix of cell types. We used estimated cell counts to account for cell variability using two reference methods specific to the umbilical cord that can be combined across studies [41,42,84]. DNA methylation was measured in the participating studies by the Illumina EPIC and 450K BeadChip arrays. We decided to meta-analyze DNA methylation levels measured in at least three studies which reduced the total number of CpG included in the analysis to less than 400,000, representing a small fraction (~1.5%) of the 28.3 million of total CpG sites in the genome [85]. Furthermore, genome-wide technologies are mainly employed to discover biomarkers that in turn require validation by other testing methods (e.g., pyrosequencing, quantitative methylation-specific polymerase chain reaction). Further validation of the identified signals is warranted. RWG was defined based on a fixed threshold (0.67 SD scores), which is well established in the literature [47]. Future studies analyzing the association of DNA methylation and the continuous increase of RWG are needed to robustly replicate these results and may provide further biological insights. Nevertheless, we performed sensitivity analyses for prenatal exposures (adjusting for delivery mode or excluding mothers affected gestational diabetes and children with non-white European ethnicity) known to affect both DNA methylation and RWG and most of the signals we identified were stable to sensitivity analyses, we cannot exclude the possibility of unmeasured confounding by genetic and other prenatal factors, including paternal factors. Finally, in the mediation analysis of single CpGs, we considered each potential mediator as independent, while genomewide techniques are available for correlated high-dimensional data [86]. However, the three CpGs identified as mediators in our analysis were not correlated between themselves.

Conclusions
In conclusion, our findings show that DNA methylation of regions of DNA, gestational age acceleration, and to a lesser extent single CpGs in cord blood were associated with RWG. The DNA methylation signatures identified showed a slight improvement in the prediction of RWG in addition to conventional risk factors. Furthermore, our mediation analysis indicated that some of the identified DNA methylation signatures were mediating the effect of gestational age on RWG. Finally, throughout the analysis incorporating RWG, we identified the AURKC DMR as predictive of overweight in childhood. Results were enriched in CpGs previously reported associated with birthweight. By increasing knowledge of the molecular mechanisms underlying RWG, our results can contribute to identifying target groups and developing prevention strategies already at birth.
Additional file 2: Table S1. Comparison of characteristics of the ALSPAC, ENVIRONAGE, EXPOsOMICS and GXXI cohorts by rapid weight growth status. Table S2. 49 CpGs associated with rapid weight growth in the meta-analysis of EWASs with P Suggestive <1e-05. Table S3. Results from sensitivity analyses adding birthweight and delivery mode, removing cell types from confounders, excluding mothers with gestational diabetes and non-white European children, for the CpGs with P Suggestive <1e-05 in the meta-analysis of EWASs of rapid weight growth. Table S4. Information available on maternal gestational diabetes, and white European ethnicity in the study population. Table S5. Results from sensitivity analyses of differentially methylated regions adding birthweight and delivery mode, removing cell types from confounders, excluding mothers with gestational diabetes and non-white European children. Results are shown only if the region has FDR-and Siddak adjusted p-value < 0.01 in DMRcate and ENmix-comb-p, respectively, in sensitivity analyses. Table S6. Count of missing CpGs for calculation of gestational age clocks and Spearman's correlation coefficients between DNA methylation and chronological gestational age. Table S7. Results from meta-analyses of total, natural direct and indirect effect of prenatal exposures on rapid weight growth via gestational age acceleration. Table S8. Results from overrepresentation analyses (ORA) of transcripts associated at P Suggestive < 1e-05 with the 44 CpGs related to RWG and the 96 CpGs belonging to the 16 DMRs related to RWG. Table S9. List of CpGs belonging to each of the 16 DMRs associated with RWG with FDR-and Siddak adjusted p-value <0.01 in DMRcate and ENmix-comb-p, respectively. Table S10. Study population characteristics of the four studies (ALSPAC, ENVIRONAGE, EXPOsOMICS and GXXI) and total population included in the meta-analysis of EWAS of childhood overweight. Table S11. Look-up analysis of childhood overweight for the 44 CpGs associated at P Suggestive <1e-05 in the meta-analysis of EWASs of rapid weight growth. Table S12. DMR that is associated with rapid growth and which is also significantly associated with childhood overweight.
Additional file 3: Figures S1-S12 [60]. Fig. S1. QQ-plot of I 2 p-values of the meta-analysis of EWASs of rapid weight growth. Fig. S2. Forrest plots and plots from leave out one cohort at time meta-analysis from EWASs of rapid weight growth for the CpGs with P Suggestive < 1e-05. Fig. S3. Quality control of cohort specific EWAS of rapid weight growth. Coefficients, standard errors and p-values distributions visualized via box plots and QQ-plot. Fig. S4. Cohort specific correlation plots between chronological and DNA methylation gestational age. Fig. S5. Forrest plot of the meta-analysis of gestational age acceleration and rapid weight growth in sensitivity analyses adding (A) delivery mode, (B) removing cell types from confounders, excluding mothers with (C) gestational diabetes and (D) non-white European children. Fig. S6. Calibration plots of Random Forest models of rapid weight growth including (A) conventional risk factors, (B) CpGs related to rapid weight growth and (C) both. Fig. S7. Calibration plots of Random Forest models of rapid weight growth including (A) conventional risk factors, (B) CpGs belonging to DMRs related to rapid weight growth and (C) both. Fig. S8. Heatmap shows Pearson's correlation between methylation levels of the 44 CpGs associated with P Suggestive <1e-05 in the meta-analysis of EWAS of rapid weight growth.
Red boxes indicate the three CpGs identified in mediation analyses. Fig.  S9. Volcano plots of the association between the 44 CpGs associated with rapid weight growth at P Suggestive < 1e-05/the 96 CpGs belonging to the 16 DMRs associated with rapid weight growth at FDR-adjusted p-value in DMRcate and Siddak p-values in ENmix-comb-p <0.01 and the entire transcriptome (A/C) and restricted to cis transcripts (B/D). Red lines represent Bonferroni-significant threshold, and black lines suggestive threshold (p-value=10-e05). For analyses restricted to cis transcripts only Bonferroni-significant threshold is represented. Fig. S10. Volcano plots of the association between the 44 CpGs associated with rapid weight growth at P Suggestive < 1e-05/the 96 CpGs belonging to the 16 DMRs associated with rapid weight growth at FDR-adjusted p-value in DMRcate and Siddak p-values in ENmix-comb-p <0.01 and the entire metabolome (A/B). Red lines represent Bonferroni-significant threshold, and black lines suggestive threshold (p-value=10-e05). Fig. S11. Forrest plot of the meta-analysis of the analysis of gestational age acceleration and childhood overweight. Fig. S12. Volcano plot from the look-up in the study population of the CpG sites associated with child anthropometrics in a previous systematic review by Alfano et al. [60].