Metabolic Dynamics and Prediction of Gestational Age and Time to Delivery in Pregnant Women

1Pregnancy is one of the most critical periods for mother and child. It involves a tremendous flow of physiological changes and metabolic adaptations week by week, and even small deviations from the norm might have detrimental consequences at different pregnancy stages. For example, approximately 20% of all pregnancies end in miscarriage (< 20 weeks), and around 10% end in preterm birth (< 37 weeks). The latter is the leading cause of global neonatal morbidity and mortality. Of 200 million annual pregnancies, 300,000 pregnancy- and birth-related maternal deaths and 7 million perinatal deaths occur worldwide. With a better understanding of how pregnancy is regulated, even small improvements in obstetric health care can enhance the well-being of many women and children.

2An accurate estimation of the timing of pregnancy and birth is important for many clinical decisions in obstetrics, including determination of preterm birth and related treatment regimens. The current clinical method of determining the gestational age and due date is based on information about the last menstruation date, which can be imprecise, or ultrasound imaging, which depends on accessibility at early pregnancy. Missing the time window is common even in developed countries: in the United States, approximately 900,000 pregnancies annually do not have a prenatal visit before the second or third trimester.

3The maternal circulatory system connects with the fetal circulatory system through the placenta, carrying bioactive molecules and biomarkers such as steroid hormones, micronutrients, and circulating nucleic acids, whose concentrations alter as gestation progresses. Recent work on cell-free RNA suggests that markers in maternal blood can be used to estimate gestational age, but sequencing can be expensive and time-consuming, and the accuracy, at present, is not ideal. Therefore, a more accurate and cost-effective method for estimating gestational age and delivery time, possibly using blood metabolites, is needed. In addition, current clinical tests often only focus on a few markers, whereas research covering more molecules often examines the profiles at one or a few time points during pregnancy. Thus, a high-resolution landscape of pregnancy-related metabolites during healthy pregnancy and the postpartum period is still poorly understood.

4Here, we use untargeted metabolomics to systematically profile blood metabolites throughout pregnancy with weekly sampling of maternal blood. The study identified a large number of pregnancy-related metabolites and metabolic pathways offering a comprehensive view of the metabolite changes during healthy pregnancy and the postpartum period. Leveraging the high-resolution datasets, we built a metabolic clock that not only predicts gestational age in high accordance with the first-trimester ultrasound, the clinical gold standard, but also recovers personal pregnancy variations undetected by ultrasound but capable of affecting delivery time.

5To capture the highly dynamic pregnancy process at high resolution, we established a multi-year single-center Danish normal pregnancy cohort and a design of high-density blood sampling. Consenting female participants submitted weekly blood draws beginning in week 5 of pregnancy and ending in the postpartum period. A total of 30 women with weekly blood samples were assigned to a discovery (N = 21) and a validation (test set 1, N = 9) cohort (Table 1; Figures 1A and S1A). The samples were analyzed in two separate years. In addition, another separate set of women (N = 8) was included as the secondary validation cohort. These samples were analyzed independently three years apart from the discovery cohort (test set 2) (Table 1).

6We randomized the 784 samples from the first 30 subjects within each cohort (Discovery and test set 1), processed them by using a standardized protocol, and analyzed them by liquid chromatography-mass spectrometry (LC-MS) for untargeted metabolomics across two separate years. After quality control, data filtering, and normalization (see STAR Methods) (N = 30), we identified 9,651 metabolic features across the different samples. Of these, 4,995 features (51.7%) were altered during pregnancy and/or the postpartum period (false discovery rate [FDR] < 0.05), suggesting extensive metabolic changes occur during pregnancy. We examined the data globally with principal component analysis (PCA), in which the samples were distributed on the basis of the first two principal components according to their gestational stages (Figure 1B; Scree plot in Figure S1B and the partial least-squares discriminant analysis [PLS-DA] results in Figure S1C), regardless of individual variation and batches (Figures S1D and S1E). Interestingly, we found that metabolites with uni-directional behaviors dominated the features, and over half of them increased across pregnancy until reaching their peaks immediately before labor (Figures S1F and S1G).

7To understand the potential function of pregnancy-related metabolites, we first annotated metabolic features by using an in-house library and a combined public spectral database (see details in STAR Methods). A total of 952 metabolic features were mapped to 687 compounds, which include plasma metabolites with important functions in humans. We then applied significance analysis for microarrays (SAM) to examine the correlation between the abundance of each compound and the reported gestational age of a woman at blood sampling. Among the 687 annotated compounds, 460 compounds were significantly associated with pregnancy (67.0%; FDR < 0.05, SAM). In addition, 264 compounds were identified with a metabolomics standards initiative (MSI) level 1 or 2 identification, among which 176 compounds (66.7%) were significantly associated with pregnancy, as determined by linear regression with gestational age (FDR < 0.05, SAM).

8Our dense sampling revealed detailed temporal patterns of molecular changes. Among the top 68 metabolites (of the 176) that changed over 50% during whole pregnancy, those that increased (N = 30) included steroid hormones estriol-16-glucuronide, estrone 3-sulfate, and tetrahydrodeoxycorticosterone (THDOC). All three increased more rapidly than the well-known steroids progesterone and 17α-hydoxyprogesterone (FDR < 0.05, SAM) (Figures 1C and 1D). By contrast, the top metabolites that decreased during pregnancy (N = 38) were mostly lipids or lipid-like molecules, such as monoacylglycerides (MGs), lysophosphatidylcholines (LPC or lysoPC), and oleoylcarnitine (Figures 1C and 1E). Hierarchical clustering of the weekly samples on the basis of the top 68 altered metabolites revealed a week-order mostly consistent with the actual progression of gestational age (Table S1; Figures 1D and 1E). Intriguingly, most of these metabolite changes rapidly returned to baseline after childbirth (postpartum) (Figures 1B and 1D and 1E). Together, these results suggest a dramatic and programmed change of human blood metabolites at a system level during pregnancy.

9To detect the functional groups of metabolites that change during pregnancy, we performed correlation analysis on the temporal intensity profiles of the top 68 pregnancy-related compounds mentioned above. In Figure S2, metabolites that were significantly increased or decreased tended to cluster together. Using existing structural and biological information, we first categorized the top changing compounds in pregnancy into seven groups. Interestingly, compounds of the same groups tended to cluster together in the correlation matrix. On the basis of the correlation relationship, we constructed a regularized partial correlation network using all pregnancy-related compounds to explore the potential regulatory relationships (Figures 2A and S2B). The topology of the network indicates that different metabolite groups occupied different positions; dense interactions occurred between both inter- and intra- metabolite groups with the densest interactions between central steroid hormones (Figure 2A). These findings highlight that even though the amount of each compound dynamically changes during pregnancy, a highly coordinated metabolite regulatory network underlies the pregnancy process.

10We next examined the main clusters that were present in the correlation analysis. Three main clusters emerged from the hierarchical clustering of metabolites (Figure S2A), with a steroid cluster (e.g., antrostane-3,17-diol, estriol-16-glucuronide, progesterone, 17α-hydroxyprogesterone, and THDOC) sitting between the large clusters of lipids and non-lipid molecules. Compared with the other steroids in this cluster that slowly but steadily increased throughout pregnancy, estriol-16-glucuronide exhibited a rapid increase before week 24, (Figures 1D and 2B). Nearly all upregulated metabolites positively correlated with this cluster of steroids, whereas all downregulated metabolites negatively correlated with this cluster (Figure S2A). This result suggests that different steroid hormones might regulate global metabolome dynamics during pregnancy.

11Within the lipid cluster, intra-correlation was relatively high. The largest cluster was composed of LysoPCs (Figures 2A and S2A), a class of phospholipids. LysoPCs gradually decreased during pregnancy and increased after childbirth in a pattern that highly correlates with the steroid dehydroepiandrosterone sulfate (DHEA-S) (Figure 2C). LysoPCs are bioactive pro-inflammatory lipids that have been linked with organismal oxidative stress and inflammation. The second-largest cluster of lipids included several free fatty acids that were highly correlated within the cluster (Figures 2A and S2A). Long-chain fatty acids showed intricate dynamics in their amounts revealed by the dense sampling. Hexadecadienoylcarnitine and tetracosahexaenoic acid (THA) decreased at the beginning of pregnancy, followed by waves of increased amounts similar to other fatty acids in the late second and third trimesters (Figure 2D). After childbirth, the amounts of most long-chain fatty acids decreased, except for hexadecadienoylcarnitine (Figure 2D).

12Within the non-lipid cluster, one sub-cluster included five highly correlated metabolites belonging to the same caffeine metabolism pathway (Figures 2A and S2A). All five metabolites were consistently elevated during pregnancy, and caffeine reached a concentration three times higher at the end of pregnancy than at the beginning (Figure 2E). This elevation might be due to a slower caffeine metabolism in pregnant women rather than an increase in coffee intake. Overall, among the 68 top-altered metabolites in pregnancy, functional metabolite groups (e.g., steroids, LysoPCs, fatty acids, and caffeine metabolites) were altered in an orchestrated manner during pregnancy, and individual compounds within each group showed inter-correlation to each other (Figure 2A).

13Next, we longitudinally examined the global pathway changes of all 687 annotated compounds during normal pregnancy. Among the 48 mapped Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, 34 showed significant changes (70.8%, adjusted FDR < 0.05, global test) through metaboanalystR (Figure 3A), suggesting large-scale pathway changes of metabolism in pregnancy. To quantify the pathway activities through gestational age, we calculated the average intensity of metabolites in the pathways (Figure 3B; see STAR Methods). Among the top altered pathways (Figure 3A), steroid hormone biosynthesis showed elevated activity precisely timed to gestation, peaking before the end of pregnancy and then declining sharply shortly after delivery (Figure 3B). Along with the essential roles of steroid hormones in maintaining pregnancy and later inducing parturition, we observed an orchestrated elevation of many components centered on progesterone, including some less well-characterized hormones (Figure S3A). Consistent with known sources of pregnancy metabolites (e.g., hormones), metabolite set enrichment analysis (MSEA) revealed that the adrenal cortex, gonad, and placenta were among the top origins of pregnancy-related metabolites (Figure S3B). The ability to recognize many well-known and less-characterized steroid hormone changes across pregnancy validates our approach.

14In addition to the steroid pathway, we observed a dynamic pattern of metabolite changes with pregnancy in other pathways, such as the arachidonic acid metabolism pathway (Figures 3A, 3B, and S3C). We observed 20-HETE amounts increased until week 34; 20-HETE is potentially linked to the regulation of blood pressure and renal function during pregnancy (Figures S3C and S3D). By contrast, 5-HETE amounts generally decreased during pregnancy, potentially associated with its regulation of the uterus (Figures S3C and S3D). Thus, beyond energy metabolism and hormones, a system-wide reconfiguration of the metabolome occurs as the mother adapts to pregnancy. In addition, based on MSEA analysis, many pregnancy-related metabolites are implicated in human disease states, including obesity and prepartum depression (Figure 3C).

15We next determined whether we can build a metabolic clock based on the high-resolution profile to predict gestational age for individual plasma samples. In the discovery cohort (samples n = 507, subjects N = 21), we applied feature selection (lasso [least absolute shrinkage and selection operator]) with all 9,651 features to build the linear regression model that shows optimal cross-validation performance for predicting a given phenotype in this cohort. We then ran the validation cohort data (test set 1, samples n = 245, subjects N = 9) through the model established in the discovery cohort to measure the independent performance of our model (Figure 4A; see STAR Methods).

16We first tested whether the metabolome change can quantitatively determine the gestational age in normal pregnant women. Feature selection in the discovery cohort yielded a linear model that included 42 metabolic features (Figure S4A; Table S2). In the cross-validation test of 507 samples in the discovery cohort, the metabolic model predicted gestational age in weeks (GAmetabolic) that correlated with gestational age estimated by the first-trimester ultrasound (GAultrasound, in compliance with the clinical standard of care) with a Pearson correlation coefficient (R) of 0.96 (R2 = 0.93, p < 1 X 10−100, root mean squared error [RMSE]= 2.49) (Figure S4B). In the independent-validation cohort, the model yielded a similar R of 0.95 (R2 = 0.91, p < 1 X 10−100, RMSE = 2.76, test set 1) (Figure S4C). This indicates metabolic features can accurately predict the gestational age on the basis of a blood sample from a pregnant woman.

17For potential clinical use, we next tested whether we can use the annotated compounds in blood to predict the gestational age in pregnant women. We performed feature selection in discovery cohort by using the 264 level 1 and level 2 compounds identified in the Human Metabolome Database (HMDB) in the discovery cohort (Table S3). This yielded a linear model including five compounds (Figures S4D and 4D) that together are highly predictive. We first evaluated the performance of the model in a 10-fold cross-validation (CV) test in the discovery cohort, in which samples were distributed into folds by subject instead of by sample to prevent person-specific information cross-over between the training folds and the test fold. In the CV test, the metabolic-clock model produced a result (GAmetabolic) that correlated with the gestational age estimated by the first-trimester ultrasound (GAultrasound) with a Pearson correlation coefficient (R) of 0.92 (R2 = 0.85, p = 8e−222, RMSE = 3.67) (Figure 4B). To avoid the hyperparametric selection bias, we further evaluated the performance of our model in two independent validation cohorts (test set 1 and test set 2). In test set 1, the model yielded an R of 0.89 (R2 = 0.80, p = 8e−93, RMSE = 4.11) (Figure 4C). The model, including four steroids and one lipid (Figure 4D), was further verified in a second independent-validation cohort of eight individuals with R of 0.91 (R2 = 0.83, RMSE = 3.05, samples n = 32, test set 2) (Table 1; Figure 4E). The compound identifications were confirmed by chemical standards (Figures 4F–4H, S4E, and S4F; Table S4; see STAR Methods). We noted that four of these five compounds are among the central steroid cluster forming a dense correlation network with one another (Figure 2A).

18As pregnancy progresses toward term, clinical classifications and decisions often need to be made based on the timing of pregnancy (e.g., < 37 weeks for preterm birth). Babies born before 37 weeks are considered preterm, those born before 20 weeks are considered a miscarriage, and those born before 24 weeks have low survival. Therefore, for clinical action it is important to accurately classify the gestational age by clinical cutoff points at weeks 20–37. As a proof-of-principle, we tested the potential use of the metabolome data to classify the normal pregnancy samples as before or after 20, 24, 28, 32, and 37 gestational weeks (Figure S5A). First, using only samples from the third-trimester (> 28 weeks of gestation), the time window where women were more susceptible to preterm delivery, we determined whether the identified maternal blood metabolites can distinguish the sample gestational age as before or after 37 weeks. Both the discovery and the validation prediction yielded an area under the receiver operating characteristics (AUROC) over or close to 0.90 (Figure S5B; see STAR Methods). Remarkably, the prediction model contained only three metabolites, and the abundance range of each individual metabolite separated the > 37 week samples from the < 37 week samples for all but one to two validation subjects (Figures S5C–S5F). Similarly, using samples across the whole pregnancy, we found that metabolites can also accurately distinguish pregnancy samples before or after other important gestational age cutoffs, such as 20, 24, 28, and 32 gestational weeks (Figures S5A and S5G–S5J).

19Next, we examined the metabolic clock prediction performance in individuals. First, we noted that for most individuals, our model produced predictions consistently aligned with the gestational age estimated by the first-trimester ultrasound (Figures 5A and 5B). In both cross validations for Discovery and test set 1, the prediction deviation (measured by RMSE) in individuals centered around 3 weeks (Figures S6A and S6B). However, in each dataset, there is a small population of individuals with higher prediction deviation. When we examined these individuals (e.g., subjects 1, 2, and 4), we found the predictions were not more randomly scattered than other individuals. Rather, in the majority of them, predictions shifted away from the actual gestational ages in a portion of the pregnancy duration (Figures 5A and 5B), suggesting effects from non-random causes.

20We hypothesized that some of these large prediction deviations might arise from biological causes, particularly from the maternal-fetal interaction. It is reported that the fetoplacental unit secretes hormones in conjunction with fetal growth and development. Indeed, we noted that the average prediction deviation strongly correlated with adjusted infant birth weight (Figure 5C, adjusted for gestational length; see Figure S6C and STAR Methods). Thus, the overall metabolic clock tends to outpace the gestational age estimation determined by first-trimester ultrasound in mothers with a heavier fetus while being delayed in the mothers with a lighter fetus. The finding suggests that fetal growth appears to be one of the inputs read by the metabolic clock.

21In addition, within the 18 women with natural labor onset (i.e., excluding women with induction before labor onset and scheduled cesarean-section), we found that the women whose overall metabolic clock of pregnancy outpaced ultrasound evaluation tended to deliver earlier, whereas a delay in metabolic clock correlated with a delayed time to child delivery compared to ultrasound estimated due date (Figure 5D). Interestingly, five out of six women (83%) with a metabolic-clock-predicted gestational age more than one week later than the ultrasound-estimated gestational age had natural labor onset after their due date (estimated by ultrasound, marked in red in Figure 5D), and four out of five women (80%) with metabolic-clock-predicted gestational age more than one week earlier than the ultrasound-estimated gestational age had natural labor onset before their due date (marked in blue in Figure 5D). These results suggest the metabolic clock of pregnancy with maternal metabolites contains information on the timing of delivery in normal pregnancy.

22We then tested whether the maternal blood metabolites can also predict the timing of a normal delivery event within a defined period (2, 4, and 8 weeks from delivery) approaching the labor events (in the third trimester). We first examined whether metabolites can predict a delivery within 2 weeks (weeks to delivery [WD] < 2w). To predict delivery triggered naturally without outside procedures (such as scheduled cesarean-section), we only included delivery events naturally triggered (subjects N = 18, samples n = 193). With just three metabolites, the metabolome accurately predicted an upcoming delivery event within 2 weeks in both discovery and validation cohorts with AUROC close to 0.9 (Figures 5E–5H, S6D, and S6E; see STAR Methods). Similarly, identified metabolites can also be used to predict the timing of a normal delivery event within 4 and 8 weeks (Figures 5E and S6F–S6I). Intriguingly, the panels of metabolites partially overlapped between the models, whereas the individual metabolites contributed differently to the models (Figure 5E). All of the metabolite markers were identified as steroids, except for phospholipid PE(P-16:0e/0:0), and most of them (three out of five in total) also appeared in the aforementioned metabolic clock for gestational age (Figure 5E; Table S4). These results demonstrate that we can precisely categorize critical pregnancy stages in normal subjects by using a small number of maternal blood metabolites, which can be further validated in larger and independent cohorts.

23In this study, we performed untargeted metabolomics profiling and identified highly dynamic temporal regulation of metabolic changes in human pregnancy: more than half of the measured metabolites and metabolic pathways changed during pregnancy. We were able to detect many of the pregnancy-associated metabolite profiles revealed in previous targeted studies (such as progesterone, 17-hydroxyprogesterone, and the linoleic acid pathway), validating our approach. At the same time, we also noted that a large portion of pregnancy-related metabolites identified in our study was less well-studied. For example, over 95% of the pregnancy-related metabolites identified in our study were not recovered from a targeted metabolic profiling study on pregnancy, demonstrating the power of unbiased and hypothesis-independent profiling. Among the changing metabolites, the major class that increased was steroids, including progesterone, which interacts with the hypothalamic-pituitary-adrenal axis (HPA axis), and estriol-16-glucuronide produced by the placenta. Here, the detailed differences in their temporal profiles were revealed by the weekly sampling design of the study (Figures 1D and 2B). In addition, we discovered less well-studied steroids in pregnancy, such as the neurosteroid THDOC, an allosteric modulator of the GABAA receptor that potentially affects stress and depression in human pregnancy. Intriguingly, many pregnancy-related metabolites that changed, including steroids, quickly returned to the maternal non-pregnant state after childbirth (Figures 1B and 1D and 1E). In addition, we also identified a wide variety of non-steroid hormones whose abundance altered during pregnancy progression.

24These metabolite changes presumably accommodate and/or reflect important maternal biological physiology during pregnancy and fetal growth. For maternal nutrient metabolism, one of the decreased carnitines, oleoylcarnitine (Figure 1C), accumulates during certain metabolic conditions, including fasting. Also, one phosphatidylcholine that functions as a micronutrient, lecithin, increased in pregnancy, suggesting a systematic change in the maternal nutritional status during gestation. Within molecules reflecting pregnancy-related physiological changes, consistent with decreased blood pressure, the antihypertensive molecule 20-HETE of the arachidonic acid metabolism pathway is elevated during pregnancy until the early third trimester, and its synthesis is regulated in a renal-specific manner. This reveals the highly dynamic temporal regulation of 20-HETE in blood pressure and kidney function during pregnancy. In contrast, compared with early pregnancy and postpartum, the amount of 5-HETE in the same pathway was generally lower in the second and third trimesters with an increasing trend right before the childbirth, consistent with previous findings that 5-HETE elevates in the uterus and amniotic fluid at the onset of human labor. In the developing fetus, changes in hexadecadienoylcarnitine amounts are associated with congenital heart defects. Here, we revealed that the amount of hexadecadienoylcarnitine in the blood decreased continuously until week 24, then steadily increased thereafter (Figure 2D). In addition, the amount of long-chain fatty acids in maternal blood samples is associated with childhood metabolic health. Here, the omega-3 fatty acid THA decreased during early pregnancy and gradually increased before childbirth (Figure 2D), suggesting gestation-related changes in the formation of docosahexaenoic acid (DHA). Our findings are robust even without a requirement for prior fasting. It will be interesting to validate these findings in cohorts that have dietary information and detailed clinical measurements, to define critical “nutritional time-zones” for micronutrient amounts and further understand the metabolite changes that are important for physiologic changes across pregnancy.

25Our high-density sampling scheme allowed us to study the temporal alteration of metabolite levels at weekly resolution. For example, even though many steroid metabolites were elevated during pregnancy, our profiling was able to show that there were at least two different behaviors: an early wave (such as progesterone and 17α-hydroxyprogesterone) and a second wave (such as estriol-16-glucuronide). These temporal changes of steroids across pregnancy and after childbirth are at least partially regulated by the fetoplacental unit, including both maternal adrenal gland and placenta and fetal adrenals and liver. Further investigation into the interaction of fetal-maternal contribution will be necessary for understanding the temporal regulation of these metabolites.

26Untargeted metabolome and high-density sampling enabled us to identify a broad set of high-resolution temporal profiles of metabolites during pregnancy. We hypothesized that this information might help us to understand the underlying metabolic clock that times the progression of pregnancy. We found that solely using the abundance of five compounds, without any other inputs from clinical features, we can precisely determine the gestational age of a healthy pregnant woman. The precision surpasses the recent cell-free RNA model by using maternal blood. Similarly, with two to three compounds, we can categorically predict many pregnancy cutoff times with high AUC: we can determine whether a woman has reached 20, 24, 28, 32, or 37 weeks (clinical cutoffs for miscarriage, age of viability, extremely preterm, very preterm, and prematurity, respectively) into her pregnancy (Figure S5A), or whether a woman will enter into labor within the next two, four, or eight weeks (Figure 5E). The proof-of-principle study suggested that metabolome bears rich quantitative information about pregnancy progression. However, our study has its limitations. The studied population consisted of healthy Caucasian pregnant women with small variations in clinical characteristics. In the future, we need to test the models in a larger cohort with diverse ethnicities and complications. Meanwhile, targeted chemical assays need to be developed on the small panels of identified metabolite markers that were discovered by untargeted metabolomics to measure the metabolite concentration independent of batches. Intriguingly, we found the metabolic clock of pregnancy to be robust in general, but small personal deviations can be observed, most likely affected by the fetal growth (Figure 5C). Lastly, we also found that the discrepancies between metabolic timing and ultrasound suggested biological significance: the women that had advanced metabolic clock tended to deliver earlier than predicted by ultrasound, whereas a delay in metabolic clock correlated with a delayed time to child delivery (Figure 5D).

27In summary, combining untargeted metabolome and high-density sampling revealed the landscape of metabolome changes during pregnancy and the postpartum period with high resolution. The data itself can serve as a resource for future research. As a proof-of-principle, we also demonstrated that the temporal abundance information of metabolome can be used to predict gestational age with high accuracy in a cohort of healthy women. There is a great need for accurate timing of pregnancy: in the US alone, 900,000 women annually missed their first-trimester ultrasound, currently the only accurate timing method for pregnancy. In low and middle-income countries, accessibility to ultrasound is even more scarce, complicating many pregnancies and fetal care down-stream (e.g., identify imminent labor, manage complications, etc.). Our study demonstrated that the development of clinical tools with a few metabolites in maternal blood to time pregnancy is promising. Testing of blood drawn from the pregnant woman would likely be limited to once or a few times to be informative and have the potential to benefit pregnant women in both developed and developing worlds.

28Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Mike Snyder (mpsnyder@stanford.edu).

29This study did not generate new unique reagents.

30Original data have been deposited to the NIH Common Fund’s National Metabolomics Data Repository (NMDR) website (supported by NIH grant U2C-DK119886), the Metabolomics Workbench https://www.metabolomicsworkbench.org, Project ID PR000918. https://doi.org/10.21228/M81H58.

31The code for the MS/MS identification pipeline used is available on github at https://jaspershen.github.io/metID/index.html https://github.com/jaspershen/metID.

32We recruited pregnant women through family doctors and advertisements (Danish IRB number H-3-2014-004). At enrollment, all women were screened to ensure that they were healthy at baseline, without chronic conditions, and without medication intake of any kind (ages 23 to 36 at giving birth). From each woman, non-fasting blood samples were collected weekly during pregnancy and one sample was collected after pregnancy (2x9 mL EDTA tube and 1xPaxGene RNA tube).

33Within Discovery and Test Set 1 cohorts, 784 normal pregnancy samples from 30 women were completely randomized and analyzed in 12 batches across two years. The 32 normal pregnancy samples in Test Set 2 were randomized and analyzed three years later. Plasma was prepared from whole blood treated with anti-clot EDTA, aliquoted, and stored at −80°C. Plasma (200 μL) was treated with four volumes (800 μL) of an acetone:acetonitrile:methanol (1:1:1, v/v) solvent mixture with internal standards (i.e., Acetyl-d3-carnitine, Phenylalanine-3,3-d2, Tiapride, Trazodone, Reserpine, Phytosphingosine, and Chlorpromazine), mixed for 15 min at 4°C, and incubated at −20°C for 2 h to allow for protein precipitation. The supernatant was collected after centrifugation at 10,000 rpm for 10 min at 4°C and evaporated under nitrogen to dryness (Biotage Turbovap). The dry extracts were reconstituted with 200 μL 1:1 methanol:water before analysis. A quality control (QC) sample was generated by pooling all the plasma samples from 10 women and injected between every 10–15 sample injections to monitor the consistency of the retention time and the signal intensity. The QC sample was also diluted two, four, and eight times to determine the linear-dilution effect of metabolic features.

34MS-grade water (7732-18-5), methanol (A456-500), acetonitrile (A9554), and acetone (67-64-1) were purchased from Fischer Scientific (Morris Plains, NJ, USA). MS-grade acetic acid (64-19-7) was purchased from Sigma Aldrich (St. Louis, MO, USA). Analytical grade chemical standards were purchased [Progesterone (Sigma-Aldrich, P-069-1ML), THDOC (Sigma-Aldrich, P2016-5MG), Estriol-16-Glucuronide (Sigma-Aldrich, E1877-10MG), DHEA-S (Sigma-Aldrich, Dehydroepiandrosterone-D5-3-sulfate (DHEAS-D5) (2,2,3,4,4,-D5) sodium salt solution, D-066-1ML), PE(P-16:0e/0:0) (Avanti Polar lipids, 852470), Androstane-3,17-diol (Sigma-Aldrich, A7755-100MG), 17α-Hydroxyprogesterone (Sigma-Aldrich, H-085-1ML)] and prepared in methanol, except PE(P-16:0e/0:0), which was prepared in chloroform/methanol (8:2).

35Metabolic extracts were analyzed by reversed-phase liquid chromatographic (RPLC)-mass spectrometry (MS) in both positive and negative ionization modes. Thermo Q Exactive Hybrid Quadrupole-Orbitrap plus and Q Exactive mass spectrometers (Xcalibur, Thermo Scientific, San Jose, CA, USA) were operated in full MS-scan mode for data acquisition (acquisition from m/z 500 to 2,000) with a scan rate of approximately 4 Hz and a resolution set at 30,000 (at m/z 400). The MS/MS spectra of the QC sample were acquired under different fragmentation energy (25 NCE and 50 NCE) of the top 10 parent ions. The resulting mass spectra were exported into Progenesis QI Software (Nonlinear Dynamics, Durham, NC, USA) for further processing.

36RPLC separation was performed using Zorbax SB columns (2.1 X 50mm, 1.8 Micron, 600 Bar; 827700-914) purchased from Agilent Technologies (Santa Clara, CA, USA). Mobile phases for RPLC consisted of 0.06% acetic acid in water (phase A) and 0.06% acetic acid in MeOH (phase B). Metabolites were eluted from the column at a flow rate of 0.6 mL/min, leading to a backpressure of 220–280 bar at 99% phase A. A linear 1%–80% phase B gradient was applied over 9–10 min. The oven temperature was set to 60°C, and the sample injection volume was 5 μL.

37Metabolomic features were extracted with a unique mass/charge ratio and retention time, then aligned and quantified with the Progenesis QI software (Nonlinear Dynamics, Durham, NC, USA, http://www.nonlinear.com/progenesis/qi/). Peak deconvolution was performed under default settings in Progenesis QI. Acquired data were processed using an analysis pipeline written in R (https://www.R-project.org). Progenesis QI output was then processed by removing all metabolites that were quantified in less than 30% of the samples or had a median intensity of less than twofold signal over the noise threshold (S/N < 2). The noise threshold was estimated by using the median signal across all the blank runs (if no quantitation was reported in any of the blank runs, the feature was also included in the analysis, as it likely had good S/N characteristics). Then the data were log-transformed and normalized. For each run, the median of all features was centered to correct for variation in the sample amount. Then for each analyte, a linear correction was applied per batch to correct for any linear decrease or increase in abundance during the acquisition of a batch. In short, for each analyte and each batch, a linear model was fitted with the log-abundance of the analyte as the dependent variable and the acquisition number [run order (randomized)] as the independent variable. The model prediction was interpreted as an underlying drift in mass spectrometric sensitivity and subtracted from the analyte level to yield within-batch normalized abundances. Finally, for each analyte, the abundances were median centered by batch to correct for sensitivity differences between batches. The positive- and negative-mode features were then concatenated for downstream analysis. In total, 9,651 features were included in the final analysis. In addition, for samples with more than 50% of the values missing, the sample was removed (one sample in total). The remaining missing values were imputed by the nearest 10 neighbors using the k-Nearest Neighbor algorithm. Note that Discovery and Test Set 1 were normalized together, while samples of Test Set 2 were normalized independently.

38We applied principal component analysis (PCA) to examine the overall distribution of the sample data (with all 9,651 features) and check the run quality. The gestational ages (based on first-trimester ultrasound measurements) were superimposed to facilitate the analysis. During the analysis, the vast majority of the samples were separated by pre- and postpartum in PCA space defined by two components, which explained the largest variations (PC1 and 2, Figure 1B), while two samples of a same subject (last two in her collection, before and after childbirth) displayed irregular behavior in PCA and unsupervised clustering analysis. The two samples were treated as outliers and excluded from further analysis. We also performed partial least-squares discriminant analysis (PLS-DA) according to the categories of gestational age (by the mixOmics package).

39Metabolite identification was performed using a two-step approach. First, to identify compounds, we used our in-house metabolite library, which contains chemical standards and a manually curated compound list based on accurate mass (m/z, ± 5 ppm), retention time and spectral patterns. Second, further metabolites were identified based on accurate mass, isotope pattern and MS/MS spectra against public databases, including HMDB, MoNA, MassBank, METLIN, and NIST.

40Specifically, tandem mass spectrometry (MS/MS) data of QC samples were acquired using a Thermo Q Exactive plus mass spectrometers. The raw MS data (.raw format) were converted to.mgf format files using ProteoWizard (Version 3.0.19095-938eda31a, http://proteowizard.sourceforge.net). Using the metabolic features table (from Waters Progenesis QI) and QC MS/MS data (.mgf format), the metabolic features and MS/MS spectra were matched according to their accurate masses (±25 ppm), and RT values (±30 s). If one metabolic feature matched multiple MS/MS spectra, then all matched MS/MS spectra were used for the identification.

41Next, the generated MS1/MS2 pairs were automatically searched in the public databases: HMDB (http://www.hmdb.ca/), MoNA (http://mona.fiehnlab.ucdavis.edu/), and MassBank (http://www.massbank.jp/). The MS/MS spectra similarity score was calculated using the forward dot-product algorithm, which considers both fragments and intensities. The similarity score cutoff was set as 0.5.

42Furthermore, the metabolic features with MS/MS spectra and not matched in download public databases were searched in the online public databases, METLIN (https://metlin.scripps.edu) and NIST (https://www.nist.gov/). Then the MS/MS spectra match was manually checked to confirm the identifications, which was considered a level 2 identification according to MSI. In addition, the metabolic peaks with MS/MS spectra that were not matched in public databases were analyzed by MetDNA and given a MSI level 4 identification.

43Finally, predictors from the machine-learning models were further confirmed with chemical standards by matching the accurate masses (±5 ppm), retention time (±30 s), and the MS/MS spectra for a MSI level 1 identification.

44In the rare cases, when a given metabolic feature was matched differently between different matching methods, we choose the matching based on the identification level: standards > MS/MS > MetDNA.

45A statistical method specialized for multi-testing, SAM (Significance Analysis of Microarrays) was applied to identify metabolic features/compounds altered significantly in metabolome-wide analysis. Specifically, we used SAM to examine the correlation between abundance of each compounds and the gestational age of each sample in Discovery and Test Set 1 cohorts. For all SAM analyses, distribution-independent ranking tests (based on the Wilcoxon test) and the sample-wise permutation (default by the samr package) were used to ascertain significance (false discovery rate, FDR < 0.05). The adjusted gestational ages were included in a number of plots to present the changes in metabolites among individuals, which were calculated by scaling all delivery event timing to 40 weeks. The populational baseline was calculated by taking the mean intensity values of all women with samples before 14 (20 out of 30 women).

46To identify top changed compounds with abundance increases or decreases more than 50% during the whole pregnancy (40 weeks), we performed a linear regression between log2 abundance and the gestational weeks of samples, and only those compounds with absolute slope larger than log2(1.5)/40 weeks = 0.015 were chosen.

47The regularized partial correlation network captures the remaining association between two nodes after controlling all other information (indirect correlations) in the network. Namely, each node represents a compound, and each edge represents the strength of partial correlation between two nodes after conditioning on all other variables in the datasets. Edge weights represent the partial correlation coefficeients. Lasso (least absolute shrinkage and selection operator) was used to shrink small association coefficient to zero and thus limit spurious correlations in the network. To perform the lasso-based regularized partial correlation, we used qgraph package in R. The tuning parameter gamma(γ), which controls the complexity of the network, was set to 0.5 as suggested. Three measures, strength (the sum of absolute edge weight connected to each metabolite), closeness (inverse of the sum of distances from one metabolite to all others), and betweenness (how often one metabolite is in the shortest paths between other metabolites), indicated how important metabolites are in the network.

48The compounds identified by the methods mentioned above were pooled together. We utilized MetaboAnalystR (https://github.com/xia-lab/MetaboAnalystR) to perform the metabolite set enrichment analysis (MSEA) as well as metabolic pathway analysis (MetPA) on all identified metabolites. For the potential location/organ analysis on metabolites, we excluded male organ/cell types for MSEA. To quantify pathway activity, we averaged the intensities of all identified metabolites for each pathway that includes no less than three identified metabolites and plotted them on the heatmap (Figure 3B). The pathway activity before 14 weeks were averaged across all available samples and subtracted from all later time points. The statistical significance of the changes in a pathway’s activity across pregnancy was evaluated by global testing, the default method used by MetaboAnalystR. The topological pathway impacts were quantified using published method, with MetaboAnalystR. Human desease states that correlated with pregnancy-related metabolites were calculated based on published metabolomics data.

49Three cohorts of data collected and run at different years but from the same center were used to establish Discovery (subjects N = 21, samples n = 507), Test Set 1 (subjects N = 9, samples n = 245), and Test Set 2 (subjects N = 8, samples n = 32) datasets, excluding non-pregnant (postpartum) samples. We applied lasso (R package: glmnet) in the Discovery dataset to select compounds/metabolic features to build the linear regression model to predict gestational age. A 10-fold cross validation was performed to choose optimal lambda (penalty for the number of features), which determines the performance of the lasso model (number of features included in the model and prediction deviations). For the practical utility of a signature in potential clinical settings, in the identified compound-prediction models, if the number of predictors exceeded five under a given optimal lambda, we increased the lambda value so that the number of predictors is no more than five in the final models. We then used two different methods to evaluate the prediction deviation of our lasso model produced under a given lambda value: 1) A 10-fold cross validation within the Discovery cohort, in which the optimal lambda was used. In the CV, samples were distributed into folds by subject instead of by samples to prevent person-specific information cross-over between the training folds and the test fold. 2) Validation tests in the separate Test Set 1 and Test Set 2 cohorts, in which independent subjects were included and samples were analyzed one year and three years from the Discovery cohort. We built the model using the optimized lambda and full discovery datasets. This model was applied to the validation cohort for prediction and verification. A linear fitting from the two above evaluations were performed, between the predicted value and the actual values, with Pearson correlation coefficient (R), R2, and RMSE reported. The contribution of each predictor (metabolite) in each prediction model is defined by:

50i: the metabolite included in the linear model

51Unlike cross-validation in the discovery dataset, validation tests are not prone to hyperparametric selection bias for lambda value. Since samples from Test Set 2 cohort were normalized independently from other samples, a scaling was done in the end. Note that since the sequential nature of the data were not used in the machine-learning methods, other statistical tools, such as recurrent neural networks (e.g., LSTM), may be explored to improve the model.

52For samples collected after 28 weeks (third-trimester samples), we started with 264 level 1 and 2 compounds, and we used a similar discovery and validation pipeline described for predicting gestational age (above) to build logistic regression models predicting the categorical labels of gestational age > 32 and 37 weeks or delivery within 2, 4, 8 weeks. The prediction models for > 20, 24, and 28 weeks were built using samples from all three trimesters. For the prediction on delivery within 2, 4, and 8 weeks, only the 18 women (out of 30) with natural labor onset were included, excluding subjects with induction before labor onset and scheduled cesarean-section (induction by oxytocin/membrane strip after the onset is allowed). To estimate the confidence interval for each AUROC, we performed bootstrapping by person (instead of by samples) for 1000 times, and calculate the 95% confidence interval for AUROC.

53We evaluated the individual correlation between the predictions made by the metabolic clock and the estimations from first-trimester ultrasound: We first examined the correlation between metabolic clock predictions and the gestational age based on first-trimester ultrasound in individual persons. Each correlation was evaluated by Pearson’s correlation. We then performed meta-analysis across the persons to generate a summary p value, using Fisher’s method, to describe the overall correlation in each cohort (cross-validation in the Discovery, Independent validation of Test Set 1).

54Previous literature and our own observations suggest that birth weight and gestational length are positively correlated; later delivery is associated with a heavier absolute birth weight of an infant. To determine whether an infant’s birth weight falls above or below the group mean, we performed a linear regression between the two parameters and took the residuals to represent the birth weight deviation adjusted for delivery timing.

55Average Δ(GAmetabolic- GAultrasound): For each person, at each time point, we examined the differences between the metabolic clock and first-trimester ultrasound estimation of gestational age. These values were averaged for each person to represent the overall relative pace of metabolic clock compared to the first-trimester ultrasound estimation. We then examined the correlation between delivery timing adjusted birth weights and average Δ(GAmetabolomic- GAultrasound) (Figure 5C).

56To examine whether an accelerated metabolic clock (compared to the first-trimester ultrasound estimation) associates with advanced delivery, we performed the correlation between average Δ(GAmetabolomic- GAultrasound) and delivery timing, only in women with a natural labor onset (Figure 5D).