Identification of potential biomarkers and immune infiltration characteristics in recurrent implantation failure using bioinformatics analysis

Introduction: Recurrent implantation failure (RIF) is a frustrating challenge because the cause is unknown. The current study aims to identify differentially expressed genes (DEGs) in the endometrium on the basis of immune cell infiltration characteristics between RIF patients and healthy controls, as well as to investigate potential prognostic markers in RIF.

Methods: GSE103465, and GSE111974 datasets from the Gene Expression Omnibus database were obtained to screen DEGs between RIF and control groups. Gene Ontology analysis, Kyoto Encyclopedia of Genes and Genomes Pathway analysis, Gene Set Enrichment Analysis, and Protein-protein interactions analysis were performed to investigate potential biological functions and signaling pathways. CIBERSORT was used to describe the level of immune infiltration in RIF, and flow cytometry was used to confirm the top two most abundant immune cells detected.

Results: 122 downregulated and 66 upregulated DEGs were obtained between RIF and control groups. Six immune-related hub genes were discovered, which were involved in Wnt/-catenin signaling and Notch signaling as a result of our research. The ROC curves revealed that three of the six identified genes (AKT1, PSMB8, and PSMD10) had potential diagnostic values for RIF. Finally, we used cMap analysis to identify potential therapeutic or induced compounds for RIF, among which fulvestrant (estrogen receptor antagonist), bisindolylmaleimide-ix (CDK and PKC inhibitor), and JNK-9L (JNK inhibitor) were thought to influence the pathogenic process of RIF. Furthermore, our findings revealed the level of immune infiltration in RIF by highlighting three signaling pathways (Wnt/-catenin signaling, Notch signaling, and immune response) and three potential diagnostic DEGs (AKT1, PSMB8, and PSMD10).

Conclusion: Importantly, our findings may contribute to the scientific basis for several potential therapeutic agents to improve endometrial receptivity.

1 Introduction

Although there is no unified definition for recurrent implantation failure (RIF), it is usually referring to women who fail to achieve a clinical pregnancy after the transfer of at least four high-quality embryos in a minimum of three fresh or frozen cycles in a woman under the age of 40 years (1, 2). As multiple failed cycles can often cause frustration and desperation for the struggling couples, the causes for these failures remain unknown. According to a recently published retrospective cohort study, the incidence of RIF was around 5% in women who had an anatomically normal uterus and underwent consecutive single euploid blastocyst transfer (3), among which 50% of affected cases remain unexplained. It is essential and urgent to determine the etiology of idiopathic RIF and propose novel therapeutic approaches for these patients.

The success of embryo implantation depends on embryonic and endometrial factors. With the innovative breakthrough in preimplantation genetic testing (PGT), abnormal embryos can be discriminated by screening and selection technologies (4). The interventions targeting endometrium have also been more matured in the past decade. For instance, such developments include optimized IVF protocols to synchronize cooperation between the endometrium and the embryo, antibiotic regimens to revert endometrial dysfunction resulting from chronic endometritis (5), and surgical management to rectify congenital uterine malformation and intrauterine adhesion. Nevertheless, endometrium is a complicated tissue composed of various cell types that are dynamically adapted during the window of implantation (WOI) to support the apposition, attachment and invasion of embryo in a synchronic way. The elements involved in this endometrial transformation may open an avenue to elucidate other potential causes and therapeutic targets of implantation failure.

In recent years, increasing accumulated evidence suggests the fundamental role of immunologic factors in embryo implantation. Women with RIF are usually manifested by abnormal immune signaling in the endometrium and peripheral blood, which is particularly characterized by altered Th1/Th2 ratio, and disrupted NK cell and macrophage numbers (6, 7). Additionally, in the endometrium of patients with RIF, the levels of numerous cytokines, chemokines, and their receptors such as IL2, IL6, IFNG, IL17A, IL23A, CCR4, CCR5, CXR3, CCL2, TLR4 are aberrant (8). Further studies are still needed to determine whether these are causative of RIF or just a secondary event. Though the contribution of immune factors in endometrium receptivity can be further explored with the application of RNA-sequencing, the sample sizes in these studies are usually very small, ranging from 9 to 48 samples (9–11), and the immune cell distribution was not analyzed comprehensively. Therefore, it is intriguing to do an in-depth review of the immune network in the endometrium by a systematic analysis of global infiltrating levels of various immune cells and potential immunologic-related markers for the direction of future studies.

Bioinformatics analysis on gene expression microarray has been recommended as a powerful weapon to explore potential biomarkers and signaling pathways in complicated diseases for understanding the pathogenesis and therapeutic development (12). By using microarray data from a publicly accessible database, we further evaluated the differently expressed genes (DEG) in the endometrium of women with RIF and control patients for processing by functional enrichment analysis. To further characterize the features of immune response enriched by aforementioned analysis, we investigated the panorama of immune infiltration by CIBERSORT, and the main findings were preliminarily verified by flow cytometry. With the identification of some immune signaling-related genes, the corresponding ROC curves were plotted for the potential diagnostic values of RIF. Moreover, potential drugs and their potential molecular targets in RIF were identified by cMAP analysis. In this study, we aimed to provide a comprehensive overview of the immune infiltration levels and signaling pathways in women with RIF, and to reveal the potential immune signaling-related biomarkers or therapeutic targets for RIF.

2 Materials and methods

2.1 RIF datasets

Four independent RIF gene expression profiles (GSE103465, GSE111974, GSE58144, and GSE183837) were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE103465, and GSE111974 were exploited as discovery datasets to identify DEGs. GSE58144 and GSE183837 were selected as the validation datasets. GSE103465 is an expression profile based on the GPL16043 platform (GeneChip^® PrimeView™ Human Gene Expression Array (with External spike-in RNAs)) and contains samples of fertile control (n=3) and RIF (n=3). GSE111974 was obtained from the microarray platform of Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381 (Probe Name version), and included 48 endometrium samples (24 control endometrium samples, and 24 RIF endometrium samples). GSE58144 was an expression profile based on the GPL15789 platform (A-UMCU-HS44K-2.0), and contains 71 control endometrium samples, and 43 RIF endometrium samples. GSE183837 was a single-cell RNA-seq (scRNA-seq) dataset of fertile control (n=3) and RIF (n=6), based on the 10X Genomics Chromium platform. The detailed information of these four datasets were listed in Table 1.

TABLE 1

Table 1 Baseline Characteristics.

2.2 Data preprocessing

The “limma” package (13) in R software (version 4.2.0; https://www.r-project.org/) was used to background correction and quantile normalization of all the raw data files, and the expression values were then obtained. The averages of the probe set of values were calculated as the expression values for the same gene with multiple probe sets (14).

2.3 Identification of DEGs

To establish the DEGs between control and RIF samples, the “limma” package of R (13) was used. Moreover, a volcano plot was generated to assess the DEGs. |log2 fold change (FC) |≥ 1, adjusted P value < 0.05 were taken as differentially expressed genes between RIF and control endometrium. Volcano map of DEGs was drawn by GraphPad Prism v8.0 software.

2.4 Functional enrichment of DEGs

The Gene Ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were carried out by using R software and the “clusterProfiler” package (15). In this analysis symbol codes were converted to Entrez ID using Human genome annotation package “org.Hs.eg.db.” The “ggplot2” (16), “pathview” (17) and “gplots” packages of R software to visualize the plots. GSEA software (version 4.1.0, http://www.gsea-msigdb.org/gsea/index.jsp) was used to carry out GSEA analysis (18). The “c2.cp.kegg.v7.5.1.symbols.gmt”, “c5.go.bp.v7.5.1.symbols.gmt”, “c5.go.cc.v7.5.1.symbols.gmt”, “c5.go.mf.v7.5.1.symbols.gmt” was chosen as the reference gene set. Gene set permutations were performed 1,000 times for each analysis. Adjusted P-value < 0.05 and Q value < 0.05 were considered significant for GO and KEGG analysis. The cutoff point of significance was |normalized enrichment score (NES)| > 1, P-value < 0.05, false discovery rate (FDR) Q value < 0.25 for GSEA.

2.5 Protein–protein interaction (PPI) network analysis

The STRING database (available online: http://string-db.org) was performed for PPI network prediction. Cytoscape (v.3.7.2) was used for visual representation and further PPI network experimental studies.

2.6 Hub gene extraction and submodule analysis

Hub genes play an important role in biological processes. Based on the PPI network, hub genes were screened according to network topology. Cytoscape software (version 3.7.2, cytoHubba and MCODE plug-ins) was used to discover the key targets or subnetworks of complex networks (19, 20). The functional annotation and pathway enrichment analysis of hub genes was used by g:Profiler (https://biit.cs.ut.ee/gprofiler/gost) (21).

2.7 Immune infiltration analysis

The different immune cell types of tissues were analyzed by CIBERSORT (22), using normalized gene expression profiles of GSE103465, and GSE111974. LM22, a leukocyte gene signature matrix, contains 547 genes that distinguish 22 human hematopoietic cell phenotypes, including seven T cell types, naïve and memory B cells, plasma cells, NK cells, and myeloid subsets. CIBERSORT results matched ground truth phenotypes in 93% of external datasets of variably purified leukocyte subsets (22). Based on the analysis of LM22, a matrix of 22 kinds of immune cells (naïve B cells, memory B cells, plasma cells, CD8⁺ T cells, CD4⁺ naive T cells, CD4⁺ resting memory T cells, CD4⁺ activated memory T cells, follicular helper T cells, regulatory T cells (Tregs), γδT cells, resting NK cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils) was obtained. CIBERSORT P < 0.05 was used to filter the samples, and the percentage of each immune cell type in the samples was calculated and displayed in a bar plot. The “box” package was used to compare the levels of 22 kinds of immune cells between the two groups. The “ggplot2” package (16) was adopted to visualize the plots.

2.8 ROC analysis

The multivariate modelling with combined selected genes used to identify biomarkers with high sensitivity and specificity for RIF diagnosis by using visualization tool (https://hiplot.com.cn/basic/roc). GSE103465, and GSE111974 datasets were analyzed for training the model, and GSE58144 dataset as validation sample iteratively. The receiver operator characteristic (ROC) curves were plotted and area under curve (AUC) was calculated separately to evaluate the performance of each model using the R packages “pROC” (23). And AUC > 0.7 indicated that the model had a good fitting effect.

2.9 Screening of potential drugs against RIF

Connectivity map (cMap) online tool (http://clue.io) was used to search for compounds that may cause genomic changes similar or opposite to the changes in RIF patients. Compounds with enrichment score >90 or score <−90 were considered significantly positive or negative compounds, respectively.

2.10 Patients’ recruitment for the cohort of validation

The protocol for this study was approved by the Human Research Ethics Committees of Obstetrics and Gynecology Hospital of Fudan University, and written informed consent was obtained from all participants. Human endometrial tissues were collected from women attending the Department of Assisted Reproduction, a dedicated research clinic at the Obstetrics & Gynecology Hospital of Fudan University. Surplus tissue from endometrial biopsies obtained for diagnostic purposes at the Department of Assisted Reproduction was used for this study. The timing of endometrial biopsy was 5 days after ovulation evaluated by ultrasonography, equivalent to LH+7 (hallmark of WOI) in a natural cycle. The definition of RIF was consistent with the studies recruited for bioinformatics. The control group recruited women who had a history of pregnancy but received the assistance of reproductive technology due to obstruction of fallopian tube or male infertility. All the samples were transported to the laboratory on ice in Dulbecco’s modified Eagle’s medium (DMEM)/F-12 (Gibco) for further study.

2.11 Sample preparation and immunocyte isolation

The procedures of sample preparation and immune cell isolation were described with details in our previous studies (24). In brief, the endometrial tissue samples were minced into 1 mm³ pieces on ice and subsequently digested with 1 mg/mL collagenase type IV (Sigma, USA). The tissue pieces were later filtered through a 70 μm cell strainer (Falcon, USA), followed by centrifugation at 400 g for 8 min for the collection of all types of cells. Later, the cells were resuspended in DMEM/F-12 containing 10% fetal bovine serum (FBS; Hyclone, Logan, UT, USA), plated on culture flasks and incubated in a humidified incubator with 5% CO₂ at 37 °C. Finally, after 24 hours, the supernatant was removed and centrifuged at 400 g at 4 °C for flow cytometry analysis.

2.12 Flow cytometry analysis (FCM)

Flow cytometry analysis (FCM) was used to validate the main changes in endometrial immunocytes between women with RIF and controls. All antibodies were from Biolegend (San Diego, CA, USA). The antibodies to detect endometrial NK cells and macrophages included: APC/Cyanine 7 (APC/Cy7)-conjugated anti-human CD45, fluorescein isothiocyanate (FITC)-conjugated anti-human CD3, brilliant violent 605-conjugated anti-human CD56, PE-Cy7-conjugated anti-human CD14, and PE-conjugated anti-human CD16. Staining was performed with the above antibodies (5 μL separately) at room temperature for 30 minutes in the dark. In addition, an isotype IgG antibody of each marker (5 μL separately) was used as the control. Human Trustain FcX (cat. No. 422301; Biolegend, Inc.) was used to block Fc receptors prior to flow cytometry. Subsequently, cells were washed twice and resuspended in PBS for flow cytometry analysis. Samples were analyzed using a CytoFLEX flow cytometer (Beckman Coulter, Inc.), and data were analyzed using FlowJo (version 10.07, FlowJo LLC.).

2.13 Statistical analysis

Statistical analysis was performed using SPSS (version 23.0, Chicago, IL) and GraphPad Prism (version 8, San Diego, CA) software. Continuous variables were expressed as mean ± standard deviation (SD), and differences between the two groups were compared using Student’s t-test for normally distributed variables. A Kruskal–Wallis H test with Dunn’s multiple-comparisons test was used for continuous variables with a non-normal distribution, and the results were presented as the median and interquartile range. The sensitivity and specificity of feature genes to distinguish RIF from controls were assessed using a receiver operating characteristic (ROC) curve. Differences were considered statistically significant when P < 0.05.

3 Results

3.1 Data preprocessing and heat map of top 50 genes in RIF and Ctrl groups

All expression values of GSE103465 and GSE111974 datasets were normalized by principal component analysis (PCA), a method to identify strong patterns in large, complex datasets, which is widely used in data reduction technique captures the essential information in high-dimensional data by identifying principal components that account for most of the variability in the dataset. The data before and after normalization are presented in Figures 1A, B. The heatmap of the top 50 DEGs between RIF and Ctrl groups is shown in Figure 1C.

FIGURE 1

Figure 1 Data Preprocessing and Heat map of top 50 genes in RIF and Ctrl groups. (A, B) PCA analysis of GSE103465, and GSE111974 expression profiles before (A) and after (B) normalization. (C) Expression of top 50 genes in RIF and Ctrl groups was presented by heatmap.

3.2 Identification, PPI network and module analysis of DEGs

The two discovered datasets contained 18301 genes with 27 RIF and 27 Ctrl samples. In total, 188 DEGs (122 genes were downregulated and 66 genes were upregulated) were identified (Figure 2A and Table 2). With PPI enrichment analysis being performed, significant modules were identified (Figures 2B, C). The four MCODE components extracted were found to be mainly associated with NOTCH signaling pathway, cell cycle, translation and ribosome biogenesis.

FIGURE 2

Figure 2 Identification, PPI Network and Module Analysis of DEGs. (A) The 188 DEGs were visualized by volcano map. (B) PPI network constructed with the DEGs from GSE103465, and GSE111974 datasets. (C) The significant module identified from the PPI network using the molecular complex detection (MCODE) method with a score of ≥ 5.0. Different nodes’ color present different functions.

TABLE 2

Table 2 Up and Down DEGs.

3.3 Functional enrichment analysis of DEGs

To further investigate the biological functions of the 188 DEGs, GO analysis, KEGG pathway analysis and GSEA were performed. GO analysis and KEGG pathway analysis demonstrated that DEGs were mainly associated with cell cycle, Wnt/β-catenin signaling, NOTCH signaling, immune system regulation, ATP metabolic process, cellular protein catabolic process, lipid biosynthetic process, and thyroid hormone signaling pathway (Figure 3A). While for the GSEA analysis, activation of immune response and positive regulation of leukocyte migration were remarkably increased in patients with RIF, and cilium movement and steroid hormone biosynthesis were suppressed in patients with RIF (Figure 3B). Besides immune system regulation, rose diagram of the validation dataset (GSE183837) showed the genes in the pathways of cell cycle, the Wnt/β-catenin signaling, NOTCH signaling, ATP metabolic process, cellular protein catabolic process, lipid biosynthetic process, and thyroid hormone signaling pathway were significantly down-regulated in patients with RIF (Figure 3C).

FIGURE 3

Figure 3 Functional Enrichment Analysis of DEGs. (A, B) GO analysis, KEGG pathway analysis (A), and GSEA enrichment analysis (B) of DEGs. (C) Up or down gene numbers related to cell cycle, signaling by Wnt/β-catenin, signaling by NOTCH, immune system regulation, ATP metabolic process, cellular protein catabolic process, lipid biosynthetic process, and thyroid hormone signaling pathway in RIF group from validation dataset (GSE183837) were presented by rose diagram.

3.4 Identification and functional enrichment analysis of hub genes

To further evaluate the significance of the genes in the PPI network of DEGs, the top 20 hub genes were obtained from PPI network by using cytoHubba and presented in Figure 4A. GO analysis, KEGG pathway analysis and REAC analysis demonstrated that hub genes were mainly associated with cell cycle, the Wnt/β-catenin signaling, NOTCH signaling, immune response, cellular protein catabolic process, and thyroid hormone signaling pathway (Figure 4B).

FIGURE 4

Figure 4 Identification and Functional Enrichment Analysis of Hub Genes. (A) The top 20 hub genes were discovered by Cytoscape software (version 3.7.2, cytoHubba plug-ins). (B) GO analysis, KEGG pathway analysis and REAC analysis of hub genes were performed by g:Profiler (https://biit.cs.ut.ee/gprofiler/gost).

3.5 Immune infiltration analysis

More importantly, functional enrichment analysis of DEGs exhibited a significant enrichment of immune response in endometrium from RIF patients. Herein, we further performed the CIBERSORT algorithm to compare the immune cells infiltration level between patients with RIF and the control group (Figure 5A). The percentage of each of the 22 types of immune cells in each sample is shown in the bar plot, and it was intuitively illustrated that most of the endometrial immune cells were NK cells, followed by macrophages, accounting for approximately 30% of total cells which is similar to previous studies (25, 26). The boxplot of the infiltration levels of various immune cells was sorted according to their median proportion in the RIF group, suggesting that NK and macrophages were the top two of the most abundant immune cells. The relative proportions of activated NK cells and macrophages were not significantly different between the two groups, but M0 macrophages were significantly increased in women with RIF (p<0.05) (Figure 5B). These findings were validated by another set of data (GSE183837). As shown in Figure 5C, we found that the proportion of endometrial immune cells in patients with RIF and Ctrl patients also changed by single-cell RNA sequencing (GSE183837). Since NK cells and macrophages were the top two abundant cells in endometrium during the WOI and their crucial roles in implantation is unneglectable, we further measured the proportions of endometrial NK cells and macrophages of patients with RIF and healthy women by FCM to validate the changes. Our data showed that the proportions of macrophages and NK cells were significantly increased in the endometrium of patients with RIF, and the proportions of CD16⁺ macrophage and CD16⁺ NK cells were also elevated (Figures 6A, B), which echoed the results of array in Figure 5 partly. These data implied that RIF was associated with abnormally elevated NK cells and macrophages numbers, as well as increased cytotoxicity.

FIGURE 5

Figure 5 The landscape of immune infiltration between RIF and control samples. (A) The relative percentage of 22 types of immune cells. (B) The difference of immune infiltration between RIF and controls. The control group was marked as blue color and RIF group was marked as red color. (C) Stacked area plot showed the relative percentage of 7 types of immune cells from validation dataset (GSE183837). Data are presented as the median and the interquartile range. Statistical significance (Kruskal–Wallis test with Dunn’s multiple-comparison test): *p < 0.05.

FIGURE 6

Figure 6 The Proportions of Macrophage and NK Cells. (A) The gating strategy for identifying macrophage and NK cells from endometrial leukocytes. (B) The quantitative analysis of the proportions of macrophage, NK cells, CD16⁺ macrophage and CD16⁺ NK cells (n=10) in endometrium from control and RIF patients. Data are presented as the mean ± standard error of the mean. Statistical significance (Student’s t-test): **p < 0.01, ***p < 0.001.

3.6 Expression levels and diagnostic significance of hub genes

Regarding to the importance of Wnt/β-catenin signaling pathway, NOTCH signaling pathway, and immune response in RIF, we identified 6 genes (AKT1, PSMB8, PSMD10, RBX1, PSMC6, and PSMD3) from the hub genes (Figure 7A). As they were all related to these key pathways, 4 of them were remarkably decreased in the endometrium of patients with RIF according to the GSE103465 and GSE111974 datasets (Figure 7B). To determine which hub genes might act as diagnostic markers for patients with RIF, ROC analysis was conducted to explore the diagnostic values of the 6 hub genes for RIF. The results showed that AKT1, PSMB8, and PSMD10 had potential predictive values for RIF in the model of training sets (AUC > 0.7) (Figure 7C). We then validated the model composed of the above 3 qualified genes (AKT1, PSMB8, and PSMD10) in the validation set, the AUC was found to be 0.673 (Figure 7D).

FIGURE 7

Figure 7 Expression Levels and Diagnose Significance of Hub Genes. (A) Venn diagram of common hub genes identified among DEGs, Wnt/β-catenin signaling pathway, NOTCH signaling pathway, immune response, and hub genes. (B) The expression levels of 6 hub genes in GSE103465, and GSE111974 datasets. (C) The GSE103465, and GSE111974 datasets were used to train the diagnostic effectiveness of the biomarkers for RIF by ROC analysis. (D) The GSE58144 dataset was used to validate the diagnostic effectiveness of the biomarkers (AKT1, PSMB8, and PSMD10) for RIF by ROC analysis. Data are presented as the median and the interquartile range. Statistical significance (Kruskal–Wallis test with Dunn’s multiple-comparison test): NS, no significant difference, **p < 0.01, ***p < 0.001.

3.7 Drug-target genes interaction network analysis for RIF

To explore potential therapeutic or induced compounds for RIF, we queried the cMap touchstone database with expression signatures from DEGs. Compounds with enrichment score >90 or <−90, which were regarded as significantly positive or negative compounds, respectively, and their target genes were shown together (Figure 8 and Table 3). Among the drugs with the highest connectivity scores, fulvestrant (estrogen receptor antagonist) was found to be positively correlated with the occurrence of RIF. Among the drugs with the lowest connectivity scores, bisindolylmaleimide-ix (CDK and PKC inhibitor) and JNK-9L (JNK inhibitor) were shown to be negatively correlated with the occurrence of RIF, indicating their potential as therapeutic compounds. These compounds could influence the pathophysiology of RIF by regulating relevant pathways, such as cell cycle, the Wnt/β-catenin signaling, NOTCH signaling, immune response, reproduction, and thyroid hormone signaling pathway. Notably, AKT1, on the hub genes list, was one of bisindolylmaleimide-ix’s targets.

FIGURE 8

Figure 8 Drug-Target genes Interaction Network Analysis. From the left to right of alluvial diagram: the potential pharmaceuticals, target genes and enriched pathways by GO and KEGG analysis.

TABLE 3

Table 3 Information of Top 8 prospective drugs with significant scores for the treatment.

4 Discussion

It is a widely-accepted consensus that the process of embryo implantation is analogous to seed germination, in which the quality of the seed (the embryo) and the receptivity of soil (the endometrium) are both indispensable to successful implantation. In order to promote the success for embryo implantation, the endometrium undergoes a series of changes in functioning signaling pathways and immune cell distribution, subsequently resulting in the maternal immunonutrition and tolerance of the semiallogenic fetus. To comprehend the molecular pathways connected with endometrial receptivity and the foundations of RIF, it is essential to investigate the disruptions in signaling pathways and immune responses in the local endometrium milieu.

The exact mechanisms underlying endometrial immunomodulation by which the embryo is shielded from the maternal immune attack remained poorly understood. As one of the predominant leukocytes in the endometrium, the effects of uNK cells on embryo implantation should not be overlooked (27). In the non-pregnant endometrium, uNK cells are mostly quiescent but they are capable to differentiate throughout the menstrual cycle in preparation for pregnancy (28). When embryo implantation is accomplished, uNK cells contribute to placentation by encouraging trophoblast invasion and spiral artery remodeling to ensure sufficient maternal-fetal material exchange (29). Aberrant uNK activity may produce undesirable outcomes during trophoblast invasion, such as vascular remodeling, local ischemia, and oxidative stress, which are detrimental to implantation (30). According to a latest systematic review and meta-analysis (27), the proportion of CD56⁺ uNK was significantly elevated in patients with RIF when compared with healthy controls. Moreover, most studies found lower expressions of inhibitory receptors and increased expressions of angiogenic factors in uNK cells of women with RIF. These studies may indicate that the assessment of uNK level/activity may be developed as diagnostic tool for RIF. However, a standard reference range would still need to be established prior for use in the clinical setting. In addition to NK cells, macrophages are another significant regulator of the embryo implantation microenvironment. They have been shown to play an important role in implantation, placentation, and embryo development (6). These cells may facilitate trophoblast invasion into spiral arteries, which is related to vascular remodeling at the maternal-fetal interface. Ideal pregnancy outcome calls for appropriate proportion of M1/M2 macrophage ratio as well. However, once distorted, several complications during pregnancy have been shown to occur, including preeclampsia, intrauterine growth restriction, and RIF (31). In line with previous studies, our data validated the higher proportions of uNK cells and macrophages in RIF, notably CD16⁺ NK cells and CD16⁺ macrophages, which revealed stronger cytotoxicity to likely result in maternal immune rejection of the fetus. Nevertheless, further experimental research is still warranted for elucidating the underlying pathophysiology of NK cells and macrophages alterations.

As depicted above, we emphasized the Wnt/β-catenin signaling as a RIF-related signaling pathway by a series of functional enrichment analysis. As shown in previous literature, Wnt/β-catenin signaling acts as a vital part involved in various biological processes ranging from development to early animal evolution and cancer (32). In particular, female reproduction in mammals was shown to be regulated by Wnt4 signaling (33). It has been reported in mice that uterine Wnt/β-catenin signaling is essential for implantation, and the suppression of Wnt4/β-catenin signaling can interferes with the process of embryo implantation (34). On the basis of animal studies, either depletion or overexpression of Wnt/β-catenin in the endometrium would led to implantation failure in mice (35). Meanwhile, Franco et al. demonstrated that the conditional knockout of Wnt4 in the uterus results in infertility and abnormal decidualization in mice (36). Similarly, Zhou et al. discovered during decidualization that EHD1 inhibited the Wnt4/-catenin pathway in patients with RIF by regulating LRP5/6 protein activity via the endocytic route. As a Wnt4 agonist could restore a defective decidualization process, this suggesting that modulation of the EHD1-Wnt4 pathway may be a viable therapeutic target for enhancing endometrial receptivity in women with RIF (37). While studies in human, clinical research has found that the expression of seven Wnt genes (Wnt 2, 3, 4, 5a, 7a, 7b and 10b) in the endometrium during pregnancy occupied a decisive position in the regulation of survival, proliferation and differentiation in endometrial stromal cell (38).

Another highlighted pathway in this current study was the Notch signaling pathway. As an evolutionarily conserved pathway, Notch signaling has been shown to regulate cell proliferation, invasion, differentiation, and apoptosis (39), which contributed significantly in endometrial remodeling. Additionally, during implantation and placentation, Notch signaling is part and parcel of maternal-fetal communication (40). As reported in previous literature, aberrant activation of canonical Notch1 signaling in the mouse uterus could decrease the level of progesterone receptor by hypermethylation that can lead to infertility (41). Likewise, hampered Notch signaling in mice is linked to the development of endometriosis which contributes to dysdecidualization through the down-regulation of FOXO1, which might lead to subsequent implantation failure (42). In a recent study involving women with endometrial disorders including RIF, endometriosis, and PCOS, the expression pattern of NOTCH pathway molecules including NOTCH1, 3, JAG1, 2, and survivin in the mid-luteal phase were distinctly expressed in the patient groups compared to controls (43).

These two pathways have also been demonstrated to regulate immune cell differentiation and activation in autoimmune diseases and cancers, which are commonly proposed to be a potential therapeutic target for immunotherapies (44, 45). Whether the downregulation in Wnt/β-catenin signaling and Notch signaling in women with RIF can cause disrupted endometrial immune cell proportions or not, it is still a promising direction to further dissect the underlying mechanisms of RIF. In this regard, Nevertheless, we identified 3 hug genes (AKT1, PSMB8 and PSMD10) that have strong correlation to the two signaling pathways and immune response when discriminating women with RIF from controls. However, the performance of the predictive model was found to be unsatisfactory. This might be due to the likelihood of small sample size. Hence, further studies with a larger sample size are still needed to explore the potential of these hug genes for elucidating their key pathological pathways of RIF for potential clinical applications. In addition to diagnostic values, compounds (fulvestrant, bisindolylmaleimide-ix and JNK-9L) with potential to interfere the pathogenesis of RIF were identified based on the DEGs between women with RIF and controls. However, the exact effects of these compounds on endometrium and embryo implantation were not validated with additional studies. Further investigations on this aspect are still warranted.

In summary, our study demonstrated a significant difference in the transcriptional profiles between RIF group and the control group (Figure 9). Based on the 122 downregulated and 66 upregulated DEGs, PPI network and module analysis and functional enrichment analysis were performed. Hub genes were further identified by STRING and CytoHubba. Upon our analysis, Wnt/β-catenin signaling pathways, NOTCH signaling pathways, and immune response were the top 3 enriched pathways involved in RIF. With our in-depth analysis, we calculated the immune infiltration level of 22 types of immune cells on validating the increase of NK cells and macrophages in women with RIF. Moreover, the DEGs were assessed to identify potential diagnostic markers and therapeutic compounds for RIF in providing strategies to tackle RIF. However, studies using a larger sample size are still needed to validate the results of the potential genes to potentially be used for diagnosis or prediction of RIF. Another limitation of this study is the lack of functional studies for the potential therapeutic compounds. Nevertheless, our current findings have helped to shed lights for potential directions toward understanding the importance of immunity during RIF.

FIGURE 9

Figure 9 Diagram of the study design. Discover datasets, GSE103465, and GSE111974, contained 18301 genes with 27 RIF and 27 Ctrl samples. In total, 188 DEGs were identified: 122 genes were downregulated and 66 genes were upregulated. Ten important pathways related to RIF were enriched, and were validated by validation dataset (GSE183837). Top 20 hub genes were identified by Cytoscape through PPI network, and three of them may have potential predictive values for RIF in the validation set (GSE58144). Endometrial macrophages’ infiltration level was increased in patients with RIF, compared to the control group, which were performed by CIBERSORT algorithm, and validated by FCM analysis. Bisindolylmaleimide-ix and JNK-9L may be potential therapeutic compounds for RIF, related to the result of cMap analysis.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/; GSE103465, GSE111974, GSE58144, and GSE183837.

Ethics statement

The studies involving human participants were reviewed and approved by Human Research Ethics Committees of Obstetrics and Gynecology Hospital of Fudan University. The patients/participants provided their written informed consent to participate in this study.

Author contributions

Z-ZL and JZ conducted all experiments, prepared the figures and the manuscript, and were responsible for bioinformatics analysis. W-JZ, J-WS, H-LY, S-LY and FX helped sample collection. J-NW helped data analysis. M-QL and TZ designed, initiated and supervised the project and edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding

We are grateful to Dr. Gene Chi Wai Man from Assisted Reproductive Technology Unit, Department of Obstetrics and Gynecology, Faculty of Medicine, Chinese University of Hong Kong for helping edit the manuscript and response latter. This study was supported by the Major Research Program of National Natural Science Foundation of China (NSFC) (No. 92057119, 31970798), the Program for Zhuoxue of Fudan University (JIF157602), the Support Project for Original Personalized Research of Fudan University (IDF157014/002).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest.

The reviewer JZ declared a shared parent affiliation with the author WJZ to the handling editor at the time of the review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Coughlan C, Ledger W, Wang Q, Liu F, Demirol A, Gurgan T, et al. Recurrent implantation failure: definition and management. Reprod BioMed Online (2014) 28(1):14–38. doi: 10.1016/j.rbmo.2013.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Polanski LT, Baumgarten MN, Quenby S, Brosens J, Campbell BK, Raine-Fenning NJ. What exactly do we mean by 'recurrent implantation failure'? a systematic review and opinion. Reprod BioMed Online (2014) 28(4):409–23. doi: 10.1016/j.rbmo.2013.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Pirtea P, De Ziegler D, Tao X, Sun L, Zhan Y, Ayoubi JM, et al. Rate of true recurrent implantation failure is low: results of three successive frozen euploid single embryo transfers. Fertil Steril (2021) 115(1):45–53. doi: 10.1016/j.fertnstert.2020.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Simpson JL, Rechitsky S. Preimplantation genetic testing to reduce preterm births in assisted reproductive technology. Int J Gynaecol Obstet (2020) 150(1):34–40. doi: 10.1002/ijgo.13192

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Devesa-Peiro A, Sebastian-Leon P, Garcia-Garcia F, Arnau V, Aleman A, Pellicer A, et al. Uterine disorders affecting female fertility: what are the molecular functions altered in endometrium? Fertil Steril (2020) 113(6):1261–74. doi: 10.1016/j.fertnstert.2020.01.025

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Mukherjee N, Sharma R, Modi D. Immune alterations in recurrent implantation failure. Am J Reprod Immunol (New York NY 1989) (2022) e13563. doi: 10.1111/aji.13563

CrossRef Full Text | Google Scholar

7. Tremellen KP, Russell P. The distribution of immune cells and macrophages in the endometrium of women with recurrent reproductive failure. II: adenomyosis and macrophages. J Reprod Immunol (2012) 93(1):58–63. doi: 10.1016/j.jri.2011.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Amjadi F, Zandieh Z, Mehdizadeh M, Aghajanpour S, Raoufi E, Aghamajidi A, et al. The uterine immunological changes may be responsible for repeated implantation failure. J Reprod Immunol (2020) 138:103080. doi: 10.1016/j.jri.2020.103080

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Babayeva G, Purut YE, Giray B, Oltulu P, Alakus R, Colakoglu MC. Endometrial CD56+ natural killer cells in women with recurrent implantation failure: An immunohistochemical study. Turk J Obstet Gynecol (2020) 17(4):236–9. doi: 10.4274/tjod.galenos.2020.90359

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Bastu E, Demiral I, Gunel T, Ulgen E, Gumusoglu E, Hosseini MK, et al. Potential marker pathways in the endometrium that may cause recurrent implantation failure. Reprod Sci (Thousand Oaks Calif) (2019) 26(7):879–90. doi: 10.1177/1933719118792104

CrossRef Full Text | Google Scholar

11. Lai Z-Z, Wang Y, Zhou W-J, Liang Z, Shi J-W, Yang H-L, et al. Single-cell transcriptome profiling of the human endometrium of patients with recurrent implantation failure. Theranostics (2022) 12(15):6527–47. doi: 10.7150/thno.74053

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Gauthier J, Vincent AT, Charette SJ, Derome N. A brief history of bioinformatics. Brief Bioinform (2019) 20(6):1981–96. doi: 10.1093/bib/bby063

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Li W, Li K, Zhao L, Zou H. Bioinformatics analysis reveals disturbance mechanism of MAPK signaling pathway and cell cycle in glioblastoma multiforme. Gene (2014) 547(2):346–50. doi: 10.1016/j.gene.2014.06.042

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Walter W, Sanchez-Cabo F, Ricote M. GOplot: an r package for visually combining expression data with functional analysis. Bioinformatics (2015) 31(17):2912–4. doi: 10.1093/bioinformatics/btv300

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics (2013) 29(14):1830–1. doi: 10.1093/bioinformatics/btt285

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinf (2003) 4:2. doi: 10.1186/1471-2105-4-2

CrossRef Full Text | Google Scholar

20. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Reimand J, Kull M, Peterson H, Hansen J, Vilo J. g:Profiler–a web-based toolset for functional profiling of gene lists from large-scale experiments. Nucleic Acids Res (2007) 35(Web Server issue):W193–200. doi: 10.1093/nar/gkm226

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: an open-source package for r and s+ to analyze and compare ROC curves. BMC Bioinf (2011) 12:77. doi: 10.1186/1471-2105-12-77

CrossRef Full Text | Google Scholar

24. Lai ZZ, Ruan LY, Wang Y, Yang HL, Shi JW, Wu JN, et al. Changes in subsets of immunocytes in endometrial hyperplasia. Am J Reprod Immunol (New York NY 1989) (2020) 84(4):e13295. doi: 10.1111/aji.13295

CrossRef Full Text | Google Scholar

25. Nagler A, Lanier LL, Cwirla S, Phillips JH. Comparative studies of human FcRIII-positive and negative natural killer cells. J Immunol (Baltimore Md 1950) (1989) 143(10):3183–91. doi: 10.4049/jimmunol.143.10.3183

CrossRef Full Text | Google Scholar

26. Tang AW, Alfirevic Z, Quenby S. Natural killer cells and pregnancy outcomes in women with recurrent miscarriage and infertility: a systematic review. Hum Reprod (Oxford England) (2011) 26(8):1971–80. doi: 10.1093/humrep/der164

CrossRef Full Text | Google Scholar

27. Woon EV, Greer O, Shah N, Nikolaou D, Johnson M, Male V. Number and function of uterine natural killer cells in recurrent miscarriage and implantation failure: a systematic review and meta-analysis. Hum Reprod Update (2022) 28(4):548–82. doi: 10.1093/humupd/dmac006

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Strunz B, Bister J, Jonsson H, Filipovic I, Crona-Guterstam Y, Kvedaraite E, et al. Continuous human uterine NK cell differentiation in response to endometrial regeneration and pregnancy. Sci Immunol (2021) 6(56):eabb7800. doi: 10.1126/sciimmunol.abb7800

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Huhn O, Zhao X, Esposito L, Moffett A, Colucci F, Sharkey AM. How do uterine natural killer and innate lymphoid cells contribute to successful pregnancy? Front Immunol (2021) 12:607669. doi: 10.3389/fimmu.2021.607669

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Donoghue JF, Paiva P, Teh WT, Cann LM, Nowell C, Rees H, et al. Endometrial uNK cell counts do not predict successful implantation in an IVF population. Hum Reprod (Oxford England) (2019) 34(12):2456–66. doi: 10.1093/humrep/dez194

CrossRef Full Text | Google Scholar

31. Tsao FY, Wu MY, Chang YL, Wu CT, Ho HN. M1 macrophages decrease in the deciduae from normal pregnancies but not from spontaneous abortions or unexplained recurrent spontaneous abortions. J Formos Med Assoc (2018) 117(3):204–11. doi: 10.1016/j.jfma.2017.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Nusse R, Clevers H. Wnt/beta-catenin signaling, disease, and emerging therapeutic modalities. Cell (2017) 169(6):985–99. doi: 10.1016/j.cell.2017.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Vainio S, Heikkila M, Kispert A, Chin N, McMahon AP. Female development in mammals is regulated by wnt-4 signalling. Nature (1999) 397(6718):405–9. doi: 10.1038/17068

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Mohamed OA, Jonnaert M, Labelle-Dumais C, Kuroda K, Clarke HJ, Dufort D. Uterine wnt/beta-catenin signaling is required for implantation. Proc Natl Acad Sci U S A (2005) 102(24):8579–84. doi: 10.1073/pnas.0500612102

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Jeong JW, Lee HS, Franco HL, Broaddus RR, Taketo MM, Tsai SY, et al. Beta-catenin mediates glandular formation and dysregulation of beta-catenin induces hyperplasia formation in the murine uterus. Oncogene (2009) 28(1):31–40. doi: 10.1038/onc.2008.363

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Franco HL, Dai D, Lee KY, Rubel CA, Roop D, Boerboom D, et al. WNT4 is a key regulator of normal postnatal uterine development and progesterone signaling during embryo implantation and decidualization in the mouse. FASEB J (2011) 25(4):1176–87. doi: 10.1096/fj.10-175349

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhou Q, Yan G, Ding L, Liu J, Yu X, Kong S, et al. EHD1 impairs decidualization by regulating the Wnt4/beta-catenin signaling pathway in recurrent implantation failure. EBioMedicine (2019) 50:343–54. doi: 10.1016/j.ebiom.2019.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Bui TD, Zhang L, Rees MC, Bicknell R, Harris AL. Expression and hormone regulation of Wnt2, 3, 4, 5a, 7a, 7b and 10b in normal human endometrium and endometrial carcinoma. Br J Cancer (1997) 75(8):1131–6. doi: 10.1038/bjc.1997.195

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Weijzen S, Rizzo P, Braid M, Vaishnav R, Jonkheer SM, Zlobin A, et al. Activation of notch-1 signaling maintains the neoplastic phenotype in human ras-transformed cells. Nat Med (2002) 8(9):979–86. doi: 10.1038/nm754

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Cuman C, Menkhorst E, Winship A, Van Sinderen M, Osianlis T, Rombauts LJ, et al. Fetal-maternal communication: the role of notch signalling in embryo implantation. Reproduction (2014) 147(3):R75–86. doi: 10.1530/REP-13-0474

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Su RW, Strug MR, Jeong JW, Miele L, Fazleabas AT. Aberrant activation of canonical Notch1 signaling in the mouse uterus decreases progesterone receptor by hypermethylation and leads to infertility. Proc Natl Acad Sci U S A (2016) 113(8):2300–5. doi: 10.1073/pnas.1520441113

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Su RW, Strug MR, Joshi NR, Jeong JW, Miele L, Lessey BA, et al. Decreased notch pathway signaling in the endometrium of women with endometriosis impairs decidualization. J Clin Endocrinol Metab (2015) 100(3):E433–42. doi: 10.1210/jc.2014-3720

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Amjadi F, Salehi E, Zandieh Z, Rashidi M, Taleahmad S, Javedani Masrour M, et al. Comparative evaluation of NOTCH signaling molecules in the endometrium of women with various gynecological diseases during the window of implantation. Iran J Basic Med Sci (2019) 22(4):426–31. doi: 10.22038/ijbms.2019.32961.7874

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Janghorban M, Xin L, Rosen JM, Zhang XHF. Notch signaling as a regulator of the tumor immune response: To target or not to target? Front Immunol (2018) 9. doi: 10.3389/fimmu.2018.01649

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Suryawanshi A, Tadagavadi RK, Swafford D, Manicassamy S. Modulation of inflammatory responses by wnt/beta-catenin signaling in dendritic cells: A novel immunotherapy target for autoimmunity and cancer. Front Immunol (2016) 7:460. doi: 10.3389/fimmu.2016.00460

PubMed Abstract | CrossRef Full Text | Google Scholar