Research Paper Volume 13, Issue 4 pp 6076—6090

Hsa_circ_0038383-mediated competitive endogenous RNA network in recurrent implantation failure

Huishan Zhao1, *, , Lili Chen1, *, , Yinghua Shan1, *, , Gang Chen2, , Yongli Chu3, , Huangguan Dai1, , Xuemei Liu1, , Hongchu Bao1, ,

  • 1 Reproductive Medicine Centre, The Affiliated Yantai Yuhuangding Hospital of Qingdao University, Yantai, China
  • 2 Department of Breast Surgery, The Affiliated Yantai Yuhuangding Hospital of Qingdao University, Yantai, China
  • 3 Department of Obstetrics and Gynecology, The Affiliated Yantai Yuhuangding Hospital of Qingdao University, Yantai, China
* Equal contribution

Received: September 5, 2020       Accepted: December 19, 2020       Published: February 20, 2021      

https://doi.org/10.18632/aging.202590
How to Cite

Copyright: © 2021 Zhao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Inadequate endometrial receptivity contributes to recurrent implantation failure (RIF) during IVF–embryo transfer. Though multiple circRNAs have been confirmed differentially expression in RIF, the potential function of novel circRNAs needed to be detected.

Results: The top ten DEcircRNAs were selected as initial candidates. A ceRNA network was conducted on the basis of circRNA–miRNA–mRNA potential interaction, consisting of 10 DEcircRNAs, 28 DEmiRNAs and 59 DEmRNAs. Three down-regulation circRNAs with high degree of connectivity were verified by RT-qPCR, and results suggested that only hsa_circ_0038383 was significantly downregulation in RIF compared with control group. Subsequently, three hub genes (HOXA3, HOXA9 and PBX1) were identified as hub genes. Ultimately, a subnetwork was determined based on one DEcircRNA (hsa_circ_0038383), two DEmiRNAs (has-miR-196b-5p and has-miR-424-5p), and three DEmRNAs (HOXA3, HOXA9 and PBX1). Following verification, hsa_circ_0038383/miR-196b-5p/HOXA9 axis may be a key pathway in affecting RIF.

Conclusion: In summary, a hsa_circ_0038383-mediated ceRNA network related to RIF was proposed. This network provided new insight into exploring potential biomarkers for diagnosis and clinical treatment of RIF.

Methods: We retrieved the expression profiles of RIF from GEO databases (circRNA, microRNA and mRNA) and constructed a competing endogenous RNAs (ceRNA) network based on predicted circRNA–miRNA and miRNA–mRNA pairs. The expression levels of three hub DEcircRNAs identified by cytoscape were validated by RT-qPCR.

Introduction

Recurrent implantation failure (RIF), defined as a clinical symptom that implanted good quality embryos experience at least three times failure during in vitro fertilization-embryo transfer (IVF–ET), remains a challenging problem for assisted reproductive technology (ART) [1, 2]. However, a precise definition of RIF is still argumentative, leading to accurately recognize patients difficulty [3, 4]. The pathogenesis of RIF is various, relating to disturbed immune system, poor quality embryos and insufficient endometrial receptivity, among which poor endometrial receptivity is believed as the decisive cause of RIF [5, 6]. Currently, the receptivity of endometrium is evaluated by morphological features, but the invasive damage to childbearing age women during laparoscopy is inevitable. Therefore, seeking noninvasive diagnostic and prognostic biomarkers of RIF is urgently necessary.

Circular RNAs (circRNAs), as novel class of endogenous non-coding RNAs (ncRNAs), are covalently closed loop structures generated by a process named back-splicing with high stability, abundance and tissue specificity [79]. The loop structure of circRNAs makes it more suitable for biomarkers than their corresponding linear transcripts [10]. Recently, numerous of studies have suggested that circRNAs regulate gene expression by competitively binding as miRNA sponges, or interacted with RNA-binding proteins (RBPs) [11]. Growing evidences revealed that circRNAs were correlated with the diseases processes, containing immune disorder, tumorigenesis, metabolic disorders, degenerative disease, and disease occurrence and progression, etc [1214]. Thus, circRNAs may become potential therapeutic target and biomarkers with diagnostic capabilities. Multiple studies have reported that the gene expression of endometrium is aberrant in RIF patients compared with healthy women, involving in coding RNA and non-coding genes [5, 1517]. Yet, reports about focused on the differential expression of circRNAs in RIF are few. Only one study reported a set of differently expressed circRNAs (DEcircRNAs) in female with RIF for the first time [18]. Hence, discovering more RIF-related DEcircRNAs and exploiting their responding mechanism to understand pathogenesis and filter diagnostic biomarkers of RIF is necessary.

The competing endogenous RNA (ceRNA) hypothesis was described as ncRNAs regulating mRNA expression through competitively binding to shared miRNAs and forming a regulatory RNA network [19]. In recent years, it has been widely used in seeking candidate gene and exploring new mechanisms of circRNAs by constructing ceRNA regulated network based on bioinformatics analysis. Lu et al. found that hsa_circ_0011385 may play a key role in carcinogenesis-related pathways of bladder cancer as a ceRNA by identifying a circRNA–miRNA–mRNA regulatory network [20]. Another study disclosed that hsa_circ_0028883 could be regarded as a potentially reliable biomarker to diagnose active tuberculosis via comprehensive bioinformatics analysis and further verification with real-time quantitative polymerase chain reaction (RT-qPCR) [21]. In this study, we discussed the distinct circRNAs and their molecular mechanisms in RIF with the assistance of public databases. Firstly, we downloaded the expression profiles of circRNA (GSE147442), microRNA (GSE71332) and mRNA (GSE58144) of endometrium tissues of RIF patients and the control group from the Gene Expression Omnibus (GEO) database. Next, we selected the top 10 DEcircRNAs (5 upregulation and 5 downregulation) of GSE147442 as the initial candidates based on t-test. Then we established a circRNA–miRNA–mRNA regulation network depending on circRNA-miRNA interaction and miRNA–mRNA interaction by applying website tools and screening. Also, aim to expound the underlying mechanism of differently expressed mRNA (DEmRNAs) in ceRNA network, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway analysis were performed. Moreover, we chose three down-regulation circRNAs with a high degree of connectivity in ceRNA network as hub-circRNAs. While only hsa_circ_0038383 was identified with differential expression. Furthermore, hub genes of mRNA were screened out by a protein-protein interaction (PPI) network with Cytoscape MCODE (Molecular Complex Detection). Eventually, a circRNA-miRNA–mRNA subnetwork axis was successfully conducted. Thus, this research provided a forceful tool for searching potential biomarkers or therapeutic targets for RIF. A workflow summarizing our procedure is presented in Figure 1.

Flow chart of the ceRNA network analysis.

Figure 1. Flow chart of the ceRNA network analysis.

Results

Identification of DEcircRNAs, differently expressed miRNAs (DEmiRNAs) and DEmRNAs in RIF endometrial tissues

The basic information of three GEO datasets (GSE147442, GSE71332 and GSE58144) used in this study was shown in Table 1. 1468 DEcircRNAs (387 upregulation and 1081 downregulation) were filtered based on predetermined threshold in GSE147442. The DEcircRNAs were visualized by a volcano plot and a heat map (Figure 2A, 2B). The top ten DEcircRNAs in Table 2 were selected as initial candidates.

Table 1. Basic information of the three microarray datasets from GEO.

Data sourcePlatformTissueSample size (RIF/control)AuthorYearRegionRNA type
GSE147442GPL21825Endometrium8/8Yan2020ChinacircRNA
GSE71332GPL18402Endometrium7/5Fan2017ChinamiRNA
GSE58144GPL15789Endometrium43/42van Hooff2015NetherlandsmRNA
Differentially expressed circRNAs in RIF patients compared the control group. (A) Heatmap of the differentially expressed circRNAs in RIF based on GSE147442. (B) Volcano map for all circRNAs in GSE147442.

Figure 2. Differentially expressed circRNAs in RIF patients compared the control group. (A) Heatmap of the differentially expressed circRNAs in RIF based on GSE147442. (B) Volcano map for all circRNAs in GSE147442.

Table 2. Basic characteristics of the top ten DEcircRNAs.

CircRNA IDFCP valueGene symbolRegulation
hsa_circ_00809683.170.0064ADAM22Up
hsa_circ_00007002.890.0026CHD9Up
hsa_circ_00009222.850.0081ZNF536Up
hsa_circ_00261342.460.0018TUBA1CUp
hsa_circ_00958572.420.0021PRDM11Up
hsa_circ_00001150.340.0003CSDE1Down
hsa_circ_00853260.340.0040EIF3EDown
hsa_circ_00754020.360.0015GNB2L1Down
hsa_circ_00910530.360.0035RPS4XDown
hsa_circ_00383830.360.0089THUMPD1Down

After batch effect normalization and further analyses using a limma (versions 3.30.0) R package, we screened out 143 upregulation miRNAs and 58 downregulation miRNAs with p<0.05 from GSE71332 (Figure 3A, 3B). According to the above top ten DEcircRNAs, circRNA–miRNA pair prediction was conducted using online tool CircInteractome and ENCORI database. 217 target miRNAs based on 5 upregulation circRNAs and 331 target miRNAs based on 5 downregulation circRNAs were obtained. Subsequently, the predicted MREs of circRNAs and DEmiRNAs were intersected to acquire 5 down- and 23 upregulation miRNAs using Venn software online (Figure 3C, 3D).

Identification of differentially expressed miRNAs. (A) Heatmap of the differentially expressed miRNAs from the GEO microarray GSE71332. (B) Volcano map for all miRNAs in GSE71332. (C, D) Overlapping between circRNA-related target miRNAs predicted by online tool and DEmiRNAs in GSE71332.

Figure 3. Identification of differentially expressed miRNAs. (A) Heatmap of the differentially expressed miRNAs from the GEO microarray GSE71332. (B) Volcano map for all miRNAs in GSE71332. (C, D) Overlapping between circRNA-related target miRNAs predicted by online tool and DEmiRNAs in GSE71332.

Using the same data processing method, 678 upregulation and 988 upregulation mRNAs in GSE58144 were selected (Figure 4A, 4B). On basis of similar methods, 59 candidate mRNAs were chosen, of which 2 were upregulated and 57 were downregulated (Figure 4C, 4D).

Identification of differentially expressed mRNAs. (A) Heatmap of the differentially expressed mRNAs in GSE58144. (B) Volcano map for all mRNAs in GSE58144. (C, D) Venn diagram analysis of DEmRNA-predicted targets and differentially expressed mRNAs in GEO.

Figure 4. Identification of differentially expressed mRNAs. (A) Heatmap of the differentially expressed mRNAs in GSE58144. (B) Volcano map for all mRNAs in GSE58144. (C, D) Venn diagram analysis of DEmRNA-predicted targets and differentially expressed mRNAs in GEO.

Construction of the ceRNA network

To better comprehend the role of circRNAs in miRNAs mediated mRNAs, a circRNA–miRNA–mRNA regulatory network was generated employing a combination of circRNA–miRNA pairs and miRNA–mRNA pairs after multiple screening steps. The ceRNA contained 33 circRNA–miRNA pairs and 63 miRNA–mRNA pairs, including 10 DEcircRNAs, 28 DEmiRNAs and 59 DEmRNAs. The subnetwork was presented using cytoscape software (Figure 5).

The ceRNA network of circRNA–miRNA–mRNA in RIF. Diamonds represent circRNAs, triangles indicate miRNAs, and rounds indicate mRNAs. Node size represent the degrees of node in the ceRNA network. The node of red and green color express upregulation and downregulation, respectively.

Figure 5. The ceRNA network of circRNA–miRNA–mRNA in RIF. Diamonds represent circRNAs, triangles indicate miRNAs, and rounds indicate mRNAs. Node size represent the degrees of node in the ceRNA network. The node of red and green color express upregulation and downregulation, respectively.

GO and KEGG enrichment analyses

To find out the biological functions and pathway of DEmRNAs in ceRNA network, GO annotation and KEGG pathway analysis was performed. GO enrichment analysis was constituted with biological processes, molecular functions and cellular components. In biological process terms, the DEmRNAs were primarily involved in ‘regulation of transcription’ and ‘gene expression’. For cellular component terms, the DEmRNAs were mainly enrichment in ‘nucleus’, ‘nucleoplasm’, ‘transcription regulator complex’ and ‘cell-cell junction’. While ‘transcription factor activity’, ‘BMP receptor activity’, ‘transmembrane receptor protein serine/threonine kinase activity’ and ‘SMAD binding’ were mainly related to molecular function (Figure 6A). As shown in Figure 6B, KEGG pathway analysis indicated that the 59 aberrant DEmRNAs in RIF were significantly enriched in ‘apoptosis’, ‘fluid shear stress and atherosclerosis’, ‘NOD-like receptor signaling pathway’, ‘Mitophagy’, ‘TGF-beta signaling pathway’ and ‘NF-kappa B signaling pathway’, etc. (Figure 6B).

GO and KEGG analyses of 59 mRNAs. (A) GO analysis. (B) KEGG pathway analysis.

Figure 6. GO and KEGG analyses of 59 mRNAs. (A) GO analysis. (B) KEGG pathway analysis.

Identification and verification of hub DEcircRNAs

CircRNAs act as hub node in biological networks. Three down-regulation circRNAs (hsa_circ_0038383, hsa_circ_0000115 and hsa_circ_0091053) with high degree of connectivity calculated by the cytoHubba plugin of cytoscape were filtered as hub-circRNAs. Their basic characteristics were expressed in Table 3. To further assess the expression of these three selected hub DEcircRNAs, 18 RIF patients and 16 healthy women were included as a validation cohort. The detailed characteristics of these patients were summarized in Table 4. Full-length sequences were got from CircInteractome website, and sequences of specific primers for each circRNA were designed by Primer-BLAST in PubMed. The primers information was shown in Table 5. It was noteworthy that only hsa_circ_0038383 was significantly downregulated in endometrial tissues of RIF patients compared with control groups (Figure 7A). Hsa_circ_0000115 and hsa_circ_0091053 showed no significant difference (Figure 7B, 7C). The basic structural patterns of the three circRNAs were displayed in Figure 7D7F. Consequently, hsa_circ_0038383 was selected as the circRNA target for future research.

Table 3. Basic characteristics of the three selected differently expressed circRNAs.

CircRNA IDPositionStrandGenomic lengthBest transcriptGene symbolRegulation
hsa_circ_0038383chr16:20744988-20749278-4290NM_017736THUMPD1Down
hsa_circ_0000115chr1:115276352-115280693-4341NM_001242892CSDE1Down
hsa_circ_0091053chrX:71492452-71496084-3632NM_001007RPS4XDown

Table 4. Clinical characteristics of women recruited in the present study.

VariablesControl groups (n=16)RIF groups (n=18)p value
Age (years)30.19±2.7930.89±3.430.52a
BMI (kg/m2)21.14±1.5722.01±2.210.20a
FSH (mIU/ml)7.53±2.027.25±1.560.66a
LH (mIU/ml)4.68±2.084.01±2.030.37a
Estradiol (pg/ml)38.60±10.7634.87±11.130.33a
Infertility duration (years)4.16±1.793.64±2.260.43a
Number of embryo transfer cycles1(14/16, 87.5%)
2(2/16, 12.5%)
3(9/18, 50.0%)
4(4/18, 22.2%)
5(3/18, 16.7%)
6(2/18, 11.1%)
<0.001b
Number of embryos per transfer1.50±0.631.61±0.610.61a
aStudent t-test; bFisher’s exact test.
Data presented as mean ± SD. P < 0.05was considered statistically significant.

Table 5. Primer sequences of circRNAs for RT-qPCR.

circRNAsForward primer (5’-3’)Reverse primer (5’-3’)
hsa_circ_0038383TACGGCCAGAAATCGAGCTTATGTGCCTGAGATGGGTAACA
hsa_circ_0000115TGCCTCAAGGAACAGTCATTAATGGGTTTCCCAGTCCGTC
hsa_circ_0091053GAGCAGTGGGTGAAATGGGTATGGACGAGGAGCAAACACA
hsa-miR-196b-5pCGCGCTAGGTAGTTTCCTGTTGTTGGGTailing Reaction, Tm (65° C)
hsa-miR-424-5pCGCGCGCAGCAGCAATTCATGTTTTGTailing Reaction, Tm (65° C)
U6CGCAGAGAAGATTAGCATGGCCCCTGTailing Reaction, Tm (65° C)
HOXA9TAAACCTGAACCGCTGTCGGCCAGTTGGCTGCTGGGTTAT
PBX1GCTGATGCATTCCCATGCTGTGGGCTCCTCGGATACTCAA
HOXA3GCGACCTACTACGACAGCTCCTCACTCAGTTCGTGTGCCT
β-actinCTGGACTTCGAGCAAGAGATGGAGTTGAAGGTAGTTTCGTGGA
Expression verification and structure of circRNAs. The expression of hsa

Figure 7. Expression verification and structure of circRNAs. The expression of hsa_circ_0038383 (A), hsa_circ_0000115 (B) and hsa_circ_0091053 (C) in RIF compared with the controls. Structural patterns of hsa_circ_0011385 these three DEcircRNAs by circRNADb.

PPI network construction and hub mRNA validation

The PPI network was established based on protein interaction provided by STRING database. 30 nodes and 134 edges were included in this PPI network (Figure 8A). Then the most significant cluster in this PPI network was uncovered applying cytoscape MCODE module. The results indicated that only one cluster as key genes, including PBX1, HOXA9 and HOXA3 (Figure 8B). Based on the above results, we constructed a circRNA–miRNA–mRNA subnetwork, consisting of one DEcircRNA (hsa_circ_0038383), two DEmiRNAs (has-miR-196b-5p and has-miR-424-5p), and three DEmRNAs (HOXA3, HOXA9 and PBX1), which provided a novel insight in the pathogenesis and treatment of RIF (Figure 8C). Aimed to confirm whether the expression level of miRNAs and mRNAs in the particular axis is in accordance with prediction, we detected their expression in endometrial tissues of RIF patients compared with the healthy women. The results showed that only has-miR-196b-5p and HOXA9 was consistent the predicted results (Figure 9A9E), which pointed that the hsa_circ_0038383/miR-196b-5p/HOXA9 axis may be a key pathway leading to RIF.

PPI network, hub genes, and circRNA–miRNA–mRNA subnetwork for hsa

Figure 8. PPI network, hub genes, and circRNA–miRNA–mRNA subnetwork for hsa_circ_0038383. (A) PPI network of 59 DEmRNA. (B) Three hub genes extracted from the PPI network based on the MCODE algorithm. (C) The circRNA-miRNA-mRNA axes of hsa_circ_0038383.

Expression verification of miRNAs and mRNAs in the subnetwork. The expression of has-miR-196b-5p (A), has-miR-424-5p (B), HOXA9 (C), PBX1 (D) and HOXA3 (E) in endometrial tissues of RIF and the control group.

Figure 9. Expression verification of miRNAs and mRNAs in the subnetwork. The expression of has-miR-196b-5p (A), has-miR-424-5p (B), HOXA9 (C), PBX1 (D) and HOXA3 (E) in endometrial tissues of RIF and the control group.

Discussion

RIF is a major cause of female infertility and a challenging problem in assisted reproduction field. In addition, the pregnancy rate after the therapy of IVF-ET still remains approximate 25% [22]. Good endometrium receptivity during the window of implantation is necessary for successful embryo transfer. Aberrant gene expression is a key factor in causing insufficient endometrium receptivity, involving in coding and non-coding gene. Wang et al. [23] found several mRNAs which may affect embryo implantation, and identified one possible target in endometrial tissues may be used for treating RIF based on GSE58144. Another study demonstrated that 15 lncRNAs abnormally expressed in endometrial tissues could be regarded as potential biomarkers for RIF through performing gene sequencing and further RT-qPCR verification [24]. Additionally, six key lncRNAs were screened out and their ceRNA subnetworks were conducted by comprehensive bioinformatics analysis, which signified lncRNAs could be developed into predictive biomarkers for RIF [25]. The high stability circRNAs with cyclic structures, exhibit specific expression in tissue or developmental stage and make this type of ncRNAs more valuable as new biomarkers for diverse diseases. Recently, a study showed that several circRNAs were existed with different expression between RIF patients and normal controls, which may furnish diagnosis molecular and therapeutic target for RIF [18].

Nowadays, more and more attention was attracted in investigating the underlying functions and mechanisms of circRNAs through establishing ceRNA network. For instance, in the research of lung adenocarcinoma, Liang et al. [26] appraised significant circRNAs and meaningful circRNA–miRNA–mRNA networks based on the guidance of bioinformatics. Also, hsa_circ_0001955/hsa_circ_0000977-mediated ceRNA sub-network was confirmed to uncovering potential target of colorectal cancer by whole-transcriptome analysis [27]. Furthermore, circUBAP2-mediated circRNA–miRNA–mRNA network modulated the development and progression of pancreatic adenocarcinoma via regulating the infiltration and function of immune cells [28]. In addition, a circRNA-miRNA-hub gene axis was constructed to discover potential therapeutic targets for spinal cord injury by way of microarray data mining and comprehensive bioinformatics analyses [29].

To further understand the potential molecular mechanism and functions of circRNAs on RIF, we constructed a circRNA–miRNA–mRNA regulatory network by differential expression analysis, intersection analysis and correlation analysis. The DEmRNAs involving in biological progress and pathway were performed by GO and KEGG analysis. Although there were few reports on the GO terms directly related to RIF, the GO terms related to endometrial receptivity have been gradually discovered. The GO analysis showed BMP signaling pathway and SMAD binding are two important terms. For example, BMP mediated endometrial receptivity and decidualization [30]. Additionally, smad dependent TGF-β signaling acted a major role in the process of embryo implantation. KEGG pathway analysis indicated that Apoptosis, TGF-beta signaling pathway and signaling pathway. Yu et al. found that the apoptosis of HIF-1α affected embryo implantation during implantation window of women with RIF [31]. Another study indicated that TGF-β1 was down-regulated in patients with RIF and played an important role in endometrial receptivity and embryo implantation [32]. Moreover, Ersahin A et al. found that the expression of endometrial NF-κB expression was disturbed in women with RIF [33]. Therefore, the GO terms or KEGG pathways provided direction for the future mechanism research. Subsequently a circRNA–miRNA–mRNA subnetwork was conducted, of which hsa_circ_0038383 may work as a ceRNA to capture hsa-miR-196b-5p in regulating the expression of HOXA9 and PBX1. Also, hsa_circ_0038383 may positively modulate expression of HOXA3 by acting as hsa-miR-424-5p sponge. These results revealed that the circRNA–miRNA–mRNA network of hub-circRNA not only be a novel candidate biomarker, but also provided evidence for regulation mechanism of circRNA.

Hsa_circ_0038383/miR-196b-5p/HOXA9 and hsa_circ_0038383/miR-196b-5p/PBX1 were two essential axes from the subnetwork. The verification results revealed that only the hsa_circ_0038383/miR-196b-5p/HOXA9 axis may be a key pathway leading to RIF. HOXA9 (Homeobox A9), a homeodomain transcription factor, was a regulator of embryonic development similar to HOX members [34]. Recently, its role in endometrium was gradually reported. The expression of HOXA9 was increased during mid-secretory phase of menstrual cycle and silence of Hoxa9 in the murine uterus reduced average implantation rates [35]. This result prompted that HOXA9 low expression could result in decreasing endometrial receptivity. In our study, the expression of HOXA9 should be downregulated in RIF with low endometrial receptivity, which was positively with hsa_circ_0038383. These consequences suggested that our research was consistent with other study. HOXA9 as the target of hsa-miR-196b-5p, the targeting relationship between them had been confirmed by multiple studies. One study reported that miR-196b directly targeted HOXA9 adjusting aggressiveness through NF-κB activity in non-small cell lung cancer cells [36]. Another study showed that in mixed lineage leukaemia (MLL)-rearranged leukaemia, miR-196b could directly target left HOXA9/MEIS1 oncogenes and FAS tumor suppressor [37]. However, the relationship between miR-196b and PBX1 was still unknown. Therefore, we speculated that hsa_circ_0038383 influenced uterine receptivity and embryo implantation by miR-196b/HOXA9 axis, and we would test this axis in the future.

Hsa_circ_0038383/hsa-miR-424-5p/HOXA3 was also included in the subnetwork. The expression of miRNAs and target genes between menstrual endometria and early pregnancy were compared, of which miR-424-5p were significantly downregulated during early pregnancy deciduas [38]. But related research about HOXA3 in endometrial receptivity and the target relationship between HOXA3 and miR-424-5p were still unknown. Moreover, the correlation of genes expression in hsa_circ_0038383/hsa-miR-424-5p/HOXA3 axis was not confirmed in our research. Whether this axis played roles in endometrium acceptability needed more experiments.

In this study, a novel circRNA–miRNA–mRNA network related with RIF was proposed and hsa_circ_0038383 imbalance was verified in RIF. This network would contribute to explore the initiation and progression of RIF and further develop potential treatment strategies for this disease. However, in view of the results are mainly based on computational biology and RT-qPCR, further biological and molecular experiments are indispensable to verify our hypothesis.

Conclusions

In summary, we constructed a potential hsa_circ_0038383-mediated ceRNA subnetwork which provided new insight into exploring molecular mechanism and offering a candidate biomarker for RIF. Especially, hsa_circ_0038383/miR-196b-5p/HOXA9 axis may be a key pathway in affecting uterine receptivity and embryo implantation, and further in-depth molecular biology experiments on this axis are necessary to verify the circRNA role in RIF. This study contributed to seeking new ideas for diagnostic biomarkers and therapeutic targets for patients with RIF.

Materials and Methods

Microarray data information

The raw data of circRNA, microRNA and mRNA were obtained from the GEO http://www.ncbi.nlm.nih.gov/geo/) dataset, which is an international public repository and recording platform employed for searching any appropriate datasets. GSE147442 is a circRNAs database of RIF and based on the GPL21825 platform, containing 8 RIF endometrial and 8 normal endometrium specimens. GSE71332 is a miRNAs database of RIF and based on the GPL18402 platform, consisting of 7 RIF endometrial and 5 normal endometrium specimens. GSE58144 is a mRNAs database of RIF and based on the GPL15789 platform, involving in 43 RIF endometrial and 42 normal endometrium samples. The fundamental information of this gene chip was shown in Table 1. DEcircRNAs were identified based on the fold-change and Student’s t testing (p<0.05, FC>1.5). While data in GSE71332 and GSE58144 were normalized and processed using bioconductor limma (versions 3.30.0) R package. DEmiRNA and DEmRNA were acknowledged with p value <0.05. The basic information of these three GEO datasets used in this study was shown in Table 1.

Prediction of circRNA–miRNA pairs

The circBase website (http://www.circbase.org/) was utilized to observing the basic information of circRNAs [39]. The target miRNAs and miRNA response elements (MREs) presented in circRNAs were predicted by The Circular RNA Interactome (CircInteractome, https://circinteractome.nia.nih.gov/) [40] and the Encyclopedia of RNA Interactomes (ENCORI database, http://starbase.sysu.edu.cn/index.php) [41]. The overlapping miRNAs were further screened as potential target miRNAs to the DEcircRNAs based on the corresponding to target miRNAs to DEcircRNAs and DEmiRNAs of GSE71332.

Prediction of miRNA–mRNA pairs

The miRNA–mRNA interaction were predicted with three online websites respectively, including TargetScan (http://www.targetscan.org/) [42], miRDB (http://www.mirdb.org/) [43] and miRTarBase [44]. Only target mRNAs recognized consistently by these three databases were selected as candidate targets and further were used to intersect with differentially expressed mRNAs of GSE58144 to screen out potential target mRNAs for the miRNAs.

Construction of the circRNA–miRNA–mRNA network

To better uncover the correspondence among circRNAs, miRNAs and mRNAs, we constructed a circRNA–miRNA–mRNA regulatory network utilizing a combination of circRNA-miRNA pairs and miRNA–mRNA pairs. And the ceRNA network was visualized using cytoscape (http://cytoscape.org/; version 3.7.1) software [45].

GO term and KEGG pathway enrichment analysis

The definitions of GO and KEGG analysis was described in our previous study [46]. To evaluate the main function pathways of RIF, DEmRNAs in ceRNA network were performed through GO annotation and KEGG pathway analyses. The results were depicted using clusterProfiler (versions 3.18.0) package in R (R-3.6.0).

Construction of PPI network and identification of hub genes

The PPI network was established for indicating the interaction among the determined DEmRNAs more intuitively, according to the Search Tool for the Tetrieval of Interacting Genes (STRING) database (http://string-db.org/, version 11.0) [47]. And then the result was shown using cytoscape software. Subsequently, the hub mRNAs were appraised by cytoscape MCODE module, which is a plug-in used to find closely connected nodes in a complex network based on topology.

Endometrial tissues

34 women (aged 24-40) treated at Reproductive Medicine Center of Yantai Yuhuangding Hospital were enrolled in our study. All of them signed informed consent forms approved by the Institutional Review Board of Yantai Yuhuangding Hospital (reference number 2019-121). The endometrial tissue samples in mid-secretory phase were obtained from RIF patients (n = 18) and control groups (n = 16), respectively. The women still did not get pregnancy after at least three IVF–ET failure cycles were assigned as RIF patients. The control women who experienced IVF–ET cycle due to tubal obstruction without hydrosalpinx achieved a clinical pregnancy after their first or second embryo transfer. Endometrium-related diseases, hydrosalpinx, polycystic ovarian syndrome (PCOS), etc. were excluded. All participants hold normal hormone level, normal endometrial thickness and morphology and a regular menstrual cycle (28-31 days).

RNA extraction and RT-qPCR

Total RNA was isolated from the tissue samples using Trizol reagent following manufacturer’s instructions (Sparkjade, Qingdao, China). Then 1 μg RNA of each sample was reverse-transcribed to obtain cDNA using SPARKscript II RT Plus Kit (With gDNA Eraser) (Sparkjade, Qingdao, China). The expression of circRNAs and mRNAs in these individual samples was performed by qRT-PCR reaction using SYBR Green qPCR Mix kit (With ROX) (Sparkjade, Qingdao, China) following: 94° C for 2 min, followed by 40 cycles of 95° C for 10 s and 60° C for 30 s. The extraction and amplification of miRNAs was performed according to the manufacturer’s instructions (SPARKeasy animal tissue/Cell microRNA kit, miRNA first strand synthesis kit and miRNA SYBR Green qPCR Mix kit, Sparkjade, Qingdao, China). All RT-qPCR were repeated three times. The relative expression of all genes was calculated using the 2-ΔΔCT method.

Statistical analysis

The continue variables were expressed with means ± standard deviation. The statistical analysis of circRNAs was performed using GraphPad Prism (GraphPad Software, Inc., La Jolla, CA). The significance between groups was tested by Student’s t test or Fisher’s exact test. When p<0.05 were considered as statistically significant.

Abbreviations

GEO: Gene Expression Omnibus; RIF: recurrent implantation failure; IVF–ET: in vitro fertilization-embryo transfer; ART: assisted reproductive technology; circRNAs: circular RNAs; DEcircRNAs: differently expressed circRNAs; DEmiRNAs: differently expressed miRNAs; DEmRNAs: differently expressed mRNAs; ceRNA: competing endogenous RNAs; ncRNAs: non-coding RNAs; RT-qPCR: real-time quantitative polymerase chain reaction; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Gene and Genome; STRING: the Search Tool for the Tetrieval of Interacting Genes database; PPI: protein-protein interaction; MCODE: Molecular Complex Detection; RBPs: RNA-binding proteins.

Author Contributions

Huishan Zhao: conception, idea, design of the study, experiments and writing the original draft; Lili Chen and Yinghua Shan: data collection and curation, experiments, review and editing the draft; Gang Chen: methodology, data analysis and visualization; Yongli Chu and Huangguan Dai: clinical case collection and clinical data analysis; Xuemei Liu and Hongchu Bao: supervision, project administration, funding acquisition and paper finalization. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Funding

This work was supported by National Natural Science Foundation of China (grant no. 81601276) and the Science and Technology Plan of Yantai (grant no. 2020MSGY089 and 2019YD020) and the Natural Science Foundation of Shandong Province (grant no. ZR2017PH047).

References

  • 1. Kupka MS, D'Hooghe T, Ferraretti AP, de Mouzon J, Erb K, Castilla JA, Calhaz-Jorge C, De Geyter Ch, Goossens V, and European IVF-Monitoring Consortium (EIM), and European Society of Human Reproduction and Embryology (ESHRE). Assisted reproductive technology in Europe, 2011: results generated from European registers by ESHRE. Hum Reprod. 2016; 31:233–48. https://doi.org/10.1093/humrep/dev319 [PubMed]
  • 2. Polanski LT, Baumgarten MN, Quenby S, Brosens J, Campbell BK, Raine-Fenning NJ. What exactly do we mean by ’recurrent implantation failure’? a systematic review and opinion. Reprod Biomed Online. 2014; 28:409–23. https://doi.org/10.1016/j.rbmo.2013.12.006 [PubMed]
  • 3. Margalioth EJ, Ben-Chetrit A, Gal M, Eldar-Geva T. Investigation and treatment of repeated implantation failure following IVF-ET. Hum Reprod. 2006; 21:3036–43. https://doi.org/10.1093/humrep/del305 [PubMed]
  • 4. Coughlan C, Ledger W, Wang Q, Liu F, Demirol A, Gurgan T, Cutting R, Ong K, Sallam H, Li TC. Recurrent implantation failure: definition and management. Reprod Biomed Online. 2014; 28:14–38. https://doi.org/10.1016/j.rbmo.2013.08.011 [PubMed]
  • 5. Kliman HJ, Frankfurter D. Clinical approach to recurrent implantation failure: evidence-based evaluation of the endometrium. Fertil Steril. 2019; 111:618–28. https://doi.org/10.1016/j.fertnstert.2019.02.011 [PubMed]
  • 6. Davari-Tanha F, Shahrokh Tehraninejad E, Ghazi M, Shahraki Z. The role of G-CSF in recurrent implantation failure: a randomized double blind placebo control trial. Int J Reprod Biomed. 2016; 14:737–42. [PubMed]
  • 7. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013; 495:333–38. https://doi.org/10.1038/nature11928 [PubMed]
  • 8. Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015; 365:141–48. https://doi.org/10.1016/j.canlet.2015.06.003 [PubMed]
  • 9. Kristensen LS, Andersen MS, Stagsted LV, Ebbesen KK, Hansen TB, Kjems J. The biogenesis, biology and characterization of circular RNAs. Nat Rev Genet. 2019; 20:675–91. https://doi.org/10.1038/s41576-019-0158-7 [PubMed]
  • 10. Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, Hanan M, Behm M, Bartok O, Ashwal-Fluss R, Herzog M, Schreyer L, Papavasileiou P, et al. Circular RNAs in the mammalian brain are highly abundant, conserved, and dynamically expressed. Mol Cell. 2015; 58:870–85. https://doi.org/10.1016/j.molcel.2015.03.027 [PubMed]
  • 11. Li X, Yang L, Chen LL. The biogenesis, functions, and challenges of circular RNAs. Mol Cell. 2018; 71:428–42. https://doi.org/10.1016/j.molcel.2018.06.034 [PubMed]
  • 12. Beermann J, Piccoli MT, Viereck J, Thum T. Non-coding RNAs in development and disease: background, mechanisms, and therapeutic approaches. Physiol Rev. 2016; 96:1297–325. https://doi.org/10.1152/physrev.00041.2015 [PubMed]
  • 13. Zhou Z, Sun B, Huang S, Zhao L. Roles of circular RNAs in immune regulation and autoimmune diseases. Cell Death Dis. 2019; 10:503. https://doi.org/10.1038/s41419-019-1744-5 [PubMed]
  • 14. An T, Zhang J, Ma Y, Lian J, Wu YX, Lv BH, Ma MH, Meng JH, Zhou YT, Zhang ZY, Liu Q, Gao SH, Jiang GJ. Relationships of non-coding RNA with diabetes and depression. Sci Rep. 2019; 9:10707. https://doi.org/10.1038/s41598-019-47077-9 [PubMed]
  • 15. Koot YE, van Hooff SR, Boomsma CM, van Leenen D, Groot Koerkamp MJ, Goddijn M, Eijkemans MJ, Fauser BC, Holstege FC, Macklon NS. An endometrial gene expression signature accurately predicts recurrent implantation failure after IVF. Sci Rep. 2016; 6:19411. https://doi.org/10.1038/srep19411 [PubMed]
  • 16. Shi C, Han HJ, Fan LJ, Guan J, Zheng XB, Chen X, Liang R, Zhang XW, Sun KK, Cui QH, Shen H. Diverse endometrial mRNA signatures during the window of implantation in patients with repeated implantation failure. Hum Fertil (Camb). 2018; 21:183–94. https://doi.org/10.1080/14647273.2017.1324180 [PubMed]
  • 17. Shi C, Shen H, Fan LJ, Guan J, Zheng XB, Chen X, Liang R, Zhang XW, Cui QH, Sun KK, Zhao ZR, Han HJ. Endometrial MicroRNA signature during the window of implantation changed in patients with repeated implantation failure. Chin Med J (Engl). 2017; 130:566–73. https://doi.org/10.4103/0366-6999.200550 [PubMed]
  • 18. Liu L, Li L, Ma X, Yue F, Wang Y, Wang L, Jin P, Zhang X. Altered circular RNA expression in patients with repeated implantation failure. Cell Physiol Biochem. 2017; 44:303–13. https://doi.org/10.1159/000484887 [PubMed]
  • 19. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011; 146:353–58. https://doi.org/10.1016/j.cell.2011.07.014 [PubMed]
  • 20. Lu HC, Yao JQ, Yang X, Han J, Wang JZ, Xu K, Zhou R, Yu H, Lv Q, Gu M. Identification of a potentially functional circRNA-miRNA-mRNA regulatory network for investigating pathogenesis and providing possible biomarkers of bladder cancer. Cancer Cell Int. 2020; 20:31. https://doi.org/10.1186/s12935-020-1108-3 [PubMed]
  • 21. Zhang X, Zhang Q, Wu Q, Tang H, Ye L, Zhang Q, Hua D, Zhang Y, Li F. Integrated analyses reveal hsa_circ_0028883 as a diagnostic biomarker in active tuberculosis. Infect Genet Evol. 2020; 83:104323. https://doi.org/10.1016/j.meegid.2020.104323 [PubMed]
  • 22. Cavalcante MB, Cavalcante CT, Sarno M, Barini R. Intrauterine perfusion immunotherapies in recurrent implantation failures: systematic review. Am J Reprod Immunol. 2020; 83:e13242. https://doi.org/10.1111/aji.13242 [PubMed]
  • 23. Wang F, Liu Y. Identification of key genes, regulatory factors, and drug target genes of recurrent implantation failure (RIF). Gynecol Endocrinol. 2020; 36:448–55. https://doi.org/10.1080/09513590.2019.1680622 [PubMed]
  • 24. Xu H, Zhou M, Cao Y, Zhang D, Han M, Gao X, Xu B, Zhang A. Genome-wide analysis of long noncoding RNAs, microRNAs, and mRNAs forming a competing endogenous RNA network in repeated implantation failure. Gene. 2019; 720:144056. https://doi.org/10.1016/j.gene.2019.144056 [PubMed]
  • 25. Feng C, Shen JM, Lv PP, Jin M, Wang LQ, Rao JP, Feng L. Construction of implantation failure related lncRNA-mRNA network and identification of lncRNA biomarkers for predicting endometrial receptivity. Int J Biol Sci. 2018; 14:1361–77. https://doi.org/10.7150/ijbs.25081 [PubMed]
  • 26. Liang L, Zhang L, Zhang J, Bai S, Fu H. Identification of circRNA-miRNA-mRNA networks for exploring the fundamental mechanism in lung adenocarcinoma. Onco Targets Ther. 2020; 13:2945–55. https://doi.org/10.2147/OTT.S235664 [PubMed]
  • 27. Ding B, Yao M, Fan W, Lou W. Whole-transcriptome analysis reveals a potential hsa_circ_0001955/hsa_circ_0000977-mediated miRNA-mRNA regulatory sub-network in colorectal cancer. Aging (Albany NY). 2020; 12:5259–79. https://doi.org/10.18632/aging.102945 [PubMed]
  • 28. Zhao R, Ni J, Lu S, Jiang S, You L, Liu H, Shou J, Zhai C, Zhang W, Shao S, Yang X, Pan H, Han W. CircUBAP2-mediated competing endogenous RNA network modulates tumorigenesis in pancreatic adenocarcinoma. Aging (Albany NY). 2019; 11:8484–501. https://doi.org/10.18632/aging.102334 [PubMed]
  • 29. Peng P, Zhang B, Huang J, Xing C, Liu W, Sun C, Guo W, Yao S, Ruan W, Ning G, Kong X, Feng S. Identification of a circRNA-miRNA-mRNA network to explore the effects of circRNAs on pathogenesis and treatment of spinal cord injury. Life Sci. 2020; 257:118039. https://doi.org/10.1016/j.lfs.2020.118039 [PubMed]
  • 30. Doherty LF, Taylor HS. Leiomyoma-derived transforming growth factor-β impairs bone morphogenetic protein-2-mediated endometrial receptivity. Fertil Steril. 2015; 103:845–52. https://doi.org/10.1016/j.fertnstert.2014.12.099 [PubMed]
  • 31. Yu X, Gao C, Dai C, Yang F, Deng X. Endometrial injury increases expression of hypoxia-inducible factor and angiogenesis in the endometrium of women with recurrent implantation failure. Reprod Biomed Online. 2019; 38:761–67. https://doi.org/10.1016/j.rbmo.2018.12.027 [PubMed]
  • 32. Guo F, Si C, Zhou M, Wang J, Zhang D, Leung PCK, Xu B, Zhang A. Decreased PECAM1-mediated TGF-β1 expression in the mid-secretory endometrium in women with recurrent implantation failure. Hum Reprod. 2018; 33:832–43. https://doi.org/10.1093/humrep/dey022 [PubMed]
  • 33. Ersahin A, Acet M, Acet T, Yavuz Y. Disturbed endometrial NF-κB expression in women with recurrent implantation failure. Eur Rev Med Pharmacol Sci. 2016; 20:5037–40. [PubMed]
  • 34. Domingo-Reines J, López-Ornelas A, Montes R, Romero T, Rodriguez-Llamas JL, Lara-Rodarte R, González-Pozas F, Ayllon V, Menendez P, Velasco I, Ramos-Mejia V. Hoxa9 and EGFP reporter expression in human embryonic stem cells (hESC) as useful tools for studying human development. Stem Cell Res. 2017; 25:286–90. https://doi.org/10.1016/j.scr.2017.08.004 [PubMed]
  • 35. Xu B, Geerts D, Bu Z, Ai J, Jin L, Li Y, Zhang H, Zhu G. Regulation of endometrial receptivity by the highly expressed HOXA9, HOXA11 and HOXD10 HOX-class homeobox genes. Hum Reprod. 2014; 29:781–90. https://doi.org/10.1093/humrep/deu004 [PubMed]
  • 36. Yu SL, Lee DC, Sohn HA, Lee SY, Jeon HS, Lee JH, Park CG, Lee HY, Yeom YI, Son JW, Yoon YS, Kang J. Homeobox A9 directly targeted by miR-196b regulates aggressiveness through nuclear factor-kappa B activity in non-small cell lung cancer cells. Mol Carcinog. 2016; 55:1915–26. https://doi.org/10.1002/mc.22439 [PubMed]
  • 37. Li Z, Huang H, Chen P, He M, Li Y, Arnovitz S, Jiang X, He C, Hyjek E, Zhang J, Zhang Z, Elkahloun A, Cao D, et al. miR-196b directly targets both HOXA9/MEIS1 oncogenes and FAS tumour suppressor in MLL-rearranged leukaemia. Nat Commun. 2012; 3:688. https://doi.org/10.1038/ncomms1681 [PubMed]
  • 38. Muqri F, Helkin A, Maier KG, Gahtan V. Thrombospondin-5 and fluvastatin promote angiogenesis and are protective against endothelial cell apoptosis. J Cell Biochem. 2020; 121:4154–65. https://doi.org/10.1002/jcb.29686 [PubMed]
  • 39. Glažar P, Papavasileiou P, Rajewsky N. circBase: a database for circular RNAs. RNA. 2014; 20:1666–70. https://doi.org/10.1261/rna.043687.113 [PubMed]
  • 40. Dudekula DB, Panda AC, Grammatikakis I, De S, Abdelmohsen K, Gorospe M. CircInteractome: a web tool for exploring circular RNAs and their interacting proteins and microRNAs. RNA Biol. 2016; 13:34–42. https://doi.org/10.1080/15476286.2015.1128065 [PubMed]
  • 41. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-seq data. Nucleic Acids Res. 2014; 42:D92–97. https://doi.org/10.1093/nar/gkt1248 [PubMed]
  • 42. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. Elife. 2015; 4:e05005. https://doi.org/10.7554/eLife.05005 [PubMed]
  • 43. Chen Y, Wang X. miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res. 2020; 48:D127–31. https://doi.org/10.1093/nar/gkz757 [PubMed]
  • 44. Chou CH, Shrestha S, Yang CD, Chang NW, Lin YL, Liao KW, Huang WC, Sun TH, Tu SJ, Lee WH, Chiew MY, Tai CS, Wei TY, et al. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018; 46:D296–302. https://doi.org/10.1093/nar/gkx1067 [PubMed]
  • 45. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
  • 46. Zhao H, Jiang A, Yu M, Bao H. Identification of biomarkers correlated with diagnosis and prognosis of endometrial cancer using bioinformatics analysis. J Cell Biochem. 2020; 121:4908–21. https://doi.org/10.1002/jcb.29819 [PubMed]
  • 47. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 47:D607–13. https://doi.org/10.1093/nar/gky1131 [PubMed]