The Effect of Vorinostat on the Development of Resistance to Doxorubicin in Neuroblastoma

Histone deacetylase (HDAC) inhibitors, especially vorinostat, are currently under investigation as potential adjuncts in the treatment of neuroblastoma. The effect of vorinostat co-treatment on the development of resistance to other chemotherapeutic agents is unknown. In the present study, we treated two human neuroblastoma cell lines [SK-N-SH and SK-N-Be(2)C] with progressively increasing doses of doxorubicin under two conditions: with and without vorinsotat co-therapy. The resultant doxorubicin-resistant (DoxR) and vorinostat-treated doxorubicin resistant (DoxR-v) cells were equally resistant to doxorubicin despite significantly lower P-glycoprotein expression in the DoxR-v cells. Whole genome analysis was performed using the Ilumina Human HT-12 v4 Expression Beadchip to identify genes with differential expression unique to the DoxR-v cells. We uncovered a number of genes whose differential expression in the DoxR-v cells might contribute to their resistant phenotype, including hypoxia inducible factor-2. Finally, we used Gene Ontology to categorize the biological functions of the differentially expressed genes unique to the DoxR-v cells and found that genes involved in cellular metabolism were especially affected.


Introduction
Neuroblastoma is the most common extra-cranial solid organ malignancy of childhood. Outcome is heterogeneous and depends on several clinical and biologic factors. Five year survival approaches 90% for infants diagnosed with neuroblastoma in the first year of life, but remains only 65% among all children older than one year-of-age and 30-60% for those with high-risk tumors [1][2][3]. Acquired resistance to commonly used chemotherapeutic agents remains a major barrier to successful therapy, especially among older children with high-risk neuroblastoma. Doxorubicin, along with cisplatin, cyclophosphamide and etoposide, are key components of modern chemotherapy protocols for intermediate and high-risk neuroblastoma [2,4]. The development of resistance to these commonly used agents can be associated with upregulation of multidrug transporter genes such as mdr1 and mrp1. Upregulation of the genes encoding these drug efflux pumps has been correlated with a poor clinical prognosis [5,6]. Novel strategies are therefore required to prevent or reverse the development of the multidrug resistant phenotype in order to improve the prognosis for children with high-risk, relapsed and recurrent neuroblastoma.
Vorinostat is a class I and II histone deacetylase (HDAC) inhibitor which has demonstrated efficacy against neuroblastoma in preclinical studies [7][8][9]. Vorinostat was well-tolerated in a phase I study of children with solid organ tumors including neuroblastoma, and demonstrated promising early results [10]. Additional phase I trials of novel protocols incorporating vorinostat are currently ongoing [4]. However, the long-term effect of vorinostat treatment on the expression of genes encoding multidrug efflux pumps and the development of resistance to commonly used chemotherapeutic agents is unknown.
In the current study, we investigated the effect of vorinostat cotreatment on the development of drug resistance using two human neuroblastoma cell lines (SK-N-SH and SK-N-Be(2)C) exposed to progressively increasing concentrations of doxorubicin. Our aims were to determine (1) the functional effect of vorinostat on the long term responsiveness of these cells to doxorubicin, and (2) the genetic changes induced by vorinostat during this process.

Cell Lines, Reagents and Proliferation Assay
Human neuroblastoma (SK-N-SH and SK-N-Be(2) C) lines were purchased from American Type Culture Collection (Rockville, MA). DMEM and fetal bovine serum (FBS) were obtained from Mediatech (Herndon, VA) and Atlanta Biologicals (Atlanta, GA), respectively. Doxorubicin, 3-(4,5-dimethyl-2-thiazolyl)2,5diphenyl tetrazolium bromide (MTT), and deferoxamine mesylate (DFO) were purchased from Sigma (St. Louis, MO). Vorinostat was purchased from ChemieTek (Indianapolis, IN). Cell proliferation was assessed by MTT (2-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) assay. Cells were seeded in 96 well plates and incubated with logarithmic concentrations of doxorubicin ranging from 10 29 to 10 25 M. MTT (10 mL of 5 mg/ml solution) was added to each well of the titration plate and incubated for 4 hours at 37uC. The cells were then solubilized by the addition of 100 mL of 10% SDS/0.01 mmol/L HCL and incubated for 15 hours at 37uC. The absorbance of each well was determined in an ELISA plated reader using an activation wavelength of 570 nm and a reference wavelength of 650 nm. Cell viability in the presence of different doses of doxorubicin was determined by comparison with untreated control cells.
Quantitative Real-time PCR 1 mg of RNA extracted from each cell lines as previous described was reverse transcribed to cDNA in a 20 ml reaction by using GeneAmpH RNA PCR Core Kit (Applied Biosystems, Foster City, CA) with random primers. Quantitative RT-PCR assays were then performed in triplicate on 1 mL of each cDNA dilution using the Power SYBR Green Master Mix (Applied Biosystems, Foster City, CA), the following thermal cycling specifications were performed on the ABI 7500 Fast Real-Time PCR system (Applied Biosystems, Foster City, CA) 60 s at 95uC and 40 cycles each for 9 s at 95uC and 60 s at 60uC. Specific intron-spanning primers were designed for EPAS1 (forward primer: CATTTGAGTCCTACCTGCTGC and reverse primer: GTAGAAGCCTGGCTCAGGTG), SST (forward primer: TTGTCCTCCCCACTTCTCT and reverse primer: AAA-GATTTACTCAATAGAAACCA) and Pon2 (forward primer: AACAGCCATAGTAGTCAC and reverse primer: CAAGGGA-CAGAAAAGAAAG). The amplification of Glyceraldehyde-3phosphate dehydrogenase (GAPDH) (forward primer: GGAGTGTGCCCGGGATGAAGAAT and reverse primer: TGGGGCTGGCAGGCTAAACA) was selected as the endogenous control to normalize all data, relative expression values were obtained by compared to the parental lines' value.

Whole Genome Expression
RNA was isolated from triplicate specimens of each cell line using the RNeasy Mini Kit (Qiagen, Venlo, Netherlands). RNA quality and concentration were determined using the 2100 Bioanylzer (Agilent Technologies, Santa Clara, CA) and NanoDropH ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE). Samples with RNA integrity numbers greater than 9.5 were used in subsequent expression analyses. All RNA samples were stored at 280uC prior to labeling and hybridization.
Amplification and labeling of the RNA was performed with the Illumina TotalPrep RNA Amplification Kit from Life Technologies (Ambion,Austin, TX) using 150 ng of RNA per sample. 750 ng of labeled cRNA were hybridized to Human HT-12 v4 Expression BeadChips for 18 hours and processed (Illumina, San Diego, CA) according to the manufacturer's instructions. Slides were scanned using the Illumina iScan System Human HT-12 v4 beadchip containing 47,231 probes.

Image and Data Analysis
Quality check and probe level processing of the Ilumina microarray data were further made with R bioconductor package, Lumi [11]. The data processing includes a normalization method [12] to reduce the obscuring variation between microarrays, which might be introduced during the processes of sample preparation, manufacturing, fluorescence labeling, hybridization and/or scanning. Hierarchical clustering and Principal Component Analysis were performed on the normalized signal data to assess the sample relationship and variability. Probes that were expressed in none of the samples as judged by the criteria of Illumina detection p-value $0.01 were filtered out. Differential gene expression between the different conditions was assessed by statistical linear model analysis using bioconductor package, limma, in which an empirical Bayes method is used to moderate the standard errors of the estimated log-fold changes of gene expression. It results in more stable inference and improved power, especially for experiments with small numbers of microarrays [13,14]. The moderated t statistic pvalues derived from the limma analysis above were further adjusted for multiple testing by Benjamini and Hochberg's method to control false discovery rate (FDR). The list of differentially expressed genes was obtained by the FDR criteria of ,10% and foldchange cutoff of .1.5, and visualized by volcanoplots. The gene ontology analysis of the lists of differentially expressed genes were performed with R bioconductor package, topGO [15]. All bioconductor packages are available at http://bioconductor.org and all computation was performed under R environment (http:// www.r-project.org).

Development of Drug Resistant Cells
The half maximal inhibitory concentration (IC 50 ) of doxorubicin was determined by MTT assay. The IC 50 of doxorubicin for the parental SK-N-SH and SK-N-Be (2) Resistance to doxorubicin was confirmed by MTT assay as shown in Figure 1 (a-b). The DoxR and DoxR-v cells were both equally resistant to doxorubicin, as well as cross-resistant to etoposide ( Figure 1c) but not cisplatin (data not shown). Vorinostat treatment alone (WT-v cell lines) did not cause doxorubicin resistance.
P-glycoprotein (P-gp) expression was assessed by western blot (Figure 1d). P-gp, the protein product of the ABCB1 or multidrug resistance mdr1 gene was significantly upregulated in both the SK-N-SH and SK-N-Be(2)C DoxR cells relative to the WT or WT-V cells. Interestingly, although the DoxR-v cells were equally resistant to doxorubicin, these cells had reduced P-gp expression compared to their DoxR counterparts. Likewise in the SK-N-Be(2)C cell line, hypoxia inducible factor-1a (HIF-1a) expression was similarly upregulated in the Dox-R but not the DoxR-v cells (Figure 1e) under hypoxia mimicking conditions.

Whole Genome Expression
Whole genome expression analysis was performed using the WT, WT-v, DoxR, and DoxR-v cells from the SK-N-SH and SK-  N-Be(2)C lines. To understand the effect of vorinostat treatment on the development of drug resistance, we first analyzed the expression of known drug resistance genes in the DoxR and DoxR-v cells relative to their parental (WT) lines (Table 1). In agreement with the aforementioned western blot findings, expression of the ABCB1 gene (which encodes the P-gp protein) was approximately 4-fold elevated in the DoxR cells, but only 2fold elevated in the DoxR-v cells. Of note, expression of the ABCC1 (mrp1) gene was not upregulated in the DoxR nor the DoxR-v cells relative to their parental line. No other known drug resistance genes had significant differential expression between the DoxR-v and DoxR cells to explain how the DoxR-v cells were equally resistant to doxorubicin despite having lower P-gp expression. Genes differentially expressed in the DoxR-v cells relative to the parental (WT) lines were analyzed in an effort to determine  Table 2. Genes with a potential role in doxorubicin-resistance following vorinostat treatment.

Gene Symbol
Gene Entrez ID Gene Name SK-N-SH Fold SK-N-Be (2) alternative pathways for doxorubicin resistance following Vorinostat co-treatment. Volcano plots (Figure 2, a-b), demonstrate a similar pattern of differentially expressed genes (DEGs) in the SK-N-SH and SK-N-Be(2)C DoxR-v lines. DEGs common to both cell lines were compared between the DoxR-v, DoxR, and WT-v lines (Figure 2c). The 405 DEGs unique to the DoxR-v versus WT comparison were selected for further analysis, to hone in on the genetic changes responsible for drug resistance following this specific treatment combination. The 51 genes with the greatest differential expression (fold change .2) are shown in Table 2, and displayed visually as heat maps in Figure 3. The complete list of all 405 DEGs (fold change .1.5, p,0.1) unique to the DoxR-v versus WT comparison are shown in Table S1. Expression of these 405 DEGs in WT and DoxR-v cells for the SK-N-SH and SK-N-Be(2)C lines are displayed visually as heat maps in Figures S1 and S2, respectively. Among these DEGs, 181 (44.7%) were significantly down-regulated in the DoxR-v cells, while the remaining 224 (55.3%) were up-regulated.
To confirm the results of the whole genome analysis, qRT-PCR was performed on selected DRGs. Expression of somatostatin (SST), paraoxanase 2 (PON2) and endothelial PAS domain protein 1 (EPAS1) in the WT-v, DoxR and DoxR-v cells relative to their parental lines is shown in Figure 4. Relative expression of these genes based on microarray analysis and qRT-PCR is compared in Table 3.

Functional Analysis
Functional analysis based on mapping of biological processes from the Gene Ontology Database was performed using the 405 DEGs unique to the DoxR-v versus WT comparison for both cell lines. For 34 different biological processes, a statistically greater than expected number of the genes annotated to that process were differentially expressed (Table 4). These altered biological processes included a predominance of those involved in cellular metabolism, cellular biosynthetic processes, and regulation of the cytoskeleton and   microtubule structure. A hierarchical flowchart of these altered biological processes is included as Figure S3.

Discussion
In the first study of the long-term effect of HDAC inhibition with vorinostat on the development of chemotherapy resistance, we herein demonstrate that vorinostat alters the mechanism but does not prevent the development of multidrug resistance. Cells treated intermittently with vorinostat each time the dose of doxorubicin was increased in generating resistance had significantly reduced expression of the mdr1 gene and its protein product, P-gp, yet still demonstrated equivalent functional resistance to doxorubicin and etoposide. Whole genome analysis did not identify any known multidrug resistance genes that were upregulated in the vorinostat-treated cells to compensate for the lower expression of mdr1. However, several specific genes and biological processes that may play a role in the alternative mechanism of resistance following vorinostat treatment were identified.
Prior studies of histone deacetylase inhibitors, including Vorinostat, in neuroblastoma have focused on their short term effects where they have been shown to induce apoptosis of human neuroblastoma cells in vitro and in vivo [16,17]. As vorinostat becomes more widely utilized in human clinical trials, it is important to understand the long-term effect of treatment with this drug on neuroblastoma responsiveness to standard cytotoxic agents [4,10,18]. Interest in the ability of HDAC inhibitors to prevent or reverse the development of drug resistance stems from findings of elevated HDAC expression in drug resistant neuroblastoma cell lines. Keshelava et al demonstrated that HDAC1 was upregulated in multidrug resistant neuroblastoma cell lines relative to drug-sensitive lines [19]. Inhibition of HDAC1 expression or activity resulted in synergistic killing of multidrug resistant cells when combined with standard cytotoxic agents. Oehme et al found elevated HDAC8 expression in primary human neuroblastoma samples from children with advanced and metastatic disease [20]. They reported that small-molecule inhibition of HDAC8 induced cellular differentiation and reduced proliferation. In both of the above studies, however, the effect of the HDAC inhibitors was studied in neuroblastoma cell lines that already displayed the multidrug resistance phenotype or other high-risk features. Our study is the first to investigate the effect of an HDAC inhibitor during induction of resistance to a commonly used cytotoxic agent by exposure to progressively increasing concentrations of the drug.
Our results provide further proof that expression of the P-gp drug efflux pump is only one component of the multi-drug resistance phenotype, which is also mediated by other drug efflux pumps and aspects of the cellular structure and tumor microenvironment [21]. We demonstrate that co-treatment with vorinostat can reduce expression of the ABCB1 (mdr1) gene and its protein product P-gp, yet these cells are equally resistant to doxorubicin and etoposide. We hypothesize that the effect of vorinostat on P-gp expression may be partially mediated by hypoxia inducible factor 21a (HIF-1a). We and others have demonstrated that histone deacetylase inhibitors, including vorinostat, can reduce HIF-1a activity in tumor cells [22,23]. In turn, HIF-1a has been shown to bind to the promotor region of the mdr1 gene to increase its expression under hypoxic conditions [24]. Therefore, long-term inhibition of HIF-1a activity by vorinostat treatment could partially explain the reduced expression of P-gp. At the same time, however, a recent prospective analysis of neuroblastoma tissue specimens suggested that the level of mdr1 expression had no prognostic significance [6]. The present work confirms that neuroblastoma cells can develop equivalent multidrug resistance despite decreased P-gp expression.
We identified 405 additional genes with differential expression unique to the DoxR-v cells. Although no drug transporter proteins were included in this list, several of the identified genes have been previously shown to play a role in drug resistance. Endothelial PAS domain protein 1 (EPAS1), otherwise known as hypoxiainducible factor 2a has been shown to play a key role in tumor responsiveness to chemotherapy. Elevation of HIF-2 (EPAS1) has been shown to suppress p53 function, and its knockdown can restore p53 function and reverse resistance to chemotherapyinduced cell death [25]. Furthermore, HIF-2 (EPAS1) has been shown to promote hypoxia response element (HRE)-driven gene expression, and shRNA knockdown increased renal cell carcinoma sensitivity to ionizing radiation [26]. If, as postulated above, HIF-1a activity is suppressed by vorinostat, which in turn leads to decreased P-gp expression, then increased HIF-2a (EPAS1) activity may mediate an alternative pathway by which the cells respond to the stress of the chemotherapeutic agents. Prior studies have shown that HIF-1a and HIF-2a activities have a reciprocal relationship: knockdown of HIF-1a induces HIF-2a (EPAS1) upregulation [27]. It is also known that HIF-2 activity is stimulated by Sirtuin 1 (sirt1), a class III HDAC that was upregulated in both our DoxR and DoxR-v cells, and which is not suppressed by vorinostat [28]. Further work is required to identify downstream targets of HIF-2a (EPAS1) which may mediate drug-resistance in vorinostat-treated cells.
A number of additional DEGs identified by our whole genome analysis have been previously postulated to play a role in chemotherapy resistance and other aspects of an aggressive tumor phenotype. Paraoxanose-2 upregulation has been found to correlate with chemotherapy resistance via an antioxidative mechanism which reduces oxidative damage to the endoplasmic reticulum [29]. Stomatin is a membrane protein that may serve as a skeleton anchor and regulate cation transport through lipid membranes, and it has been found upregulated in doxorubicin resistant tumors [30]. Finally, EBF3 has been identified as a tumor suppressor gene which mediates cell cycle arrest and apoptosis and whose expression is down regulated in a number of malignancies [31].
The genes with unique differential expression in the vorinostat co-treated doxorubicin-resistant cells are known to play key roles in a number of biological processes which have been implicated in the drug resistance phenotype. Of particular note, the DoxR-v cells were found to have statistically significant differential expression of genes involved in cellular metabolism, including ketone metabolism (GO:0042180) and lipid metabolism (GO:0006629). Alterations in the metabolic needs of tumor cells are increasingly recognized as a mechanism by which cells survive the stress of chemotherapy. Metabolic alterations in lipid metabolism have been proposed as an advantageous mechanism which helps tumor cells escape chemotherapy-induced apoptosis [32]. In drug-resistant colorectal cancer cells, alterations in lipid metabolism have been demonstrated and are postulated to offer a growth advantage under stress conditions [33]. Alterations in lipid metabolism may also alter drug uptake via changes in the properties of the cell's lipid bilayer [21].
In conclusion, we have shown that intermittent therapy with Vorinostat has long-term effects on the development of resistance to doxorubicin in neuroblastoma cell lines. Vorinostat reduces expression of P-glycoprotein, the protein product of the mdr1 gene, but does not prevent the functional development of multidrug resistance. Whole genome analysis identified several genes, including HIF-2a (EPAS1), whose differential expression might partially explain the equivalent multidrug resistance in the vorinostat treated cells. Further work is needed to investigate these differentially expressed genes, their downstream targets, and the biological process involved to fully delineate the implications of Table 4. Altered biological processes based on the genes with differential expression (fold change.1.5; adjusted p,0.1) unique to the DoxR-v cell lines as mapped by Gene Ontology. vorinostat on the long-term development of multidrug resistance in neuroblastoma.