Introduction

Miscarriage is one of the most prevalent complications of pregnancy, occurring in 15.3% of all clinically recognized pregnancies. Recurrent miscarriage (RM) is frustrating for clinicians and heartbreaking for patients. Approximately 1.9% (1.8−2.1%) and 0.7% (0.5−0.8%) of women experience two and three or more miscarriages, respectively.¹ Miscarriage is known to be associated with several etiological factors, including chromosomal abnormalities in couples or embryos, anatomical factors in the reproductive organs, endocrine abnormalities, endometrial dysfunction, infectious factors, and thrombophilia. Nevertheless, the etiology of this disease remains an enigma in up to 50% of cases² and requires in-depth etiopathogenesis studies.

Given that successful pregnancy involves multiple maternal physiological changes and the intricate interactions between the fetus and mother, researchers have concentrated on the relationship between miscarriage and chromosomal abnormalities or maternal-fetal interface imbalance. Decidualization refers to the optimal transformation of endometrial stromal cells (ESCs) after embryo implantation, which is vital for establishing pregnancy in humans. Throughout pregnancy, immune responses are carefully regulated to balance antimicrobial protection for both the mother and fetus with the necessity for immune tolerance to prevent fetal rejection.³ Accordingly, the maternal-fetal interface, represented by the uterine decidua, harbors diverse populations of immune cells capable of performing various functions throughout the gestational period.

Gut microbiota is largely accountable for balancing host defense and immune tolerance.⁴ Butyrate is a specific short-chain fatty acid (SCFA) that is readily generated by microbial fermentation from fiber-rich diets, which plays a crucial role in sustaining intestinal homeostasis and holds a significant immunomodulatory property.⁵ Reportedly, elevated systolic and diastolic blood pressures in early pregnancy are entailed by the gut microbiota and butyrate.⁶ A prior study revealed that administration of butyrate to pregnant mice enhanced the resistance of newborn mice to inflammation and bile duct injury and improved their survival.⁷ Nonetheless, it remains ambiguous whether butyrate mediates the immune environment at the maternal-fetal interface. Butyrate metabolism-related genes (BMRGs) have been unveiled as potential prognostic biomarkers for renal clear cell carcinoma.⁸ In addition, another integrative bioinformatics study further revealed three biomarkers strongly associated with butyrate metabolism, which provides new insights into understanding the role of butyrate metabolism in pathological mechanisms such as ulcerative colitis (UC) and Metabolic-associated steatohepatitis (MASH), and signals that butyrate metabolism may become a highly valuable diagnostic biomarker.⁹ However, while butyrate metabolism has demonstrated importance in a variety of diseases, there is a research gap regarding the potential biological role of BMRGs at the maternal-fetal interface in patients with recurrent miscarriages (RM).

Given the impact of butyrate regulation on barrier function and the immune response, we hypothesized that imbalanced BMRGs may be associated with immune dysfunction or abnormal decidual function at the maternal-fetal interface during pregnancy loss. Therefore, this study aimed to identify butyrate-metabolism-related biomarkers in patients with RM, analyze the biological pathways they are involved in, perform immune profiling, and predict their associated drugs to define their mechanisms in RM further.

Materials and Methods

Data Collection

The RM-related datasets (GSE165004, GSE111974, GSE73025, and GSE179996) were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/). The GSE165004 and GSE111974 datasets were derived from endometrial tissues. The training set was the GSE165004 dataset that contained 24 RM and 24 control samples.¹⁰ The GSE111974 dataset containing 24 RM and 24 control samples was used as the validation set.¹¹ miRNA sequencing data (miRNA-seq) were downloaded from GSE73025 and included five control samples and 5 RM samples. The samples of GSE73025 were derived from human chorionic villi. lncRNA-sequencing data (lncRNA-seq) were downloaded from GSE179996, which contained five control samples and five RM samples. And the samples of GSE179996 were derived from peripheral blood. Additionally, 393 butyrate metabolism-related gene sets were collected from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/index.jsp), and 11,726 butyrate metabolism-related genes (BMRGs) were obtained.

Screening of Differential Expression BMRGs (DE-BMRGs) and Enrichment Analysis

Initially, the mRNA expression levels between RM and control samples in the GSE165004 dataset were compared via “limma” package (v 3.52.4) with adj.P < 0.05 and |log2FC| > 0.5,¹² to acquire differentially expressed genes (DEGs). DE-BMRGs were acquired by taking the intersection of the DEGs and BMRGs. Subsequently, the “clusterProfiler” R package (v 4.4.4) was used to carry out enrichment analysis on DE-BMRGs, based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (P < 0.05).¹³

Screening of the Biomarkers

STRING was used to create a protein-protein interaction (PPI) network of DE-BMRGs (confidence score = 0.4). Module analysis of the PPI network was performed using MCODE software. Characteristic genes were identified using the least absolute shrinkage and selection operator (LASSO), support vector recursive feature elimination (SVM-RFE), and random forest (RF), respectively. The top 5 genes in terms of importance were identified as characteristic genes using RF. The characteristic genes acquired by the three algorithms were intersected to obtain biomarkers. In addition, the correlation between biomarkers was analyzed using Spearman’s algorithm. The functional similarity among the biomarkers was assessed via the “GOSemSim” R package. Finally, the diagnostic value of these biomarkers for RM was assessed using the receiver operating characteristic (ROC) curves in GSE165004 and GSE111974. The expression of biomarkers in the RM and control groups was analyzed using GSE165004 and GSE111974.

Generation of Nomogram

The nomogram was generated to predict the probability of developing RM using the “RMS” R package (v 6.0–1). Calibration curves, decision curve analysis (DCA), and ROC curves were used to examine the reliability of this nomogram.

GSEA of Biomarkers

To investigate the signaling pathway associated with these biomarkers, we accomplished Gene Set EnrichmentAnalysis (GSEA)via “clusterProfiler” R package (v 4.4.4) and org.Hs.eg.db. The correlation between the biomarkers and all genes in the RM samples was computed using Pearson’s test, and genes were ranked in accordance with their relevance. The threshold value was adjusted to P < 0.05.

Immune Analysis

The infiltrating abundance of 16 immune cells and the activity of 13 immune-related gene sets were assessed via single-sample gene set enrichment analysis (ssGSEA) (v 1.0.13)¹⁴ in the RM and control samples. Differences between RM and control samples were analyzed using the Wilcoxon test. The relationship between the differential immune cell/immune-related gene sets and biomarkers was evaluated using the Spearman algorithm.

Construction of Transcription Factors (TFs)-mRNA Network

We generated a TF mRNA network to further reveal the mechanisms involved in RM biomarkers. Biomarker TFs were predicted using CHEA3 (https://maayanlab.cloud/chea3/). The connectivity of the nodes in the network was calculated, and the network was constructed by selecting the top 100 nodes in order of connectivity.

Construction of the Competing Endogenous RNA (ceRNA) Network

Firstly, the differentially expressed miRNAs (DE-miRNAs), and differentially expressed lncRNAs (DE-lncRNAs) between RM and control samples were acquired via “limma” R package (v 3.52.4) (|log2FC| > 0.5 and P < 0.05). The miRNAs regulating these biomarkers were identified using the TarBase Database (http://mirtarbase.mbc.nctu.edu.tw/index.html). miRNAs were identified by considering the intersection of the predicted miRNAs and DE miRNAs. Subsequently, the LncBase was used to predict lncRNAs interacting with miRNAs (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex). The predicted lncRNAs were intersected with DE-lncRNAs to obtain the lncRNAs. Eventually, the ceRNA network was demonstrated using Cytoscape software (v. 3.9.0).¹⁵

Drug Prediction

Each biomarker was used as a keyword to search for drugs that interacted with biomarkers in the Drug-Gene Interaction database (DGIdb, https://dgidb.genome.wustl.edu/). The Cytoscape software (v. 3.9.0) was used to present the biomarker-drug network.

Expression Level Verification of Biomarkers

To verify the role of the biomarkers in RM, their expression levels were validated by reverse transcription quantitative polymerase chain reaction (RT-qPCR). This study collected 10 clinical tissue samples from 5 patients with unexplained RM and 5 controls from the Wenzhou Hospital of Integrated Traditional Chinese and Western Medicine. This study was approved by the Ethics Committee of Wenzhou Hospital of Integrated Traditional Chinese and Western Medicine (approval number: 2024-K001). Total RNA was prepared from tissue samples using TRIzol reagent. Reverse transcription was performed using the SureScript First-Strand cDNA Synthesis Kit to generate cDNAs. RT-qPCR was performed as follows: 40 cycles of 95°C for 1 min, 95°C for 20s, 55°C for 20s, and 72°C for 30s. GAPDH was used as an internal reference gene. RT-qPCR primers used are listed in Table 1. The relative expression levels of biomarkers were calculated using the 2^−ΔΔCT method.

Table 1 Primers for Biomarkers and Intrinsic Reference Genes

Statistical Analysis

Statistical analysis was performed using the R software. Differences between the groups were analyzed using the Wilcoxon test. Statistical significance was set at P < 0.05.

Results

Acquisition and Enrichment Analysis of DE-BMRGs

A total of 595 DEGs between the RM and control groups were identified, including 220 upregulated and 375 downregulated genes (Figure 1a and b). A total of 213 DE-BMRGs were obtained by taking the intersection of DEGs and BMRGs (Figure 1c). The DE-BMRGs were enriched in 607 GO-BP items, 57 GO-CC items, 51 GO-MF items (Table S1), and 20 KEGG pathways (Table S2), including homeostasis of the number of cells, viral process, response to parathyroid hormone, heterotrimeric G-protein complex, GTPase complex, focal adhesion, cell adhesion mediator activity, cell-cell adhesion mediator activity, alpha-mannosidase activity, long-term depression, Wnt signaling pathway and Fc gamma R-mediated phagocytosis (Figure 1d and e).

Figure 1 Identification of differentially expressed genes between RM and control groups. (a) The volcano map of DEGs. (b) The heatmap of DEGs. (c) The Venn diagram of 213 DE-BMRGs was obtained by overlapping DEGs and BMRGs. (d) GO enrichment analysis of DE-BMRGs. (e) KEGG enrichment analysis of DE-BMRGs. Abbreviations: RM, recurrent miscarriage; DEGs, differentially expressed genes; DE-BMRGs, differentially expressed BMRGs; BMRGs, butyrate metabolism-related genes; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes.

Acquisition of the Biomarkers

The PPI network of DE-BMRGs contained 156 nodes and 250 edges, indicating that GNG2 interacts with various proteins, including CRHR2, GIPR, and GNA11 (Figure 2a). One key module contained 11 genes: CFL1, ACTR2, TPM4, PFN1, ANXA2, IRF2, IFIT3, OAS1, IFI27, MX1, and IRF6 (Figure 2b). For LASSO, we identified eight characteristic genes: ACTR2, ANXA2, IFI27, IRF2, IRF6, OAS1, PFN1, and TPM4 (Figure 2c). Six characteristic genes were identified using the SVM-RFE algorithm: ACTR2, ANXA2, IFI27, PFN1, IRF6, and OAS1 (Figure 2d). For the RF algorithm, the top 5 genes in terms of importance were ACTR2, ANXA2, IFI2, PFN1, and OAS1 (Figure 2e). Subsequently, the intersection of feature genes from Lasso, SVM-RFE and RF machine learning screening was taken, and a total of 4 overlapping genes were obtained, namely ACTR2, ANXA2, PFN1 and OAS1, which were the biomarkers obtained from the screening (Figure 2f). Correlation analysis showed that all four biomarkers were significantly correlated (Figure 2g). PFN1 was negatively correlated with OAS1, ANXA2, and ACTR2, whereas OAS1, ANXA2, and ACTR2 were positively associated with each other. The functional similarity results demonstrated that the functions of all four biomarkers were similar (Figure 2h).

Figure 2 Identification of biomarkers for RM. (a) The protein-protein interaction network of DE-BMRGs. (b) Network diagram of key modules obtained via the MCODE plug-in. (c) Adjustment of feature selection in the LASSO model. (d) Plot of error rate versus number of features via SVM-RFE algorithm. (e) Presentation of the results of the RF analysis. A random forest decision tree diagram of 11 different DE-BMRGs and a lollipop diagram representing the importance of each gene. (f) Venn diagram for three algorithms screening biomarkers. (g) Results of correlation analysis of biomarkers. ** P < 0.01, *** P < 0.001. (h) Violin plots for functional similarity analysis of biomarkers. Abbreviations: RM, recurrent miscarriage; DE-BMRGs, differentially expressed BMRGs; MCODE, molecular complex detection; LASSO, least absolute shrinkage, and selection operator; SVM-RFE, support vector machine recursive feature elimination; RF, Random forest.

The AUC values of all biomarkers in the training set were 0.986 (ACTR2), 0.899 (ANXA2), 0.896 (OAS1), and 0.887 (PFN1) (Figure 3a). The AUC values of all biomarkers in the validation set were 0.97 (ACTR2), 0.757 (ANXA2), 0.785 (OAS1), and 0.552 (PFN1), respectively (Figure 3b). Moreover, the AUC value of the logistic model was 1 for both GSE165004 and GSE111974 datasets (Figure 3c and d). These results suggest that all four biomarkers have a diagnostic value for RM.

Figure 3 Diagnostic efficacy of the biomarkers in the prediction of RM. (a and b) The ROC curves estimate the diagnostic performance of the biomarkers (ACTR2, ANXA2, OAS1, PFN1) in the prediction of the RM in both GSE165004 (a) and GSE111974 (b) datasets. (c and d) ROC curves for the logistic model of biomarkers in both GSE113079 (c) and GSE166780 (d) datasets. Abbreviations: RM, recurrent miscarriage; ROC, receiver operating characteristic.

Construction and Assessing of the Nomogram

The nomogram is shown in Figure 4a. The calibration curve suggests that the predicted prevalence was very close to the actual prevalence (Figure 4b). DCA showed that the nomogram demonstrated the best net benefit (Figure 4c). The area under the curve of the nomogram was 1 (Figure 4d). These results indicated that the nomogram had a high predictive accuracy for the prevalence of RM.

Figure 4 Construction and validation of the RM diagnostic nomogram model. (a) A nomogram was used to predict the occurrence of RM. (b) The calibration curve assessed the predictive power of the nomogram model. (c) DCA curves were used to assess the clinical value of the nomogram model. (d) The ROC curve assessed the predictive accuracy of the nomogram. AUC, the area under the curve. Abbreviations: RM, recurrent miscarriage; DCA, decision curve analysis; ROC, receiver operating characteristic.

Enrichment Analysis of Biomarkers

To explore the biological functions of these biomarkers, we performed a GSEA. The top 5 GO items and KEGG pathways are shown. According to the GO results, ANXA2 and OAS1 were negatively linked to nuclear chromosomes and nuclear division (Figure 5a and b). ACTR2 was positively correlated with mRNA binding, ribonucleoprotein complex biogenesis, and cytosolic ribosomes (Figure 5c). PFN1 was actively correlated with the mitochondrial inner membrane, actin cytoskeleton, mitochondrial protein-containing complex, etc. (Figure 5d).

Figure 5 Gene set enrichment analysis of four biomarkers. GO terms (a-d) and KEGG pathways (e-h) enriched in ANXA2 (a and e), OAS1 (b and f), ACTR2 (c and g), and PFN1 (d and h). Abbreviations: GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes.

KEGG results showed that ANXA2 and OAS1 were positively linked to the MAPK signaling pathway, Rap1 signaling pathway, and NOD-like receptor signaling pathways (Figure 5e and f). ACTR2 was positively related to ribosomes, herpes simplex virus 1 infection, and oxidative phosphorylation (Figure 5g). PFN1 was negatively correlated with herpes simplex virus 1 infection and actively correlated with the pathways of neurodegeneration, multiple diseases, and diabetic cardiomyopathy (Figure 5h).

Immune Analysis Between RM and Control Groups

The RM and control groups differed in two immune cell and three immune-related gene sets, namely T helper cells, regulatory T cells (Treg), MHC class I, parainflammation, and type I IFN responses (Figure 6a). PFN1 was negatively correlated with Treg and type I IFN responses. ACTR2, ANXA2, and OAS1 were positively correlated with MHC class I, parainflammation, Treg, and type I IFN responses, while the correlation between ANXA2 and type I IFN responses was not significant (Figure 6b). In addition, PFN1 was significantly positively correlated with T helper cells. However, ACTR2, ANXA2, and OAS1 showed opposite correlations with T helper cells (Figure 6b).

Figure 6 Correlation between the biomarkers and immunity. (a) Comparison of the ssGSEA scores of the immune cells and immune-related gene sets between the RM and control groups. (b) Heatmap of correlations between biomarkers and differential immune cells/immune-related gene sets. Redder colors indicate more positive correlations, bluer colors indicate more negative correlations, and the magnitude of the value indicates the correlation coefficient. * P < 0.05, ** P < 0.01, *** P < 0.001, ns: P > 0.05. Abbreviations: RM, recurrent miscarriage; ssGSEA, single sample gene set enrichment analysis.

The TF-mRNA and lncRNA-miRNA-mRNA Networks

A total of 937 TF mRNA pairs were obtained, including AHR-ACTR2, AHR-ANXA2, AHR-OAS1, and ARID2-PFN1 (Table S3). The TF-mRNA network based on the top 100 nodes ranked in order of connectivity is shown in Figure 7a, including 100 nodes and 253 edges. For example, the TF mRNA pairs included MTERF2-ACTR2, NKX23-PFN1, STAT1-OAS1, and SP100-ANXA2. Subsequently, 129 DE-miRNAs and 199 DE-lncRNAs were identified in the RM and control groups, respectively. A total of 236 mRNA-miRNA pairs containing four mRNAs and 193 miRNAs were predicted. Finally, 13 miRNAs were identified (Figure 7b). A total of 5341 lncRNA-miRNA pairs were predicted, including 11 miRNAs and 4115 lncRNAs. Twelve lncRNAs were identified (Figure 7c). An lncRNA-miRNA-mRNA network was generated that contained two mRNAs, five miRNAs, and 12 lncRNAs (Figure 7d). For example, the regulated relationship pairs include PPP4R1L-hsa-miR-124-3p-ANXA2 and AC012073.1-AS1-hsa-miR-941-PFN1.

Figure 7 Construction of the regulatory networks and drug prediction. (a) TFs-mRNAs network based on the first 100 nodes in order of connectivity, red indicates mRNAs, and cyan indicates TFs. (b) Venn diagram of the intersection of miRNAs predicted by the TarBase database with differentially expressed miRNAs in the dataset. (c) Venn diagram of the intersection of lncRNAs predicted by the LncBase database with differentially expressed LncRNAs in the dataset. (d) Construction of the lncRNAs-miRNAs-hub regulation network, cyan, Orange, and red correspond to lncRNAs, miRNAs, and mRNAs, respectively. (e) Drugs prediction network diagram, red represents biomarkers, and yellow represents drugs. Abbreviations: TFs, transcription factors; mRNAs, messenger ribonucleic acids; miRNAs, microRNAs; lncRNAs, long non-coding RNAs.

Drug Prediction

Using DGIdb, we compiled a drug-gene interaction network. The results demonstrated that three drugs (withaferin A, N-ethylmaleimide, and etoposide) interacted with two biomarkers (ANXA2 and ACTR2) (Figure 7e). We found that withaferin A and N-ethylmaleimide are the target drugs of ANXA2. Etoposide is a target drug of ACTR2.

Validation of Biomarker Expression Levels

By visualizing the data from the training set, ACTR2, ANXA2, and OAS1 were found to be notably downregulated, whereas PFN1 was markedly overexpressed in the RM group (Figure 8a). According to the RT-qPCR results, the expression trends of ACTR2, ANXA2, OAS1 and PFN1 were consistent with the results in a public database (P<0.05) (Figure 8b). Between the RM group and the control group, there was no significant difference in the expression of ACTR2 (P>0.05); however, the expression of ACTR2 in the RM group was still slightly lower than that in the control group (Figure 8b). This might be due to the small sample size; further validation will be performed in the future.

Figure 8 Expression validation of four biomarkers. (a) Boxplot of expression levels of four biomarkers in the GSE165004 dataset. (b) Expression validation of four biomarkers between RM and control groups by RT-qPCR. * P < 0.05, **** P < 0.0001, ns: P > 0.05. Abbreviations: RM, recurrent miscarriage; RT-qPCR, quantitative reverse transcription polymerase chain reaction.

Discussion

Complex factors are involved in successful embryo implantation, as well as pregnancy establishment and maintenance. Only 50% of RM cases have a definite etiology, whereas the remaining cases are of unexplained etiology and may have a direct or indirect potential link to endometrial factors. Recent research has unraveled that there may be impaired decidualization prior to the onset of symptoms such as RM, preeclampsia, and placental abruption, indicating that disorders such as these are linked to abnormal decidualization.¹⁶ As early as 1984, attention was paid to the relationship between the gut state and pregnancy outcomes in women.¹⁷ Butyrate, one of the most important SCFAs, has sophisticated regulatory mechanisms as reflected by its multiple functions in a wide range of diseases. Nevertheless, few studies have analyzed butyrate metabolism in RM to date.

Our GO enrichment analysis showed that 213 BMRGs associated with RM were all enriched in cell-cell signaling by Wnt, focal adhesion, G protein-coupled receptor binding, and cell adhesion mediator activity. According to KEGG enrichment analysis, these were predominantly associated with the Wnt signaling pathway, the NOD-like receptor signaling pathway, etc. These are all closely related to uterine receptivity and decidualization and have been shown to play a vital role in reproductive tissues and early pregnancy events.^18,19

In the present study, three machine learning algorithms (LASSO, SVM-RFE, and RF) were utilized to further determine four butyrate metabolism-related biomarkers (ACTR2, ANXA2, PFN1, and OAS1) that could be used for the diagnosis of RM. Most of these biomarkers are linked to decidualization or the progression of RM.

Decidualization is a process requiring the reorganization of the actin cytoskeleton by motor proteins. ACTR2 is one of the two essential actin-associated proteins (ACTR2 and ACTR3) in the ARP2/3 complex,²⁰ which is a major actin nucleator and a central regulator of the actin cytoskeleton.²¹ Importantly, ACTR2 has been widely demonstrated to drive cell migration and invasion in various cancers.²² Given the similar invasive process between trophoblast and tumor cells, we speculated that ACTR2 downregulation may compromise normal decidualization and trophoblast cell migration. In a recently published study,²³ ACTR2 was significantly lowly expressed in RM patients, corroborating the results of our training set (GSE165004). Past research has pointed out that the decreased expression of ACTR2, a key hub gene in recurrent pregnancy loss (RPL), probably interferes with the cellular functions and signaling pathways closely related to ACTR2, which was related to the occurrence and development of RPL. More in-depth studies are expected to unravel the specific mechanism of ACTR2’s role in RPL. They may open up new avenues for RPL therapeutic strategies and provide potential points of targeted intervention.²⁴ As a pivotal member of the annexin family, ANXA2 is expressed on the surface of endothelial cells, macrophages, mononuclear cells, and various types of cancer cells. ANXA2 has a crucial role in many biological processes, such as endocytosis, exocytosis, autophagy, cell-cell communication, and biochemical plasminogen activation.²⁵ ANXA2 is highly expressed in endometria at mid- and late-secretory phases and is primarily distributed in the luminal epithelium,²⁶ suggesting that the endometrial epithelium is the primary epithelium of ANXA2. As reported, ANXA2 is an intracellular molecule abundantly expressed in the human receptive endometrium, which is critical for embryo attachment.²⁷ ANXA2 has been identified as a momentous marker in the human receptive endometrium.²⁸ In another study, ANXA2 was identified as having a key role in maintaining normal blood flow and preventing thrombosis, especially during pregnancy, when it is essential for maintaining the health of the placenta and fetal blood vessels. An immune response or dysfunction against ANXA2 may be important in RM and other obstetric complications.²⁹ In our study, both the training and validation sets showed decreased ANXA2 expression in the RM group. Endometrial defects in ANXA2 expression interfere with the decidualization of endometrial stromal cells, as well as the microenvironment that promotes embryo implantation and placentation in vivo,³⁰ suggesting that ANXA2 may provide a tool for novel diagnostic and/or therapeutic interventions in RM. The 2’-5’-oligoadenylate synthetases (OAS), including OAS1, OAS2, and OAS3, are components of interferon-induced genes recognized for their antiviral properties and have additional cellular functions, such as the induction of apoptosis, enhancement of IFN-αβ signaling, modulation of immune cell receptors, and autophagy.³¹ The functions of OAS1 include antiviral response, cellular growth, and differentiation control, ensuring successful implantation and pregnancy.³² In our training and validation sets, OAS1 expression was reduced in RM patients, illustrating the importance of OAS1 in uterine receptivity and embryonic survival. PFN1 is an actin-binding protein that contributes to the development of endothelial and vascular diseases,^33–35 which is now considered a novel predictor of hypertension and a potential pro-atherosclerotic factor. This gene activates the RhoA/ROCK signaling to produce ROS. Wu et al found that the PFN1/RhoA/ROCK pathway was tightly linked to endothelial damage in preeclampsia.³⁶ Several factors, including oxidant-induced endothelial damage, impaired placental vascularization, and immune dysfunction.³⁷ As a consequence, we postulate that abnormal upregulation of PFN1 may provoke placental dysfunction by promoting cytoskeletal rearrangement, increasing stress fibers, activating inflammation-related pathways, and triggering endothelial cell dysfunction and oxidative stress, thereby culminating in miscarriage and RPL. However, PFN1 expression was slightly increased in the validation set but was not statistically significant. This may have been caused by differences in the overall experimental design. Our validation process involved collecting five normal and five unexplained RM samples further to demonstrate the clinical diagnostic value of these 4 genes. These four genes were expressed in a similar manner to GSE165004, confirming the diagnostic properties of RM-related biomarkers and focusing on butyrate metabolism. It was determined that all biomarkers had potential diagnostic value in the management of RM, and a high correlation was found among these four biomarkers, with a positive correlation among ACTR2, ANXA2, and OAS1 and a negative correlation between PFN1 and other genes. These findings further support the notion that BMRGs linked to RM, as identified in this study, have promising potential as targets for experimental laboratory designs to elucidate the molecular mechanisms involved in RM advancement.

Currently, several scientific studies have unveiled the application value of biomarkers in the diagnosis, prediction, and treatment of RM. For instance, Xie et al screened four key biomarkers (PTPN6, GJA1, CPT1A, and CREB3L1) associated with oxidative stress and ferroptosis, establishing a solid theoretical foundation for the development of diagnostic and therapeutic strategies of RM.³⁸ The study by Firatligil et al demonstrated that CYR61 might be a potential biomarker for predicting the risk of recurrent pregnancy failure, which not only deepens the understanding of molecular mechanisms underpinning recurrent pregnancy failure but also unlocks novel avenues for the development of potential preventive strategies for this condition.³⁹ Additionally, the strong association between abnormal miRNA expression and RM has garnered widespread attention. For example, miRNAs circulating in peripheral blood, particularly miR-100-5p and let-7c, hold the promise of becoming novel non-invasive biomarkers and a robust tool for diagnosis and prognostic assessment of RM.⁴⁰ On this basis, our study achieved a breakthrough by identifying PFN1 and OAS1 as novel biomarkers for RM for the first time, which not only enriches the biomarker library in this domain but also offers a precious resource for exploring new therapeutic targets, further advancing the development of novel therapies for RM.

Subsequently, the biological functions of the four biomarkers were further analyzed with GSEA. KEGG enrichment analysis results displayed that the four biomarkers were highly enriched in the MAPK pathway, Rap1 pathway, NOD-like receptor pathway, and oxidative phosphorylation. MAPK is implicated in pathological pregnancy by modulating follicular implantation, trophoblast cell proliferation and invasion, immunomodulation, and apoptosis.^41,42 An abnormal MAPK signaling pathway attenuates trophoblast invasion and increases TNF-α levels. Overexpression of TNF-α promotes trophoblast apoptosis, stimulates the release of prostaglandin E2 (PGE2), stimulates uterine smooth muscle cells, and induces miscarriage.⁴³ In addition, the targeted deletion of c-Jun and JunB in AP-1, a target protein downstream of MAPK, is associated with early embryonic death.^44,45 The Rap1 pathway performs crucial biological functions of regulating cell adhesion, junctions, and proliferation.^46–48 Kusama et al observed that Rap1 was a crucial factor involved in endometrial decidualization in pregnant rats.⁴⁹ The exchange protein activated by the cAMP (Epac)/Rap1 signaling pathway is involved in cAMP-mediated decidualization of human endometrial stromal cells (ESCs).⁵⁰ Luo et al discovered that Rap1 could act as a key pathway that could coordinate with other pathways to mediate the proliferation, invasion, migration, and fusion of trophoblast cells.⁵¹ Oxidative phosphorylation is a principal catabolic and cellular energy-producing pathway in eukaryotic cells, which has been recently demonstrated to selectively coordinate tissue macrophage homeostasis.⁵² Decidua is the second-largest population of immune cells and macrophages,⁵³ which are essential for pregnancy establishment and maintenance because of their unique phenotypes and heterogeneity.

Conventionally, women are in a state of adaptive immunosuppression during pregnancy to avert a semi-allogeneic immune response to the fetus. Of note, there is an emerging view that fetal development from embryo implantation to birth is regulated by T helper (Th) 1/Th2 cells and their cytokines. Additionally, IL-17-producing Th (Th17) cells are elevated in patients with pregnancy complications including miscarriage,^54,55 with decreased and dysfunctional Treg cells.⁵⁶ Consequently, the Th1/Th2 model is insufficient to explain the mechanisms of maternal-fetal immune tolerance and needs to be extended to the Th1/Th2/Th17/Treg model.⁵⁶ Immune-mediated placental rejection results from inappropriate expression of trophoblasts and MHC class I antigens. The regulation of MHC class I genes is crucial for immunological acceptance in allogeneic conception⁵⁷ because cellular immune imbalance may lead to RM. Accordingly, an immune cell infiltration analysis was performed in our study to further ascertain the immunoregulatory mechanism of RM. The findings disclosed that RM patients presented with higher levels of Th cells and lower levels of Tregs, MHC class I, para-inflammation, and type I IFN responses. Further correlation analysis results revealed that ACTR2, ANXA2, and OAS1 levels were positively correlated with significantly different immune cells and that PFN1 levels were correlated strongly positively with Th cells and negatively with Treg and type I IFN responses. Overall, these four biomarkers may evoke the occurrence of RM by controlling Th1/Th2/Th17/Treg cells and immune imbalance at the mother-fetus interface.

Recently, the number of lncRNAs identified in the placenta has also increased. According to current research,^58–60 lncRNAs bind to miRNAs that regulate the proliferation, migration, and invasion of placental trophoblast cells through various signaling pathways and may also be involved in the development of RM and other diseases by promoting the inflammatory response, autophagy, and affecting glycolysis. However, studies on the mechanism underlying the role of lncRNA in RM are still at an early stage, and most molecular regulatory mechanisms are not yet clear. In this study, we constructed a lncRNA-miRNA-mRNA network. Finally, we identified 12 lncRNAs, five miRNAs (hsa-mir-183-5p, has-mir-124-3p, hsa-mir-361-3p, hsa-mir-941, has-let-7e-5p) and two mRNAs (ANXA2 and PFN1) that play important roles in RM. In the lncRNA-miRNA-mRNA network, hsa-mir-183-5p, hsa-mir-361-3p, hsa-mir-941, and has-let-7e-5p had a regulatory effect on PFN1, and has-mir-124-3p had a regulatory effect on ANXA2. To date, no direct studies have addressed the roles of these five miRNAs in RM. We found that eight lncRNAs served as ceRNAs to regulate the expression of PFN1 in RM by sponging has-let-7e-5p. LncRNAs 00841, AC012073.1, and Linc AC016876.1 can act as ceRNAs for hsa-miR-183-5p, hsa-miR-941, and hsa-miR-36-3p, respectively, which in turn regulate the expression of PFN1. Additionally, PPP4R1L can act as a ceRNA that regulates the expression of ANXA2 through the sponge hsa-miR-24-3p. During pregnancy, invasive migration of trophoblast cells into the uterine decidua, alteration of the uterine spiral arteries, and immunological immunity of the mother to the embryo are similar to tumorigenesis. Given that lncRNAs perform tumor prevention and provide therapeutic outcomes, lncRNA-micRNA-mRNAs associated with BMGEs may provide new insights into the pathogenesis of RM. To assess the potential of BMRGs as drug targets for treating RM, we constructed a drug-gene interaction network using DGIdb. These results indicated that the three drugs could interact with two biomarkers (ANXA2 and ACTR2), and a number of properties have been demonstrated for Withaferin A, including anticancer, antiangiogenic, anti-invasive, anti-inflammatory, and proapoptotic properties.⁶¹ Many studies suggest that women with pathological pregnancies, such as RM, exhibit signs of exaggerated inflammatory immune responses before and during pregnancy, as well as signs of tolerance to autoantigens and fetal antigens.⁶² Recent research has indicated that withaferin A may act as a selective glucocorticoid receptor modulator with anti-inflammatory effects,⁶³ suggesting a potential role in treating RM.

Nonetheless, this study also faced some limitations. First, the GSE111974 dataset was used as the validation set. Although GSE111974 is highly relevant to the field of RPL, it is strictly a sample set centered on cases of repeat implantation failure, which limits the representativeness of the dataset to some extent. Secondly, limited by the sample size, the findings may not present a comprehensive picture of the overall characteristics of patients with RM, and the small sample size may lead to insufficient statistical validity, increasing the uncertainty of the results. To address these limitations, we plan to expand the sample size in the future to weaken the sampling error and enhance the statistical validity, thereby elevating the representativeness and reliability of the results. Additionally, more in-depth clinical studies will be conducted to validate our findings, therefore more precisely revealing the role of the screened biomarkers in the pathogenesis of RM and providing a stronger foundation for the clinical diagnosis of RM and the development of therapeutic strategies.

Conclusion

In conclusion, we identified and validated four biomarkers using BMRGs that can be used as diagnostic markers for RM. These genes may influence RM progression through the regulation of decidualization.