Placental epigenetic clocks: estimating gestational age using placental DNA methylation levels

1Gestational age (GA) of the fetus is used to forecast the date of delivery, optimize prenatal care, and monitor the growth and development of the fetus relative to other pregnancies. Short GA at delivery impacts neonatal morbidity and mortality, as well as brain development. Thus, accurate classification of the fetus may help predict neonatal risk. In this regard, the World Health Organization defined extremely preterm (<28 weeks of gestation), very preterm (28-32 weeks of gestation) and moderate or late preterm (32-37 weeks of gestation) birth to reflect the newborn’s developmental stage.

2Traditional methods for estimating GA include early obstetric ultrasound measures or calculations based on the last menstrual period (LMP). The early ultrasound method estimates GA based on the visible fetal size (e.g., crown-rump length during the first trimester or biparietal diameter after the second trimester). The LMP method calculates GA based on the time elapsed since the known first day of the LMP. The early ultrasound method is widely accepted as the gold standard due to its higher accuracy but is not routinely available in low and middle-income countries. More accurate classification of GA at birth may help predict neonatal risk for adverse outcomes and measure GA more accurately than through the assessment of physical and neurological features of the newborn, especially when early ultrasound measures are lacking, or the infant is growth-restricted but not preterm.

3Here, we aim to develop a new molecular estimator of GA based on placental tissue samples that is more accurate than the previous clock. Earlier studies have revealed profound molecular changes in placental chorionic villi, the placental structures that project into maternal decidua and are bathed in maternal blood, during gestation. We focus on placental DNA methylation (DNAm) data, because prior work demonstrated that accurate estimators of chronological age (epigenetic clocks) can be developed based on DNAm levels from a variety of tissues, that one can estimate GA based on DNAm data derived from umbilical cord blood samples, and most pertinently that one can estimate GA based on placental methylation data (Mayne et al. 2017). Our study provides more accurate placenta-based GA estimators (i.e. placental epigenetic clocks) than those developed previously, because we use a substantially larger sample for our training set (more than six times larger than that of Mayne et al. 2017). We aim to develop three different placental epigenetic clocks: 1) a “robust placental clock” (RPC) that is largely unaffected by pregnancy conditions (e.g., preeclampsia, gestational diabetes, and trisomy), 2) a “control placental clock” (CPC), tailor-made for measuring GA in normal pregnancies, and 3) a “refined RPC”, trained for uncomplicated term (GA>36) pregnancies. For the RPC, we purposely included placental samples from a variety of pregnancy complications in the training data (e.g., hypertension or diabetes) as well as congenital abnormalities (e.g., trisomy 13, 18 and 21).

4We downloaded publicly available DNAm data from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/; Table 1) that assessed DNAm levels in placental tissues. Eighteen datasets used the Illumina HumanMethylation 450K BeadChip (450K) platform and one used the more recent Illumina MethylationEPIC BeadChip (EPIC) array. Our analyses focused on the 441,870 autosomal CpG probes that are shared between the two Illumina platforms such that the resulting GA estimators (RPC and CPC) would be applicable to data from both platforms.

5An overview of our analysis is presented in Figure 1. We developed the RPC using several placental DNAm datasets (training n=1,102, Table 1, Figure 1). We regressed GA (dependent variable) on DNAm levels of CpG sites using a penalized regression model (elastic net regression). The elastic net regression model automatically selected 558 CpG sites for the RPC model (Supplementary File 1). Predicted GA is a weighted average of DNAm levels at these 558 CpGs, where the weights are the regression coefficients.

6Figure 2 shows the results of a comparison between the RPC and the placental clock from Mayne et al. (2017) in independent test data (test n=187, Table 1). The predictive accuracy of the placental clocks was quantified using the median absolute error (MAE, defined as the median absolute deviation between predicted GA and observed GA), and the degree of the linear association between predicted GA and observed GA was measured using the Pearson correlation coefficient (r). According to both measures, the RPC (MAE=0.96 weeks; 95% confidence interval (CI) [0.88, 1.19], r=0.99; 95% CI [0.98, 0.99]) outperformed Mayne’s placental clock (MAE=2.63 weeks; 95% CI [2.17, 3.01], r=0.94; 95% CI [0.92, 0.96]). Note that Mayne’s placental clock underestimated GA in two data sets: GSE73375 (green dots) and GSE75196 (blue dots), and overestimated GA in two other data sets: GSE66210 (black) and GSE70453 (red).

7The advantage of the RPC is particularly pronounced in later gestation, e.g., when restricting the analysis to placental samples with GA > 25 weeks, the RPC (MAE=0.89 [0.73, 1.02], r=0.82 [0.76, 0.87]) greatly outperforms Mayne's clock (MAE=2.25 [1.9, 2.63], r=0.61 [0.05,0.71], Figure 2C and 2D).

8As expected by its construction, the RPC predicted GA accurately even in placental samples with adverse pregnancy conditions such as preeclampsia, gestational diabetes, and trisomy 13, 18 or 21 (Supplementary Figure S1). However, Mayne’s placental clock underestimated GA in placental samples from preeclampsia cases and overestimated GA in cases of gestational diabetes and trisomy (Supplementary Figure S1). In case of trisomy, the RPC (MAE=2.26 [1.63, 2.88], r=0.12 [-0.25, 0.46]) was more accurate than Mayne’s clock (MAE=3.99 [3.35, 5.4], r=0.02 [-0.34, 0.39]) but still showed a slight overestimation. The RPC's GA estimate was not associated with fetal sex (Supplementary Figure S2). We could not evaluate the effect of ethnicity because our test data did not include ethnic information except for GSE73375 (n=36, Supplementary Figure S3).

9The training data used in the construction of the RPC employed several different normalization methods: functional normalization (funNorm,), subset-quantiles within arrays (SWAN,) and the normal-exponential out-of-band (noob,) approach. This lack of uniformity in normalization methods in the training data has a statistical advantage: it makes it more likely that the RPC will be robust with respect to different normalization methods. In support of this, we found that the RPC validated in test data that were normalized using various methods: beta-mixture quantile dilation (BMIQ,), quantile normalization (quanNorm,), data-driven separate normalization (dasen,) as detailed in Table 1.

10We trained the CPC on placental samples (training n=963, Table 1) that had been designated as "control" samples. Hence, placental samples with higher GA were probably from relatively normal pregnancies. However, placental samples with lower GA might contain samples that would be considered abnormal (i.e. premature rupture of membranes, spontaneous premature labor) but minimal placental pathology relative to preeclampsia cases. The analysis flow was identical as for the RPC, except for the composition of the training and test sets (Supplementary Figure S4). The elastic net regression model used for the CPC automatically selected 546 CpG sites (Supplementary File 1).

11To assess whether adverse pregnancy conditions influence the epigenetic GA estimate, we applied the CPC to placental samples associated with chromosomal abnormalities (confined placental mosaicism, diandric triploidy, trisomy 13, 16, 18 and 21), neural tube defects (anencephaly and spinal bifida), intrauterine growth restriction, maternal complications (gestational diabetes and preeclampsia), and chorioamnionitis (test, n=326). Interestingly, the CPC accurately predicted the GA of fetuses with the above-mentioned conditions (MAE=1.02, r=0.98, Figure 3A) even though the CPC was constructed using unaffected control samples only.

12To test whether pregnancy conditions are associated with faster/slower epigenetic aging, we used epigenetic measures of GA acceleration that were formally defined as raw residuals resulting from regressing the DNAm GA estimate on observed GA. By definition, this residual-based measure of GA acceleration is not correlated with true GA (r=0). GA acceleration did not significantly deviate from zero for any pregnancy conditions mentioned above (Figure 3B), but we acknowledge the small sample sizes for diandric triploidy (n=3) and trisomy 16 (n=3). When restricting the analysis to placental samples from the first trimester (weeks 1 to 12), we found the CPC’s GA estimates to be slightly inaccurate, which was due to the small training set (only n=7 fetuses with GA < 12 weeks).

13For researchers who are particularly interested in uncomplicated term pregnancies, we also developed a second version of the RPC using placental samples from "uncomplicated term” pregnancies (defined as GA > 36 weeks) without any known pregnancy condition.

14Toward this end, we selected "uncomplicated" term placental samples (n=733) from the training set used for the original RPC. Further, we restricted the penalized regression model analysis to the 558 CpGs that make up the original RPC. The penalized regression model automatically selected 395 CpG sites out of the 558 sites (Supplementary File 1). We find that the "refined" RPC for uncomplicated term pregnancies leads to highly accurate GA estimates (MAE=1.49, r=0.98, Figure 4A) in the RPC’s test set (n=187).

15Using RPC’s test data (n=187), we found that previously published epigenetic clocks derived from cord blood samples or other tissues do not apply to the estimation of GA based on placental samples.

16No significant correlation between GA and predicted DNAm age could be observed for clocks by Hannum (2013), Horvath (2013), Levine (2018), and Horvath (2018) (Supplementary Figure S5). However, the DNAm age estimate is close to zero for Horvath’s pan-tissue clock and the more recently developed Skin & Blood clock. Similarly, GA estimators for cord blood (Bohlin's cord blood clock, Knight’s cord blood clock) failed to accurately predict GA in placental samples (Supplementary Figure S6). Overall, these studies demonstrate that the placenta is quite distinct from other tissues regarding the development and application of DNAm based age estimators.

17Several GEO datasets did not report fetal sex (e.g., GSE70453, GSE73375 and GSE76641) and CpGs present on sex chromosomes. Therefore, we developed a fetal sex-classifier based RPC’s training data (n=1,102) using CpGs that are present on autosomes. Toward this end, we regressed fetal sex (binary outcome) on 441,870 autosomal CpG sites using an elastic net implemented in the glmnet R package. The elastic net automatically selected 220 autosomal CpG sites. The classification accuracy was 100% for the placental test data from GSE75196 (n=24). Interestingly, the placental sex classifier turns out to be highly accurate, when applied to blood-based DNAm data from adults (e.g., an accuracy of 96% in the data from the Framingham Heart Study, n=2,356). As the sex of the fetus is typically identical to the sex of its placenta (except for rare cases of chimerism or sex-chromosome mosaicism), the sex-classifier was used to impute fetal sex in GSE66210, GSE70453, GSE73375 and GSE76641.

18We briefly report the results from an epigenome wide association study (EWAS) of GA to demonstrate the profound effect of GA on placental DNAm levels. To protect against confounding by preeclampsia, we conducted EWAS in two separate strata: first, for placental samples from control pregnancies (n=831); second, for placental samples from pregnancies with preeclampsia (n=70). We combined the summary statistics from the two EWAS using Stouffer's method for meta-analysis. The two EWAS summary statistics presented consistent DNAm-GA correlations across 441,870 autosomal CpGs (Figure 5A).

19Strikingly, 10,827 CpG sites exhibit a genome-wide significant correlation with GA (P<1E-07; Figure 5C, Supplementary File 2). Among these sites, 5,940 were in CpGs islands, 262 were in the north shelf, 1,165 in the north shore, 241 in the south shelf, and 902 in the south shore. The top four genes with the largest number of significant CpG sites were MAD1L1 (17 CpGs), BRD2 (13 CpGs), INPP5A (12 CpGs) and RPTOR (9 CpGs). The top 25 CpG sites and their nearest gene(s) are reported in Table 2.

20The RPC had 36 epigenome-wide significant (P<1E-07) CpG sites, the CPC had 39, and the refined RPC had 32.

21Using the largest placental training set to date (n=1,102), we developed highly robust molecular estimators of GA. The robust placental epigenetic clock (RPC) is expected to perform well, even when applied to cases with adverse fetal outcomes or pregnancy complications. We developed this clock using a placenta-based training set that included several adverse conditions, including chromosomal abnormalities (trisomy and triploidy), neural tube defects (anencephaly and spinal bifida), intrauterine growth restriction, maternal complications (gestational diabetes and preeclampsia), and chorioamnionitis.

22In contrast, the only other published placental clock by Mayne and colleagues was trained on a small training set (n=170). In our independent test set (n=187), Mayne’s clock under/overestimated GA according to pregnancy conditions. (Supplementary Figure S1). These systematic deviations from Mayne's clock might reflect interesting biological effects or technical artifacts (batch effects, normalization methods). Another potential limitation of Mayne’s clock is that the authors limited the eligible CpG sites to the approximately 18,437 autosomal sites on the 27K and 450K bead chips. This might explain why the Mayne’s clock uses only 62 CpG sites, whereas our RPC uses 558.

23To infer biological processes under the 558 and 546 CpG sites, we conducted functional gene enrichment analyses using the Genomic Regions Enrichment of Annotation Tool (GREAT, v.3.0,). However, we did not find any significant biological annotations associated with fetal aging. Elastic net regressions automatically select predictive CpG sites of gestational age (GA), but these CpG sites are not always biologically meaningful.

24Our study had several limitations. First, the "observed" GA used for building these epigenetic clocks were estimated either by early pregnancy ultrasound or the LMP method. Although early pregnancy ultrasound based on fetal growth is the gold standard in a clinical setting, it is susceptible to variations in fetal size and leads to a systematic underestimation of GA in smaller fetuses.

25There is also a concern that some of the training sets might be subject to systematic confounding due to adverse pregnancy conditions, as is the case for preeclampsia (Figure 5B). GA tends to be overestimated for placentas linked to preeclampsia, which is consistent with the associated pathology of advanced villous maturation, as well as previous reports of molecular signs of advanced aging. In this hypothetical example, preeclampsia confounds the association between placental DNAm and GA (Figure 5B,). However, this type of confounding probably does not affect our placental clocks for the following reasons. First, the CPC for control samples and the refined RPC for uncomplicated term samples also accurately predicted GA even in pregnancies with known complications. Second, our EWAS of GA reveals profound associations between GA and DNA methylation levels even after stratifying the analysis by preeclampsia.

26Moreover, it is possible that the RPC and the CPC might not perform well in case of non-live births was because the proportion of non-live births was extremely small amongst the third trimester samples in the training datasets, while unavoidably all first and second trimester samples are non-live births. In addition, it has been suggested that gravidity or parity may change placental physiology (e.g., higher placental weight associated with higher parity) and therefore might modify the relationship between the placental epigenome and GA.

27The clinical application of the RPC might be limited, because obtaining placental samples during pregnancy is highly invasive (e.g., chorionic villus sampling). However, the existence of a predictive placental clock – the RPC – opens the possibility to develop another epigenetic clock based on cell-free fetal DNA (cffDNA). cffDNA is fragmented from placenta trophoblasts, and circulates in maternal blood during pregnancy. If the development of a cffDNA clock is successful, clinicians readily estimate GA simply by collecting and analyzing maternal blood anytime during pregnancy.

28We collected publicly available data from Gene Expression Omnibus (GEO) using the GEOparse Python package (Python 3.6.5: Anaconda, Inc.). Table 1 details each dataset. GSE71678 examined the correlation between placental DNAm and arsenic exposures in the New Hampshire Birth Cohort Study. GSE75248 examined placental DNAm in relation to newborns’ neurobehavioral outcomes. GSE71719 studied the association between DNA hydroxymethylation and gene expression using placental samples. The Robinson laboratory (RL) at the University of British Columbia (Vancouver, BC, Canada) transferred placental DNAm data that are publicly available in the GEO database. GSE100197 and GSE98224 were studies that aimed to find placental DNAm profiles for preeclampsia and intrauterine growth restriction in women recruited at the University of British Columbia Women’s and Children’s Hospital (Vancouver, Canada) and at Mount Sinai Hospital (Toronto, Canada), respectively. GSE108567 investigated batch effects in DNAm micro array data. GSE69502 explored DNAm patterns in multi-tissue samples (placental chorionic villi, kidney, spinal cord, brain, and muscle) from fetuses that were aborted due to neural tube defects. GSE74738 aimed to identify differentially-methylated imprinted regions using a genome-wide approach. GSE115508 compared DNAm patterns in cases of placental inflammation (acute chorioamnionitis) with those in unaffected controls. GSE44667 studied the association between placental DNAm in gene enhancer regions and early-onset preeclampsia. GSE49343 investigated placental DNAm with trisomy and preeclampsia. GSE42409 enhanced probe annotation of Illumina HumanMethylation 450K BeadChip to facilitate biologically meaningful data interpretation. GSE120250 examined the impact of assisted reproductive technology on the placental DNA methylome. GSE70453 conducted epigenome-wide and transcriptome-wide analyses of gestational diabetes. GSE73375 examined DNAm in the preeclamptic placenta in relation to the transforming growth factor beta pathway. GSE75196 studied different DNAm patterns in patients with preeclampsia and unaffected controls. GSE76641 studied the transcriptional and DNAm trajectory of 21 organs during fetal development.

29Either the Illumina Infinium HumanMethylation450 BeadChip or the Infinium MethylationEPIC BeadChip was used to measure DNAm level at each CpG site. The DNAm level (-value) was the ratio of two fluorescence signals (methylated and unmethylated). The minfi R package was used to preprocess all the DNAm datasets except for GS2E115508 and GSE120250 (preprocessed by Illumina’s proprietary software, Genome Studio). The preprocessing methods and probe exclusion criteria differed across studies. For example, Marsit and colleagues, the largest GEO submitter, used the funNorm, whereas Robinson and colleagues mostly used the funNorm or the SWAN (Table 1). Other GEO submitters used the BMIQ, funNorm, quanNorm, dasen, or noob. Most GEO submitters excluded probes on sex-chromosomes, near single nucleotide polymorphisms, with cross-hybridization or with a detection p-value > 0.01.

30We ensured that all samples were included only one time in our training data. Some GEO datasets re-used the same samples or included technical replicates. For example, 154 samples were re-used in GSE100197, GSE108567, GSE69502, GSE74738, GSE44667, GSE49343 and RL data; and 15 technical replicates were found in GSE100197 and RL data. The sample size (N) in Table 1 refers to the counts after excluding the re-used samples and replicates.

31We detected and removed outliers using the following steps: 1) we defined a gold standard DNAm profile as the inter-sample median value. For each CpG, we computed the median beta value across all placental samples. 2) The gold standard was correlated with each placental sample to calculate the Pearson correlation coefficient. 3) Placental samples were excluded if their correlation with the gold standard profile was lower than 0.9. Overall, only four putative outliers were removed from the analysis.

32Missing DNAm levels were imputed with the gold standard DNAm levels. Thus, if the beta value of a CpG was missing, the missing value was imputed with the interpersonal median value across all samples. These imputations were only implemented in the training data.

33We fit a penalized regression model using the “glmnet” R package. GA was regressed on 441,870 CpG sites that are shared between the 450K and the EPIC array. The glmnet mixing parameter alpha was set to 0.5 (specifying elastic net regression), and the shrinkage parameter, lambda resulting in the minimum mean square error, was chosen using 10-fold cross-validation in the training data. The RPC automatically selected 558 CpG sites (lambda=0.0936), the CPC did 546 CpG sites (lambda=0.0892), the refined RPC did 395 CpG sites (lambda=0.0116), and the fetal sex-classifier did 220 CpGs (lambda=0.0073). The number of overlapping CpGs between the RPC and CPC was 199. Supplementary File 1 includes CpG sites and their corresponding coefficients for the RPC, CPC, refined RPC and fetal sex-classifier.

34We used the R function "standardScreeningNumericTrait" from the weighted gene co-expression network analysis R package (WGCNA;) to carry out a robust correlation test (based on the biweight midcorrelation) between each CpG and GA. We conducted two separate EWAS of GA: one in control placental samples (n=831) and the other in placental samples from preeclampsia cases (n=70). We computed biweight midcorrelations between DNAm levels and GA, and the corresponding Z statistics and p-values in each stratum. The Z statistics of the two sets of EWAS were combined using the weighted Stouffer’s method as:, where is the square root of the sample size in the th stratum. The corresponding p-values were computed as. The EWAS was limited on the 411,870 autosomal probes available on both the 450K and the EPIC array platform.

35The coefficient values of the placental clocks and the fetal sex classifier can be found in Supplementary File 1.