Pregnant women were enrolled between September 2013 and May 2016 at Stanford University, California (California cohort) and between December 2013 and February 2016 at the University of Alabama, Alabama (Alabama cohort). Urine samples were collected weekly from week 9 until delivery for women in the California cohort and week 15 until delivery for subjects in the Alabama cohort, and one sample collected postpartum. Dates of sampling and delivery were documented. Clinical information was collected. Ultrasound measurements at the first trimester were recorded. Women with normal, full-term pregnancies were included in the study cohorts. A normal, full-term pregnancy refers to a pregnancy at 37–41 weeks and without known complications. The study was approved by ethics committees at Stanford University and the University of Alabama, and written informed consents were obtained from all participants.

The study was conducted in two phases: 1) modeling to devise a metabolite-based estimation of GA during normal, full-term pregnancies; and 2) modeling to devise a metabolic panel predictive of time to delivery. In this study, the “gold” standard of GA was the ultrasound measurement based on the crown-rump length at the first trimester (Robinson, 1973). Metabolic concentrations in each urine sample were measured by global mass spectrometry (MS) analysis. Models that estimated GA and predicted time to delivery were developed using the California cohort and validated using the Alabama cohort. All statistical analyses were done in R software.

The urinary metabolites were extracted using a protein precipitation-based approach. Briefly, 10 µL of urine sample was extracted with 100 µL of methanol containing 5 μg/ml of ¹³C₅, ¹⁵N-l-proline^{, 13}C₆-l-arginine, and D₅-l-glutamine. These externally spiked exogenous isotope-labelled metabolites were used as references for sample preparation and extraction efficiency. The extract was then vortexed for 1 minute and centrifuged at 12,000 × g for 5 min 90 µL of supernatant was then collected for the global metabolomics analysis.

Following sample preparation, 10 µl of extract was injected onto a ZIC-HILIC column (2.1 mm × 100 mm × 3.5 μm) from EMD Millipore (Burlington, MA). The chromatographic separation was carried out using 10-mM ammonium acetate in water (A) and 10-mM ammonium acetate in acetonitrile (B) as mobile phases. The urine metabolites were eluted off the column using a linear gradient of 10–90% phase B over 7.5 min at 0.4 ml/min, and the column was re-equilibrated at 10% phase B from 7.6 to 10 min at 0.5 ml/min before the next injection. The column oven and autosampler temperatures were maintained at 40 and 4°C, respectively, throughout the analysis.

The data was acquired by a Q Exactive plus mass spectrometer in both electrospray positive and negative modes via a data-dependent manner. The source conditions were set at 3.5/−3.5 kV for Spray Voltages, 300°C for Vaporizer Temp, 300°C for Capillary Temp, 55.0% for S-Lens, 40 for Sheath Gas, 10 for Auxiliary Gas, and 0 for Sweep Gas. The spectra were acquired using following parameters: automatic gain control (AGC, MS1) = 1 × 10⁶, AGC (MS2) = 1 × 10⁵, Injection Time (MS1) = 100 m, Injection Time (MS2) = 50 m, Mass Range (MS) = 70–1,000 Da (Da), Resolution (MS1) = 70,000, Resolution (MS2) = 17,500, data dependent top-5 experiment (MS2) data acquisision, Isolation Window (MS2) = 1.0 Da, Dynamic Exclusion = 14 s, and Collision Energy (CE) = 25 eV. Between batch injections, the instrument was calibrated using external standards in both positive and negative modes to ensure a mass accuracy of less than one part per million (ppm).

A set of urine samples were collected longitudinally from a healthy volunteer across gestation, and the longitudinal urine samples were pooled to generate a urine sample for quality control (QC) purposes. Prior to the batch analysis, five replicates of QC urine extract were injected into the system for proper conditioning. Along the batch analysis, QC urine samples were analyzed repetitively at a frequency of one QC injection per 10 testing samples to allow the systematic assessment of the data quality. In addition, the signal responses of spiked internal standards (ISs) were checked batchwise to evaluate for the quality of metabolite extraction and sample preparation.

The detected metabolic features were extracted, aligned, integrated, grouped, and annotated using the XCMS package in R to obtain a matrix of metabolite intensities versus samples. The obtained features were normalized by a QC-based robust locally estimated scatterplot smoothing (LOESS) signal correction approach, and each feature was independently corrected by fitting a LOESS curve to the signal response measured in QC replicates injected repeatedly along the batch (Dunn et al., 2011). Features with coefficient variation (CV) ≤ 20% in QC and missing values in ≤30% of samples were selected for downstream analysis. After missing value imputation, the normalized intensities of qualified metabolites were divided by the creatinine values in a sample-wise manner to correct for the differential dilutions of urine samples.

Partial least-squares discriminant analysis (PLS-DA) was performed to characterize metabolic profiling of samples collected at the first, second, third trimesters and the postpartum period. The correlation coefficients of detected metabolites to GA were determined by the Pearson approach. A global false discovery rate (FDR) of 5% was applied to correct for the errors from multiple hypothesis testing. Correlations between the abundances of metabolites that were significantly associated with GA (Peasron’s r > 0.3 or r < −0.3) were presented using a two-dimensional self-organizing map (SOM), where metabolites were categorized into clusters with similar abundance patterns. Prior to the modeling process, data preprocess (average with 3 weeks window; 40 weeks). All statistical analyses were preformed using R packages (Heinemann, 2019).

The metabolic features were identified by matching their accurate precursor masses against a local database compiling a variety of public databases such as HMDB and MS-DIAL by using a tolerance of ± 5 ppm. Upon matching, the assigned metabolite identities were further validated by the retention time and fragmentation pattern from our local library containing 600 + authentic standards that were analyzed under the identical condition. The significant metabolites with validated identities were mapped to known Homo sapiens pathways in Kyoto Encyclopedia of Genes and Genomes (KEGG) database to characterize the pertaining pathway network and identify the meaningful pathways in the course of pregnancy development. Metabolic pathway-based pregnancy modeling was implemented as previously described (Sylvester et al., 2020) to estimate GA and time to delivery. Models were derived with the samples from the California cohort, and validated with the samples from the Alabama cohort. A tree based gradient boosting algorithm was utilized to construct the models (Chen and Guestrin, 2016). Model performance was assessed by Pearson’s correlation between model estimation and ultrasound results (for the GA dating model), and between model prediction and observed values (for the time-to-delivery predictive model). Errors were calculated with samples collected at the second trimester, the third trimester, and the combined period, respectively, and aggregated by subjects. Profiles of specific metabolites in urine were compared with those in sera.

Metabolite biomarker identification was performed as a Tier 1 or 2 identification with chemical standards according to MSI (Viant et al., 2017). With tandem mass spectrometry (MS/MS, Thermo Q Exactive plus) data of urine samples and manual review confirmation, the generated MS1/MS2 pairs were searched in the public databases: HMDB (http://www.hmdb.ca/), MoNA (http://mona.fiehnlab.ucdavis.edu/), MassBank (http://www.massbank.jp/), METLIN (https://metlin.scripps.edu), and NIST (https://www.nist.gov/). The metabolites of interest were procured and subjected to a Tier one identification comparing the retention time, MS1 and MS2 patterns with the biomarker candidates, using the same LCMS/MS protocol with the sample analysis.

To allow direct insight into the model, we implemented an “explainer”, based on Shapley values (Lundberg et al., 2020) from game theory, to explain the association between the input and predicted output. Our tree-based machine learning models such as random forests, decision trees, and gradient boosted trees are popular non-linear predictive models, yet the interpretability of tree-based models is low. Shapley value approach combines many high-quality local explanations to allow us to represent global structure while retaining local faithfulness to the original model. Shapley values used for explanations of the model were calculated with the SHAP package (https://github.com/slundberg/shap).

To investigate model performance between different machine learning approaches, we conducted comparative modelling with XGBoost, random forest, and elastic net classifiers for all combinations of predictive features and cohorts.

As shown in Figure 1, the California cohort consisted of 19 full-term pregnancies with 478 urine samples, and the Alabama cohort consisted of 10 full-term pregnancies with 171 urine samples. Each subject had 9 to 32 samples collected prior to delivery and one sample collected postpartum. Demographics of the two cohorts are shown in Supplementary Table S1 and S2. All the subjects in the California cohort were white while all the subjects in the Alabama cohort were black. Compared with subjects in the California cohort, subjects in the Alabama cohort were younger (p < 0.001: 32.2 vs. 26). The average gestational duration of California subjects is 39.5 weeks, longer than the Alabama moms’ (38.2 weeks). Compared with subjects in the California cohort, subjects in the Alabama cohort had higher pre-pregnancy body mass index (BMI; p = 0.01). The average BMI of California subjects is 21.8, ranging from 20.2 to 24.7, while that of Alabama subjects is 29.8, ranging from 26.6 to 32.5. All the Alabama subjects had previous pregnancies, while 47% (9 of 19) of the California subjects were in the first pregnancies.

A total of 5,473 metabolic features were identified by LC-MS-based global metabolomic profiling of samples in the California cohort (Supplymentary Figure S1). Of these, 1,716 features were selected with a missing value percentage less than 30% across all samples. These features were examined globally with PLS-DA (Figure 2), revealing a clustering pattern of ordered pregnancy progression as a function of the GAs in the first/second/third trimesters and the postpartum period.

A process of metabolite feature reduction was performed using the workflow in Supplymentary Figure S1. 885 features were found to be significantly correlated to gestations with a FDR [Benjamin and Hochberg method (Hochberg and Benjamini, 1990)] < 0.05, suggesting extensive urine metabolic changes occurred during the pregnancy progression. Filtered by an additional GA correlation criterion (Pearson’s r > 0.3 or < −0.3), a total of 119 metabolites were identified by LCMS profiling. 37/119 compound identities (Figure 3A) were determined by LC-MS/MS and compound library annotation analyses.

We analyzed the correlation in metabolite abundance along GA by characterizing and clustering the top 37 annotated metabolites. Based upon the existing structural and biological annotations, these 37 metabolites were categorized/clustered into five groups (Figure 3B): amino acid, carbohydrate, organic acid, quaternary amine, and steroid metabolisms.

The largest cluster was composed of various carnitine metabolites. These metabolites are indispensable for the transport of activated long-chain fatty acids from the cytosol to the mitochondrial matrix, where β-oxidation occurs, and the transfer of products of peroxisomal β-oxidation to the mitochondria for oxidation in the citrate cycle, the modulation of the acyl-coenzyme A (CoA)/CoA-ratio and the storage of energy as acetylcarnitine (McGarry and Brown, 1997; Rebouche and Seim, 1998; Steiber et al., 2004). We found that their abundance decreased along GAs (Figure 3C). The gestational patterns of the carnitine metabolites in maternal urine were similar to those in maternal plasma (Bargen-Lockner et al., 1981; Keller et al., 2009), which markedly declined to about half of the concentrations of non-pregnant women.

Positively-correlated patterns were observed between carnitine and basic amino acid clusters (Figure 3B). Such findings are in line with the fact that the endogenous biosynthesis of carnitine involves a complex series of reactions with lysine providing the carbon backbone (Keller et al., 2009). In the amino acid cluster, deoxypyridinoline had a positive correlation with GA, while taurine had a negative correlation with GA (Figure 3C).

Positively correlated patterns were also observed between steroid and carbohydrate metabolisms (Figure 3B). Metabolites associated with these two clusters, including 5-alpha-pregnan-3beta, 20alpha-diol disulfate, dehydroepiandrosterone sulfate (DHEA-S), estriol-3-glucuronide, pregnanediol 3-O-glucuronide, 3′-sialyllactose, and α-lactose, were all positively correlated with GA (Figure 3C).

The 37 annotated metabolites were mapped to five pathways: bile secretion, steroid hormone biosynthesis, pantothenate and CoA biosyntheses, benzoate degradation, and phenylpropanoid biosynthesis. Metabolic-pathway-based GA dating models was constructed based on the five pathways, using XGBoost and random forest approaches (Figure 4A). XGBoost outperforming random forest (Figure 4A), bile secretion was the most important feature in the XGBoost modeling (Figure 4B). Model performance was assessed with each cohort (California cohort: Pearson’s R = 0.96; Alabama cohort: Pearson’s R = 0.76; Supplymentary Figure S2) and with each subject [California cohort: Pearson’s R = 0.97 (0.96, 0.98); Alabama cohort: Pearson’s R = 0.84 (0.77, 0.87); Figure 4C]. Summary plot of Shapley values computed for each patient individually in the test partition. The pathways are sorted top-down based on their global contribution which is in line with Figure 4D. The distance of a dot representing a sample from the vertical line indicates its contribution. The color of a dot indicates feature value for that sample. Blue and pink color represent extreme values of the feature. Shapley values on the right side of vertical axes “push” predictions towards the high GA correlation and those on the left side towards the low correlation.

Six metabolic marker compounds of pregnancy GA dating were procured and subjected to a Tier 1 identification comparing the retention time, MS1 and MS2 patterns with the biomarker candidates, using the same LCMS/MS protocol with the urine sample analysis: acetylcholine, estriol-3-glucuronide, DHEA-S, α-lactose, hydroxyhexanoycarnitine, and l-carnitine (Figure 5). GA dating models were constructed based on the six metabolic markers, using XGBoost, random forest and elastic net approaches (Figure 6A). XGBoost outperforming random forest (Figure 6A, R, R ², and RMSE analysis), DHEA-S and α-lactose were two most important features in XGBoost modelling (Figure 6B). Model performance was assessed with each cohort (California cohort: Pearson’s R = 0.95; Alabama cohort: Pearson’s R = 0.79; Figure 6C) and with each subject [California cohort: Pearson’s R = 0.97 (0.96, 0.98); Alabama cohort: Pearson’s r = 0.87 (0.82, 0.92); Figure 6D]. Differences between the model estimation and ultrasound results were calculated (Figure 6C). There were 73.7% (14 of 19 subjects) of the estimates in California cohort and 40.0% (4 of 10 subjects) of the estimates in Alabama cohort within ±1 week of the ultrasound results. Summary plot (Figure 6E) of Shapley values computed for each patient individually in the test partition. The metabolite biomarkers are sorted top-down based on their global contribution which is in line with Figure 6B. The distance of a dot representing a sample from the vertical line indicates its contribution. The color of a dot indicates feature value for that sample. Blue and pink color represent extreme values of the feature. Shapley values on the right side of vertical axes “push” predictions towards the high GA correlation and those on the left side towards the low correlation.

We developed and validated a separate model predictive of time to delivery using the 6 GA dating markers. Like the GA dating model, DHEA-S and α-lactose were top features in predicting time to delivery (Figure 7A). In the second trimester, 52.6% (10 of 19) of the prediction in the California cohort and 30.0% (3 of 10) of the prediction in the Alabama cohort had errors of less than 1 week, which was similar to the performance of the expectation by the first-trimester ultrasound results (California: p = 0.9; Alabama: p = 0.2; Figure 7B).