Molecular subtype of recurrent implantation failure reveals distinct endometrial etiology of female infertility.

To comprehensively characterize the endometrial transcriptome, we analyzed gene expression data from four publicly available RIF gene expression datasets and one independent in-house dataset, including three microarray datasets and two RNA sequencing datasets. Subsequently, three datasets with large sample sizes— GSE111974 , GSE58144 , and GSE106602 —were integrated to identify robust differentially expressed genes (DEGs). These datasets contained 10,383 common genes from 108 controls and 85 RIF samples, while the GSE71331 dataset served as the validation cohort. Additionally, an independent cohort of 33 endometrial biopsy samples—comprising 21 controls and 12 RIF cases—was prospectively collected, sequenced, and analyzed following predefined inclusion criteria to enable independent validation and clinical characterization of the study findings. Table  1 provides a comprehensive summary of the characteristics of all datasets utilized in this study, encompassing 134 controls and 104 RIF samples. Supplementary Table 3 provided detailed demographic information for each dataset. Table 1 Datasets used in identifying RIF subtypes GSE GPL Stage RIF Control GSE111974 GPL17077 (Agilent) mid-luteal 24 24 GSE58144 GPL15789 (Agilent) mid-luteal 42 68 GSE106602 GPL16791 (RNA-seq) mid-luteal 19 16 GSE71331 GPL19072 (Agilent) mid-luteal 7 5 Private Data Mars-seq2.0(RNA-seq) mid-luteal 12 21 ALL 104 134 Datasets used in identifying RIF subtypes Initial principal component analysis (PCA) revealed significant separation between the RIF and normal groups in each dataset (Supplementary Fig.  1 A and B), underscoring the substantial contribution of the endometrial factors in the RIF etiology. Subsequent differential gene expression analysis identified 5,792 ( GSE106602 ), 6,912 ( GSE111974 ), 352 (GSE58144_b1), and 1,584 (GSE58144_b2) DEGs between the normal and RIF groups across the batches ( p  < 0.05, limma). Assessment of the consistency of these DEGs across the training datasets revealed that only 2.23% of the identified genes were differentially expressed in at least three datasets (Fig.  2 A). Fig. 2 Identification and functional characterization of recurrent implantation failure (RIF)-related genes. (A) Venn diagram shows the overlap of gene expression data among four different datasets: GSE58144_B1 (red), GSE58144_B2 (light purple), GSE111974 (yellow), and GSE106602 (green). (B) Upset and pie charts illustrating the intersection of five previously published RIF gene signatures: Dec_score_2020 (purple), Shi_2018 (light green), Koot_2016 (blue), Bastu_2019 (pink), and Pathare_2017 (orange). (C) Heatmap showing expression patterns of RIF-related genes across training samples. (D) KEGG pathway enrichment analysis displaying significantly upregulated (red) and downregulated (blue) pathways associated with RIF-related genes in the training dataset Identification and functional characterization of recurrent implantation failure (RIF)-related genes. (A) Venn diagram shows the overlap of gene expression data among four different datasets: GSE58144_B1 (red), GSE58144_B2 (light purple), GSE111974 (yellow), and GSE106602 (green). (B) Upset and pie charts illustrating the intersection of five previously published RIF gene signatures: Dec_score_2020 (purple), Shi_2018 (light green), Koot_2016 (blue), Bastu_2019 (pink), and Pathare_2017 (orange). (C) Heatmap showing expression patterns of RIF-related genes across training samples. (D) KEGG pathway enrichment analysis displaying significantly upregulated (red) and downregulated (blue) pathways associated with RIF-related genes in the training dataset Furthermore, overlap evaluation among the five previously published sets of RIF-related genes [ 18 , 20 , 21 , 41 , 42 ] demonstrated a limited concordance, with a duplication rate of 15% (Fig.  2 B, Supplementary Table 4 ). These observations indicate that earlier analyses may have been significantly impacted by dataset-specific biases attributable to batch effects and inconsistent application of the inclusion and exclusion criteria. Consequently, reliable identification of robust molecular pathological changes associated with RIF requires the standardization of data acquisition protocols and the integration of multi-platform datasets. A total of 1,776 RIF-associated DEGs (871 upregulated and 905 downregulated) were identified in the integrated dataset following batch effect correction (Fig.  2 C). Notable genes among the significantly upregulated genes were: ADAMTS8 ( p  = 1 × 10 − 20 , FC = 4.79), a gene that inhibits embryo implantation through anti-angiogenic mechanisms [ 43 ]; PER1 ( p  = 1 × 10 − 20 , FC = 4.64), a core clock gene that plays an essential role in generating circadian rhythms and senses hormone dynamics necessary for homeostatic regulation of physiological functions [ 44 ]; and LPIN1 ( p  = 6.74 × 10 − 6 , FC = 4.27), a gene potentially implicated in infertility cases through its influence on lipid metabolism [ 45 ]. Notable genes among the significantly downregulated genes were: SPTLC1 ( p  = 1 × 10 − 20 , FC =-5.04), a gene essential for the process of decidualization [ 46 ]; EML4 ( p  = 1 × 10 − 20 , FC = -5.27), a key gene that is involved in the regulation of cellular processes, including microtubule dynamics; CA12 ( p  = 1 × 10 − 20 , FC = -4.77), whose expression is regulated by hypoxia and estrogen receptors and is significantly responsible for vascular protection. All these genes are implicated in the RIF etiology. Functional pathways enrichment analysis revealed distinct patterns in RIF (Fig.  2 D). Notably, the upregulated pathways in RIF included mineral absorption ( p  = 7.87 × 10 − 6 ) and nucleotide metabolism ( p  = 9.03 × 10 − 4 ), while the downregulated pathways were primarily associated with cell growth, signal transduction, and immune response. Notably, endometrial decidualization is significantly associated with cell senescence [ 47 ] ( p  = 1.56 × 10 − 2 ). These findings indicate that the RIF-associated DEGs potentially play significant regulatory roles during the process of embryo implantation. To identify genes that are significantly associated with endometrial dysfunction, we employed a forward-backward feature selection strategy combined with the KNN algorithm. We ultimately identified 54 key RIF- DEGs (Supplementary Fig.  2 A, Supplementary Table 5 ). Subsequently, given the substantial effects of hormones on endometrial receptivity, we analyzed the relationships between these RIF-associated DEGs and hormone receptor genes to elucidate the underlying regulatory patterns (Fig.  3 A). Notably, our results indicated a strong, significant positive correlation between estrogen receptor 1 ( ESR1 ) and CPE ( p  = 4.7 × 10⁻¹², R  = 0.47). Additionally, UNC119 , GRHL2 , and CDC14B showed significant correlations with progesterone receptors, PGR ( p  < 2 × 10⁻¹⁶, R  = 0.5), PGRMC1 ( p  = 7.8 × 10⁻¹⁰, R  = 0.43), and PGRMC2 ( p  < 2 × 10⁻¹⁶, R  = 0.72), respectively. Among prostaglandin receptors, MCM3AP-AS1 , PDSS2 , and SNX30 exhibited significant displayed correlations with PTGER1 ( p  = 2.4 × 10⁻⁹, R  = 0.41), PTGER2 ( p  = 7.9 × 10⁻⁸, R  = 0.38), and PTGIR ( p  = 2.4 × 10⁻¹², R = -0.47). Additionally, hydroxysteroid 17-β dehydrogenase 2 ( HSD17B2 ) exhibited a significant inverse correlation with PEPD ( p  = 1.8 × 10⁻¹³, R = -0.5). The androgen receptor ( AR ) exhibited a strong correlation with MEX3B ( p  = 8.7 × 10⁻¹⁵, R  = 0.52) (Supplementary Fig.  2 B). These results highlight a complex interplay between RIF-associated DEGs and hormone receptor genes. Fig. 3 Clinical and hormonal heterogeneity among RIF patients. (A) Heatmap illustrating the correlation between important RIF-DEGs and hormone receptors. *: p  < 0.05. (B) Heatmap showing the correlation between clinical characteristics and important RIF-DEGs. * p  < 0.1 and ** p  < 0.01. For categorical variables (batch and smoke), ANOVA was used for correlation analysis, while Spearman correlation was applied for other continuous variables. (C) Scatterplots showing representative correlations between individual clinical variables and the expression levels of important RIF-DEGs Clinical and hormonal heterogeneity among RIF patients. (A) Heatmap illustrating the correlation between important RIF-DEGs and hormone receptors. *: p  < 0.05. (B) Heatmap showing the correlation between clinical characteristics and important RIF-DEGs. * p  < 0.1 and ** p  < 0.01. For categorical variables (batch and smoke), ANOVA was used for correlation analysis, while Spearman correlation was applied for other continuous variables. (C) Scatterplots showing representative correlations between individual clinical variables and the expression levels of important RIF-DEGs The correlations between the identified RIF-associated DEGs and various clinical characteristics were investigated, with the result revealing several significant associations (Fig.  3 B). The top five genes that were strongly correlated with the specific clinical characteristics are presented in Fig.  3 C. Notably, the upregulated gene ADAMTS8 showed a significant negative correlation with the duration of the luteal phase (LH + days; p  = 8.1 × 10 − 11 , R = -0.49), a finding that was consistent with previous research [ 48 ]. We hypothesized that this gene potentially influences embryo implantation through its anti-angiogenic properties. The down-regulated gene ADAM17 displayed a significant negative correlation with BMI ( p  = 9.4 × 10  − 3 , R = -0.31), suggesting that its dysregulation may be associated with the dysfunction linked to obesity, which is an established risk factor for infertility. Furthermore, the upregulated gene NTRK3 was negatively significantly associated with smoking status ( p  = 3.11 × 10  − 3 , R = -0.28), another result consistent with observations in recent research [ 49 ], highlighting its potential role in smoking-related reproductive dysfunction. Additionally, the HUS1 gene demonstrated a significant positive correlation with age ( p  = 4.9 × 10  − 2 , R  = 0.17), with its effects potentially associated with embryo implantation failure through its role in inducing endometrial senescence. The study also revealed that DOK5 exhibited a borderline significant correlation with endometrial thickness ( p  = 0.07, R  = 0.27), which potentially contributes to endometrial abnormalities by inhibiting cell proliferation and differentiation. Collectively, these findings reveal the multifactorial etiological complexity of RIF at the genetic level. Given the heterogeneity of RIF patients, coupled with the clinical evidence that various endometrial diseases can cause infertility, we speculated that endometrial dysfunction in RIF involves multiple molecular pathological alterations. Subsequently, we conducted ConsensusClusterPlus analysis on the training dataset of 85 RIF patients using 1,776 RIF-associated DEGs to identify the molecular subtypes. Using the analysis from the CDF curve and Delta area, we identified k = 2 as the number of unique and non-overlapping subtypes (Fig.  4 A, Supplementary Fig.  3 A and B), and designated them as RIF-I and RIF-M. The PCA results showed significant differences between RIF-I, RIF-M, and the normal group, independent of batch effects (Supplementary Fig.  3 C). Fig. 4 Molecular classification of RIF subtypes based on differentially expressed genes. (A) Heatmap showing Consensus clustering heatmap of RIF samples using the “ConsensusClusterPlus” algorithm, identifying two robust molecular subtypes (K = 2). (B) Rank order plots displaying differentially expressed genes in RIF-I (left) and RIF-M (right) compared to controls. The top ten upregulated and downregulated genes in each subtype are labeled. Y-axis represents log₂ fold-change relative to the normal group. (C) Scatterplots showing the top two clinical variables associated with gene dysregulation in the RIF-I subtype: LH + days (left) and endometrial thickness (right). (D) Scatterplots showing the top two clinical variables associated with gene dysregulation in the RIF-M subtype: LH + days (left) and smoking status (right). (E) Gene set enrichment analysis (GSEA) for RIF-I (left) and RIF-M (right). Butterfly plots depict the top 10 significantly enriched genesets among upregulated (red bars) and downregulated (blue bars) pathways for each subtype Molecular classification of RIF subtypes based on differentially expressed genes. (A) Heatmap showing Consensus clustering heatmap of RIF samples using the “ConsensusClusterPlus” algorithm, identifying two robust molecular subtypes (K = 2). (B) Rank order plots displaying differentially expressed genes in RIF-I (left) and RIF-M (right) compared to controls. The top ten upregulated and downregulated genes in each subtype are labeled. Y-axis represents log₂ fold-change relative to the normal group. (C) Scatterplots showing the top two clinical variables associated with gene dysregulation in the RIF-I subtype: LH + days (left) and endometrial thickness (right). (D) Scatterplots showing the top two clinical variables associated with gene dysregulation in the RIF-M subtype: LH + days (left) and smoking status (right). (E) Gene set enrichment analysis (GSEA) for RIF-I (left) and RIF-M (right). Butterfly plots depict the top 10 significantly enriched genesets among upregulated (red bars) and downregulated (blue bars) pathways for each subtype We explored the differential expression of the 1,776 RIF-associated DEGs within the identified RIF-I and RIF-M subtypes. Notably, within the RIF-I subtype, 366 RIF-associated DEGs exhibited significant dysregulation, whereas 1,571 dysregulated RIF-associated DEGs were identified in the RIF-M subtype (Supplementary Table 6 ). Notably, the specific RIF-associated DEGs that were dysregulated significantly differed between the two subtypes. In the RIF-I subtype, key upregulated genes included the following: DNTTIP2 ( p  = 6.42 × 10 − 7 , FC = 4.89), responsible for encoding a nuclear protein that enhances estrogen receptor transcription; TCF7L2 ( p  = 1.51 × 10 − 6 , FC = 4.19), a key transcription factor in the WNT signaling pathway, plays a critical in endometrial progression [ 50 ]; RAB12 ( p  = 8.76 × 10 − 5 , FC = 3.75), associated with macrophage activation and Th1 response [ 51 ]; and STK3 ( p  = 1.04 × 10 − 5 , FC = 3.70), a key gene in the Hippo signaling pathway; NDUFS4 ( p  = 1.32 × 10 − 5 , FC = 3.64), which plays a significant role in mitochondrial function and oxidative phosphorylation and is upregulated in polycystic ovary syndrome and congenital uterine diseases [ 52 ] (Fig.  4 B). In the RIF-M subtype, prominently upregulated genes included the following: NFIA ( p  = 1.00 × 10 − 20 , FC = 7.93), a key transcriptional regulator associated with endometrial receptivity [ 53 ]; TM7SF2 ( p  = 1.00 × 10 − 20 , FC = 7.46), which play a significant role in lipid reprogramming and cholesterol synthesis was dysregulated in thin endometrium [ 54 ]; PEA15 ( p  = 1.00 × 10 − 20 , FC = 7.33), a negative regulator of apoptosis; and potassium ion channel-related genes TMEM175 ( p  = 1.00 × 10 − 20 , FC = 7.22) and LRRC26 [ 55 ] ( p  = 1.00 × 10 − 20 , FC = 7.21) (Fig.  4 B). Furthermore, across the identified subtypes, these dysregulated genes were potentially correlated with various clinical characteristics. In the RIF-I subtype, the dysregulated genes were primarily associated with LH + days and endometrial thickness: TCF7L2 exhibited the strongest positive correlation with LH + days ( p  = 0.025, R  = 0.27), while DNTTIP2 showed a significant negative correlation with endometrial thickness ( p  = 0.027, R = -0.33) (Fig.  4 C). In the RIF-M subtype, the dysregulated genes were primarily linked associated with LH + days and smoking status: TM7SF2 exhibited the strongest negative correlation with LH + days ( p  = 2.65 × 10 − 5 , R = -0.35), and NFIA demonstrated a highly significant negative correlation with smoking status ( p  = 1.69 × 10 − 6 , R = -0.18) (Fig.  4 D). To further elucidate the transcriptomic characteristics of the identified RIF subtypes, Gene Set Enrichment Analysis (GSEA) was conducted to compare the differences in the pathways involving the RIF-I, RIF-M, and the normal group. The results revealed that RIF-I patients exhibited significant activation of immune-related pathways, including the IL-17 ( p  = 2.73 × 10 − 3 ) and TNF ( p  = 5.86 × 10 − 3 ) pathways, as well as signal transduction pathways like Ras ( p  = 2.49 × 10 − 3 ) and MAPK ( p  = 3.05 × 10 − 3 ) signaling pathway (Fig.  4 E). Additionally, RIF-I patients exhibited inhibition of pathways related to fatty acid ( p  = 1.46 × 10 − 5 ) and carbon metabolism ( p  = 1.61 × 10 − 5 ). Conversely, RIF-M patients exhibited significant activation of the metabolism-related pathways, including oxidative phosphorylation ( p  = 9.69 × 10 − 8 ) and fatty acid metabolism ( p  = 4.07 × 10 − 4 ), along with hormone-related pathways such as insulin secretion ( p  = 4.86 × 10 − 3 ). Additionally, in the RIF-M subtype, pathways associated with cell growth, apoptosis, and immune-related pathways were significantly suppressed (Fig.  4 E). To characterize the immunological features of the two RIF subtypes, we performed CIBERSORT analysis of bulk transcript data based on single cell RNA sequence data of endometrium. RIF-I patients exhibited a higher proportion of lymphocyte infiltration compared to RIF-M patients (Supplementary Fig.  4 A). Furthermore, we utilized the characteristic expression profiles of 28 distinct immune cell types [ 56 ] to quantify their infiltration scores in each sample. Analysis revealed that multiple immune cell types including plasmacytoid dendritic cells (pDCs), macrophages, TH17 cells, and CD56 dim NK cells—were significantly enriched in RIF-I compared to RIF-M (Fig.  5 A and B). For example, pDCs exhibited the most substantial increase in infiltration in RIF-I, with an FC of 218.63 compared to 7.30 in RIF-M. Macrophages were upregulated by 15.44 FC in RIF-I while downregulated by -14.87 FC in RIF-M. TH17 and CD56 dim NK cells also showed marked infiltration in RIF-I but were depleted in RIF-M. Fig. 5 Distinct immune functions and metabolic features of RIF subtypes compared to controls. (A) Network illustrating immune cell interactions in the RIF-I subtype. (B) Immune cell interaction network for the RIF-M subtype. (C) Violin plots showing different immune functionalities, including chemokine activity, APC co-stimulation, inflammation, and cytotoxicity, across normal, RIF-I, and RIF-M groups. (D) Radar chart depicting metabolic profiles across the three groups. (E) Boxplots showing different levels of steroid and estrone metabolites in RIF-I and RIF-M patients compared to controls, with statistical significance indicated. (F) Representative immunohistochemical staining of T-bet and GATA3 across control, RIF-I, and RIF-M endometrial tissues Distinct immune functions and metabolic features of RIF subtypes compared to controls. (A) Network illustrating immune cell interactions in the RIF-I subtype. (B) Immune cell interaction network for the RIF-M subtype. (C) Violin plots showing different immune functionalities, including chemokine activity, APC co-stimulation, inflammation, and cytotoxicity, across normal, RIF-I, and RIF-M groups. (D) Radar chart depicting metabolic profiles across the three groups. (E) Boxplots showing different levels of steroid and estrone metabolites in RIF-I and RIF-M patients compared to controls, with statistical significance indicated. (F) Representative immunohistochemical staining of T-bet and GATA3 across control, RIF-I, and RIF-M endometrial tissues Subsequent correlation analyses revealed that RIF-I patients exhibited elevated infiltration of Th1 (FC = 3.01), Th17 (FC = 6.79), TFH (FC = 2.40) cells, and activated CD8 + T cells (FC = 0.68), suggesting coordinated immune activation. In contrast, RIF-M patients exhibited weaker immune networks, with downregulation of TFH cells (FC = -2.69), macrophages (FC = -14.87), and myeloid-derived suppressor cells (MDSCs) (FC = -4.03). Functional pathway analysis supported this distinction: RIF-I was associated with higher expression of genes involved in chemotaxis, cytotoxicity, (Enrichment score = 0.22, p  < 0.001), and increased cytotoxic activities (Enrichment score = 0.33, p  < 0.001), whereas RIF-M was characterized by reduced antigen presentation capacity and a suppressed inflammatory profile (Enrichment score = -0.33, p  < 0.001) (Fig.  5 C). Metabolic pathway analysis further distinguished the subtypes. RIF-I samples exhibited downregulation of the tricarboxylic acid (TCA) cycle (Enrichment score = -0.24, p  < 0.001), with compensatory upregulation of arginine biosynthesis (Enrichment score = 0.13, p  < 0.05) (Fig.  5 D and E). In contrast, RIF-M samples demonstrated increased activity in multiple metabolic domains, including the TCA cycle (Enrichment score: 0.17, p  < 0.05), oxidative phosphorylation (Enrichment score = 0.19, p  < 0.01). lipid metabolism (Enrichment score = 0.07), particularly steroid (Enrichment Score = 0.20, p  < 0.01) and fatty acid (Enrichment score = 0.28, p  < 0.01). RIF-M also showed elevated purine metabolizing activity (Enrichment score = 0.09, p  < 0.01). MetaFlux analysis revealed abnormally high fluxes through steroidogenesis and estrogen metabolism in RIF-M patients (Supplementary Fig.  4 B), suggesting altered hormonal processing. To validate the molecular subtypes at the protein level, we performed immunohistochemistry for T-bet (TBX21) and GATA3 in mid-luteal endometrial samples. In RIF-I samples, T-bet showed markedly higher expression compared to normal and RIF-M, while GATA3 was predominantly expressed in RIF-M. Notably, the T-bet/GATA3 ratio was elevated in RIF-I and reduced in RIF-M, consistent with transcriptomic findings (Fig.  5 F). Notably, the clinical application of the identified large number of RIF-DEGs is a major challenge. Therefore, we utilized the forward-backward feature selection method to identify a refined set of signature genes, selecting 64 RIF-I DEGs (Supplementary Table 7 ) and 19 RIF-M DEGs (Supplementary Table 8 ). Next, we established several robust machine-learning models based on 64 combinations of algorithms within the training dataset. The predictive power of each model was evaluated by calculating the F-score across all cohorts (Supplementary Fig.  5 ). MetaRIF, which represents an ensemble classifier combining the Nnet (neural network) and RF (random forest) algorithms, had the highest average F-score among the tested models. Furthermore, principal component analysis provided evidence that the selected genes could accurately predict the RIF status and help to distinguish among the identified RIF subtypes (Supplementary Fig.  6 A). The MetaRIF algorithm utilizes the hierarchical series of decision rules to classify the patient’s RIF status based on gene expression profiles (Fig.  6 A). Initially, it determines whether the patient’s normal scores exceed RIF-I and RIF-M’s. If this condition is satisfied, the endometrial status is classified as “Normal”. If not, the MetaRIF evaluates the RIF-I score, if it exceeds normal levels and exceeds the RIF-M score by at least 0.1, the classification is designated as “RIF-I”. In cases where neither Normal nor RIF-I criteria are met, the RIF-M score is calculated and analyzed. If the RIF-M score exceeds normal levels and is greater than the RIF-I score by 0.1, the status is classified as “RIF-M”. If none of these classifications apply, the patient is categorized as “Uncertain”. This structured methodology enables precise RIF status classification using gene expression profiles, providing a more robust foundation for clinical decision-making. Fig. 6 Development of the MetaRIF classifier, clinical characterization of RIF subtypes, and identification of candidate therapeutics. (A) Schematic workflow of the MetaRIF classification system. The process includes gene filtering to produce RIF-I and RIF-M arrays, which are analyzed to classify samples into subtypes of Normal, RIF-I, RIF-M, or Uncertain based on specific thresholds. (B) Bar plot showing the top 10 genes contributing most significantly to the prediction of RIF-I and RIF-M subtypes within the MetaRIF model. (C) Receiver operating characteristic (ROC) curves demonstrating the performance of the MetaRIF classifier in two independent test datasets: GSE71331 and an in-house cohort. (D) Comparative ROC analysis showing that MetaRIF outperforms two published gene signatures (Koot_sig, Wang_sig and OSR_score) in subtype classification accuracy. (E) Boxplots comparing clinical parameters—body mass index (BMI) and antinuclear antibody (ANA) levels—among Normal, RIF-I, and RIF-M groups. (F) Scatterplot of Connectivity Map (CMAP) scores identifying candidate compounds associated with RIF-I and RIF-M subtypes, highlighting distinct drug response profiles Development of the MetaRIF classifier, clinical characterization of RIF subtypes, and identification of candidate therapeutics. (A) Schematic workflow of the MetaRIF classification system. The process includes gene filtering to produce RIF-I and RIF-M arrays, which are analyzed to classify samples into subtypes of Normal, RIF-I, RIF-M, or Uncertain based on specific thresholds. (B) Bar plot showing the top 10 genes contributing most significantly to the prediction of RIF-I and RIF-M subtypes within the MetaRIF model. (C) Receiver operating characteristic (ROC) curves demonstrating the performance of the MetaRIF classifier in two independent test datasets: GSE71331 and an in-house cohort. (D) Comparative ROC analysis showing that MetaRIF outperforms two published gene signatures (Koot_sig, Wang_sig and OSR_score) in subtype classification accuracy. (E) Boxplots comparing clinical parameters—body mass index (BMI) and antinuclear antibody (ANA) levels—among Normal, RIF-I, and RIF-M groups. (F) Scatterplot of Connectivity Map (CMAP) scores identifying candidate compounds associated with RIF-I and RIF-M subtypes, highlighting distinct drug response profiles Finally, the main predictive genes for RIF-I decision-making were LIPG , PA2G4 , MUTYH , and TRAPPC1 . In contrast, the main predictive genes for RIF-M were NFIA , RNF10 , and GRHL2 (Fig.  6 B). In addition, we tested the classification performance of the MetaRIF algorithm in two independent validation datasets. In GSE71331 , the predicted area under the receiver operating characteristic curve (AUC) for RIF was 0.94, while in the in-house dataset, it was 0.85 (Fig.  6 C). Subsequently, we compared the MetaRIF’s performance with two currently published gene signatures for predicting RIF in the integrated validation dataset. The results indicated that the MetaRIF significantly outperformed the koot_sig [ 21 ], Wang_sig [ 25 ] and Oxidative Stress-Related genes (OSR_score) [ 57 ] in terms of classification accuracy (AUC: MetaRIF = 0.88, Koot_sig = 0.48, Wang’s_sig = 0.54, OSR_score = 0.72) (Fig.  6 D). To explore clinical distinctions between the molecular subtypes, we analyzed key patient characteristics including body mass index (BMI), antinuclear antibody (ANA) levels, and endometrial aging markers. RIF-M patients exhibited significantly higher BMI compared to RIF-I patients, consistent with a metabolism-associated phenotype. Conversely, RIF-I patients, characterized by heightened immune activity, showed significantly elevated ANA levels, reflecting systemic immune activation (Fig.  6 E). We also observed a trend toward accelerated endometrial aging in RIF-M patients relative to controls ( p  = 2.7 × 10⁻⁸; Supplementary Fig.  6 B), further supporting the hypothesis that metabolic stress contributes to impaired endometrial function in this subtype. To identify potential therapeutic compounds targeting each RIF subtype, we performed drug-response signature matching using the CMap database. Several compounds were positively associated with the RIF-I transcriptomic profile, including sirolimus, CAY-10,618, alimemazine, tyrphostin AG-1478, and phloretin. For RIF-M, top-ranking candidate compounds included prostaglandins, PJ-34, apoptosis activator II, amonafide, and helveticoside (Fig.  6 F). Notably, sirolimus, the leading candidate drug for RIF-I, is an immunosuppressive agent that inhibits interleukin-2–mediated T cell proliferation. Its therapeutic potential in RIF was supported by a previous double-blind study showing benefit in patients with elevated Th17/Treg ratios [ 58 ]. For RIF-M, prostaglandins emerged as the most promising candidate. Prostaglandins are essential for endometrial remodeling and vascular permeability, and their deficiency has been linked to reduced implantation success and infertility [ 59 ].

Microarray expression datasets for RIF were retrieved from the Gene Expression Omnibus (GEO) database ( GSE111974 , GSE71331 , GSE58144 , and GSE106602 ). The GSE111974 dataset, based on the GPL17077 microarray platform, includes gene expression profiles of 24 RIF cases and 24 healthy controls in the natural menstrual cycle. The GSE71331 dataset, based on the GPL9072 microarray platform, included gene expression profiles of 7 RIF cases and 5 healthy controls in the natural menstrual cycle. The GSE58144 dataset, based on the GPL15789 microarray platform, includes gene expression profiles of 42 RIF cases and 68 healthy controls in the natural menstrual cycle. Lastly, the GSE106602 dataset, derived from the GPL16791 microarray platform, features gene profiles of 19 RIF cases and 16 healthy controls in the natural menstrual cycle. All patients who failed due to endometrial cavity abnormalities or endometrium disorders were excluded (Supplementary Table 1 ). Additionally, endometrial biopsy samples were collected from 33 women at Shenzhen Zhongshan Obstetrics & Gynecology Hospital. The participants were categorized into two groups: women diagnosed with RIF ( n  = 12) who did not achieve clinical pregnancy following three or more embryo transfers, encompassing a cumulative total of ten or more embryos; women with tubal factor infertility ( n  = 21) who successfully achieved a clinical pregnancy after the initial single embryo transfer and served as the Normal group. Notably, all participants were recruited for this study, with all the samples obtained—derived from human endometrial tissue—meeting the inclusion criteria described below. The following inclusion criteria was employed on all the eligible individuals: age between 18 and 38 years; a BMI within the range of 18 to 25 kg/m²; no documented case of pregnancy within past 18 months; a history of regular menstrual cycles within 25 to 35 days; evidence of normal fallopian tube patency; evidence of normal ovarian function; and absence of any hormonal treatments over the past three months before endometrial biopsy. The exclusion criteria used were as follows: a history of intrauterine pathologies, including congenital uterine anomalies, uterine fibroids, endometrial polyps, and intrauterine adhesions; presence of hydrosalpinx; a diagnosis of polycystic ovary syndrome (PCOS); confirmed diagnosis of endometriosis or adenomyosis; a history of abnormal chromosomal karyotypes in either partner or previous abortion tissues; serological evidence of lupus anticoagulant or anticardiolipin antibodies; presence of active infectious diseases; pre-existing endocrine disorders, defined as abnormal blood glucose levels or thyroid dysfunction; and current use of hormonal contraception, specifically contraceptive pills or intrauterine devices over the past six months. Additionally, all RIF samples included in our study were evaluated for the presence of CD138 + plasma cells and chronic endometritis, with such cases being excluded from the analysis. The assessment of the criteria was conducted using a multifaceted approach encompassing vaginal ultrasonography, hysteroscopy, laparoscopy, and karyotyping analysis, supplemented by relevant hormonal and immunological assays. Endometrial tissue specimens were collected during the mid-secretory phase of the menstrual cycle, specifically between 5 and 8 days after the peak of luteinizing hormone. The precise timing of endometrial sampling was further corroborated using histological evaluation based on Noyes’ criteria [ 26 ]. All tissue samples were immediately cryopreserved and stored at − 80 °C for subsequent analysis. The tissue samples were immersed and rinsed with plain RPMI − 1640 (Gibico) to remove blood and mucus. Total RNA was isolated using Qiagen RNeasy Mini Kits (Qiagen) following the manufacturer’s protocol, followed by the preparation of the transcriptome libraries using the Massively Parallel Single-Cell RNA-seq method (MARS-seq). Notably, all mRNA (1–50 ng) was barcoded, reverse transcribed into cDNA, and pooled together. Subsequently, the pooled single-strand cDNA was converted into double-strand DNA, followed by linear amplification into RNA using in vitro transcription. The resulting RNA was fragmented and ligated using Illumina sequences for sequencing and barcoding (DNA). Finally, the RNA-DNA chimera was reverse transcribed into DNA and amplified using the Polymerase Chain Reaction (PCR) sequencing technique. The quality of the constructed library was evaluated using real-time PCR, while its concentration level was quantified using Qubit 2.0. Endometrial tissues were fixed in formalin and embedded in paraffin. Tissues were sectioned on a rotary microtome at 4 μm thickness and mounted on glass slides. Slides from the embedded tissues were baked for 60 min at 60℃, deparaffinized in xylene, and rehydrated in graded concentrations of ethanol in water. T-bet (BD Pharmingen, 561262) and GATA3 (R&D systems, 634913) immunostaining was performed on a Leica Biosystems BOND ® RX instrument. Counterstained with Hematoxylin. Heat mediated antigen retrieval with Tris-EDTA buffer (pH 9.0, epitope retrieval solution 2) for 10 min. Slides were scanned with Olympus VS200 and images were analyzed with the HALO Image Analysis System (Indica Labs, Albuquerque NM, USA). Processed datasets corresponding to data generated on the Agilent platform were directly downloaded. For the datasets generated using the Illumina platform, sequencing reads were mapped to the human reference genome hg38 using HISAT [ 27 ] (version 0.1.6). Notably, reads with multiple mapping positions were excluded. Reads were associated with genes if they were mapped to an exon using Homo sapiens Ensembl90 for reference. Additionally, reads counts were raw count values, and the fragments per kilobase of transcript per million mapped reads normalized count values were extracted. The expression data used in our study were quality-controlled and quantile-normalized using the Affy package [ 28 ]. The random effect method was then employed to integrate the datasets at the mRNA level [ 29 ], this approach offers flexibility and robustness comparable to the classical procedure while effectively reducing batch effect (Supplementary Fig.  1 B). Subsequently, upregulated and downregulated expressed genes were identified using the MetaDE [ 30 ] (version 2.2.3) package in R ( p  < 0.05). The CIBERSORT [ 31 ], a gene expression-based deconvolution algorithm, was initially used to quantify cell fractions from gene expression profiles in bulk tissues. Subsequently, single-cell RNA sequencing data [ 32 , 33 ] was leveraged to construct an endometrial special reference signature gene file, with 1000 permutations integrated to quantify the relative proportions of the 14 endometrial cell types (Supplementary Table 2 ). Additionally, we collected 28 genes characteristically expressed by immune cells to refine the infiltration degree of immune cell types [ 34 ]. We then used the Gene Set Variation Analysis (GSVA) [ 35 ] method to determine the infiltration scores of these immune cells in each sample separately. Based on the RIF-DEGs, consensus clustering was performed in the training cohort using the ConsensusClusterPlus package [ 36 ]. We set the parameters as the Euclidean distance-based Pam algorithm with 1000 iterations, with 80% of the samples drawn randomly at each iteration. The cumulative distribution function (CDF) was used to select an optimal number of clusters, which typically ranged from 2 to 8 clusters. Notably, two clusters were identified: RIF-I and RIF-M. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the clusterProfiler [ 37 ] package in R to further examine the biological function of RIF-related genes. Reference gene sets were sourced from KEGG and public research databases [ 38 ]. The GSVA was employed to compare the functional module activity between the RIF-I and RIF-M groups. Metabolic fluxes were determined and normalized from gene expression data using the METAFlux [ 39 ]. A p  < 0.05 was considered statistically significant. We evaluated our dataset using eight machine-learning algorithms encompassing 64 model configurations to design a consensus signature optimized for accuracy and reliability. Due to the wide range of applications and the established nature of these methods, we selected eight algorithms: Random Forest (RF), Lasso, Naive Bayes, K-Nearest Neighbors (KNN), Neural Network (Nnet), Light Gradient-Boosting Machine (LightGBM), AdaBoost.M1, and Support Vector Machine (SVM). Initially, genes associated with RIF-I and RIF-M were incorporated into 64 algorithmic combinations to train predictive models within the training cohort, utilizing leave-one-out cross-validation (LOOCV) for model evaluation. Subsequently, all predictive models were validated using the GSE71331 dataset and an independent in-house dataset. Ultimately, the F-score was computed for all training and validation datasets for each model, and the model with the greatest average F-score was selected as ideal. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:F-score=\frac{2*sensitivity*specificity}{sensitivity+specificity}$$\end{document} The Connectivity Map (CMap) [ 40 ], which integrates large perturbational datasets, advances our understanding of pathological conditions, thereby accelerating the development of novel therapeutic interventions. The query tool of a cloud-based software platform, CMap and LINCS Unified Environment (CLUE), was used to conduct drug prediction, mode of action analysis (MOA), and drug-target analysis. Briefly, a connectivity score was obtained by comparing the query gene set to the reference data from the database [ 40 ]. Notably, a positive score indicates the similarity between a given perturbagen’s signature and the query’s gene, whereas a negative score indicates dissimilarity between the signatures. Overlapping DEGs were then subjected to the query tool with specific parameters configured as “Gene expression (L1000)”. Notably, RIF-I and RIF-M scores were determined as the inverse of the connectivity score, with these scores used to classify the samples. The top five small-molecule compounds with the highest positive RIF-I or RIF-M score were visualized using the ggplot2 package. All statistical analyses were performed using R software (version 4.1.1). The Wilcoxon test was used to compare data between the two groups, while the Spearman correlation analysis was employed to determine the correlation between the immune infiltration scores or genes and disease subtypes. The diagnostic efficacy of the model was calculated using receiver operator characteristic (ROC) curves and area under the curve (AUC) using the pROC (version 1.18.5) package in R. A p-value < 0.05 was considered statistically significant, unless stated otherwise.

Recurrent implantation failure (RIF) remains a persistent barrier in routine ART clinics, with patients and clinicians often left without definitive answers despite the transfer of morphologically normal embryos. Previous studies of endometrial factors in RIF have typically focused on “endometrial receptivity”, where the concept of WOI is often confused with receptivity itself. However, WOI only reflects a temporal window, whereas endometrial receptivity presents a broader range of molecular and structural properties for embryo acceptance. This conceptual confusion has led to misdiagnoses and ineffective, one-size-fits-all treatment strategies. In this study, we integrate multi-cohort transcriptomic data with clinical profiling to identify two distinct molecular subtypes of RIF—an immune-driven subtype (RIF-I) and a metabolic-driven subtype (RIF-M). These two subtypes are associated with unique molecular signatures, supporting the hypothesis that RIF is not a uniform clinical entity but rather results from diverse pathogenic mechanisms. The RIF-I subtype is characterized by significant activation of immune signaling pathways, including IL-17 and TNF, alongside elevated infiltration of effector immune cells such as Th1, Th17, TFH, and cytotoxic CD56 dim NK cells, and elevated expression of immune-related genes such as, such as LIPG [ 60 ], PA2G4 [ 61 ], and TRAPPC [ 62 ]. Notably, PA2G4 and DNTTIP2 both displayed a negative correlation with endometrial thickness, suggesting that RIF-I may be linked to insufficient decidual transformation or structural endometrial inadequacy [ 63 ]. Elevated Th17/Treg ratios and cytotoxic NK activity have long been associated with implantation failure and adverse pregnancy outcomes [ 6 , 64 , 65 ]. Interestingly, we did not observe marked changes in CD56 bright NK cells, which have previously been implicated in many RIF patients [ 66 , 67 ]. This discrepancy may suggest that the cytotoxic CD56 dim subset, rather than the more tolerogenic CD56 bright NK cells, may be more relevant in RIF-I pathogenesis. Macrophages were also significantly elevated in RIF-I samples, consistent with their involvement in other inflammatory gynecological disorders such as endometriosis and adenomyosis [ 68 , 69 ]. In these conditions, macrophages contribute to chronic inflammation, tissue remodeling, and disease progression through robust inflammatory responses and overexpression of macrophage-related markers, such as CD68 and CD163 [ 70 ]. Our previous study also showed that absent expression of CD68 is a risk factor for embryo implantation failure [ 71 ]. However, we were unable to characterize macrophage polarization (M1 vs. M2) or spatial localization in this study. Future integration of single-cell data will be essential to clarify these functional roles. Patients classified under the RIF-I group also exhibited elevated levels of antinuclear antibodies (ANA), supporting the presence of an autoimmune contribution. These findings align with earlier reports associating autoimmunity with implantation failure and provide additional rationale for considering immunomodulatory therapy in this subtype [ 72 – 74 ]. These patterns raise the possibility that RIF-I may be more prevalent among patients with abnormal endometrial development or insufficient decidualization, thus emphasizing the need for a comprehensive evaluation incorporating molecular, immunologic, and morphologic assessments of the endometrium. In contrast, RIF-M is characterized by a distinct molecular phenotype manifesting as enhanced oxidative phosphorylation, fatty acid metabolism, and steroid hormone biosynthesis. These findings agree with recent research emphasizing the significance of metabolic-immune interactions in reproductive health [ 60 ]. Predictive genes for this subtype, such as NFIA [ 53 ], GRHL2 [ 75 ], and AASDHPPT [ 76 ] are linked to endometrial tolerance and metabolic processes. Notably, NFIA and AASDHPPT exhibited a strong negative correlation with smoking status, suggesting that environmental factors such as tobacco exposure may exacerbate metabolic dysfunction in the endometrial tissues [ 77 ]. We also identified dysregulated expression of the circadian rhythm gene PER1, which is known to regulate hormone-responsive gene expression [ 44 , 78 ]. Dysregulated PER1 expression may lead to altered hormone-responsive gene expression and a shift in the window of implantation. This hypothesis is supported by clinical observations from endometrial receptivity assays (ERA), suggesting that some patients with RIF have a temporal asynchrony in endometrial development [ 79 , 80 ]. Mechanistically, circadian rhythm disturbances may not only reflect intrinsic metabolic dysregulation but may also be influenced by external factors, including medication use [ 44 ]. Given that some RIF patients may receive long-term glucocorticoids or use hormone replacement therapy to prepare the endometrium during repeated ART cycles, such iatrogenic effects on clock genes may lead to impaired receptivity in susceptible individuals [ 81 ]. However, this possibility warrants further investigation. Importantly, our findings also underscore the value of lifestyle modification as a complementary strategy for RIF-M patients [ 7 ], including obesity, type-2 diabetes, insulin resistance, and polycystic ovary syndrome (PCOS) [ 82 , 83 ]. Elevated BMI, disrupted energy metabolism, and evidence of premature endometrial aging suggest that weight management, dietary intervention, and metabolic support may play a meaningful role in improving endometrial receptivity [ 84 ]. Interventions targeting insulin resistance, oxidative stress, and circadian rhythm alignment—such as structured physical activity and micronutrient support—could be particularly relevant for this group. Furthermore, we extended our molecular findings to the protein level by validating key markers using immunohistochemistry. We selected T-bet ( TBX21 ) and GATA3, based on their transcriptomic specificity and established biological relevance. The T-bet/GATA3 ratio has been developed of Th1/Th2 immune balance and is clinically associated with excessive inflammation and poor reproductive outcomes [ 85 ]. Our IHC results confirmed elevated T-bet expression in RIF-I samples and higher GATA3 expression in RIF-M, consistent with the immune-metabolic dichotomy observed at the transcriptomic level. In addition, we explored differential genes shared between RIF-I and RIF-M in an attempt to discover shared molecular alterations between RIFs. Through our analysis, we identified 1 up-regulated and 13 down-regulated genes (Supplementary Table 9 ). The results showed limited shared genes between the different subtypes, mainly involving disruption of functions such as cell cycle regulation and immune microenvironment homeostasis. This suggests to us that these alterations in shared genes and pathways may be secondary consequences of implantation failure. In order to make it possible to translate these findings into clinical practice, we developed the MetaRIF classifier, successfully distinguishing these molecular subtypes in independent validation datasets. Our results suggest that MetaRIF outperforms existing genetic signatures and provides a practical tool for stratifying patients and guiding individualized treatment [ 21 , 57 ], although these require further optimization and prospective validation. For RIF-I, which is driven by immune activation and autoimmune features, immunosuppressive therapy (particularly with rapamycin, an mTOR inhibitor known to modulate T-cell responses) may help restore endometrial immune tolerance and improve implantation potential. In contrast, RIF-M is characterized by metabolic and circadian dysregulation, often associated with higher BMI, smoking, and features of the metabolic syndrome. This prediction is supported by multiple lines of evidence. For example, a case-control study in women with previous implantation failure demonstrated reduced endometrial expression of cyclooxygenase-2 and terminal prostaglandin synthases compared to fertile controls [ 59 ]. Additionally, experimental data in mice have shown that early embryonic signals trigger PGE₂ release via epithelial sodium channel (ENaC) activation, which in turn supports decidualization—a critical step for implantation—through EP2/EP4-mediated cAMP and CREB signaling [ 86 ]. These findings align with our computational flux analysis uncovering impaired steroid-to-prostaglandin conversion in RIF-M. Taken together, this evidence provides strong biological plausibility for targeting prostaglandin pathways—especially PGE₂ analogues or EP receptor agonists—as a therapeutic strategy for metabolism-dominant RIF. Further preclinical and prospective clinical studies are warranted to test this hypothesis. Despite these promising findings, several limitations should be acknowledged. First, our analysis was primarily computational and transcriptomic, necessitating further experimental validation through protein-level assessments, functional immunological assays, and targeted metabolic profiling. While bioinformatics approaches are powerful in identifying patterns and generating hypotheses, they are inherently dependent on the quality and scope of the available datasets, limiting the generalizability of the findings. Second, the current binary classification of RIF into RIF-I and RIF-M may overlook the full biological complexity of recurrent implantation failure. Each subtype may contain additional subgroups with different pathogenic drivers. For example, the immune-driven group may contain both infection-related immune activation and autoimmune-driven inflammation. Similarly, the metabolic subtype can be further divided into hormone resistance-dominated and circadian rhythm disturbance-dominated forms of dysfunction. Additionally, Finally, although we rigorously excluded chronic endometritis (i.e., CD138⁺ plasma cell infiltration) in our in-house cohort through histological assessment, such confirmation was not uniformly available across public datasets. This introduces potential confounding, particularly regarding the immune activation observed in RIF-I samples. While T-bet and GATA3 validation supports the immune signature of RIF subtypes, broader protein-level characterization and functional assays will be necessary to deepen mechanistic insight. Finally, our analyses focused on mid-luteal phase biopsies; dynamic changes across the menstrual cycle and embryo transfer timing were not addressed and warrant future investigation.

In conclusion, by identifying immune-driven and metabolic-driven subtypes, this study provides an initial approach to classifying recurrent implantation failure based on endometrial molecular features. Further refinement of the MetaRIF classifier may contribute to more personalized treatment decisions in the future, but prospective validation and clearer subtype definitions are still needed before clinical application.

Recurrent implantation failure (RIF), defined as failure to achieve clinical pregnancy after three or more good-quality embryo transfers, remains a major challenge in assisted reproductive technology (ART) [ 1 – 3 ]. RIF is not caused by a single factor but is the result of a combination of potential causes. This biological and clinical heterogeneity significantly complicates both diagnosis and treatment, contributing to significant physical, emotional, and financial burdens on patients and their families [ 4 , 5 ]. While embryo-related abnormalities such as aneuploidy have been widely recognized, increasing attention has shifted toward maternal factors—particularly endometrial dysfunction—as critical contributors to RIF [ 6 ]. These include immunological dysregulation, inflammatory responses, and impaired endometrial receptivity. Clinical conditions such as chronic endometritis (CE), intrauterine adhesions, and elevated body mass index (BMI) have also been identified as important risk factors [ 7 ]. However, the mechanisms by which these diverse factors impair implantation remain poorly understood. Successful embryo implantation depends on finely tuned communication between the embryo and the endometrium during the window of implantation (WOI). This process is regulated by a complex network of hormones, immune cells, and molecular signaling pathways [ 8 ]. Disruptions in these regulatory systems—such as abnormal infiltration of natural killer (NK) cells or macrophages, or cytokine imbalances (e.g., Th1/Th2, Th17/Treg)—can interfere with implantation. Although advances in imaging and diagnostics have improved the identification of gross anatomical abnormalities, many RIF cases remain unexplained due to the limited understanding of molecular alterations in the endometrial microenvironment [ 9 , 10 ]. This diagnostic uncertainty often results in empirical and inconsistent treatment strategies. Current clinical therapeutic approaches for managing RIF-related endometrial dysfunction mainly focus on hormonal regulation, improvement of endometrial receptivity, immunomodulation, and antithrombotic therapy [ 8 ]. These empirical treatments aim to improve the local microenvironment of the endometrium; however, their efficacy varies significantly among individuals, and clinical benefits are often inconsistent, with a lack of universally accepted standardized protocols [ 8 , 11 , 12 ]. For example, granulocyte colony-stimulating factor (G-CSF) has been proposed to enhance endometrial receptivity and improve reproductive outcomes [ 13 ]. However, a randomized trial found that in RIF patients with normal endometrial morphology, administration of G-CSF did not significantly improve clinical pregnancy or live birth rates compared to placebo [ 14 ]. Such inconsistencies underscore the importance of understanding the molecular basis of RIF and the need for more detailed bioinformatics studies to reduce diagnostic and therapeutic uncertainty caused by endometrial heterogeneity. The rapid advancement of high-throughput sequencing technologies has provided unprecedented opportunities to systematically investigate the complex molecular landscape of the endometrium in RIF patients [ 15 ]. Previous studies have used techniques such as microarrays and RNA sequencing to preliminarily reveal the potential links between abnormal signaling pathways and RIF [ 16 – 21 ]. However, the results of these studies often lack consistency due to limitations such as single-center designs, small sample sizes, heterogeneous technical platforms, and variations in sampling standards [ 22 , 23 ]. Therefore, integrating multicenter, multi-platform datasets to objectively reveal key molecular characteristics and subtypes of RIF-related endometrial factors is crucial in overcoming current research limitations and advancing the field. Recent transcriptomic efforts have identified potentially meaningful RIF subtypes for example, an “immune-activated” phenotype marked by elevated cytokine signatures, and a “mitochondrial-deficient” phenotype associated with reduced oxidative phosphorylation [ 24 , 25 ]. While these studies offer important mechanistic insights, they did not yield clinically deployable classifiers, nor did they validate subtype-specific features across multiple biological levels. Moreover, they focused on individual biological axes in isolation, without addressing the potential crosstalk between immune and metabolic dysfunction. In this study, we innovatively integrate multiple independent datasets to define and validate molecular subtypes of endometrial-related RIF. We identify two major subtypes of endometrial-related RIF: immune activation type (RIF-I) and metabolic disorder type (RIF-M), and establish a corresponding molecular classification system. Additionally, based on an in-depth analysis of the subtype-specific molecular characteristics, we preliminarily screened potential targeted intervention drugs, including rapamycin for RIF-I and prostaglandin for RIF-M. To our knowledge, this is the first study to propose a molecularly defined subtyping system for RIF based on endometrial features. By characterizing these subtypes and linking them to personalized treatment approaches, we aim to reduce diagnostic ambiguity and therapeutic inefficacy caused by RIF heterogeneity. A schematic of the study process is shown in Fig.  1 . Fig. 1 Flowchart for depicting the development, validation, and application of MetaRIF Flowchart for depicting the development, validation, and application of MetaRIF

Below is the link to the electronic supplementary material. Supplementary Material 1: The inclusion and exclusion criteria within RIF in public data (Table S1 ) Supplementary Material 1: The inclusion and exclusion criteria within RIF in public data (Table S1 ) Supplementary Material 2: The list of endometrium cell types (Table S2 ) Supplementary Material 2: The list of endometrium cell types (Table S2 ) Supplementary Material 3: Demographic characteristics of the women included in this study (Table S3 ) Supplementary Material 3: Demographic characteristics of the women included in this study (Table S3 ) Supplementary Material 4: The list of published RIF-related marker genes (Table S4 ) Supplementary Material 4: The list of published RIF-related marker genes (Table S4 ) Supplementary Material 5: The filtered RIF-related gene set (Table S5 ) Supplementary Material 5: The filtered RIF-related gene set (Table S5 ) Supplementary Material 6: Dysregulated RIF-related genes in distinct RIF subtypes (Table S6 ) Supplementary Material 6: Dysregulated RIF-related genes in distinct RIF subtypes (Table S6 ) Supplementary Material 7: The list of RIF-I related differential genes (Table S7 ) Supplementary Material 7: The list of RIF-I related differential genes (Table S7 ) Supplementary Material 8: The list of RIF-M related differential genes (Table S8 ) Supplementary Material 8: The list of RIF-M related differential genes (Table S8 ) Supplementary Material 9: The validated shared differential genes among both RIF classes (Table S9 ) Supplementary Material 9: The validated shared differential genes among both RIF classes (Table S9 ) Supplementary Material 10: The supplementary figures Supplementary Material 10: The supplementary figures