Prediction of potential prognostic biomarkers in metastatic prostate cancer based on a competing endogenous RNA regulatory network


 Background

Recently, studies on competing endogenous RNA (ceRNA) networks were becoming prevalent, and circular RNAs (circRNA) have crucial implications for the development and progression of carcinoma. However, studies relevant to metastatic prostate cancer (mPCa) are scanty.
Methods

RNA-Seq data was obtained from MiOncoCirc database and Gene Expression Omnibus (GEO). Differential expression patterns of RNAs were examined using R packages. Circular RNA Interactome, miRTarBase, miRDB and TargetScan were applied to predicting the corresponding relation between circRNAs, miRNAs and mRNAs. The Gene Ontology (GO) annotations were performed to present related GO terms, and Gene Set Enrichment Analysis (GSEA) tools were applied for pathway annotations. Moreover, survival analysis was conducted for hub genes.
Results

We found 820 circRNAs, 81 miRNAs and 179 mRNAs, which were distinguishingly expressed between primary prostate cancer (PCa) and mPCa samples. A ceRNA network including 45 circRNAs, 24 miRNAs and 56 mRNAs was constructed. Besides, the protein-protein interaction (PPI) network was built up, and 10 hub genes were selected by using CytoHubba application. Among the 10 hub genes, survival analysis showed ITGA1, LMOD1, MYH11, MYLK, SORBS1 and TGFBR3 were significantly connected with disease-free survival (DFS).
Conclusion

The circRNA-mediated ceRNA network provides potential prognostic biomarkers for metastatic prostate cancer .


Abstract Background
Recently, studies on competing endogenous RNA (ceRNA) networks were becoming prevalent, and circular RNAs (circRNA) have crucial implications for the development and progression of carcinoma.
However, studies relevant to metastatic prostate cancer (mPCa) are scanty.
Methods RNA-Seq data was obtained from MiOncoCirc database and Gene Expression Omnibus (GEO).
Differential expression patterns of RNAs were examined using R packages. Circular RNA Interactome, miRTarBase, miRDB and TargetScan were applied to predicting the corresponding relation between circRNAs, miRNAs and mRNAs. The Gene Ontology (GO) annotations were performed to present related GO terms, and Gene Set Enrichment Analysis (GSEA) tools were applied for pathway annotations. Moreover, survival analysis was conducted for hub genes.

Results
We found 820 circRNAs, 81 miRNAs and 179 mRNAs, which were distinguishingly expressed between primary prostate cancer (PCa) and mPCa samples. A ceRNA network including 45 circRNAs, 24 miRNAs and 56 mRNAs was constructed. Besides, the protein-protein interaction (PPI) network was built up, and 10 hub genes were selected by using CytoHubba application. Among the 10 hub genes, survival analysis showed ITGA1, LMOD1, MYH11, MYLK, SORBS1 and TGFBR3 were significantly connected with disease-free survival (DFS).

Conclusion
The circRNA-mediated ceRNA network provides potential prognostic biomarkers for metastatic prostate cancer .

Background
In men, prostate cancer (PCa) is a type of malignant tumor of the urogenital system with high morbidity. Patients with de novo metastatic PCa (mPCa) have a high risk of disease-specific mortality [1]. Previous studies have reported prostate-specific antigen (PSA) testing could contribute to reducing the number of patients with metastatic cancer. Despite screening for earlier diagnosis, the number of men who develop metastases is still large [2]. Many studies have reported ceRNA networks were crucial for PCa biology and therapeutics [3,4]. Nevertheless, there is no enough understanding of the molecular mechanism that may cause mPCa. Accordingly, it is essential to explore effective biomarkers, which could improve the perception of development and progression of mPCa.
CircRNA belongs to the non-coding RNA formed from a non-canonical splicing event called backsplicing, and their covalently closed ring structure prevents them from exonuclease-mediated degradation resulting in high stability [5]. Relevant study showed that circRNAs acted as miRNA sponges, functioned through interactions with proteins, and participated in protein synthesis [4].
Accumulating evidence suggested circRNAs had impacts on physiological development and multifarious diseases, such as brain development and diseases, cardiovascular diseases and cancers [6].
A ceRNA hypothesis was raised by Pandolfi et al. indicating that all types of RNA transcript harboring miRNA response elements (MREs) could serve as one type of ceRNA, which may influence the tumorigenesis of malignant tumors, such as mRNAs, pseudogenes, circRNAs and lncRNAs [7]. Certain study showed ceRNA interactions had implications for gastric cancer, mammary cancer, melanoma, spongioblastoma, prostate cancer and so on [8].
In this study, RNA-seq data was obtained from MiOncoCirc and GEO database. We identified the circRNAs that acted as miRNA sponges and found target genes for miRNA respectively. Following, a ceRNA network was built up, including 45 circRNAs, 24 miRNAs and 56 mRNAs. GO and GSEA were applied for differentially expressed mRNAs. Finally, we established a PPI network to better comprehend the molecular mechanism and pathogenesis of mPCa. This study aims to figure out the relations between ceRNA and the carcinogenesis and progression of mPCa, and provide further insight into the improvement for treatment and prognosis of mPCa. The flow diagram of the whole ceRNA analysis was represented in Fig. 1.

MiOncoCirc Compendium
MiOncoCirc Compendium is an accessible compendium of cancer-focused circRNAs for the scientific community [9]. In this study, the circRNA profiles of 88 mPCa tissues were retrieved from MiOncoCirc (https://MiOncoCirc.github.io/). A complete description of the donor age range, sample location and cell lines was obtained from relevant annotation documents.

GEO database
The circRNA data of 144 primary PCa tissues was downloaded from GSE113120 dataset of GEO database. The expression patterns of miRNA and mRNA were acquired from GSE21036 (99 mPCa tissues and 14 primary PCa tissues) and GSE21034 (131 mPCa tissues and 19 primary PCa tissues) respectively.
Screening of differentially expressed genes (DEGs) R packages were employed to identify DEGs in mPCa tissues and primary PCa tissues, and the differences were described by fold-change (FC) and false discovery rate (FDR). |Log2FC|>5 and FDR < 0.05 were deemed significant for circRNAs. The differential expression levels of mRNAs and miRNAs were analyzed with cutoff values of |Log2FC|>1 and FDR < 0.05. Then, we use pheatmap packages to draw heatmaps for obtained DEGs.

Construction of ceRNA network
Circular RNA Interactome was applied to predicting circRNA-miRNA interactions. MiRDB, miRTarBase and TargetScan databases were employed to explore -corresponding mRNAs for miRNAs. At last, a ceRNA network was constructed by selected circRNAs, miRNAs and mRNAs.

PPI network construction
We employed the String database to analyze the interactions between proteins of DEmRNAs. Then the network was visualized by using Cytoscape v3.7.2. Cytohubba was used to obtain 10 hub genes ranked by degree.

GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis
ClusterProfiler packages were employed to identify GO terms. GO consists of three aspects of genes: biological process, cellular component and molecular function. GSEA was applied to analyzing KEGG pathways at a significant level which was set at p < 0.05.

Survival analysis
For the purpose of identifying specific correlations of the expression levels of 10 hub-genes with patient prognosis, Gene Expression Profiling Interactive Analysis (GEPIA) was applied for survival analysis. We set group cutoff as median, calculated the hazards ratio based on the Cox PH Model, and added the 95% confidence interval as a dotted line. Finally, the corresponding DFS was obtained.

Identification of differentially expressed circRNAs (DEcircRNAs)
We examined the expression level of circRNAs between 144 primary PCa and 88 mPCa samples.
After intersecting with DEmiRNAs, 45 of 62 DEcircRNAs were found to target 24 of 81 DEmiRNAs.
Then, miRDB, miRTarBase and TargetScan were adapted for predicting target mRNAs of 24 circRNAtargeted miRNAs, and only mRNAs identified by at least two of these databases were selected. The selected target mRNAs then were overlapped with DEmRNAs. Finally, we got a ceRNA regulatory network constructed by 45 circRNAs, 24 miRNAs and 56 mRNAs, and made it visualize by using Cytoscape v3.7.2 (Fig. 5).

PPI network construction
We identified the mutual effect between proteins of DEmRNAs through String online database, and constructed a PPI network including 21 nodes and 31 edges (Fig. 6a). Then, we evaluated the degree between proteins to find hub genes (Fig. 6b). Based on the rank method of Cytohubba, we got 10 hub genes, involving ITGA1, MYH11, MYLK, SORBS1, CALD1, LMOD1, TGFB3, DMD, F5, TGFBR3.

Functional assessment
To probe the biological process of target mRNAs, ClusterProfiler packages of R were used to make GO analysis (Table 4). We found the up-regulated mRNAs were related to interleukin-7-mediated signaling pathway, chromatin silencing at rDNA, cellular response to interleukin-7 and protein heterooligomerization. In addition, we found base excision repair, cell cycle, one carbon pool by folate, DNA replication, mismatch repair and homologous recombination were related to genes that up-regulated in mPCa samples by KEGG-GSEA (Fig. 8).

Survival analysis of hub genes
We analyzed DFS of the 10 genes by GEPIA with a cutoff value of p < 0.05. Results indicated the expression levels of ITGA1, LMOD1, MYH11, MYLK, SORBS1 and TGFBR3 were positively associated with DFS in patients with PCa (Fig. 9).

Discussion
Recently, on account of high attention paid to the ceRNA hypothesis, the development of the ceRNA network is greatly fast, contributing to a deep perception of the function and mechanism of ceRNAs.
Increasing evidence showed that ceRNAs played remarkable parts in tumor development and advancement via changing the expression of pivotal oncogenic or tumor-suppressive genes [10][11][12]. CircRNAs often act as effective clinical biomarkers owing to the resistance for exonucleases degradation [9]. For example, circRNAs were extremely steady in mammalian cells, and hsa_circ_002059 could be a fresh type of diagnostic biomarker for stomach cancer [15].
Hsa_circ_0003221 enhanced the cell multiplication and migration and appeared to be an oncogenic circRNA with the potential for therapeutic target for bladder cancer [16]. In mPCa, relevant researches reported the down-regulated Circ-ITCH was connected with a high chance for lymphatic metastasis and poor outcome, and circZNF609 could up-regulate miR-186-5p to inhibit development, migration, and invasion of mPCa [17,18]. However, the details of circRNAs in mPCa are not entirely illustrated.
The purpose of our research is to explore a novel way which could promote comprehensive recognition of the underlying molecular mechanism of circRNAs in mPCa by identifying the relations between circRNAs, miRNAs and mRNAs.
Increasing evidence has proved the important effect of miRNAs in cancer cell proliferation, metastasis and immune evasion, which could be used as possible diagnostic and therapeutic methods in various cancers [19][20][21]. In the present research, we found 24 DEmiRNAs in the ceRNA network. Certain studies have demonstrated that certain DEmiRNAs were necessary for tumor development. MiR-296-5p was considered to be a valuable miRNA, and a report showed the overexpression of miR-296-5p restricted tumor development and aggression of HCC via inactivating NRG1-Fra-2 signaling [22]. In colorectal cancer (CRC), cellular experiments have proved downregulated miR-548c-5p restrained cancer cell growth and reproduction of inflammatory cytokines via directly targeting PGK1 [23].
Downregulation of miR-548c-5p could predict poor survival in CRC [24]. Similarly, miR-665 could predict poor survival and facilitate metastasis of breast cancer via the NR4A3-MEK signaling pathway [25]. CircABCC2 had the potential to bond miR-665 to manage ABCC2 expression, influencing the development and progression of HCC [26]. Thus, miRNA plays an essential part in cancer progression, and the 24 DEmiRNA involved in the ceRNA network may have inseparable relations with mPCa.
Among the 56 DEmRNAs in ceRNA regulatory network, the multivariable analysis indicated MEIS2 expression correlated with reduction for metastatic progression of PCa [27]. Yan et al. demonstrated CNTN1, a type of cell adhesion protein in the neural system, could promote PCa cell invasion, play a part in tumor formation and enhance metastasis [28]. The oncoprotein STMN1 was able to boost PCa cell growth and aggression. And miR-34a downregulated STMN1 to inhibit cancer proliferation and colony formation [29]. DDR2 in PCa may stimulate bone resorption by up-regulating RANKL, and promote the adhesion between cancer cells and type I collagen, so it was able to induce bone metastasis and could be a new target for PCa treatment [30]. The flow diagram of ceRNA analysis Figure 2 The heatmaps of DEcirRNAs Figure 3 The heatmaps of DEmiRNAs The heatmaps of DEmRNAs  The PPI network and corresponding hub genes. a. The network includes 21 nodes and 31 edges. b. 10 hub genes extracted from a. The node color changes gradually from yellow to red indicating the ascending order of the degree of the genes.