Maternal characteristics during pregnancy and offspring DNAm GA at birth

Figures 2, 3, and 4 show the associations between maternal characteristics during pregnancy with the raw DNAm GA difference and the DNAm GA residual of the offspring. The regression models are adjusted for cellular heterogeneity and population stratification. When based on the raw DNAm GA difference, GAA was associated with a maternal age of above 40 years at delivery, pre-eclampsia in a previous pregnancy, fetal demise in a previous pregnancy, and having three or more of the pre-pregnancy risk factors for pre-eclampsia and intrauterine growth restriction (Fig. 2). GAA was also associated with the presence, onset, and severity of maternal pre-eclampsia in the index pregnancy and maternal treatment with betamethasone in the index pregnancy, particularly if the treatment was started a maximum of 30 days before delivery (Fig. 3). GAD was associated with insulin-treated GDM in a previous pregnancy (Fig. 2).

When based on the DNAm GA residual, GAA was associated with a maternal age of above 40 years at delivery, and GAD with insulin-treated GDM in a previous pregnancy and maternal Sjögren’s syndrome (Fig. 4).

Additional file 2: Table S1 shows the unstandardized regression coefficients and 95% confidence intervals for the associations depicted in Figs. 2, 3, and 4 and for the associations between the other tested maternal characteristics during pregnancy and offspring DNAm GA at birth. Additional file 2: Table S2 shows that all of the significant associations remained significant when additionally adjusted for the birth weight SD score based on Finnish national growth references [23].

Additional file 2: Table S3 shows the associations between maternal characteristics and the offspring’s Horvath epigenetic age at birth.

Offspring characteristics and DNAm GA at birth

GAA, based on the raw DNAm GA difference, was associated with lower birth weight, birth length, ponderal index at birth, birth head circumference, placental weight (Fig. 5), being a lower birth weight for GA (continuous and being small-for-gestational-age, <−2 SD), a lower 1-min Apgar score, and female sex (Fig. 6). All models were adjusted for cellular heterogeneity, population stratification, and additionally for sex in the analyses of the offspring birth anthropometry.

When based on the DNAm GA residual, GAA was associated with a lower 1-min Apgar score and female sex (Fig. 7). Additional file 2: Table S4 shows the unstandardized regression coefficients and 95% Confidence Intervals for the associations depicted in Figs. 5, 6, and 7 and for the associations between the other tested offspring characteristics and offspring DNAm GA at birth.

Additional file 2: Table S5 shows the associations between offspring characteristics and the offspring Horvath epigenetic age at birth.

Strengths

The main strength of our study is the use of a well-characterized prospective, ethnically homogenous cohort with data on pre-pregnancy risk factors of pre-eclampsia and intrauterine growth restriction, pregnancy disorders validated by a clinical jury consisting of two qualified physicians and a nurse, and data on other maternal and neonatal characteristics extracted from both medical records and the Finnish Medical Birth Register (MBR) [35]. Our sample was enriched for women with known risk factors for pre-eclampsia and intrauterine growth restriction. This resulted in greater variation of the pre- and neonatal characteristics, thus increasing the possibility of being able to detect their effects on the DNAm GA predictor. Furthermore, clinical GA estimation was based on an ultrasound scan conducted in all women between 12 + 0 and 13 + 6 weeks + days of gestation. DNAm GA was estimated from fetal umbilical cord blood, and we applied novel methods to account for any contamination of the samples by maternal blood. A number of studies have shown that cellular heterogeneity [36], and genetic variation in the population structure [37, 38], can influence epigenetic profiles. Therefore, we removed any potential effects of cell type heterogeneity using bioinformatics methods and corrected for population structure using principal components derived from genome-wide genotypes.

Limitations

Several strengths of this study are also limitations. The ethnic homogeneity of our sample may preclude generalizations for other ethnic groups. Our inclusion criteria, which resulted in enrichment of pregnancy disorders in the study population and increased statistical power to detect their effects, also limit generalizability of these findings to women without such risk factors. Finally, our findings are limited to one tissue type; therefore, we could not test cross-tissue correlations. It is also important to note that for some prenatal characteristics, e.g., for insulin-treated diabetes in a previous pregnancy, maternal Sjögren’s syndrome, less than ten pairs of women and neonates were present in the risk group. Therefore, although we observed significant associations, they should be interpreted with caution and need to be replicated. Finally, while maternal pre-eclampsia in the index pregnancy, maternal treatment with betamethasone, lower birth weight and length, and lower placental weight remained associated with GAA after Bonferroni correction for multiple testing, this correction may be too conservative as the tested associations were not independent.

Study population

Data were taken from the PREDO study, which is a longitudinal multicenter pregnancy cohort study of Finnish women and their singleton children born alive between 2006 and 2010 [39]. We recruited 1079 pregnant women, of whom 969 had 1 or more, and 110 had none of the known risk factors for pre-eclampsia and intrauterine growth restriction (Table 1). The recruitment took place in arrival order when these women attended the first ultrasound screening at 12 + 0–13 + 6 weeks + days of gestation in 1 of the 10 hospital maternity clinics participating in the study. The cohort profile [39] contains details of the study design, inclusion criteria, and all the data that are available. Additional file 1: Figure S1 shows a flowchart of the 814 mother-offspring pairs with data available for the current study.

Offspring DNA methylation, GA, and DNAm GA at birth

Fetal cord blood samples were collected according to standard procedures. DNA was extracted at the National Institute for Health and Welfare, Helsinki, Finland, and the Department of Medical and Clinical Genetics, University of Helsinki, Finland. Methylation analyses were performed at the Max Planck Institute of Psychiatry in Munich, Germany. DNA was bisulphite-converted using the EZ-96 DNA Methylation kit (Zymo Research, Irvine, CA). Genome-wide methylation status of over 485,000 CpG sites was measured using the Infinium Human Methylation 450 BeadChip (Illumina Inc., San Diego, USA) according to the manufacturer’s protocol. The arrays were scanned using the iScan System (Illumina Inc., San Diego, USA). The quality control (QC) pipeline was set up using the R-package minfi (http://bioconductor.org/packages/release/bioc/html/minfi.html). The samples were excluded if they were duplicates, outliers in the median intensities, and because of sex discrepancy. Furthermore, any probes on chromosome X or Y, cross-hybridizing probes as well as probes containing SNPs, and CpGs with a detection P value > 0.01 in at least 50% of the samples, or maternal blood contamination were excluded. Maternal blood contamination was tested using DNAm data at 10 CpGs independently identified as differentially methylated between cord and adult blood and indicative of maternal blood contamination (paper under review). The samples with DNAm values above the previously identified thresholds at five or more of these CpGs were considered contaminated and removed from all future analyses. The final dataset contained 428,619 CpGs. Additional file 1: Figure S1 shows that of the 876 samples available for these analyses, 51 were excluded. Methylation beta-values were normalized using the funnorm function and incorporating the first ten principal components from the internal control probes. To check for batch effects, principal components were computed on these beta values. Two batches, i.e., slide and well, were significantly associated to the main principal components and were removed iteratively using the combat package.

DNAm GA was calculated using the method published by Knight et al. [21] and was based on the methylation profile of 148 selected CpGs.

We calculated a raw DNAm GA difference by subtracting the chronological GA assessed at the first ultrasound screening conducted at 12 + 0–13 + 6 weeks + days of gestation from the predicted DNAm GA. DNAm GA residual was extracted from a linear regression of predicted DNAm GA on ultrasound-based GA.

Offspring cord blood cell counts at birth

Seven cell types (nucleated red blood cells, granulocytes, monocytes, natural killer cells, B cells, CD4(+)T cells, and CD8(+)T cells) were estimated from cord blood methylation using the method of Bakulski et al. [40] which was also incorporated in the R-package minfi.

Offspring genome-wide genotyping and multi-dimensional scaling analysis

Genotyping was performed on Illumina Human Omni Express Exome Arrays containing 964,193 SNPs. Only markers with a call rate of at least 98%, a minor allele frequency of 1%, and a P value for deviation from Hardy-Weinberg-Equilibrium (P > 1.0 e-06) were kept in the analysis. After QC, 587,290 SNPs were available. Any sample pair with IBD estimates >0.125 was checked for relatedness. For most pairs, high IBD estimates could be explained due to partly African origin. As we corrected for admixture in our analyses, these samples were kept except for one pair. For this pair, the high IBD estimate could not be resolved and the other one of this pair was excluded from further analysis. The samples showing discrepancies between phenotypic and genotypic sex were removed. We also checked for heterozygosity outliers, but found none. In total, we genotyped 996 samples of which 13 were excluded. Of the samples with DNA methylation data, eight were excluded based on the above criteria (Additional file 1: Figure S1).

We performed multi-dimensional scaling (MDS) analysis on the IBS matrix of quality-controlled genotypes. Apart of those samples with African admixture, no outliers were detected. The first two components depict this admixture and were included as covariates in the regression analysis.

Maternal characteristics during pregnancy

Pre-pregnancy risk factors of pre-eclampsia and intrauterine growth restriction were derived from medical records screened for by trained study nurses or clinic personnel at maternity clinics of study hospitals at enrolment into the study. These pre-pregnancy risk factors are listed in Table 1.

Pregnancy disorders, derived from hospital records and further verified by a clinical jury [39, 41], included gestational diabetes, which was defined as fasting, 1- or 2-h plasma glucose during a 75-g oral glucose tolerance test ≥5.1, 10.0 or 8.5 mmol/L, respectively, that emerged or was first identified during pregnancy, and further categorized according to treatment as diet or insulin treated; gestational hypertension, which was defined as systolic blood pressure ≥140 mmHg and/or diastolic blood pressure ≥90 mmHg on ≥2 occasions at least 4 h apart in a woman who was normotensive before 20 weeks of gestation; pre-eclampsia, which was defined as systolic blood pressure ≥140 mmHg and/or diastolic blood pressure ≥90 mmHg on ≥2 occasions at least 4 h apart with proteinuria ≥300 mg/24 h. Pre-eclampsia diagnosis was further divided into early (diagnosis before 34 weeks of gestation) and late pre-eclampsia (diagnosis 34 weeks of gestation or later), and also into severe (blood pressure ≥160 mmHg systolic and/or ≥110 mmHg diastolic and/or proteinuria ≥5 g/24 h) and non-severe pre-eclampsia (blood pressure 140–159.9 mmHg systolic and/or 90–109.9 mmHg diastolic and/or proteinuria 0.3–4.9 g/24 h); chronic hypertension was defined as systolic/diastolic blood pressure ≥140/90 mmHg or antihypertensive medication before 20 weeks of gestation (in 24 out of 135 women chronic hypertension was diagnosed during pregnancy). In addition to pre-pregnancy obesity, data on maternal pre-pregnancy BMI calculated from measured weight and height at the first antenatal clinic visit at 8 + 4 (SD 1 + 3) weeks + days were derived from the Finnish Medical Birth Register (MBR) [35] and data on weight change during pregnancy from medical records.

Information on the antenatal corticosteroid (betamethasone) treatment was derived from hospital records (one woman received half a standard dose, i.e., 12 mg/24 h, 22 women received one standard dose of 24 mg/24 h, and one woman received two standard doses totalling 48 mg/24 h), and timing of exposure was further defined as the number of days before birth (over 30 vs 30 days or less before birth).

Maternal smoking during pregnancy (non-smoker, quit during first trimester, smoked throughout), parity (primiparous vs multiparous) and mode of delivery (vaginal delivery vs caesarian section) were derived from the MBR, and alcohol use (yes vs no) and education level (lower secondary or less, upper secondary, tertiary) were reported at 12 + 0–13 + 6 gestational weeks + days.

Offspring characteristics at birth

Weight (kg), length (cm), head circumference (cm), fetal cord blood venous and arterial pH, and 1-min Apgar score were measured at birth, and the birth ponderal index (kg/m³) was calculated. We further divided birth weight and length into normal and small (≤ − 2SD) for GA using Finnish national growth references [23].

Statistical analysis

We tested the associations between maternal and offspring characteristics with the raw DNAm GA difference, DNAm GA residual, and Horvath epigenetic age with linear regressions. All models were adjusted for cell-type composition and population stratification estimated with two multi-dimensional scaling components based on genome-wide data. Maternal characteristic data were further adjusted for birth weight SD score and neonatal anthropometric data for child’s sex. Unstandardized regression coefficients and 95% confidence intervals (CI) represent effect sizes in weeks (raw DNAm GA difference and Horvath epigenetic age) and SD units with a mean of 0 and SD 1 (DNAm GA residual). Nominal 2-tailed p values are given in the tables and Bonferroni-corrected p value threshold reaching a level of p < 0.05 in footnotes. All statistical analyses were performed using SAS 9.4 (SAS Institute, Inc., Cary, NC, USA).