Artificial Neural Network Model Using Immune-infiltration Modules for Endometrial Receptivity Assessment of Implantation Failure

doi:10.21203/rs.3.rs-860927/v1

Download PDF

Research

Artificial Neural Network Model Using Immune-infiltration Modules for Endometrial Receptivity Assessment of Implantation Failure

https://doi.org/10.21203/rs.3.rs-860927/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Objectives

This study was anchored on the state of local immune-infiltration in the endometrium, which acts as critical factors affecting embryonic implantation, and aimed at establishing novel approaches to assess endometrial receptivity for patients with IVF failure.

Methods

Immune-infiltration levels in the GSE58144 dataset (n=115) from GEO were analyzed by digital deconvolution and validated by immunofluorescence (n=30), illustrating that dysregulation of the ratio of Mf1 to Mf2 is an important factor contributing to implantation failure. Then, modules most associated with M1/M2 macrophages (Mfs) and their hub genes were then selected by weighted gene co-expression network and univariate analyses, then validated by GSE5099 macrophage dataset, qPCR analysis (n=16), and western blot. It revealed that closely related gene modules dominated three biological processes in macrophages: antigen presentation, interleukin−1−mediated signalling pathway, and phagosome acidification, respectively. Their hub genes were significantly altered in patients and related with ribosomal, lysosome, and proteasomal pathways. Finally, the artificial neural network (ANN) and nomogram models were established from hub genes, of which efficacy was compared and validated in the GSE165004 dataset (n=72). Models established by the selected hub genes exhibited excellent predictive values in both datasets, and ANN performed best with an accuracy of 98.3% and an AUC of 0.975 (95% CI 0.945-1).

Conclusions

Macrophages, proven to be essential for endometrial receptivity, were regulated by gene modules dominating antigen presentation, interleukin−1−mediated signalling pathway, and phagosome acidification. Selected hub genes can effectively assess endometrial dysfunction receptivity for IVF outcomes by the ANN approach.

Sexual & Reproductive Medicine

Cancer Biology

Immune infiltration

Endometrial receptivity

Implantation failure

Deep machine learning

Artificial neural network

Recurrent implantation failure (RIF) is one of the most frustrating and difficult areas in reproductive medicine because the etiology is often unknown and there are few evidence-based diagnostic and treatment strategies. Defective endometrial receptivity is currently becoming a critical theory in the study of this disease. Endometrial receptivity, is a complex process that enables embryonic attachment, invasion, and development. For healthy females, during the secretory phase, the window of implantation (WOI) lasts from 3 to 6 days. In certain inflammatory or anatomical cases, this window can be narrowed or displaced to inhibit normal implantation, resulting in infertility or loss of pregnancy (1). Therefore, a receptive endometrium is a prerequisite for successful embryonic implantation. Defective endometrial receptivity is often associated with RIF and unexplained infertility (UI). However, a lack of understanding of defective endometrial receptivity has led to poor diagnosis and treatment. The current diagnostic tests for endometrial receptivity anomalies include sonography (endometrial thickness, blood flow, and morphology), histopathology (integrin αvβ3, leukemia inhibitory factor, vascular endothelial growth factor, and uterine natural killer cells), as well as electron microscopy of endometrial cell morphology. However, these methods have a limited clinical guidance (2). An emerging assay, Endometrial Receptivity Array (ERA), guides clinical practice by assessing endometrial receptivity status through microarray analysis of 238 genes (3). However, the clinical significance of the test has not been fully established (4). The high number of genes that are required for analysis and strict requirements for sample storage conditions lead to high costs and difficulties in clinical practice.

Immune status of the endometrium is closely associated with normal reproductive functions(5). Trophoblast cells of the implanted blastocyst invade the endometrium of the maternal uterus by forming the placenta. Endometrial mesenchymal cells undergo a decidual response (called decidualization) and establish an environment that is conducive for trophoblast invasion. Trophoblast implantation and placenta formation require appropriate maternal immune tolerance to the hemi allogeneic foetus. After embryonic implantation and decidual development, the number of leukocytes in the uterus change significantly (6). Natural killer (NK) cells and macrophages predominate during the early gestation period. Macrophages, which account for approximately 20–30% of total infiltrating leukocytes, play an essential role in fetal tolerance, trophoblast invasion, and tissue as well as vascular remodelling (7). The immune environment of the endometrium is in a dynamic balance due to hormonal regulation and stress responses to clinical operation and dysbacteriosis. Therefore, endometrial receptivity assessment based on immune-related modules can predict implantation outcomes and predict pregnancy timing after a uterine cavity procedure (induced abortion, adhesiolysis, and polypectomy) or other therapies (antibiotics, hormones etc). However, few studies have evaluated its value in assessing endometrial receptivity.

Due to genetic heterogeneity, epistatic interactions, and environmental factors, identification of a single gene or pathway underlying the complex traits is difficult (8). Integration of gene expression is a key method for solving this problem. Network methods have been used to identify and characterize various biological interactions and have helped in predicting gene functions (9). In this study, we used weighted gene co-expression network analysis (WGCNA) to establish gene modules and evaluated their association with integrated clinical traits(10). The predictors selected in this approach have representative biological structures and functions. Previous evaluation models were established by regression equations. However, measurement errors across platforms and patients' individualization often make it difficult to guarantee the consistency of results, which in turn affects repeatability of the method. The artificial intelligence approach can be used to guarantee the consistency of results across platforms and populations by training the programs with specific patterns. The artificial neural network (ANN) is an artificial intelligence (machine learning) method that works similarly to the human brain(11). Its incorporated feature variables, also known as predictors, input variables and covariates, are the input signals that inform pattern recognition. Each characteristic variable is weighted according to its clinical significance. The task is accomplished by dendrites in the biological nervous system. An activation function sums the weighted signals (12). We performed WGCNA to identify immunological factors that were most associated with defective endometrial receptivity and used them to develop a fertility prognostic model through ANN. The results were validated in different platforms and populations to obtain a more practical and reliable approach for endometrial receptivity assessment.

2.1 Datasets and patient selection

All datasets were selected from chip microarrays. The GSE58144 dataset, containing 43 RIF patients and 72 controls, was used for immune infiltration analysis, hub gene selection, and machine deep learning prognostic model establishment. The GSE5099 dataset is a three replicate measure dataset of gene expression matrices of Mφ1 and Mφ2 that was used to validate the association between hub genes and macrophage polarization. The GSE165004 dataset includes 48 implantation failure (IF) patients and 24 controls, and was used to validate hub gene expression differences and machine deep learning prognostic model. Details of the microarrays are available in supplementary Table 1.

Validation of Mφ1/Mφ2, mRNA and protein levels were performed using clinical samples. Ethical approval for this study was obtained from the Research Ethics Committee of the Beijing Obstetrics and Gynecology Hospital. Experiments were performed (under protocol number 2017-KY-082-02) based on the Helsinki Declaration of 1975 (revised in 2013). Patients eligible for hysteroscopy were required to sign an informed consent before surgery. Samples from endometrial biopsies were acquired at Beijing Obstetrics and Gynecology Hospital. The inclusion criteria for patients were: i. Those aged less than 40 years; ii. Those whose sex hormone levels, including follicle-stimulating hormone, luteinizing hormone, testosterone, estradiol and prolactin, were within normal ranges and iii. Those without endometriosis, fibroids, active or a history of pelvic inflammatory disease or other medical comorbidities (hyperprolactinemia, thyroid disease etc). Study participants in the control group (n = 14) had at least one live birth while others (n = 16) were examined for unexplained infertility or implantation failure. Basic demographic characteristics for each group are presented in supplementary Table 2. Specific flow of the study design is shown in Fig. 1.

2.2 Processing of primary datasets

Pre-processing and normalization of microarray datasets based on raw data of the Affymetrix platform (GSE5099) were performed using the affy package (in R, version 3.6.2) (13) with the following methods: i. Robust multi-array average (RMA, for background correction) (14); ii. Quantile (for normalization) (15); iii. pmonly (perfect match correction) (16) and iv. Median polish (as a summary method) (17). For the microarray datasets that were based on the Agilent platform (GSE58144 and GSE165004), Biobase and limma packages (in R, version 3.6.2) were used for pre-processing and normalization after the data had been converted to log (base 2). The RMA and normalizeBetweenArrays methods were used for background correction and normalization, respectively(18). Annotation files for different microarray platforms were downloaded from the NCBI GEO database (19).

2.3 Digital deconvolution of bulk tissues

Cell-type deconvolution was performed using CIBERSORTx (http://cibersortx.stanford.edu), which is an analytical tool developed by Newman et al. (20) to impute gene expression profiles and provide estimations of the abundances of immune cell infiltration levels in mixed cell populations, using gene expression data. We used the LM22 gene signature matrix for 22 immune cell types. CIBERSORTx was run with batch correction and 100 permutations. Barplot and vioplot were established using the plot function (in R, version 3.6.2).

2.4 Immunofluorescence

From each sample, 5 µm sections were prepared and dewaxed in xylene, dehydrated using graded alcohol and rinsed in distilled water. For antigen retrieval, sections were boiled in citric saline (10 mmol/L, pH 6.0) for half an h. Then, samples were treated with 3% hydrogen peroxide solution for 25 min to block the activity of endogenous peroxidase, blocked using 3% bovine serum albumin (BSA, Servicebio, Wuhan, China) for 30 min at room temperature, after which they were incubated at 37°C for 1 h with primary antibodies, which included mouse anti-CD68 (ab201973, Abcam, dilution 1:200), rabbit anti-CD86 (13395-1-AP, Proteintech, dilution 1:200) for Mφ1 as well as mouse anti-CD68, rabbit anti-CD86 (13395-1-AP, Proteintech, dilution 1:200) for Mφ2. Next, sections were rinsed 3 times in phosphate-buffered saline (PBS) and stained using anti-rabbit-Alexa Fluor® 488 (ab150073, abcam) and anti-mouse-Alexa Fluor® 594 (ab150064, abcam) for 1 h at room temperature (both Invitrogen). Slides were mounted in the SlowfadeGold reagent containing DAPI (Thermofisher, Landsmeer, The Netherlands) and examined using a NIKON CORPORA microscope (Nikon, Japan). The percentage of Mφ1 or Mφ2 was obtained by counting fluorescence positive cells and dividing them by DAPI signal points using Image J (Version 1.50b), respectively.

2.5 Weighted gene co-expression network

The network module was established using the WGCNA package (in the R environment, version 3.6.2) in GSE58144 dataset (21). To minimize noise in the gene expression dataset, data was filtered as follows. Pearson correlation analysis was used to rank all genes according to their association with Mφ1/Mφ2. Correlations between genes and Mφ1/Mφ2 were established at a cut-off of p ≤ 0.05. To reduce the computational burden and enhance signals in our data, we used 2,185 of the 5,531 genes with the greatest variability ranked by the variance in our initial network construction with a cut-off value of var > 0.05 (21). By definition, module genes are highly connected (i.e., module genes tend to have relatively high connectivity). Therefore, for module detection, restricting analysis to the most connected genes should not lead to a major loss of information. Then, we performed cluster analysis of 2,185 genes in these 115 patients. The theory of network construction algorithm has been previously described (22). Briefly, for co-expression module identification, Pearson correlation matrices were first generated (average linkage method) for all pairwise genes. An adjacency matrix was then constructed using a "soft" power adjacency function, \({a}_{ij}={\left|cor\right({x}_{i} {y}_{j}\left)\right|}^{\beta }\). Based on scale-free topology criteria (R² = 0.85), we selected a power of β = 10. In WGCNA, a soft threshold parameter, beta, of the power function was used to ensure that the co-expression network (adjacency matrix) best approximates scale-free topology. This adjacency matrix was then transformed into a topological overlap matrix to measure relative gene interconnectedness and proximity. Finally, gene co-expression modules corresponded to branches of the resulting hierarchical clustering tree (dendrogram). To ensure that genes in the analyzed network exhibited sufficient correlation, we set the weight threshold of the co-expression network to 0.03. Uniform Manifold Approximation and Projection (UMAP) analysis and visualization of the modules was performed using the umap package (in R, version 3.6.2). Visual network diagram was constructed using Cytoscape (version 3.4.0).

2.6 Enrichment analysis of functional categories

The STRING v11.5 online tool (https://string-db.org/) was used for functional enrichment analysis of the gene module that was most associated with endometriosis, identified after WGCNA analysis. In the enrichment analysis, Gene Ontology (GO) terms (including Biological Process BP, Cellular Component CC, and Molecular Function MF) as well as the Kyoto Encyclopedia of Genes and Genome (KEGG) were used to evaluate functional categories and pathways for genes involved in the module. Correlations with Mφ1/Mφ2 for all genes was calculated using the cor function in R environment, version 3.6.2. Then, the GOplot package (in the R environment, version 3.6.2) was used to visualize GO enrichment analysis results and their regulatory conditions. The formula for calculating Z-score is \(z-score=\frac{(up-down)}{\sqrt{count}}\),which is a value that indicates whether GO terms are more likely to be decreased (negative value) or increased (positive value). Gene-set enrichment analysis (GSEA) was performed on the GSE5099 dataset using the “ClusterProfiler” (24) package in R. The Broad Molecular Signature Database (MSigDB v7.0) dataset in the Kyoto Encyclopedia of Genes and Genomes (KEGG) (c2.cp.kegg.v7.0.symbols) was used. This database summarizes and presents specifically well-defined biological states and pathway processes. For statistical significance estimation, the GSEA program was run with 1,000 permutations, and correlations between selected genes and other genes were used to rank all genes.

2.7 Quantitative Real-Time PCR Analysis (qRT-PCR)

Total RNA was extracted from each sample with RNAiso Plus (Takara Bio Inc., Shiga, Japan) and quantified with a NanoDrop™ One Spectrophotometer (Thermo Fisher Scientific Inc., Massachusetts, USA). The First-Strand cDNA Synthesis SuperMix Kit (AT301-3, EasyScript, China) was used to synthesize cDNA from 1 µg of total RNA per sample. The primers used in this study were designed by Sangon Biotech Co., Ltd. Shanghai, China. The sequences are presented in supplementary Table 3. PCRs were performed on an LightCycler 480 PCR System (Roche, Germany) with the protocol for the SYBR Premix Ex TaqTM II (RR820A, Takara). The reaction began at 95°C for 30 seconds for initial denaturation, followed by 35 cycles of 5 seconds at 95°C and 34 seconds at 60°C(23). The measurements were repeated three times, and the relative quantification was performed by the comparative CT (2^-ΔΔCT) method.

2.8 Western blot

Western blot analysis was performed as previously described. After sample preparation, equal amounts of a sample protein (50 µg) were loaded onto an SDS-PAGE gel. After electrophoresis, the samples were transferred onto a nitrocellulose membrane. Then, the membrane was blocked for 2 h at room temperature and incubated overnight at 4°C with the following primary antibodies: anti-MARF1(KIAA0430) (1:500, proteintech, China), anti-dUTPase (DUT) (1:10000, Abcam, USA), anti-RPS9 (1:500, proteintech, China), anti-GADPH (1:1000, LABLEAD, China). Appropriate secondary antibodies (1:5000, ABclonal, China) were selected to incubate with the membrane for 1 h at room temperature. Then, blot bands were visualized with an ECL reagent (Santa Cruz, USA). Non-saturated bands were selected to perform densitometry quantification using Image J software (Version 1.50b).

2.9 prognostic model establishment and statistics

The ANN prognostic model was implemented using the nnet package in R environment, version 3.6.2 (12). To determine the number of units in the hidden layer, the GSE58144 dataset was used: \(Average Accuracy=\sum _{i=1}^{l}\frac{{TP}_{i}+{TN}_{i}}{{TP}_{i}+{FN}_{i}+{FP}_{i}+{TN}_{i}}/l\). Moreover, to ensure maximum optimization of the prognostic model, NNET was run with 100 permutations in the GSE165004 dataset for verification.

For prognostic model efficacy testing, Z tests were used to determine the significance of the area under the receiver operating characteristic (ROC) curve (AUC) using the pROC package in R, version 3.6.2. Sensitivity, specificity, Youden index (YI), positive predictive value (PPV), and negative predictive value (NPV) were used to assess prognostic values of hub genes and the machine learning model.

Comparisons of cell proportions and single gene expression levels between the two groups were performed by Wilcox test. Normally distributed continuous variables were analysed by the Students t-test and paired t-test was used for paired analysis (in R, version 3.6.2).

3.1 Immune infiltration levels on endometrial receptivity

Immune infiltration levels of 22 immune cells in endometrial mixed tissue samples from the GSE58144 dataset were determined using the CIBERSORTx platform. These findings are shown in Fig. 2A. In different groups, there were no significant differences in the percentage of immune cells, except for eosinophils, which were present in deficient levels in the tissues (Fig. 2B). Considering that macrophage polarization may be an important factor in endometrial receptivity, we determined the ratio between Mφ1 and Mφ2 and found significant differences in Mφ1/Mφ2 between IF and normal groups (p = 0.019) (Fig. 2C). We further obtained endometrial tissues from 30 patients within the mid-secretion phase for immunofluorescence detection, and found that the balance between Mφ1 and Mφ2 in the IF group was significantly altered compared to the controls (p = 0.043), as shown in Fig. 2D&E.

3.2 Macrophage polarization-related gene module functions

As shown in Fig. 3A, genes were clustered into different groups, referred to as modules. The GSE58144 gene set had 7 different gene modules with a high topological overlap. To distinguish the modules, we allocated a colour to each module (including brown, black, Green, turquoise, red, blue, and yellow).

Then, we evaluated the pathological correlation for each module by examining the overall correlation of module genes with clinical traits of immune infiltration. The measure of gene significance was defined by the absolute value of the correlation between clinical factors and gene expression levels. Average genetic significance of a particular module is considered module significance (MS). As shown in Fig. 3B, three modules (blue, turquoise, and green) exhibited excellent correlations with Mφ1/Mφ2, that is, -(0.42 (p = 3×10^− 6), -0.35 (p = 10^− 4), and 0.32 (p = 5×10^− 4), respectively. Moreover, there were significant correlations between module membership of the genes within the module and gene significance between these genes and Mφ1/Mφ2 (P_blue=1.6×10^− 39, p _turquoise =6.6×10^− 15, and p _green=0.0093), as shown in Fig. 3C. In Fig. 3D, the integrated molecular profiles of aforementioned three modules are visualized using two-dimensional maps generated by the dimension reduction technique UMAP. In the resulting plot, the samples in the dataset can be grossly divided into two subgroups by the disease, revealing that implantation failure patients have unique profile based on Mφ1/Mφ2 related modules.

Then, we performed enrichment analysis of the three modules. We found that the main biological processes of the three modules were the three aspects of macrophage functions (antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP − dependent, FDR = 4.5×10^− 10; phagosome acidification, FDR = 6.2×10^− 4; interleukin-1-mediated signaling pathway, FDR = 1.13×10^− 13), while the corresponding GO-CC and GO-MF terms of the three were also different. KEGG enrichment analysis revealed that the three modules were enriched in ribosomal, lysosome, and proteasomal pathways, respectively (Fig. 4A).

3.3 Selection and verification of hub genes associated with Mφ1/Mφ2

Due to differences in biological roles of the three gene modules and high-levels of consistency in expression among genes within the modules, we selected three hub genes in each of the modules. First, we selected genes that were most associated with modules and Mφ1/Mφ2 based on the median of the membership in the module and gene significance for Mφ1/Mφ2 (Fig. 3C). Then, we screened the hub genes in the network based on the number of connections between nodes. CytoHubba was used to calculate the top 40 hub genes (Fig. 4B). Intersections between the two sets of hub genes were established. Intersections for the genes in blue, turquoise, and green modules were 15, 15, and 19, respectively (Fig. 4C). Finally, univariate analysis revealed that RPS9, DUT, and KIAA0430 genes were significantly associated with implantation failure in the blue, turquoise, and green modules (Fig. 4D).

3.4 Lab and external verification of hub genes expression and correlation with Mφ1/ Mφ2

We first validated the altered mRNA expression levels of DUT, RPS9, and KIAA0430 in endometrial tissues. Compared with the control group, the expression of all three genes (P = 0.029, 0.028, and 0.0006 for DUT, RPS9, and KIAA0430, respectively) in patients with IVF failure was significantly downregulated (Fig. 5A). Furthermore, we verified that the three mRNA-translated proteins had the same trend of alteration at protein levels (Fig. 5B).

Then, we performed correlation analysis by macrophages count and qPCR results to verify the association between the above screened genes and Mφ1/Mφ2. As shown in Fig. 6C, RPS9, DUT, and KIAA0430 genes exhibited significant correlation with Mφ1/ Mφ2, as 0.52 (P = 0.058), 0.50 (P = 0.048), and 0.55(P = 0.027), respectively (Fig. 5C).

Then, we selected Mφs datasets in GSE5099 to verify the association between the above screened genes and macrophage polarization. As shown in Fig. 5D, RPS9 (P = 0.038), DUT (P < 0.0001), and KIAA0430 (P < 0.0001) genes exhibited significant differences between Mφ1 and Mφ2. Subsequently, GSEA results for the three genes were consistent with our previous findings, that is ribosomal (RPS9, DUT and KIAA0430), lysosome (DUT), and proteasomal (RPS9 and KIAA0430), respectively (Fig. 5E).

3.5 Establishment and validation of an artificial neural network prognostic model

Then, we used a deep machine learning algorithm (ANN) to evaluate the predictive power of selected hub genes for defective endometrial receptivity. First, we set the number of units in the hidden layer of the GSE58144 dataset. The model achieved the highest prediction accuracy of 98.3% when the hidden layer was set to 24 (Fig. 6A).

Further, we validated this prognostic model using the GSE165004 dataset. Aforementioned hub genes were incorporated into the ANN prognostic model. Based on an earlier finding, we set the number of units in the hidden layer to 24 and performed 100 cycles of simulation to improve the prognostic accuracy (Fig. 6B). As a further verification that the ANN model is the optimal model for implantation failure prediction, we also established a prediction model based on logistic regression. In Fig. 6C, expression levels of the three hub genes were found to be significantly different between the two groups, consistent with the results of the GSE58144 dataset. Also, the risk score, based on the logistic regression, showed the ability for defective endometrial receptivity identification. The Nomogram based on logistic regression is shown in Fig. 6D. Prognostic results were averaged and plotted as ROC, and were found to compare to those by modules and hub genes. As result, the AUC of the ANN model was 0.975 (95% CI 0.945-1), significantly better than that of nomogram model, DUT, RPS9, and KIAA0430 (p = 0.0439, 9.27×10^− 6, 1.15×10^− 6, and 5.33×10^− 4), as shown in Fig. 7A. Sensitivity, specificity, YI, PPV, and NPV of the ANN model were 89.58 (95% CI 77.3–96.5), 95.83 (95% CI 78.9–99.9), 0.854, 97.7 (95% CI 86.3–99.7), and 82.1 (95% CI 66.6–91.4), superior to those of the nomogram model and other hub genes, as shown in Fig. 7B.

The altered immune microenvironment in the uterine cavity, an essential aspect of endometrial receptivity, is a crucial factor for a successful IVF(24). However, immune-related factors have been seriously underestimated when used to assess endometrial receptivity. In this study, we found that the balance of the ratio between Mφ1 and Mφ2 is an important factor that affects endometrial receptivity for patients. Then, we screened different biological functional gene modules and obtained hub genes that may associate with Mφ1/Mφ2 and developed a machine deep learning prognostic model with an excellent predictive performance through the artificial neural network approach.

Macrophages can be classified into the activated Mφ1 and the alternatively activated Mφ2 based on phenotypes and function. The balance of macrophage ratios is important for WOI and embryonic implantation in the endometrium(6). However, the type of macrophage dominating the mid-secretory phase has not yet been established. Russell et al. (25) reported a significant increase in Mφ2 during the luteal phase, which are thought to contribute to the establishment of maternal immune tolerance to fetal antigens at the onset of implantation. Conversely, Diao et al. (24) reported that the abundance of Mφ2 was significantly lower in control subjects than in patients with failed embryonic implantation. Some studies have also suggested that, during implantation, activated Mφ1 produce inflammatory cytokines and mediators, such as IL-6, IL-1β, TNF-α, and nitric oxide, inducing a pro-inflammatory response and promoting embryonic attachment to the decidua. Consistent with these observations, in this study, immune infiltration analysis and immunofluorescence validation revealed Mφ1/Mφ2 alterations in patients in the reduced endometrial receptivity group compared to the control group. This finding confirms the importance of the balance of macrophage polarization for embryonic implantation.

Factors that induce macrophage infiltration into the endometrium in mid-secretory phase, including chemokines, colony-stimulating factor (CSF)-1, and granulocyte macrophage-colony-stimulating factor (GM-CSF) are abundantly secreted by endometrial stromal cells in response to hormonal stimulation. Vascular endothelial growth factor receptor-1 (VEGFR-1) performs an important role in macrophage recruitment and angiogenesis at the implantation site. Macrophages have a high potential for plasticity and can modify their functions depending on the changing microenvironment in tissues and are, therefore, involved in different physiological functions and disease development. This study identified three highly correlated modules (antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP − dependent, phagosome acidification, and interleukin-1-mediated signaling pathway) in macrophage polarization in the endometrium that affect macrophage functions and their effects on endometrial receptivity through different aspects. The main pathways of the modules were ribosome, lysosome and proteasome pathways. Then, we selected representative genes (DUT, RPS9 and KIAA0430) in each gene module. These genes and their closely related are important in the homeostasis of macrophage polarization and maintenance of favourable endometrial receptivity (Fig. 8).

In blue module, the principal biological functions performed by the gene cluster is antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP − dependent. MHC class I plays an intermediate role in the regulation of macrophage phagocytosis. As part of local immunity, macrophages recognize antigens in MHC class I and present them to T cells which recognize the MHC-antigen complex through their T cell receptors, which requires additional costimulatory and cytokine signals. Macrophages can thus mediate and provide the costimulatory signals and cytokine secretion required for effective T-cell activation (26). Mφ1 secrete various cytokines, including tumor necrosis factor alpha (TNF-α), interleukin (IL)-1α, and IL-6, which may be closely associated with endometrial decidualization and embryonic implantation process. Mφ2 maturation can be induced by various cytokines, including IL-4, IL-13, glucocorticoids, as well as M-CSF/CSF1, and together with Tregs, they promote tissue remodeling and regeneration, wound healing, and anti-inflammation in the endometrial tissue. Altered expression levels of some cytokines, such as TNF-α and IL-1β, can affect the polarization shift between Mφ1 and Mφ2, impacting reproductive outcomes. Meanwhile, class I expression may be critical for avoidance of immunological rejection(27). Its hub gene RPS9, a component of the ribosomal 40s subunit, is an important factor in macrophage activation and polarization through its protein translation function. Various cytokines, such as IFN-γ, IL-10a, and IL-6 have been shown to regulate ribosomes through receptors, thereby affecting macrophage metabolism as well as mRNA translation (28). Silvia et al. (29) evaluated endometrial receptivity-related pathways and found that the ribosomal pathway is the most relevant for endometrial fertility. Furthermore, the result of GSEA showed a negative correlation between RPS9 and proteasome. According to Fan et al. (30), molecules such as ribosomes in macrophages can be hydrolyzed by the proteasome through ubiquitination. In addition, inhibition of this pathway has also been shown to have positive implications in protection of macrophages as well as the periphery from hypoxia-reoxygenation injury. In this study, RPS9 was significantly differentially expressed in IF and control groups for both GSE58144 and GSE165004 datasets. These findings imply that module of antigen processing and presentation of exogenous peptide antigen via MHC class I in macrophage polarization play important roles in endometrial receptivity.

As for the turquoise module, the genes are enriched in the biological process of phagosome acidification, which is directly linked to the phagocytic function of macrophages. DUT, as the hub gene, were enriched in lysosome and ribosome pathways by its related genes in the GSE5099 dataset, consistent with its module. Phagocytosis and antigen presentation by macrophages is dependent on phagosome and lysosomal activity. Phagosome-like compartments containing antigen at some stages fused with lysosomes to form a phagolysosome. Correspondingly, according to Ariza et al. (31), the regulation of DUT to macrophages depends on their phagocytosis. DUT (dUTPase), an essential enzyme during nucleotide metabolism, hydrolyzes dUTP to dUMP and pyrophosphate. Alternative splicing of this gene leads to different isoforms that are localized in the mitochondria or nucleus. dUTPase modulates innate immunity in human primary monocyte-derived macrophages through toll-like receptor (TLR) 2 leading to NF-kB activation and the production of pro-inflammatory cytokines. This process is achieved via macrophage phagocytosis of exosomes containing dUTPase. Then, interferon (IFN)-γ, lipopolysaccharide (LPS), or granulocyte macrophage-colony-stimulating factor induce macrophage maturation towards the Mφ1, which activates Toll-like receptor signaling pathways, thereby playing a crucial role in clearance of residual fibers and tissue debris, and in the synthesis of pro-inflammatory cytokines and growth factors. In addition, Wang et al(31). reported that DUT activation led to the production of large amounts of mtDNA, which bound the mitochondrial ribosomal proteins to co-synthesize mitochondria-associated proteins, further demonstrating the association between DUT and the ribosomal pathway. Therefore, enrichment analysis of the turquoise gene module in endometrial tissues and the DUT in macrophage GSE5099 dataset revealed a close association among DUT, phagosome acidification, and the aforementioned two pathways, which affect the Mφ1/Mφ2 ratio in patients who were in the mid-secretory phase.

The green module's biological function was enriched in the interleukin-1-mediated signaling pathway, which is vital for macrophage functions and polarization. Elevated expression levels of cytokines such as TNF-α in the early stages promote macrophages polarization, as well as the synthesis and secretion of IL-1 by macrophages. In addition, IL-1 in the early stages can act on interleukin-1-mediated signaling pathway to promote other cytokines’ synthesis and secretion,is intimately involved in the regulation of the ribosome pathway (32). Its hub gene, KIAA0430, also known as meiosis regulator and mRNA stability factor 1 (MARF1), encodes a putative peroxisomal protein, which can silence targeted mRNA and inhibit gene expression (33). It can regulate the translational function of ribosomes in the manner above, consistent with the result that it was enriched in the ribosome pathway in the GSE5099 dataset. Meanwhile, we found that gene enrichment in the green module, where KIAA0430 is located, was closely related to negative regulation of the proteasome pathway. The negative correlation between KIAA0430 and the proteasome pathway was verified in the macrophage GSE5099 dataset. Proteasomes are protein disruption devices that are involved in many essential cellular functions, such as cell cycle regulation, cell differentiation, signal transduction pathways, antigen processing for proper immune responses, stress signaling, inflammatory responses, and apoptosis. Moreover, they are involved in macrophage polarization. Cytokines such as IFN-γ and TNF-α have been shown to regulate macrophage functions by moderating the proteasome pathway. Han et al. (34) reported that proteasomes in endometrial stromal cells can mediate diminished protein stabilities of HOXA10, a histone important for promoting endometrial decidualization, leading to a defective endometrial receptivity. In contrast, negative regulation of the proteasome pathway by KIAA0430 facilitates the maintenance of favourable endometrial receptivity.

Although macrophage polarization is of great importance in pregnancy, few are available to diagnose reproductive prognosis, mainly because of subjectivity of cell counts by immunohistochemistry and difficulty with standardization. Moreover, immune cell infiltration dynamics leads to the possibility of some systematic errors. Therefore, in this study, we used macrophage polarization-related factors as predictors for reproductive prognostic outcomes. To ensure that genes included in the prognostic model are representative of certain molecular functions or structures, we performed gene selection through the WGCNA approach. Then, we used ANN and logistic regression approaches to establish a prognostic model. Traditional models were often constructed through the regression approach. Model optimization requires the inclusion of factors that are both normally distributed and independent. For factors with small values, the weight of their impact on the outcome is often ignored because of their low coefficients. The ANN model, however, can weigh each feature variable according to its importance and then perform the summation of activity function, thus, it requires fewer factors to provide a more accurate classification. Also, it was certainly demonstrated in our study that the ANN was superior to the traditional logistic regression-based Nomogram in terms of IVF failure prediction. Although we found an imbalance between Mφ1 and Mφ2, we performed the diagnosis by expression of factors in different gene modules instead of direct microscopic counts to improve reproducibility of the results and reduce measurement errors associated with subjectivity. The predictive accuracy of the model was up to 98.3% in the GSE58144 dataset, while its AUC was 0.975 in the validation GSE165004 dataset, significantly better than the predictive method with a single factor and regression approach. Its performance was also superior to the conventional endometrial receptivity examinations, such as ultrasound for endometrial thickness (> 7 mm) with a sensitivity of 99% and a specificity of 3% (2). The significance of this predictive method is that it effectively determines the receptive status of the endometrium. Jena et al. (6) reported that there is a dynamic balance of macrophage polarization patterns in the endometrium in response to the menstrual cycle, suggesting that this model is appropriate for assessing endometrial receptivity in patients with IVF failure due to WOI changes. Importantly, uterine cavity procedures (induced abortion, adhesiolysis, or polypectomy), inflammatory responses and transient immune dysregulation (dysbacteriosis or endometritis) due to stress can lead to embryonic implantation failure (35, 36). However, this model is highly correlated with the immune environment and is appropriate for informing optimal timing of natural or IVF pregnancies in such patients. For patients, assessment of endometrial receptivity using this tool is beneficial in identifying the cause and predicting the immediate outcome of IVF.

This study has its advantages. Firstly, the effects of 22 immune cell infiltrations on embryo implantation were explored for the first time, screening factors of significance. Secondly, instead of applying the commonly used PCA and other dimensionality reduction methods, we performed a hub gene selection in the module by WGCNA, which ensured the representativeness of the characteristics data and reduced the number of items to be examined. This dramatically saves tests cost in clinical application. Finally, a non-regression model was used to establish predicted reproductive outcomes, reducing the influence of data characteristics and weights on prediction. This experiment preliminary validated changes in the ratio between Mφ1 and Mφ2 and established a prognostic model, however, there are some limitations. First, other cheaper assays, such as qRT-PCR or ELISA, were not performed for the test of the model, which needs to be validated through further large sample experiments. Second, for the mechanistic study, only enrichment analysis was performed, which are all based on mRNA expression levels. Therefore, it was not possible to observe the effects of changes in translational or post-translational protein levels for mechanism prediction, which should be validated further. Importantly, Lab experiments results and cross-validation between multiple data sets guarantee reliability of the mechanism as well as the ANN model. In short, based on our results, the consistency of validation, as analysed by different methods and datasets, was mutually confirmed.

In summary, the balance between Mφ1 and Mφ2 is essential for the pregnancy process. Gene modules associated with biological processes of antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP − dependent, phagosome acidification, and interleukin-1-mediated signaling pathway can impact macrophage polarization that ameliorate endometrial receptivity. Furthermore, the established ANN model based on hub genes can effectively assess endometrial receptivity to inform patients' reproductive outcomes and individualized clinical management strategies.

WOI	window of implantation
(R)IF	(repeated) implantation failure
NK	Natural killer
WGCNA	weighted gene co-expression network analysis
ANN	artificial neural network
DER	decreased endometrial receptivity
RMA	Robust multi-array average
BSA	bovine serum albumin
PBS	phosphate-buffered saline
GO	Gene Ontology
BP	Biological Process
CC	Cellular Component
MF	Molecular Function
KEGG	Kyoto Encyclopedia of Genes and Genome
GSEA	Gene-set enrichment analysis
qRT-PCR	Quantitative Real-Time PCR Analysis
ROC	receiver operating characteristic
PPV	positive predictive value
NPV	negative predictive value
MS	module significance
CSF-1	colony-stimulating factor
GM-CSF	granulocyte macrophage-colony-stimulating factor
VEGFR-1	Vascular endothelial growth factor receptor-1
IFN-γ	Interferon
LPS	lipopolysaccharide
TNF-α	tumor necrosis factor alpha
IL-1α	interleukin.

Competing interests

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.

Author Contributions

BL, Conception and design, Collection and assembly of data, Manuscript writing, Final approval of manuscript, Administrative support; HD, Conception and design, Manuscript writing, Administrative support, Final approval of manuscript; SW, Provision of study materials or patients, Manuscript writing, Final approval of manuscript; JW and YL Manuscript writing, Final approval of manuscript.

Funding

This study was financially supported by the National Key Research and Development Program of China (2018YFC1004803).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics approval and consent to participate

Ethical approval for this study was obtained from the Research Ethics Committee of the Beijing Obstetrics and Gynecology Hospital. Experiments were performed (under protocol number 2017-KY-082-02) based on the Helsinki Declaration of 1975 (revised in 2013).

Acknowledgments

This study was financially supported by the National Key Research and Development Program of China (2018YFC1004803). The paper is edited by Home for researchers editorial team. We also thank the researchers of datasets GSE58144, GSE5099 and GSE165004 for their selfless sharing.

Consent for publication

All presentations have consent for publication.

1. Lessey BA, Young SL. What exactly is endometrial receptivity? Fertility and sterility. 2019;111(4):611-7.

2. Craciunas L, Gallos I, Chu J, Bourne T, Quenby S, Brosens JJ, et al. Conventional and modern markers of endometrial receptivity: a systematic review and meta-analysis. Hum Reprod Update. 2019;25(2):202-23.

3. Ruiz-Alonso M, Blesa D, Díaz-Gimeno P, Gómez E, Fernández-Sánchez M, Carranza F, et al. The endometrial receptivity array for diagnosis and personalized embryo transfer as a treatment for patients with repeated implantation failure. Fertility and sterility. 2013;100(3):818-24.

4. Ben Rafael Z. Endometrial Receptivity Analysis (ERA) test: an unproven technology. Hum Reprod Open. 2021;2021(2):hoab010.

5. Zhang Y-H, He M, Wang Y, Liao A-H. Modulators of the Balance between M1 and M2 Macrophages during Pregnancy. Front Immunol. 2017;8:120.

6. Jena MK, Nayak N, Chen K, Nayak NR. Role of Macrophages in Pregnancy and Related Complications. Archivum immunologiae et therapiae experimentalis. 2019;67(5):295-309.

7. Li Z-H, Wang L-L, Liu H, Muyayalo KP, Huang X-B, Mor G, et al. Galectin-9 Alleviates LPS-Induced Preeclampsia-Like Impairment in Rats via Switching Decidual Macrophage Polarization to M2 Subtype. Front Immunol. 2018;9:3142.

8. Schork NJ. Genetics of complex disease: approaches, problems, and solutions. Am J Respir Crit Care Med. 1997;156(4 Pt 2):S103-S9.

9. Lan H, Chen M, Flowers JB, Yandell BS, Stapleton DS, Mata CM, et al. Combined expression trait correlations and expression quantitative trait locus mapping. PLoS Genet. 2006;2(1):e6.

10. Li B, Wang S, Duan H, Wang Y, Guo Z. Discovery of gene module acting on ubiquitin-mediated proteolysis pathway by co-expression network analysis for endometriosis. Reprod Biomed Online. 2021;42(2):429-41.

11. Hou Q, Bing ZT, Hu C, Li MY, Yang KH, Mo Z, et al. RankProd Combined with Genetic Algorithm Optimized Artificial Neural Network Establishes a Diagnostic and Prognostic Prediction Model that Revealed C1QTNF3 as a Biomarker for Prostate Cancer. EBioMedicine. 2018;32:234-44.

12. Zhang Z. A gentle introduction to artificial neural networks. Ann Transl Med. 2016;4(19):370.

13. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004;20(3):307-15.

14. Li Y, Bai W, Zhang L. The Overexpression of CD80 and ISG15 Are Associated with the Progression and Metastasis of Breast Cancer by a Meta-Analysis Integrating Three Microarray Datasets. Pathol Oncol Res. 2020;26(1):443-52.

15. Yoneya T, Miyazawa T. Integration of pre-normalized microarray data using quantile correction. Bioinformation. 2011;5(9):382-5.

16. Shakya K, Ruskin HJ, Kerr G, Crane M, Becker J. Comparison of microarray preprocessing methods. Adv Exp Med Biol. 2010;680:139-47.

17. Forero DA, Guio-Vega GP, González-Giraldo Y. A comprehensive regional analysis of genome-wide expression profiles for major depressive disorder. J Affect Disord. 2017;218:86-92.

18. Li B, Duan H, Wang S, Li Y. Ferroptosis resistance mechanisms in endometriosis for diagnostic model establishment. Reprod Biomed Online. 2021.

19. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. NCBI Resource Coordinators, 2015;43(Database issue):D6-17.

20. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37(7):773-82.

21. Ghazalpour A, Doss S, Zhang B, Wang S, Plaisier C, Castellanos R, et al. Integrating genetic and network analysis to characterize genes related to mouse weight. PLoS Genet. 2006;2(8):e130.

22. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

23. Wang S, Li B, Duan H, Wang Y, Shen X, Dong Q. Abnormal expression of connective tissue growth factor and its correlation with fibrogenesis in adenomyosis. Reprod Biomed Online. 2021;42(3):651-60.

24. Diao L, Cai S, Huang C, Li L, Yu S, Wang L, et al. New endometrial immune cell-based score (EI-score) for the prediction of implantation success for patients undergoing IVF/ICSI. Placenta. 2020;99:180-8.

25. Russell P, Anderson L, Lieberman D, Tremellen K, Yilmaz H, Cheerala B, et al. The distribution of immune cells and macrophages in the endometrium of women with recurrent reproductive failure I: Techniques. Journal of reproductive immunology. 2011;91(1-2):90-102.

26. Guerriero JL. Macrophages: Their Untold Story in T Cell Activation and Function. Int Rev Cell Mol Biol. 2019;342:73-93.

27. Davies CJ, Fisher PJ, Schlafer DH. Temporal and regional regulation of major histocompatibility complex class I expression at the bovine uterine/placental interface. Placenta. 2000;21(2-3):194-202.

28. Su X, Yu Y, Zhong Y, Giannopoulou EG, Hu X, Liu H, et al. Interferon-γ regulates cellular metabolism and mRNA translation to potentiate macrophage activation. Nat Immunol. 2015;16(8):838-49.

29. Pérez-Debén S, Bellver J, Alamá P, Salsano S, Quiñonero A, Sebastian-Leon P, et al. iTRAQ comparison of proteomic profiles of endometrial receptivity. J Proteomics. 2019;203:103381.

30. Fan T, Huang Z, Wang W, Zhang B, Xu Y, Mao Z, et al. Proteasome inhibition promotes autophagy and protects from endoplasmic reticulum stress in rat alveolar macrophages exposed to hypoxia-reoxygenation injury. J Cell Physiol. 2018;233(10):6748-58.

31. Wang Y-P, Sharda A, Xu S-N, van Gastel N, Man CH, Choi U, et al. Malic enzyme 2 connects the Krebs cycle intermediate fumarate to mitochondrial biogenesis. Cell Metab. 2021;33(5).

32. Li B, Duan H, Wang S, Wang Y, Chang Y, Guo Z, et al. Hierarchical cluster analysis in the study of the effect of cytokine expression patterns on endometrial repair and receptivity after hysteroscopic adhesiolysis. Ann Transl Med. 2021;9(9):746.

33. Nishimura T, Fakim H, Brandmann T, Youn J-Y, Gingras A-C, Jinek M, et al. Human MARF1 is an endoribonuclease that interacts with the DCP1:2 decapping complex and degrades target mRNAs. Nucleic Acids Res. 2018;46(22):12008-21.

34. Han K, Zhou Q. Skp2 Deteriorates the Uterine Receptivity by Interacting with HOXA10 and Promoting its Degradation. Reprod Sci. 2021;28(4):1069-78.

35. Ke L, Lin W, Liu Y, Ou W, Lin Z. Association of induced abortion with preterm birth risk in first-time mothers. Sci Rep. 2018;8(1):5353-.

36. Zhao X, Zhao Q, Zhu X, Huang H, Wan X, Guo R, et al. Study on the correlation among dysbacteriosis, imbalance of cytokine and the formation of intrauterine adhesion. Ann Transl Med. 2020;8(4):52.

SupplementaryTable.docx

Download PDF

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Artificial Neural Network Model Using Immune-infiltration Modules for Endometrial Receptivity Assessment of Implantation Failure

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methods

2.1 Datasets and patient selection

2.2 Processing of primary datasets

2.3 Digital deconvolution of bulk tissues

2.4 Immunofluorescence

2.5 Weighted gene co-expression network

2.6 Enrichment analysis of functional categories

2.7 Quantitative Real-Time PCR Analysis (qRT-PCR)

2.8 Western blot

2.9 prognostic model establishment and statistics

3. Results

3.1 Immune infiltration levels on endometrial receptivity

3.2 Macrophage polarization-related gene module functions

3.3 Selection and verification of hub genes associated with Mφ1/Mφ2

3.4 Lab and external verification of hub genes expression and correlation with Mφ1/ Mφ2

3.5 Establishment and validation of an artificial neural network prognostic model

4. Discussion

Abbreviations

Declarations

Competing interests

Author Contributions

Funding

Data Availability Statement

Ethics approval and consent to participate

Acknowledgments

Consent for publication

References

Supplementary Files

Status:

Version 1