Analysis of the genomic response of human prostate cancer cells to histone deacetylase inhibitors

Abstract

Histone deacetylases (HDACs) have emerged as important targets for cancer treatment. HDAC-inhibitors (HDACis) are well tolerated in patients and have been approved for the treatment of patients with cutaneous T-cell lymphoma (CTCL). To improve the clinical benefit of HDACis in solid tumors, combination strategies with HDACis could be employed. In this study, we applied Analysis of Functional Annotation (AFA) to provide a comprehensive list of genes and pathways affected upon HDACi-treatment in prostate cancer cells. This approach provides an unbiased and objective approach to high throughput data mining. By performing AFA on gene expression data from prostate cancer cell lines DU-145 (an HDACi-sensitive cell line) and PC3 (a relatively HDACi-resistant cell line) treated with HDACis valproic acid or vorinostat, we identified biological processes that are affected by HDACis and are therefore potential treatment targets for combination therapy. Our analysis revealed that HDAC-inhibition resulted among others in upregulation of major histocompatibility complex (MHC) genes and deregulation of the mitotic spindle checkpoint by downregulation of genes involved in mitosis. These findings were confirmed by AFA on publicly available data sets from HDACi-treated prostate cancer cells. In total, we analyzed 375 microarrays with HDACi treated and non-treated (control) prostate cancer cells. All results from this extensive analysis are provided as an online research source (available at the journal’s website and at http://luigimarchionni.org/HDACIs.html). By publishing this data, we aim to enhance our understanding of the cellular changes after HDAC-inhibition, and to identify novel potential combination strategies with HDACis for the treatment of prostate cancer patients.

Keywords: :

• analysis of functional annotation
• HDACis
• prostate cancer
• mitotic spindle checkpoint
• major histocompatibility complex
• valproic acid
• vorinostat
• gene expression analysis

Introduction

An important mechanism of cells to epigenetically regulate gene expression is by acetylating and deacetylating histones.¹ Histone deacetylases (HDACs) are a class of enzymes that deacetylate lysine residues in the N-terminal tails of histones, thereby blocking gene transcription.¹ HDACs are frequently overexpressed in cancer; their overexpression leads among others to epigenetic silencing of tumor suppressor genes.¹ Therefore, various HDAC-inhibitors (HDACis) have been developed for cancer therapy, of which vorinostat (SAHA) and Romidepsin are approved by the United States Food and Drug Administration (US FDA) for the treatment of cutaneous T-cell lymphomas (CTCL). HDACis arrest cells in G0/G1 or G2/M phase dependent on the dose of HDACi and/or cell type used.² Despite pre-clinical data showing great promise and their success in liquid tumors, the potential of HDACis as single agents against solid tumors, specifically prostate cancer (PCa), seems to be limited in clinical studies.²

It seems that improving DNA accessibility with HDACis is merely the first step in cancer treatment. Recent studies have therefore focused on combination strategies involving HDACis, with success. Valproic acid (VPA) in combination with epirubicin/FEC (5-fluorouracil, epirubicin, cyclophosphamide) resulted in an objective response in 64% of patients with solid advanced malignancies.³ Combination therapy with the HDACi magnesium valproate and DNA demethylating agent hydralazine resensitized 80% of cancer patients to chemotherapy on which they had previously progressed.⁴ This combination was successfully added to doxorubicin and cyclophosphamide therapy in breast cancer patients as well.⁵ The addition of vorinostat to the mammalian target of rapamycin (mTOR) inhibitor temsirolimus improved anti-cancer activity against renal cell carcinoma in vitro and in vivo.⁶ Other recent preclinical studies indicated that HDACis such as VPA may sensitize cancer cells, among others PCa cells, to radiotherapy.^7,8 In non-small cell lung cancer studies it was found that cells may be sensitized for radiotherapy through acetyl p53-mediated downregulation of c-myc.⁹ The rationale for such combination studies with HDACis was that HDACis may reverse epigenetic changes made by the tumor, downregulate gene expression involved in DNA damage repair and/or upregulate apoptosis in cancer cells.

In this study, we apply analysis of functional annotation (AFA) to HDACi-treated PCa cells, thereby providing a rationale for novel combination strategies with HDACis. AFA is a high-throughput bioinformatics approach to identify sets of genes that are differentially expressed between conditions, such as cancer cells pre- and post-treatment. It is conceptually similar to gene set enrichment analysis (GSEA).^10-14 This unbiased method enables the interpretation of large amounts of gene expression data generated by microarray analysis through superimposition, selection, analysis and visualization of information encompassing distinct biological concepts, such as cellular signaling pathways, protein-protein interaction (PPI) networks, gene ontology (GO), gene expression regulation by transcription factors and microRNA targets. In our study AFA was used to detect cellular processes that are affected by HDACis in PCa cell lines. We analyzed data from an elaborate microarray experiment with an HDACi-sensitive (DU-145) and a (relatively) HDACi-insensitive PCa cell line (PC3) treated with the HDACis VPA and vorinostat (SAHA).¹⁵ We further complemented this analysis with gene expression profiles from HDACi-treated LNCaP and PC3 cells using all available data from the public domain.^16,17 Overall we analyzed 375 distinct microarray experiments involving HDACi-treated PCa cells, which to our best knowledge is the largest analysis of this kind to date. By creating an encyclopedia of expression changes induced by HDACis in PCa cell lines, we have implemented a resource publicly available to the research community for use in future research, both for improving our understanding of the cellular effects of HDACis on PCa cells, and for developing hypotheses to test drug combinations with HDACis in preclinical studies. Finally, based on the AFA results, we present a rationale for the combination of HDACis with mitotic spindle checkpoint inhibitors and immunotherapy. Further (pre-)clinical evaluation is necessary to validate these rationales.

Results

Analysis of functional annotation (AFA) after treatment with HDACis

We previously used a “multiple-loop, double-cube” design to identify genes differentially expressed in PCa cell lines upon HDAC-inhibition by VPA or SAHA.¹⁵ In the current study, we applied AFA, as previously described, to enable the interpretation of these results in the context of relevant cancer biology.^10,11,18-20 To this end, we selected Functional Gene Sets (FGS) from distinct databases, recapitulating different and complementary biological concepts: (1) cellular signaling pathways from Pathway Commons, (2) PPI networks from the National Center for Biotechnology Information (NCBI) Entrez Gene database, (3) downstream transcriptional responses, and (4) gene expression regulatory networks orchestrated by transcription factors and microRNA targets. These collections included, among others, the Human Protein Reference Database (HPRD), GO, Kyoto Encyclopedia of Genes and Genomes (KEGG), the Molecular Signature Database (MSigDB), the Pathway Commons and NCBI Entrez Gene databases.^14,21-30 FGS whose expression altered most significantly upon HDAC-inhibition in DU-145 and PC3 cells across all conditions after correction for multiple testing (adjusted p-value < 0.05, top five FGS or more in case of ties), were retrieved and are displayed in Figure 1. This approach enabled the identification of biological themes that are differentially expressed upon HDAC-inhibition across all conditions irrespective of the gene expression direction change (Fig. 1A), as well as biological processes that are selectively up- or downregulated across all conditions (Fig. 1B and C, respectively). The complete lists with all AFA results, including FGS that were differentially expressed between cell lines, treatments and/or time-points, with all relevant statistics (e.g., fold-change, p-values, false discovery rates (FDR), log-odds), are provided as supplementary data available online at the journal’s website and at http://luigimarchionni.org/HDACIs.html.

Figure 1. Heat maps visualizing differentially expressed functional gene sets (FGS) as determined by analysis of functional annotation (AFA) performed on DU-145 and PC3 cells treated with HDAC-inhibitors (HDACis) VPA or vorinostat for 48 or 96 h. Each row represents a distinct FGS, while each column represents a distinct coefficient from our previous linear model analysis. The FGS that were most significantly differentially expressed across all comparisons performed are shown in the figure (top 5 FGS showing an adjusted p-value ≤ 5%, or more in case of ties). Color scales correspond to the absolute adjusted p-values obtained from our analysis after base 10 logarithmic transformations (i.e., the number on the color scale increases with decreasing adjusted p-values). Differentially expressed FGS were selected from different collections in order to encompass distinct biological concepts, as shown by the color bar on the left of each heat map. Cell signaling FGS are highlighted in red and yellow (Pathway Commons Reactome and NCI pathways, respectively), signaling pathway target gene sets in green (Human Protein Reference Database, HPRD), protein-protein-interaction networks in cyan (PPI, as compiled in the NCBI Entrez Gene database), FGS for shared transcriptional factor binding sites (TFBS) in blue, and microRNA (MIR) targets gene sets in pink (both from the Broad Institute Molecular Signature Database collections). For complete tables and heat maps visualizing all differentially expressed FGS, see supplementary data (available at the journal’s website and at http://luigimarchionni.org/HDACIs.html). (A) Differentially expressed FGS after ordering the genes based on absolute moderated t-statistics, thus irrespective of the direction of gene expression modulation upon HDAC-inhibition. (B) FGS differentially expressed due to upregulated genes, as determined by signed t-statistics. (C) FGS differentially expressed due to downregulated genes, as determined by signed t-statistics.

We further explored the AFA results by analyzing the relationship among identified FGS in terms of gene overlap, identifying common and distinct sets of genes driving the enrichment. To this end we assembled “gene to FGS” membership matrices for enriched FGS presented in Figure 1, and performed hierarchical clustering and social network analysis to group them based on gene overlap (for details see Figs. S1-S6 and the supplementary methods section).³¹ This revealed both FGS groups whose enrichment was driven by similar, if not identical, genes, as well as FGS groups with minimal gene overlap, whose enrichment was hence due to separate and distinctive genes. Even more interesting, such divergent FGS groups pointed to complementary, yet distinct, biological concepts. For instance, a number of different FGS corresponding to specific signaling pathways and to their downstream targets were enriched, with minimal overlap of the underlying genes, providing a stronger rationale for pathway modulation upon HDACi-treatment. An example of such findings is the epidermal growth factor signaling, for which the FGS corresponding to both the pathway machinery (e.g., “ErbB1 receptor signaling network” from National Cancer Institute (NCI)-Pathway Commons) and the downstream targets (e.g., “EGFR Signaling Pathway” targets from HPRD) proved to be enriched (Fig. 1B; Figs. S2 and S5).

Overall, treatment with HDACis altered gene expression programs related to biological processes that are well known to be involved in tumor development and progression, such as cell adhesion and cell cycle related FGS. In addition, well-established potential drug targets were modulated, such as the c-myc, the androgen receptor signaling, and the Epidermal Growth Factor Receptor 1 (EGFR1) signaling pathways. Furthermore, the interpretation of the AFA results attracted attention to two groups of functionally connected FGS that were up- or downregulated after HDAC-inhibition: major histocompatibility complexes (MHC) genes and mitotic spindle checkpoint genes, respectively (Fig. 1). These FGS will be further discussed below.

Treatment of PCa cells with HDACis results in upregulation of major histocompatibility complex (MHC) genes

Interpretation of the AFA results indicated that HDACis may play a role in immuno-editing, as it transcriptionally modulated genes involved in this biological process. In fact, the “Signaling in Immune system” Reactome FGS proved to be most significantly upregulated across all studied conditions along with FGS related to IL-4 and Interferon signaling pathways (Fig. 1B). Beta-2 microglobulin (B2M), Human Leukocyte Antigen (HLA) class I molecules, and Interferon induced proteins were among the most upregulated genes driving this upregulation (Table 1). Upregulation of FGS related to the immune system was observed in other collections as well: for example PPI networks with LILRB1 (Leukocyte immunoglobulin-like receptor B1) and LILRB2 were also upregulated upon HDAC-inhibition (Fig. 1B). Furthermore, the cAMP Response Element Binding protein (CREB) and the interferon-stimulated response element (ISRE) site gene sets were differentially expressed upon HDACi-treatment (see supplementary data online available at the journal’s website and at http://luigimarchionni.org/HDACIs.html). Interestingly, CREB is part of an enhanceosome that increases transcription of MHC-I, B2M and MHC-II genes, and interferon is a major regulator of MHC-I transcription.^32,33

Table 1. Differentially expressed genes from the “Signaling in Immune system” FGS

			DU-145								PC3
			SAHA				VPA				SAHA				VPA
			48 h		96 h		48 h		96 h		48 h		96 h		48 h		96 h
SYMBOL	GENENAME	ENTREZID	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR
B2M	β-2-microglobulin	567	0.74	0.000021	0.80	0.000007	0.96	0.000002	0.92	0.000002	0.61	0.000138	0.68	0.000046	0.60	0.000172	0.56	0.000217
C3	complement component 3	718	1.19	0.000025	–	–	1.12	0.000025	0.85	0.000165	–	–	–	–	–	–	–	–
CD44	CD44 molecule (Indian blood group)	960	0.76	0.000151	0.74	0.000125	–	–	–	–	–	–	–	–	-0.78	0.000211	–	–
CDH1	cadherin 1, type 1, E-cadherin (epithelial)	999	–	–	–	–	1.07	0.000007	1.24	0.000002	–	–	–	–	–	–	–	–
CFD	complement factor D (adipsin)	1675	–	–	–	–	0.54	0.000005	0.45	0.000013	–	–	–	–	–	–	–	–
CSF2	colony stimulating factor 2 (granulocyte-macrophage)	1437	–	–	–	–	–	–	–	–	–	–	–	–	–	–	0.96	0.000588
DDX58	DEAD (Asp-Glu-Ala-Asp) box polypeptide 58	23586	–	–	–	–	–	–	0.63	0.000172	–	–	0.69	0.000193	–	–	–	–
EGR1	early growth response 1	1958	–	–	–	–	–	–	1.68	0.000011	–	–	–	–	–	–	–	–
FYN	FYN oncogene related to SRC, FGR, YES	2534	–	–	–	–	0.67	0.000022	0.78	0.000004	–	–	–	–	–	–	0.51	0.000334
GAB2	GRB2-associated binding protein 2	9846	–	–	–	–	–	–	0.43	0.000148	–	–	–	–	–	–	0.42	0.000398
GBP2	guanylate binding protein 2, interferon-inducible	2634	–	–	–	–	0.52	0.000019	0.64	0.000003	–	–	–	–	–	–	–	–
HLA-A	major histocompatibility complex, class I, A	3105	0.71	0.000221	0.66	0.000263	1.01	0.000010	1.23	0.000002	0.67	0.000533	–	–	–	–	0.69	0.000365
HLA-B	major histocompatibility complex, class I, B	3106	0.85	0.000105	0.90	0.000046	0.69	0.000043	1.01	0.000002	0.90	0.000115	0.71	0.000069	–	–	0.94	0.000065
HLA-C	major histocompatibility complex, class I, C	3107	0.88	0.000010	0.98	0.000003	1.02	0.000002	1.33	0.000000	0.82	0.000022	0.92	0.000008	0.78	0.000034	1.09	0.000003
HLA-F	major histocompatibility complex, class I, F	3134	0.57	0.000295	0.62	0.000121	0.85	0.000010	1.13	0.000001	0.61	0.000294	0.66	0.000135	–	–	0.73	0.000065
HLA-G	major histocompatibility complex, class I, G	3135	–	–	–	–	0.81	0.000243	1.19	0.000008	0.79	0.000658	0.79	0.000497	–	–	0.87	0.000268
ICAM1	intercellular adhesion molecule 1	3383	0.72	0.000009	0.81	0.000002	0.52	0.000052	0.65	0.000007	–	–	–	–	–	–	–	–
IFI6	interferon, α-inducible protein 6	2537	0.55	0.000101	0.74	0.000007	0.61	0.000029	0.90	0.000001	–	–	–	–	–	–	–	–
IFIT1	interferon-induced protein with tetratricopeptide repeats 1	3434	–	–	–	–	–	–	1.03	0.000027	–	–	–	–	–	–	–	–
IFIT2	interferon-induced protein with tetratricopeptide repeats 2	3433	–	–	–	–	–	–	–	–	0.47	0.000199	0.52	0.000075	–	–	–	–
IFNGR2	interferon gamma receptor 2 (interferon gamma transducer 1)	3460	–	–	0.52	0.000070	–	–	–	–	–	–	–	–	–	–	–	–
IL18	interleukin 18 (interferon-gamma-inducing factor)	3606	1.16	0.000013	1.02	0.000022	0.76	0.000199	0.77	0.000131	0.93	0.000102	0.79	0.000251	–	–	0.71	0.000595
IL1B	interleukin 1, β	3553	–	–	–	–	–	–	–	–	–	–	–	–	1.59	0.000352	2.11	0.000034
IL1R2	interleukin 1 receptor, type II	7850	–	–	–	–	–	–	–	–	–	–	1.65	0.000433	–	–	1.91	0.000165
IL6	interleukin 6 (interferon, β 2)	3569	–	–	–	–	–	–	–	–	–	–	–	–	–	–	1.31	0.000332
IL6ST	interleukin 6 signal transducer (gp130, oncostatin M receptor)	3572	–	–	0.40	0.000286	–	–	–	–	–	–	–	–	–	–	–	–
IP6K2	inositol hexakisphosphate kinase 2	51447	–	–	0.86	0.000041	–	–	–	–	–	–	–	–	–	–	–	–
IRF5	interferon regulatory factor 5	3663	–	–	–	–	0.51	0.000033	0.61	0.000006	–	–	–	–	–	–	–	–
ISG20	interferon stimulated exonuclease gene 20kDa	3669	–	–	–	–	–	–	–	–	0.71	0.000539	–	–	–	–	–	–
ITGB1	integrin, β 1 (fibronectin receptor, β polypeptide, antigen CD29 includes MDF2, MSK12)	3688	0.28	0.000252	0.58	0.000001	0.39	0.000015	0.33	0.000042	–	–	0.44	0.000011	0.46	0.000012	0.40	0.000025
JUN	jun proto-oncogene	3725	–	–	0.80	0.000048	–	–	0.86	0.000019	–	–	0.71	0.000196	–	–	–	–
KLRC1	killer cell lectin-like receptor subfamily C, member 1	3821	0.76	0.000073	0.66	0.000150	–	–	–	–	–	–	–	–	–	–	–	–
LY96	lymphocyte antigen 96	23643	–	–	0.56	0.000019	0.42	0.000160	0.68	0.000003	–	–	0.44	0.000193	–	–	0.43	0.000280
MAP3K7	mitogen-activated protein kinase kinase kinase 7	6885	–	–	–	–	0.31	0.000213	–	–	–	–	–	–	–	–	–	–
MEF2C	myocyte enhancer factor 2C	4208	–	–	–	–	–	–	0.32	0.000182	–	–	0.33	0.000328	–	–	–	–
OAS3	2'-5′-oligoadenylate synthetase 3, 100kDa	4940	–	–	–	–	–	–	–	–	–	–	–	–	–	–	0.44	0.000255
OASL	2'-5′-oligoadenylate synthetase-like	8638	–	–	–	–	0.28	0.000211	0.42	0.000006	–	–	–	–	–	–	–	–
PAG1	phosphoprotein associated with glycosphingolipid microdomains 1	55824	–	–	–	–	0.58	0.000081	0.66	0.000022	–	–	–	–	0.51	0.000486	0.59	0.000137
PELI3	pellino E3 ubiquitin protein ligase family member 3	246330	–	–	–	–	–	–	–	–	–	–	–	–	–	–	0.42	0.000345
PIK3C3	phosphoinositide-3-kinase, class 3	5289	0.63	0.000103	0.68	0.000038	–	–	0.59	0.000083	–	–	–	–	–	–	–	–
PML	promyelocytic leukemia	5371	–	–	–	–	0.70	0.000005	0.75	0.000002	–	–	–	–	–	–	–	–
PTPN1	protein tyrosine phosphatase, non-receptor type 1	5770	0.42	0.000357	0.51	0.000064	–	–	–	–	–	–	–	–	–	–	–	–
PTPN2	protein tyrosine phosphatase, non-receptor type 2	5771	–	–	–	–	–	–	–	–	–	–	–	–	–	–	-0.71	0.000545
PVRL2	poliovirus receptor-related 2 (herpesvirus entry mediator B)	5819	–	–	–	–	–	–	0.57	0.000022	–	–	–	–	–	–	0.56	0.000065
RAP1A	RAP1A, member of RAS oncogene family	5906	0.37	0.000421	0.40	0.000187	0.41	0.000139	0.42	0.000096	–	–	–	–	–	–	–	–
RIPK2	receptor-interacting serine-threonine kinase 2	8767	–	–	–	–	0.54	0.000051	0.60	0.000016	–	–	–	–	–	–	–	–
SHC1	SHC (Src homology 2 domain containing) transforming protein 1	6464	0.58	0.000170	0.60	0.000092	–	–	–	–	–	–	–	–	–	–	–	–
SOCS3	suppressor of cytokine signaling 3	9021	–	–	0.37	0.000153	–	–	–	–	–	–	–	–	–	–	–	–
SP100	SP100 nuclear antigen	6672	0.46	0.000154	0.63	0.000010	–	–	–	–	–	–	–	–	–	–	–	–
STAT1	signal transducer and activator of transcription 1, 91kDa	6772	–	–	–	–	0.32	0.000086	–	–	–	–	–	–	–	–	–	–
TICAM1	toll-like receptor adaptor molecule 1	148022	0.46	0.000026	0.71	0.000001	0.37	0.000080	0.58	0.000002	0.48	0.000032	0.72	0.000001	0.57	0.000010	0.78	0.000001
TRAF2	TNF receptor-associated factor 2	7186	-0.39	0.000232	–	–	-0.37	0.000242	–	–	–	–	–	–	–	–	–	–
YWHAB	tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, β polypeptide	7529	–	–	0.41	0.000118	–	–	–	–	–	–	0.39	0.000320	–	–	–	–
YWHAZ	tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide	7534	–	–	–	–	–	–	–	–	0.51	0.000133	0.50	0.000112	0.44	0.000389	–	–

Table displaying genes from the “Signaling in Immune system” FGS (Pathway Commons Reactome database) that were significantly differentially expressed (Log Odds > 3, FDR < 0.007) in DU-145 and PC3 cell lines upon HDACi-treatment with VPA or vorinostat (SAHA). In both cell lines the FGS alteration is driven by gene upregulation. P-values per gene can be found in the supplementary data (available in the Supplemental Material and at http://luigimarchionni.org/HDACIs.html).

We validated our data by performing differential gene expression analysis and AFA on an additional 353 microarray experiments from three publicly available data sets (GSE8645 and GSE31620, and Connectivity Map) involving PCa cells treated with HDACis.^16,17 We assessed the overall agreement of differential gene expression across all analyzed data sets using Correspondence At the Top (CAT) curves as described previously, which overall proved to be greater than expected by chance (Figs. S7-S16).^10,18,34 We also compared and validated AFA results identifying the biological themes and processes which were consistently differentially expressed across the data sets, treatment conditions, and data points analyzed (Figs. S17-S19). Overall, FGS related to the immune system were significantly upregulated across most of the analyzed conditions, spanning distinct HDACis, cell lines and treatment time points (Fig S18; the complete results in the Supplemental Material and at http://luigimarchionni.org/HDACIs.html).

Based on the finding that FGS related to the immune system are differentially expressed upon HDACi-treatment, we hypothesized that PCa cells may alter the expression of genes related to the immune system, a change that is reverted by HDACis. To further explore this hypothesis, we assessed whether the expression of MHC proteins in human PCa tissue is deregulated compared with normal prostate tissue. For this purpose, a PCa Tissue Microarray (TMA) was performed using specimens from 71 patients who had undergone a radical prostatectomy. Protein expression of HLA-A, HLA-B, HLA-DRA and B2M was scored in PCa tissue and normal surrounding tissue (Fig. 2). Repeated measurement analysis using the panel of all four proteins revealed a significant combined change in MHC protein expression in human PCa tissue compared with normal surrounding tissue (p < 0.001). B2M and HLA-B were primarily responsible for this deregulation. B2M protein expression was downregulated compared with normal surrounding tissue (p = 2.5·10⁻⁷), while HLA-B was upregulated (p = 0.009), as determined by Student’s t-tests (Fig. 2B). There was no eminent change in HLA-A and HLA-DRA protein expression. As a final note, B2M expression of non-cancerous tissue tended to be decreased in patients with PSA recurrence (p = 0.06), possibly indicating that deregulation of B2M is correlated with PCa progression (data not shown).

Figure 2. IHC staining of human prostate tissue. A TMA performed with PCa tissue and normal surrounding tissue from 71 patients with primary PCa indicated deregulated MHC protein expression in human PCa tissue (p < 0.001). (A) Representative images of human PCa and normal prostate tissue stained for B2M (top) and HLA-B (bottom) via immunohistochemistry show decreased expression of B2M and increased expression of HLA-B in PCa tissue compared with normal surrounding tissue in a PCa patient. (B) B2M, HLA-A, HLA-B and HLA-DRA protein expression was scored by measuring the staining intensity of these proteins in immunohistochemistry samples from 71 PCa patients by two investigators independently (0, no staining; 1, light staining; 2, medium staining; 3, strong staining). The mean staining intensity of 71 cancerous and non-cancerous tissue samples is displayed in the bar graph. B2M and HLA-B protein expression was significantly different between PCa and normal surrounding tissue (p = 2.5·10⁻⁷ and p = 0.009, respectively), while there was no eminent change in HLA-A and -DRA protein expression (p > 0.05).

Treatment of PCa cells with HDACis results in downregulation of mitotic spindle checkpoint genes

Our AFA results indicated that HDACis changed the expression of genes relevant to mitosis as well (Fig. 1; Figs. S17 and S19). Table 2 summarizes the genes contributing to altered expression of the “cell cycle, mitotic” Reactome FGS in DU-145 and PC3 cells treated with VPA or vorinostat. Treatment of DU-145 cells with 1 mM VPA or 1 μM vorinostat resulted in a transcriptional downregulation of genes coding for kinetochore and other mitotic proteins (adjusted p-value < 10⁻⁹). In contrast, the more resistant PC3 cells showed less downregulation (adjusted p-value < 10⁻⁹ only after 96h treatment with VPA) (Fig. 1 and Table 2).

Table 2. Differentially expressed genes from the “Cell Cycle, Mitotic” FGS

			DU-145								PC3
			SAHA				VPA				SAHA				VPA
			48 h		96 h		48 h		96 h		48 h		96 h		48 h		96 h
SYMBOL	GENENAME	ENTREZID	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR	Log2 Fold Change	FDR
ACTR1A	ARP1 actin-related protein 1 homolog A, centractin α (yeast)	10121	–	–	–	–	–	–	–	–	0.36	0.000540	0.39	0.000235	–	–	–	–
ANAPC1	anaphase promoting complex subunit 1	64682	-0.59	0.000293	-0.73	0.000042	-0.57	0.000263	-0.58	0.000180	–	–	–	–	–	–	–	–
AURKA	aurora kinase A	6790	-0.74	0.000005	-0.87	0.000001	-0.60	0.000013	-0.70	0.000003	–	–	–	–	–	–	–	–
AURKB	aurora kinase B	9212	–	–	–	–	-0.37	0.000118	–	–	–	–	–	–	–	–	–	–
BIRC5	baculoviral IAP repeat containing 5	332	-0.67	0.000294	-0.84	0.000039	–	–	–	–	–	–	–	–	–	–	–	–
BUB1	budding uninhibited by benzimidazoles 1 homolog (yeast)	699	-0.70	0.000009	-0.65	0.000009	-0.47	0.000108	-0.64	0.000007	–	–	–	–	–	–	–	–
BUB1B	budding uninhibited by benzimidazoles 1 homolog β (yeast)	701	-0.41	0.000125	-0.41	0.000088	-0.41	0.000084	-0.52	0.000011	–	–	–	–	–	–	–	–
CCDC99	coiled-coil domain containing 99	54908	–	–	-0.30	0.000246	–	–	–	–	–	–	–	–	–	–	–	–
CCNA1	cyclin A1	8900	–	–	–	–	0.56	0.000001	0.54	0.000001	–	–	–	–	–	–	–	–
CCNA2	cyclin A2	890	-0.58	0.000020	-0.56	0.000016	-0.62	0.000002	-0.70	0.000001	–	–	–	–	–	–	-0.36	0.000277
CCNB1	cyclin B1	891	-0.85	0.000295	–	–	–	–	-0.82	0.000227	–	–	–	–	–	–	–	–
CCNB2	cyclin B2	9133	-0.75	0.000028	-0.88	0.000006	–	–	-0.54	0.000165	–	–	–	–	–	–	–	–
CCND1	cyclin D1	595	0.88	0.000053	0.93	0.000024	–	–	0.74	0.000095	–	–	–	–	–	–	–	–
CCNE1	cyclin E1	898	–	–	1.24	0.000091	–	–	–	–	–	–	–	–	–	–	–	–
CDC20	cell division cycle 20 homolog (S. cerevisiae)	991	-0.82	0.000002	-0.87	0.000001	-0.63	0.000007	-0.54	0.000018	–	–	–	–	–	–	–	–
CDC25B	cell division cycle 25 homolog B (S. pombe)	994	–	–	-0.44	0.000121	–	–	–	–	–	–	–	–	–	–	–	–
CDC45	cell division cycle 45 homolog (S. cerevisiae)	8318	-0.84	0.000010	-0.90	0.000003	-0.79	0.000008	-0.77	0.000007	–	–	–	–	–	–	–	–
CDCA8	cell division cycle associated 8	55143	-0.78	0.000027	-0.86	0.000008	-0.65	0.000071	-0.66	0.000047	–	–	–	–	–	–	–	–
CDK5RAP2	CDK5 regulatory subunit associated protein 2	55755	-0.35	0.000159	-0.41	0.000033	–	–	–	–	–	–	–	–	–	–	–	–
CDKN1A	cyclin-dependent kinase inhibitor 1A (p21, Cip1)	1026	–	–	1.80	0.000033	–	–	1.36	0.000220	–	–	–	–	–	–	–	–
CDKN2D	cyclin-dependent kinase inhibitor 2D (p19, inhibits CDK4)	1032	–	–	–	–	0.81	0.000068	1.00	0.000010	–	–	0.69	0.000429	0.99	0.000039	1.18	0.000009
CENPA	centromere protein A	1058	-0.62	0.000132	-0.70	0.000036	–	–	–	–	–	–	–	–	–	–	–	–
CENPF	centromere protein F, 350/400kDa (mitosin)	1063	–	–	-1.10	0.000029	–	–	–	–	–	–	–	–	–	–	–	–
CENPH	centromere protein H	64946	-0.66	0.000301	-0.87	0.000029	–	–	–	–	–	–	–	–	–	–	–	–
CENPM	centromere protein M	79019	-0.86	0.000002	-0.75	0.000003	-0.71	0.000005	-0.64	0.000006	-0.48	0.000223	–	–	–	–	–	–
CENPN	centromere protein N	55839	-0.96	0.000013	-1.06	0.000004	-0.67	0.000129	-0.84	0.000016	–	–	–	–	–	–	–	–
CEP41	centrosomal protein 41kDa	95681	–	–	–	–	0.38	0.000027	0.34	0.000041	–	–	–	–	–	–	–	–
CEP72	centrosomal protein 72kDa	55722	–	–	-0.39	0.000138	–	–	–	–	–	–	–	–	–	–	–	–
CETN2	centrin, EF-hand protein, 2	1069	–	–	0.60	0.000195	–	–	–	–	–	–	0.62	0.000244	–	–	–	–
CKS1B	CDC28 protein kinase regulatory subunit 1B	1163	-0.97	0.000084	-1.17	0.000014	–	–	–	–	–	–	–	–	–	–	–	–
DBF4	DBF4 homolog (S. cerevisiae)	10926	-0.49	0.000063	-0.66	0.000007	–	–	-0.52	0.000018	–	–	–	–	–	–	–	–
DHFR	dihydrofolate reductase	1719	-0.30	0.000275	–	–	–	–	–	–	–	–	–	–	–	–	–	–
DYNC1H1	dynein, cytoplasmic 1, heavy chain 1	1778	0.62	0.000004	0.55	0.000006	–	–	–	–	–	–	–	–	–	–	–	–
E2F1	E2F transcription factor 1	1869	-0.52	0.000238	-0.51	0.000201	-0.54	0.000128	–	–	–	–	–	–	–	–	–	–
FGFR1OP	FGFR1 oncogene partner	11116	-0.64	0.000048	-0.60	0.000049	-0.51	0.000170	-0.50	0.000151	–	–	–	–	–	–	–	–
HAUS2	HAUS augmin-like complex, subunit 2	55142	–	–	–	–	–	–	-0.62	0.000030	–	–	–	–	–	–	–	–
KIF20A	kinesin family member 20A	10112	-0.64	0.000006	-0.62	0.000005	-0.58	0.000008	-0.67	0.000002	–	–	–	–	–	–	–	–
KIF23	kinesin family member 23	9493	–	–	-0.91	0.000046	–	–	–	–	–	–	–	–	–	–	–	–
LIG1	ligase I, DNA, ATP-dependent	3978	-0.60	0.000004	-0.59	0.000003	-0.40	0.000041	-0.34	0.000125	–	–	–	–	–	–	–	–
MAD2L1	MAD2 mitotic arrest deficient-like 1 (yeast)	4085	–	–	-0.82	0.000086	–	–	–	–	–	–	–	–	–	–	–	–
MCM5	minichromosome maintenance complex component 5	4174	–	–	–	–	–	–	–	–	–	–	–	–	-0.51	0.000412	–	–
MCM7	minichromosome maintenance complex component 7	4176	-0.55	0.000239	-0.60	0.000087	-0.70	0.000024	-0.61	0.000054	–	–	–	–	-0.51	0.000614	–	–
MNAT1	menage a trois homolog 1, cyclin H assembly factor (Xenopus laevis)	4331	0.34	0.000106	0.28	0.000304	0.43	0.000011	–	–	0.30	0.000370	–	–	0.46	0.000017	–	–
NDC80	NDC80 kinetochore complex component homolog (S. cerevisiae)	10403	-0.81	0.000021	-0.80	0.000015	–	–	-0.59	0.000112	–	–	–	–	–	–	–	–
NDE1	nudE nuclear distribution E homolog 1 (A. nidulans)	54820	–	–	-0.42	0.000050	–	–	–	–	–	–	–	–	–	–	–	–
NEK2	NIMA (never in mitosis gene a)-related kinase 2	4751	-0.36	0.000132	-0.33	0.000188	–	–	–	–	–	–	–	–	–	–	–	–
NSL1	NSL1, MIND kinetochore complex component, homolog (S. cerevisiae)	25936	–	–	0.47	0.000175	–	–	–	–	–	–	–	–	–	–	–	–
NUF2	NUF2, NDC80 kinetochore complex component, homolog (S. cerevisiae)	83540	-0.78	0.000002	-0.80	0.000001	–	–	-0.50	0.000027	–	–	–	–	–	–	–	–
NUP133	nucleoporin 133kDa	55746	–	–	–	–	–	–	–	–	–	–	–	–	-0.59	0.000444	–	–
NUP160	nucleoporin 160kDa	23279	–	–	–	–	–	–	–	–	–	–	–	–	–	–	-0.35	0.000069
NUP37	nucleoporin 37kDa	79023	–	–	-0.60	0.000042	–	–	–	–	–	–	–	–	–	–	–	–
NUP85	nucleoporin 85kDa	79902	–	–	-0.44	0.000063	–	–	-0.40	0.000090	–	–	–	–	–	–	–	–
ORC3	origin recognition complex, subunit 3	23595	-0.41	0.000029	-0.56	0.000002	-0.33	0.000095	-0.43	0.000010	–	–	–	–	–	–	–	–
ORC6	origin recognition complex, subunit 6	23594	-0.69	0.000014	-0.55	0.000050	–	–	–	–	–	–	–	–	–	–	–	–
PAFAH1B1	platelet-activating factor acetylhydrolase 1b, regulatory subunit 1 (45kDa)	5048	-0.50	0.000266	-0.59	0.000049	-0.66	0.000021	-0.84	0.000002	–	–	–	–	–	–	–	–
PER2	period homolog 2 (Drosophila)	8864	–	–	–	–	–	–	–	–	-0.75	0.000033	-0.72	0.000035	–	–	–	–
PKMYT1	protein kinase, membrane associated tyrosine/threonine 1	9088	-0.56	0.000106	–	–	-0.50	0.000172	–	–	–	–	–	–	–	–	–	–
PLK4	polo-like kinase 4	10733	–	–	–	–	–	–	-0.55	0.000078	–	–	–	–	–	–	–	–
PMF1	polyamine-modulated factor 1	11243	-0.44	0.000063	-0.43	0.000048	-0.55	0.000007	-0.59	0.000003	-0.40	0.000179	-0.40	0.000149	-0.36	0.000454	-0.40	0.000170
POLE	polymerase (DNA directed), epsilon	5426	-0.24	0.000280	-0.28	0.000055	-0.26	0.000088	-0.26	0.000084	–	–	–	–	–	–	–	–
PRIM2	primase, DNA, polypeptide 2 (58kDa)	5558	–	–	0.46	0.000272	–	–	–	–	–	–	–	–	–	–	–	–
PRKAR2B	protein kinase, cAMP-dependent, regulatory, type II, β	5577	–	–	0.57	0.000191	0.93	0.000005	0.94	0.000003	–	–	–	–	–	–	–	–
PSMB10	proteasome (prosome, macropain) subunit, β type, 10	5699	-0.88	0.000368	-0.93	0.000171	–	–	-1.02	0.000061	–	–	–	–	–	–	–	–
PSMC1	proteasome (prosome, macropain) 26S subunit, ATPase, 1	5700	–	–	-0.71	0.000221	–	–	–	–	–	–	–	–	–	–	–	–
PSMC2	proteasome (prosome, macropain) 26S subunit, ATPase, 2	5701	–	–	–	–	-0.31	0.000257	-0.43	0.000014	–	–	–	–	–	–	–	–
PSMC6	proteasome (prosome, macropain) 26S subunit, ATPase, 6	5706	–	–	–	–	–	–	-0.45	0.000076	–	–	–	–	–	–	–	–
PSMD14	proteasome (prosome, macropain) 26S subunit, non-ATPase, 14	10213	–	–	-0.33	0.000083	–	–	-0.31	0.000097	–	–	–	–	–	–	–	–
PSMD2	proteasome (prosome, macropain) 26S subunit, non-ATPase, 2	5708	–	–	1.02	0.000028	–	–	–	–	–	–	–	–	–	–	–	–
PTTG1	pituitary tumor-transforming 1	9232	-0.64	0.000264	-0.80	0.000035	–	–	–	–	–	–	–	–	–	–	–	–
RAD21	RAD21 homolog (S. pombe)	5885	-0.37	0.000336	-0.38	0.000199	–	–	–	–	–	–	–	–	–	–	–	–
RPA2	replication protein A2, 32kDa	6118	-0.44	0.000217	-0.52	0.000046	–	–	–	–	–	–	–	–	–	–	–	–
SEH1L	SEH1-like (S. cerevisiae)	81929	-0.69	0.000006	-0.64	0.000006	-0.70	0.000003	-0.70	0.000002	-0.68	0.000008	-0.63	0.000014	-0.86	0.000003	-0.86	0.000002
SGOL2	shugoshin-like 2 (S. pombe)	151246	–	–	-0.39	0.000175	–	–	–	–	–	–	–	–	–	–	–	–
SKP2	S-phase kinase-associated protein 2, E3 ubiquitin protein ligase	6502	-0.34	0.000178	-0.32	0.000171	-0.45	0.000015	-0.45	0.000010	–	–	–	–	–	–	–	–
SPC25	SPC25, NDC80 kinetochore complex component, homolog (S. cerevisiae)	57405	-0.62	0.000155	–	–	–	–	-0.56	0.000172	–	–	–	–	–	–	–	–
TUBB	tubulin, β class I	203068	–	–	-0.63	0.000191	–	–	–	–	–	–	–	–	–	–	–	–
TUBB4A	tubulin, β 4A class IVa	10382	-0.55	0.000235	-0.62	0.000069	–	–	–	–	–	–	–	–	–	–	–	–
TUBB4B	tubulin, β 4B class IVb	10383	-0.47	0.000010	-0.44	0.000010	–	–	–	–	–	–	–	–	–	–	–	–
TUBGCP4	tubulin, gamma complex associated protein 4	27229	-0.35	0.000151	-0.39	0.000049	-0.41	0.000033	-0.50	0.000005	–	–	-0.30	0.000537	–	–	–	–
TYMS	thymidylate synthetase	7298	-1.16	0.000001	-1.22	0.000001	-0.90	0.000004	-0.92	0.000002	–	–	–	–	–	–	–	–
UBE2C	ubiquitin-conjugating enzyme E2C	11065	-0.45	0.000256	-0.59	0.000055	–	–	–	–	–	–	0.43	0.000178	–	–	–	–
ZWINT	ZW10 interactor	11130	–	–	-0.45	0.000216	–	–	–	–	–	–	–	–	–	–	–	–

Table displaying genes from the “Cell Cycle, Mitotic” FGS (Pathway Commons Reactome database) that were significantly differentially expressed (Log-odds > 3, FDR < 0.007) in DU-145 and PC3 cell lines upon HDACi-treatment with VPA or vorinostat (SAHA). In DU-145 cells most gene changes involved a downregulation of the genes; in PC3 cells only few genes proved to be differentially expressed and less FGS related to mitosis were altered. P-values per gene can be found in the Supplemental Material and at http://luigimarchionni.org/HDACIs.html.

Without proper kinetochore formation, chromosomes will not divide equally across the two daughter cells. To prevent aneuploidy, cells activate the spindle assembly checkpoint, which prevents progression through mitosis if metaphase is not completed properly.³⁵ We hypothesized that treatment with HDACis VPA or vorinostat could result in a functional spindle checkpoint defect by downregulation of kinetochore proteins and other proteins involved in this checkpoint. In PC3 cells we expected the spindle checkpoint to be intact, as there is less transcriptional downregulation of mitotic spindle checkpoint and kinetochore genes. To test this hypothesis, DU-145 and PC3 cells were synchronized in S-phase by a double thymidine block (Fig. 3), released and subsequently treated with colcemid—a spindle poison that depolymerizes microtubules and limits microtubule formation in mitosis. Cells with an intact spindle checkpoint will recognize a defect in microtubule attachment through kinetochore signaling and be arrested in mitosis. In the absence of a functioning kinetochore, the spindle checkpoint is not intact and cells progress through mitosis in the absence of proper chromosome segregation by microtubules. Both PC3 and DU-145 cells accumulated in G2/M-phase after 18h of treatment with colcemid; a small sub-G0 population was measured as well. Therefore, it can be assumed that these cells have an intact spindle assembly checkpoint (Fig. 4A). Separate doses of colcemid were chosen for DU-145 and PC3 cells since the standard dose of 0.1µg/ml led to toxicity in DU-145 cells (Fig. 4A). Next, we treated synchronized PC3 and DU-145 cells after release from S-phase with SAHA or simultaneously with HDACis and colcemid. During the 20 h after administration of 9 µM SAHA, most PC3 cells were initially in the G2/M phase, but later in G1-phase as well. Treatment of DU-145 cells with SAHA resulted in an increased G1- and (later) sub-G0 population. After combination treatment, PC3 cells mainly accumulated in G2/M-phase after ≥ 9 h of treatment, most likely these cells are mitotically arrested by activation of the spindle checkpoint (Fig. 4B; Fig. S20, Supplemental Material and http://luigimarchionni.org/HDACIs.html). On the other hand, treatment of DU-145 cells for ≥ 9 h resulted in a time-dependent increase in a large sub-G0 population apart from a population in G2/M phase, the latter population getting smaller after 20 h of treatment (Fig. 4C; Fig. S21 and http://luigimarchionni.org/HDACIs.html). This may indicate that part of the HDACi-sensitive cells (DU-145 cells) override the spindle checkpoint, leading to a G0/G1 arrest and eventually apoptosis, while cells that are relatively resistant to HDACis (PC3 cells) are arrested in mitosis by activation of the spindle assembly checkpoint.

Figure 3. Summary of the flow cytometry protocol. Cells were synchronized with a double thymidine block, which arrests cells at the beginning of the S-phase. For this purpose, cells were plated, and at 30–50% confluency incubated for 12 h in DMEM-10 media containing 2 mM thymidine. Subsequently cells were washed twice in PBS and incubated in DMEM-10. After 16 h, the media was replaced with DMEM-10/thymidine and cells were incubated for 12–14 h. Then cells were released from S-phase by removing the thymidine containing media. After three hours cells were treated with HDACis and/or colcemid, and for 17 h control and treated cells were harvested (†) every 2–3 h to monitor cell cycle progression by flow cytometry.

Figure 4. HDACis and spindle checkpoint activation in PCa cells. Simultaneous inhibition of PCa cells with HDACis and colcemid leads to accumulation of cells in mitosis in PC3 cells, but to a time-dependent increase of a sub-G0 population in DU-145 cells. (A) DU-145 and PC3 cells were treated for 18 h with colcemid at a final concentration of 0.1 µg/ml after being released from S-phase (M3). This resulted in mitotic accumulation of PC3 cells (M2), but DU-145 cells exhibited an increased sub-G0 population (M4), indicating cytotoxicity. DU-145 cells were successfully arrested in mitosis with low toxicity at a final concentration of 0.04 µg/ml colcemid. (B) PC3 cells were treated for variable periods with either 9 µM vorinostat alone or a combination of 0.1 µg/ml colcemid with 4 µM vorinostat or 4 mM VPA. Combining HDACis with colcemid resulted in mitotic accumulation of PC3 cells. (C) DU-145 cells were treated for variable periods with either 9 µM vorinostat alone or a combination of 0.04 µg/ml colcemid with 4 µM vorinostat or 4 mM VPA. Combining HDACis with colcemid resulted in a time-dependent increase of a sub-G0 population in DU-145 cells.

Discussion

In the past decade, a growing number of scientists have concluded that single agents, often targeting a single pathway, are insufficiently effective for cancer treatment, as many complex signaling pathways are involved in cancer development and progression. Therefore clinical trials have shifted focus to combination therapy to either target different pathways or target the same major pathway from different angles. However, many recent phase III studies have concluded that clinical benefit of these combination therapies is limited, while toxicity is increased.^36-38 To speed up the implementation of combination therapies in clinic and to reduce costs of failing expensive clinical trials, there is a strong need to develop preclinical tools to identify combination strategies that will be successful in clinic.

In this study, we show that AFA allows unbiased hypothesis generation for combination therapies. This systems biology approach successfully unraveled the mechanisms of action of HDACis, without the bias of assumptions based on previous literature or the preference of the researcher. We show major FGS expression changes after HDACi-treatment, induced either directly through regulation of histones or indirectly by regulating (multiple) genes in pathways. A major advantage of AFA is that it is able to identify and visualize modest gene expression changes when occurring across a predefined category, without prior arbitrary cut-offs for fold change, p-value, and/or false discovery rates (FDR). This is particularly relevant since previous studies have indicated that fluctuations in such basal gene expression patterns at levels less than what is considered significant (i.e., FDR < 0.05) do have a functional relationship.³⁹ Furthermore, AFA combines biological information from different knowledge domains and different conditions and is therefore equal to GSEA done over multiple contrasts. FGS representing similar biological concepts proved to be in agreement with each other, supporting the strength of AFA. In summary, AFA enables us to identify specific biological and molecular processes affected by chemotherapeutic agents in greater detail, as subtle but consistent changes can be identified, providing a framework for further investigation.

We hence performed differential gene expression analysis and applied AFA to four distinct data sets, ultimately analyzing over 42,475 distinct FGS over a total of 375 different samples. This analysis, which we made publically available online, is the largest performed on PCa cells treated with HDACis. It therefore constitutes an invaluable resource for the identification of FGS that are targeted by HDACis. Finally, we have investigated two of the identified molecular concepts for which multiple FGS proved affected upon HDAC-inhibition more into detail, namely, MHC genes and mitotic spindle checkpoint genes.

The immune system plays a prominent role in the development, progression and treatment of cancers, among others PCa.⁴⁰ Recognition of cancer cells by the immune system leads to eradication of these cells, therefore cancer cells need to escape this immune response (immune escape). This is accomplished by both genetic mutations and epigenetic modifications resulting in the loss of MHC class I and II molecules in some cancer cells (immuno-editing); these cells then get selected for further proliferation.⁴⁰ Indeed, we show in this study that in human primary PCa tissue samples MHC-related gene expression is deregulated. This is in line with earlier reports: for instance Sharpe et al. concluded that expression of MHC class I and II determinants (such as B2M) was downregulated in tissue samples from PCa patients.⁴¹ While B2M expression is downregulated, HLA-A expression does not differ significantly in our results. Tamura et al. concluded that expression of HLA-A is downregulated when PCa cells progress from hormone therapy sensitive to hormone resistant PCa.⁴² Downregulation of HLA-A may not have occurred in our patients’ samples yet, as most cancers are in an early stage when a radical prostatectomy is performed.

Immunotherapy has become a common approach to treat cancer. It can be used to target specific pathways, such as the use of trastuzumab to target HER2 positive breast cancer cells,⁴³ or to modulate the immune system against cancer cells.The role of immunotherapy in PCa treatment is less well established. Although it is generally accepted that the immune system plays a role in PCa, the exact mechanism by which the immune system is involved in PCa initiation, progression and treatment remains unclear.^41,44 Sipuleucel-T has been approved by the US FDA as a cellular immunotherapy that improves overall survival in asymptomatic or minimally symptomatic metastatic castrate-resistant PCa patients, and other immunotherapies such as ipilimumab are in advanced stages of clinical development.^45,46 However, some have doubted its efficacy.^47,48 Therefore, immunotherapy is not yet widely used in PCa treatment.

As VPA and vorinostat upregulate FGS involved in the immune system in DU-145 and PC3 cells, treatment by HDACis may prompt a natural or induced immune response against tumor cells by the immune system or immunotherapy, respectively. Studies in other cancer types have also found that HDACis enhance the expression of immunologically important molecules and induce immune responses, strengthening our findings.⁴⁹ As there currently is no tissue available from patients after treatment with HDACis, it is not possible to assess whether such immune responses occur in PCa patients. However, our results indicate that HDACis may revert the acquired B2M deficiency found in cancer cells resistant to immunotherapy.⁵⁰ Considering the increasing importance of immunotherapy in PCa, it is imperative that future studies test the hypothesis that pretreatment with HDACis may sensitize PCa cells for immunotherapy. The clinical response to sipuleucel-T, ipilimumab and/or other immunotherapies may be increased through reversal of (epi)genetic modifications that lead to immune escape of the cancer cell.

Results from our study further indicate that posttranslational changes may not be solely responsible for the regulatory effects of HDACis during mitosis, but translational downregulation of genes involved in mitotic regulation may also play a role in creating mitotic defects. Treatment with high doses of SAHA seems to result in progression through mitosis in the absence of a mitotic spindle in HDACi-sensitive cells. This may lead to aneuploidy, a G0/G1 phase arrest and eventually to cell death.⁵¹ Similar results have been reported by Noh et al.,⁵² in which treatment of HeLa cells with Trichostatin A resulted in a transcription-dependent defective mitosis with loss of checkpoint control. This effect may be increased when treating PCa cells with agents that deregulate the spindle assembly checkpoint. Combination therapy of HDACis with one such spindle checkpoint deregulator, aurora kinase inhibitors, has recently been tested in lymphoma cells in vitro and in vivo, with success.⁵³ Both results from this study and our results encourage testing combinations of HDACis with aurora kinase inhibitors or other agents that deregulate the spindle checkpoint, such as polo-like kinase 1 (Plk1-)inhibitors, in PCa and other solid tumors. Preliminary results of preclinical studies by our group involving combination therapies with HDACis (VPA/vorinostat) and an aurora kinase inhibitor (AMG 900) or Plk1-inhibitors (BI 2536/volasertib [BI 6727]) in PCa indicate that these combinations result in increased anti-tumor activity with acceptable toxicities.^54,55

Our AFA results identified other FGS differentially expressed upon HDACi-treatment that may be of clinical importance too, such as upregulation of the androgen receptor signaling pathway. Increased sensitivity to androgen therapy after VPA treatment has been reported with other PCa cell lines.⁵⁶ Our data suggest that DU-145 and PC3 cells may have increased sensitivity to androgen therapy after HDACi-treatment as well, potentially by upregulation of androgen receptor expression.⁸ Further studies need to be performed to assess increased androgen sensitivity after HDACi-treatment in PCa.

In summary, in this study we successfully applied AFA. This created an expansive and complete library of biological processes modulated by HDACis in PCa cell lines, and is therefore a valuable resource for researchers involved in research with HDACis. We further explored two pathways that were differentially expressed upon HDACi-treatment that may be of interest for future combination therapy. Novel studies combining HDACis with inhibitors of identified pathways are granted to assess clinical benefit for PCa patients. Furthermore, this study indicates that AFA can be used as an unbiased rational approach to identify novel combination strategies against cancer. Performing such analyses generates a comprehensive list of FGS differentially expressed upon cancer cell treatment.

Materials and Methods

Ethics statement

Anonymized human PCa specimens were obtained from the University Medical Center Utrecht (UMCU), after it was ensured that enough material remained to serve the patient’s and family’s needs, and in accordance to the Dutch code of conduct for the use of leftover body material.⁵⁷ All investigations involving human samples were performed in strict adherence to the Declaration of Helsinki.

Cell lines and chemicals

PCa cells were obtained from the American Type Culture Collection (ATCC) and maintained as described by our group previously.¹⁵ Stock solutions of 10 mM vorinostat (AtonPharma/Merck) dissolved in DMSO were used for all experiments; vorinostat was further diluted in complete RPMI-1640 media (Invitrogen). VPA (Sigma-Aldrich) was stored as a salt, and freshly dissolved in complete RPMI-1640 media (RPMI-1640 media supplemented with 10% fetal bovine serum [FBS]) before each experiment. Colcemid (Sigma-Aldrich) was acquired as a 10 µg/ml solution in Hank’s balanced salt solution (HBSS) and stored at 4°C.

Microarray experiment

Microarray data were obtained from time course experiments in which PC3 and DU-145 cells were treated with VPA or vorinostat for 48 h, and are available from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (series GSE34452). Expression data were processed using the R-Bioconductor library limma.^58-63 Differential gene expression data were acquired using a generalized linear model approach (“multiple-loop, double-cube” design) as described by Kortenhorst et al.¹⁵ A detailed explanation of all procedures and methods used for microarray data pre-processing, and differential gene expression analysis and detection has been described in detail previously.^15,64 Briefly, in this previous study we applied a complex design based on a “multiple-loop, double-cube” cDNA microarray experiment to enhance our understanding of the molecular underpinnings of HDACi resistance in PCa cells. To this end, we analyzed gene expression data from a total of 22 dual-color microarray hybridizations comparing DU-145 and PC3 cells treated with vorinostat and VPA after an incubation of 48 and 96 h. Data were normalized using within-array “loess” normalization, and an ANOVA coupled with empirical Bayes standard errors shrinkage was used to identify differentially expressed genes between treated cells and controls. HDAC inhibition resulted in differential expression of 2.8% to 10% of genes (adjusted P-value < 0.001) across all conditions. Between 51% to 72% of such genes were upregulated, while 28% to 49% were downregulated. Larger differential gene expression was observed in DU-145 cells compared with PC3 cells, with VPA treatment compared with vorinostat, and with longer drug exposures.

Other public domain data were obtained from GEO and the Connectivity Map, publicly available from the NCBI GEO database, along with MIAME (minimum information about a microarray experiment) compliant information. Details regarding the source of gene expression data sets used, the platforms and annotation packages employed, the pre-processing procedures adopted for each data set, as well as the statistical details of the analysis performed, are described extensively in the Supplemental Material and at http://luigimarchionni.org/HDACIs.html.

Analysis of Functional Annotation (AFA)

We applied Analysis of Functional Annotation (AFA), a gene set analysis approach, on differential gene expression data obtained from distinct comparisons and across different studies to identify biological processes and signaling pathways modulated by HDACis in PCa cells.^{10,11,18,19,65} Overall, this methodology extends gene set analysis procedures, such as gene set enrichment analysis (GSEA) or parametric analysis of gene set enrichment (PAGE), by investigating biological processes enrichment over multiple experimental conditions as briefly summarized below.^12-14

FGS, recapitulating distinct and complementary biological concepts such as cellular signaling pathways, PPI networks, downstream transcriptional responses, gene expression regulatory networks orchestrated by transcription factors and microRNA targets, were retrieved in the form of gene lists from various publicly available databases (see Table S1 and http://luigimarchionni.org/HDACIs.html). These collections included the Reactome, the HPRD, GO, KEGG, the MSigDB and NCBI Entrez Gene databases.^14,21-29

A one-sided Wilcoxon rank sum test was separately applied across all investigated comparisons to test whether any given FGS was differentially expressed, upregulated or downregulated, using the absolute and signed t-statistics to order genes (for details on the linear model analysis see the Supplemental Material and Kortenhorst et al.¹⁵). The enrichment analysis was performed on all non-redundant genes present on the microarray, according to the NCBI Entrez Gene database annotation.⁶⁶ Filtering of redundant microarray features (i.e., probes mapping to the same NCBI Entrez Gene identifier) was achieved by retaining only the probes with the largest absolute t-statistics for further analysis. Correction for multiple hypothesis testing was obtained separately for each FGS collection by applying the Benjamini and Hochberg method.⁶⁷ Differentially expressed FGS were visualized using heatmaps; an adjusted p-value < 0.05 was considered significant.

To validate the results from the AFA on our microarray data, we further performed differential gene expression analysis and AFA on publicly available data sets of HDACi-treated PCa cells. Three data sets were available (GSE8645, GSE31620, and Connectivity Map). In one data set (GSE8645), LNCaP cells had been treated with either 7.5 µM CG-1521 or 5 µM Trichostatin A (both HDACis) for 24 h, after which a microarray was performed.⁶⁸ In data set GSE31620, LNCaP cells were treated for 36 h with Trichostatin A and/or the DNA-methylating agent 5-Azacytidine after which microRNA microarrays were performed.¹⁶ In the Connectivity Map PC3 cells had been treated with various HDACis at various dosages for 6 h.¹⁷ A detailed description of these data sets is available in the Supplemental Material and http://luigimarchionni.org/HDACIs.html.

Flow cytometry

DU-145 and PC3 cells were synchronized in S-phase by a double thymidine block (Fig. 3). Cells were plated in 100 mm dishes; at 30–50% confluency, cells were incubated in DMEM-10/thymidine media (DMEM [Invitrogen] supplemented with 10% FBS and 2 mM thymidine [Sigma-Aldrich]) for 12h. Subsequently cells were washed twice in PBS and incubated in DMEM-10 media. After 16h, the media was replaced with DMEM-10/thymidine and cells were incubated for 12–14h. Cells were released from S-phase by removing the thymidine-containing media; three hours later DMEM-10 media containing an HDACi and/or colcemid was added. Cells were harvested immediately after the S-phase release, at the start of treatment, and 3, 6, 9, 11, 14 or 17h after treatment initiation.

Cells were prepared for flow cytometry analysis using the method of Vindelov.⁶⁹ Flow cytometry was performed on a BD-LSR II flow cytometer (BD Biosciences); data were interpreted using CellQuest (BD Biosciences).

Tissue microarrays (TMAs) and immunohistochemistry (IHC)

For TMAs, human PCa samples were obtained from archives of the Pathology Department of the UMCU. Samples from 71 patients who had undergone radical prostatectomy were collected; information concerning survival and PSA recurrence of these patients was collected. Tumors had been formalin-fixed and paraffin-embedded after resection.

From paraffin blocks four 1 mm cores were taken per specimen. These cores were transferred to a recipient composite paraffin block using a Beecher MT1 manual arrayer (Beecher Instruments). Four micron slides were prepared from each TMA. Spots that contained necrotic/dead tissue or fibrin instead of cancer cells were excluded.

IHC was performed with 4 µm paraffin sections on the BOND-MAX (Leica) using the Bond Polymer Refine Detection™ (Leica) kit. With the exception of B2M, antigen retrieval was performed in citrate buffer (pH 6.0) for 20 min at 100 °C. Sections were then incubated with anti-HLA-A antibody HCA2 (HCE kweek, UMCU) 1:200, anti-HLA-B antibody HC10 (HCE kweek) 1:400, anti-B2M (Dako) 1:600, anti-HLA-DR-α (Dako) 1:60 or anti-HLA-DR-β (Dako) 1:800, consecutively followed by incubation with Bond Polymer Refine Detection kit items (Leica) and incubation with DAB solution (Bond Polymer Refine Detection, Leica). Slides were counterstained with hematoxylin (Bond Polymer Refine Detection, Leica), and blindly examined and scored for protein expression by two investigators (MSQK, PJD). Appropriate negative and positive controls were used.

Abbreviations:
AFA	=	analysis of functional annotation
ANOVA	=	analysis of variance
ATCC	=	American type culture collection
B2M	=	β-2 microglobulin
CAT	=	correspondence at the top
cDNA	=	complementary deoxyribonucleic acid
CREB	=	cAMP response element binding protein
CTCL	=	cutaneous T-cell lymphoma
DAB	=	3,3′-Diaminobenzidine
DMEM	=	Dulbecco’s modified eagle medium
DMSO	=	dimethyl sulfoxide
DNA	=	deoxyribonucleic acid
EGFR1	=	epidermal growth factor receptor 1
FBS	=	fetal bovine serum
FDR	=	false discovery rate
FGS	=	functional gene set
GEO	=	gene expression omnibus
GO	=	Gene Ontology
GSEA	=	Gene Set enrichment analysis
HBSS	=	Hank’s balanced salt solution
HDAC	=	histone deacetylase
HDACi	=	histone deacetylase inhibitor
HER2	=	human epidermal growth factor receptor 2
HLA	=	human leukocyte antigen
HPRD	=	human protein reference database
IHC	=	immunohistochemistry
ISRE	=	interferon-stimulated response element
KEGG	=	Kyoto Encyclopedia of Genes and Genomes
LILRB	=	Leukocyte immunoglobulin-like receptor B
MHC	=	major histocompatibility complex
MIAME	=	minimum information about a microarray experiment
MSigDB	=	molecular signatures database
mTOR	=	mammalian target of rapamycin
NCBI	=	National Center for Biotechnology Information
PCa	=	prostate cancer
Plk1	=	polo-like kinase 1
PPI	=	protein-protein interaction
PSA	=	prostate-specific antigen
RPMI	=	Roswell Park Memorial Institute medium
SAHA	=	suberoylanilide hydroxamic acid
TMA	=	tissue microarray
UMCU	=	University Medical Center Utrecht
US FDA	=	United States Food and Drug Administration
VPA	=	valproic acid

Acknowledgments

The authors would like to thank Marianna Zahurak for her contribution to the statistical analysis for this project and René Medema for his support and discussion. This study was supported by the Prostate Cancer Foundation, the Flight Attendant Medical Research Institute (SK), NIH (National Institute of Health) - NCI (National Cancer Institute) SPORE (Specialized Programs of Research Excellence) in Prostate Cancer P50CA058236 (MAC), by NIH-NCI 5P30CA006973 (MAC, LM, SK), and Aegon International Research Fellowship in Oncology (MAC, SK). MSQK was also partially supported by a fellowship from the Catharina van Tussebroek Fund. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Supplemental Materials

Supplemental materials may be found here: http://www.landesbioscience.com/journals/epigenetics/article/25574