Identication of Androgen Receptor Variant 7-related RNAs Affecting Abiraterone Ecacy in Castration-resistant Prostate Cancer Treatment by RNA-sequencing

Background: This study aimed to identify androgen receptor variant 7 (AR-V7)-related RNAs affecting Abiraterone treatment of castration-resistant prostate cancer (CRPC) using RNA-sequencing. Methods: To identify AR-V7-related RNAs affecting Abiraterone treatment of CRPC, a series of in vitro experiments were employed, including cell proliferation assay, cell apoptosis assay, Western blot analysis, and Real-time quantitative reverse transcription PCR (qRT-PCR). After RNA-sequencing, the differentially expressed (DE) mRNAs and DE long non-coding RNAs (lncRNAs) were screened and the lncRNA-mRNA pairs were identied. In addition, enrichment and Protein-protein interaction (PPI) network analyses were performed. Finally, the mRNA-miRNA-lncRNA competing endogenous RNAs (ceRNA) network was built, along with survival analysis. Results: A total of 1,387 AR-V7-related RNAs affecting Abiraterone treatment of CRPC were identied. The enrichment analysis showed that the target genes of DEmRNAs and DElncRNAs were primarily involved in cancer-related pathways, including the ErbB signaling pathway, pathways in cancer, and basal cell carcinoma. In addition, notch receptor 1 (NOTCH1), ionotropic NMDA glutamate receptor type subunit 1 (GRIN1), U-box containing protein 1 (STUB1), and mitogen-activated protein kinase 7 (MAP2K7) with high degrees in the PPI network. Moreover, MAP2K7 was regulated by hsa-miR-6825-5p, which in turn was regulated by MAFG-AS1 lncRNA in the ceRNA network. The survival analysis revealed that a total of four lncRNAs, including MAFG-AS1, and 17 mRNAs, including high muscle blind-like splicing regulator 2 (MBNL2), were associated with disease-free survival (RFS). Among them, only MBNL2 overexpression correlated with good survival outcome. The NOTCH1, GRIN1, STUB1, MBNL2, and ErbB signaling pathways are likely related to the ecacy of Abiraterone In this study, we rst investigated the effect of Abiraterone administration on the AR-V7 expression level of androgen-dependent cell lines in vitro. Moreover, AR-V7-related RNAs affecting Abiraterone treatment of CRPC were identied through bioinformatics analysis based on RNA sequencing. Here, we explored the mechanism of Abiraterone effecting AR-V7 and AR-V7 affecting therapy ecacy of Abiraterone, and provide new options for Abiraterone treatment of CRPC.

Prostate cancer (PCa), the most common malignancy in men, had a global incidence rate that was second only to lung cancer in 2018 [1] . With changes in lifestyle and diet, as well as an aging population, the incidence of PCa in China is increasing every year [2] . Androgen deprivation therapy (ADT) achieves prostate-speci c antigen (PSA) rate decline or disease control in patients with advanced PCa [3] . However, almost all PCa patients progress to castration-resistant prostate cancer (CRPC) after a period of treatment. Thus, it is urgent to understand the speci c molecular mechanism of CRPC progression and search for potential therapeutic targets for CRPC.
The androgen receptor (AR) plays a key role in the progression of CRPC. CRPC invariably develops due to the aberrant reactivation of the androgen/AR signaling axis [4] . In addition, the expression levels of AR splice variants (AR-Vs) in CRPC were generally higher than those in hormone-naive PCa, especially AR-V7 (also known as AR3) and AR v567es [5] . Numerous studies show that AR-V7 is associated with the development of CRPC [6][7][8] . Observations of primary tumor tissues before and after castration resistance indicate that AR-V7 expression increases with the outgrowth of CRPC tumors [9,10] . Abiraterone, a secondgeneration anti-androgen drug, is approved for CRPC treatment. Resistance to Abiraterone occurs frequently, although it is initially effective for CRPC treatment [11] . Evidence from experimental and clinical studies illustrate that AR-V7 plays a vital role in the promotion and progression of CRPC during ADT and the induction of Abiraterone resistance [12,13] , and in some CRPC patients started with Abiraterone didn't respond, but after several months continuous treatment they could respond to it [12] . However, the mechanism of AR-V7 affecting Abiraterone e cacy and Abiraterone affecting AR-V7 in CRPC treatment is unknown. Notably, noncoding RNA transcripts, such as miRNAs and lncRNAs, play signi cant roles in the molecular mechanism of cancers [14,15] . In this study, we rst investigated the effect of Abiraterone administration on the AR-V7 expression level of androgen-dependent cell lines in vitro. Moreover, AR-V7related RNAs affecting Abiraterone treatment of CRPC were identi ed through bioinformatics analysis based on RNA sequencing. Here, we explored the mechanism of Abiraterone effecting AR-V7 and AR-V7 affecting therapy e cacy of Abiraterone, and provide new options for Abiraterone treatment of CRPC.
The cell apoptosis assay The 22RV1 cells (7 10 5 ) and LNCaP cells (2 10 5 ) were treated with 175 µM Abiraterone, and cell proliferation studies were performed at 24, 48 and 72 h using an annexin V-uorescein isothiocyanate (FITC)-propidium iodide (PI) cell apoptosis kit (Invitrogen, Thermo Fisher Scienti c, Inc.). The cells were washed twice with PBS buffer (pH 7.4) then resuspended in the binding buffer. Moreover, 5 µl PI and 5 µl Annexin V-FITC were mixed with the cells. The mixtures were incubated at 25°C for 15 min in the dark, and then the combinations were measured utilizing FACScan ow cytometry (FACSCalibur, BD).
The FASTQ sequence data were acquired from the RNA-sequencing data using Base Calling version 0.11.4. In addition, the FASTQ data quality control was performed utilizing the READ QC tool in FastQC version 0.11.4. Cutadapt version 1.16 was used for raw data trimming. Also, to acquire clean reads, the low-quality reads, including sequences with a quality score < 15, adaptor sequences, and sequences with an N number of raw reads > 1 were removed.

Expression level and Pearson correlation coe cient analysis of samples
The cor function in R language was employed to calculate the Pearson correlation coe cient (P) between each pair of samples. The prcomp algorithm in R language was utilized to conduct PCa of samples. Finally, the PCa diagram was performed using the Ggfortify package version 0.4.5.

Differential expression analysis
The voom function in the limma package [19] was employed to standardize the raw counts. The differentially expressed (DE)mRNAs and DE Long non-coding RNAs (lncRNAs) between the JY-48 h group and the control group, or the JY-72 h group and the control group were screened utilizing the limma package with the threshold of adj.P < 0.05 and a |log fold change (FC)| > 1. The Benjamini and Hochberg (BH) method was used to adjust the P-value for multiple comparisons. In addition, to obtain the AR-V7related RNAs affecting Abiraterone treatment of CRPC, the DEmRNAs and DElncRNAs screened from the JY-48 h group, and the control group were intersected with DEmRNAs and DElncRNAs screened from the JY72 h group and the control group, and the overlapping DEmRNAs and DElncRNAs were obtained. Then, the overlapping DEmRNAs and DElncRNAs were removed from the DEmRNAs, and DElncRNAs screened from the JY-48 h group and the control group, and the remaining DEmRNAs and DElncRNAs were used for subsequent analyses. Moreover, the featured lncRNAs recorded in the LncBook database [20] were intersected with the other DElncRNAs, and the overlapping lncRNAs were obtained for subsequent analyses.
DEmRNA-DElncRNA co-expression network construction The P of DEmRNAs and DElncRNAs was calculated based on mRNA and lncRNA data-matched samples. In addition, the BH method was performed to adjust the P values, and the DEmRNAs-DElncRNAs pairs were screened with the cutoff value of r > 0.95 and adj.P < 0.01. Then, the DEmRNA-DElncRNA coexpression networks were constructed utilizing Cytoscape software [21] . Moreover, "prostate cancer" was used as a keyword to search disease-related lncRNAs based on the LncBook database, and the diseaserelated lncRNAs were marked in the DEmRNA-DElncRNA co-expression networks. Also, "prostate neoplasms, castration-resistant" were used as keywords to search disease-related RNAs based on the CTD database [22] (http://ctdbase.org/) and the GeneCards [23] (https://www.genecards.org/) database, and the disease-related RNAs were marked in the DEmRNA-DElncRNA co-expression networks.

Functional enrichment analysis
Both Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted on the target genes of lncRNAs utilizing the Metascape online tool [24] based on the DEmRNAs-DElncRNAs pairs. The parameters were set to min overlap = 3, P = 0.05, min enrichment = 1.5. Based on the mRNAs of the DEmRNAs-DElncRNAs pairs, the DAVID tool [25] was used to perform GO and KEGG pathway enrichment analyses with a cutoff value of P < 0.05 and count ≥ 2.

Protein-protein interaction network building
The Search Tool for the Retrieval of Interacting Genes (STRING) database [26] was performed to analyze the interactions between protein and protein encoded by mRNAs in the DEmRNA-DElncRNA co-expression network. The PPI score was set as 0.4 (referred to as the median con dence). Afterward, the PPI network was built using Cytoscape software. The CytoNCA plugin [27] was used to analyze the degree of nodes in the network, and the parameter was set to "without weight." lncRNA-miRNA-mRNA network construction According to the mRNAs in the DEmRNA-DElncRNA co-expression network, the miRTarBase database [28] was utilized to predict microRNAs (miRNAs), and the miRNA-mRNA pairs were obtained. The DIANA-LncBase database [29] was employed to predict miRNAs based on lncRNAs in the DEmRNA-DElncRNA coexpression network, and the lncRNA-miRNA pairs were obtained. The lncRNA-miRNA-mRNA pairs were further identi ed according to the positive co-expression relationships between mRNA and lncRNA (r > 0.95), and the competing endogenous RNAs (ceRNA) network was built.

Survival analysis
Based on the lncRNAs and mRNAs in the ceRNA network, the survival analysis was performed utilizing the GEPIA2 online tool [30] (http://gepia2.cancer-pku.cn/#survival), and the statistical signi cance level was de ned as P < 0.05.

Statistical Analysis
The data from each assay were shown as the mean ± standard deviation. One-way analysis of variance (ANOVA), t-tests, and Newman-Keuls multiple comparison tests were performed to assess the statistical signi cance between groups. All data were analyzed using SPSS 17.0 software and GraphPad Prism. A P < 0.05 was regarded as signi cant.

Results
Abiraterone inhibited cell growth Drug resistance and sensitivity are conventionally quanti ed using half-inhibitory concentration values (IC50). To explore the effects of Abiraterone on survival, the cells were separately treated with different concentrations of Abiraterone. As shown in Figure 1A, four cell lines were sensitive to Abiraterone and the inhibition rate increased with increasing Abiraterone concentrations. The IC50 of all cell lines were unchanged at 24 h, while the IC50 of PC-3 and LNCaP cell lines decreased to half of the initial IC50 at 48 h ( Figure 1B). The IC50 for LNCaP cells was 344.9 mM. Therefore, half of the LNCaP cell IC50 (172.45 mM) was applied in subsequent experiments.
AR-V7 protein expression in 22RV1 cells is higher than VCaP, LNCaP, or PC-3 cells AR-V7 has been revealed to be a potential clinical marker for patients with CRPC [31] . To identify the suitable cell line for subsequent analysis, AR-V7 protein expression levels in untreated cell lines (LNCaP, PC-3, VCaP, or 22RV1) were analyzed using western blot analysis. As illustrated in Figure 2, the AR-V7 protein expression level in 22RV1 cells was higher than in VCaP, LNCaP, and PC-3 cells. Therefore, LNCaP and 22RV1 cells were chosen for subsequent analysis because VCaP cells grew slowly.

Abiraterone effectively reduced AR-V7 expression in treated cells
Though some research indicated that AR-V7 status is related to the therapy e cacy of Abiraterone [12,13] , but few papers detected the in uence of Abiraterone to AR-V7 expression. In our study, to assess the effect of Abiraterone on AR-V7 expression, qRT-PCR and western blot analyses were conducted. The expression level of AR-V7 mRNA in the Abiraterone-treated 22RV1 cell group was lower than in the control group at 48 h (P < 0.01). In comparison, there was no signi cant difference between the Abiraterone and the control groups at 72 h (P > 0.05). AR-V7 mRNA expression in LNCaP cells following Abiraterone treatment was decreased compared with that in the control group at 48 and 72 h (P < 0.01) ( Figure 3A). In addition, the results of western blot analysis further indicated that Abiraterone administration reduced AR-V7 expression levels ( Figure 3B).

Abiraterone resistance in 22RV1 cells 48 and 72 h post-treatment
As shown in Figure 4A, the results showed that cell apoptosis signi cantly increased in 22RV1 cells following 24 h Abiraterone treatment (P < 0.01), while the Abiraterone group and the control group showed no signi cant differences (P > 0.05) after 48 or 72 h. The proportion of apoptotic LNCaP cells in the Abiraterone group was signi cantly increased compared with that in the control group at 24, 48, and 72 h (P < 0.01) ( Figure 4B). The results revealed that resistance to Abiraterone occurred in 22RV1 cells at 48 and 72 h, although Abiraterone is effective at 24 h. Thus, 22RV1 cells were used for RNA-sequencing.

Expression level analysis and PCA of samples
After quality control and ltering, a total of 24,043 RNAs, including 16,457 mRNAs and 7,586 lncRNAs, were identi ed in 22RV1 cells. The cor function in R language was employed to calculate the Pearson correlation coe cient (P) between each pair of samples. The correlation heatmap of each pair of samples is shown in Figure 5A. In addition, the control, JY-48 h, and JY-72 h samples were entirely separated, meaning that the expression patterns of these samples could be used to ultimately distinguish between the control, JY-48 h, and JY-72 h groups ( Figure 5B).

Functional enrichment analyses of lncRNA target genes and mRNAs
To further elucidate the functional role of the DElncRNAs and DEmRNAs in CRPC, we used the Metascape online tool and the DAVID tool to perform functional enrichment and KEGG pathway analyses of the target genes of lncRNAs and DEmRNAs. Responses to hypoxia, decreased oxygen levels, and oxygen levels were signi cantly enriched GO terms of the DElncRNAs target genes ( Figure 8A). The ErbB signaling pathway (hsa04012), pathways in cancer (hsa05200), and endocrine resistance (hsa01522) were the primary enriched KEGG pathways of DElncRNA target genes ( Figure 8B). The upregulated DEmRNAs were mainly enriched in 4 GO terms [including lactation, protein autophosphorylation, protein phosphorylation, and embryonic pattern speci cation] and the 0 KEGG pathway ( Figure 8C). In addition, the downregulated DEmRNAs were mainly involved in 27 GO terms and 2 KEGG pathways (hsa04012: ErbB signaling pathway, hsa05217: Basal cell carcinoma) ( Figure 8D).

PPI network construction
To analyze the interactions between protein and protein encoded by mRNAs in the DEmRNA-DElncRNA co-expression network, the PPI network was built. In total, 379 PPI pairs were acquired from the PPI network ( Figure 9). In this PPI network, the notch receptor 1 (NOTCH1), the ionotropic NMDA glutamate receptor type subunit 1 (GRIN1), U-box containing protein 1 (STUB1), and mitogen-activated protein kinase 7 (MAP2K7) with high degrees from a local network and may be the key proteins in the development of CRPC (Supplementary le 2). In addition, the results showed that NOTCH1, STUB1, and the CCAAT enhancer-binding protein alpha (CEBPA) were related to CRPC.

Discussion
Previous studies have established that AR-V7 plays critical roles in the progression of CRPC. For instance, Chen et al. suggested that HoxB13 might serve as a therapeutic target for AR-V7-driven prostate tumor treatment [32] . Scher et al. validated that nuclear expression of AR-V7 protein in circulating tumor cells in men with metastatic CRPC could serve as a treatment-speci c biomarker [6] . However, the regulatory mechanism of AR-V7 affecting Abiraterone treatment of CRPC is unknown. In this study, we found that Abiraterone administration signi cantly reduced AR-V7 expression levels in vitro, this may explain why some CRPC patients started with Abiraterone didn't respond, but after several months continuous treatment they could respond to it for AR-V7 expression is reduced by lone term Abiraterone treatment.
Therefore, AR-V7-related RNAs affecting Abiraterone treatment of CRPC were identi ed using bioinformatics analysis. The target genes of DEmRNAs and DElncRNAs were primarily involved in cancerrelated pathways, including the ErbB signaling pathway, pathways in cancer, and basal cell carcinoma. In addition, NOTCH1, GRIN1, STUB1, and MAP2K7 were found with high degree in the PPI network.
Moreover, MAP2K7 was regulated by hsa-miR-6825-5p, and hsa-miR-6825-5p was regulated by lncRNA MAFG-AS1 in the ceRNA network. The survival analysis showed that a total of four lncRNAs (including MAFG-AS1) and 17 mRNAs (including MBNL2) were associated with RFS. Among them, only MBNL2 overexpression correlated with a good survival outcome.
The Notch family of receptors determines the fate of cells in many organ systems, including the prostate [33] . Rice et al. have suggested that the loss of NOTCH1 in aggressive PCa cells decreases proliferation, invasion, and tumorsphere formation [34] . Kwon et al. shown that Notch1 promotes cell survival and proliferation in rat luminal prostate cells through the activation of the pro-survival NF-κB pathway [35] . STUB1 is a type of accessory chaperone protein and an E3 ubiquitin ligase. STUB1 connects the polypeptide binding activity of HSP70 to the ubiquitin-proteasome system. Liu et al. considered that inhibiting the HSP70 targeting protease stabilization pathway may be a valuable strategy for overcoming anti androgen resistance in the next generation of CRPC patients [36] .
Nevertheless, the role of GRIN1 in PCa development has not been explored. Herein, NOTCH1, GRIN1, STUB1 exhibited high degrees in the PPI network, and the survival analysis showed that MBNL2 overexpression revealed a good survival outcome. Bii et al. illustrated that MBNL2 was identi ed as a candidate gene in PCa progression [37] . Zhu et al. demonstrated that miR-636 promotes cell invasion and migration, which might promote bone metastasis of PCa by targeting MBNL2 [38] . Taken together, we speculate that NOTCH1, GRIN1, STUB1, and MBNL2 contribute to the progression of CRPC.
The lncRNA MAFG-AS1 is reportedly related to the progression of cancer. For instance, MAFG-AS1 lncRNA is considered a novel oncogenic lncRNA in the progression of colorectal cancer through its regulation of miR-147b/NDUFA4 [39] . Jia et al. demonstrated that overexpression of MAFG-AS1 lncRNA primarily downregulated miR-339-5p levels in non-small-cell carcinoma cells, and MMP15 is a miR-339-5p target [40] . In addition, Hu et al. illustrated that MAP2K7 might be a promising prognostic molecular biomarker in PCa [41] . Nevertheless, the role of hsa-miR-6825-5p in PCa development has not been explored. Moreover, Wang et al. implied that N-cadherin promotes cancer stem cell-like traits and epithelial-mesenchymal transition through the ErbB signaling pathway in PCa cells [42] . In this study, MAP2K7 was regulated by hsa-miR-6825-5p, and hsa-miR-6825-5p was regulated by MAFG-AS1 lncRNA in the ceRNA network.
Moreover, the enrichment analysis found that MAP2K7 is involved in the ErbB signaling pathway. Altogether, this study revealed that the MAFG-AS1-hsa-miR-6825-5p-MAP2K7 axis could be related to the progression of CRPC treatment.
Although we explored the potential molecular mechanisms of Abiraterone treatment in CRPC occurrence and development using a bioinformatics approach, there still exist some limitations in our study. For instance, relevant experiments, including cell biology assays, and animal and clinical studies need to be performed to verify the multiple target candidates and signaling pathways identi ed by our bioinformatics analyses.

Conclusions
In summary, NOTCH1, GRIN1, STUB1, MBNL2, and the ErbB signaling pathway were likely related to the progression of CRPC treatment. Moreover, the MAFG-AS1-hsa-miR-6825-5p-MAP2K7 axis might be a therapeutic target for Abiraterone in CRPC treatment. This study may provide new strategies for further studies of Abiraterone in CRPC treatment. Not applicate, this research did not include in vivo and human experiment, so Ethics approval and consent to participate is not included.

Consent for publication:
Written informed consent for publication was obtained from all participants; Availability of data and material: All experiment were nished in our Lab, all data used in this article are promising true and available, anyone who is interested in this research can contact me by E-mail: dupeng9000@126.com, I can sent the original data by E-mail; Competing interests: All authors declared no competing interests; Funding: The research was nancially supported by The Capital Clinical Characteristics Application Research.