The present study enrolled 192 participants, comprising 122 in the control group and 70 in the PL group. The participants were recruited from Lanzhou University Second Hospital between February 2023 and September 2023. The inclusion criteria for the PL group included women aged 18 to 42 years who had experienced at least one PL in the past six months, in accordance with the ESHRE diagnostic criteria (1). The exclusion criteria for the PL group were as follows: presence of endocrine disorders, infections, or immunological diseases; or a history of the most recent pregnancy ending in an ectopic pregnancy, hydatidiform mole, or congenital defects. The inclusion criteria for the control group were women with no history of PL and at least one successful full-term pregnancy. The exclusion criteria for the control group included the presence of endocrine disorders, infections, or immunological diseases. The present study was approved by the Institutional Review Board of Lanzhou University Second Hospital (ethical approval number: 2023A-553). All patients provided written informed consent.

Fasting blood samples were collected from participants using EDTA tubes and centrifuged at 1,500×g for 10 minutes at 4°C to obtain plasma, which was then stored at -80°C. For metabolomics analysis, 100 μL of plasma was placed in Eppendorf tubes, mixed with 80% prechilled methanol, and vortexed thoroughly. The samples were incubated on ice for 5 minutes, followed by centrifugation at 15,000×g for 20 minutes at 4°C. A portion of the supernatant was diluted with LC–MS-grade water to reach a 53% methanol concentration. This mixture was transferred to fresh Eppendorf tubes and centrifuged again at 15,000×g for 20 minutes at 4°C. The final supernatant was then injected into the LC–MS/MS system for analysis via an Orbitrap Q Exactive™ HF–X mass spectrometer (19, 20).

UHPLC–MS/MS analysis was performed using a Vanquish UHPLC system coupled with an Orbitrap Q Exactive™ HF-X mass spectrometer (Thermo Fisher, Germany) at Novogene Co., Ltd. Samples were injected onto a Hypersil Gold column (100 × 2.1 mm, 1.9 μm) and separated over a 12-minute linear gradient at 0.2 mL/min. The eluents for positive polarity mode consisted of 0.1% formic acid in water (eluent A) and methanol (eluent B), whereas those for negative polarity mode included 5 mM ammonium acetate (pH 9.0, eluent A) and methanol (eluent B). The solvent gradient was programmed as follows: initial 2% B for 1.5 minutes, increased to 85% B over 3 minutes, increased to 100% B over the next 10 minutes, decreased to 2% B over 0.1 minutes, and maintained at 2% B for 12 minutes. The mass spectrometer was operated in both positive and negative modes with a spray voltage of 3.5 kV, capillary temperature of 320°C, sheath gas flow rate of 35 psi, auxiliary gas flow rate of 10 L/min, S-lens RF level of 60, and auxiliary gas heater temperature of 350°C. High-energy collision dissociation (HCD) was used as the fragmentation method, with an isolation window scanning the m/z range of 100–1,500. Collision energy was applied in three steps, namely, 20 V, 40 V, and 60 V. The quality control (QC) samples were interspersed among the study samples during analysis. The QC samples were used to monitor and assess the data quality by evaluating the correlation between QC injections, ensuring the reliability and consistency of the analytical results. All samples were analyzed in a single batch during the analytical run.

The raw data files from UHPLC–MS/MS were processed using Compound Discoverer 3.3 for peak alignment, peak picking, and metabolite quantitation, following the protocols established in previous studies (21, 22). First, retention times and mass-to-charge ratios (m/z) were aligned across samples with a mass tolerance of 5 ppm, a signal intensity tolerance of 30%, and a defined intensity threshold. Metabolite identification was performed by matching accurate m/z values, adduct ions, and MS/MS fragmentation patterns against the mzCloud (https://www.mzcloud.org/), mzVault, and MassList databases. In databases containing fragmentation data, the acquired MS spectra were compared with reference fragment ions and collision energies, ensuring higher confidence annotations. Next, peak intensities were normalized to reduce run-to-run variations and enable consistent comparative analyses. Quality control (QC) samples were injected to monitor instrument stability, and only features with a coefficient of variation (CV) < 30% in the QC samples were retained. Finally, the normalized peak areas of the confidently identified metabolites were used for subsequent statistical analyses.

The identified metabolites were annotated using the KEGG and HMDB databases (https://www.metaboanalyst.ca/). Since our research focus is not on exogenous metabolites, we have excluded them from further analysis. A detailed list of these excluded exogenous metabolites is provided in Supplementary Table 1 . To visualize metabolic differences and identify potential metabolic signatures, orthogonal partial least squares discriminant analysis (OPLS-DA) was performed using the ‘ropls’ package in R software. Statistical significance was assessed using the t test, with differentially abundant metabolites selected on the basis of a variable importance in projection (VIP) score > 1 and a p value < 0.05. Heatmaps and volcano plots were generated using the ‘pheatmap’ and ‘ggplot2’ packages in R. Pathway analysis was conducted using the MetaboAnalyst platform (version 6.0), a widely used web-based tool for metabolomics data interpretation. The list of significantly altered metabolites was uploaded using HMDB (Human Metabolome Database) IDs, and these metabolites were assigned on the basis of MS/MS spectral matching and database searches. The KEGG database was used for pathway enrichment analysis. Enrichment analysis was performed using the hypergeometric test to identify overrepresented pathways compared with a reference metabolome; pathways with a p value < 0.05 were considered statistically significant. To assess the biological importance of metabolites within pathways, the relative-betweenness centrality algorithm was applied during pathway topology analysis. The results were visualized via a pathway impact plot, which integrates enrichment and topology analyses to highlight the most significantly impacted pathways.

LASSO regression was employed via the ‘glmnet’ package in R to identify significant metabolites associated with PL (23). The data were divided into a training set (60% of the data) and a validation set (40% of the data). Within the training set, 10-fold cross-validation was applied using the ‘caret’ package in R to optimize the model and avoid overfitting (24). The selected features from the LASSO model were then validated using the validation set. ROC curve, sensitivity, specificity, and confusion matrix analyses were conducted using the ‘pROC’ package in R. Additionally, ROC curve analysis on clinical data, including age, BMI, and WHR, was conducted to compare their predictive performance with that of the differentially metabolite abundance-based model.

Pearson correlation analysis was performed using the ‘Hmisc’ package in R to investigate the relationships between differentially abundant metabolites and clinical parameters, such as age, BMI, and WHR. Heatmaps were created with the ‘ComplexHeatmap’ package, and significant correlations were visualized using scatter plots generated with the ‘ggplot2’ package, offering insights into the metabolic alterations linked to clinical features in the control and PL group. Additionally, a permutation test with 200 iterations was employed to assess the statistical significance of the correlation coefficient differences between the control group and PL group.

Table 1 presents the baseline clinical characteristics of the study participants, which included 122 women in the control group and 70 women in the PL group. The median age in the PL group was significantly lower than that in the control group (30 years vs. 33 years, p < 0.001). The PL group had a higher median BMI (22.6 kg/m² [Q1-Q3: 21.5-24.7]) than the control group (21.8 kg/m² [Q1-Q3: 20.2-23.4], p = 0.006). Furthermore, the PL group presented a significantly greater WHR (0.9 [Q1-Q3: 0.8-0.9]) than the control group (0.8 [Q1-Q3: 0.8-0.8], p < 0.001). In terms of race, the majority of participants in both groups were Han Chinese, with no significant difference observed between the two groups (p = 0.726). However, education levels significantly differed between the two groups. A significantly greater proportion of women in the control group had a university education or above (91% vs. 67.1%, p < 0.001), whereas the PL group had a greater percentage of participants with a high school education or below (32.9% vs. 9%, p < 0.001). With respect to the number of pregnancy losses, none of the women in the control group experienced a loss, whereas the PL group included participants with varying numbers of losses, including 17.1% with one loss, 50% with two losses, and 32.9% with three or more losses (p < 0.001).

Figure 1 illustrates the OPLS-DA model and its evaluation through permutation tests. The OPLS-DA plot ( Figure 1A ) revealed a separation between the control and PL group, indicating distinct metabolic profiles. The x-axis (t1) explained 7.0% of the variance between the two groups, whereas the y-axis (to1) explained 8.0% of the variance within each group. A permutation test with 200 iterations ( Figure 1B ) was used to evaluate the OPLS-DA model. The histogram shows the frequency distribution of the permuted R²Y and Q² values, with an R²Y value of 0.932 and an Q² value of 0.899, which were both significantly higher than those obtained from the permuted models (p = 0.005), confirming the reliability and predictive performance of the model.

Figure 2 shows an in-depth analysis of the differentially abundant metabolites between the control and PL group. In total, 57 differentially abundant metabolites were identified on the basis of VIP> 1 and p values<0.05 ( Supplementary Table 2 ). The heatmap shown in Figure 2A displays the differentially abundant metabolites, revealing distinct clustering patterns between the two groups. In addition, the volcano plot shown in Supplementary Figure S1A illustrates the differentially abundant metabolites using log2-fold changes and -log10 p values. Testosterone glucuronide, 6-hydroxymelatonin, and (S)-leucic acid demonstrated statistical significance and fold changes, with testosterone glucuronide exhibiting the greatest significance and fold change. Among the 57 significantly altered metabolites, 54 were successfully matched to known metabolic pathways in the KEGG database using the MetaboAnalyst platform. The pathway impact plot revealed several significantly impacted pathways, including caffeine metabolism, tryptophan metabolism, and riboflavin metabolism ( Figure 2B ). The bubble plot shows these pathways, with larger and more vividly colored bubbles indicating pathways with greater impact and significance.

Figure 3 shows the results of LASSO regression and the subsequent analysis of selected metabolites. In Figure 3A , the coefficients of the metabolites are plotted against the log of the regularization parameter (λ), demonstrating the selection process for the most important metabolites in distinguishing between the control and PL group. Figure 3B shows the three most important metabolites selected by the LASSO model, namely, testosterone glucuronide, 6-hydroxymelatonin and (S)-leucic acid.

Figure 4A compares the abundance levels of these metabolites between the control and PL group. Significant differences in abundance were observed for all three metabolites, with p values indicating significant differences (p<0.0001). These findings suggested that these metabolites play crucial roles in the metabolic alterations associated with PL. An ROC curve of the LASSO model revealed that the combined AUC value was 0.993, indicating high predictive accuracy for the selected metabolites in the validation set ( Figure 4B ). The performance of the optimal metabolite model was evaluated using a confusion matrix ( Supplementary Figure S1B ) and classification performance indicators ( Supplementary Figure S1C ). The confusion matrix revealed that the model achieved a sensitivity of 1.0000 and a specificity of 0.9643, indicating high diagnostic accuracy. In addition, the model had an accuracy of 0.9868. The positive predictive value (PPV) was 0.9796, and the negative predictive value (NPV) was 1.0000, reflecting the robust ability of the model to distinguish between the control and PL group.

The diagnostic performance of the clinical parameters and metabolites was evaluated via receiver operating characteristic (ROC) curve analysis. Among the clinical parameters, WHR had the highest AUC (0.724), followed by age (AUC = 0.696) and BMI (AUC = 0.435). The combined clinical model achieved an AUC value of 0.775 ( Figure 5A ). Furthermore, the combined model, which integrated both key metabolites and clinical parameters, achieved an AUC value of 1.000 ( Figure 5B ).

Figure 6 shows the Pearson correlations between differentially abundant metabolites and clinical parameters within the PL group. The heatmap shown in Figure 6A depicts the Pearson correlation coefficients between core differentially abundant metabolites and clinical phenotypes. The scatter plot shown in Figure 6B further illustrates the negative correlation between WHR and testosterone glucuronide, which indicated that higher WHR values were associated with lower levels of testosterone glucuronide (r = -0.291, p = 0.0146), whereas (S)-leucic acid was significantly positively correlated with WHR (r = 0.248, p = 0.0381). Additionally, the correlation of WHR with the levels of testosterone glucuronide or (S)-leucic acid were evaluated in the control group. The scatter plot shown in Supplementary Figure S1D indicates that both correlations were nonsignificant in the control group (WHR vs. testosterone glucuronide: r = –0.0294, p = 0.748; WHR vs. (S)-leucic acid: r = 0.126, p = 0.168), suggesting no evident associations ( Supplementary Figures S1D, E ). Moreover, a permutation test comparing the correlation coefficients between the control and PL group revealed a statistically significant difference (p = 0.04) in the correlation between WHR and testosterone glucuronide, whereas the correlation between WHR and (S)-leucic acid did not differ significantly (p = 0.28).

Despite the promising findings, this study has several limitations. Firstly, the sample size is relatively small, which may limit the generalizability of the results and prevent detailed stratification analyses, such as differentiating metabolic variations among PL cases with different recurrence histories. Secondly, the study design is cross-sectional, which prevents us from establishing causal relationships between the identified metabolic alterations and PL. Although we applied OPLS-DA and LASSO regression to prioritize key metabolites and minimize irrelevant variation and multicollinearity, these methods cannot fully account for confounding factors such as age, BMI, and WHR. We employed a permutation test to further address this issue, but future studies incorporating matched or statistically adjusted designs are warranted. Additionally, we did not collect detailed dietary information or account for variations in education level, which may influence metabolic profiles. Finally, the study focuses solely on plasma metabolites and includes limited clinical data. Future studies should aim to integrate multi-omics approaches and more comprehensive clinical datasets to provide a deeper understanding of the molecular mechanisms underlying PL.