Expression Profile Analysis of Circular RNAs in Ovarian Endometriosis by Microarray and Bioinformatics.

BACKGROUND Endometriosis is a common gynecologic disorder with enigmatic etiopathogenesis and is characterized by tumor-like biological behaviors. Recently, circular RNAs (circRNAs) have attracted considerable attention because they exert very important functions in the progression of human cancers. However, little is known about the functions and molecular mechanism of circRNAs in endometriosis. MATERIAL AND METHODS A total of 20 patients with ovarian endometriosis and 4 normal endometrium from women free of endometriosis were included in this study. Ectopic endometrium tissues and paired eutopic endometrium tissues were collected from ovarian endometriosis patients. We assessed the expression profiles of circRNAs in endometriosis by microarray analysis. Expression of selected circRNAs in those tissues was detected by quantitative real-time PCR (qRT-PCR). Based on the target prediction, we constructed a circRNA-miRNA-mRNA competing endogenous RNA (ceRNA) network and elucidated circRNAs through Kyoto Encyclopedia of Genes and Genomes and Gene Ontology analyses. RESULTS We detected 2237 circRNAs, differentially expressed among 3 groups, and then found 8 circRNAs that may be involved in the epithelial-mesenchymal transition process. The qRT-PCR validation suggested that circ_103470 and circ_101102 matched the microarray results. The functional analysis revealed 17 pathways, such as the mTOR signaling pathway, the Hippo signaling pathway, and the HIF-1 signaling pathway, which may be associated with the pathogenesis and development of endometriosis. CONCLUSIONS In general, our results suggest that 2 downregulated circRNAs (circ_103470 and circ_101102) may regulated epithelial-mesenchymal transition in endometriosis via miR-141-5p, which may be a promising therapeutic target in the future.

Background

Endometriosis is classified by the presence of functional ectopic endometrial glands and stroma beyond the uterine cavity at several sites, primarily within the ovaries [1]. It is a leading cause of disability that affects 6–10% of women of reproductive age, which is approximately 176 million women worldwide, and for women who are enduring pelvic pain and/or infertility it may affect 35–50% [2]. However, a diagnosis and treatment delay of 5 to 10 years is common, due to non-specific symptomatology and a lack of noninvasive methods for detecting endometriosis [3]. Taking into account both direct and indirect costs, endometriosis severely impairs the quality of life for women and imposes a considerable national economic burden [4].

Endometriosis is a heritable, complex chronic disease, for which many of the causal pathogenic mechanisms remain unknown. Based on twin studies, the heritability of the condition is estimated at approximately 50% [5], thus emphasizing the importance of genetic or epigenetic contribution to the disease etiology and pathogenesis. Although it is not a malignant condition, its biological behaviors, including metastasis and implantation, are similar to that of tumors [6,7]. There exists a large amount of molecular aberrations among ectopic endometrium (EC) (endometriotic lesions), paired eutopic endometrium (EU), and normal endometrium from women free of endometriosis, which might explain the mechanism of the abnormal growth of EU outside the uterus with endometriosis. From accumulating studies, it was found that epithelial-mesenchymal transition of endometrial cells may be a prerequisite for the establishment of endometriotic lesions [8]. During epithelial-mesenchymal transition, epithelial cells lose their cell polarity and cell-cell adhesion when gaining migratory and invasive properties. Thus, it is necessary to identify potential molecular biomarkers for diagnosis and regulatory factors, underlying the progress of endometriosis, especially the cases that may be induced by epithelial-mesenchymal transition.

Circular RNAs (circRNAs) belong to non-coding RNAs (ncRNAs) family that regulate transcriptional and post-transcriptional gene expression [9]. Unlike traditional linear RNAs, their expression is highly cell-type-specific, tissue-specific, and developmental-specific. Increasing evidence demonstrates that intracellular circRNAs with competing endogenous RNAs (ceRNAs) may act as miRNA sponges that sequester miRNAs by binding miRNAs with miRNA response elements [10]. A study by Conn et al. showed that hundreds of circRNAs were regulated during human epithelial-mesenchymal transition, and these findings demonstrated that circRNAs are both purposefully synthesized and regulated by cell-type specific mechanisms [11]. Therefore, circRNAs are consider important biological regulators for understanding the molecular mechanisms of disease and identifying effective diagnostic biomarkers or therapeutic targets.

Recent studies have shown and emphasized the importance of circRNAs in regulating cancer-related signaling pathways [12]. To our knowledge, however, little is known about the expression and function of circRNAs in endometriosis. In our study, the expression profiles of circRNAs in endometriosis (EC and EU) and normal endometrium from women free of endometriosis were analyzed, with a focus on circRNAs that play a big part in endometriosis, especially that associated with the epithelial-mesenchymal transition process. In particular, using public databases, the clinical significance of circRNAs was identified, and their expression in endometriosis patients was validated by quantitative real-time PCR (qRT-PCR).

Material and Methods

This study was approved by institutional ethics review board of the Second Xiangya Hospital, Central South University (#2016243). All tissue samples were obtained with full and informed patient consent, and all experiments were followed the ethical principles outlined by the 1964 Helsinki Declaration and its later amendments or comparable ethical standards.

Specimens

In this study, snap-frozen cyst walls of ovarian endometriosis and matched EU of the uterus, from the same patient who underwent laparoscopy, were simultaneously collected for qRT-PCR analysis (26 to 36 years of age, proliferative phase, n=20, revised classification of American Fertility Society, r-AFS III–IV). Paired specimens of 4 women, selected randomly, and 4 normal endometrium from women free of endometriosis (control group) (27 to 39 years of age, proliferative phase) were obtained for microarray analysis (Table 1). All the patients had regular menstrual cycles and did not receive hormonotherapy from 6 months before the specimen collection. Furthermore, patients in the control group had no evidence of tumor in the endometrium and did not have histological diagnosis of adenomyosis. The specimen diagnosis was determined by histopathology and the menstrual cycle was confirmed by both last menstrual period and histological ascertain.

Table 1

Detailed information of patients in microarray analysis.

Group	Sample ID	Diagnosis	G/P	rAFS score/stage
Test1(EC)	N1-B-160808	O-EMs (Bilateral)	G0P0	56/IV
Test2(EU)	N1-G-160808
Test1(EC)	N1-B-160901	O-EMs (Bilateral)	G4P1	80/IV
Test2(EU)	N1-G-160901
Test1(EC)	Z1-B-160907	O-EMs (left)	G0P0	24/III
Test2(EU)	Z1-G-160907
Test1(EC)	J1-B-161221	O-EMs (right)	G1P0	24/III
Test2(EU)	J1-G-161221
Control	W1-DZG-161025	CIN III	G3P2
Control	W2-DZG-161107	CIN II	G2P2
Control	W1-DZG-161110	CIN III	G6P1
Control	J2-DZG-161216	CIN III	G3P1

G/P – gestation/production; O-EMs – ovarian endometriosis; rAFS – revised American Fertility Society stage: I: 1–5 score; II: 6–15 score; III: 16–40 score; IV: >40 score. CIN – cervical intraepithelial neoplasias.

Total RNA isolation and quality control

Total RNA was extracted from frozen tissues of 3 groups using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s instructions. The RNA from each sample was quantified using the NanoDrop ND-1000 (Nano Drop Thermo, DE, USA), and the RNA integrity was verified by using the 1% formaldehyde denaturing gel electrophoresis.

CircRNA microarray analysis

The labeling and hybridization of all specimen was conducted following the manufacturer’s protocol (Arraystar, MD, USA). Briefly, total RNAs were digested with Rnase R (Epicentre, Inc.) to remove linear RNAs and enrich circular RNAs. Then, the enriched circular RNAs were amplified and transcribed into fluorescent cRNA utilizing a random priming method (Arraystar Super RNA Labeling Kit; Arraystar). The labeled cRNAs were hybridized onto the Arraystar Human circRNA Array V2 (8×15K, Arraystar), which contains 15 000 probes for 13 617 human circRNAs [10].

After the arrays were washed and scanned in the Agilent Scanner G2505C, the subsequent data were exported by using Agilent Feature Extraction software (version 11.0.1.1) and quantile normalization were performed by using the R software limma package. Differentially expressed circRNAs with statistical significance between 2 groups were identified through Volcano Plot filtering and Fold Change (FC) filtering, and subjected to statistical analyses of type, chromosome distribution, and length distribution. The confirmation of the distinguishable circRNAs expression pattern among samples were performed by Hierarchical Clustering.

Construction of circRNA/miRNA interaction network

Raw data were analyzed by the Student t-tests and the false discovery rate method. The circRNAs with expression change >2 folds and P values <0.05 were considered as statistically significant differential circRNAs. The parental genes of differentially expressed circRNAs were obtained from circBASE database (). To further elucidate the interactions between circRNAs and miRNAs, 2 types of software, Miranda (v3.3a) and TargetScan (Version: 7.0), were used to predict the targets miRNA of each circRNA. In cases of abnormally high predicted miRNA binding site numbers, a circRNA could be considered as a super-sponge regulating gene expression for miRNAs. The top 5 selected miRNA targets for each circRNA were used to describes the circRNA/miRNA interactions, and narrow the scope of circRNA in subsequent research. The circRNA-miRNA networks were visualized by Cytoscape (version 3.6.1).

Validation for the differential expression level of circRNAs

According to the analysis of the raw data and bioinformatic prediction, 3 circRNAs (circ_102362, circ_103470, and circ_101102) aberrantly expressed in EC-EU were selected for further research.

Beta-actin was used as a reference. Total RNA was isolated by using TRIzol reagent (Invitrogen) from 20 paired EC and EU, then reversed into cDNA by using superscript TM III reverse transcriptase (Invitrogen). We evaluated the expression level of the circRNAs by qPCR using SYBR Green assay (Arraystar, MD, USA). Divergent primers were designed and optimized for circRNAs transcripts were as follows: circ_102362 forward: 5′ CTGAATGACATCGAGAAGGTCCT 3′ reverse: 5′ CAGAGAGTTGTCTTTGCTGCCAT 3′; circ_103470 forward: 5′CTTA ACACGCTGGCCCATT3′, reverse: 5′CCAGGTCCCGAAGGATGTAG3′. circ_101102 forward: 5′TCTAAGCATTGAACCAGTGGGATT3′, reverse: 5′TCCTGTCCATTGTGGGCTCAT3′. The reaction program was setting as follows: 50°C, 120 seconds; 95°C, 600 seconds; (95°C, 10 seconds; 60°C, 60 seconds) 40 amplification cycles. All targets and reference were run in triplicate wells. The relative amounts of the selected circRNAs were calculated by double standard curve method, and were shown as the mean ±SEM. The comparisons were analyzed by the Student t tests, and P value of less than 0.05 were considered statistically significant.

Delineation of circRNA-miRNA-mRNA associations and pathway analysis

As demonstrated by lots of researches previously, circRNAs were functioning as miRNA “sponges” and influence the activities of miRNAs in regulating gene expressions. A circRNA-miRNA-mRNA network, based on the 2 verified circRNAs, was constructed to investigate their interactions by using Arraystar target prediction software. The verified circRNAs, the predicted miRNAs and the overlapped mRNAs were included in the network.

Based on the overlapped mRNAs reflected in the circRNA-miRNA-mRNA network, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses were carried out to annotate the gene function and pathways, respectively. Two analyses denote the significance of the pathway correlated to the target genes which mean the selected circRNAs are likely to influence these pathways in the pathological process of endometriosis.

Results

CircRNA microarray expression profiling in endometriosis

Our circRNA profiles show that 2237 circRNAs were differentially expressed among 3 groups (Figure 1A). According to the circRNA microarray results, 1056 circRNAs were differentially expressed between EC and EU, of which 528 circRNAs were upregulated and another 528 circRNAs were downregulated. 1465 circRNAs were differentially expressed in the EC group and the control group, including 842 circRNAs that were upregulated and 623 circRNAs were downregulated. What’s more, 372 circRNAs were found differentially expressed between the EU group and the control group in gross, of which 193 circRNAs of them were upregulated and 179 circRNAs of them were downregulated. The types and the chromosome distribution of these circRNAs are presented in Figure 1B and 1C.

Figure 1

(A) Hierarchical clustering of differentially expressed circular RNAs. (B) Crocs plot showing differently expressed circRNAs on human chromosomes in endometriosis. (C) The typical distribution of differentially expressed circRNAs.

Furthermore, we screened the target circRNAs as the protocol that the P-value <0.05 and FC ≥2 are satisfied at the same time. According to the 2 terms, between EC and EU, 41 circRNAs were upregulated and 107 circRNAs were downregulated (Figure 2A–2D). Between the EC group and the control group, 98 circRNAs were upregulated and 224 circRNAs were downregulated (Figure 2E–2H). Between the EU group and the control group, 49 circRNAs were upregulated and no circRNAs were downregulated (Figure 2I–2L).

Figure 2

(A) Hierarchical clustering of differentially expressed circular RNAs between ectopic endometrium (EC) and eutopic endometrium (EU). The color scale reflects the log2 signal intensity and runs from green (low intensity) to red (strong intensity). (B) Scatter plots show the difference in the expression of circRNAs between EC and EU. (C) Volcano plots illustrated differential expression between 2 different conditions. (D) The cluster shaped bar chart demonstrates the distribution of upregulated and downregulated circRNAs between EC and EU. (E) Hierarchical clustering of differentially expressed circular RNAs between EC and Control. (F) Scatter plots show the difference in the expression of circRNAs between EC and Control. (G) Volcano plots illustrated differential expression between 2 different conditions. (H) The cluster shaped bar chart demonstrates the distribution of upregulated and downregulated circRNAs between EC and Control. (I) Hierarchical clustering of differentially expressed circular RNAs between EU and Control. (J) Scatter plots show the difference in the expression of circRNAs between EU and Control. (K) Volcano plots illustrated differential expression between 2 different conditions. (L) The cluster shaped bar chart demonstrates the distribution of upregulated and downregulated circRNAs between EU and Control. (A, E, I) Each column represents a sample and each row represents a circRNA. The color scale reflects the log2 signal intensity and runs from green (low intensity) to red (strong intensity). (B, F, J) The outermost layer of the circos plot is a chromosome map of the human genome. The increased or decreased circRNAs have been marked in red or green bars, respectively. The second outermost circle represents all target circRNAs detected by sequencing, and the larger inner circle indicates the significantly by whole transcriptome differentially upregulated circRNAs with FC ≥2.0, P<0.05. The smaller inner circle indicates the significantly differentially downregulated circRNAs with FC ≥2.0, P<0.05. (C, G, K) The green vertical lines represent 2.0-FC (log2 scaled) up and down, respectively, and the horizontal green line manifests a P-value of 0.05 (−log10 scaled). The red points in plot denote the differentially expressed circRNAs with statistical significance.

In order to elucidate the circRNAs for further research, we predicted the target miRNAs of these differentially expressed circRNAs. According to the analysis, we selected 8 differentially expressed circRNAs with FC ≥2, which may interact with microRNAs in the miR-200s family, including miR-200a/b/c, miR-141 and miR-429, among the 3 groups (Table 2). A circRNA/miRNA interaction network of selected circRNAs was constructed and is shown in Figure 3A.

Table 2

Related Information of Selected circRNAs.

Group	circRNA ID	CircBase	P-value	FC	Regulation	Type	Length
EC-EU	circ_102362	circ_0008545	1.84E-02	2.27	Up	Exonic	6993
	circ_103470	circ_0067301	3.51E-02	5.66	Down	Exonic	104
	circ_101102	circ_0000419	7.05E-03	2.45	Down	Exonic	1513
EC-C	circ_100332	circ_0014130	1.18E-02	3.95	Up	Exonic	5843
	circ_103237	circ_0063526	2.48E-03	2.35	Up	Exonic	26775
	circ_103362	circ_0065284	6.22E-04	2.15	Down	Exonic	7252
	circ_103130	circ_0001193	7.754E-03	2.01	Down	Exonic	790
EU-C	circ_104700	circ_0005273	3.46E-02	2.86	Up	Exonic	5315

circRNA ID – The circRNA ID is in circRNA microarray profile; CircBase – the circRNA corresponds to that in the circBase (); FC – Fold Change, the absolute ratio (no log scale) of normalized intensities between two conditions; p-value – p-value calculated from paired t-test.

Figure 3

(A) The circRNA/miRNA interaction network analysis of 8 selected circRNAs. The network comprising 4 upregulated (red nodes) and 4 downregulated circRNAs (green nodes) and their target miRNAs and mRNAs is presented. (B) The microarray and qRT-PCR result of selected circRNAs. * P<0.05; ** P<0.01. (C) The circRNA/miRNA/mRNA network analysis of 2 downregulated (green nodes) verified circRNAs, and their target miRNAs and mRNAs are presented. (D, E) GO analysis and KEGG pathway analysis for the 2 circRNAs gathering genes. The vertical axis shows the annotated functions of the target genes. The horizontal axes show the enrichment score (−log2 transformed p-value) and the gene number of each cluster respectively. Only the considerably enriched clusters are listed. (F) Mapping of mTOR signaling pathway. Red pentagram mark nodes are associated with deferentially expressed circRNAs. GO – Gene Ontology; KEGG – Kyoto Encyclopedia of Genes and Genomes (KEGG).

QRT-PCR validation of selected circRNAs

To validate the microarray data, we use qRT-PCR to detect 3 selected circRNAs expression in 20 endometriosis patients, which belonged to the circRNA/miRNA interaction network and differently expressed in EC and EU. Through our confirmation, circ_102293 and circ_103470 had different expressions in endometriosis, but circ_102362 lost significance (Figure 3B).

We constructed a ceRNA-network for the 2 verified circRNAs (Figure 3C), and there are 122 miRNAs and 137 mRNA in the network. In order to explore potential pathways relevant to endometriosis, we use the mRNAs in ceRNA-network for further bioinformatic prediction. In this procedure, KEGG and GO analyses were conducted, and 17 pathways were revealed, such as the mTOR signaling pathway, the Hippo signaling pathway, the HIF-1 signaling pathway, and the PI3K-Akt signaling pathway. which may be involved in the pathogenesis and development of endometriosis. (Figure 3D, 3E). Among these, we focused on the mTOR pathway which we believed to be the most likely to be involved in the pathogenesis of endometriosis (Figure 3F).

Discussion

For over 100 years, endometriosis, as a chronic, estrogen-dependent, inflammatory, heritable disease affecting millions of women across the world with enigmatic etiopathogenesis [1-5]. It is necessary to identify accurate diagnostic biomarkers and effective therapeutic targets. One hypothesis is that the endometriotic cell has undergone epigenetic changes, as DNA methylation and ncRNAs, and those specific changes determine the development into endometriosis [13]. CircRNAs, a type of noncoding RNA, display tissue and cell specificity and are highly conserved. Accumulating evidence has demonstrated that circRNAs are involved in the development of cancers and could be used as diagnostic and prognostic biomarkers [25]. To better understand their intricate roles in endometriosis, circRNAs profiles analysis of endometrium tissue with endometriosis and normal endometrium tissue without endometriosis were conducted and identified several circRNAs that are aberrantly expressed. The results of qRT-PCR were consistent with the circRNAs profiles analysis results, both indicated that the specific circRNAs might promote epithelial-mesenchymal transition and reduce endometriosis. This study highlighted the importance of circRNAs in ovarian endometriosis and provided unique insights into the molecular basis concerning the pathogenesis of endometriosis.

Epithelial-mesenchymal transition of endometrial cells may be a prerequisite for the establishment of endometriotic lesions [8]. MiR-200s family members seem to be strongly associated with pathologic epithelial-mesenchymal transition and to have a metastasis suppressive role [14]. Various studies on miRNAs in endometriosis have identified their cardinal role in the pathogenesis of the disease, such as the miR-200 family [13], and have taken them as potential biomarkers in endometriosis. CircRNAs contain the miRNA targeting site, which can be bind to miRNAs. Previous studies revealed that circVMA21 acted as a sponge of miR-200c and functioned in NP cells through targeting miR-200c and XIAP [15]. Our array-based circRNA profiling bioinformatics analysis identified 5 miRNAs that each had circRNA most likely to interact (i.e., containing complementary base pair sites), and it was found that 8 circRNAs might interact with the miR-200 family. In the current research reports, more and more significant targeted miRNAs have been associated with endometriosis. MiR-503 was reported to target circRNA_101102, induce apoptosis and cell arrest, and inhibit cell proliferation, angiogenesis, and contractility of endometriotic cyst stromal cells [16]. In addition, miR-34, targeted by circRNA_103237, was found to suppress cell proliferation, migration, and cell-cycle maintenance in the early secretory endometrium of women with endometriosis [17].

In an effort to further investigate the specific mechanisms, 2 downregulated circRNAs (circ_103470 and circ_101102) were validated to be significant in endometriosis patients by qRT-PCR. Circ_101102 is derived from RAB3IP which is reported to suppress autophagy in mammalian cells [18]. It has been reported that autophagy is induced in normal endometrial tissue throughout the menstrual cycle, and aberrant autophagy occurs in the eutopic endometrium and ectopic endometriotic foci with endometriosis [19]. Hence, it is speculated that the regulation of this circRNA might contribute to the pathogenesis of endometriosis by inducing autophagy, and this may uncover a potential new treatment for endometriosis.

CircRNAs have been proven to function as ceRNA to interact with miRNAs and influence the expression of miRNA target mRNAs [10,11]. Analysis of the circRNA-miRNA-mRNA network provided potential targets of circRNAs. In the current study, for the first time, we performed functional enrichment analyses and predicted to the target mRNAs of circRNA. In addition, 17 important signaling pathways were delineated, when most of them are consistent with the current knowledge of endometriosis, such as the mTOR signaling pathway, the Hippo signaling pathway, the HIF-1 signaling pathway, and the PI3K-Akt signaling pathway. Endometriosis lesions exist in a unique microenvironment characterized by increased concentrations of hormones, inflammation, and oxidative stress. This environment influences cell survival through the binding of membrane receptors and the subsequent cascading activation of intracellular kinases that stimulate a cellular response. Previous findings indicate that progesterone treatment of endometriotic cells suppresses AKT and ERK1/2 activity, thereby in turn inhibiting mTOR, inducing autophagy, and promoting apoptosis [20,21]. The findings of Song et al. highlighted that dysregulation by genetic inactivation of the Hippo/Yap pathway components or amplification of its downstream effector YAP results in increased endometrial stromal cell proliferation and decreased cell apoptosis [22]. The hypoxia microenvironment stimulates endometrial stromal cell to produce excessive TGF-β1 and activates the TGF-β1/Smad signaling pathway, thus enhancing integral expression and the adhesion ability of endometrial stromal cell [23]. The PI3K-Akt pathway, which integrates a variety of extracellular signals and regulates various cellular functions, including cell growth, differentiation, transformation and survival, is being investigated as a therapeutic target in other diseases and thus may also represent a target for endometrial stromal cell treatment [24]. These signaling pathways can affect the progression of epithelial-mesenchymal transition in cancer-related diseases but are rarely reported in endometriosis-related epithelial-mesenchymal transition mechanisms. In addition, our enrichment analysis revealed some novel dysregulated signaling pathways, such as the insulin signaling pathway and the cGMP-PKG signaling pathway, and their functions should be further investigated.

Overall, our research had several limitations. First and foremost, the sample size was small, which cause random errors to our results. In further research, more samples should be involved. Second, the diagnostic value of circRNA levels needs to be further assessed in other samples, such as blood and ascitic fluid, from women with or without endometriosis. This will broaden the application of the result since this might function not only as a blocking target of endometriosis but also is important as a diagnostic biomarker in the progression of the disease. Last, their pathophysiological relevance to the regulation of circRNAs requires further study, which means experiments should be carried out at the cell, tissue, and individual levels.

To conclude, we believe that this study, dealing with endometriosis, an extremely frequent and debilitating gynecological disease, is a significant step toward understanding epigenetic mechanisms, which is emerging as a very promising and stimulating field of future research.

Conclusions

In general, our results suggested that 2 downregulated circRNAs (circ_103470 and circ_101102) might regulated epithelial-mesenchymal transition in endometriosis via miR-141-5p, which might be a new target for diagnosis or therapy in the future. However, further research is still needed to confirm our conclusion.