Construction of circRNA‐miRNA‐mRNA network in the pathogenesis of recurrent implantation failure using integrated bioinformatics study

Abstract This research attempted to elucidate the molecular components are involved in the pathogenesis of recurrent implantation failure (RIF). We initially identified that 386 mRNAs, 144 miRNAs and 2548 circRNAs were differentially expressed (DE) in RIF and then investigated the genetic cause of the observed abnormal expression by constructing a circRNA‐miRNA‐mRNA network considering the competing endogenous RNA theory. We further analysed the upstream transcription factors and related kinases of DEmRNAs (DEMs) and demonstrated that SUZ12, AR, TP63, NANOG, and TCF3 were the top five TFs binding to these DEMs. Besides, protein‐protein interaction analysis disclosed that ACTB, CXCL10, PTGS2, CXCL12, GNG4, AGT, CXCL11, SST, PENK, and FOXM1 were the top 10 hub genes in the acquired network. Finally, we performed the functional enrichment analysis and found that arrhythmogenic right ventricular cardiomyopathy (ARVC), hypertrophic cardiomyopathy (HCM), pathways in cancer, TNF signalling pathway and steroid hormone biosynthesis were the potentially disrupted pathways in RIF patients. Optimistically, our findings may deepen our apprehensions about the underlying molecular and biological causes of RIF and provide vital clues for future laboratory and clinical experiments that will ultimately bring a better outcome for patients with RIF.


| INTRODUC TI ON
Recurrent implantation failure (RIF) has been defined as an abnormal condition in which the transfer of more than ten good-quality embryos frequently fails to implant following at least two cycles of in vitro fertilization (IVF) or intracytoplasmic sperm injection (ICSI). 1,2 Multiple aetiologies have been proposed for RIF occurrence including, advanced maternal age, increased body mass index, smoking status, stress levels, infectious disease, genetics and immunological factors. 3 Plenty of literature has highlighted insufficient endometrial receptivity and constitutive endometrial dysfunction as a leading contributor to this condition. 4 Moreover, many studies also have supported the idea that endometrial gene expression profiles may be altered in patients experiencing RIF. 5 However, the exact mechanisms are not yet fully apprehended and need further elucidation. Therefore, it is crucial to investigate the possible genetic cause of abnormal endometrial gene expression profiles in RIF's patients to find appropriate and effective targets for diagnostic and treatment purposes.
Non-coding RNAs (ncRNAs) are a type of gene without any protein production that modulates cell function through many different mechanisms. Functionally, ncRNAs have been classified into two categories: housekeeping and regulatory ncRNAs. The first ones (eg tRNAs, rRNAs and snRNAs) play crucial roles in many cellular processes. The later ncRNAs (eg miRNAs, siRNAs and lncRNAs) are tissue-specific and work in response to various internal and environmental stimuli. 6 Among ncRNAs, circular RNAs (CircRNAs) are newly discovered and have a specific structure recognized by a covalently closed-loop without 5'-end caps, or 3'-end poly (A) tails. 7 Several research pieces have indicated that circRNAs exert their biological tasks by regulating transcription, maturation of RNAs, methylation of histones, and function as microRNA (miRNA) sponges. 8 These shreds of evidence demonstrated the potential roles of circRNAs in physiological and pathological cellular conditions. 9,10 Interestingly, L Liu et al reported that the expressions of certain circRNAs were misregulated in RIF patients versus healthy controls. 11 In 2011, Salmena et al presented the competing endogenous RNA (CeRNA) theory, in which ncRNAs control the expression of corresponding mRNAs by targeting the shared miRNAs through miRNA response elements. 12,13 Up to now, CeRNA networks have been constructed for many diseases and osteoarthritis, 14 Diabetes Mellitus, 15 coronary artery disease, 16 lung cancer 17 and gastric cancer. 18 Here, we assumed that the establishment of the circRNA-miRNA-mRNA network could help us to boost our apprehension about the involved processes in the RIF.
In this research, first, we gathered the expression profiles of mRNAs, miRNAs and circRNAs in RIF patients from the Gene Expression Omnibus (GEO) database. Next, using GEO2R analysis, the Differentially Expressed mRNAs (DEMs), miRNAs (DEMIs) and circRNAs (DECs) were determined. Furthermore, CircRNA-miRNA and miRNA-mRNA axes were established. Then, by integrating them, a circRNA-miRNA-mRNA network was established. Moreover, an assessment of the upstream transcription factor (TFs) and Protein Kinases (PKs) of DEMs was performed. Besides, commonly used enrichment analysis approaches such as protein-protein interactions (PPI), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were involved in the prediction of the probable mechanisms and tasks of dysregulated genes and pathways in RIF.
The present work could advance our perception about the underlying molecular processes of RIF, providing novel diagnostic biomarkers, therapeutic targets and remarkable hints for further studies.

| Expression profile data
First, a total of three microarray expression profile data sets were selected from the National Center of Biotechnology Information (NCBI) GEO (https://www.ncbi.nlm.nih.gov/geo/) (   were ruled out. We converted the multiple expression values referred to the same gene to a single value using the following approach: Genes with the same direction of expression were averaged and discarded if they had the reversed expression tendency.

| Assessment of the reliability of GEO data sets
In order to check the reliability of GEO data sets included in our re-

| Construction of the miRNA-mRNA network
The miRDIP v4.1 online tool (http://ophid.utoro nto.ca/mirDI P/), as an integrated database of human microRNA-target predictions across from 30 independent resources with an integrative score, was used to identify the possible interactions between DEMs and DEMIs via bidirectional mode search under a high confidence filter (score class 5% = high). 32 It hastwo major approaches for computational prediction of the appropriate miRNA-mRNA pairs: identification of downstream gene targets for selected miRNA or upstream miRNA regulators for query gene. The second is to implement validation of the miRNA-gene correlation analysis using previous experiments. As miRNAs tend to control the expression of their desired mRNAs, even to a small extent, only the pair of miRNA-mRNAs with reverse expression patterns were included for miRNA-mRNA network construction. 33,34

| Establishment of the circRNA-miRNA network
Targeting circRNAs of DEMIs in miRNA-mRNA network was predicted through the circBank (http://www.circb ank.cn/) database. 35 Then, we screened these circRNAs by DECs in GSE14 7442. As the expression levels of circRNAs typically do not alter the expression of corresponding miRNAs, 36 the circRNA-miRNA pairs were selected without the restriction that the correlation between circRNAs and miRNAs must be negative.

| Reconstruction of circRNA-miRNA-mRNA network
We combined the pairs of circRNA-miRNA and miRNA-mRNA to reconstruct the circRNA-miRNA-mRNA network. The nodes that could not reach a circRNA-miRNA-mRNA axis were eliminated. The circRNA-miRNA-mRNA network was visualized utilizing Cytoscape 3.8.0. 37

| TFs and PKs related to DEGs
The upstream regulators and protein kinases related to 386 DEMs in GSE11 1974 were recognized by submitting their official gene names to Expression2Kinases (X2K) tool (https://amp.pharm.mssm.edu/ X2K/). 38 The top 10 most enriched TFs and PKs were listed regarding the integrated (p-value and z-score) score.

| Searching for a potential connection between DEMs, DEMIs and DECs
We studied the possible regulatory relationship between DEMIs in GSE71332 and DEMs in GSE11 1974 through the miRDIP v4.1 database. As depicted in Figure S1, the results presented 966 miRNA-mRNA pairs in which 87 DEMIs presumably sponged 170 DEMs.

| Construction of the circRNA-miRNA-mRNA network
We established a circRNA-miRNA-mRNA network by merging the circRNAs-miRNAs and miRNA-mRNA pairs and excluding the nodes that were not linked to the central network. As exhibited in Figure 2, the final network was comprised of 437 nodes and 1375 edges.

| Analysis of upstream regulator of DEMs
We recognized the enriched upstream TF, mediatory protein, and related PKs to elucidate the molecular processes of 386 DEMs Transcription Factor 3 (TCF3) were the top five TFs binding to DEMs ( Figure 3A). We also found that these TFs were related to 58 intermediate proteins ( Figure 3B). Our investigations then uncovered the prominent protein kinases such as mitogen-activated protein kinase 1 (MAPK1), cyclin-dependent kinase 4 (CDK4), homeodomain-interacting protein kinase 2 (HIPK2), Casein Kinase 2 Alpha 1 (CSNK2A1) and mitogen-activated protein kinase 3 (MAPK3) for these DEMs ( Figure 3C) had relations with a plenty number of mediatory proteins and TFs ( Figure 3D).

| Identifying hub genes in PPI network construction
PPI networks were constructed using the STRING database. The PPI network showed 279 nodes and 667 edges (Figure 4). The top 10 identified hub genes were ranked as follows ( Table 2)

| Enrichr database analysis
The DAVID v6.8 database was utilized for functional enrich-  Based on the CeRNA hypothesis, we sought to obtain the potential regulatory between the retrieved differently expressed genes. To the best of our knowledge, the present research was the first attempt to build a circRNA-miRNA-mRNA network and identify the most potential corresponding pathways in the pathogenesis of RIF. However, our study had some limitations. Firstly, due to the limited number of GEO data sets, we only have utilized the data obtained from three data sets deposited in the GEO database. Additionally, the sample size of these data sets was few.

| D ISCUSS I ON
Therefore, it was hard to prevent random errors. Besides, our analyses' findings were not confirmed with reliable laboratory analysis, and we just focused on in-silico assessments. Depending solely on bioinformatics investigations may trigger deviations. Hence, in vitro experiments remain to be carried out to validate the results of the present report.

| CON CLUS ION
In conclusion, the present research performed various bioinformatics analyses on the DEGs in the endometrium of RIF patients. These investigations demonstrated that the hub genes in RIF directly or indirectly affect cell cycle, apoptosis, angiogenesis and immune response-related pathways that modulate the endometrial receptivity and embryo implantation. Hopefully, our study's findings may bear an increased knowledge about the underlying molecular and biological aetiologies of RIF and valuable clues for additional experiments that will ultimately lead to improved capacity for accelerated and accurate diagnosis and treatment of patients with RIF.

ACK N OWLED G EM ENT
This study was supported by a grant (980305) from the Hormozgan University of Medical Sciences.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data generated\examined for the present research are presented in the manuscript context and the supplementary file.