Five Circular RNAs in Metabolism Pathways Related to Prostate Cancer

Prostate cancer (PCa) is the most common malignant tumor in men, and its incidence increases with age. Serum prostate-specific antigen and tissue biopsy remain the standard for diagnosis of suspected PCa. However, these clinical indicators may lead to aggressive overtreatment in patients who have been treated sufficiently with active surveillance. Circular RNAs (circRNAs) have been recently recognized as a new type of regulatory RNA that is not easily degraded by RNases and other exonucleases because of their covalent closed cyclic structure. Thus, we utilized high-throughput sequencing data and bioinformatics analysis to identify specifically expressed circRNAs in PCa and filtered out five specific circRNAs for further analysis—hsa_circ_0006410, hsa_circ_0003970, hsa_circ_0006754, hsa_circ_0005848, and a novel circRNA, hsa_circ_AKAP7. We constructed a circRNA-miRNA regulatory network and used miRNA and differentially expressed mRNA interactions to predict the function of the selected circRNAs. Furthermore, survival analysis of their cognate genes and PCR verification of these five circRNAs revealed that they are closely related to well-known PCa pathways such as the MAPK signaling pathway, P53 pathway, androgen receptor signaling pathway, cell cycle, hormone-mediated signaling pathway, and cellular lipid metabolic process. By understanding the related metabolism of circRNAs, these circRNAs could act as metabolic biomarkers, and monitoring their levels could help diagnose PCa. Meanwhile, the exact regulatory mechanism for AR-related regulation in PCa is still unclear. The circRNAs we found can provide new solutions for research in this field.


INTRODUCTION
Prostate cancer (PCa) is a slow-growing malignant tumor, the incidence of which increases with age (Carroll and Mohler, 2018;Etzioni and Nyame, 2020). At the beginning of diagnosis, most patients are asymptomatic; however, it is still among the top three causes of cancer-related deaths in men (Siegel et al., 2017). Patients with a high risk of PCa must undergo periodic testing for serum prostate-specific antigen (PSA). Tissue biopsy remains the care standard for diagnosis for suspected PCa (Litwin and Tan, 2017). After the tumor is confirmed by biopsy, the next step is to determine the invasiveness of tumor cells. The Gleason score is the most commonly used scale to assess the grade of PCa. When the grade is high, tumors tend to spread. Most Gleason scores used to evaluate prostate biopsy samples range from 6 to 10. A score of 6 indicates low-risk PCa; a score of 7 indicates intermediate-risk PCa; a score of 8 to 10 indicates high-risk PCa (2019). However, researchers have found that the "normal" PSA level of 0-4 ng/mL does not guarantee cancer-free status; in approximately 25% of men with a PSA below 4 ng/mL, a biopsy still reveals PCa (Kitagawa et al., 2014). Thus, these clinical indexes cannot guarantee the reliability of diagnosis. New assistant biomarkers need to be developed for the diagnosis of PCa.
With the development of sequencing technology, circRNAs are recognized as a new type of regulatory RNA. They were first identified by analysis using next generation RNA sequencing (RNA-seq) in a study of pediatric acute lymphoblastic leukemia (Salzman et al., 2012). Most circRNAs are composed of protein-coding exons; thus, the expression of these circRNAs competes with the production of pre-mRNAs. These events also lead to the expression of circRNAs being higher than that of their cognate linear RNAs under certain conditions (Jeck et al., 2013;Jeck and Sharpless, 2014). CircRNAs can function by directly regulating gene expression or by acting as miRNA sponges (Tay et al., 2014). They are similar to competitive endogenous RNAs (ceRNAs) and contain miRNA response elements (MREs). Therefore, they can function by competing with mRNAs to bind miRNAs (Hansen et al., 2013). For example, dysregulation of circRNA-0001946 contributes to tumor cell proliferation and metastasis in colorectal cancer by targeting microRNA-135a-5p (Deng et al., 2020). CircRNAs are not easily degraded by RNases and other exonucleases due to a covalent closed cyclic structure without free 5 or 3 ends. They have a longer half-life (>48 h) than linear RNAs (Suzuki et al., 2006;Jeck and Sharpless, 2014). CircRNAs have been known to be rich in tumors (Salzman et al., 2012). Therefore, compared with other RNAs, circRNAs have more advantages as novel biomarkers of cancer and other diseases (Arnaiz et al., 2019).
Studies have shown that circRNAs are functional in PCa. Overexpression of circ0005276 and its host gene X-linked inhibitor of apoptosis protein (XIAP) can promote cell proliferation, migration, and epithelial-mesenchymal transition in PCa tissues compared with that in normal tissues (Feng et al., 2019). CircRNAs can act as oncogenes in the progression of PCa and are differentially expressed between cancer tissues and normal tissues (Feng et al., 2019). It is also reported that circRNAs can act as therapeutic targets. For example, the overexpression of circRNA cir-ITCH significantly inhibits the proliferation, migration, and invasion of PCa cells. By targeting miR-17 in PC-3 and LNCaP cell lines, circRNAs could act as therapeutic targets in PCa, especially in castration-resistant prostate cancer (CRPC) . CircRNAs also affect carbohydrate, lipid, and amino acid metabolism in cancer. By regulating transcription factors, circRNAs can modulate glycolysis (Yu et al., 2019). Thus, it is important to identify differentially expressed circRNAs in PCa and explore their potential as diagnostic and therapeutic targets in cancer.
In this study, we compared circRNAs between four PCa tissues and two adjacent normal tissues of two PCa patients by sequencing six sets of RNA-seq. We selected five circRNAs that were highly expressed in tumor tissues and found that the fold change in expression of these five circRNAs was significantly higher than that of their cognate linear RNAs. We verified these circRNAs by PCR in PCa cell lines and used the circRNA-miRNA-mRNA method to predict biological pathways regulated by these circRNAs. Some well-known pathways in PCa were enriched, such as the p53 signaling pathway, MAPK signaling pathway, hormone-mediated signaling pathway, and cellular lipid metabolic process. These pathways also confirmed the high reliability of the five circRNAs that participated in the regulation of PCa.

Patients and Samples
For sequencing samples, two pairs of PCa tissues and adjacent tissues were derived from surgical samples. Sections from normal and malignant tissues were examined after staining with hematoxylin and eosin. The tumor specimen comprised >80% malignant cells, and the benign specimen comprised an approximately equal admixture of normal epithelial and stromal cells. The pathology of the prostate tumor was checked by a pathologist and established as a combined Gleason Score of 6 (3 + 3), stage T2a, with focal involvement of the surgical margin. RNA was purified from minced frozen tissue using Trizol reagent (Life Technologies, Inc., Rockville, MD, United States). Total RNA was briefly treated with DNase I. For each sample, an RNA library was constructed using 3 µg total RNA. Ribo-Zero Gold Kits were used to remove rRNA. According to the instructions of the NEB-Next Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ispawich, United States), different index tags were selected. The constructed libraries were sequenced using Illumina; the sequencing strategy used was PE150.
For qRT-PCR, 20 pairs of PCa tissues and adjacent normal tissues were collected from the Department of Pathology of Beijing Hospital with Gleason scores of 6 (14 cases) and 7 (6 cases). All tissues were fixed in phosphate-buffered formalin, dehydrated with ethanol, and embedded in paraffin. The malignant status and Gleason score were obtained for these samples by histological analysis. The work was approved by the Beijing Hospital Ethics Committee.
None of these patients had undergone hormonal therapy prior to surgery.

Quality Control and Mapping of Sequencing Data
The quality of the fastq data of RNA-seq was evaluated using fastQC (Andrews, 2010). We found some reads mixed with adapters, and then Trimmomatic (Bolger et al., 2014) was used to filter these sequences. The reads with lengths of less than 28 were dropped. The read quality was filtered through a four-base sliding window with an average quality threshold of 15. After the reads passed the sequence quality tests, the filtered reads were mapped to the human hg38 genome (Lander et al., 2001) using the aligner software STAR (Dobin et al., 2013) with parameter "-chimSegmentMin 10."

Differential Expressed circRNAs Filtered
CIRCexplorer2  and CLEAR/CIRCexplorer3 (Ma et al., 2019) were mainly used to obtain circRNAs in our research. For CIRCexplorer2, we used the parse module to analyze circRNA fusion junction reads and annotate modules for circRNA gene information. The expression of circRNAs was quantified by CIRCscore in CLEAR/CIRCexplorer3, which indicates the circRNA expression level by linear RNA expression level adjustment. The expression of fold change between tumor circRNAs and normal adjacent prostatic tissue circRNAs was calculated to filter tumor-specific expressed circRNAs.

Bioinformatics Analysis of CircRNAs
All of the interaction binding sites between circRNAs and miRNAs were downloaded from circBank . For the novel circRNAs, we used miRNDB (Chen and Wang, 2020) to predict the related miRNAs. We identified the function of the predicted miRNAs by manual literature mining. Then, Cytoscape (Shannon et al., 2003) was used to build a network between circRNAs and miRNAs. For differentially expressed mRNAs, we chose featureCounts (Liao et al., 2014) to quantify read counts for each gene. Based on paired-end data, "requireBothEndsMapped = TRUE" and "isPairedEnd = TRUE" were set additionally. Then, we calculated the normalized expression levels in fragments per kilobase per million mapped reads (FPKM) by using the DGEList and rpkm function from edgeR (Robinson et al., 2010).
To predict circRNA-related pathways, we regarded miRNAs as a middleman to find circRNA-related mRNAs. The interactions between miRNAs and mRNAs were obtained from miRDB (Chen and Wang, 2020). Meanwhile, it provided interaction scores to assess accuracy; only scores higher than 90 and differentially expressed mRNAs were considered for further analysis.

Functional Enrichment Analysis and Survival Analysis
The circRNA-related mRNA list was analyzed using the functional enrichment tool GOseq (Young et al., 2010). Compared with other tools, GOseq can alleviate selection bias more effectively. The pathways were drawn using ggplot2 (Wickham, 2016) in the R language. For survival analysis, TCGA PCa (PRAD) data were selected in UCSC Xena (Goldman et al., 2020) with progression free intervals to draw Kaplan-Meier plots of circRNA cognate genes. The expression level of the gene was used for survival analysis.

RNA Extraction and Real-Time Quantitative PCR (qPCR)
A total of 20 PCa tissues and 20 adjacent normal tissue samples were prepared. RNA was extracted from three 10µm FPE sections per sample. Paraffin was removed by xylene extraction followed by washing with ethanol. RNA was isolated from the sectioned tissue blocks using the purification kit, total RNA was extracted, and RNA was subjected to DNase I (Invitrogen, AM2222) treatment. The qRT-PCR was performed using the TransScript II Green One-Step qRT-PCR SuperMix kit (TransGen Biotech, AQ311-01) with 100 ng RNA as template in a 20 µL reaction volume on an ABI 7,500 real-time cycler (Qiagen).  The red points represent differentially expressed circRNAs in PCa with fold-change >2 or fold-change <0.5. The paired t-test was used to obtain p-value. (B) The fold-change of circRNAs and its cognate mRNA expression between tumor tissue and adjacent normal tissue. The fold-change of circRNAs was much higher than its mRNAs especially in hsa_circ_0006410. In both tumor as well as normal samples, the expression of hsa_circ_0003970 cognate mRNA was 0. Thus, the mRNA fold change of hsa_circ_0003970 was 0. (C) The exon composition of hsa_circ_0006410, hsa_circ_0003970, hsa_circ_0006754, hsa_circ_0005848, and hsa_circ_AKAP7. (D) Heatmap of the expression of five specific circRNAs. These five circRNAs are highly expressed in tumor tissues than in normal tissues. The value used in this figure is the expression of circRNAs, and we used "scale = row" to make the graph more consecutive.
Frontiers in Genetics | www.frontiersin.org PCR cycling was performed as follows: one cycle at 95 • C for 10 min, 95 • C for 20 s, and 40 cycles at 60 • C for 45 s. The threshold cycle for a given amplification curve during RT-PCR occurs at the point where the fluorescent signal grows beyond a specified fluorescence threshold setting.
FIGURE 3 | The network between five specific circRNAs and their predicted interactions of miRNAs. All interactions between circRNAs and miRNAs were obtained and extracted for hsa_circ_0006410, hsa_circ_0003970, hsa_circ_0006754, hsa_circ_0005848 and hsa_circ_AKAP7 to build the network.
Due to the small sample size in our sequencing control set, the p-value was not accurate for initial screening. The differentially expressed circRNAs included 32 upregulated circRNAs and 57 downregulated circRNAs (Figure 2A). According to this result, we selected the top five differentially expressed circRNAs for further analysis-included hsa_circ_0006410, hsa_circ_0003970, hsa_circ_0006754, hsa_circ_0005848, and a novel circRNA that we named hsa_circ_AKAP7 ( Figure 2C and Table 1). Because of the circRNA back-splicing feature, an mRNA may correspond to multiple circRNAs. Figure 2C shows the precise corresponding exon, enabling the identification of its component. These five circRNAs were highly expressed in tumor tissues compared with their cognate mRNAs ( Figure 2B). The Pearson correlation coefficients of circRNAs and mRNAs were consistent with those found in previous studies; circRNAs and mRNAs tend to have a high correlation . The produced circRNAs and their cognate mRNAs usually inhibit each other but Figure 2B did not show an opposite trend of circRNA and mRNA expression. This indicated that these five circRNAs have key functions in tumor tissue instead of their cognate mRNAs. The expression of these five circRNAs is depicted in Figure 2D. The prominent higher expression of these five specific circRNAs led us to further research. Therefore, we regarded circRNAs as miRNA sponges to explore the functions of these five circRNAs.

CircRNA-miRNA Network and miRNA-Related mRNAs
After obtaining circRNA and miRNA interactions from circBank and miRNDB, we built five circRNA and miRNA interaction networks (Figure 3). We found 215 circRNA-miRNA interactions in the network, and each circRNA had an average of 43 miRNA interactions. Many of these miRNAs have been reported to be associated with PCa (Table 2). Hsa_circ_0003970 and hsa_circ_0005848 interacted with miRNA-204-5p, which is a tumor suppressor that promotes apoptosis by targeting BCL2 in PCa cells (Lin et al., 2017). Hsa_circ_0005848 and hsa_circ_0006754 related miR-3160-5p is a PCa cell proliferation suppressor that targets the F-box protein (Lin et al., 2018). MiR-548 acts as an anti-oncogenic factor that inhibits the phosphoinositide three-kinase (PI3K)/AKT signaling pathway in lung cancer and is associated with high-risk Gleason scores in PCa (Shi et al., 2015). The PI3K/AKT pathway is involved in tumor immunological surveillance and immune suppression (Dituri et al., 2011). Hsa-miR-181b-2-3p and hsa-miR-96-5p were associated with the androgen receptor and Gleason score (Mekhail et al., 2014). These results further highlight the contribution of this study.
The miRNA target prediction based on the short seed sequence provided many false positive results. Thus, we filtered differentially expressed mRNAs in PCa to analyze miRNA and mRNA interactions. Only scores >90 interactions were selected to predict circRNAs function. Then, we used functional enrichment analysis to explore the function of the five circRNAs.

CircRNA-Related Pathways in Metabolism Pathways
Analysis of these five circRNA-related mRNAs showed that they were all enriched in many well-known PCa pathways ( Figure 4A). Hsa_circ_0006410, hsa_circ_0003970, hsa_circ_AKAP7, hsa_circ_0006754, and hsa_circ_0005848 were all related to the MAPK signaling pathway. MAPK signaling is an important regulator of cancer, especially PCa. It includes three cross-signaling pathways: p38, JNK, and ERK (Dhillon et al., 2007). Each pathway comprises several levels of kinases. The p38-MAPK pathway is important for the production of inflammatory cytokines and IFN-γ. It can also positively regulate Th1 differentiation instead of Th2 (Martinez et al., 2009). The JNK-MAPK pathway plays pro-inflammatory roles in macrophages, inducing M1 differentiation. Activation of the ERK-MAPK pathway favors cell differentiation into CD4 lineage and is critical for CD4 T cell polarization of Th2 because it is required for IL-4 receptor function (Alessandro et al., 2019). These regulators are significant in PCa.
Other significant pathways were the hormone-mediated signaling pathway and cellular lipid-related process, which were associated with hsa_circ_0006410, hsa_circ_0003970, hsa_circ_AKAP7, hsa_circ_0006754, and hsa_circ_0005848 (Figures 4B-F). Steroid androgen hormones play key roles in the progression and treatment of PCa. Androgen deprivation therapy (ADT) is the first-line treatment used to control cancer growth (Munkley et al., 2016). It functions by inhibiting the production of male hormone testosterone and preventing it from reaching PCa cells. ADT can cause apoptosis of PCa cells and can make them grow slowly. Studies have indicated that dietary fat intake is related to PCa development, suggesting that lipid metabolism plays a role in the carcinogenesis and progression of PCa (Tamura et al., 2009). Dysregulation of metabolism of lipids, especially sphingolipid, is a hallmark of the malignant phenotype. Increased lipid accumulation leading to changes in levels of lipid metabolic enzymes has been verified in various tumors, including PCa (Wu et al., 2014). Castration-resistant PCa (CRPC) is considered to utilize de novo lipid synthesis to produce fatty acids to obtain energy (Eidelman et al., 2017). The five circRNAs were related to both the hormone-mediated signaling pathway and the lipid-related process, indicating that they are involved in PCa regulation. Several pathways found to be closely related to PCa include the chemokine pathway, cell cycle, p53 signaling    pathway, apoptosis and transcriptional misregulation in cancer ( Figure 4A).

Survival Analysis of CircRNA Cognate Genes and qRT-PCR Validation of Five CircRNAs
To observe the effect of the five circRNAs in PCa patients, survival analysis of circRNA cognate genes in TCGA data was performed ( Figure 5). TUSC3 (hsa_circ_0006410 cognate gene), AKAP7 (hsa_circ_AKAP7 cognate gene), and RBM39 (hsa_circ_0005848 cognate gene) were all significantly associated with progressionfree survival of PCa patients, as shown by the Kaplan-Meier plot (P-value < 0.05) ( Figure 5A). Patients with high expression of TUSC3 and AKAP7 showed better overall survival. This is consistent with the fact that TUSC3 is a tumor suppressor gene (Yu et al., 2017). Meanwhile, low expression of RBM39 was found to be associated with low overall survival. RBM39 is associated with precursor messenger RNA (pre-mRNA) splicing factors, and inactivation of RBM39 causes aberrant pre-mRNA splicing. Previous studies have shown that several single amino acid substitutions in RBM39 confer resistance to the toxic effects of indisulam in cultured cancer cells and in mice with tumor xenografts (Han et al., 2017). Since the direction of differential expression varied among the five mRNAs, we think that circRNAs might act as oncogenes or tumor suppressor genes in PCa. The direction of different functional circRNAs is different in its cognate mRNAs (Figure 5). However, further experiments are required for confirming this. We used qRT-PCR to validate the five circRNAs in 20 PCa and normal samples (Figure 6). In our results, four circRNAs were significantly validated-hsa_circ_0006410, hsa_circ_0003970, hsa_circ_AKAP7, and hsa_circ_0006754. The relative expression of circRNAs indicated that these circRNAs were highly expressed in tumor tissues compared to normal tissues and validated our analysis.

DISCUSSION
Based on sequencing data of PCa tissues and adjacent normal tissues, we identified differentially expressed circRNAs in PCa. We filtered five specific highly expressed circRNAs that had never been studied in PCa before for further analysis. We found that the miRNAs and mRNA pathways related to these circRNAs were related to known metabolic pathways, such as PI3K-Akt signaling pathway, MAPK signaling pathway, and lipid metabolic process. This also confirmed the reliability of our findings. Through bioinformatics analysis, we analyzed expression levels of circRNA through linear RNA expression level adjustment using CIRCexplorer3. For these five circRNAs, the expression of circRNAs and mRNAs in tumor tissue was highly correlated, which is consistent with the results of previous studies. However, the fold changes of circRNAs expression were notably larger than those of their cognate mRNAs, suggesting that circRNAs play a role in tumor tissues. The underlying mechanism, however, is still unknown and requires further research. We predicted circRNAs function by using circRNA-miRNA-mRNA interactions and showed that they were all significantly enriched in the lipid metabolism pathway. The link between PCa development and lipid metabolism is well established, with AR intimately involved in a number of lipogenic processes. Altered lipid signatures may offer insights into metabolic reprogramming. Lipid pathway deregulation in advanced PCa is a hot research field to identify a therapeutic pathway. Several therapeutic agents, such as warfarin, atostatin, and orlistat, are known to block key processes in lipid metabolism and negatively influence PCa progression. Lipid metabolism is also activated by the PI3K-Akt signaling pathway by sterol regulatory element-binding protein 1 (SREBP1) (Edlind and Hsieh, 2014). The hsa_circ_0006410, hsa_circ_0003970, hsa_circ_0006754, hsa_circ_0005848, and hsa_circ_AKAP7 were all enriched in lipid related pathways, which shows their potential as targets.
The PI3K-Akt signaling pathway is deregulated in 42% of localized disease and 100% of advanced-stage disease in PCa. This implies that the alteration of this pathway is another factor in the development of CRPC. Gene amplifications, mutations, and changes in mRNA expression of PI3K signaling pathway are highly correlated with PCa patients (Edlind and Hsieh, 2014). Hsa_circ_0006410, hsa_circ_0003970, hsa_circ_0005848, and hsa_circ_AKAP7 were all related with this pathway. CircRNAs have been reported to activate the PI3K/Akt signaling pathway by regulating gene expression in PCa . Although the exact mechanism affecting PI3K/Akt signaling pathway is unclear, the circRNAs identified in this study also provided support for this field.
We used survival analysis and qRT-PCR to validate our findings. Survival analysis is a good indicator to assess the function of genes. Three cognate genes of these five circRNAs were significantly identified in survival analysis, alluding that their cognate genes were key genes in regulating tumor progression. Meanwhile, four circRNAs were well verified by qRT-PCR, except for hsa_circ_0005848. We inferred that this may be due to the space structure or the false positive expression of hsa_circ_0005848.
Our research was based on the bioinformatics analysis of RNA-seq between prostate tumor tissues and adjacent normal tissues. We found five specific circRNAs that were highly related to the AR signaling pathway, MAPK signaling pathway, hormone-mediated signaling pathway, and cellular lipid metabolic process. Furthermore, survival analysis and qRT-PCR validation also verified that the circRNAs were closely related to tumor progression of PCa. These five circRNAs can provide new solutions for research in this field.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Genome Sequence Archive for Human (https://bigd.big.ac.cn/bioproject/browse/ PRJCA003890) (accession: PRJCA003890).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Beijing Hospital Ethics Committee. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
LZ and WZ conceived the study and wrote the manuscript. LZ, FS, and YZ designed the detail analysis pipeline. LZ and FS did the bioinformatics analysis. WZ, HL, XT, and SX performed the experiments. MW and LW participated in revising the manuscript. All authors read and approved the final manuscript.