Integrative Multi-Omics Analysis Unveils Biomarkers Linking the Gut Microbiota, Blood Metabolites, and Recurrent Pregnancy Loss

Introduction

Recurrent pregnancy loss (RPL) is a problem that still troubles women of childbearing age to date. According to the newest European Society of Human Reproduction and Embryology (ESHRE) guidelines, it is defined as the loss of two or more consecutive pregnancies (excluding ectopic pregnancy and molar pregnancy) and affects approximately 1–2% of couples worldwide.¹ It has been confirmed that many factors can result in RPL, such as embryonic chromosomal abnormalities, autoimmune disorders, endometrial dysfunction and infections.^2,3 Nonetheless, the cause of more than half of RPL remains to be clarified.^2,3

The microbial community that colonizes within the human digestive system is known as the gut microbiota, and it is vital to human health.^4,5 Studies have demonstrated that patients with pregnancy loss exhibit significantly reduced gut microbial diversity, along with decreased relative abundance of Prevotellaceae and Selenomonas.⁶ In patients with RPL, those positive for antiphospholipid antibodies and antinuclear antibodies show higher gut microbial richness and diversity, as well as increased proportions of Megasphaera and Enterococcus.⁷ Thus, gut microbial dysbiosis—regardless of whether it presents as decreased or abnormally elevated diversity—adversely affects pregnancy outcomes and increases the risk of RPL. During pregnancy, the embryo acts as a semi-allograft, and the establishment and maintenance of maternal-fetal immune tolerance serve as the prerequisite for successful pregnancy.^8,9 In this process, the gut microbiota contributes to the maintenance of normal pregnancy by systemically regulating immune homeostasis and maternal-fetal immune tolerance^10,11 through a gut-placenta immune axis, mainly by inducing myeloid-derived suppressor cells (MDSCs) and gut-derived RORγt+ regulatory T cells (Tregs) to suppress excessive IFN-γ+ and IL-17A+ T cell responses at the maternal-fetal interface.⁹ Gut microbial dysbiosis disrupts immune balance, triggers a systemic pro-inflammatory state, and indirectly impairs embryo implantation, decidualization, angiogenesis, spiral artery remodeling, and placental development, thereby representing a crucial potential mechanism underlying RPL.¹²

Additionally, the embryo may alter maternal metabolic pathways through a series of complex regulatory mechanisms during pregnancy.^6,8,13 Notably, a metabolomic profiling has shown significant differences between RPL patients and the control group.¹⁴ It also remains unclear whether the gut microbiota is involved in and interacts with these metabolic alterations to affect RPL. Hence, clarifying the role of gut microbiota and blood metabolites in RPL is worth exploring. However, conducting clinical trials to address these questions is rather difficult due to the constraints of sample size and the interference of many confounders. By contrast, Mendelian Randomization (MR) has emerged as a powerful research methodology for assessing the potential associations between exposures and specific diseases as it can minimize the interference of confounders and the bias caused by reverse causality while providing reliable and robust estimations.^15,16

Herein, we explored the potential causal associations of gut microbiota and blood metabolites with RPL, and the potential interaction/mediating relationship between them through bidirectional two-sample MR and mediation analyses. Then, biomarkers linking the gut microbiota-blood metabolites network in RPL, as well as their roles in the RPL immune microenvironment were further identified via integrative bioinformatics analysis. This study aimed to improve the understanding of RPL pathogenesis and offer fresh perspectives and strategies for its clinical treatment.

Materials and Methods

MR Analysis Exploring the Associations of Gut Microbiota and Blood Metabolites with RPL

Study Design

This MR study was carried out in accordance with the STROBE-MR statement (Table S1).¹⁷ First and foremost, to eliminate bias and guarantee the reliability of causal inference, the genetic variants that were utilized as instrumental variables (IVs) in the MR study must comply with three fundamental assumptions: (1) there should be a significant association between IVs and the exposures under investigation; (2) there is no association between IVs and confounders of the exposures and outcomes; (3) the effect of IVs on the outcomes should be fully mediated by the exposures (Figure 1). Based on the above, the two-sample MR study was initially implemented to explore the potential causal associations of gut microbiota and blood metabolites with RPL. Subsequently, a two-step mediation MR analysis was implemented to further explore the effect of gut microbiota on RPL via blood metabolites, as well as that of blood metabolites on RPL via gut microbiota. Finally, reverse MR analysis was implemented to gauge the effect of RPL on both gut microbiota and blood metabolites. The steps involved and the results were accompanied by a series of rigorous statistical analyses and verifications. The flowchart of the overall study design was shown in Figure 2.

Figure 1 Fundamental assumptions and data of genetic variants for the bidirectional two-sample MR study. (A) MR study for gut microbiota and RPL; (B) MR study for blood metabolites and RPL. The black X symbol indicates the exclusion of potential confounding pathways, fulfilling the core exclusion assumptions of the MR study. The bold font represents important factors in the analysis: exposure and outcome. The label “bidirectional” indicates a bidirectional two-sample MR design, which simultaneously tests the potential causal effects in both directions (exposure factor → outcome and outcome → exposure factor) to separate the direction of causality. Abbreviations: MR, Mendelian randomization; RPL, recurrent pregnancy loss; GWAS, genome-wide association study.

Figure 2 Flowchart of the overall design for the bidirectional two-sample MR study and two-step mediation analysis. The dark black arrow indicates the sequential analysis process of the research, showing the direction of data processing and analysis from left to right. The bold font indicates important factors, analysis methods, and critical thresholds. Abbreviations: SNPs, single nucleotide polymorphisms; LD, linkage disequilibrium; IVW, inverse variance weighted; MR-PRESSO, MR pleiotropy residual sum and outlier; NA, no association.

Data Sources

Data regarding gut microbiota, blood metabolites, and RPL were obtained from the NHGRI-EBI GWAS Catalog database (https://www.ebi.ac.uk/gwas/). The human gut microbiota data (GCST90032172-GCST90032644) used in this study were derived from the largest GWAS published to date, which were based on fecal samples from 5959 Finnish individuals, involving 7,967,866 single nucleotide polymorphisms (SNPs) and covering 473 taxonomic units, including 11 phyla, 19 classes, 24 orders, 62 families, 146 genera, and 209 species.¹⁸ The metabolite data utilized in this analysis were from one of the most comprehensive metabolite studies, which encompassed 1091 blood metabolites and 309 metabolite ratios, deriving from the analysis of 8299 European participants in the Canadian Longitudinal Study on Aging cohort and involving approximately 15.4 million SNPs.¹⁹ The full GWAS summary statistics for these 1400 blood biomarkers were accessible to the public in the GWAS Catalog (GCST90199621-GCST90201020). The RPL dataset was identified as GCST011887 and consisted of 150,965 samples, including 750 RPL cases and 150,215 control samples, all of which were of European descent and encompassed 23,125,249 SNPs.²⁰

Selection of IVs

To fulfill the three assumptions of MR analysis (Figure 1), SNPs used as IVs need to adhere to the following criteria: Firstly, a relatively lenient statistical threshold of P < 1×10^−521,22 was employed to recognize SNPs associated with 473 gut microbiota and 1400 blood metabolites, to avoid excessive filtering and retain sufficient statistical power. For RPL, a threshold of P < 5×10⁻⁶ was utilized. Secondly, to avoid biased results caused by LD in this study, the “ieugwasr” package (v 1.0.0) was utilized to remove inapplicable IVs with LD, with a threshold of r² < 0.001 and a clustering distance of 10,000 kb. Thirdly, the F-statistic was applied to measure the strength of the IVs, those IVs exhibiting the F-statistic below 10 were deemed weak instruments and thus removed using the formula F = (N-k-1)/k × [R²/(1-R²)], where N represents the sample size of the exposure GWAS, k is the number of IVs, and R² is the proportion of variance explained.²³ Fourthly, the IVs strongly associated with the outcome GWAS traits (P < 1×10⁻⁵) were excluded to minimize potential confounding factors (Figure 2). Subsequently, to ensure the unidirectionality of causalities, the Steiger directionality test was used to further screen IVs.²⁴ Finally, the SNPs containing palindromic sequences were eliminated after matching the results to ensure the effects of SNPs on the exposure consistent with the same allele as those on the outcome.

Bidirectional Two-Sample MR Study and Mediation Analysis

To explore the impact of gut microbiota and blood metabolites on RPL, the TwoSampleMR package (v 0.6.0) and MRPRESSO package (v 1.0) were utilized.^25,26 The inverse variance weighted (IVW) method was selected as the primary method for MR analysis as it could integrate information from all IVs to estimate the causal effect through a weighted average. On this basis, to further confirm the robustness of the MR analysis, weighted median, MR Egger, weighted mode, and simple mode were employed to assess the potential causal effects.^27–29 Moreover, reverse MR analysis was also conducted, in which RPL served as the exposure while gut microbiota or blood metabolites as the outcome, to assess whether RPL exerted a causal effect on gut microbiota or blood metabolites.

To illustrate whether the mediator mediated the link between exposure and outcome, two-step mediation analyses were further executed. The effect estimates of gut microbiota (blood metabolites) on RPL (β1), the effect estimates of gut microbiota (blood metabolites) on blood metabolites (gut microbiota) (β2), and the effect estimates of blood metabolites (gut microbiota) on RPL (β3) were calculated. The mediation effect was quantified as (β2×β3), while the direct effect of exposure on outcome was represented as (β1-β2×β3). The mediation proportion was calculated as (β2×β3/β1) (Figure 2). Standard errors for the mediation effect, direct effect, and mediation proportion were estimated using the delta method.

Integrative Transcriptome Analysis of the Mechanisms Underlying Mediation Effects

Data Source

The training (GSE165004), validation (GSE26787), and scRNA-seq (GSE214607) datasets were acquired from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/). Dataset selection was based on the following criteria: the species was restricted to Homo sapiens, all samples were human tissue specimens, subjects were clearly divided into RPL patients and healthy controls, and all datasets were authentic and reliable. GSE165004 dataset (GPL16699) comprised endometrial samples from 24 RPL patients and 4 controls. GSE26787 (GPL570) dataset comprised 5 RPL and 5 control endometrial samples. GSE214607 dataset (GPL24676) comprised 6 RPL samples and 10 control decidual and villous samples. The transcriptomic data were obtained as preprocessed and standardized expression matrices. No further upstream processing or data merging was performed, and all analyses were conducted independently in each dataset. Additionally, 94 metabolic pathway-related genes (KEGG_GLUTATHIONE_METABOLISM and KEGG_LYSINE_DEGRADATION) were gathered via the MSigdb database (http://www.gsea-MSigdb.org/gsea/msigdb).

Differential Expression and Enrichment Analyses

In the GSE165004 dataset, genes corresponding to blood metabolites with mediating effects were subjected to the Wilcoxon test between the RPL and control groups, with genes having P-values less than 0.05 being recognized as key genes. Subsequently, the clusterProfiler package (v 4.6.2)³⁰ was implemented to execute Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses on these key genes (P < 0.05). Finally, the STRING database was employed to build protein-protein interaction (PPI) network, setting the low-confidence threshold at 0.150.

Identification of Biomarkers via Machine Learning

Candidate biomarkers were identified using three machine learning algorithms: least absolute shrinkage and selection operator (LASSO), Boruta, and support vector machine-recursive feature elimination (SVM-RFE), which were implemented using the glmnet (v 4.1–8),³¹ Boruta (v 8.0.0),³² and e1071 (v 1.7–14)³³ packages, respectively. LASSO excels in handling high-dimensional data and aids in variable selection and model interpretation. Boruta can select any feature associated with the response variable without lowering the cost function of the model. To optimize the margin between classes and facilitate efficient classification and regression, SVM-RFE finds the best hyperplane. Both Lasso and SVM-RFE analyses were performed with 10-fold cross-validation, while maxRuns = 500 was set in the Boruta analysis. Candidate biomarkers were determined by identifying overlapping genes from the results of these three algorithms. In the GSE165004 and GSE26787 datasets, genes that showed significant differences and consistent expression trends between the RPL and controls were designated as the final biomarkers.

Construction of the Nomogram

Based on the expression levels of biomarkers, the nomogram was established in the GSE165004 dataset via the rms package (v 6.7–0)³⁴ to predict RPL risk. The predictive accuracy of the nomogram was gauged by generating a receiver operating characteristic (ROC) curve, and its reliability was validated using a calibration curve.

Gene Set Enrichment Analysis (GSEA)

Within the GSE165004 dataset, Spearman analysis was executed between each biomarker and all genes using correlation coefficients as the ranking metric. Then GSEA was conducted using the clusterProfiler package (|NES| > 1 and P < 0.05), basing on the “KEGG.v2024.1.Hs.symbols.gmt” file downloaded from the MSigdb database.

Immune Infiltration Analysis

The CIBERSORT algorithm, in conjunction with the LM22 gene signature,³⁵ was implemented to analyze the proportions of 22 immune cells in RPL and controls from the GSE165004 dataset. Wilcoxon test was implemented to contrast variations in immune infiltration between the groups. Additionally, Pearson analysis was conducted to evaluate the relationships between biomarkers and immune cells.

Creation of the Regulatory Network

The TFs targeting the biomarkers were predicted by the ChEA3 database (https://maayanlab.cloud/chea3/), and those supported by ChIP-seq data from the ENCODE (https://www.encodeproject.org/) were selected. The microRNAs (miRNAs) regulating biomarkers were predicted in the miRDB (https://mirdb.org/) (score > 70). Subsequently, lncRNAs interacting with predicted miRNAs were predicted by the starbase database (https://rnasysu.com/encori/) (clipExpNum > 10). Finally, the networks of biomarker-TF and biomarker-miRNA-lncRNA were constructed using the Cytoscape software.

scRNA-Seq Analysis

The Seurat package (v 5.1.0)³⁶ was implemented to process the GSE214607 dataset. As part of the initial quality control, cells with fewer than 200 genes and genes present in fewer than three cells were excluded. Cells with mitochondrial gene ratios exceeding 25%, as well as those with gene numbers ≤ 200 and ≥ 8000, and total counts ≤ 200 and ≥ 100,000, were also removed. To identify the top 2000 highly variable genes, the vst method within the FindVariableFeatures function was applied. Data were standardized using ScaleData and subjected to principal component analysis with RunPCA, selecting the top 30 principal components for subsequent analysis. Cells were clustered using the FindNeighbors and FindClusters functions (resolution = 0.2), and the clustering outcomes were showed using UMAP. The expression of each cell type’s marker genes as documented in the literatures was used to annotate the various cell types.^37–41 DoubletFinder (v 2.0.4)⁴² was used to identify and remove potential doublets. The proportions of different cell types in RPL and controls were analyzed, and the distribution of biomarkers within these cells, as well as expression differences between RPL and controls, were explored. Finally, the CellChat package (v 1.6.1)⁴³ was employed to visualize cell communication networks, revealing ligand-receptor interactions and signaling patterns to elucidate communication mechanisms among various cells in the RPL microenvironment.

Statistical Analysis

To confirm the robustness of the causal effects, sensitivity test was executed, including heterogeneity test, horizontal pleiotropy test, and leave-one-out (LOO) test. The mr_heterogeneity function was used to inspect the heterogeneity level. When heterogeneity was present (P < 0.05), the IVW-random effects model was selected; otherwise, the IVW-fixed effects model was implemented. To gauge potential horizontal pleiotropy of the results, MR-Egger regression intercept tests and MR-PRESSO global tests were conducted. Any result demonstrating substantial pleiotropy (P < 0.05) was omitted from further analysis. In addition, the effect of each SNP on the results was measured by the LOO test. This involved sequentially excluding each SNP to evaluate its effect on the overall causal estimation.⁴⁴ The Wilcoxon test was employed to gauge group differences in bioinformatics analysis. The R software was employed for all statistical analyses. A statistically significant P-value was less than 0.05.

Results

Potential Associations of Gut Microbiota and Blood Metabolites with RPL

For the IVs corresponding to each exposure and outcome, the F-statistic ranged from 19.56 to 41.22 for gut microbial SNPs, from 19.52 to 999.88 for blood metabolites/metabolite ratios, and from 21.07 to 97.06 for RPL. All F-statistics exceeded the conventional threshold of 10, indicating the absence of weak instrument bias and ensuring the reliability of subsequent causal inference. Of the 473 gut microbiota, 28 gut microbiota were found to be significantly associated with RPL. Among them, the abundances of Brevibacillaceae, Faecalicatena sp001517425, Photobacterium, Staphylococcus aureus, Syntrophomonadia, Alloprevotella, etc were connected with an increasing risk of RPL [odds ratio (OR) > 1.000, P < 0.05]. The others such as Fibrobacterales, Gluconobacter, Pararhizobium, Ruminococcus, CAG-495, and UBA2658 sp002841545 were identified as potential protective factors for RPL (OR < 1.000, P < 0.05) (Figure S1).

Among 1091 blood metabolites and 309 metabolite ratios examined, suggestive causal associations with RPL were identified for 82 traits at a nominal threshold (P < 0.05). Among them, 37 blood metabolites and 12 metabolite ratios were identified as potential protective factors for RPL, such as levels of trimethylamine n-oxide, 3-hydroxyisobutyrate, phosphate, caprate (10:0), 12,13-DiHOME, and cysteine-glutathione disulfide (OR < 1.000, P < 0.05); while 27 blood metabolites and 6 metabolite ratios were recognized as risk factors for RPL, such as levels of 3-hydroxylaurate, 4-ethylphenylsulfate, 3-amino-2-piperidone, N-acetylvaline, X-23587, and cytidine to N-acetylneuraminate ratio (OR > 1.000, P < 0.05) (Figure S2). No heterogeneity or horizontal pleiotropy was observed in the findings (Table S2).

Potential Mediation Effect of Blood Metabolites in the Association Between Gut Microbiota and RPL

Based on the above results, an investigation was conducted to further explore the potential association between the 28 gut microbiota and 82 blood metabolites/metabolite ratios. A total of 114 significant associations were screened, comprising 60 potential risk factors and 54 potential protective factors (Figure S3). Subsequently, mediation analysis revealed that 3-amino-2-piperidone levels (GCST90200188) had a suggestive mediation effect (β = 0.017, P = 0.0478) on Photobacterium abundance (GCST90032511) and RPL with a mediation proportion of 14.4% (Figure 3A and Table 1). This suggestive mediation supports a potential mechanistic pathway linking Photobacterium abundance to RPL risk.

Table 1 Mediation Effect of 3-Amino-2-Piperidone on the Causality Between Photobacterium Abundance and RPL

Figure 3 The illustrative diagram of mediation analysis. (A) Mediation effect of 3-amino-2-piperidone on the causality between Photobacterium abundance and RPL. (B) Mediation effect of CAG-495 on the causality between cysteine-glutathione disulfide levels and RPL. Arrows indicate potential causal effect pathways: β2 (exposure → mediator), β3 (mediator → outcome), and β1 (exposure → outcome, total effect). The bold font indicates exposure factors, mediating variables, and outcomes.

Potential Mediation Effect of Gut Microbiota in the Association Between Blood Metabolites and RPL

The 82 blood metabolites/metabolite ratios were utilized as exposures, and 28 gut microbiota were designated as outcomes. A total of 120 nominally significant associations were identified, of which 63 showed potential protective associations with gut microbiota and 57 showed potential risk associations (Figure S4). Mediation analysis suggested that CAG-495 abundance (GCST90032290) exhibited a suggestive mediating role in the pathway between cysteine-glutathione disulfide levels (GCST90199784) and RPL (β = 0.003, P = 0.0497), with a mediation proportion of 15.5% (Figure 3B and Table 2).

Table 2 Mediation Effect of CAG-495 on the Causality Between Cysteine-Glutathione Disulfide Levels and RPL

Reverse MR Analysis

Reverse MR analysis was implemented to study the potential reverse associations between RPL and gut microbiota as well as blood metabolites. The analysis showed no reverse associations between them (Figures S5 and S6).

Transcriptome Analysis Identified Key Genes Linked to Metabolites in RPL

To further elucidate the roles of 3-amino-2-piperidone and cysteine-glutathione disulfide in RPL, integrative transcriptome analysis was conducted. Multiple families of secondary metabolites originate from basic amino acids, such as lysine and ornithine, as well as from intermediates or derivatives of their biosynthetic pathways.⁴⁵ Considering that 3-amino-2-piperidone is a product of ornithine metabolism and that ornithine shares high structural similarity with lysine, differing only by one fewer carbon atom in its side chain,^46,47 we investigated 94 genes associated with the glutathione metabolism and lysine degradation pathways. Among these, 20 genes showed significant differences between RPL and control samples and were identified as the key genes (Table S3). Enrichment analysis unveiled that these key genes were associated with 241 GO terms, including glutathione metabolic process, cellular modified amino acid metabolic process, and antioxidant activity (Figure S7A). KEGG analysis identified 21 pathways, such as glutathione metabolism, lysine degradation, and arginine and proline metabolism, further supporting the metabolic roles suggested by the mediation analysis at the transcriptomic level (Figure S7A). Additionally, the PPI network of these key genes uncovered 20 nodes and 78 interaction pairs, with ALDH9A1, G6PD, and GCLM showing more interactions with other proteins, suggesting their potential key roles (Figure S7B).

Determination and Functional Analysis of Biomarkers (ASH1L, G6PD, SETDB1, and LAP3)

Three machine learning algorithms were integrated to simplify the most critical feature variables. LASSO identified 10 feature genes (Lambda.min = 0.006) including G6PD, GCLM, GSTO2, GSTP1, etc (Figure 4A). Boruta analysis determined 11 of the most important feature genes, such as G6PD, GCLM, GSTO1, etc (Figure 4B). The SVM-RFE algorithm identified 16 genes when the error rate was lowest at 0.0672, including SETDB1, G6PD, GSTP1, etc (Figure 4C). By cross-referencing the results from the three algorithms, nine overlapping genes (G6PD, GCLM, GSTO2, GSTP1, LAP3, ASH1L, EHHADH, SETD1A, and SETDB1) were identified (Figure 4D). In GSE165004 and GSE26787, ASH1L, G6PD, and SETDB1 were significantly upregulated in RPL, while LAP3 was downregulated, showing consistent expression trends and serving as biomarkers (Figure 4E). Results from GSEA further indicated that ribosome and complement and coagulation cascades were commonly enriched among these biomarkers (Figure 5A–D and Table S4). Additionally, FC gamma R mediated phagocytosis and B-cell receptor signaling pathways were enriched in ASH1L, G6PD, and SETDB1, suggesting that these biomarkers may participate in the RPL process by regulating immune responses and cellular homeostasis.

Figure 4 Machine learning algorithms identified biomarkers for RPL. (A) Feature gene selection using the LASSO algorithm. (B) Feature gene selection using the Boruta algorithm. (C) Feature gene selection using the SVM-RFE algorithm. (D) Cross-referencing the results from the three algorithms yielded 9 candidate biomarkers. (E) Validation of biomarker expression in independent datasets. Genes labeled in red were upregulated in the RPL group, and genes labeled in blue were downregulated in the RPL group. *P < 0.05; **P < 0.01; ***P < 0.001. Abbreviations: LASSO, least absolute shrinkage and selection operator; SVM-RFE, support vector machine-recursive feature elimination.

Figure 5 GSEA of the biomarkers. (A-D) GSEA of ASH1L (A), G6PD (B), LAP3 (C), and SETDB1 (D), respectively. Abbreviation: GSEA, gene set enrichment analysis.

Development of the Nomogram Based on Biomarkers

The nomogram was developed based on the four biomarkers to predict the risk of RPL. Each feature variable in this nomogram was assigned a unique score, and the total of all feature scores in each sample indicated the likelihood of RPL occurrence (Figure 6A). ROC analysis demonstrated that the nomogram had an area under the curve (AUC) of 0.972, indicating good diagnostic value (Figure 6B). The calibration curve unveiled that the nomogram ensured a high consistency between the predictions and actual observations (Figure 6C). These findings suggested that the nomogram based on these biomarkers could serve as an effective tool for predicting RPL risk.

Figure 6 Development and validation of the nomogram for RPL risk prediction. (A) The nomogram was developed based on 4 biomarkers to predict the risk of RPL occurrence. The blue area shows the overall distribution of each variable. The gray dots are marked on the black median line, precisely indicating the specific value of the median. The red solid points represent the scores corresponding to the expression levels. The red quadrilateral and arrow at the bottom indicate the total score and its corresponding prediction probability. (B) ROC curve of the nomogram. (C) Calibration curve of the nomogram. Abbreviations: ROC, receiver operating characteristic; AUC, area under the curve.

Immune Cell Infiltration Patterns in RPL

To understand the immune microenvironment characteristics of RPL, CIBERSORT was implemented to analyze the infiltration patterns of 22 immune cell types (Figure 7A). Notably, regulatory T cells (Tregs) and M1 macrophages exhibited markedly variations between RPL and controls (Figure 7B). Correlation analysis further showed that SETDB1 had the strongest positive correlation with Tregs (cor = 0.473, P < 0.001), while ASH1L had the strongest negative correlation with M2 macrophages (cor = −0.457, P < 0.01), providing clues regarding the connection between biomarkers and immune regulation (Figure 7C).

Figure 7 Immune cell infiltration landscape and its correlation with biomarkers in RPL. (A) Composition of 22 immune cell types in RPL and control samples. (B) Differences in immune cell infiltration between RPL and control samples. Green represents the control group and Orange represents the RPL group. Red text indicates cell subsets with significant differences in proportion between the two groups. ns, no significance; **P < 0.01. (C) Correlation between biomarkers and immune cells. *P < 0.05; **P < 0.01; ***P < 0.001.

Regulatory Landscape of Biomarkers

To scout the regulatory mechanisms of the biomarkers, TFs were predicted via the ChEA3 database, and a biomarker-TF regulatory network comprising 89 nodes and 116 relationships was constructed (Figure S8A). GABPA was predicted to regulate G6PD, SETDB1, and LAP3; while CTCF, EBF1, and EGR1 were predicted to jointly regulate SETDB1 and LAP3. Additionally, 54 miRNAs corresponding to 4 biomarkers and 125 miRNA-lncRNA interactions were predicted, and a biomarker-miRNA-lncRNA network involving 4 biomarkers, 19 lncRNAs, and 54 miRNAs was constructed (Figure S8B). This network revealed relationships such as ASH1L-hsa-miR-139-5p-AC084082.1, LAP3-hsa-miR-1297-MALAT1, and SETDB1-hsa-miR-1296-5p-LZTS1-AS1, depicting the complex regulatory landscape of the biomarkers.

Single-Cell Level Analysis of RPL and Revelation of Complex Communication

scRNA-seq was performed to further scout the features of RPL at the single-cell level. After quality control, 29,299 genes and 125,973 cells were included (Figure S9A). Subsequently, all cells were classified into 20 clusters using UMAP, identifying 14 distinct cell types [B cells, T cells, dendritic cells, decidual macrophages (dM), mast cells, endothelial cells, epithelial cells, syncytiotrophoblast (SCT), villous cytotrophoblast (VCT), extravillous trophoblast (EVT), perivascular (PV), decidual stromal cells (DSC), decidual natural killer cells (dNK), and red blood cells (RBC)] (Figures S9B, C and 8A). Moreover, after removing 9448 (7.5%) high-confidence doublets, the cell boundaries became clearer, reducing the interference of technical noise with biological signals (Figure S9D). Among the 14 annotated cell types, dM and dNK had higher proportions in RPL (Figure 8B). Additionally, ASH1L was distributed across almost all cell types, while LAP3 was mainly found in dM (Figure 8C). Importantly, in dM, SCT, VCT, EVT, DSC, dNK, and RBC, significant expression differences of these biomarkers were observed between RPL and control groups (Figure 8D).

Figure 8 Single-cell analysis in RPL. (A) Annotation of cell types. (B) Proportion of cells in the RPL and control sample. (C) Distribution of biomarkers across cell types. (D) Expression of biomarkers across cell types in RPL vs. controls. ns, no significance; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Cell communication patterns unveiled that in comparison with the control group, interactions and interaction frequencies between dM and VCT increased in RPL (Figure 9A). Overall, RPL exhibited stronger signal input, output, and total communication patterns, which might lead to immune tolerance imbalance at the maternal-fetal interface (Figure 9B).

Figure 9 Cell communication analysis. (A) Cell communication across cell types. (B) RPL exhibited stronger signal input, output, and total communication patterns compared to the controls.

Discussion

Pregnancy is a long and arduous task for women, which involves a series of complex immune and metabolic regulatory mechanisms and is characterized by major shifts in maternal biology.¹³ Currently, studies regarding RPL, especially unexplained RPL, are more focused on the failure of maternal-fetal crosstalk caused by immune factors.^6,8,16 Here, from the perspectives of microbiota and metabolome, our MR analysis suggested potential links within the gut microbiota–blood metabolites axis in RPL and further identified four biomarkers (ASH1L, G6PD, SETDB1, and LAP3) at the transcriptomic level, providing an entero-metabolic axis perspective to understand the pathogenesis of RPL.

Short-chain fatty acids and other metabolites synthesized by gut microbiota can enter the bloodstream and regulate the host’s immune system together with endogenous metabolites,⁴⁸ while changes in metabolite abundance can also affect the composition and function of gut microbiota.¹³ This interaction may play a unique and more significant role in the occurrence of RPL. Surprisingly, our mediation analysis suggested that Photobacterium might putatively increase the risk of RPL through the mediation of 3-amino-2-piperidone (mediation proportion = 14.4%, P = 0.0478). Photobacterium is a Gram-negative facultative anaerobic coccobacillus that is widely distributed in marine habitats and can infect marine humans through the consumption of fish, causing various primary diseases.^49,50 3-amino-2-piperidone, an ornithine cycle-related metabolite which is probably involved in lysine metabolic pathways, has been found to promote the release of inflammatory factors.^47,51,52 Based on data from European populations, particularly fish-consuming Finns, our findings suggest that Photobacterium may enter the human body through diet and putatively influence RPL via immune regulation mediated by this metabolite. On the other hand, cysteine-glutathione disulfide, a non-endogenous metabolite formed under oxidative stress, has demonstrated protective properties in non-alcoholic fatty liver disease by influencing lipid metabolism.⁵³ Similarly, as an oxidized form of glutathione, it may putatively reduce RPL risk by lowering oxidative stress and inflammation,^54,55 although this protection was suggested to be partially weakened by the mediation of CAG-495 (mediation proportion = 15.5%, P = 0.0497). Studies have found that the synthesis of the disulfide is regulated by gut microbiota in mice with adenomyosis,⁵⁶ and CAG-495 was found to be increased during the remission phase of ulcerative colitis,⁵⁷ suggesting its putative regulatory role under certain pathological conditions.

To further explore the biological mechanisms of blood metabolites with mediating effects, we focused on the lysine degradation and glutathione metabolism pathways associated with these metabolites. Women with recurrent miscarriage have been reported to exhibit markedly reduced erythrocyte glutathione levels, particularly in those with autoimmune, unexplained, or luteal phase defect etiologies, indicating impaired antioxidant defense and elevated oxidative stress.⁵⁸ Glutathione depletion may promote lysine acetylation,⁵⁹ whereas lysine uptake enhances glutathione metabolism and reduces oxidative stress.⁶⁰ Accordingly, upregulation of lysine and significantly reduced levels of glutathione in placental tissues have been observed in the RPL group,^61,62 suggesting that lysine and glutathione metabolism probably play an indispensable role in RPL.

Subsequently, four biomarkers—ASH1L, G6PD, SETDB1, and LAP3 were identified. ASH1L is essential for controlling transcription and chromatin remodeling, as well as for promoting the methylation of certain histone lysine residues.⁶³ Recent studies have emphasized the pathogenic role of ASH1L in congenital malformations of the female genital tract.⁶⁴ In mouse models, ASH1L mutation leads to partial postnatal death, while the surviving mutant mice exhibit growth dysfunction and infertility due to defects in epididymis and uterus development.⁶⁵ In addition, ASH1L regulates the expression of p63 and p-CHK2 during early meiosis in mice, thereby protecting oocyte genome integrity and eliminating oocytes with severe DNA damage.⁶⁶ G6PD is the rate-limiting enzyme of the pentose phosphate pathway and is a key molecule for cells to resist oxidative damage.^67,68 Maternal G6PD deficiency has been proven to be fatal to embryos and causes severe placental abnormalities.⁶⁹ Additionally, children with G6PD deficiency showed increased oxidative damage to embryonic DNA, fetal mortality, and birth abnormalities when treated with the anticonvulsant medication phenytoin.⁷⁰ SETDB1 was originally thought to H3K9 in the nucleus, where it regulates chromatin function.⁷¹ It is crucial for maintaining embryonic stem cell pluripotency and inhibiting trophoblast over-differentiation.⁷² Maternal SETDB1 deficiency leads to preimplantation developmental arrest of embryos, accompanied by cell cycle and chromosome segregation abnormalities.⁷³ LAP3 is participated in the processing of bioactive peptides and the presentation of MHC-I antigens in mammals.⁷⁴ Its expression level significantly affects the development of sheep embryonic myoblasts.⁷⁵ In conclusion, these findings link specific biomarkers to metabolic disorders and RPL pathology, providing a gene-level understanding of the gut microbiota-metabolic axis mechanism of this disease.

GSEA revealed that the four biomarkers were commonly enriched in the complement and coagulation cascade pathways. Over-activation of this pathway has been reported to be related to RPL.⁷⁶ By expressing various regulatory proteins, human trophoblasts facilitate controlled complement activation, which is beneficial for spiral artery remodeling and the clearance of cell debris.⁷⁷ Crosstalk between the complement and coagulation cascades can prevent maternal rejection of the embryo and is favorable for maintaining normal pregnancy. However, excessive activation triggers an intrinsic immune feedback loop that simultaneously induces a compensatory anti-inflammatory response, which rapidly amplifies other targeted responses, thereby causing inflammation at the maternal-fetal interface or systemically and ultimately leading to miscarriage.⁷⁸ These results suggested that abnormal activation of the complement and coagulation cascades might significantly contribute to the pathogenesis of RPL.

Moreover, the interconnection between the coagulation system and immune system is particularly prominent under pathological conditions, further suggesting that immunity may be involved. During pregnancy, the maternal immune system undergoes significant changes to maintain maternal-fetal tolerance.⁷⁹ Our immune infiltration analysis unveiled marked variations in Tregs and M1 macrophages between RPL and controls. Notably, at the single-cell level, the dM showed more complex cell communication patterns in RPL. Abnormal polarization of dM is known to be associated with RPL,⁸⁰ and M1/M2 macrophage imbalance is considered one of the important causes of spontaneous abortion.⁸¹ M1 macrophage-dominated pro-inflammatory responses (such as high expression of TNF-α) are markedly connected with the occurrence of RPL.^82,83 On the other hand, Tregs multiply following exposure to fetal antigens and are crucial for sustaining maternal-fetal immunological tolerance.⁸⁴ They often proliferate at the maternal-fetal contact and in peripheral circulation during parturition.⁸⁵ Certain studies indicate that both elevated and diminished levels of Tregs may increase the incidence of miscarriage, exhibiting a U-shaped impact curve.⁸⁶ These findings highlight the immune imbalance in RPL. It is noteworthy that the gut microbiota and blood metabolites screened in this study may participate in the pathogenesis of RPL by modulating the immune microenvironment.

The advantages and innovations of this study are reflected in the following aspects. Firstly, the effect of microbiota and metabolism on human reproductive health is currently a widely-concerned topic, and to the best of our knowledge, this is the first MR study examining microbiota and metabolome to investigate potential causes of RPL. In this study, we identified dozens of gut microbiota and blood metabolites as risk or protective factors for RPL. Secondly, mediation analysis provided support for potential mediating effects involving gut microbiota and blood metabolites. Lastly, by combining mediation MR analysis with transcriptomics, we not only identified potential biomarkers but also clarified that these molecules might mediate their effects through immunomodulatory mechanisms, thereby establishing a comprehensive gut microbiota-metabolites-immune axis in RPL.

However, our study has certain limitations. Firstly, although we utilized the largest GWAS datasets that currently available, all the data were from individuals of European ancestry. This reduced heterogeneity and restricted the generalizability of our findings to other populations at the same time. Further validation in diverse groups is needed to confirm their broader generalizability. Secondly, due to the limited number of SNPs reaching genome-wide significance, we adopted a relatively lenient P-value threshold. Although this approach is well-established in similar studies^53,87 and was supplemented multiple validation methods to ensure reliability, further large-scale randomized controlled trials are still warranted to enhance the validity of the findings. Thirdly, as with all MR investigations, our study cannot completely exclude horizontal pleiotropy or weak instrument bias, which may affect the stability of causal estimates. Additionally, the nomogram had potential overfitting risk due to the small control sample size in the training set. The scarcity of public RPL datasets limits further sample expansion and external validation, requiring cautious interpretation and future large-scale independent validation. Finally, the specific mechanisms of the obtained biomarkers also need to be functionally verified through in vivo and in vitro experiments.

Conclusion

MR and mediation analyses unveiled the potential associations of gut microbiota and blood metabolites with RPL and two mediation pathways: Photobacterium-3-amino-2-piperidone-RPL and cysteine-glutathione disulfide-CAG-495-RPL. Integrative transcriptome analysis further identified four candidate key biomarkers involved in these pathways: ASH1L, G6PD, SETDB1, and LAP3. Additionally, our study revealed that these biomarkers and metabolic pathways converged to regulate immune homeostasis. Notably, this study is hypothesis-generating rather than definitive, and all findings should be considered preliminary and suggestive. In summary, our study identifies candidate microbial, metabolic, and candidate biomarkers with putative diagnostic and therapeutic potential for RPL, supporting the role of the gut microbiota–metabolites–immune axis. These hypothesis-generating findings may enhance the translational relevance of RPL research but require validation in future studies.