INTRODUCTION
As substrates and products of metabolism, metabolites drive essential cellular functions, including energy production and storage, signal transduction and apoptosis
1. The heritability of circulating metabolite levels ranges from 10.5% to 93.2%, with a median value of 48.8%
2. With advancements in whole-genome sequencing and metabolomics, genome-wide association studies for metabolites (mGWASs) have successfully identified numerous loci associated with metabolic traits. These findings provide novel insights into the biology of metabolism and offer opportunities to identify pathways that underlie the links between genetic determinants and health-related phenotypes
3-8.
Pregnancy involves a wide variety of physiological changes and remarkable metabolic shifts
9. Core pathways, such as bile acid biosynthesis, arachidonic acid metabolism, and pyrimidine, purine, glutamate, and tyrosine metabolism, exhibit significant changes during pregnancy, influencing traits in pregnant women and their children
10. While the impact of diet and lifestyle on the maternal metabolome has been partially investigated, the role of genetics in shaping maternal metabolome during pregnancy has not yet been systematically elucidated. Moreover, the use of assisted reproductive technology (ART) has steadily increased worldwide, with ART-conceived babies now accounting for approximately 3% of total births in China
11,12. ART-conceived pregnancies have been associated with increased risks of pregnancy complications and adverse outcomes, including gestational diabetes mellitus (GDM), hypertensive disorders of pregnancy (HDP) and preterm birth, many of which are linked to metabolic disturbance
13. These findings raise the important question: do naturally conceived (NC) pregnancies and ART pregnancies exhibit differences in the genetic regulation of the maternal metabolome?
Thus, we performed mGWASs of 748 qualified plasma metabolites across early, mid and late pregnancy in 2,342 pregnant individuals from the Jiangsu Birth Cohort (JBC) Study, including 1,525 with NC pregnancies and 817 with ART pregnancies (Supplementary Fig. S1). This study systematically characterizes the genetic regulation of the maternal metabolome across gestation, identifies candidate pregnancy context-dependent genetic effects, and evaluates the extent to which genetic architecture is shared between NC and ART pregnancies. By providing a large-scale genomic resource for maternal metabolism during pregnancy, our findings lay the groundwork for future mechanistic studies and precision medicine strategies in maternal and reproductive health.
RESULTS
mGWASs of pregnant women on common and low-frequency variants identify 1,080 associations for 385 metabolites
Utilizing whole-genome sequencing data and untargeted metabolomic data in plasma samples collected in the first, second and third trimesters, we constructed a genomic atlas of plasma metabolite levels throughout pregnancy. The characteristics of women with NC pregnancies and ART pregnancies are shown in Supplementary Table S1. Details on the quality control of whole-genome sequencing data and untargeted metabolomic data are provided in the Materials and Methods section. A total of 748 metabolites that passed quality control were included in the subsequent analyses and classified into 9 superclasses and 107 subclasses based on the Human Metabolome Database (HMDB)
14 (Supplementary Table S2). We performed mGWASs for the trimester-specific maternal metabolome across more than eight million common (minor allele frequency (MAF) > 5%) and low-frequency (1% < MAF ≤ 5%) variants in the NC and ART pregnancies, respectively. Notably, early-pregnancy samples were collected from NC women at gestational weeks (GW) 8–14, while ART women provided samples at GW 4–6. To evaluate the potential influence of this difference in early pregnancy sampling time, we compared the effect sizes and
P values of significant associations both with and without adjustment for gestational week. The high consistency (
r > 0.99) observed suggests that the discrepancy in early-pregnancy sample collection is unlikely to have affected the results significantly. Results from the two groups were then combined using inverse variance meta-analysis. The magnitude and direction of most significant associations in meta-analysis were highly consistent across NC and ART groups, with effect sizes correlating at
r = 0.92 and 99.81% of associations showing concordant direction (Materials and Methods; Supplementary Fig. S2). Single-variant association tests in the NC and ART datasets, as well as in the meta-analysis, identified a total of 1,025 significant and independent variant-metabolite associations (Bonferroni-corrected genome-wide significance
P < 5.56 × 10
−10, linkage disequilibrium (LD)
r2 < 0.1) for 385 metabolites across 232 loci (Materials and Methods; Fig. 1; Supplementary Table S3). Inflation was well controlled (genomic control inflation factor median = 1.00, range = 0.96–1.04; Supplementary Table S4). Conditional analyses were also conducted within metabolite-specific regions, yielding 55 secondary signals (Materials and Methods; Supplementary Table S5). On a global scale, the genetic regulation of metabolite levels was generally consistent across three trimesters, with the estimated effect sizes of the total 1,080 associations (primary plus secondary, corresponding to 805 index variants ) showing high correlations across trimesters (1
st trimester vs 2
nd trimester,
r = 0.90; 2
nd trimester vs 3
rd trimester,
r = 0.94; 1
st trimester vs 3
rd trimester,
r = 0.90; Fig. 1; Supplementary Fig. S3).
Overall, approximately 40–60% of the analyzed metabolites exhibited independent variant-metabolite associations, with the exception of xenobiotics (32.0%) and energy (14.3%) (Fig. 2a). An inverse relationship was observed between effect size and MAF, with low-frequency variants (1% < MAF ≤ 5%) generally exhibiting larger effect sizes. Specifically, 327 (30.3%) independent associations had large absolute effect sizes (β) (> 0.5 standard deviations per allele), of which 230 (70.3%) were low-frequency variants (Fig. 2b; Supplementary Table S3). Functional annotation using the Ensembl Variant Effect Predictor (VEP)
15 indicated that over 60.0% of the independent variant-metabolite associations involved index variants located in intronic or intergenic regions. Notably, 59 index variants (11.1% of the associations) had a direct functional consequence on the corresponding transcript, such as loss-of-function or missense variants (Fig. 2c).
We then systematically assessed protein-coding genes within 1-MB windows centered on index variants to align gene functions with relevant metabolites in the 232 loci corresponding to the 1,080 independent associations identified in this study. This analysis identified 141 loci with plausible or established biochemical links, implicating 166 unique putative causal genes (Supplementary Table S6), the majority of which encode enzymes (66.9%) and transport proteins (22.3%) (Fig. 2d). Approximately 60% (n = 140) of these loci were linked to a single metabolic trait (Fig. 2e). In contrast, the locus on chromosome 11, involving the fatty acid desaturase gene family, was associated with 45 lipid metabolites, exhibiting notable pleiotropy. Extensive pleiotropy was also evident in other loci. For example, within-class pleiotropy was observed at the locus on chromosome 15, which contains the Lipase C coding gene (LIPC) and is associated with 18 lipid metabolites. Across-class pleiotropy was identified at the locus on chromosome 12, which contains the Solute Carrier Organic Anion Transporter Family genes (SLCO1A2, SLCO1B1, SLCO1B3, and SLCO1B7). This locus was linked to metabolites in both the cofactors and vitamins superclass and the lipids superclass (Fig. 2e; Supplementary Table S3). In addition, 242 out of the 385 metabolites (62.9%) were associated with a single locus (Fig. 2f). The remaining metabolites demonstrated polygenicity, with 38 out of 385 metabolites (9.9%) being influenced by more than two loci. O-acetylhomoserine was associated with the greatest number of loci, although no putative causal genes were identified for these associations (Fig. 2f).
To evaluate the novelty of the variant-metabolite associations identified in our study, we performed cross-study comparisons with 12 previously published mGWASs based on untargeted metabolomic profiling using the Metabolon platform, each with a minimum sample size of approximately 1,000
2, 7, 8, 16-24 (Materials and Methods; Supplementary Table S7). Among the 1,080 associations identified in the present study, we excluded replicated signals, resulting in 321 novel associations involving 279 index variants (Materials and Methods; Supplementary Table S8). Using a purely regional definition approach, we identified 82 additional loci (Materials and Methods; Supplementary Fig. S2), further expanding our understanding of the genetic landscape of metabolite regulation (Supplementary Table S8).
The 12 prior studies collectively reported 8,500 unique significant genetic associations involving 1,737 metabolites (including both known and unknown metabolites) (Supplementary Table S9). Of these, 601 metabolites were assayed and analyzed in the present study, corresponding to 3,258 associations and 2,042 variants. We then validated these associations and found that 2,587 previously reported associations (79.4%) were replicated in our dataset, with concordant effect directions and
P < 0.05 (Supplementary Table S10). Additionally, a recently published genome-wide association study (GWAS) by Liu
et al. investigated the associations between variants and 60 maternal metabolites during pregnancy using a targeted approach
25 (Supplementary Table S11). Fourteen of these metabolites were detected in our study, corresponding to 20 associations. Among these, 18 associations reached statistical significance (
P < 0.05) and demonstrated consistent effect directions with our findings (Supplementary Table S11). These replicated findings further underscore the reliability of our dataset.
Ancestry-specific analyses identify new associations in the East Asian population
Given that most published mGWASs have been conducted in individuals of European ancestry, we assessed whether the 321 novel associations were primarily identified due to ancestry specificity or pregnancy status specificity. To do this, we compared MAFs of the 279 index genetic variants in East Asian populations with those in other ancestry groups using data from the 1000 Genomes Project Phase 3 (1KGP3) (Supplementary Table S12). Our analysis identified 68 variants (corresponding to 78 associations) that exhibited substantial differences in genetic architecture. Of these, 34 variants (43 associations) were found exclusively in the East Asian population, while 34 variants (35 associations) displayed substantially higher MAFs (> 10-fold) in East Asians compared with other ancestry groups.
Furthermore, we assessed the novel associations in 844 participants with available pre-pregnancy metabolomic data (Supplementary Fig. S1). Of the 78 associations, 14 reached statistical significance in pre-pregnancy samples and showed no heterogeneity in effects compared with analyses performed in pregnancy samples (Materials and Methods; Supplementary Fig. S2). This ruled out the potential influence of pregnancy status and confirmed ancestry specificity (Fig. 3a; Supplementary Table S13). As an example, a missense variant in the coding region of
SLC10A1, Chr14:70245193:G > A (rs2296651, c.800 C > T/p.Ser267Phe), exclusively present in East Asians (MAF = 0.07), was associated with increased levels of six bile acids (Fig. 3b; Supplementary Table S13). The ancestry-specific associations between Chr14:70245193:G > A and bile acid levels were also recently demonstrated in the Born in Guangzhou Cohort Study in China, where the bile acid levels were measured through clinical serum biochemical tests
26.
Novel associations highlight candidate pregnancy context-dependent associations
To determine whether the identification of the 321 novel associations was due to specific or intensified genetic regulation in metabolism in pregnancy, we applied interaction analyses to assess whether the relationship between genetic variants and metabolites changed by pre-pregnancy/pregnancy status, using a significance threshold for interaction at P < 0.1 (Materials and Methods). Significant interactions between variant and pregnancy/non-pregnancy status were observed for 174 signals. Furthermore, 65 associations were identified as potential pregnancy-specific or intensified associations, as they exhibited significant effects during pregnancy, while the associations were attenuated or nullified before pregnancy. To ensure rigor, we further tested the heterogeneity between the effects of the 65 variant-metabolite associations observed pre-pregnancy and during pregnancy. Of these, estriol 3-sulfate showed a pregnancy-specific association, while the remaining 13 associations (involving 11 variants) exhibited significant heterogeneity and had sufficient statistical power (power > 0.9) to detect pre-pregnancy significance but showed no such association before pregnancy. These 14 associations (involving 12 variants) showed significant heterogeneity and were ultimately identified as candidate pregnancy context-dependent associations (Supplementary Fig. S4; Table S14).
The 14 candidate pregnancy context-dependent associations implicated 11 unique putative causal genes, including seven encoding enzymes and two encoding transport proteins (Supplementary Table S14). One notable association was between estriol 3-sulfate and the missense variant Chr12:21331549:T > C (rs4149056, c.521T > C /p.Val174Ala) in SLCO1B1. This association was null preconceptionally (β = −0.24, P = 0.425, effect allele: C) but became significant during early pregnancy (β = 0.42, P = 2.69 × 10−20, effect allele: C), middle pregnancy (β = 0.60, P = 3.92 × 10−43, effect allele: C), and late pregnancy (β = 0.56, P = 2.32 × 10−39, effect allele: C; Fig. 4a). Maternal circulating estriol 3-sulfate levels varied by rs4149056 genotype, with participants carrying the C allele showing significantly higher estriol 3-sulfate levels, compared with those with T allele, in both NC and ART populations (Fig. 4b). In the majority of pre-pregnancy samples (94.1%), estriol 3-sulfate was below the limit of detection, while it was specifically detected in maternal plasma collected during pregnancy (Fig. 4c). This finding is consistent with the fact that estriol 3-sulfate is synthesized by the fetal-placental unit. Putative causal genes in the region, including SLCO1A2, SLCO1B1, SLCO1B3, and SLCO1B7, encode members of the organic anion transporter family, which may be involved in the transport of estriol 3-sulfate (Fig. 4d). Leveraging the follow-up data from this cohort, we found that elevated plasma levels of estriol 3-sulfate were in relation to fetal growth, with the top tertile of estriol 3-sulfate levels in mid- and late pregnancy being associated with 0.52- and 0.93-fold increased risk of large-for-gestational-age (LGA) offspring (2nd trimester: OR = 1.52, 95% CI: 1.11–2.10; 3rd trimester: OR = 1.93, 95% CI: 1.40–2.67, Fig. 4e).
The metabolites involved in the remaining 13 associations were detected in both non-pregnant and pregnant individuals, but the metabolite-variant associations were observed only during pregnancy (Supplementary Table S14). For instance, the intronic variant Chr2:158975570:C > T (rs12987618) in
UPP2 (encoding Uridine Phosphorylase 2), identified as an expression quantitative trait locus (eQTL) for
UPP2 in liver tissues from the Genotype-Tissue Expression (GTEx) database, was significantly associated with decreased levels of 5-methyluridine throughout pregnancy (early pregnancy: β = –0.46,
P = 4.31 × 10
−25; middle pregnancy: β = –0.43,
P = 3.40 × 10
−20; late pregnancy: β = –0.43,
P = 5.64 × 10
−21, effect allele: T). However, the effect size of this association was considerably attenuated and did not reach significance before pregnancy (β = –0.22,
P = 0.002, effect allele: T) (Fig. 5a, b; Supplementary Table S14). Our data demonstrated enhanced 5-methyluridine degradation during pregnancy (Fig. 5c). 5-methyluridine, also known as ribothymidine, is a nucleoside and a modified form of RNA. Given that recent evidence showed that the degradation product ribose-1-phosphate supports energy supply in rapidly proliferating cells
27, 28, it is biologically plausible that the degradation process of 5-methyluridine might be crucial for the rapid cellular proliferation necessary for pregnancy establishment and maintenance. To test this hypothesis, we conducted a nested case-control study within our cohort, which included 201 cases of biochemical pregnancy loss and 781 controls with live births (Materials and Methods; Fig. 5d). We observed that the top-tertile levels of 5-methyluridine in early pregnancy were significantly associated with an increased risk of biochemical pregnancy loss after the adjustment for conventional covariates (OR = 1.58, 95% CI: 1.07–2.34; Fig. 5e). After adjustment for the kidney and liver function indexes measured in the clinical serum biochemical tests (including blood creatinine level, blood urea nitrogen level, and levels of alanine aminotransferase and aspartate aminotransferase), the associations remained unchanged.
Genetic determinants of metabolite dynamic changes across gestation
Given the substantial metabolic changes that occur during pregnancy, identifying the genetic determinants influencing these metabolic trajectories can provide valuable insights into the genetic regulation of maternal metabolism. To assess metabolic changes, we calculated trimester-to-trimester ratios within the same individual, yielding three ratios for each metabolite (i.e., 2nd trimester/1st trimester, 3rd trimester/2nd trimester, 3rd trimester/1st trimester). Genome-wide association analyses of these ratios revealed 56 independent associations at a multiple-testing-adjusted genome-wide significance threshold of P < 1.67 × 10−9 (Fig. 6a). Assessment of genomic inflation factors indicated no excessive test statistic inflation or population stratification (genomic control inflation factor median = 1.00, range = 0.87–1.10; Supplementary Table S15). Among these 56 associations, we identified six significant signals involving five hormone-related metabolite trajectories, including two androgenic steroids, two estrogenic steroids, and one progestin steroid, indicating that gestational changes in hormone-related metabolites are under genetic regulation (Supplementary Table S16). To further validate these 56 signals, we applied a mixed-effects regression model incorporating interactions between genetic variants and gestational weeks to assess their temporal effects during pregnancy. Our analysis revealed that 47 out of 56 signals exhibited significant interactions with gestational weeks of sampling (Materials and Methods; Supplementary Fig. S2). The annotation of these signals identified 28 out of 47 signals as being linked to ten putative causal genes, including seven metabolic enzymes and two transporters (Supplementary Table S16).
We then visualized the trajectories of 39 metabolites implicated in the 47 signals across gestation by genotypes. Notably, the
UGT3A1 intronic variant (Chr5:35983283:CA > C, rs78247641) exerted a concerted regulatory effect on the dynamic changes of six metabolites throughout pregnancy. These metabolites included two bile acids (i.e., glycochenodeoxycholate glucuronide [GCDCA glucuronide] and glycodeoxycholate glucuronide [GDCA glucuronide]), and three dicarboxylate fatty acids (i.e., octadecenedioate [C18:1-DC], hexadecanedioate [C16-DC] and tetradecanedioate [C14-DC], Fig. 6b; Supplementary Table S16). Prior studies have shown that bile acids and fatty acids can act as substrates for UGT3A1
29. Further analyses using our follow-up data (
n = 2,070) indicated that elevated levels of C14-DC, GCDCA glucuronide and GDCA glucuronide were associated with an increased risk of preterm birth (Fig. 6c). These results suggest that
UGT3A1 locus pleiotropically regulates multiple metabolites, which may play a role in pregnancy progression.
Different associations by obstetric conditions suggest heterogeneous regulatory effects
Obstetric factors, such as maternal age and pre-pregnancy body mass index (BMI) are known to shape maternal metabolic profiles and may modify genetic associations with circulating metabolites. To assess whether metabolite-variant associations vary across various obstetric conditions, we conducted stratified analyses by maternal age, pre-pregnancy BMI, parity and GDM status (Materials and Methods; Supplementary Fig. S2). Our findings revealed heterogeneity in three metabolite-variant associations across BMI subgroups. Specifically, the variant Chr9:6665010:C > T (rs1658972) was significantly associated with reduced mid-pregnancy levels of 3-methylglutarylcarnitine (β = –0.83, P = 4.42 × 10–45, effect allele: T) in individuals with a BMI < 24 (underweight or normal). However, this effect was substantially attenuated in overweight or obese women (β = –0.34, P = 2.56 × 10–3 in BMI ≥ 24, effect allele: T; Supplementary Table S17). This heterogeneity was further supported by the significant interaction between BMI and rs1658972 (P for interaction = 1.45 ×10–4). We further identified significant heterogeneity in three metabolite-variant associations between nulliparous and multiparous women, with stronger genetic effects on metabolite concentrations observed in multiparous women. In addition, two metabolite-variant associations showed significant heterogeneity according to GDM status (Supplementary Table S17). However, no significant interactions were found between genetic variants and maternal age in association with maternal metabolite levels. Given the low incidences of HDP and preterm birth in our cohort, sensitivity analyses were conducted by excluding these samples, which showed consistent results with the overall analysis.
For the associations that reached study-wide significance in either the NC or ART dataset, we assessed heterogeneity in the effects between the two conception types. While the associations between genetic variants and metabolites are generally consistent in the NC and ART women (Fig. 7a), we identified 55 associations exhibiting significant heterogeneity across the datasets (P for heterogeneity test after correction for false discovery rate (FDR) < 0.05; Fig. 7b; Supplementary Table S17). This heterogeneity is likely partly due to the procedures and medications used in ART (Supplementary Table S18). For example, the index variant at the NAALAD2 locus (Chr11:89899487:A > G, rs7951089), an eQTL of NAALAD2 across multiple tissues (GTEx database), was associated with decreased beta-citrylglutamate levels in early pregnancy in the NC population (β = –0.38, P = 2.96 × 10−23, effect allele: G), whereas this association was considerably attenuated in the ART women (β = –0.10, P = 0.039, effect allele: G; Fig. 7c, d). This discrepancy may be attributed to gonadotrophin-releasing hormone antagonists (GnRHA) used for ovulation induction in ART. After GnRHA treatment, beta-citrylglutamate level in individuals with the AA genotype dropped to match those of the individuals with the GG genotype, thereby nullifying the association (β = –0.02, P = 0.789, effect allele: G, Fig. 7e; Supplementary Table S18). Follow-up data indicated that elevated beta-citrylglutamate was associated with an increased risk of GDM (OR = 1.57, 95% CI: 1.23–1.99; Fig. 7f). Thus, individuals with the A allele at Chr11:89899487 appear to carry a higher risk of GDM, with GnRHA treatment potentially mitigating this risk, specifically in women with the AA genotype (Fig. 7g). The differences in associations between the NC and ART populations may also stem from distinct physical conditions, which warrants further investigation in large sample sizes.
Rare variants affecting metabolite levels
To better understand the contribution of rare variants to maternal metabolite levels, we performed a gene-based burden analysis focusing on coding sequence variations (Materials and Methods; Supplementary Fig. S2). This analysis identified 268 statistically significant metabolite-gene associations, involving 211 unique genes and 191 metabolites, at a significance threshold of P < 3.44 × 10−6 (Bonferroni correction for 4,847 genes and three trimesters, Supplementary Table S19). We also found overlaps between the burden test results and single variant analysis. For example, a common variant (rs147919763) tagging the ACY1 gene has been linked to N-acetylmethionine levels in our single variant analysis. Consistently, our gene-burden analysis also showed the association of ACY1 with maternal N-acetylmethionine levels (P = 2.00 ×10−29), further validating the relationship between ACY1 and N-acetylmethionine. ACY1 encodes an enzyme involved in the catabolism and salvage of acylated amino acids. Interestingly, the burden test revealed new associations between ACY1 and seven additional N-acyl metabolites, including N-acetylserine, N-acetylalanine, N-acetylglutamate, N-acetylglutamine, N-acetylglycine, N-acetylvaline, and N-acetylthreonine (Supplementary Table S19). In addition, we identified significant gene-metabolite associations for 85 metabolites that did not show significant associations in single variant analysis. For example, the TDO2 gene was associated with tryptophan levels in the burden test (P = 7.94 × 10−7), though it was not significantly associated with any single variant. These findings add to our knowledge of the genetic determinants of maternal metabolism during pregnancy.
Genetic-predicted metabolite levels are associated with pregnancy outcomes
To explore the relationship between metabolite levels and pregnancy outcomes, we imputed polygenic score (PGS)-predicted metabolite levels in 7,385 pregnant women using weighted genetic scores. We then evaluated the associations of these predicted metabolite levels with seven pregnancy-related phenotypes, including GDM, preeclampsia, birthweight, gestational duration, preterm birth, LGA, and small for gestational age (SGA) (Materials and Methods; Supplementary Fig. S2). Although no associations reached statistical significance after the correction for FDR, we identified 120 associations between 51 genetically predicted metabolite levels and phenotypes that showed nominal significance (P < 0.05; Supplementary Table S20), suggesting potential causal relationships between metabolites and pregnancy outcomes.
For example, consistent with previous studies, our analysis found that PGS-predicted maternal retinol levels were associated with a decreased risk of GDM (OR = 0.44, 95% CI: 0.29–0.66,
P = 1.08 ×10
−4; Supplementary Table S20). This supports the proposed causal role of retinol in reducing GDM risk
30. In addition, PGS-predicted dimethylglycine levels across pregnancy were associated with a reduced risk of preterm birth, while PGS-predicted succinate levels were associated with an increased risk of preterm birth (Supplementary Table S20). These results align with findings from previous observational studies and reinforce the potential causal effects of metabolites on preterm birth
31.
DISCUSSION
Pregnancy is a unique and critical period characterized by significant metabolic adaptation and regulation. Despite the importance, the genetic determinants of maternal metabolome have not been systematically studied. In this study, we combined whole-genome sequencing with untargeted metabolomic profiling on 2,342 women throughout their pregnancy to create a comprehensive genomic atlas of the cross-gestation metabolome. Our analysis identified 1,025 genome-wide significant associations and 55 secondary associations, involving 385 circulating metabolites and 805 index variants. Notably, 321 of these variant-metabolite associations were newly identified, among which 14 mQTL were determined as ancestry-specific. Specifically, this study identified pregnancy context-dependent associations, suggesting potential new or intensified genetic regulation on metabolism during pregnancy. Leveraging the maternal metabolome data across pregnancy, we additionally identified 47 mQTL that show associations with metabolite-changing trajectories throughout pregnancy. The changing intensity of metabolites, rather than cross-sectional level, may also have effects on pregnancy progression and development of complications. Furthermore, we found that variant-metabolite associations varied by obstetric conditions, partially attributable to physical conditions or clinical interventions, offering insights into future strategies in personalized prenatal care.
The genetic architecture of the associations in the present study revealed that 14 novel associations were attributable to ancestry-specificity. This suggests that mGWASs conducted in diverse populations can expand our understanding of ancestry-specific variation in metabolic traits. In addition, discovering ancestry-specific mQTL might provide new insights for biomedical and pharmaceutical research. Intrahepatic cholestasis of pregnancy (ICP) is a significant obstetric complication, of which the clinical feature is elevated serum total bile acid levels
32. In East Asians, we found that rs2296651 in
SLC10A1 is significantly associated with bile acid levels, providing a new potential drug target for ICP and other bile acid metabolism disorders.
In addition to the documentation of genetic associations with broader metabolic traits, our study identified novel associations that highlight the specificity of metabolic regulation during pregnancy. Given that naturally conceived participants were recruited in early pregnancy, pre-pregnancy metabolomic profiles were only available for ART women, creating a data imbalance that hinders validating pregnancy context-dependent genetic effects and raises the risk of spurious signals due to limited pre-pregnancy sample size. We integrated interaction analysis and statistical power evaluation to filter credible signals and finally identified 14 candidate pregnancy context-dependent associations. One such novel association involves estriol 3-sulfate, a metabolite derived from estriol synthesized by the fetal-placental unit
33 and detected only during pregnancy. The SLCO gene family, whose members encode organic anion transporters, is implicated in these associations. The genetic variants may impair the transport of estriol 3-sulfate metabolites, resulting in the accumulation of estriol 3-sulfate in maternal plasma.
Some variant-metabolite associations were presented in both pregnant and non-pregnant status but notably amplified during pregnancy, pointing to context-dependent genetic regulation shaped by gestational physiological conditions, which warrants further validation in future investigations. The
UPP2 locus and its association with 5-methyluridine exemplify this phenomenon. While we cannot completely rule out the presence of this association in non-pregnant individuals, the
UPP2 locus showed a much stronger effect on 5-methyluridine during pregnancy, implying altered regulation in this context.
UPP2 is minimally expressed in adult organs and embryonic tissues
34, highlighting its pregnancy-specific activity. 5-methyluridine, a pyrimidine nucleoside structurally similar to uridine, can be salvaged to meet energy requirements through its conversion into ribose-1-phosphate and subsequently into fructose-6-phosphate and glyceraldehyde-3-phosphate via the non-oxidative branch of the pentose phosphate pathway. These intermediates are then utilized in glycolysis to fuel ATP production, biosynthesis, and gluconeogenesis
28. Recent evidence suggests that in conditions of glucose limitation, alternative nutrients like uridine may serve as a fuel in proliferative cells, such as cancer cells
27. During pregnancy, glucose is preferentially allocated for fetal growth and development, limiting its supply to the placenta
35. It is biologically plausible that 5-methyluridine serves as an alternative fuel for placental development when glucose is scarce. Disruption in nutrient allocation or metabolism within placenta is often linked to adverse pregnancy outcomes
36. In an independent case-control study, we found that elevated levels of 5-methyluridine were associated with a 58% increased risk of biochemical pregnancy loss. These findings underscore the critical role of 5-methyluridine degradation in the early establishment and maintenance of pregnancy.
The unique design of the JBC allowed for the identification of associations showing heterogeneous genetic effects between NC and ART subcohorts. While the associations between genetic variants and metabolites are generally consistent in the NC and ART women, we observed 55 associations that displayed considerable heterogeneity. In particular, rs7951089, an eQTL of
NAALAD2, was associated with beta-citrylglutamate levels in the NC women but not in the ART women. This heterogeneity was partially attributable to clinical interventions in the ART process. The treatment with GnRHA in ART women carrying the A allele of rs7951089 resulted in a reduction in beta-citrylglutamate level to a point comparable to that of those with the G allele. Previous research suggested that the NAALADase family members are beta-citrylglutamate hydrolases
37. Thus, the functional variant of
NAALAD2 and GnRHA may interact to affect circulating beta-citrylglutamate, a known glutamate metabolite. Glutamate, which modulates the function and viability of endocrine cells in pancreatic islets, plays a role in the regulation of glucose homeostasis in diabetes
38, 39. Our follow-up data suggest that GnRHA treatment specifically lowers beta-citrylglutamate levels in women carrying the A allele of rs7951089, which is consequently related to a reduced risk of GDM. Such findings have clinical significance, as they could provide opportunities for precision intervention strategy in the ART process based on genotype-specific considerations.
The mGWASs primarily focus on common variants
16, 21, 40, which are mostly non-coding and have undetermined functions. Protein-coding variants, often with very low frequencies, tend to have direct biological effects. We additionally performed gene-centric rare variant analyses to investigate whether rare variants, in aggregate, affect metabolite regulation. Notably, 85 metabolites did not show significant associations with single variants while were associated with certain genes through the effects of rare variants, thereby adding additional knowledge to the genetic regulation of metabolome.
The present study investigated the genetic regulation of metabolome across pregnancy and identified genetic determinants of maternal metabolites that might be implicated in maternal and child health. However, several limitations should be noted. First, the untargeted plasma metabolome approach limits the accurate quantification of metabolites. Second, potential selection bias may exist in our study population. ART women were diagnosed with infertility causes, whereas NC women were naturally fertile, which suggests potential differences in their genetic background and metabolic characteristics. While benefiting from this design, we reported signals consistently exist in both ART and NC women, and also presented associations exhibiting significant heterogeneity. However, due to the limited case number of each infertility cause, we are not able to identify specific signals stemming from distinct physical conditions, which warrants further investigation in large sample sizes. Third, similar to other GWASs, the reported signals with rigorous P values were subject to inflation in the effect sizes. Further studies are warranted to validate the identified effect sizes. Fourth, pre-pregnancy metabolomic data were only available in ART women but not in the NC group due to the fact that they were enrolled in the cohort during early pregnancy, resulting in limited statistical power to detect pregnancy context-dependent genetic effects and restricted ability to validate corresponding interaction effects in the naturally conceived population. Future studies incorporating additional pre-pregnancy metabolomic data are thus warranted to address this gap. Lastly, the findings on novel associations should be considered as hypothesis-generating and warrant replications and further investigations in future studies to determine the causality and uncover the underlying mechanisms.
Taken together, our findings advance the mechanistic understanding of genetic regulation of the maternal metabolome throughout human pregnancy, deepening insights into the biological mechanisms underlying maternal metabolic adaptation. Specifically, we demonstrate that maternal metabolism is governed by common and rare genetic variants, and undergoes candidate pregnancy context-dependent associations, through which we identified key metabolites and genes that may be critically involved in pregnancy maintenance and fetal growth (Supplementary Fig. S5). These findings may provide valuable implications for the development of improved prenatal care and personalized prevention strategies for pregnant women.
MATERIALS AND METHODS
Study design and participants
The present study included 1,525 women with natural conceptions and 850 women conceived via assisted reproductive techniques from the JBC Study (Supplementary Fig. S1). These participants were randomly selected from eligible fully followed-up participants while preserving the original proportion of spontaneous and ART conception populations within the cohort. The JBC Study was originally designed to investigate the disparities between ART and NC in relation to perinatal outcomes and child health
41. All participants in the present study underwent genome-wide genotyping, and their circulating metabolic profiles were measured in early, mid or late pregnancy. The study protocol was approved by the Human Research Ethics Committee of Nanjing Medical University (NJMUIRB(2014) 248), and written informed consent was obtained from all eligible participants at recruitment.
DNA sample extraction and whole-genome sequencing
Genomic DNA was extracted from peripheral blood samples using the QIAamp DNA Mini Kit (Qiagen, 51306) following the manufacturer’s protocol. DNA concentration was measured using a Qubit 4.0 fluorometer (Thermo Fisher Scientific) for accurate quantification, and DNA integrity was assessed through 1% agarose gel electrophoresis under standard conditions. Only samples meeting high-quality criteria for both concentration and integrity were included in subsequent library preparation steps. Library preparation was carried out using the KAPA HyperPlus Library Preparation Kit (Roche, KK8514) following the manufacturer’s instructions. Genomic DNA was first fragmented to an average size of approximately 350 bp and then subjected to end-repair and A-tailing before adapter ligation, which included dual-index barcodes for multiplexing. After adapter ligation, the libraries were enriched by PCR amplification and purified using AMPure XP beads (Beckman Coulter) to generate final libraries. Whole-genome sequencing was performed on Illumina NovaSeq 6000 platforms (Illumina, San Diego, CA, USA) to generate 150 bp paired-end reads, with a target coverage of approximately 30× per sample.
Alignment and variant calling
The FastQC package (
bioinformatics.babraham.ac.uk/projects/fastqc) was used to assess the quality-score distribution of sequencing reads. Read sequences were mapped to the human reference genome (GRCh37) using the Burrows-Wheeler Aligner (BWA-MEM v0.7.15-r1140
)42 with the default parameters. The quality metrics (including coverage, median insert size, and percentage of chimeric reads) and contamination estimation of “Bam” files were obtained by Picard (v1.70) (
broadinstitute.github.io/picard) and VerifyBamID2 (v1.04) (
github.com/Griffan/VerifyBamID), respectively.
Single-nucleotide variants (SNVs) and short insertion/deletion variants (indels) were jointly called using the Genome Analysis Toolkit (GATK v3.8.1)
43 following GATK Best Practices for germline SNVs and indels. Briefly, variant calling was performed on individual samples using HaplotypeCaller in gVCF mode with default reads filtering parameters, including local realignment. The individual gVCF files were then jointly genotyped to identify high confidence alleles using GenotypeGVCFs tool for all autosomes and the X chromosome. Further details on alignment and variant calling pipeline can be found in a previous publication
44.
Quality control of variants
We implemented a rigorous variant quality control (QC) process that combined hard filters and a Random Forest (RF) model
45. The initial hard filters excluded variants based on the following criteria: (1) excess of heterozygosity, indicated by an inbreeding coefficient of < –0.3; (2) spanning deletion; (3) base changes exceeding 50 bp; and (4) call rates below 95%. The RF model was optimized and trained using the Python auto-sklearn package
46 to distinguish true variants from potential artifacts. This model operated at the allele level and encompassed both SNVs and indels. A total of 61 features were included in the model, categorized as follows: (1) site-level annotations from GATK HaplotypeCaller (e.g., genotype quality, inbreeding coefficient, and ReadPosRankSum); (2) allele-level annotations derived from in-house scripts (e.g., variant allele fraction, allele-level quality by depth); and (3) region-specific annotations (e.g., GC content, presence in the simple repeats sequences). The training sets for the RF model included both positive and negative examples constructed from whole-genome sequencing data. The positive training set included alleles previously genotyped or confidently discovered in three databases: Omni 2.5, Mills, and 1000 Genomes Project. To better identify rare genetic variants, singletons from unrelated samples that were transmitted in one of the trios were also included. The negative training set consisted of alleles that failed traditional GATK hard filters (QD < 2 or FS > 60 or MQ < 30). Additionally, chromosome 20 was excluded for further evaluation of its model performance.
Untargeted UPLC-MS/MS assay and quality control
Fasting blood samples were collected from participants in the morning after an overnight fast of at least 8 h at pre-pregnancy and all three gestational trimesters. This standardized fasting protocol was used to reduce dietary-related fluctuations in circulating metabolites and improve the reliability of metabolic phenotyping and genetic association analyses. Plasma and blood cells were separated and stored at –80 °C. A blind randomization approach was implemented using participant IDs, ensuring that samples from the same participant (at different gestational ages) were processed in the same batch.
Each plasma sample was extracted using a methanol-based extraction solution and shaken vigorously for two minutes. Proteins were denatured and removed by centrifugation, and the resulting supernatant containing the extracted metabolites was divided into four fractions for distinct LC-MS assays, including two reverse-phase LC-MS methods (different chromatography conditions) under positive electrospray ionization (ESI) mode, one reverse-phase LC-MS method under negative ESI mode, and one hydrophilic interaction liquid chromatography (HILIC) LC-MS method under negative ESI mode. The LC-MS instruments were the ACQUITY 2D UPLC system (Waters, Milford, MA, USA) coupled with the Q Exactive Orbitrap mass spectrometer (Thermo Fisher Scientific, San Jose, USA). To maximize analytical throughput, a two-column parallel configuration was used. During the analysis of one column, the other column was subjected to high-organic strong washing followed by full re-equilibration to initial gradient conditions.
The reverse-phase C18 column was BEH C18, 2.1 × 100 mm, 1.7 μm (Waters, Milford, MA, USA), and the HILIC column was BEH Amide, 2.1 × 150 mm, 1.7 μm (Waters, Milford, MA, USA). The mobile solutions for the first positive ESI LC-MS method were water and methanol containing 0.05% perfluoropentanoic acid (PFPA) and 0.1% FA, with a gradient of the methanol mobile phase linearly increasing from 5% to 95% within 3.4 min, and the flow rate was 0.35 mL/min. The column was then washed with high organic solvent and re-equilibrated to initial conditions for subsequent injection. The mobile solutions for the second positive ESI LC-MS method were optimized for more hydrophobic metabolites: mobile phase A was water with 0.05% PFPA and 0.1% FA; mobile phase B was methanol/acetonitrile/water (50:45:5, v/v/v) with 0.05% PFPA and 0.1% FA. The gradient program for acetonitrile mobile phase was initiated at 40%, then linearly increased to 99% in 1 min, held for 2.4 min, resulting in a total run time of 3.4 min with a flow rate of 0.5 mL/min. After each analysis, the column was washed and re-equilibrated to the starting conditions. For the negative ESI C18 LC-MS method, the mobile phase consisted of water and methanol containing 6.5 mM ammonium bicarbonate. The gradient for methanol increased linearly from 1% to 99% over 5.5 min, with a flow rate of 0.35 mL/min. For the negative ESI HILIC LC-MS method, the mobile solutions were water and acetonitrile containing 10 mM ammonium formate. The gradient for acetonitrile mobile solution started at 80%, decreased to 20% within 5 mins, and then linearly decreased to 5% in the final 0.5 min, resulting in a total run time of 5.5 min with a flow rate of 0.5 mL/min. After each run, the column was re-equilibrated to ensure reproducible retention times.
The mass spectrometer analysis was conducted using a QE mass spectrometer operated at 35,000 mass scan resolution. The instrument alternated between MS and data-dependent MS2 scans, using dynamic exclusion. The scan range was 70–1,000 m/z. Key operation conditions included a capillary temperature of 350 °C, a sheath gas flow rate of 40, and aux gas flow rate of 5. Given the high column efficiency of UPLC (2–3 times higher than conventional HPLC), the relatively short run times did not compromise chromatographic separation or increase co-elution. High mass accuracy (< 10 ppm) of the Q Exactive Orbitrap mass spectrometer further discriminated co-eluting metabolites with different m/z values. In addition, isomeric metabolites (e.g., isoleucine/leucine, mannose/glucose) were chromatographically separated where possible to ensure reliable annotation.
Samples from NC and ART women were analyzed in two separate batches. Pooled quality control (PQC) sample, generated by combining a small volume of all experimental samples, was included as a technical replicate and analyzed every 10 experimental samples. Process blanks consisting of extracted water samples were also analyzed. Internal standards were added to every sample to monitor instrument performance and facilitate chromatographic alignment.
Metabolite identification and quantification
Raw data processing, peak detection and annotation were performed using an in-house software
47. To ensure high-quality data, peaks were excluded if they met the following criteria: apex intensity less than 3,000, fewer than 7 scan points, retention time width less than 0.02 min, and a signal-to-noise ratio < 5. Peaks approved by the software, along with their associated spectra underwent further manual inspection. Metabolite peak areas were calculated using area-under-the-curve method. Metabolites were identified by comparing experimental data with an in-house reference standard library containing thousands of entries generated from purified metabolite standards analyzed through the same experimental LC-MS platform. This approach ensured that metabolite identification adhered to level 1 standards as defined by the Chemical Analysis Working Group of the Metabolomics Standards Initiative (MSI)
48, 49. The identification process was based on three criteria: narrow window of retention time (less than 0.1 min during one runday), accurate mass with variation less than 10 ppm, and MS/MS spectra with above 75% forward and reverse matching scores of the experimental spectrum to the reference entries in the library. Metabolites identified using public databases or previous publications by matching MS and MS/MS data but not validated with reference standards on the current platform were denoted with asterisks (not MSI level 1 identification).
Metabolomics data processing
The samples from NC and ART women were assayed separately in two batches, resulting in the identification of 869 metabolites across the two populations. To minimize batch effects, the samples were fully randomized and divided into several run-day blocks. The raw peak areas were corrected in run-day blocks by aligning the medians across experimental samples to equal one. Data points were normalized proportionately to ensure comparability
50. Only metabolites with a detection rate ≥ 50% across the total study population, at least 3-fold higher than the corresponding blank signals, and RSD ≤ 30% across pooled QC samples were retained for GWAS analyses, yielding a final panel of 748 eligible metabolites. Metabolite levels were then log
2-transformed, subjected to outlier trimming at ± 5 standard deviations from the mean, and standardized to a mean of 0 and standard deviation of 1. To assess changes in metabolite levels during pregnancy, trimester-to-trimester ratios were computed for each metabolite using batch-normalized values across trimesters within the same individual. The metabolite ratios were then transformed, trimmed, and standardized using the same procedures for the raw metabolite data.
Statistical analysis
Genome-wide association study of common and low frequency variants
After data processing and quality control, the variants that passed QC were further filtered using the following to remove those with an MAF < 0.01 and Hardy-Weinberg equilibrium (HWE)
P value < 1.00 × 10
–6 in single variant analysis. GWASs were performed within both the NC and ART populations. Linear regression of the metabolites and metabolite ratios were performed, with adjustment for age, GW of sampling, and the first five genetic principal components (fastGWA tool from GCTA version 1.94.1)
51. Associations identified in the NC and ART women were combined using inverse variance metaanalysis based on effect size estimates, and the heterogeneity between datasets was tested by Cochran’s Q test (which is equivalent to the McNemar test here as the number of data sets is two). All above analyses were carried out using Metal software
52.
Genomic inflation
The genomic inflation factor for each GWAS result was calculated to assess potential population stratification or systematic biases. This was done by dividing the median of the observed Chi-squared test statistics by the median of the expected Chi-squared test statistics for each metabolite and its trimester-to-trimester change.
Independent genome-wide significant associations
Independent genome-wide associations for each metabolite were identified using the LD-based clumping procedure implemented in PLINK 1.9
53. This approach was designed to filter out variants in LD while retaining statistically significant variants. Variants with a
P value < 5.56 × 10
–10 (Bonferroni-corrected for multiple testing of different trimesters and different populations) from GWAS conducted in the NC population, ART population, or from the GWAS meta-analysis were considered, and then clumps were formed around these 'index' variants based on LD threshold. The PLINK1.9 parameters used were: --clump-p1 5.56 × 10
-10, --clump-r2 0.1, and --clump-kb 500.
Identification of secondary associations
To uncover secondary signals within metabolite-specific region with only one signal, conditional analysis was conducted using GCTA v1.91.4. This approach identifies additional variant-metabolite associations after accounting for the primary signals. Briefly, the analysis began by conditioning on the most strongly associated regional variant identified in the marginal analyses, and the association of each additional regional variant was estimated independently within this conditional model. The regional variant with the lowest P value from the conditional analysis was then added to the model. All other regional variants were re-estimated using the updated conditional model. This process was repeated iteratively until no additional regional variants within the region met the significance threshold (P < 1.0 × 10–6).
Variant functional annotation
All variants were annotated using the VEP
15 (version 108; VEP documentation). The annotation process incorporated the “-pick_order” option to assign each variant using a single transcript, prioritizing transcripts in the following order: transcript support level, transcript biotype, APPRIS isoform annotation, deleteriousness of annotation as estimated by Ensembl, transcript CCDS status, canonical status of the transcript, and transcript length. Additionally, to enhance biological interpretation, we expand variant annotations to include potential regulatory elements of our findings. Element-gene predictions were obtained from Genehancer predictions from UCSC table browser (genome.ucsc.edu/cgi-bin/hgTables, table geneHancerRegElements, build hg19)
54, and from engreitzlab.org/resources/ (all element-gene connections with ABC scores ≥ 0.015)
55. Ultimately, a total of 207 variants were predicted to reside within regulatory elements, and these variants were linked to 1,648 element-gene pairs (details in Supplementary Table S3).
Putative causal gene nomination
To identify putative causal genes for independent genome-wide significant associations, a multi-step evaluation was conducted, leveraging genomic, transcriptomic and biological data. To identify genes with metabolite-associated loci, we first retrieved protein-coding genes located within or overlapping the 1 Mb region of metabolite-associated loci from the human GENCODE resource (www.gencode-genes.org/) using bedtools
56. These genes were then evaluated based on the following criteria:
(1) Biological relevance to associated metabolites
To determine whether the genes within the 1 Mb region were involved in biological processes related to the associated metabolites, we investigated their roles in enzymatic reactions, transportation, and other relevant processes. This evaluation was based on data from the HMDB
14, KEGG pathway database
57, and PubChem Chemical Co-occurrences in Literature database
58. This step allowed us to identify genes with biological relevance to the metabolites.
(2) Differential expression during pregnancy progression
Given that the study focuses on a pregnant population, we assessed these genes for their differential expression throughout pregnancy. Differentially expressed genes were annotated based on findings from Knight
et al.
59.
(3) Transcriptional regulation by associated variants
To evaluate whether the expressions of these genes were influenced by the associated variants, we first identified genes within or overlapping with metabolite-associated loci whose expressions were affected by these variants and their high linkage disequilibrium (
r2 > 0.8) counterparts. This evaluation was conducted by querying multi-tissue gene expression data from the GTEx project
60 using v7.signif_pairs.txt files. We included all variant-gene pairs that met statistical significance thresholds as determined by GTEx's permutation approach in any tissue. To further examine whether the same genetic variants were driving the associations with metabolites and eGenes, we performed a colocalization analysis of GWAS and eQTL signals using a stringent Bayesian approach implemented in the coloc R package (version 5.2.3)
61. Variants within ± 500 kb of the lead SNP at each locus were included, and eGenes within this region were analyzed. A PPH4 > 0.8 was used as the cutoff for colocalization. Additionally, we investigated whether these genes were targets of regulatory elements associated with the identified variants. This step helped us pinpoint expression-related genes.
Genes meeting both expression-related and biological relevance criteria were classified as first-tier causal genes, whereas those deemed relevant solely based on biological evidence were classified as second-tier candidates. In total, we identified first-tier causal genes for 470 metabolite-variant associations and second-tier candidate genes for 469 metabolite-variant associations.
Detection of novel associations and novel loci
To assess the novelty of the variant-metabolite associations identified in our study, we compiled a summary of 12 published mGWASs that utilized Metabolon's metabolomics profiling
2, 16-24 (see Supplementary Table S7). Collectively, these studies reported 8,500 unique significant genetic associations involving 1,737 metabolites (including both identified and unidentified metabolites). Before identifying novel variant-metabolite associations, we first examined the replication of signals reported in previous studies using our dataset. After confirming the robustness of our dataset through these replicated associations, we defined novel associations for each metabolite as those having no prior known associations near the index variant, determined using the following PLINK 1.9 parameters: --ld-window 1000000 and --ld-window-r2 0.01. To assess the novelty of loci, we evaluated whether any genetic variants within a locus region were associated with the same metabolites or different metabolites. A locus was defined as “known” if it contained genetic variants associated with the same or other metabolites reported in earlier studies. Conversely, loci that did not meet this criterion were classified as “novel”.
Identification of ancestry-specific associations
We further investigated the 321 novel associations in the pre-pregnancy period by conducting additional mGWASs in 844 participants, using metabolic profiles measured from pre-conception plasma samples. Furthermore, we compared the MAFs of the 279 index genetic variants involved in these 321 associations between East Asian populations and other ancestry groups, using data from 1KGP3. Associations were classified as ancestry-specific if the variants were either exclusive to East Asian populations or had MAFs in East Asians at least 10-fold higher than those in other ancestries. Additionally, to qualify as ancestry-specific, the associations had to meet a stringent statistical significance (P < 0.05/321 = 1.55 × 10–4) in the pre-pregnancy analysis and exhibit no significant heterogeneity in effects during pregnancy.
Identification of pregnancy context-dependent associations
For the 321 novel metabolite-variant associations not previously reported in non-pregnant populations, we applied a single model incorporating interactions between genetic variants and gestational stage to evaluate whether the relationship between genetic variants and metabolites changes across pre-pregnancy period and different stages of pregnancy. Families with data available for both pre-pregnancy and pregnancy were included in this analysis. Given the limited sample size, we set a significance threshold for interaction at P < 0.1. A total of 174 metabolite-variant associations met this threshold. Among these, 65 associations showed no heterogeneity between natural and assisted pregnancies and were not significant in the pre-pregnancy analysis (P > 0.05/321 = 1.55 × 10–4).
To further validate these pregnancy-specific or intensified associations, we conducted heterogeneity tests. First, an mGWAS was performed across all pre-pregnancy participants to compare the effects observed pre-pregnancy with those during pregnancy (meta-analyzed across natural and assisted pregnancies). Associations that were not significant pre-pregnancy and exhibited significant heterogeneity between pre-pregnant and pregnant women were selected. Among these, 14 associations that achieved a power of 0.9 and were significant in the interaction analysis were identified as candidate pregnancy context-dependent associations.
Logistic regression analysis of 5-methyluridine levels and biochemical pregnancy loss risk
To investigate the relationship between 5-methyluridine levels and early pregnancy outcomes, we conducted a nested case-control study within the JBC. A total of 201 cases of biochemical pregnancy loss and 781 controls who delivered live births (confirmed by hospital records) were included, all of whom had undergone metabolic profiling. Plasma 5-methyluridine levels were quantified using the Metabolon HD4 platform and batch-normalized with other metabolites. Logistic regression models were fitted to evaluate the associations between standardized plasma 5-methyluridine levels and biochemical pregnancy loss, with adjustment for age, pre-pregnancy BMI, ovulation induction protocols, and cycle type.
Validation of signals identified by GWASs of trimester-to-trimester ratios
To further validate potential genetic determinants of metabolite dynamics across gestation, we employed a mixed-effects regression model incorporating interaction terms between genetic variants and gestational weeks. This approach allows us to assess the temporal effects of genetic variants during pregnancy. Women were treated as random effects in the model to account for individual variability. A statistically significant interaction term (P < 0.05) was considered indicative of a consistent genetic influence on metabolite variation throughout gestation.
Stratified analyses of obstetric conditions
To evaluate whether associations differ across various health-related conditions, we conducted stratified analyses based on age, BMI, GDM status, parity, and conception type. Heterogeneity tests were performed for all associations that were significant in at least one group (P < 5.56 × 10–10) and involved variants with MAF greater than 0.05. The heterogeneity P values were adjusted using FDR correction, with significant results identified at an FDR threshold of 0.05. For obstetric conditions with relatively low incidence, such as preterm birth and HDP, we conducted sensitivity analyses excluding participants with these conditions. This approach allowed us to assess whether these conditions influenced the genomic atlas of the metabolome across gestation.
We accessed participants’ health-related information through electronic medical records (EMRs) maintained at the maternal and childcare hospitals where they were recruited and received antenatal care. These EMRs provided comprehensive details on maternal complications during pregnancy, delivery outcomes, and other relevant clinical data. At GW 24–28, pregnant women underwent universal screening for GDM using the 75-g 2-h oral glucose tolerance test (OGTT) after a 12-h fast. Additionally, fasting glucose tests were performed during early (GW 10–14) and late (GW 30–34) pregnancy. GDM was diagnosed if one or more plasma glucose values during the OGTT met or exceeded the following thresholds: fasting blood glucose (FBG) 5.1–6.9 mM, 1-h plasma glucose (1-h PG) ≥ 10.0 mM, or 2-h plasma glucose (2-h PG) 8.5–11.0 mM. Alternatively, GDM was also diagnosed if an FBG of 5.1–6.9 mM was detected independently or during the OGTT, following the WHO 2014 criteria. Pregnancies with FBG ≥ 7.0 mM or 2-h PG ≥ 11.1 mM were classified as diabetes in pregnancy and excluded from GDM-related analyses.
HDP were defined to include chronic hypertension, gestational hypertension, and preeclampsia. Chronic hypertension was diagnosed as systolic blood pressure (SBP) ≥ 140 mmHg or diastolic blood pressure (DBP) ≥ 90 mmHg documented before 20 weeks of gestation or prior to pregnancy, with elevated readings confirmed on more than one occasion. Gestational hypertension was diagnosed if SBP ≥ 140 mmHg or DBP ≥ 90 mmHg occurred after 20 weeks of gestation in the absence of proteinuria or prior hypertension. Preeclampsia was defined as SBP ≥ 140 mmHg or DBP ≥ 90 mmHg after 20 weeks of gestation, accompanied by one or more of the following: proteinuria (≥ 0.3 g/24 h), urine protein-to-creatinine ratio ≥ 0.3, positive random urine protein test, or evidence of organ or system involvement (e.g., cardiovascular, pulmonary, hepatic, renal, hematologic, gastrointestinal, or neurological dysfunction) or placental-fetal compromise. Detailed birth records, including gestational duration and birth weight, were documented by trained healthcare professionals attending all deliveries. Gestational duration was estimated using ultrasonography performed by licensed practitioners, and all participants had complete birth records. Preterm birth was defined as a live birth occurring before 37 complete weeks of gestation.
Phenome-wide associations of metabolite levels
To conduct phenome-wide association studies (PheWAS) for metabolites, we imputed plasma metabolite levels in 7,385 pregnant women using in-house whole-genome sequencing data from the JBC. To reduce the potential impact of horizontal pleiotropy, we limited our analysis to genetic variants associated with fewer than five metabolites, including both primary and secondary signals. The strength of these genetic instruments was evaluated by calculating the proportion of phenotypic variance explained by the instruments for each metabolite. Based on these evaluations, we generated weighted summed scores representing the genetic load for metabolite levels. The weight was assigned according to the marginal effect of the variants. These genetic scores were computed for each pregnancy trimester and applied to metabolites with at least two associated variants and > 5% instrument-explained phenotypic variance in either the NC or ART cohort. These genetic scores served as exposure variables to test associations with seven pregnancy-related phenotypes: GDM, preeclampsia (PE), birthweight, gestational duration, preterm birth, LGA, and SGA. For GDM and PE, we restricted the analysis to predicted metabolite levels from the first trimester. To address the issue of multiple testing in our PheWAS analysis, we adjusted P values using FDR correction and reported significant results at an FDR threshold of 0.05.
Gene-centric rare variant analyses
To investigate rare variant associations, we applied an MAF filter of 1% and a missingness filter of 5%. Rare variant burden associations were tested on a gene-by-gene basis using the default burden test implemented in Regenie
62. The analysis focused on coding variants annotated as loss-of-function (LoF) variants, including stop gain, frameshift, and splice donor/acceptor variants, as identified by VEP. Following our previously established analysis strategy
44, we included a total of 4,847 genes with at least three LoF variants in their coding sequences. Covariates such as age, gestational weeks at sampling, and the first five genetic principal components were incorporated into the models to account for potential confounding. Analyses were performed separately for natural and assisted conception groups. These were subsequently combined using meta-analyses. Statistical significance was defined as meta-
P < 3.44 × 10
–6, applying Bonferroni correction for 14,541 tests (4,847 genes × 3 trimesters).
DATA AVAILABILITY
Metabolomics data have been deposited in MetaboLights and are publicly available as of the date of publication. Accession number: MTBLS10849. The summary data for variant-metabolite associations can be queried on the website: medomicscgp.com/mQTL_Query/.
The GWAS was performed using GCTA-fastGWA (v1.94.1)
51. PLINK v1.9
53 was used to identify LD-independent variants, and regional association plots were generated using LocusZoom
63. All other data analyses were performed using R (version 4.2.2), with the following R packages utilized for analysis and plotting: dplyr (1.1.4), data.table (1.15.0), tidyverse (2.0.0), dplyr (4.4.0), ggbeeswarm(0.7.2), ggplot2 (3.4.4), circlize (0.4.16)
64, ComplexHeatmap (2.15.4)
65, RColorBrewer (1.1.3), ggpubr (0.6.0) and ggbreak (0.1.2).
The Author(s) 2026. Published by Higher Education Press. This is an Open Access article distributed under the terms of the CC BY license (https://creativecommons.org/licenses/by/4.0/).
