Integrative transcriptome analysis identifies genes and pathways associated with enzalutamide resistance of prostate cancer

Abstract

Background: Enzalutamide, a novel androgen receptor (AR) signaling inhibitor, has been widely used to increase survival in patients with castration-resistant prostate cancer. However, resistance to enzalutamide invariably develops.

Methods: To understand the underlying mechanisms of resistance to enzalutamide, we performed integrative analysis on multiple transcriptome datasets to identify those genes constantly up- or down-regulated in response to enzalutamide treatment.

Results: There were 703 and 581 differentially expressed genes derived from enzalutamide-sensitive and -resistant cell lines, respectively. Functional enrichment analysis on these genes demonstrated that biological processes of cell proliferation and ubiquitin mediated proteolysis pathway are specifically disturbed in sensitive cell lines but not resistant ones. Such divergence explained why enzalutamide ineffective for resistant prostate cancer.

Conclusions: Taken together, the present study revealed a set of critical genes, which can provide etiologic clues as to enzalutamide-resistant prostate cancer and guide novel therapeutic approaches.

Keywords:

• Prostate cancer
• enzalutamide resistance
• transcriptome analysis
• gene expression omnibus
• gene ontology
• KEGG pathway

Introduction

Prostate cancer is a common cancer among men [1], in which malignant cells form in prostate tissue and may spread to other parts of the body. Prostate cancer is highly sensitive to hormones, especially androgen [2,3]. To take advantage of this unique property of prostate cancer, androgen deprivation therapy (ADT) with surgical or medical castration is widely adopted as an effective systemic treatment for hormone-sensitive prostate cancer (HSPC) [4,5]. However, the emergence of castration-resistant prostate cancer (CRPC) is almost inevitable within 2–3 years of initiation of ADT [6].

In recent years, a series of new agents are developed to treat CRPC. Among them, enzalutamide and abiraterone acetate have provided additional survival benefit and have been approved by the US Food and Drug Administration (FDA) [7,8]. Enzalutamide is a new generation of androgen-receptor signaling inhibitor, with significantly higher affinity for AR than prior anti-androgens. Furthermore, enzalutamide prevents nuclear translocation and coactivator recruitment of the ligand-receptor complex [9]. And compared with androgen synthesis inhibitors such as abiraterone, a potential advantage of enzalutamide is that the co-administration of steroids is not required [10]. In spite of multiple advantages, unfortunately, response to enzalutamide is temporary for the majority of treated patients due to acquired resistance [11]. As a result, the overall survival rate for CRPC patients is only modestly increased by enzalutamide [12,13]. Understanding the mechanisms of resistance to enzalutamide is critical to future research into targeted therapies.

Microarray is a powerful tool in basic and translational cancer research [14]. Effective analysis of microarray data is highly important for understanding the biological mechanisms of carcinogenesis and the downstream effects of various anticancer agents [15]. Since microarray produces enormous data from each specimen, it offers not only great opportunity for discovery but also great risk of error and misleading output [16]. In this case, combining the results of different independent studies allows more reliable discovery of gene expression patterns due to larger sample size and increased statistical power [17]. The extraction and integration of public microarray data have been extraordinarily valuable in studying prostate cancer. For instance, one study recently identified a novel gene signature that can discriminate between prostate tumor tissue and benign tissue [18]. And another inter-study validated the consistent dysregulation of polyamine and purine biosynthesis pathways in prostate cancer [19]. Experience can be learned from the above research.

The current study was conducted to investigate the underlying mechanisms of resistance to enzalutamide. At first, we implemented an integrative analysis of several genomic expression datasets of enzalutamide-sensitive and enzalutamide-resistant prostate cancer cell lines. A set of genes were retrieved for being consistently up- or down-regulated across datasets upon enzalutamide treatment. Then, these genes were further interrogated in functional analysis, so as to identify the association between key genes and enzalutamide resistance. The findings of our study provided novel gene signatures and mechanistic insights into the acquired resistance to enzalutamide, which might lead to innovation of therapeutic interventions against prostate cancer.

Materials and methods

Raw data collection and review

We searched Gene Expression Omibus (GEO) for all enzalutamide-related mRNA expression studies that were performed on prostate cancer cell lines. We selected and downloaded four profile datasets with GEO accession numbers GSE44905, GSE78201, GSE44927, and GSE69896. As these four datasets included multiple matched treatment and control groups, the datasets were further curated into six custom datasets (Table 1). Samples in each dataset were categorized based on enzalutamide treatment. The objective was to compare global gene expression profiles in enzalutamide-treated cells to vehicle controls.

Table 1. Gene expression datasets extracted from GEO.

Enzalutamide-sensitive
Dataset 1	Enzalutamide treatment	GSM1093805, GSM1093806, GSM1093807
	Vehicle control	GSM1093796, GSM1093797, GSM1093798
Dataset 2	Enzalutamide treatment	GSM2069522, GSM2069523, GSM2069524, GSM2069525
	Vehicle control	GSM2069526, GSM2069527, GSM2069528, GSM2069529
Dataset 3	Enzalutamide treatment	GSM1094017, GSM1094019, GSM1094021, GSM1094037, GSM1094039, GSM1094041
	Vehicle control	GSM1094016, GSM1094018, GSM1094020, GSM1094036, GSM1094038, GSM1094040
Enzalutamide-resistant
Dataset 4	Enzalutamide treatment	GSM1712132, GSM1712138, GSM1712144
	Vehicle control	GSM1712130, GSM1712136, GSM1712142
Dataset 5	Enzalutamide treatment	GSM1712118, GSM1712119, GSM1712120
	Vehicle control	GSM1712121, GSM1712122, GSM1712123
Dataset 6	Enzalutamide treatment	GSM1094023, GSM1094031, GSM1094033, GSM1094035, GSM1094043, GSM1094051, GSM1094053, GSM1094055
	Vehicle control	GSM1094022, GSM1094030, GSM1094032, GSM1094034, GSM1094042, GSM1094050, GSM1094052, GSM1094054

Microarray data analysis

To remove systemic error occurred during the experimentation, each expression dataset was individually pre-processed and normalized with optimal normalization methods (MAS5.0 for Affymetrix and Quantile for Illumina). Samples of each data set were divided into two groups, the enzalutamide-treated and the vehicle control group. The Affymetrix and Illumina probes were annotated with corresponding gene symbols and Entrez IDs by using the SOURCE (http://source-search.princeton.edu/) and g:Profiler (http://biit.cs.ut.ee/gprofiler/gconvert.cgi) [20] online tools, respectively. The Student’s t-test was performed and fold-change was calculated to identify differentially expressed genes (p < .05 and the same direction of fold-change) upon enzalutamide treatment. The global differentially expressed genes across datasets were chosen following the vote counting procedure [21], only those genes consistently up- or down-regulated in the majority of expression datasets were selected.

Gene ontology, pathway, and protein-protein interaction network analysis

In order to interpret the biological difference between ‘sensitive’ differentially expressed genes (sDEGs) and ‘resistant’ differentially expressed genes (rDEGs), we further performed enrichment analyses on biological processes of Gene Ontology (GO) and pathways of Kyoto Encyclopedia of Genes and Genomes (KEGG). The gene sets of interest were investigated by using the WebGestalt online tool (http://www.webgestalt.org/). This tool runs hypergeometric test and evaluates which functional categories are significantly enriched in a gene list as compared with a reference gene set (usually the genome background). Pathways and GO terms with p values <.05 were selected as significantly enriched category.

Results

Gene expression patterns in enzalutamide-resistant and -sensitive cells

Raw data were retrieved from GEO datasets GSE44905, GSE78201, GSE44927, and GSE69896, which comprised enzalutamide-treated and vehicle control samples. From these datasets, we further curated three custom datasets of enzalutamide-resistant cells and three datasets of enzalutamide-sensitive cells (Table 1). For each dataset, the enzalutamide-treated samples were compared to the vehicle controls by performing Student’s t-test test and calculating expression fold-change.

With regard to enzalutamide-sensitive cells, a total of 703 sDEGs were found significant (p < .05) in at least two of three datasets and having the same direction of the fold change (Supplement S1). On the other hand, 581 rDEGs were identified from enzalutamide-resistant cells (Supplement S2). Of note, the sDEGs and rDEGs were highly distinct from each other with only 33 genes in common (Figure 1(A)). In addition, we noticed that most sDEGs (418 out of 703) were down-regulated upon enzalutamide treatment while most rDEGs (446 out of 581) were up-regulated (Figure 1(B)). Such divergence (Fisher’s exact test p = 5.26 × 10⁻⁴⁰) suggested that sensitive and resistant cells could activate highly different mechanisms in response to enzalutamide treatment.

Figure 1. Divergence between sDEGs and rDEGs. (A) sDEGs and rDEGs had only 33 genes in common. (B) Nearly 60% of sDEGs were down-regulated while over 76% of rDEGs were up-regulated.

Gene ontology annotation

The differentially expressed genes served as a starting point to explain the resistance to enzalutamide. To gain insights into the biological difference between sDEGs and rDEGs, we used the WebGestalt online server [22] to perform enrichment analysis with regard to GO [23] biological process. The rDEGs and sDEGs were tested against the background set of all genes in human genome. Hypergeometric test with p < .05 was adopted as the criteria for statistical significance. We found that the sDEGs were significantly enriched in GO terms related to cell proliferation, including chromosome segregation (GO:0007059), DNA replication (GO:0006260), DNA repair (GO:0006281), cell cycle phase transition (GO:0044770), cell cycle checkpoint (GO:0000075), mitotic nuclear division (GO:0007067), telomere organization (GO:0032200), anatomical structure homeostasis (GO:0060249), and DNA strand elongation (GO:0022616) (Figure 2(A)) (Supplement S3). On the other hand, the rDEGs were not significantly enriched in this category (Figure 2(B)) (Supplement S4). Such difference indicated that certain biological processes of cell proliferation are specifically disturbed in sensitive cell lines but not resistant ones.

Figure 2. The top 10 enriched GO terms and the numbers of differentially expressed genes. (A) Biological processes for sDEGs. (B) Biological processes for rDEGs.

KEGG pathway enrichment

To further evaluate the biological significance for the differentially expressed genes, we also performed pathway enrichment analysis based on the KEGG database [24]. We found that the sDEGs were significantly enriched in pathways such as ubiquitin mediated proteolysis (hsa04120), steroid biosynthesis (hsa00100), antifolate resistance (hsa01523), lysine degradation (hsa00310), DNA replication (hsa03030), pathways in cancer (hsa05200), fanconi anemia pathway (hsa03460), one carbon pool by folate (hsa00670), mismatch repair (hsa03430), histidine metabolism (hsa00340) (Figure 3(A)) (Supplement S5). On the other hand, the sDEGs were significantly enriched in small cell lung cancer (hsa05222), toxoplasmosis (hsa05145), Fc gamma R-mediated phagocytosis (hsa04666), circadian entrainment (hsa04713), other types of O-glycan biosynthesis (hsa00514), sphingolipid signaling pathway (hsa04071), non-small cell lung cancer (hsa05223), cholinergic synapse (hsa04725), axon guidance (hsa04360), AMPK signaling pathway (hsa04152), and other pathways (Figure 3(B)) (Supplement S6).

Figure 3. The top 10 enriched KEGG pathways and the numbers of differentially expressed genes. (A) Pathways for sDEGs. (B) Pathways for rDEGs.

A number of the above pathways can be further explored. For instance, ubiquitin mediated proteolysis (hsa04120) was the top pathway enriched with sDEGs, but was not significantly enriched with rDEGs. A series of genes related to this pathway were generally up-regulated (e.g. PML, SOCS3, MID1, UBE2H, FBXO2, RNF7, VHL, HERC4, CUL5, UBE3B, and UBE2K) or down-regulated (e.g. UBE3C, PIAS2, UBE2D2, and UBE2I) upon enzalutamide treatment only in sensitive cell lines but not resistant ones (Figure 4). Such differences suggested that enzalutamide may exert antineoplastic effect in sensitive cell lines through the ubiquitin-proteasome system, which is abnormal in resistant cell lines and leads to enzalutamide resistance.

Figure 4. Genes involved in the ubiquitin mediated proteolysis pathway were disturbed by enzalutamide in sensitive rather than resistant cell lines.

Discussion

Prostate cancer is recognized as the second leading cause of cancer-related deaths in men [25,26]. Despite advances in diagnosis and treatment, morbidity from prostate cancer remains high. Enzalutamide has been placed great expectations to be effective in increasing survival rate. However, the application of enzalutamide therapy is complicated by the development of resistance. The molecular response to enzalutamide treatment in patients with prostate cancer is largely undefined. Consequently, there is an urgent need for transcriptomic biomarkers to guide personalized therapy selection and elucidate molecular mechanisms of enzalutamide resistance.

In this study, we combined differentially expressed genes from multiple transcriptomic datasets to highlight genes that were consistently up- or down-regulated. In total, 703 sDEGs and 581 rDEGs were identified from enzalutamide-sensitive and enzalutamide-resistant cells, respectively. While most sDEGs were down-regulated, the great majority of rDEGs were up-regulated. Moreover, we found that sDEGs and rDEGs were highly distinct from each other with little overlap. Such divergence indicated that the response to enzalutamide treatment may be different in sensitive and resistant prostate cancer cells.

In order to uncover the biological differences between enzalutamide-sensitive and enzalutamide-resistant cells, GO enrichment analysis was performed on the differentially expressed genes. We found that sDEGs, but not rDEGs, were highly enriched in a series of GO terms related to cell proliferation, such as chromosome segregation. It has been proved by many studies that chromosome segregation errors can promote chromosomal instability and foster multi-drug resistance in vitro and in vivo [27,28]. Our results suggested that chromosome mis-segregation may also induce enzalutamide resistance in prostate cancer.

Furthermore, the differentially expressed genes were analyzed for KEGG pathway enrichment. As shown in Figure 4, a number of genes involved in the ubiquitin mediated proteolysis pathway were disturbed by enzalutamide only in sensitive cell lines but not resistant ones. For example, SOCS3 (suppressor of cytokine signaling 3) is a tumor suppressor gene [29], whose loss of expression is an early event in carcinogenesis [30]. It has been found that inactivation of SOCS3 in prostate cancer can reduced the efficacy of enzalutamide [31]. In line with previous finding, our data showed that enzalutamide can restore the expression of SOCS3 in only sensitive cells but not resistant ones. As another example, PML (promyelocytic leukemia) protein expression has been found to be reduced or abolished in prostate cancer and cancers of other histologic origins [32]. Our results suggested that PML was up-regulated in most enzalutamide-sensitive cell lines, which could suppress prostate cancer cell growth by enhancing p53 activity and inhibiting AR transactivation [33]. Therefore, the absence of elevated PML expression in resistant cells may partially explain the inefficacy of enzalutamide. The above difference in ubiquitin-proteasome system between sensitive and resistant cell lines provided new clues as to the molecular mechanisms of enzalutamide resistance in prostate cancer.

The present study also has some limitations that should be taken into account. First, the results could be distorted by heterogeneity between datasets and confounding factors in experiments. Although we conducted global normalization for different datasets, the heterogeneity of cell culture protocols and various microarray platforms used in different studies can hardly be removed. Second, the current data were derived from in vitro system, which may not fully represent in vivo pathology and pharmacology. Third, the technical limitations of microarray and RNA-seq platforms, such as low accuracy regarding transcripts with low expression, cannot be overcome in this data analysis. However, our integrative analysis on multiple datasets enabled us to detect differentially expressed genes that may otherwise be missed in single-dataset analysis. Despite these limitations, our discovery has important implications in the molecular mechanisms of resistance to enzalutamide. Further experimental research is still required to confirm and extend our results.

In summary, our integrative transcriptome analysis maximized the value of microarray expression data, so as to correlate differentially expressed genes with enzalutamide resistance in prostate cancer. The results pointed to a set of biological processes and signaling pathways that may play substantial role in this issue. These insights suggested several new avenues for mechanistic research into cancer drug resistance and potential novel therapies.

Disclosure statement

No conflict of interest was reported by the authors.

Additional information

Funding

This study is supported by the Science and Technology Commission of Shanghai Municipality [No. 15ZR1427600].