A systematic method introduced a common lncRNA-miRNA-mRNA network in the different stages of prostate cancer

Introduction The present study aimed to investigate the interaction of the common lncRNA-miRNA-mRNA network involved in signaling pathways in different stages of prostate cancer (PCa) by using bioinformatics and experimental methods. Methods Seventy subjects included sixty PCa patients in Local, Locally Advanced, Biochemical Relapse, Metastatic, and Benign stages, and ten healthy subjects were entered into the current study. The mRNAs with significant expression differences were first found using the GEO database. The candidate hub genes were then identified by analyzing Cytohubba and MCODE software. Cytoscape, GO Term, and KEGG software determined hub genes and critical pathways. The expression of candidate lncRNAs, miRNAs, and mRNAs was then assessed using Real-Time PCR and ELISA techniques. Results 4 lncRNAs, 5 miRNAs, and 15 common target genes were detected in PCa patients compared with the healthy group. Unlike the tumor suppressors, the expression levels of common onco-lncRNAs, oncomiRNAs, and oncogenes showed a considerable increase in patients with advanced stages; Biochemical Relapse and Metastatic, in comparison to the primary stages; Local and Locally Advanced. Additionally, their expression levels significantly increased with a higher Gleason score than a lower one. Conclusion Identifying a common lncRNA-miRNA-mRNA network associated with prostate cancer may be clinically valuable as potential predictive biomarkers. They can also serve as novel therapeutic targets for PCa patients.


Introduction
Prostate cancer (PCa) is the second most diagnosed cancer in males worldwide and causes more than 350,000 deaths annually (1). Rectal examinations and prostate-specific antigen (PSA) levels are the first steps to check for PCa. PSA evaluation is beneficial for early diagnosis, but insufficient for disease grouping. Because a high value might be present in a person without cancer, and a low value can be present in someone with cancer, PSA values are difficult to interpret (2,3). Prostate biopsies can also cause considerable infectious complications, psychological damage, and financial costs. So, the development of novel molecular biomarkers for PCa detection has allowed for better therapeutic target evaluations and prognostic assessments (4). Indeed, it is essential to explore molecular mechanisms, such as non-coding RNAs (lncRNA and miRNA), and identify novel targets for developing predictive, prognostic, and therapeutic goals for PCa.
Non-coding RNAs are specific, accurate, and accessible noninvasively, which makes them appealing for use as biomarkers for identifying tumor presence and subtypes. They control numerous biological processes, including cell proliferation, differentiation, and apoptosis (5). They are also recognized as critical regulators in various networks as both tumor repressors and oncogenic. In this case, lncRNA and miRNA interaction by integrating mRNA can provide opportunities for further experimental studies and improved biomarker predictions for developing novel diagnostic approaches. For example, miR-21 is increased in PCa (6,7), and its overexpression can activate TGF-b and Hedgehog signaling pathways, promoting invasion and metastasis (8). Therefore, detecting non-coding RNAs, actively released from tumor cells into body fluids, makes them suitable as diagnostic and therapeutic biomarkers (9). Herein, we aim to investigate the interactions between common lncRNA-miRNA-mRNA networks in patients with prostate cancer. We hypothesize that a better understanding of the underlying mechanisms involved in this network and identifying commonly expressed biomarkers at different stages of PCa can lead to developing predictive, prognostic, and therapeutic goals.

The data collection of PCa datasets
Expression profiles of mRNAs (GSE16120, GSE119195, GSE62872, GSE30994, GSE69223, GSE55945, GSE68555, GSE51066, GSE134499, GSE36135, GSE3325, GSE116918, GSE70770, and GSE35988) from PCa with different tumor stages (local, locally-advanced, biochemical relapse, and metastatic), normal samples, and benign tumor samples were extracted from t h e G e n e E x p r e s s i o n O m n i b u s ( G E O ) ( h t t p s : / / www.ncbi.nlm.nih.gov/geo/) database. Gene expression information for Agilent-012097 Human 1A Microarray, Affymetrix Human Gene 1.0 ST Array, Agilent-014850 Whole Human Genome Microarray, Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix Human Genome U95A Array, Affymetrix Genome-Wide Human SNP 6.0 Array, Almac Diagnostics Prostate Disease-Specific Array, Agilent-039494 SurePrint G3 Human Gene Expression v2 8x60K Microarray, and Agilent-014698 Human Genome CGH Microarray 105A were obtained from prostate cancer patients with different tumor stages. GEPIA2 and cBioPortal analyze differentially expressed genes (DEGs) from the TCGA (10). The expression profile pattern of PCa patients was compared to that of healthy individuals to identify DEGs. The GSEs' data were downloaded using the GEOquery R package. We analyzed the selected lncRNAs, miRNAs, and genes with P-value < 0.05 and |LogFC| > 1 in the datasets as differentially expressed genes (DEGs), differentially expressed miRNAs (DEMs), and differentially expressed lncRNAs (DELs) (11).

Enrichment analysis for DEGs
Cytoscape software uses the ClueGo tool to identify gene ontologies (11). It is accepted that p < 0.05 is statistically significant in GO analysis. Moreover, the Reactome and WikiPathway databases were used to identify signaling pathways (10). The Enrichr software performed pathway enrichment analyses of the DEGs.

Experimental design and sampling
Eighty PCa patients and ten healthy subjects were entered into the current study from September 2019 to September 2020. The eligible cases were categorized into six groups: Local (L, N = 22); those whose tumor has not spread outside the prostate capsule, Locally Advanced (LA, N = 9); those whose tumor has invaded the seminal vesicles and other pelvic organs, Biochemical Relapse (BR, N = 11); those in which PSA increased again after the treatment period, Metastatic (MET, N = 8); those whose cancer tissue has spread in secondary areas such as lungs and bones, and Benign Prostate Hypertrophy (BPH, N = 10) patients, and healthy samples based on pathologic reports.
Biosystems; Thermo Fisher Scientific, Inc.) (10). For normalization, U6 and B-actin were used as internal controls. The fold changes of candidate lncRNAs, miRNAs, and mRNAs were calculated by the equation −DCT (11). The Supplementary File (S1) shows the primers' list.

Statistical tests
The statistical tests were conducted with GraphPad Prism 7.04 (San Diego, CA). The K-S test estimates data normality. ANOVA and t-tests were applied to analyze multiple and two groups. The Mann-Whitney U test evaluated nonparametric data. REST analyzed the PCR data. Regression analysis was performed to avoid bias in confounding variables. P < 0.05 was considered as a level of significance.

Identification of hub genes at different stages of prostate cancer
A total of 23 genes (9 up-regulation and 14 down-regulation) in the local group, 30 genes (17 up-regulation and 13 downregulation) in the locally-advanced group, 68 genes (17 upregulation and 51 down-regulation) in the metastatic group, 44 genes (9 up-regulation and 35 down-regulation) in the biochemical relapse group, and 228 genes (53 up-regulation and 175 down-regulation) in the BPH group (Figures 2A, B) were analyzed from GEO and TCGA datasets to identify differentially expressed genes in prostate cancer and healthy samples (Supplementary S2). A Venn diagram and Upset plot of DEGs at different stages of the disease were drawn by calculating and drawing custom Venn diagrams using the online tools and FunRich_3.1.3 software and upset plot using Up Set R package ( Figure 2). Furthermore, differentially expressed miRNAs (Table 1) and lncRNAs ( Table 2) were determined at various stages of prostate cancer.

Identified hub genes by GO term
The GO analysis revealed the enriched pathways identified in the DEGs (Figure 3). Moreover, the signaling pathways involved in the hub genes from the Reactome, KEGG, and Wikipathway databases were listed in Figure 4. Databases' common pathways involve focal adhesion and cell proliferation, such as PI3K/AKT/ mTOR ( Figure 4).

Validation of hub genes by PPI network
STRING software was used to construct and predict PPI network information ( Figure 5A). Candidate hub genes were identified by analyzing Cytohubba and MCODE software ( Figure 5B). Our results showed 4 lncRNAs, 5 miRNAs, and 15 mRNAs from a regulatory lncRNA-miRNA-mRNA network in PCa patients ( Figure 5B). Besides, a heatmap map of candidate A B D C FIGURE 2 A Venn diagram and UpSet plot of the differently expressed mRNAs between the GEO and TCGA datasets (A, B) 23 genes (9 up-regulation and 14 down-regulation) in the local group, 30 genes (17 up-regulation and 13 down-regulation) in the locally advanced group, 68 genes (17 up-regulation and 51 down-regulation) in the metastatic group, 44 genes (9 up-regulation and 35 down-regulation) in the biochemical relapse group, and 228 genes (53 up-regulation and 175 down-regulation) in the BPH group were differentially expressed compared to healthy samples. (C, D) the UpSet plot of overlapped differential expression genes among candidate groups.

Patients' characteristics
The clinicopathological characteristics of the participants are summarized in Table 3. There were no significant differences in age, BMI, and prostate volume in the different groups. The most considerable tumor volume (cm 3 ) belonged to the metastatic group (58.3 ± 7.7), and most prostate cancer patients had adenocarcinoma in pathology assessment. Except for the biochemical relapse group (due to surgery and complete prostate resection), serum PSA levels increased significantly with disease progression in other groups. All patients in the biochemical relapse group were given hormone therapy and a radical prostatectomy. All metastatic patients received hormone therapy and palliative radiotherapy. In the other groups, most received hormone therapy alone.

The correlation of clinicopathological status with commonly expressed biomarkers
The expression level of onco-lncRNAs, including PCAT19, MALAT1, and NEAT1, and the tumor suppressor lncRNA, CASC2, was increased and decreased, respectively, in PCa patients compared to the healthy group (P < 0.05) ( Table 4). Interestingly, the altered expression levels of lncRNAs were significantly associated with increased tumor stages (Table 4) and Gleason score (Table 5).
Except for miRNA-195, the expression levels of the oncomiRs, including miRNA-491, miRNA-182, miRNA-149, and miRNA-200c, showed a significant increase in PCa patients compared to the healthy group (P < 0.05) ( Table 6). Interestingly, except for miRNA-195, the expression levels of the mentioned other miRNAs were significantly associated with increased tumor stages ( Table 6) and Gleason score ( Table 7).
A heatmap of the expression of candidate lncRNAs, miRNAs, and mRNAs with clinicopathological Gleason scores was also generated by CIMminer. We showed that the high expression of the onco-lncRNAs, oncomiRNAs, and oncogenes, and the low expression of the tumor suppressors have led to increases in Gleason scores (Figure 7).
Moreover, we compared candidate biomarkers' profiles according to age. Our analysis showed no significant difference between lncRNA, miRNA, and mRNA expression with participants' age.

Discussion
The present study aimed to investigate the underlying pathophysiological pathways involved in different stages of prostate cancer tumors. Our results revealed a significant reduction in tumor suppressor lncRNAs, miRNAs, and mRNAs alongside a considerable increase in onco-lncRNAs, oncomiRNAs, and oncogenes at different stages of prostate cancer. Likewise, these biomarkers with a higher Gleason score were significantly increased compared to tumors with lower Gleason scores. Furthermore, increased PSA levels were associated with advanced disease stages in prostate cancer. Thus, the interactions between lncRNAs, miRNAs, and mRNAs could contribute to finding novel biomarkers for assessing treatment response in PCa patients. Therefore, a better understanding of miRNA-lncRNA-mRNA signaling pathways may identify novel therapeutic targets for PCa patients.  Many patients suffering from advanced PCa will ultimately relapse after initial chemotherapy and hormone therapy-based regimens (12). Numerous genetic alterations concerning cell death and survival signaling processes may be responsible for PCa induction in this setting. Over the past decade, several investigations have been conducted concerning prostate cancer-related proteinencoding genes (13). However, only 2% of human transcriptomes are protein-coding RNAs (14). In order to better understand the underlying mechanisms of PCa initiation and progression, non-coding RNAs should be investigated (15,16). Previous studies demonstrated that ceRNAs play pivotal roles in Pca progression and metastasis and even resistance to treatment (17). Recently, Liu et al. (18) showed that some miRNAs, such as miRNA-141, could reduce cancer growth and metastasis through several mechanisms (18). They also reported that miRNA-34a could directly prevent Pca regeneration and metastasis by suppressing the CD44 gene (19).  Pathway-based enrichment analysis of high-expressed hub genes and low-expressed hub genes at different stages of prostate cancer The top 10 functional pathways were associated with DEGs through WikiPathway, Reactome, and KEGG analyses with a p-value of less than 0.05. tissues, while miRNAs -31, -96, and -205 exhibited a significant correlation with Gleason score (20). Furthermore, Borkowetz et al. (21) revealed that PSA's diagnostic potential will significantly increase in combination with miRNA-16 and miRNA-195 as diagnostic biomarkers (21). Similarly, our results showed that the expression levels of miRNAs -491, -182, -149, and -200c were significantly increased in PCa patients compared to the healthy group which was significantly associated with advanced tumor stages.
Moreover, Aiello et al. (22) demonstrated that lncRNAs, such as HOTAIR and MALAT1, were associated with both ERa/ERb steroid receptors in prostate cancer patients. Inhibition of MALAT1 led to better treatment response (22). Another study indicated that MALAT1 overexpression was associated with poor prognostic markers, including high Gleason score, metastasis stage, and serum PSA > 20 ng/ml. In concordance with our results, its expression was significantly increased in patients who experienced disease relapse (23). Ramnarine et al. (2018) revealed that several lncRNAs were involved in developing PCa (24). Zhang et al. (25) indicated that cell proliferation in hormone-refractory PCa was stimulated by lncRNAs (25). Similarly, our results demonstrated a significant increase in onco-lncRNAs, such as PCAT19, MALAT1, and NEAT1, and a decreased expression level of CASC2 in PCa patients. Additionally, several studies investigated the interactions between lncRNAs, miRNAs, and mRNAs in various malignancies. A study indicated that elevated miRNA-182-5p was associated with cancer cells' mitotic arrest through MALAT1 expression (26). Wang et al. (27) revealed an inverse correlation of expression between miRNA-195-5p and FGF in human cancer tissues (27). Regarding the ceRNA hypothesis in PCa, Xu et al. (28) revealed a ceRNA network of 63 prostate cancer-specific lncRNAs, 13 miRNAs, and 18 mRNAs. They also established a predictive model for overall survival using HOXB5, GPC2, PGA5, and AMBN (28). Besides, Jiang et al. (29) constructed a ceRNA network that included 23 lncRNAs, 6 miRNAs, and 2 mRNAs differentially expressed in PCa. However, only 3 lncRNAs were significantly associated with overall survival (29).
Basically, androgen/androgen receptor (AR) signaling is crucial for both normal prostate development and tumorigenesis which regulates a series of cancer biological processes including cell proliferation and metastasis (30, 31). Interestingly, several studies confirmed that AR-bound ceRNAs are required for AR interaction (32, 33). AR could directly target miRNAs and lncRNAs in PCa cells, which were recognized by microarray or RNA sequencing (34,35). These studies strongly suggested that several ceRNAs containing specific sequences could bind to AR and other steroid receptors to stop their transcriptional activities (36). In this case, androgen deprivation therapy (ADT) is the standard management for patients with recurrences after primary treatment (37). Moreover, most primary cancers acquire ADT resistance and progress to castration-resistant PCa (38). Therefore, there is a critical need to clarify the molecular mechanisms contributing to the development of AR-dependent patients for emerging substitute  therapeutic opportunities in advanced PCa. AR's role and interaction with the network are poorly described. However, recent studies provide substantial evidence for this crosstalk in PCa. Several lncRNA-miRNA-mRNA networks are aberrantly expressed with significant contributions to PCa cell initiation and progression (39, 40). Application of their inhibitors can limit offtarget or non-specific effects by avoiding targeting ceRNA complexes. Ostling et al. (2011) showed 78 ceRNA inhibitors to significantly modulate AR reporter activity (41). Some of them may activate or repress mRNA translations (42,43). Further studies are required to disentangle these conflicting outcomes.

Clinical implications and limitations
Herein, we investigated commonly expressed lncRNA, miRNA, and target genes at different stages of prostate cancer. We hypothesized that a better understanding of the lncRNA-miRNA-mRNA network might develop novel therapeutic options for PCa patients. Finding therapeutic targets expressed at all prostate cancer stages has several advantages. However, there are some advantages and disadvantages. We used a large number of data from different sources to minimize off-target findings. Our patient-derived data showed network overlap. Cancer patients at any stage of the disease will benefit from their medicinal properties. While there is an association between disease stage and biomarker expression levels, these may be considered predictive biomarkers for treatment response. Since we demonstrated a significant difference in the expression levels of indicated lncRNAs and miRNAs in normal Tumor suppressor lncRNAs CASC2 -1.5 ± 0.7 -9.8 ± 1.9 * -17.9 ± 0.9 # Data were calculated based on mean ± SD and -DCT. GS, Gleason score. *P <0.05 compared to the Healthy group. # P <0.05 compared to the GS ≤ 6.      tissue and benign prostate hyperplasia, these biomarkers can act as a potential diagnostic biomarker in combination with PSA to increase diagnostic accuracy. Taken together, we introduced several promising diagnostic biomarkers. Further studies should identify this issue as a second-level test for patients. Some limitations merit consideration. First, due to the limited sample size, we could not evaluate the diagnostic accuracy of indicated lncRNAs and miRNAs as predictive tools. Therefore, future studies are warranted in this regard. Second, bioinformatic analysis predicted an interaction between lncRNAs and miRNAs. In vitro assessments directly evaluate the potential interactions between each lncRNA and its associated miRNA. Third, assessing the prognostic significance of the indicated lncRNAs and miRNAs regarding other clinical outcomes of patients and survival is warranted.

Conclusion
We have identified a common lncRNA-miRNA-mRNA network associated with prostate cancer, which may be clinically valuable as potential predictive biomarkers. These biomarkers might also serve as novel therapeutic targets or potential prognostic indicators for PCa patients.

Data availability statement
The data analyzed in this study is subject to the following licenses/restrictions: The data supporting this study's findings are available from the corresponding author, AA, upon reasonable request. Requests to access these datasets should be directed to AA, aalizadeh@sina.tums.ac.ir.

Ethics statement
The studies involving human participants were reviewed and approved by Iran University of Medical Sciences (NO: IR.IUMS.REC.1399.1348). The patients/participants provided their written informed consent to participate in this study.

Author contributions
GV: study concept and design and manuscript revision. SK: sample processing, data analysis, and manuscript preparation.
AMA: study concept and design, data analysis, and manuscript revision. MY: sample collection. MM: manuscript preparation. TR and MaB: sample processing. MR: manuscript preparation. EE, MoB, and AK: sample collection. All authors contributed to the article and approved the submitted version.

Funding
Elite Researcher Grant Committee supported the research reported in this publication under award number (Grant Number: 995719) from the National Institute for Medical Research and Development (NIMAD), Tehran, Iran. This study was also funded by the Iran University of Medical Sciences (Grant Number: 19375). None of the funding sources had any role in the study design, data collection, analysis, interpretation, or the decision to submit the article for publication.