Transcriptome analysis of differentially expressed genes and pathways associated with mitoxantrone treatment prostate cancer

Abstract The global physiological function of specifically expressed genes of mitoxantrone (MTX)‐resistant prostate cancer (PCa) is unclear. In this study, gene expression pattern from microarray data was investigated for identifying differentially expressed genes (DEGs) in MTX‐resistant PCa xenografts. Human PCa cell lines DU145 and PC3 were cultured in vitro and xenografted into severe combined immunodeficiency (SCID) mice, treated with MTX intragastrically, three times a week until all mice relapsed. Gene expression profiles of the xenografts from castrated mice were performed with Affymetrix human whole genomic oligonucleotide microarray. The Cytoscape software was used to investigate the relationship between proteins and the signalling transduction network. A total of 355 overlapping genes were differentially expressed in MTX‐resistant DU145R and PC3R xenografts. Of these, 16 genes were selected to be validated by quantitative real‐time PCR (qRT‐PCR) in these xenografts, and further tested in a set of formalin‐fixed, paraffin‐embedded and optimal cutting temperature (OCT) clinical tumour samples. Functional and pathway enrichment analyses revealed that these DEGs were closely related to cellular activity, androgen synthesis, DNA damage and repair, also involved in the ERK/MAPK, PI3K/serine‐threonine protein kinase, also known as protein kinase B, PKB (AKT) and apoptosis signalling pathways. This exploratory analysis provides information about potential candidate genes and may bring new insights into the molecular cascade involvement in MTX‐resistant PCa.

inevitable. Mitoxantrone (MTX), a synthetic anthracenedione, has been routinely used for the treatment PCa for its palliative benefit which enhances clinical remission of the PCa patients. However, despite their initial response and survival benefits, the majority of patients eventually develop resistance to these therapies. 4 The antineoplastic activity of MTX is believed to be related to its ability to bind DNA and inhibit DNA topoisomerase II, an essential enzyme in DNA synthesis and meiotic division which is highly expressed in cancer cells. 5 Damage to DNA is a notable inducer of both transient and permanent alterations in cellular phenotypes. The accumulation of DNA lesions leads to genomic instability through chromosomes breaks, amplification of oncogenes and inactivation of tumour suppression genes, driving to the acquisition of a malignant cancer phenotype. 6 However, cancer cells can overcome DNA damage by induction of a DNA damage secretory program such as proliferation, invasion, metastasis, especially treatment resistance can develop through a variety of signal pathways, including base excision repair, nucleotide excision repair, mismatch repair, direct repair and recombinational repair. 7,8 Comparative genomic hybridization can help to identify relevant genes involved in tumour chemotherapy-resistance and to predict response and cancer prognosis.

| Tumour inoculation and treatment
The animal study was carried out in a specific pathogen-free room and was approved by the Medical Ethics Committee of the Zhengzhou University in accordance with the Guide for the Care and Use of Laboratory Animals (NIH publication no. 80- 23, revised 1996).
Four to six weeks old CB-17 male SCID mice were used in the experiment. Cells (1 × 10 6 cells) were injected subcutaneously into both flanks resulting in two tumours per mouse to test the MTX sensitivity. Once tumours became palpable, the mice were randomly divided into four treatment groups (six mice per group). In the first three groups, MTX was administered three times a week at 0.35 mg/ kg, 1 mg/kg and 3.5 mg/kg respectively. The fourth group was treated with physiological saline (control) at the same time-points. In another set of experiment, animals with palpable tumours were also assigned into four groups: MTX (3.5 mg/kg), castration, MTX (3.5 mg/kg) in combination with castration and control. Surgical castration was performed after tumours have developed. MTX and saline were administered intragastriclly in a 100 µL volume three times a week in all experiments. The diameter of subcutaneously growing tumours was measured with a calliper twice a week until the animals were killed after 6 weeks of treatment. Tumour weight was calculated by the formula: Tumour weight (mg) = (length×width 2 )/2.  Images were autogridded and the chemiluminescent signals were quantified, corrected for background and spot and spatially normalized. Differentially expressed genes (DEGs) were identified through filtering the dataset using P-value <0.01 and a signal-tonoise ratio>2 for use in ANOVA statistical analysis.

| Data preprocessing
The analysis was carried out using R package oligo (version 1. 38

| Pathway enrichment and network construction
Two hundred genes with significant differences were intercepted from differential gene expression data and statistically analyzed by the GeneSpring GX software package. The target gene expression data were analysed by Cytoscape and the signal pathway is derived from tumour-related candidate genes in microarray data. In order to find out the DEGs closely related with signal pathway, additional filtering (minimum 3-fold change) was applied to extract the most significant of these genes which were further analysed using Cytoscape software. Those genes with known gene symbols and their corresponding expression values were uploaded into the software. Networks of these genes were algorithmically generated based on their connectivity. Amersham). Densitometric quantification of band intensities was performed with Kodak one-dimensional image analysis software.

| Gene validation by qRT-PCR
Quantitative real-time PCR was performed with QPK-201 SYBR Green master mix (Toyobo, Osaka, Japan) and the ABI 7300 system from Applied Biosystems. The primers used in the study were obtained from Invitrogen (Beijing, China). Thermocycling parameters included a RT step at 50°C for 20 minutes, followed by a DNA polymerase activation step at 95°C for 2 minutes and 50 PCR cycles (95°C for 20 seconds, 60°C for 30 seconds). All reactions were conducted in triplicate. The fold-change in differential expression for each gene was calculated using the comparative C T method.

| Statistical analysis
Results are presented as the mean ± SE of the mean. To determine whether differences between groups were statistically significant, Wilcoxon rank-sum test of variance was performed, and P < 0.05 was considered to indicate a statistically significant difference. SPSS 12.0 software was used for statistical analyses. and dose-dependent manner. As shown in Figure 1A,B, the highest cytotoxicity of MTX was at 72 hours, and IC50 is 0.1 mg/mL for the

| DEGs in MTX-resistant xenografts
To identify the DEGs between MTX-resistant xenografts and their controls, threshold |logFC| >1 and P-value <0.05 were used in comparative analysis. A total of 1849, 2123 genes were extracted from the DU145R and PC3R respectively ( Figure 3A,B). Among them, 986, 851 down-regulated genes, and 863, 1272 up-regulated genes were screened in the DU145R and PC3R xenografts respectively.
We selected top 20 up-regulated and down-regulated genes according to the log ratio expression values (Tables 1 & 2), and among these, the ones in the MTX-resistant groups whose expression levels were changed by more than 3-fold compared with their control groups (P < 0.01). Upon comparison of the DEGs between both MTX-resistant xenografts, 355 genes are overlapped, comprising 131 co-downregulated genes and 224 co-upregulated genes ( Figure 3C, D). After that, the overlapping DEGs were clustered which can well differentiate the MTX-treatment samples from the controls. The heatmap of the overlapping DEGs is shown in Figure 3E.

| Functional and pathway enrichment analysis of overlapping DEGs
Gene Ontology enrichment analysis revealed 355 overlapping genes that are involved in a number of BP including response to hypoxia, transforming growth factor β receptor signalling pathway, signal transduction and chemotaxis ( Figure 4A). In terms of CC, DEGS were mostly enriched in the extracellular exosome, cell surface and lateral plasma membrane ( Figure 4B). Molecular functions analysis indicated that the overlapping DEGs were mainly associated with protein binding, heparin binding and transcription factor binding ( Figure 4C). Subsequential KEGG pathway enrichment analysis revealed that the common down-regulated DEGs primarily enriched in the hippo signalling pathway, pathways in cancer, proteoglycans in cancer and cell adhesion molecules (CAMs) ( Figure 4D). The top five enriched terms were presented in Table 3. These significantly enriched GO function and KEGG pathways could aid further understanding of the roles of these DEGs, involved in the development of MTX-resistant CRPC.

| Identification of candidate MTX-resistant CRPC markers
To further clarify the core genes of DEGs identified in the microarray analysis, protein-protein interaction (PPI) networks were gener-  Table 4.
After performing edge percolated component and shortest path analysis, the most significant modules composed of 10 nodes were screened out from the PPI networks and the hub genes in the networks with a connectivity degree >16 were identified ( Figure 5C,D).
By comparing the hub genes between both MTX-resistant xenografts, PARP1, IL8 and CDH1 are overlapped.

Western blotting and qRT-PCR
The expression patterns of four DEGs, PARP1, IL1B, CDH1 and PLAUR were evaluated by Western blot ( Figure 6A) and quantitative

| DISCUSSION
This study revealed global pathways and networks of DEGs involved in DU145R and PC3R MTX-resistant PCa xenografts. Prostate cancer is a heterogeneous disease and many molecular methods have been used in the search to determine the mechanism behind its development, and find new therapeutic and prognostic targets. 11 Microarray technology for gene expression profiling has proven to be successful in a variety of experimental settings having the potential to discover the diversified and dynamic molecules during tumour progression. 12 Network analysis helps us obtain global and integrated molecular information about interactions among the significant DEGs. One important network was identified around the AKT2 genes. AKT2 is a putative oncogene encoding a protein belonging to the serine-threonine protein kinase, also known as protein kinase B, PKB (AKT) subfamily of serine/threonine-protein kinases, as well as a key node on the phosphatidylinositol 3-kinase (PI3K/AKT) pathway which is recognized as a key pathway in carcinogenesis occurring commonly in diverse human cancer cells. 26,27 Previous studies have revealed that aberrant activation of PI3K/AKT pathway is also closely associated with the process of cancer metastasis. 28  In our study, Ataxia-telangiectasia mutated kinase (ATM) was found to be highly up-regulated in MTX-resistant xenografts, and observed as a hub gene in the network. This gene is also highly expressed in non-small cell lung cancer (NSCL) exposed to carbon ion irradiation. 32 Ataxia-telangiectasia mutated kinase, a member of PI3K/AKT family protein, its main function is to control the cell cycle progression following DNA damage, particularly DSBs. 33 Reports indicated that MTX produce DNA cross-links and DNA replication defects in tumour cells, once DNA damage occurs, ATM pathway for HR repair is activated. 34 DNA repair is an essential prerequisite for

A B
F I G U R E 6 Expression profiles of PARP1, ILB1, CDH1 and PLAUR were evaluated by Western blot and qRT-PCR. A, Western blot. B, Results expressed as western blotting band intensity. C, qRT-PCR. Means ± SEM (n = 4) the maintenance of genomic integrity and cellular viability. 35 It is demonstrated that DNA-damaging drugs (including MTX) can trigger an ATM-dependent DNA damage response, leading to increased cytokine secretion and resistance to chemotherapy-induced apoptosis. 36 Report also revealed that ATM together with phosphotyrosine binding domain and a leucine zipper motif (APPL), a regulator of ATM phosphorylation, modulates DNA damage repair and consequently raises the survival of pancreatic carcinoma cells. 37 In addition, ATM has been implicated in the control of RNA splicing that may play a role in genomic stability. 38,39 DNA-damage response activates a secretory program that comprises a diverse spectrum of proteases, growth factors and, cytokines that have been shown to contribute to wound healing and altered immune responses. 40,41 Interleukin-6 (IL-6), at high concentrations, is found to protect cancer cells from therapy-induced DNA damage, oxidative stress and apoptosis by facilitating the repair and induction of antioxidant response. 42 In our study, IL6 is up-regulated and screened as a hub gene in androgen-independent, MTX-resistant xenografts. Evidence support DNA damage and stress induce sustained IL-6 and IL1B production in PCa through P2Y11 receptor-p38 MAPK-NF-κB signalling pathway, 6,43 which may play an essential role in promoting cancer cells proliferation, survival, invasiveness and metastasis. Therefore, inhibition of IL-6 or in combination with conventional anticancer therapies may be a potential therapeutic strategy for the treatment of MTX-resistant PCa.  F I G U R E 7 qRT-PCR analysis of differentially expressed genes identified in the microarray. A, DU145R vs DU145, B, PC3R vs PC3, C, optimal cutting temperature (OCT) and D, formalin-fixed, paraffin-embedded (FFPE) tumour samples from patients with metastatic castrationresistant prostate cancer. Expression data are represented by a log ratio calculated by comparing ΔCq from the xenograft with ΔCq from the controls. ΔCq was calculated as the difference between Cq of the targeted genes and Cq of the endogenous control gene ACTB Sixteen DEGs were selected to be validated by qRT-PCR according to their function and relative expression level in MTX-resistant xenografts. These DEGs were further tested in MTX-resistant CRPC patient samples (seven OCT and six FFPE samples).
Results confirmed that BLCAP, EML1, FOSL2, ADNP and FRMD3 were deregulated in the same way among the xenografts, FFPE and OCT tumour samples. However, discordant results were observed in the expression of other genes such as PLD1, which was down-expressed in MTX-resistant tumour samples but was significantly overexpressed in MTX-resistant xenografts. These conflicting results may be due to cellular heterogeneity between the xenografts from cell lines and the tumour sample from patients. Therefore, further clinical validation of these results is needed in a large cohort of patients.
In analysing, often individual genes were found in multiple categories of functions related to cancer development including cell signalling, cell death, cellular growth and proliferation. Besides, it reminds us there are certain limits in the analysis as there are many different gene interactions resulting from various cellular/experimental conditions. Nevertheless, this exploratory analysis may be still useful to bestow a theranostic perspective to the current trend of research in PCa or to develop targeted therapies to overcome MTX chemotherapy resistance.

DATA AVAILABILITY
The data of the study have been deposited into the Research Data Deposit (http://www.researchdata.org.cn), with the Approval Number as RDDB2018000423.