Meta-analysis of high-throughput transcriptomics uncovers aberrant endometrial microRNA-messengerRNA networks in recurrent implantation failure

Data acquisition

The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) was systematically interrogated to identify mRNA and miRNA expression datasets pertinent to RIF. The search strategy utilized the following key terms: “Recurrent implantation failure,” “Homo sapiens [porgn: txid9609],” and the methodology filters “Expression profiling by array,” “Non-coding RNA profiling by array,” and “Non-coding RNA profiling by high throughput sequencing.”

Two mRNA expression datasets, GSE111974 and GSE26787, were obtained from the GEO database. The GSE111974 dataset, comprising 48 endometrial tissue samples (24 RIF and 24 controls), was generated using Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray. The second dataset (GSE26787) was generated using Affymetrix Human Genome U133 Plus 2.0 Array and included 10 endometrial tissue samples (5 RIF and 5 controls).

For miRNA, two datasets (GSE71332 and GSE108966) were considered. The GSE71332 dataset was generated on 12 endometrial tissue samples (7 RIF and 5 controls) with Agilent-046064 Unrestricted_Human_miRNA_V19.0_ Microarray. The second miRNA dataset, GSE108966, was generated on Illumina HiSeq 2500 (miRNA V21.0) using 110 endometrium tissue samples (36 RIF and 74 controls). The summary of the datasets is presented in Table 1.

RNA sequencing data

The R statistical programming language (version 4.5.0, 2025-04-11; Bioconductor version 3.2.1) was used for data processing and integration. For mRNA analysis, expression data from GSE111974 and GSE26787 were log₂-transformed and quantile-normalized using the limma package.^[8] Probes were aggregated to gene-level expression through median values based on platform-specific annotation files.^[9] Genes with variance below the 25^th percentile across samples were filtered out to reduce noise.^[10,11] Batch effects (non-biological differences) were corrected using the ComBat function from the surrogate variable analysis package.^[12]

Pre-processing of the GSE71332 miRNA dataset was performed similarly to the mRNA datasets. For the second miRNA dataset (GSE108966), raw count data were retrieved, and DESeqDataSet (a central data structure of the DESeq2 package) was created. The ComBat-Seq function was used for batch correction.^[13]

To identify differentially expressed genes (DEGs) between RIF and controls, differential gene expression analysis was performed separately for each dataset using limma and DESeq2, as applicable.

Meta-analysis of differentially expressed RNAs

Weighted gene co-expression network analysis (WGCNA)

WGCNA, version 1.73 (https://rdocumentation.org/packages/WGCNA/versions/1.73) package in the R language was used to analyze gene modules and genes that correlated highly with RIF. The combined, batch-corrected mRNA expression matrix from GSE111974 and GSE26787 was used for WGCNA. Quality control was performed using the good Samples Genes function from the WGCNA package. A soft-thresholding power was selected based on scale-free topology fit (R² > 0.85), with a power of 9.^[18] An unsigned adjacency matrix was computed, followed by the topological overlap matrix calculation. Hierarchical clustering with dynamic tree cutting identified initial modules, which were merged based on eigengene dissimilarity (cutHeight = 0.25). Module– trait correlations were calculated for RIF and control traits, and the top 10% most connected genes per module were identified.^[19]

miRNA–mRNA target integration and network construction

Validated miRNA–mRNA interactions were retrieved using the multiMiR package, restricting targets to high-confidence mRNAs and miRNAs.^[20] Only interactions with consistent directionality (negative correlation between miRNA and mRNA logFC) were retained. A network of these miRNA–mRNA interactions was visualized in Cytoscape (version 3.10.2) (https://cytoscape.org/). Further, the overlapping hub genes from WGCNA modules and the miRNA–mRNA interactions were taken into consideration.^[18]

Gene ontology (GO) analysis

GO enrichment analysis was performed using the org.Hs.eg. db annotation package from Bioconductor. GO terms were mapped based on the GO database (release: February 6, 2025) and Entrez Gene annotations (sourced: February 22, 2025). The GO file was obtained from the GO consortium (http://current.geneontology.org/ontology/go-basic.obo). All annotations and mappings were performed using R version 4.5.0 and Bioconductor version 3.2.1.

Pathway enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was performed in R using clusterProfiler (v4.16.0), converting gene symbols to ENTREZ IDs through org.Hs.eg.db (v3.21.0).^[21] The enrichKEGG() function was applied with organism= “has,” P-value cutoff = 0.05, and q-value cutoff = 0.2, utilizing the most current KEGG database (2024) accessed through the KEGG REST API.^[22]

Figure 1 represents summary of the study (1a) and the detailed workflow (1b).

Close

Figure 1

(a): Graphical abstract of the present study. Created in BioRender. Paul, P. (2025) https://BioRender.com/vam5ze8

Export to PPT

Close

Figure 1

(b): Overall workflow of the study. (mRNA: messenger RNA, miRNA: microRNA, RMA: Robust Multi-array Average, DE: Differentially expressed, FDR: False discovery rate, WGCNA: Weighted Gene Co-expression Network Analysis, KEGG: Kyoto Encyclopedia of Genes and Genomes).

Export to PPT

Modules and genes screened by WGCNA

A total of 7,118 highly correlated genes were grouped into eight distinct modules using WGCNA. Each module was represented by a unique color: black, blue, brown, dark gray, dark red, gray, light green, and purple. Among these, four modules (black, brown, dark gray, and light green) were identified as significant, with adjusted P > 0.05. A gene dendrogram illustrating the clustering of modules is shown in Figure 4.

Close

Figure 4:

WGCNA results. (a) Cluster dendrogram with dynamic tree cut and merged dynamic modules, showing gene clustering. (b) Module–trait relationships heatmap with BH-adjusted P-values, indicating correlations between modules (e.g., MEblack, MEblue) and recurrent implantation failure/controls traits (red: Positive, blue: Negative). WGCNA: Weighted gene co-expression network analysis, ME: Module Eigengene.

Export to PPT

miRNA–mRNA network construction

A total of 2,641 interactions could be observed in the miRNAs-target network (14 miRNAs and 1,086 targets/ DEGs). Next, we selected the top 500 unique targets based on a combination of score-based sorting and constructed a representative subnetwork of 728 interactions [Figure 5]. Twelve hub genes (LNX2, DUOX1, MFSD6, RANBP17, PEA15, CD38, SLC7A7, UBAP1L, TYROBP, LAPTM5, SRGN, and GPR65) were identified from the intersection between WGCNA-derived genes (n = 70) and 1,086 unique miRNA targets. The complete list of miRNAs and hub genes is provided in Tables 2 and 3. As visualized from the network [Figure 5a and b], hsa-miR-101-3p, hsa-miR-27b-3p, hsa-miR-30e-5p, hsa-miR-301a-3p, hsa-miR-19b-3p, hsa-miR-455-3p, hsa-miR-27a-3p, hsa-miR-449a, hsa-miR-424-5p, hsa-miR-362-3p, hsa-miR-143-3p, hsa-miR-141-3p, hsa-miR-15a-5p, and hsa-miR-6724-5p were found to be associated with GPR65, CD38, SLC7A7 (upregulated hub genes) and LNX2, DUOX1, MFSD6, RANBP17, PEA15 (downregulated hub genes). However, according to the combination score-derived ranking, the remaining four genes (UBAP1L, TYROBP, LAPTM5, and SRGN) do not appear within those 500 interactions.

Close

Figure 5:

miRNA-mRNA interaction network. (a) Down-regulated miRNAs (red diamonds) interacting with up regulated mRNAs (yellow ovals), along with the overlapping WGCNA hub genes/mRNA (blue arrowheads). (b) miRNA-mRNA interaction network. Up-regulated miRNAs (red diamonds) interacting with down regulated mRNAs (yellow ovals), along with the overlapping WGCNA hub genes/mRNA (green arrowheads). mRNA: Messenger RNA, miRNA: MicroRNA.

Export to PPT

GO analysis

GO enrichment analysis of the DEGs revealed significant clustering in cellular components such as the early endosome, apical part of the cell, cell cortex, and various synaptic and membrane-associated structures, suggesting roles in intracellular trafficking and cell polarity [Figure 6a]. For molecular function, cell adhesion and cell–cell adhesion mediator activity were notably enriched, highlighting the importance of adhesion mechanisms in endometrial biology [Figure 6b]. In terms of biological processes, the most prominent enrichments included wound healing, cell-substrate adhesion, regulation of actin filament-based processes, and several coagulation-related pathways, collectively indicating that disruptions in cellular organization, adhesion, and tissue repair may contribute to the molecular pathology of RIF [Figure 6c].

Close

Figure 6:

GO enrichment associated with the DEGs. (a) Cellular components (b) Molecular functions (c) Biological processes (d) KEGG pathway; pathways with gene ratio (number of input genes associated with the term/total input genes) and p.adjust; circle size represents gene count, color intensity reflects p.adjust values (red: Lower, blue: Higher). GO: Gene ontology, KEGG: Kyoto encyclopedia of genes and genomes, TNF: Tumor Necrosis Factor, TOR: Target of Rapamycin, GTP: Guanosine triphosphate.

Export to PPT

Identification of pathways significantly associated with differentially expressed mRNAs

KEGG pathway enrichment analysis revealed considerable involvement of key signaling pathways [Figure 6d]. The analysis showed significant enrichment of 15 pathways (adjusted P < 0.05), with the PI3K-Akt signaling pathway emerging as the most prominently enriched pathway (GeneRatio = 0.07, adjusted P < 0.01, gene count = 35). Other significantly enriched pathways included human papillomavirus infection, apelin signaling, focal adhesion, insulin resistance, Ras, and mTOR signaling. Notably, leukocyte transendothelial migration and TNF signaling pathways, directly relevant to endometrial biology and implantation, could be identified. It is well established that these pathways are associated with immune modulation and inflammatory responses critical for embryo-endometrial interactions. The Hippo signaling pathway and focal adhesion pathway were also found to be enriched, suggesting involvement in cell adhesion and tissue remodeling processes, essential for endometrial receptivity. These results collectively highlight the central role of cell survival, proliferation, and vascularization pathways in the pathogenesis of RIF, with the PI3K-Akt pathway representing a key regulatory node in RIF-associated molecular dysfunction.

Limitations

This study is associated with a few limitations. First, identification of potential biomarkers is primarily driven by statistical analysis; validation in a patient cohort could not be performed. Second, the number of relevant studies publicly available is limited, and the sample size within these studies is relatively small. This constraint could have affected the statistical power and generalizability of our findings. These limitations underline the importance of future studies with a large, well-matched patient population to better elucidate the complex molecular mechanisms underlying RIF.