Profiling and bioinformatics analyses of differential circular RNA expression in prostate cancer cells

Aim: There is little knowledge about the expression profile and function of circular RNAs (circRNAs) in prostate cancer (PCa). Methods: The expression profiles of circRNAs in RWPE-1, 22RV1 and PC3 cells were explored via high-throughput circRNAs sequencing and validated by real-time qPCR. The roles of differentially expressed circRNAs were evaluated by bioinformatics analyses. Results: Altogether 9545 circRNAs were identified and hundreds of differentially expressed circRNAs were recognized. CircRNA–miRNA networks analysis showed that many circRNAs, including circSLC7A6, circGUCY1A2 and circZFP57 could cross-talk with tumor-related miRNAs such as miR-21, miR-143 and miR-200 family. Conclusion: The results of our bioinformatics analyses suggested that circRNAs should play critical roles in the development and progression of PCa.


RNA-seq data analysis
The FASTQ reads were aligned to the human reference genome. Transcriptome data were obtained from the UCSC genome database by using the BWA-MEM software [21], circRNAs were identified by using CIRI software as described [22] and the identified circRNAs were then annotated with circDB database, a manually curated database that collected all published circular RNAs and was developed by Cloudseq Inc. (For each sample, the original backspliced junction reads of each circRNA were converted to tags per million reads.) Then, statistical hypothesis based on binomial distribution was used to identify the differentially expressed circRNAs between two sample groups by edgeR analysis method. p-value < 0.05 was considered as significantly differential expression [23,24]. Hierarchical clusterings were performed by using the euclidean distance matrix of the heatmap.2 function in the gplots package of the R environment to analyze the expression pattern of significantly up and downregulated circRNAs in each group with the threshold of fold change 2.0 and p-value < 0.05. The function of circRNAs was predicted by Gene Ontology (GO) term and KEGG pathway enrichment analyses of the corresponding host genes with DAVID database. The results of GO term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were shown with threshold p < 0.05 and arrangement in order based on false discovery rate (FDR) value.

circRNA-miRNA regulatory network
The significantly differentially expressed circRNAs were validated with real time qPCR, and then subjected to analysis. Differentially expressed circRNAs were used to predict the potential miRNA response elements and the binding sites of miRNAs with CloudSeq's home-made software based on miRanda and TargetScan (CloudSeq Inc.). The network between circRNAs and miRNAs was also constructed by Cytoscape based on the binding sites of the differentially expressed circRNAs and miRNAs. Different nodes represent circRNAs and miRNAs, respectively. Solid lines between two nodes represent potential binding.
Real-time qPCR cDNA was synthesized from 500 ng total RNA with the PrimeScript RT Master Mix (Takara), according to manufacturer's instructions. The real-time PCR analyses of the expression level of the circRNAs were performed by using SYBR Premix Ex Taq II (Takara). RT-qPCR was performed in 20 μl reaction volume, including 1 μl cDNA, 10 μl 2 × Master Mix, 0.3 μl Forward Primer (10 μM), 0.3 μl Reverse Primer (10 μM) and 8.4 μl double distilled water. The primer sequences are presented in Supplementary Table 3. The reaction was set at 95 • C for 10 min, and then repeated for 40 cycles at 95 • C for 10 s and 60 • C for 60 s in real time PCR System (ABI, CA, USA), using β-actin as a reference. All samples were amplified in triplicate wells, and the relative level was calculated with 2 -Ct method.

Agarose gel electrophoresis
The products of real-time qPCR or PCR were validated in agarose (Biowest, Spain) with a concentration of 1.5% in TAE buffer (made up of Tris base, acetic acid, and EDTA; Songon Biotech, Shanghai, PR China) by use of an electrophoresis system (BIO-RAD, CA, USA). The results were detected by using an imaging system (Syngene, Cambridge, UK).

Sanger sequencing
The PCR amplified products were Sanger sequenced by standard methods in the Beijing Genomics Institute (Beijing, PR China).

Overview of circRNA profiles in PCa cell lines
To explore the profile of circRNAs in prostate carcinogenesis, biological triplicates, RNase R digested and rRNA depleted, circRNA-specific sequencing was performed in RWPE-1, 22RV1 and PC3 cells. More than 20 million

Figure 2. Profiles of differentially expressed circular RNAs in prostate cancer cell lines. (A-C)
The distribution of circRNAs from different catalogs in RWPE-1, 22RV1 and PC3 cell lines. (D-F) Hierarchical clusterings shown as the analysis of differential circRNAs in group 22RV1 versus RWPE-1, PC3 versus RWPE-1, PC3 versus 22RV1, where the red strips represent high relative expression, the green strips represent low relative expression and the dendrograms show the relationships between the samples and differential circRNAs, with the threshold of fold change 2.0 and p-value < 0.05. (G-I) The visualization of circRNAs between two conditions in group 22RV1 versus RWPE-1, PC3 versus RWPE-1, PC3 versus 22RV1, where the red rectangles represent differential expression of circRNAs with the threshold of fold change 2.0 and p-value < 0.05. raw reads on average were detected in each cell sample and more than 90% of reads were mapped (Supplementary Table 1). The number of circRNAs decreased with the increase of the mean junction reads in each cell sample (Supplementary Figure 2). A total of 9545 circRNAs were identified in the three cell lines, and 845 of them were coexpressed ( Figure 1A). In these circRNAs, 53.9% (5142/9545) were novel, and the others (4403/9545, 46.1%) were annotated from Circbase database, Zhang 2014 [25] and Guojunjie 2014 [4] reports ( Figure 1B)  and others (1737/9545), of which exon-derived circRNAs were the most (63.3%) ( Figure 1C). Meanwhile, they ranged in length from 17 to 499,794 nt, and the longest length group (>10,000 nt) accounted for the largest proportion (3745/9545, 39.2%) ( Figure 1D). The distribution of host genes in chromosomes was also significantly different, most (963/9545, 10.4%) in chromosome 1 and least (30/9545, 0.31%) in chromosome Y ( Figure 1E). In addition, many host genes could produce a diversity of circRNAs as a result of variable splicing, and the maximum number of variable splices was as great as 34 (Supplementary Figure 3).

Identification of differentially expressed circRNA profiles
To identify the differentially expressed circRNAs, we used statistical analysis in the groups of 22RV1 versus RWPE-1, PC3 versus RWPE-1, PC3 versus 22RV1. Venn diagrams showed the distribution of the three cell lines in each catalog and these circRNAs were mainly exon-derived and antisense-derived (Figure 2A-C). The expression pattern of hierarchical clusterings was distinguishable between different cell samples, and the dendrograms showed relationships between the samples and differential circRNAs ( Figure 2D-F). Volcano plots were used to find out the target circRNAs in each group with p-value < 0.05, fold change >2.0 ( Figure 2G Table 1). The distribution of the differentially expressed circRNAs derived from different chromosomes and catalogs are shown in Table 2. Most host genes of differentially expressed circRNAs were still exon-derived and in chromosome 1. circRNAs were treated with RNase R. As expected, the circRNAs were resistant to RNase R treatment in contrast to linear control ( Figure 3C).

Prediction of circRNA functions in PCa
To further explore the functions of circRNAs in the biological progression of PCa, we performed GO term and KEGG pathway analyses on the host genes of differentially expressed circRNAs.
All comparison data of GO term analyses including biological process, molecular function and cellular component are shown in Supplementary Figure 6, 7 and 8. Our biological process analyses showed that some host genes could regulate the molecular metabolism process, catabolic process, cell cycle G2/M phase transition, small GTPase mediated signal transduction and the characteristics of stem cells. Some other host genes were involved in regulation of translation, protein ubiquitination, kinase activity and signaling pathways, such as integrin-mediated signaling, TOR signaling pathway, Hippo signaling pathway and Ras protein signal transduction. Molecular function analyses showed that host genes had the ability to bind many kinds of molecules, including SH3 domains, p53, enzymes, DNA, RNA and ATP ( Figure 4A-F).
KEGG pathway revealed that some host genes were associated with cancer progression by participating in proteoglycan and transcriptional misregulation. In addition, the analyses also showed that some genes had enrichment on a variety of tumor-related signaling pathways (including Hippo, Rap1, ErbB, Wnt, PI3K-Akt, HIF-1, Notch and TGF-β signaling pathways). All comparison data are shown in Supplementary Figure 9. The significantly important pathways with high enrichment scores were 'Hippo signaling pathway' in the group of PC3 versus RWPE-1 (downregulation) and 'Rap1 signaling pathway' in the group of 22RV1 versus RWPE-1 (downregulation) ( Figure 4G-J).

Construction of circRNA-miRNA interaction networks
To illustrate the complex regulation mechanism of circRNAs in PCa cells, we constructed two circRNA-miRNA networks. First, we constructed a circRNA-miRNA network with differentially expressed circRNAs listed above and the five miRNAs with the stronger binding ability to them ( Figure 5A). To find out the target circRNAs associated with PCa more effectively, we selected 46 miRNAs closely related to the development of PCa (including miR-21, miR-221/222 and the miR-200 family) to draw a network diagram for showing the relationship between them and the corresponding circRNAs ( Figure 5B).

Discussion
CircRNAs are natural endogenous RNAs, the roles of which, in some human cancers, have been gradually unraveled in recent years [1]. However, the role of circRNAs in PCa is not clear yet. In the present study, we reported the circRNA expression profiles in PC3, 22RV1 and RWPE-1 cell lines and explored the differentially expressed circRNAs to predict the potential function of circRNAs in PCa by bioinformatics analyses. Unlike mRNA participating in signaling of the biological process by translating into protein, circRNAs can regulate some future science group www.future-science.com PC3 RWPE-1 important biological activities via different mechanisms, such as regulating the content of miRNAs as a reservoir of miRNAs [2], combining with protein to form a circRNA-protein complex to disrupt the effects of protein on the targets [26], regulating host gene transcription by binding with RNA polymerase II or miRNA [27] and affecting the formation of gene splicing protein [28]. They can also be translated into protein directly in a few cases [6,7]. The most important mechanism of circRNAs is known to function as miRNA sponges like competitive endogenous RNA molecules [2]. miRNAs are important regulators in gene expression and play a crucial role in cancer progression. In this study, we constructed circRNA-miRNA networks by using bioinformatics analyses to explore the relationship between differentially expressed circRNAs and miRNAs. Some important miRNAs were found to bind with the significantly differential circRNAs listed above, such as miR-143 (binding with circSLC7A6, hsa circ 0039943), miR200c and miR-200b (binding with circGUCY1A2, hsa circ 0008602) and miR-21 (binding with circZFP57). According to previous studies, miR-143 was associated with bone metastasis of PCa and involved in the regulation of epithelial-to-mesenchymal transition [29]; the miR-200 family could inhibit PCa progression [30]; and miR-21 could promote PCa development [31], implying that circGUCY1A2 and circZFP57 may affect the development of tumor formation and circSLC7A6 may promote bone metastasis through these miRNAs. The selected 46 miRNAs, which were reported to be closely associated with PCa progression [30][31][32][33][34], or with AR signaling pathway [35][36][37][38][39][40] and epithetlial-to-mesenchymal transition [29,[41][42][43][44][45], could also be found to have corresponding circRNAs. These results indicate that circRNAs may have a close relationship with the development of PCa.
Another important mechanism of circRNAs is that they can affect the expression or function of host genes. GO term analyses of this study showed that host genes of the differentially expressed circRNAs functioned in some important biological processes and molecular mechanisms, such as binding of the SH3 domain and participating in integrin-mediated signaling pathway. DOCK1, for example, has a SH3-SH2-domain in LNCap cells, and ROBO1 In tr a c e ll u la r In tr a c e ll u la r p a r t In tr a c e ll u la r o r g a n e ll e In tr a c e ll u la r o r g a n e ll e p a r t O r g a n e ll e p a r t T r a n s fe r a s e a c ti v it y , tr a n s .. .

Cellular component Biological process Molecular function
In tr a c e ll u la r p a r t In tr a c e ll u la r In tr a c e ll u la r o r g a n e ll e p a r t O r g a n e ll e p a r t In tr a c e ll u la r o r g a n e ll e     interacts with DOCK1 to negatively regulate the Rac activation at this domain [46]. FUT8 provides a versatile Nglycosylated protein for the formation of N-glycan branches. N-glycan is a basic and generic protein modification in mammals and plays a key role in various physiological and pathological events including cancer progression [47]. Branch variations are often found in cancer cells and are deeply involved in cancer growth, invasion and metastasis.
In addition, hundreds of host genes such as DOCK1, FUT8, MPP6, GSK3B, SETD3 are capable of binding with proteins, and therefore further study is warranted to see whether the corresponding circRNAs of these host genes would function by binding with protein.
future science group www.future-science.com KEGG pathway analysis of this study demonstrated that many host genes were involved in Hippo, Rap1, ErbB, Wnt, PI3K-Akt, HIF-1, Notch and TGF-β signaling pathways, all of which are known as important mechanisms for cancer progression. The Hippo signaling pathway, also known as the Salvador/Warts/Hippo pathway, is an evolutionarily conserved pathway that controls cell proliferation, apoptosis and stem cell self-renewal [48]. Its main function is phosphorylating the main modulation proteins, YAP and TAZ, in this pathway [49]. Once dephosphorylated, the coactivator will translocate to the nucleus and promote the expression of genes that contribute toward cell proliferation and survival, which in turn accelerates cancer progression [50]. Through Hippo signaling, Wnt signaling promotes androgen-independent prostate cancer cell proliferation [51], α3β1 integrin suppresses prostate cancer metastasis [52], and phosphodiesterase 5/protein kinase G signal governs stemness of prostate cancer stem cells [53]. Rap1, a small GTPase in the Ras-related protein family, is regulated by binding to GTP or GDP and functions as a switch during cellular signaling transduction, which has many important effects in tumor cell migration, invasion and metastasis [54]. In PCa, Rap1 signaling mainly regulates PCa cell adhesion and promotes PCa metastasis [55,56].
In the present study, we detected the differential expression of circRNAs in PCa cells. CDR1as (hsa circ 0001946) was validated to highly express in PC3. CDR1as is a circRNA containing 74 binding sites of miR-7, and its typical 'miRNA Sponge' mechanism has been studied in depth [2]. The downstream genes of miR-7 include KLF4, RAF1, PAK1, IRS2, IRS1, AKT, ACK1, FAK, IGF1R, HNF4, SNCA, mTOR, NOTCH1 and PIK3CD [57]. All of these are indispensable tumor suppressors or oncogenes in tumorigenesis of various cancers. In the sequencing results, the expression of CDR1as in PC3 was nearly 200 times higher than that in normal prostate epithelial cells  and prostate carcinoma epithelial cells, suggesting that CDR1as may substantially contribute to bone metastatic progress of PCa. CircBAGE2 (hsa circ 0061259) is significantly upregulated in 22RV1 compared with RWPE-1. The binding miRNA miR-103a can suppress tumor cell proliferation by targeting PDCD10 in PCa [58], suggesting that circBAGE2 may have a close relevance with the progression of PCa.

Conclusion
From the data above, we conclude that circular RNAs may have a close relationship with the development and progression of PCa. CircRNAs such as circGUCY1A2, circZFP57 and circBAGE2 may play an important role in pathogenesis, and CDR1as and circSLC7A6 may be associated with bone metastasis. All in all, the present study revealed a comprehensive expression and functional profile of differentially expressed circRNAs in PCa cell lines, which may provide a new direction for the study of PCa.
future science group www.future-science.com

Future perspective
With deeper research into this field, the function of circRNAs has been paid more and more attention. The current study shows that the expression of circRNAs in PCa cells is definitely abundant and bioinformatics analyses also predict that circRNAs may play an important role in PCa. However, the profile of differentially expressed circRNAs in the specimens of PCa patients is unknown, and the real tumor biology functions of circRNAs in PCa remain elusive. These should be explored in the future. The interactions between mRNAs, miRNAs and circRNAs constitutes a huge competitive endogenous RNA system network, the regulatory mechanism of which also needs to be further explored.
Future Sci. OA (2018) 4(9) future science group • Circular RNAs (CircRNAs) are involved in many important biological processes, as well as the pathogenesis of many cancers. The role of circRNAs in prostate cancer (PCa) is unknown. • The current study aimed to explore the profile of circRNAs in PCa cell lines and to predict their potential function.
• In total, 9545 circRNAs and hundreds differentially expressed circRNAs were identified in PCa cell lines.
• CircGUCY1A2, circZFP57 and circBAGE2 should be important in the pathogenesis of PCa.
• CDR1as and circSLC7A6 should be associated with bone metastasis in PCa.
• The profiling and differentially expressed circRNAs in the specimens of PCa patients should be explored in the future. financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Financial & competing interests disclosure
No writing assistance was utilized in the production of this manuscript.
Authors' contributions

Ethical conduct of research
The authors state that they have obtained appropriate institutional review board approval or have followed the principles outlined in the Declaration of Helsinki for all human or animal experimental investigations. In addition, for investigations involving human subjects, informed consent has been obtained from the participants involved.

Open access
This work is licensed under the Creative Commons Attribution 4.0 License. To view a copy of this license, visit http://creativecomm ons.org/licenses/by/4.0/