Xiaoxuan Zhao
Yang Zhao
Yuepeng Jiang
Qin Zhang- 1Department of Traditional Chinese Medicine (TCM) Gynecology, Hangzhou Hospital of Traditional Chinese Medicine Affiliated to Zhejiang Chinese Medical University, Hangzhou, China
- 2College of Basic Medicine, Hebei College of Traditional Chinese Medicine, Shijiazhuang, China
- 3College of Pharmacy, Zhejiang Chinese Medical University, Hangzhou, China
Recurrent implantation failure (RIF) is an extremely thorny issue in in-vitro fertilization (IVF)-embryo transfer (ET). However, its intricate etiology and pathological mechanisms are still unclear. Nowadays, there has been extensive interest in cellular senescence in RIF, and its involvement in endometrial immune characteristics during the window of implantation (WOI) has captured scholars’ growing concerns. Therefore, this study aims to probe into the pathological mechanism of RIF from cellular senescence and investigate the correlation between cellular senescence and endometrial immune characteristics during WOI based on bioinformatics combined with machine learning strategy, so as to elucidate the underlying pathological mechanisms of RIF and to explore novel treatment strategies for RIF. Firstly, the gene sets of GSE26787 and GSE111974 from the Gene Expression Omnibus (GEO) database were included for the weighted gene correlation network analysis (WGCNA), from which we concluded that the genes of the core module were closely related to cell fate decision and immune regulation. Subsequently, we identified 25 cellular senescence-associated differentially expressed genes (DEGs) in RIF by intersecting DEGs with cellular senescence-associated genes from the Cell Senescence (CellAge) database. Moreover, functional enrichment analysis was conducted to further reveal the specific molecular mechanisms by which these molecules regulate cellular senescence and immune pathways. Then, eight signature genes were determined by the machine learning method of support vector machine-recursive feature elimination (SVM-RFE), random forest (RF), and artificial neural network (ANN), comprising LATS1, EHF, DUSP16, ADCK5, PATZ1, DEK, MAP2K1, and ETS2, which were also validated in the testing gene set (GSE106602). Furthermore, distinct immune microenvironment abnormalities in the RIF endometrium during WOI were comprehensively explored and validated in GSE106602, including infiltrating immunocytes, immune function, and the expression profiling of human leukocyte antigen (HLA) genes and immune checkpoint genes. Moreover, the correlation between the eight signature genes with the endometrial immune landscape of RIF was also evaluated. After that, two distinct subtypes with significantly distinct immune infiltration characteristics were identified by consensus clustering analysis based on the eight signature genes. Finally, a “KEGG pathway–RIF signature genes–immune landscape” association network was constructed to intuitively uncover their connection. In conclusion, this study demonstrated that cellular senescence might play a pushing role in the pathological mechanism of RIF, which might be closely related to its impact on the immune microenvironment during the WOI phase. The exploration of the molecular mechanism of cellular senescence in RIF is expected to bring new breakthroughs for disease diagnosis and treatment strategies.
Introduction
Recurrent implantation failure (RIF) is usually defined as implantation failure in at least three consecutive attempts of in-vitro fertilization (IVF), in which one to three high-quality embryos are transferred in each cycle (1, 2). Despite tremendous advances in reproductive medicine, it remains an ongoing conundrum. About 15% of patients seeking IVF treatment are reported to be afflicted by RIF, which is a huge distress and frustration for both patients and clinicians (3). Statistics have shown that approximately 40% of euploid blastocysts fail to implant in transfers, which suggests the desynchronization between the embryo and the endometrium in the window of implantation (WOI) as a potential cause of RIF (4, 5). During the WOI period, the endometrium becomes poised to transition to a pregnant state, featured as differentiation of endometrial stromal cells (EnSCs) into decidual cells (DCs) and change of the influx and function of immunocytes in the local endometrium (6). Perturbations of these key cellular and molecular biological events tend to induce the breakdown of the feto–maternal interface and RIF (7, 8), but the underlying mechanisms are poorly understood. Accumulated evidence has suggested a conspicuous link between cellular senescence in peri-implantation endometrium and RIF (9, 10). Cellular senescence is a state of permanent cell cycle arrest and manifests with a prominent secretion of various bioactive molecules, including reactive oxygen species, pro-inflammatory cytokines, chemokines, and growth factors, called senescence-associated secretory phenotype (SASP) (11, 12). SASP secreted by senescent cells creates a long-lasting and highly disordered pro-inflammatory response, which also attracts multiple immunocytes, thus jointly facilitating an unfavorable microenvironment for embryo implantation (13, 14). Single-cell transcriptomics revealed that decidualized assembloids harbored senescent subpopulations, and senescence in the stroma calibrated the emergence of anti-inflammatory decidual cells and pro-inflammatory senescent decidual cells, which controlled endometrial fate decisions at implantation and was closely correlated with RIF (15). Moreover, intricate immune regulation is one of the most crucial aspects for the successful implantation of the hemiallogenic embryo (16, 17). Conversely, inappropriate immunocytes as well as immune function are implicated in RIF as evidenced by increasing high-quality studies (18). Furthermore, to date, compelling evidence has also indicated that senescent cells are capable of regulating the immune microenvironment and contributing to advances in pathological mechanisms in various diseases (19–21). However, the understanding of cellular senescence, immune infiltration landscape, and the linkages between the two in the context of RIF disease is still close to blank.
Nowadays, the revolutionary development of microarray technology and bioinformatics greatly facilitates the development of biomedicine. A tremendous amount of high-throughput data is piling up in public databases, which greatly helps with uncovering the potential etiopathogenesis and identifying candidate targets for drug design (22). Machine learning has recently been widely applied to learn the representation of high-dimensional features derived from gene expression data on account of its powerful capabilities in classification (23, 24). The ingenious combination of bioinformatics analysis and machine learning is a creative and crucial way to establish novel diagnostic models and understand pathological mechanisms at the molecular level, which is in line with the latest research trends. However, the application of this method in RIF is still blank, and much of the potential valuable information remains to be uncovered.
In this study, RIF raw microarray data were acquired from the Gene Expression Omnibus (GEO) database, and cellular senescence-related genes were obtained from the Cell Senescence (CellAge) database for systematic analysis. Firstly, the enrichment analysis of key module genes obtained by the weighted gene co-expression network analysis (WGCNA) was carried out. Subsequently, differentially expressed genes (DEGs) related to cellular senescence were screened. Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and gene-set enrichment analysis (GSEA) were conducted to further explore the biological mechanisms that the cellular senescence-associated DEGs were involved in. Moreover, we screened and validated the signature genes by machine learning algorithms of support vector machine-recursive feature elimination (SVM-RFE), random forest (RF), and artificial neural network (ANN). After that, we analyzed the difference in immune landscape between the RIF group and the normal group by assessing infiltrating immunocytes, immune function, and the expression profiling of HLA (human leukocyte antigen) genes and immune checkpoint genes. Moreover, we evaluated the correlation between the signature genes and immune landscape in RIF. Last but not least, we clustered RIF samples into two distinct subtypes according to the signature genes and compared the endometrial immunity between subtypes. The workflow chart is shown in Figure 1. To the best of our knowledge, this study is the first to explore the pathological mechanisms of RIF in terms of cellular senescence affecting the endometrial immune landscape during the WOI stage by integrated bioinformatics analysis and machine learning, which provides new insights into the treatment strategy for RIF and offers valuable insulation and foundation for further innovative studies on RIF.
Figure 1 The flow diagram of the study.
Materials and methods
Data collection
The microarray datasets were systematically extracted from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (25) with the keyword “Recurrent implantation failure.” Datasets that met the following inclusion criteria were included: 1) Homo sapiens; 2) expression profiling by array; 3) the experiment included patients with RIF and fertile controls; 4) the sample size was at least 10 people, with at least 5 patients in each group; and 5) the sample was from the endometrium during WOI. Finally, three datasets were included (GSE26787, GSE111974, and GSE106602), of which GSE26787 and GSE111974 were taken as the training set, and GSE106602 was taken as the testing set. The details of the three datasets are shown in Table 1. Furthermore, 279 cellular senescence-associated genes were downloaded from the CellAge database (https://genomics.senescence.info/cells/) (26).
Table 1 Basic information of the included dataset.
Data preprocessing and normalization
The raw data were downloaded from the GEO database, and then preprocessed and normalized by R statistical software (version 4.1.2, https://www.r-project.org/) and Bioconductor analysis tools (http://www.bioconductor.org/). The “affy” R language package was applied to conduct RMA background correction, complete log2 transformation, quantile normalization, and median polish algorithm summarization. Probes without matching gene symbols were excluded. For multiple probes mapped to the same gene, the mean value was taken as the final expression value. The results are shown in box plots and three-dimensional PCA cluster diagrams before and after normalization (Supplementary Figures 1A–D).
Construction of the co-expression network and hub module identification by WGCNA
The WGCNA is a widely used method to uncover critical interacted genetic modules and key genes by linking gene networks to clinical traits. In this study, the WGCNA co-expression system was established by using the “WGCNA” package in R software (27) with DEGs from the GSE26787 and GSE111974 datasets, with normal control and RIF as clinical features. First, genes with variation higher than 25% across samples in the combined dataset were selected as the input dataset for the subsequent WGCNA. Then, the outlier cases were removed by hierarchical clustering analysis with the “goodSamplesGenes” function. After that, the appropriate soft threshold was determined by using the pickSoftThreshold function and validated by the correlation between k and p(k). Subsequently, the correlation matrix was converted into an adjacency matrix, which was further processed into a topological overlap matrix (TOM). The dynamic tree cutting approach was performed to identify various modules. The relationship between these modules and RIF was investigated. Finally, the module with the greatest Pearson correlation coefficient was picked for further investigation.
Identification of cellular senescence-associated DEGs
We used the processed data to filter DEGs by using the “limma” packages of R software (version 4.1.2) with the screening criteria of adjusted P-value <0.05 and |log2 fold change (FC)| >0.5 (28). The Venn online mapping tool was used to screen for intersections of DEG and cellular senescence-associated genes, namely, cellular senescence-associated DEGs, which were visualized with heatmaps and volcano plots by the “ggplot2” package in R software (29).
Function enrichment analysis of cellular senescence-associated DEGs
The cellular senescence-associated DEGs were imported into the WebGestalt website (http://www.webgestalt.org/) (30) for GO and KEGG pathway enrichment analysis. The species was selected as “Homo sapiens,” and the reference set was selected as genome protein-coding. Items with a P-value <0.05 were displayed in a bubble plot by the “ggplot2” package in R software. Moreover, GSEA (version 4.1.0) on cellular senescence-associated DEGs was conducted by using the “c2.cp.kegg.v7.0.symbols.gmt” and “c5.go.v7.4.symbols” gene set derived from the Molecular Signatures Database (MSigDB; version 7.1) as reference (31). The threshold for significant terms was adjusted q value <0.05.
Screening cellular senescence-associated signature genes of RIF by machine learning
The SVM-RFE algorithm and RF were utilized to screen the cellular senescence-associated signature genes in RIF. SVM-RFE is a sequence backward selection algorithm based on the maximum margin principle of SVM, which has superior classification performance for high-dimensional datasets (32). It was first proposed in gene selection (33) and was applied to gene selection in cancer classification. After that, the SVM-RFE algorithm was further improved to improve its performance and efficiency, and it has been widely used in gene expression data analysis, protein function prediction, image detection, and other fields (34). In this study, the SVM-RFE algorithm was implemented by using the package of “e1071,” “kernlab,” and “caret” in R software for feature dimensionality reduction. The RF algorithm is an ensemble method that combines many decision trees and makes a single decision on behalf of the ensemble by combing the results of multiple classifiers together (35). Each decision tree in the forest is built by using the bootstrap technique to select various samples from the original dataset and then training it with a feature set chosen by the random bagging mechanism (36). Decisions made by a large number of distinct individual trees are then voted on, and the class with the most votes as a result of the voting is assigned as the class prediction (37). Here, we used the RF algorithm to predict RIF with the input of cellular senescence-associated signature genes by the “randomForest” package in R software. In addition, we constructed an ANN model for the feature genes obtained from the above method according to the gene score by using the packages of “neuralnet” and “neuralnettools” in R software. The artificial neural network can simulate the structure and function of the brain neural network and deduce a set of classification rules from a set of disordered and irregular data, so as to realize the correct classification and construct a high-accuracy diagnosis model (38). Furthermore, the receiver operating characteristic (ROC) curve was utilized to evaluate the accuracy of the ANN model in the training and testing sets.
Finally, feature genes screened by the above machine learning methods were regarded as signature genes of RIF from the aspect of cellular senescence, namely, RIF signature genes. We observed the expression of these eight genes in the training set (GSE26787 and GSE111974) and the testing set (GSE106602). Furthermore, the ROC curve was constructed by the “pROC” package in R software to evaluate the prediction accuracy of these feature genes.
The tissue localization and function of RIF signature genes
The information of tissue localization of RIF signature genes was obtained through the BioGPS (building your own mash-up of gene annotations and expression profiles) website (http://biogps.org) (39) and the Human Protein Atlas (https://www.proteinatlas.org/) (40). Moreover, to further perceive the functions of these genes, we searched the databases of BioGPS, GeneCards (https://www.genecards.org/) (41), Alliance of Genome Resources (https://www.alliancegenome.org/) (42), and UniProt (https://www.uniprot.org/) (43).
Immune infiltration landscape analysis of RIF and its correlation with cellular senescence-associated signature genes
Single-sample gene-set enrichment analysis (ssGSEA) is an extension of the GSEA method, which allows to analyze the pathways enriched by genes in each sample, thereby analyzing the activation degree of specific pathways (44, 45). Here, ssGSEA was performed by using the “GSVA” package in R software to explore the immune cell-related pathways in RIF patients and healthy controls. Moreover, the expression of HLA molecules and immune checkpoints was quantified and compared between groups. Furthermore, correlation analysis between cellular senescence-associated signature genes and immunocytes, immune functions, HLA, and immune checkpoints in RIF was evaluated by Spearman correlation analysis with the “corrplot” package in R software. P <0.05 was considered statistically significant.
Consensus clustering analysis
Consensus clustering analysis was carried out by using the “ConsensusClusterPlus” package of R software on the basis of the expression profiling of cellular senescence-associated signature genes (46). The above step was implemented for 1,000 iterations for guaranteeing the robustness of classification. The Euclidean distance between specimens was determined. Consensus cluster analysis was conducted with a maximum number of clusters of 9, 50 repeats, a sample proportion of 0.8, and a ratio of features to samples of 1. The pam cluster algorithm was selected, and cluster analysis was conducted by using Euclidean as a distance function. After the completion of cluster analysis, PCA analysis was conducted by using the “limma” package, and classification dot plots were drawn to verify the cluster results. In addition, the expression of RIF-specific cell senescence-related genes, immunocytes, immune functions, HLA, and immune checkpoints in the different subtypes were analyzed.
Construction of the “KEGG pathway–RIF signature genes–immune landscape” association network
In order to intuitively understand the molecular and biological processes involved and communicated by RIF signature genes, we constructed the “KEGG pathway–RIF signature genes–immune landscape” association network. First, the cross-talk relationship between the KEGG pathways enriched by RIF signature genes was obtained through the XTalkDB website (http://www.xtalkdb.org/contactus) (47). In addition, we sorted out the correspondence between the RIF signature genes and the significantly correlated immunocytes, immune function, HLA moleculars, and immune checkpoints. Next, we used Cytoscape 3.7.1 software to visualize the “KEGG pathway–RIF signature genes–immune landscape” association network.
Results
Establishment of a co-expression network and hub module identification
Module detection was performed by hierarchical clustering and dynamic tree cut functions (Supplementary Figures 2A–C). The soft thresholding power was set at 4 according to the scale independence and mean connectivity values. A total of 17 modules were divided by WGCNA and were identified by different colors. Among these modules, the green module showed the most significant correlation with RIF (r = −0.67, P = 0.003) (Supplementary Figure 2D). Therefore, we focused on the biological functions of the genes in the green module through GO and KEGG analysis. The results of GO annotation suggested that these genes were mostly enriched in processes involved in cellular response to interleukin-4, regulation of T-cell chemotaxis, positive regulation of lymphocyte chemotaxis, etc. As for the KEGG pathway, the majority of genes were enriched in the Wnt signaling pathway, p53 signaling pathway, and longevity regulating pathway, which were closely related to the regulation of cell fate decision. Based on the above results, it intimately implied that the pathological mechanism of RIF was bound up with immune abnormalities and the dysregulation of cell fate (Supplementary Figures 2E, F).
Identification and integrative analysis of DEGs
After preprocessing and normalization of the included datasets, differential expression analysis was performed. A total of 1,919 DEGs were obtained with the screening conditions of P-value <0.05 and |log2 FC| >0.5. Then, 25 overlapped genes were obtained by the Venn diagram after intersecting the 1,919 DEGs with 279 cellular senescence-associated genes derived from the CellAge database (Figure 2A). Moreover, these 25 genes were taken as a new gene set for further study, namely, cellular senescence-associated DEGs. Next, we observed the difference in the mRNA expression profiling of these 25 genes between the RIF and control groups. Among the 25 cellular senescence-associated DEGs, a total of 12 genes were upregulated and 13 genes were downregulated (Figures 2B–D). Subsequently, we explored the correlations between cellular senescence-associated DEGs through the Spearman correlation test. The results corroborated that CKB was significantly positively correlated with GNG11 (r = 0.81, P < 0.05), and SLC16A7 was significantly negatively correlated with CKB (r = −0.77, P < 0.05) (Figure 2E).
Figure 2 Identification of differentially expressed genes related to cellular senescence. (A) Venn diagram of differentially expressed genes and genes related to cellular senescence. (B, C) Volcano plot and heatmap visualized the differentially expressed genes (DEGs) related to cellular senescence. In the volcano plot, each dot represents a gene. The red plot points represent upregulated genes, and the blue plot points represent downregulated genes, and in the heatmap, each row represents a DEG related to cellular senescence, and each column represents a sample. (D) The difference in the mRNA expression profiling of cellular senescence-associated DEGs between the recurrent implantation failure (RIF) and control groups. (E) The correlations between cellular senescence-associated DEGs in the RIF group. * represents P <0.05 compared with the control group, ** represents P <0.01 compared with the control group, and *** represents P <0.001 compared with the control group.
Functional enrichment analysis
To comprehend the biological processes these cellular senescence-associated DEGs were involved in, we performed the GO and KEGG enrichment analysis and GSEA method. As shown in the GO enrichment analysis, those genes were mainly enriched in the biological process of cell growth, regulation of growth, negative regulation of canonical Wnt signaling pathway, etc. As for the cellular component, these genes were mainly enriched in insulin-like growth factor binding protein complex, growth factor complex, and cell–substrate junction. Furthermore, as for the molecular function, these genes were mainly enriched in insulin-like growth factor binding, protein serine/threonine/tyrosine kinase activity, MAP kinase activity, etc. (Supplementary Figure 3A).
The results of the KEGG enrichment analysis signified that the 25 cellular senescence-associated DEGs were significantly activated in aging-related pathways, such as mTOR signaling pathway, cellular senescence, etc., and immune inflammation-related pathways, such as Toll-like receptor signaling pathway, TNF signaling pathway, etc. The intricate mapping relationship between genes and pathways was also displayed in the Sankey diagram in which MAP2K1 corresponded to the most KEGG pathways, indicating a crucial role in communicating the complicated cross-talk between different pathways (Supplementary Figure 3B). Furthermore, we applied the GSEA method to appraise the differences in activation of GO enrichment and KEGG pathway between the RIF and fertile control samples. Our results revealed that GO terms of immune response, leukocyte migration, and T-cell activation were significantly activated in fertile control samples, while cytoplasmic pattern recognition receptor signaling pathway was activated in RIF samples (Supplementary Figures 3C, D). Meanwhile, the KEGG pathways of chemokine signaling pathway, cytokine–cytokine receptor interaction, and natural killer cell-mediated cytotoxicity pathway exhibited increased activation in fertile control samples, while abc transporters and steroid biosynthesis pathways were activated in RIF samples (Supplementary Figures 3E, F).
Screening cellular senescence-associated signature genes of RIF by machine learning
To clarify the diagnostic value of cellular senescence-associated DEGs, we performed machine learning of SVM-RFE algorithms and RF algorithms and obtained two and eight feature genes, respectively (Supplementary Figures 4A–C). After combining the results of the two algorithms, a total of eight cellular senescence-associated signature genes of RIF were obtained, namely, LATS1, EHF, DUSP16, ADCK5, PATZ1, DEK, MAP2K1, and ETS2. Moreover, these eight genes were used for constructing the neural network, and the results showed that the eight genes could distinguish the control samples from the RIF samples well, with an accuracy rate of 100%. Furthermore, ROC was constructed, and the results showed that the AUC of the training set and the testing set was 1.000 and 0.678, respectively, which demonstrated the high accuracy of the ANN model (Figures 3A, B). Thus, the eight genes were named as RIF signature genes. Among these genes, LATS1, EHF, ADCK5, and PATZ1 were highly expressed in RIF samples (P < 0.001), and DUSP16, DEK, MAP2K1, and ETS2 declined significantly in RIF samples (P < 0.001) (Figure 3C). To further evaluate the diagnosis value of the eight feature genes screened by the above machine learning methods, we constructed the ROC curves of these eight RIF signature genes and calculated the AUC. The results elucidated that the AUC ranged from 0.793 to 0.994 (Figures 3D–K), indicating a good diagnostic performance of these genes. In addition, we performed external validation with the testing dataset of GSE106602. The mRNA expression profiling of the eight RIF signature genes was evaluated, and the result indicated that the mRNA levels of ADCK5, DUSP16, EHF, ETS2, MAP2K1, and PATZ1 were significantly different between the RIF and fertile control groups (ADCK5, EHF: P < 0.05; DUSP16, ETS2, MAP2K1, PATZ1: P < 0.001) (Figures 4A–H). Then, we constructed the ROC curves and calculated the AUC of these eight genes. The AUC varied from 0.543 to 0.954 (Figures 4I–P), which also confirmed the excellent diagnostic performance of these genes.
Figure 3 Diagnostic value of cellular senescence-associated DEGs in RIF. (A) The neural network model: I1–I8 are the input layers (the score and weight of eight RIF signature genes), H1–H5 are the hidden layers, and O1–O2 are the output layers (sample attributes). (B) The receiver operating characteristic (ROC) curves for evaluating the diagnostic efficacy of the neural network model in the GSE26787 and GSE111974 (training set) and GSE106602 (testing set). (C) The mRNA expression profiling analysis of the eight RIF signature genes. (D–K) The ROC curves of the RIF signature genes. *** represents P <0.001 compared with the control group.
Figure 4 External validation to further test the diagnostic performance of the RIF signature genes. (A–H) Expression differences of the RIF signature genes among different groups in the testing set. (I–P) The ROC curves of the RIF signature genes in the testing set.
The tissue localization and function of RIF signature genes
To further perceive the tissue localization and function of these genes, we searched the databases of BioGPS, Human Protein Atlas, GeneCards, Alliance of Genome Resources, and UniProt to acquire comprehensive information, as shown in Table 2. From this table, we could discover that multiple genes were expressed in the uterine tissue and immunocytes. This was the premise that they could participate in cellular senescence, regulate immunization activities, and thus, affect endometrial receptivity.
Table 2 Immune localization of RIF signature genes.
Immune infiltration landscape of RIF and its correlation with RIF signature genes
To fully grasp the immune landscape of RIF, we analyzed differences in immune cells, immune function, and expression profiling of the HLA gene set and immune checkpoints between patients with RIF and fertile controls. Regarding the types of immunocytes, the levels of CD8+ T cells, immature dendritic cells (iDCs), macrophages, neutrophils, plasmacytoid dendritic cells (pDCs), Th1 cells, Th2 cells, and Treg in the RIF group were significantly lower than those in the control group (P < 0.05). Meanwhile, the T helper cells were significantly increased in the RIF group (P = 0.047) (Figure 5A). Furthermore, the changes of macrophages, pDCs, neutrophils, Th1 cells, and Treg cells in RIF in the testing set were consistent with the training set (Supplementary Figure 5A). The above data indicated the heterogeneity in immunocytes between control and RIF samples. In addition, Spearman correlation analysis implied that Th1 cells were significantly positively correlated with Th2 cells in RIF samples (r = 0.66, P < 0.05), and B cells were significantly negatively correlated with Th2 cells (r = −0.32, P < 0.05) (Figure 5B). Moreover, we further investigated the interactions of RIF signature genes with infiltrated immunocyte in RIF samples. We obtained the most correlated pairs of the RIF signature gene–immunocyte, including ADCK5 and neutrophils (r = −0.40, P < 0.05), DEK and B cells (r = −0.49, P < 0.01), DUSP16 and B cells (r = −0.45, P < 0.05), EHF and Th2_cells (r = −0.40, P < 0.05), ETS2 and follicular helper T cells (Tfh) (r = 0.59), Treg (r = 0.62, P < 0.001), LATS1 and Treg (r = −0.51, P < 0.01), MAP2K1 and Tfh (r = 0.54, P < 0.01), and PATZ1 and Tfh (r = −0.49, P < 0.01) (Figures 5C–J).
Figure 5 Immune cell infiltration landscape of RIF and its correlation with the RIF signature genes. (A) The distribution of immunocytes between the RIF and control samples. (B) The correlation heatmap showed the correlation between different immunocytes in RIF samples. (C–J) The lollipop chart showed the correlation between RIF signature genes and immunocytes in RIF samples.
When it came to the comparison of immune functions between groups, the results substantiated that the functions of APC co-stimulation, CCR, checkpoint, cytolytic activity, inflammation-promoting, MHC class I, parainflammation, T-cell co-inhibition, T-cell co-stimulation, and type I IFN response in the RIF group were significantly lower than those in the control group (P < 0.05) (Figure 6A). Moreover, the changes of APC co-stimulation, CCR, checkpoint, parainflammation, and T-cell co-stimulation in RIF in the testing set were consistent with the training set (Supplementary Figure 5A). The above data corroborated the heterogeneity in immune function between the control and RIF samples. In addition, Spearman correlation analysis revealed that the function of checkpoint was positively correlated with CCR (r = 0.87, P < 0.05) and inflammation-promoting (r = 0.87, P < 0.05) in RIF samples (Figure 6B). Furthermore, we obtained the immune functions that were most associated with the RIF signature genes by Spearman correlation analysis. The most strongly associated pairs included ADCK5 and T_cell_co-stimulation (r = 0.53, P < 0.01), EHF and parainflammation (r = −0.55, P < 0.01), ETS2 and CCR (r = 0.60, P < 0.001), APC co-stimulation (r = 0.63, P < 0.001), LATS1 and parainflammation (r = −0.40, P < 0.05), MAP2K1 and type II IFN response (r = 0.50, P < 0.01), and PATZ1 and APC co-stimulation (r = −0.59, P < 0.001) (Figures 6C–J).
Figure 6 Immune function landscape of RIF and its correlation with RIF signature genes. (A) The distribution of immune functions between RIF and control samples. (B) Correlation heatmap showed the correlation between different immune functions in RIF samples. (C–J) The lollipop chart showed the correlation between RIF signature genes and immune functions in RIF samples.
When it came to the HLA gene set expression profile, the expression levels of HLA-A, HLA-G, HLA-J (P < 0.01), HLA-DPB1 (P < 0.05), HLA-DPA1, HLA-DRA, and HLA-F (P < 0.001) in the RIF group were all significantly decreased (Figure 7A). Moreover, the expression profile changes of HLA-A and HLA-J in RIF in the testing set were consistent with the training set (Supplementary Figure 5B). The results of the Spearman correlation analysis illuminated that ETS2 was significantly positively correlated with HLA-A (r = 0.63, P < 0.05), and EHF was significantly negatively correlated with HLA-DPB1 (r = −0.64, P < 0.05) in the RIF group (Figure 7B).
Figure 7 Human leukocyte antigen (HLA) and immune checkpoint of RIF and its correlation with RIF signature genes. (A, C) The distribution of HLA and immune checkpoint between RIF and control samples. (B, D) The correlation heatmap showed the correlation between RIF signature genes, HLA genes, and immune checkpoints in RIF samples. * represents P < 0.05 compared with the control group,**represents P < 0.01 compared with the control group, ***represents P < 0.001 compared with the control group.
As for the immune checkpoints, the mRNA expression profiles of ICOS (P < 0.01), CD70 (P < 0.01), LAIR1 (P < 0.001), TNFRSF8 (P < 0.001), TNFSF15 (P < 0.001), CD40 (P < 0.001), TNFRSF4 (P < 0.05), TNFSF14 (P < 0.01), HAVCR2 (P < 0.01), BTLA (P < 0.01), CD28 (P < 0.05), CD48 (P < 0.001), CD40LG (P < 0.001), TNFRSF18 (P < 0.01), and NRP1 (P < 0.01) in the RIF group were significantly lower than those in the control group (P < 0.05), and the mRNA expression of CD276 (P < 0.01), VTCN1 (P < 0.01), and BTNL2 (P < 0.05) was significantly increased in the RIF group when compared with the control group (P < 0.05) (Figure 7C). Furthermore, the expression profile changes of ICOS, LAIR1, TNFSF15, CD40, TNFSF14, and CD40LG in RIF in the testing set were consistent with the training set (Supplementary Figure 5C). The results of the Spearman correlation analysis displayed that DEK was significantly positively correlated with NRP1, ETS2, and CD200 (r = 0.70, P < 0.05) in the RIF group, and EHF was significantly negatively correlated with TNFRSF8 and CD244 (r = −0.64, P < 0.05) in the RIF group (Figure 7D).
Immune infiltration landscape of distinct subtypes
Through consensus clustering analysis, we clustered RIF samples into two subtypes based on the mRNA expression profiling of the eight RIF signature genes, namely, subtype A (n = 6) and subtype B (n = 23) (Supplementary Figures 6A–E). Except for LATS1, ADCK5, and DEK, the expression differences of the remaining five RIF signature genes among the groups were statistically significant (P < 0.05) (Supplementary Figures 6F, G). In order to figure out the differences of immune infiltration landscape between the two subtypes, the ssGSEA method was utilized to quantify the infiltrating level of immunocytes and immune function, and the expression profiling of the HLA and immune checkpoint gene sets was also evaluated. The results showed that various infiltrating immunocytes in cluster B were significantly reduced when compared with cluster A, such as anchorage-dependent cells (aDCs), CD8+ T cells, neutrophils, pDCs, T helper cells, etc. (P < 0.05) (Figure 8A). The immune functions in cluster B, such as APC co-inhibition, APC co-stimulation, CCR, checkpoint, etc., were significantly decreased (P < 0.05) (Figure 8B). In addition, the mRNA expression profiles of HLA-A, HLA-DPA1, HLA-DPB1, etc. in cluster B were significantly decreased, and HLA-DOB was significantly increased when compared with cluster A (P < 0.05) (Figure 8C). Furthermore, various immune checkpoints in cluster B were also significantly decreased, such as TNFRSF9, CD70, LAIR, etc. (Figure 8D). All the above results indicated that different expression patterns of RIF signature genes exerted a distinct impact on the immune infiltration landscape in the endometrium during WOI.
Figure 8 Immune-related infiltration landscape in RIF by unsupervised clustering based on the eight RIF signature genes. (A) The abundance differences of different immunocytes between the two subtypes. (B) The activity differences of different immune functions between the two subtypes. (C) The abundance differences of different HLA between the two subtypes. (D) The abundance differences of different immune checkpoints between the two subtypes. * represents P <0.05 compared with cluster A, ** represents P <0.01 compared with cluster A, and *** represents P <0.001 compared with cluster A.
Discussion
Embryo implantation remains the rate-limiting step for IVF, and endometrial receptivity during the WOI has received extensive attention, which has been considered as a crucial impetus for deciphering the pathological mechanism of RIF (48, 49). Currently, in addition to the recognized immune factors, accumulated evidence has also highlighted the significant influence of cellular senescence on endometrial receptivity (50, 51). For instance, Deryabin et al. uncovered that premature senescence of EnSC could alter the “meta-signature” of human endometrial receptivity and affect embryo invasion based on an in-vitro implantation model (52). Moreover, studies by Chen et al. corroborated that premature aging of the endometrium existed in young women with RIF and was closely related to implantation failures (53). Furthermore, cellular senescence has been demonstrated to be closely related to the immune microenvironment (54). Senescence subjects the immune system to constant immune stressors and inflammatory assaults, which in turn promotes immune senescence and is intimately tied to diverse pathological changes (55). Now, the association of cellular senescence and immune regulation has attracted growing attention in various diseases. Relevant findings suggest that SASP factors contribute to disease progression by regulating immunity and doing harm to the tissue homeostasis (56, 57). However, the effect of cellular senescence on the immune system of RIF has not been elucidated. Studies have confirmed that multiple SASP factors promote the activation of molecular signaling cascades that are associated with embryo rejection (58–60). These discoveries provide us with a sufficient basis to demonstrate the potential linkage between cellular senescence and immune regulation in RIF.
Nowadays, the current diagnosis of the disease is retrospectively established based on repeated failed attempts in IVF-ET, which pose great challenges to clinical strategy formulation for precise prevention and treatment of RIF. In this study, we interpreted the pathological mechanism of RIF from a new perspective and explored a new strategy for the diagnosis of RIF based on bioinformatics analysis combining machine learning. Overall, we probed the involvement of cellular senescence in the pathological mechanism of RIF, revealed the immune infiltration landscape in the endometrium during WOI, and analyzed the correlation between cellular senescence and immune dysregulation. After a series of explorations, we obtained the following significant findings: 1) through the WGCNA method, the genes in the most relevant modules of RIF mainly affected immune functions (T-cell chemotaxis, positive regulation of lymphocyte chemotaxis) and cell fate decision (Wnt signaling pathway, p53 signaling pathway, longevity regulating pathway). 2) Through differential expression analysis, we obtained 25 cellular senescence-associated DEGs of RIF and verified their specific molecular mechanisms in regulating cellular senescence and immunity through functional enrichment analysis. 3) A total of eight RIF signature genes (LATS1, EHF, DUSP16, ADCK5, PATZ1, DEK, MAP2K1, and ETS2) were screened by machine learning of SVM-RFE, RF, and ANN. In addition, we accessed the BioGPS, Human Protein Atlas, GeneCards, Alliance of Genome Resources, and UniProt databases to analyze the localization and function of these eight RIF signature genes and discovered that they were all expressed in uterine tissue and immunocytes. Through findings 1–3, we preliminarily concluded that the pathological mechanism of RIF was related to cellular senescence and abnormal immune regulation. Moreover, senescence-associated DEGs exhibited the potential as RIF biomarkers, which also suggested the close relationship between cellular senescence and RIF. 4) We then analyzed the endometrial immune landscape during WOI in patients with RIF and fertile controls, including infiltrating immunocytes, immune function, HLA gene sets, and immune checkpoints, as well as evaluated their correlation with RIF signature genes. 5) In addition, we also performed consensus clustering analysis according to the RIF signature genes, and the differences in the endometrial immune landscape between subtypes also suggested the correlation between cellular senescence and immune regulation. Through findings 4–5, we could further conclude that cellular senescence was closely associated with the abnormal endometrial immunoregulation during WOI in RIF.
By integrating the above results, the following key indications deserved to be explored in depth. First, we found that MAP2K1 was one of the significantly downregulated cellular senescence genes in RIF both in the training set and the testing set and was involved in the most signal transduction, such as MAPK signaling pathway, Hippo signaling pathway, etc. Although MAP2K1 has not attracted enough attention in RIF yet, the detection of MAP2K1 has been regarded as a vital biomarker for evaluating cellular senescence in various diseases (61, 62). Therefore, this study suggested that we could focus on MAP2K1, which might bring new discoveries for the diagnosis and treatment of RIF from the aspect of cellular senescence in the future. In addition, among the various pathways regulated by MAP2K1, we take the MAPK signaling pathway as an example. Studies have confirmed that MAP2K1 encodes a dual specificity protein kinase which lies upstream of MAP kinases and stimulated the enzymatic activity of MAP kinases upon a wide variety of extra- and intracellular signals. Moreover, the MAP kinase/ERK cascade has been verified to be inactive in senescent cells, which is capable of significantly reducing many proinflammatory components of the SASP (63, 64). Fernandes et al. also demonstrated that MAP kinase/ERK rewired in senescent cells rendering them phenotypically different (65). Furthermore, studies implied that MAPK influenced cell survival in the endometrium (66). From the aforementioned analysis, we hypothesized that the decreased mRNA profiling of MAP2K1 in RIF tended to restrain the activation of the MAPK signaling pathway, thus procuring cellular senescence and impairing endometrial receptivity. In addition, to intuitively comprehend the involvement of the eight signature genes in the biologic process, the mapping relationship between genes and signaling pathways was visualized in Figure 9. The six signaling pathways communicated by RIF signature genes interacted with each other and formed a complex biological network, which provided abundant clues for follow-up research.
Figure 9 Construction of the “KEGG pathway–RIF signature genes–immune landscape” association network. The network consists of 86 nodes and 160 edges. Brown nodes represent the KEGG pathway, blue nodes represent RIF signature genes, red nodes represent immune cell type, orange nodes represent immune function, purple nodes represent HLA-related genes, and yellow nodes represent immune checkpoint-related genes. The black lines represent the relationships between nodes. The red lines represent facilitation effects between pathways, and the green lines represent inhibition between pathways.
When it comes to immunity and pregnancy, it is universally acknowledged that immune tolerance to the semi-allogeneic fetus is significantly rudimentary for a successful pregnancy. Abnormalities of endometrial immunity are involved in the pathogenesis of RIF (67, 68). At present, research on the immune microenvironment in the endometrium of RIF during WOI is in full swing. However, due to the complex composition of immunocytes and the complicated cross-talk of immune functions, traditional research methods fail to systematically reflect the immune landscape in the endometrium. This time, we utilized ssGSEA to analyze 15 types of infiltrated immunocytes and 13 types of immune functions in RIF and control groups and also compared the expression of HLA molecules and immune checkpoints between groups, both in the training and validation sets, respectively. The above analysis revealed that the immune landscape of RIF was significantly different from fertile controls. Take the Treg lymphocyte as an example. It is a pivotal class of immunosuppressive cells, expresses forkhead box P3 (Foxp3), and secretes the anti-inflammatory cytokines, such as TGF-β1 and IL-10, which can suppress the excessive immune response induced by the fetus, playing a pushing role in maintaining an immune tolerance microenvironment during pregnancy (69, 70). Moreover, higher levels of endometrial Treg cells have been confirmed to be a positive prognostic factor for improving pregnancy outcome in IVF (71). Our results manifested that Treg cells were reduced in patients with RIF when compared with the control group, which was consistent with multiple studies (72–74). However, the reason for the reduction of Treg cells in RIF was still perplexing.
We speculated that the disorder of immunocytes and immune function was related to cellular senescence. Continuing with the example of reduced Treg cells in RIF, our study illuminated that Treg cells were significantly negatively correlated with LATS1 and positively correlated with ETS2. LATS1 is the core kinases of Hippo effector Yes-associated protein 1 (YAP1), which can lead to the phosphorylation and inhibition of YAP1 (75), thus regulating cell proliferation, differentiation, and boosting cellular senescence (76–78). As for ETS2, it is a member of the ETS family of DNA-binding transcription factors. Studies have suggested that ETS2 can bind to the proximal promoter of human telomerase reverse transcriptase (hTERT) and positively regulate hTERT, thus inhibiting cellular senescence (79). Through the above analysis, we conjectured that the increase of LATS1 and the decrease of ETS2 in RIF may reduce the number of local endometrial Treg cells by promoting cellular senescence. Given that the above gene expression is not tissue specific, we cannot determine in which cells senescence occurs. According to published studies, we assumed that there were three possibilities: 1) the increase of senescent decidual cells leads to an excess of IL-6 and TGF-β, the most common SASP factors, which can induce the differentiation of CD4+ T cells to Th17 and reduce the level of Treg cells (80, 81). 2) The clearance mechanism of senescent cells by immunocytes is abnormal. Senescent cells are subject to immune surveillance by multiple components of the immune system (82), especially immunocytes with the function of phagocytosis. Moreover, impaired immune surveillance can lead to senescent cell accumulation (83, 84). Studies have affirmed that immunocytes exert a pivotal role in endometrial fate decisions at implantation (85). One of the principal functions of macrophages is phagocytic clearance of senescent cells (86, 87). However, we observed a decreased level of endometrial macrophages in RIF patients both in training and testing sets, which may account for the accumulation of senescent cells. Moreover, senescent cells attract and activate immune cells and serve as highly immunogenic targets for immune clearance, the microenvironment of which may prompt the Treg/Th17 balance to tilt toward the Th17 cell. 3) The senescence of CD4+ T lymphocytes themselves leads to the reduction of Treg cells. As a cellular counterpart, CD4+ T lymphocytes can also undergo cellular senescence, which could also contribute to an immune response biased toward Th17 lymphocytes from Treg lymphocytes (88). Now, senescent T cells are corroborated to be accumulated in aging, chronic viral infections, and autoimmune disorders, which are also considered as potential targets for disease treatment (89–91). However, the current research on immunosenescence in RIF disease is still blank. To systematically and intuitively grasp the correlation between immune landscape and cellular senescence in RIF, we visualized the association of immune cells, immune function, HLA molecules, and immune checkpoints with the eight cellular senescence-related signature genes of RIF in Figure 9. Taken together, our findings undoubtedly threw innovative light on the mechanism of an abnormal immune microenvironment in RIF, which also contributed to uncovering potential intervention targets for drug design.
Conclusion
In conclusion, our study is the first to explore the involvement of cellular senescence in the pathological mechanism of RIF at the molecular level based on bioinformatics combined with machine learning strategy. This study signified that cellular senescence was a potential pathological mechanism of RIF, and cellular senescence-associated genes were expected to serve as novel diagnostic biomarkers for RIF. Moreover, cellular senescence was involved in the regulation of the endometrial immune microenvironment in RIF during WOI, which provided a groundbreaking direction for the exploration of the pathogenesis of abnormal immunity in RIF. At present, our understanding of cellular senescence is only at the tip of the iceberg, and knowledge about its correlation with the pathological mechanism of RIF is still quite finite. Combining cellular senescence with the regulation of the immune microenvironment landscape to reveal the pathogenic mechanism of RIF is revolutionary for the research of RIF, which makes up for the blank in the field of immune mechanism in RIF to a great extent. This study will also encourage more researchers to carry out senescence-related research in the field of RIF. However, this study also has some limitations, which we must admit. Firstly, the results of this study are based on bioinformatics analysis. Although verified by external testing datasets, its accuracy needs to be further substantiated by experiments. Secondly, immune infiltration analysis in this study uses the most widely used ssGSEA method to quantify the number of immunocytes, but single-cell sequencing is still necessary to obtain the most accurate information. Last but not least, this study only found the correlation between cellular senescence and abnormal immune microenvironment in the endometrium of RIF during WOI, but the cause and effect between the two calls for further exploration. In conclusion, all our findings confirm the involvement of cellular senescence in RIF and its close correlation with the immune characteristics of RIF, which provide new insights into the pathogenesis, diagnosis, and treatment of RIF.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Author contributions
XZ and QZ: conceptualization. XZ and YZ: data curation. YJ and QZ: data analysis. XZ and YZ: writing—original draft. QZ: reviewed the manuscript. All authors read and approved the final article.
Funding
This study was supported by the Research on the Secondary Development of Jianpi Antai Mixture (2021C03080).
Acknowledgments
Special thanks to Yuepeng Jiang from the Zhejiang Chinese Medical University for his contribution to the data analysis.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.952708/full#supplementary-material
Supplementary Figure 1 | Data preprocessing and normalization. The box plots of the GSE26787 and GSE111974 datasets before (A) and after (B) normalization. Three- dimensional PCA cluster plot of the GSE26787 and GSE111974 datasets before (C) and after (D) normalization.Green nodes represent GSE111974 and blue nodes represent GSE26787.
Supplementary Figure 2 | Construction of WGCNA modules. (A) Sample dendrogram of control and RIF group. (B) Scale independence and mean connectivity under diverse soft-thresholding powers. (C) Co-expression modules for each gene under the hierarchical clustering tree are assigned diverse colors. The dynamic tree cutcorresponds to the original modules, and merged dynamic corresponds to the merged modules finally identified. (D)Heatmap of the relationships between co-expression modules and clinical traits. The number indicates the correlation coefficients between co-expression modules and clinical traits, and the number in parentheses indicates the corresponding p-values. (E)GO enrichment analysis of genes in the green module. (F)KEGG enrichment analysis of genes in the green module.
Supplementary Figure 3 | Function enrichment of cellular senescence associated DEGs. (A, B) Sankey-bubble diagram showing the GO and KEGG enrichment analysis of cellular senescence associated DEGs. (C) GSEA results for the activation of GO enrichment in control samples. (D) GSEA results for the activation of GO enrichment in RIF samples. (E) GSEA results for the activation of KEGG enrichment in control samples. (F) GSEA results for the activation of KEGG enrichment in RIF samples.
Supplementary Figure 4 | Machine learning of SVM-RFE and RF algorithms for clarifying the diagnostic value of cellular senescence associated DEGs. (A)Random forest of cellular senescence associated DEGs. X-axis represents the number of trees,Y-axis represents the error of cross validation, green curve represents the error of control group, red curve represents the error of RIF group, and black curve represents the error of all samples. (B)MeanDecreaseGini coefficient diagram: X-axis is MeanDecreaseGini value, which is used to judge the classification of the model; Y-axis is gene name. In this study, genes corresponding to MeanDecreaseGini≥1 were selected. (C)SVM-REFscreening feature genes: X-axis represents the change of gene number, Y-axis represents the error of cross validation.
Supplementary Figure 5 | Immune related infiltration landscape in testing set. (A)The distribution of immunocytes and immune functions between RIF and control samples in testing set. (B)The abundance differences of different HLA between RIF and control samples in testing set. (C)The abundance differences of different immune checkpoint between RIF and control samples in testing set. * represents P<0.05 compared with the control group,**represents P<0.01 compared with the control group, ***represents P<0.001 compared with the control group.
Supplementary Figure 6 | Identifying 2 distinct subtypes in RIF by unsupervised clustering based on 8 RIF signature genes. (A) The consensus clustering cumulative distribution function (CDF) for k=2-9. (B) The relative change in area under the CDF curve for k=2-9. (C) Item tracking plot showing the consensus cluster of items (in column) at each k (in row). (D) Consensus matrix plots depicting consensus values on a white to blue color scale ordered by consensus clustering when k= 2. (E) t-SNE plots confirming the classification accuracy of two distinct subtypes across RIF samples. (F–G) Box diagram and heatmap showed the expression level of 8 RIF signature genes in the two subtypes. * represents P<0.05 compared with the cluster A,**represents P<0.01 compared with the cluster A, ***represents P<0.001 compared with the cluster A.
Abbreviations
AUC, area under the ROC; CDF, cumulative distribution function; DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; IVF, in-vitro fertilization; RF, random forest; RIF, recurrent implantation failure; ROC, receiver operating characteristic curve; ssGSEA, single-sample gene-set enrichment analysis; SVM-RFE, support vector machine-recursive feature elimination; WGCNA, weighted gene co-expression network analysis.
References
1. Parvanov D, Ganeva R, Vidolova N, Stamenov G. Decreased number of p16-positive senescent cells in human endometrium as a marker of miscarriage. J Assist Reprod Genet (2021) 38:2087–95. doi: 10.1007/s10815-021-02182-5
2. 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:14–38. doi: 10.1016/j.rbmo.2013.08.011
3. Pantos K, Grigoriadis S, Maziotis E, Pistola K, Xystra P, Pantou A, et al. The role of interleukins in recurrent implantation failure: A comprehensive review of the literature. Int J Mol Sci (2022) 23:2198. doi: 10.3390/ijms23042198
4. Cimadomo D, Capalbo A, Dovere L, Tacconi L, Soscia D, Giancani A, et al. Leave the past behind: Women's reproductive history shows no association with blastocysts' euploidy and limited association with live birth rates after euploid embryo transfers. Hum Reprod (2021) 36:929–40. doi: 10.1093/humrep/deab014
5. Teh WT, McBain J, Rogers P. What is the contribution of embryo-endometrial asynchrony to implantation failure? J Assist Reprod Genet (2016) 33:1419–30. doi: 10.1007/s10815-016-0773-6
6. Gellersen B, Brosens JJ. Cyclic decidualization of the human endometrium in reproductive health and failure. Endocr Rev (2014) 35:851–905. doi: 10.1210/er.2014-1045
7. Ohara Y, Matsubayashi H, Suzuki Y, Takaya Y, Yamaguchi K, Doshida M, et al. Clinical relevance of a newly developed endometrial receptivity test for patients with recurrent implantation failure in Japan. Reprod Med Biol (2022) 21:e12444. doi: 10.1002/rmb2.12444
8. Al-Lamee H, Ellison A, Drury J, Hill CJ, Drakeley AJ, Hapangama DK, et al. Altered endometrial oestrogen-responsiveness and recurrent reproductive failure. Reprod Fertil (2022) 3:30–8. doi: 10.1530/RAF-21-0093
9. Zhai J, Ma L, Chang Z, Yu T. Increased expression of prokineticin 2 and its receptor in endometrium of recurrent implantation failure patients decreased the expression of MMP9 important for decidualization. Reprod Biol Endocrinol (2022) 20:76. doi: 10.1186/s12958-022-00947-w
10. Hapangama DK, Turner MA, Drury JA, Martin-Ruiz C, Von Zglinicki T, Farquharson RG, et al. Endometrial telomerase shows specific expression patterns in different types of reproductive failure. Reprod BioMed Online (2008) 17:416–24. doi: 10.1016/S1472-6483(10)60227-1
11. Roy AL, Sierra F, Howcroft K, Singer DS, Sharpless N, Hodes RJ, et al. A blueprint for characterizing senescence. Cell (2020) 183:1143–6. doi: 10.1016/j.cell.2020.10.032
12. van Deursen JM. The role of senescent cells in ageing. Nature (2014) 509:439–46. doi: 10.1038/nature13193
13. Peter DR, Aberkane A, Polanski L, Maruyama Y, Baumgarten M, Lucas ES, et al. Deregulation of the endometrial stromal cell secretome precedes embryo implantation failure. Mol Hum Reprod (2017) 23:478–87. doi: 10.1093/molehr/gax023
14. Deryabin P, Griukova A, Nikolsky N, Borodkina A. The link between endometrial stromal cell senescence and decidualization in female fertility: The art of balance. Cell Mol Life Sci (2020) 77:1357–70. doi: 10.1007/s00018-019-03374-0
15. Rawlings TM, Makwana K, Taylor DM, Molè MA, Fishwick KJ, Tryfonos M, et al. Modelling the impact of decidual senescence on embryo implantation in human endometrial assembloids. Elife (2021) 10:e69603. doi: 10.7554/eLife.69603
16. Bert S, Ward EJ, Nadkarni S. Neutrophils in pregnancy: New insights into innate and adaptive immune regulation. Immunology (2021) 164:665–76. doi: 10.1111/imm.13392
17. Köstlin-Gille N, Dietz S, Schwarz J, Spring B, Pauluschke-Fröhlich J, Poets CF, et al. HIF-1α-Deficiency in myeloid cells leads to a disturbed accumulation of myeloid derived suppressor cells (MDSC) during pregnancy and to an increased abortion rate in mice. Front Immunol (2019) 10:161. doi: 10.3389/fimmu.2019.00161
18. 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:548–82. doi: 10.1093/humrep/deab130.425
19. Wang H, Fan Y, Chen W, Lv Z, Wu S, Xuan Y, et al. Loss of CMTM6 promotes DNA damage-induced cellular senescence and antitumor immunity. Oncoimmunology (2022) 11:2011673. doi: 10.1080/2162402X.2021.2011673
20. Manakanatas C, Ghadge SK, Agic A, Sarigol F, Fichtinger P, Fischer I, et al. Endothelial and systemic upregulation of miR-34a-5p fine-tunes senescence in progeria. Aging (Albany NY) (2022) 14:195–224. doi: 10.18632/aging.203820
21. Marquez-Exposito L, Tejedor-Santamaria L, Santos-Sanchez L, Valentijn FA, Cantero-Navarro E, Rayego-Mateos S, et al. Acute kidney injury is aggravated in aged mice by the exacerbation of proinflammatory processes. Front Pharmacol (2021) 12:662020. doi: 10.3389/fphar.2021.662020
22. Abbott TR, Dhamdhere G, Liu Y, Lin X, Goudy L, Zeng L, et al. Development of CRISPR as an antiviral strategy to combat SARS-CoV-2 and influenza. Cell (2020) 181:865–76. doi: 10.1016/j.cell.2020.04.020
23. Sun P, Wu Y, Yin C, Jiang H, Xu Y, Sun H. Molecular subtyping of cancer based on distinguishing Co-expression modules and machine learning. Front Genet (2022) 13:866005. doi: 10.3389/fgene.2022.866005
24. Šín P, Hokynková A, Marie N, Andrea P, Krč R, Podroužek J. Machine learning-based pressure ulcer prediction in modular critical care data. Diagnostics (Basel) (2022) 12:850. doi: 10.3390/diagnostics12040850
25. Edgar R, Domrachev M, Lash AE. Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res (2002) 30:207–10. doi: 10.1093/nar/30.1.207
26. Avelar RA, Ortega JG, Tacutu R, Tyler EJ, Bennett D, Binetti P, et al. A multidimensional systems biology analysis of cellular senescence in aging and disease. Genome Biol (2020) 21:91. doi: 10.1186/s13059-020-01990-9
27. Langfelder P, Horvath S. WGCNA: an r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559
28. Li G, Sun J, Zhang J, Lv Y, Liu D, Zhu X, et al. Identification of inflammation-related biomarkers in diabetes of the exocrine pancreas with the use of weighted gene Co-expression network analysis. Front Endocrinol (Lausanne) (2022) 13:839865. doi: 10.3389/fendo.2022.839865
29. Ito K, Murphy D. Application of ggplot2 to pharmacometric graphics. CPT Pharmacometrics Syst Pharmacol (2013) 2:e79. doi: 10.1038/psp.2013.56
30. Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res (2019) 47:W199–205. doi: 10.1093/nar/gkz401
31. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004
32. Huang ML, Hung YH, Lee WM, Li RK, Jiang BR. SVM-RFE based feature selection and taguchi parameters optimization for multiclass SVM classifier. ScientificWorldJournal (2014) 2014:795624. doi: 10.1155/2014/795624
33. Sanz H, Valim C, Vegas E, Oller JM, Reverter F. SVM-RFE: selection and visualization of the most relevant features through non-linear kernels. BMC Bioinf (2018) 19:432. doi: 10.1186/s12859-018-2451-4
34. Mundra PA, Rajapakse JC. SVM-RFE with MRMR filter for gene selection. IEEE Trans Nanobioscience (2010) 9:31–7. doi: 10.1109/TNB.2009.2035284
35. Biau G, Scornet E. A random forest guided tour. TEST-SPAIN (2015) 25:1–31. doi: 10.1007/s11749-016-0481-7
36. Lee KS, Park H. Machine learning on thyroid disease: a review. Front Biosci (Landmark Ed) (2022) 27:101. doi: 10.31083/j.fbl2703101
37. Tang J, Henderson A, Gardner P. Exploring AdaBoost and random forests machine learning approaches for infrared pathology on unbalanced data sets. Analyst (2021) 146:5880–91. doi: 10.1039/D0AN02155E
38. Shi Y, Ying X, Yang J. Deep unsupervised domain adaptation with time series sensor data: A survey. Sensors (Basel) (2022) 22:5507. doi: 10.3390/s22155507
39. Wu C, Jin X, Tsueng G, Afrasiabi C, Su AI. BioGPS: building your own mash-up of gene annotations and expression profiles. Nucleic Acids Res (2016) 44:D313–6. doi: 10.1093/nar/gkv1104
40. Sjöstedt E, Zhong W, Fagerberg L, Karlsson M, Mitsios N, Adori C, et al. An atlas of the protein-coding genes in the human, pig, and mouse brain. Science (2020) 367:eaay5947. doi: 10.1126/science.aay5947
41. Safran M, Rosen N, Twik M, BarShir R, Stein TI, Dahary D, et al. The GeneCards suite. In: Abugessaisa I, Kasukawa T, editors. Practical guide to life science databases. Singapore: Springer Nature Singapore (2021). p. 27–56.
42. Alliance of Genome Resources Consortium. Harmonizing model organism data in the alliance of genome resources. Genetics (2022) 220:iyac022. doi: 10.1093/genetics/iyac022
43. UniProt Consortium. UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Res (2021) 49:D480–9. doi: 10.1093/nar/gkaa1100
44. Le T, Aronow RA, Kirshtein A, Shahriyari L. A review of digital cytometry methods: estimating the relative abundance of cell types in a bulk of cells. Brief Bioinform (2021) 22:bbaa219. doi: 10.1093/bib/bbaa219
45. Peng J, Han L, Shang X. A novel method for predicting cell abundance based on single-cell RNA-seq data. BMC Bioinf (2021) 22:281. doi: 10.1186/s12859-021-04187-4
46. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170
47. Tegge AN, Sharp N, Murali TM. Xtalk: a path-based approach for identifying crosstalk between signaling pathways. Bioinformatics (2016) 32:242–51. doi: 10.1093/bioinformatics/btv549
48. Jia Y, Dong YJ, Sha YL, Cai SC, Diao LH, Qiu Z, et al. Effectiveness comparison between endometrial receptivity array, immune profiling and the combination in treating patients with multiple implantation failure. Am J Reprod Immunol (2022) 87:e13513. doi: 10.1111/aji.13513
49. Liu Z, Liu X, Wang M, Zhao H, He S, Lai S, et al. The clinical efficacy of personalized embryo transfer guided by the endometrial receptivity array/analysis on IVF/ICSI outcomes: A systematic review and meta-analysis. Front Physiol (2022) 13:841437. doi: 10.3389/fphys.2022.841437
50. Kusama K, Yamauchi N, Yoshida K, Azumi M, Yoshie M, Tamura K. Senolytic treatment modulates decidualization in human endometrial stromal cells. Biochem Biophys Res Commun (2021) 571:174–80. doi: 10.1016/j.bbrc.2021.07.075
51. Cheloufi M, Kazhalawi A, Pinton A, Rahmati M, Chevrier L, Prat-Ellenberg L, et al. The endometrial immune profiling may positively affect the management of recurrent pregnancy loss. Front Immunol (2021) 12:656701. doi: 10.3389/fimmu.2021.656701
52. Deryabin PI, Borodkina AV. Stromal cell senescence contributes to impaired endometrial decidualization and defective interaction with trophoblast cells. Hum Reprod (2022) 37:1505–24. doi: 10.1093/humrep/deac112
53. Chen P, Yang M, Wang Y, Guo Y, Liu Y, Fang C, et al. Aging endometrium in young women: Molecular classification of endometrial aging-based markers in women younger than 35 years with recurrent implantation failure. J Assist Reprod Genet (2022). doi: 10.1007/s10815-022-02578-x
54. Lunin SM, Novoselova EG, Glushkova OV, Parfenyuk SB, Novoselova TV, Khrenov MO. Cell senescence and central regulators of immune response. Int J Mol Sci (2022) 23:4109. doi: 10.3390/ijms23084109
55. Shive C, Pandiyan P. Inflammation, immune senescence, and dysregulated immune regulation in the elderly. Front Aging (2022) 3:840827. doi: 10.3389/fragi.2022.840827
56. Childs BG, Baker DJ, Wijshake T, Conover CA, Campisi J, van Deursen JM. Senescent intimal foam cells are deleterious at all stages of atherosclerosis. Science (2016) 354:472–7. doi: 10.1126/science.aaf6659
57. Mikawa R, Suzuki Y, Baskoro H, Kanayama K, Sugimoto K, Sato T, et al. Elimination of p19(ARF) -expressing cells protects against pulmonary emphysema in mice. Aging Cell (2018) 17:e12827. doi: 10.1111/acel.12827
58. Wang Z. Regulation of cell cycle progression by growth factor-induced cell signaling. Cells-Basel (2021) 10:3327. doi: 10.3390/cells10123327
59. Childs BG, Baker DJ, Kirkland JL, Campisi J, van Deursen JM. Senescence and apoptosis: Dueling or complementary cell fates? EMBO Rep (2014) 15:1139–53. doi: 10.15252/embr.201439245
60. Menon R, Behnia F, Polettini J, Saade GR, Campisi J, Velarde M. Placental membrane aging and HMGB1 signaling associated with human parturition. Aging (Albany NY) (2016) 8:216–30. doi: 10.18632/aging.100891
61. Liu JF, Wu Y, Yang YH, Wu SF, Liu S, Xu P, et al. Phosphoproteome profiling of mouse liver during normal aging. Proteome Sci (2022) 20:12. doi: 10.1186/s12953-022-00194-2
62. Wang Y, Chen J, Chen X, Jiang F, Sun Y, Pan Y, et al. MiR-34a suppresses HNSCC growth through modulating cell cycle arrest and senescence. Neoplasma (2017) 64:543–53. doi: 10.4149/neo_2017_408
63. Sun Y, Liu WZ, Liu T, Feng X, Yang N, Zhou HF. Signaling pathway of MAPK/ERK in cell proliferation, differentiation, migration, senescence and apoptosis. J Recept Signal Transduct Res (2015) 35:600–4. doi: 10.3109/10799893.2015.1030412
64. Freund A, Patil CK, Campisi J. p38MAPK is a novel DNA damage response-independent regulator of the senescence-associated secretory phenotype. EMBO J (2011) 30:1536–48. doi: 10.1038/emboj.2011.69
65. Fernandes SE, Saini DK. The ERK-p38MAPK-STAT3 signalling axis regulates iNOS expression and salmonella infection in senescent cells. Front Cell Infect Microbiol (2021) 11:744013. doi: 10.3389/fcimb.2021.744013
66. Bora G, Yaba A. The role of mitogen-activated protein kinase signaling pathway in endometriosis. J Obstet Gynaecol Res (2021) 47:1610–23. doi: 10.1111/jog.14710
67. Kuroda K, Matsumura Y, Ikemoto Y, Segawa T, Hashimoto T, Fukuda J, et al. Analysis of the risk factors and treatment for repeated implantation failure: OPtimization of thyroid function, IMmunity, and uterine milieu (OPTIMUM) treatment strategy. Am J Reprod Immunol (2021) 85:e13376. doi: 10.1111/aji.13376
68. Yamaguchi K. Tacrolimus treatment for infertility related to maternal-fetal immune interactions. Am J Reprod Immunol (2019) 81:e13097. doi: 10.1111/aji.13097
69. Krop J, Heidt S, Claas F, Eikmans M. Regulatory T cells in pregnancy: It is not all about FoxP3. Front Immunol (2020) 11:1182. doi: 10.3389/fimmu.2020.01182
70. Tsuda S, Nakashima A, Morita K, Shima T, Yoneda S, Kishi H, et al. The role of decidual regulatory T cells in the induction and maintenance of fetal antigen-specific tolerance: Imbalance between regulatory and cytotoxic T cells in pregnancy complications. Hum Immunol (2021) 82:346–52. doi: 10.1016/j.humimm.2021.01.019
71. Sauerbrun-Cutler MT, Huber WJ, Krueger PM, Sung CJ, Has P, Sharma S. Do endometrial natural killer and regulatory T cells differ in infertile and clinical pregnancy patients? an analysis in patients undergoing frozen embryo transfer cycles. Am J Reprod Immunol (2021) 85:e13393. doi: 10.1111/aji.13393
72. Fan L, Sha M, Li W, Kang Q, Wu J, Chen S, et al. Intrauterine administration of peripheral blood mononuclear cells (PBMCs) improves embryo implantation in mice by regulating local Treg/Th17 cell balance. J Reprod Dev (2021) 67:359–68. doi: 10.1262/jrd.2021-006
73. Pourmoghadam Z, Soltani-Zangbar MS, Sheikhansari G, Azizi R, Eghbal-Fard S, Mohammadi H, et al. Intrauterine administration of autologous hCG- activated peripheral blood mononuclear cells improves pregnancy outcomes in patients with recurrent implantation failure; a double-blind, randomized control trial study. J Reprod Immunol (2020) 142:103182. doi: 10.1016/j.jri.2020.103182
74. Sadeghpour S, Ghasemnejad BM, Nazarian H, Ghasemnejad T, Nematollahi MH, Abroon S, et al. Effects of treatment with hydroxychloroquine on the modulation of Th17/Treg ratio and pregnancy outcomes in women with recurrent implantation failure: clinical trial. Immunopharmacol Immunotoxicol (2020) 42:632–42. doi: 10.1080/08923973.2020.1835951
75. Avruch J, Zhou D, Fitamant J, Bardeesy N, Mou F, Barrufet LR. Protein kinases of the hippo pathway: regulation and substrates. Semin Cell Dev Biol (2012) 23:770–84. doi: 10.1016/j.semcdb.2012.07.002
76. Xu X, Shen X, Wang J, Feng W, Wang M, Miao X, et al. YAP prevents premature senescence of astrocytes and cognitive decline of alzheimer's disease through regulating CDK6 signaling. Aging Cell (2021) 20:e13465. doi: 10.1111/acel.13465
77. Rosado-Olivieri EA, Anderson K, Kenty JH, Melton DA. YAP inhibition enhances the differentiation of functional stem cell-derived insulin-producing β cells. Nat Commun (2019) 10:1464. doi: 10.1038/s41467-019-09404-6
78. Jin H, Lian N, Zhang F, Bian M, Chen X, Zhang C, et al. Inhibition of YAP signaling contributes to senescence of hepatic stellate cells induced by tetramethylpyrazine. Eur J Pharm Sci (2017) 96:323–33. doi: 10.1016/j.ejps.2016.10.002
79. Flavin P, Redmond A, McBryan J, Cocchiglia S, Tibbitts P, Fahy-Browne P, et al. RuvBl2 cooperates with Ets2 to transcriptionally regulate hTERT in colon cancer. FEBS Lett (2011) 585:2537–44. doi: 10.1016/j.febslet.2011.07.005
80. Atreya R, Mudter J, Finotto S, Müllberg J, Jostock T, Wirtz S, et al. Blockade of interleukin 6 trans signaling suppresses T-cell resistance against apoptosis in chronic intestinal inflammation: Evidence in crohn disease and experimental colitis in vivo. Nat Med (2000) 6:583–8. doi: 10.1038/75068
81. Li X, Cannon AR, Hammer AM, Morris NL, Choudhry MA. IL-23 restoration of Th17 effector function is independent of IL-6 and TGF-β in a mouse model of alcohol and burn injury. J Leukoc Biol (2017) 102:915–23. doi: 10.1189/jlb.3A1216-527R
82. Xue W, Zender L, Miething C, Dickins RA, Hernando E, Krizhanovsky V, et al. Senescence and tumour clearance is triggered by p53 restoration in murine liver carcinomas. Nature (2007) 445:656–60. doi: 10.1038/nature05529
83. Ovadya Y, Landsberger T, Leins H, Vadai E, Gal H, Biran A, et al. Impaired immune surveillance accelerates accumulation of senescent cells and aging. Nat Commun (2018) 9:5435. doi: 10.1038/s41467-018-07825-3
84. Kong CS, Ordoñez AA, Turner S, Tremaine T, Muter J, Lucas ES, et al. Embryo biosensing by uterine natural killer cells determines endometrial fate decisions at implantation. FASEB J (2021) 35:e21336. doi: 10.1096/fj.202002217R
85. Brighton PJ, Maruyama Y, Fishwick K, Vrljicak P, Tewary S, Fujihara R, et al. Clearance of senescent decidual cells by uterine natural killer cells in cycling human endometrium. Elife (2017) 6:e31274. doi: 10.7554/eLife.31274
86. Recalcati S, Cairo G. Macrophages and iron: A special relationship. Biomedicines (2021) 9:1585. doi: 10.3390/biomedicines9111585
87. Behmoaras J, Gil J. Similarities and interplay between senescent cells and macrophages. J Cell Biol (2021) 220:e202010162. doi: 10.1083/jcb.202010162
88. González-Osuna L, Sierra-Cristancho A, Cafferata EA, Melgar-Rodríguez S, Rojas C, Carvajal P, et al. Senescent CD4(+)CD28(-) T lymphocytes as a potential driver of Th17/Treg imbalance and alveolar bone resorption during periodontitis. Int J Mol Sci (2022) 23:2543. doi: 10.3390/ijms23052543
89. Salminen A. Immunosuppressive network promotes immunosenescence associated with aging and chronic inflammatory conditions. J Mol Med (Berl) (2021) 99:1553–69. doi: 10.1007/s00109-021-02123-w
90. Zhang J, He T, Xue L, Guo H. Senescent T cells: a potential biomarker and target for cancer therapy. Ebiomedicine (2021) 68:103409. doi: 10.1016/j.ebiom.2021.103409
Keywords: recurrent implantation failure, immune landscape, cellular senescence, bioinformatics, machine learning
Citation: Zhao X, Zhao Y, Jiang Y and Zhang Q (2022) Deciphering the endometrial immune landscape of RIF during the window of implantation from cellular senescence by integrated bioinformatics analysis and machine learning. Front. Immunol. 13:952708. doi: 10.3389/fimmu.2022.952708
Received: 25 May 2022; Accepted: 17 August 2022;
Published: 05 September 2022.
Edited by:
Maria Foti, University of Milano-Bicocca, ItalyReviewed by:
Qiang Fu, Binzhou Medical University, ChinaPhilippa Saunders, University of Edinburgh, United Kingdom
Copyright © 2022 Zhao, Zhao, Jiang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qin Zhang, WmhhcWluMDFAMTYzLmNvbQ==