Distinguishing between cancer cell differentiation and resistance induced by all-trans retinoic acid using transcriptional profiles and functional pathway analysis

All-trans retinoic acid (ATRA) induces differentiation in various cell types and has been investigated extensively for its effective use in cancer prevention and treatment. Relapsed or refractory disease that is resistant to ATRA is a clinically significant problem. To identify the molecular mechanism that bridges ATRA differentiation and resistance in cancer, we selected the multidrug-resistant leukemia cell line HL-60[R] by exposing it to ATRA, followed by sequential increases of one-half log concentration. A cytotoxicity analysis revealed that HL-60[R] cells were highly resistant to ATRA, doxorubicin, and etoposide. A comparative genome hybridization analysis of HL-60[R] cells identified gains of 4q34, 9q12, and 19q13 and a loss of Yq12 compared with in the parental HL-60 cell line. Transcriptional profiles and functional pathway analyses further demonstrated that 7 genes (FEN1, RFC5, EXO1, XRCC5, PARP1, POLR2F, and GTF2H3) that were relatively up-regulated in HL-60[R] cells and repressed in cells with ATRA-induced differentiation were related to mismatch repair in eukaryotes, DNA double-strand break repair, and nucleotide excision repair pathways. Our results suggest that transcriptional time series profiles and a functional pathway analysis of drug resistance and ATRA-induced cell differentiation will be useful for identifying promyelocytic leukemia patients who are eligible for new therapeutic strategies.

selected for 1 month, 3 months, and 6 months, and analyzed their transcriptional time series profiles, which would be undertaken to comprehensively identify defects acquired in ATRA resistance and the ATRA-induced cell differentiation response by identifying pathways that are commonly deregulated during drug resistance.  Table 1). ATRA had significant effects on HL-60 cells, and these changes in cell growth also determined the percentage of cells that were arrested at the G1/G0 phase of the cell cycle ( Fig. 1F and G). However, HL-60[R] cells were not arrested at G1/G0 after 1.0 mM ATRA treatment (Fig. 1H) because of their resistance.

Results
CGH studies of chromosomal gains, losses, and amplifications in HL-60[R] cells. The CGH results of MDR HL-60[R] and HL-60 cells, including 1 example of a CGH ratio profile, are shown in Fig. 2; genetic imbalances were present in HL-60 and MDR HL-60[R] cells. We identified 4 genetic changes in the MDR HL-60[R] cell line, including 1 loss (Yq12) and 3 gains (4q34, 9q12, and 19q13), compared with in its parental HL-60 cell line. We also confirmed that ZNF230, NUP62, E2-EPF, and CNOT3, which are located in 19q13, were up-regulated in MDR HL-60[R] cells compared to in their parental HL-60 cells (Table S1).
Hierarchical cluster analysis of differentially expressed genes regulated by ATRA treatment in resistant HL-60[R] and parental HL-60 cells. Using a hierarchical clustering analysis, we compared the similarities in the expression patterns of the 210 differentially expressed genes, which were differentially expressed in more than 3 samples at 7 time points (HL-60 cells treated with 1.0 mM ATRA for 2 days, 4 days, or 6 days and HL-60[R] cells maintained with 1.0 mM ATRA for 1 month, 3 months, or 6 months, compared to their parental HL-60 cells). A more than 2.0fold change in the transcription level was used as the cut-off value for identifying the differentially expressed genes. The dendrogram, part of which is shown in Figure 3, demonstrates the relationships between genes, as calculated by the clustering algorithm. A gradual change over the 7 time points was observed after treatment with ATRA. The genes could be divided into 4 clusters. Cluster A included those genes for which mRNA levels peaked at day 2 of treatment ( Fig. 3, Cluster A), including 10 ribosomal protein genes (RPL23A, RPL38, RPS28, RPL3, RPLP1, RPL23, RPL11, RPL19, RPL13A, and RPS24) and 2 eukaryotic translation elongation factors (EEF1G and EEF2) that were up-regulated by ATRA treatment (80.0% of genes in Cluster A). These genes were involved in protein synthesis and metabolism. Cluster B included genes that were expressed maximally at day 4 or day 6 ( Fig. 3, Cluster B); these were biomarkers of differentiation, such as known ATRA differentiation molecules (JWA, ARPC3, TMSB10, TXNIP, and DEFA1) and CD antigens (CD4 and CD63) 4 . Cluster C included genes for which mRNA levels increased after ATRA treatment at 1 month, 3 months, or 6 months (OARS, TNRC1, HNRPA1, SNRPB, PAICS, MLC1SA, HNRPK, E2-EPF, CTSC, NUP62, POH1, COX4I1, ZNF230, SF3B2, ERP70, CNOT3, and MAP2K4). Cluster D included genes for which mRNA levels increased gradually over the treatment time course and were expressed maximally at 6 months (PLA2G1B,  ARHGD1B, ACAT1, NCL, UOCRH, PPIB, GPI, DACH, PTDSS1,  PSMD8, ODC1, TAGLN2, HRMT1L2, and PKBP1A) (Fig. 3, Cluster D). These genes are related to DNA repair, cell survival, metabolism, and drug resistance and were over-expressed in MDR HL-60[R] cells.
Gene expression profiles in MDR HL-60[R] and HL-60 cells with ATRA-induced differentiation. The gene expression profiles of MDR HL-60R cells were compared to those of parental HL-60 cells. A more than 2.0-fold change in the transcription level was used as the cut-off value to identify the differentially expressed genes (p # 0.05). We found that compared to their parental HL-60 cells, 104 genes were relatively up-regulated in ATRA-resistant HL-60[R] cells and repressed in cells with ATRA-induced differentiation; 69 genes were relatively repressed in HL-60[R] and up-regulated in differentiated HL-60 cells; 6 genes (IFNGR1, TGIF, SPAG9, CSF2RB, LRRFIP1, and SLC21A3) were up-regulated in both resistant and differentiated cells; and 20 genes were repressed in both resistant and differentiated cells (Table S1). The over-expressed genes in MDR HL-60[R] are related to DNA repair, stress response, drug resistance, the ubiquitin-proteasome pathway, and protein synthesis and metabolism, including anti-oxidation, oxidative phosphorylation, and the mitochondrial pathway.
Functional networks and pathways of ATRA-induced drug resistance were analyzed by the IPA. The genetic networks and cellular pathways were derived using the IPA program by analyzing 104 genes that were up-regulated in ATRA-resistant HL-60[R] cells. A more comprehensive network and pathway analysis of all deregulated genes revealed their association with 3 important network functions and 5 critical canonical pathways, all of which are relevant to the development of ATRA-resistant cancer. The differently expressed genes constituted about half the total molecules involved and the network-associated cellular functions and include those related to RNA post-transcriptional modification; DNA replication, recombination, and repair; cell death and survival; lipid metabolism; molecular transport; developmental disorders; hereditary disorders; and metabolic disease in resistant HL-60[R] cells (Table 2). These genes are expected to be affected in HL-60[R] cells. They belong to 5 canonical signaling pathways that are commonly deregulated in ATRA resistance (Table 3). Although only 10 (FEN1, RFC5, EXO1, XRCC5, PARP1, NME1, SET, TSTA3, POLR2F, and GTF2H3) were up-regulated in HL-60[R] cells and repressed in ATRA-induced differentiated cells with the signaling pathways (eukaryote mismatch repair, DNA double-strand break repair, granzyme A signaling, GDP-L-fucose biosynthesis I, and the nucleotide excision repair pathway), each has documented functions in controlling cell growth, and the DNA repair pathway has been implicated to play roles in ATRA-resistant HL-60[R] cells.

Discussion
We challenged the MDR cancer cells in our study with retinoic acid, doxorubicin and found that they involved the activation of different mechanisms of drug metabolism and were dependent on the bioactivities of certain cancer cell lines. ATRA is known to induce the in vitro and in vivo differentiation of APL cells and favor their release from the bone marrow into the blood at the initiation of therapy. In the presence of ATRA, HL-60 human promyelocytic cells underwent myeloid differentiation. MDR HL-60[R] cells, which were highly resistant to ATRA, did not undergo cell differentiation. We demonstrated that MDR HL-60[R] cells were more than 122-, 12-, and 25fold resistant to ATRA, doxorubicin, and VP-16, respectively, compared to parental HL-60 cells (Table 1). In addition, cell cycle www.nature.com/scientificreports SCIENTIFIC REPORTS | 4 : 5577 | DOI: 10.1038/srep05577 arrest was not induced at the G1/G0 phase of the cell cycle in HL-60[R] cells but was induced in parental HL-60 cells with 1.0 mM ATRA treatment. Next, a CGH analysis of MDR HL-60[R] cells identified gains of 4q34, 9q12, and 19q13 and a loss of Yq12 as the most prominent alterations compared with in parental HL-60 cells. It is well recognized that resistant cells with genetic alterations possess a growth or survival advantage that leads to clonal expansion. To understand the molecular mechanism by which MDR is induced by ATRA in HL-60[R] cells, we used a DNA microarray to monitor the changeable gene expression profiles, from the cell differentiation induced by ATRA to the development of MDR; we also analyzed the functional networks and pathways of ATRA-induced cell differentiation and drug resistance.
Chromosomal 19q13 has been reported to contain several genes that are important in DNA repair-specifically, nucleotide excision repair and apoptosis mechanisms such as XPD, ERCC1, and RAI. The 19q13 chromosomal region is important in cancer 13 . Our microarray analysis revealed that genes in the 19q13 chromosomal region-ZNF230, NUP62, E2-EPF, and CNOT3-were constitutively   (Table S1 and Fig. 3C).
NUP62 is an essential component of the nuclear pore complex and plays a novel role in centrosome integrity. Knockdown of NUP62 induced G2/M phase arrest, mitotic cell death, and aberrant centrosome and centriole formation 14 . E2-EPF, named UBE2S (ubiquitinconjugating enzyme E2S), accepts ubiquitin from the E1 complex and catalyzes its covalent attachment to other proteins before elongating ubiquitin chains on APC/C substrates to promote mitotic exit 15 . High expression of E2-EPF was indicative of poor overall survival in a large-scale co-expression analysis of breast cancer 16 . It is interesting that CNOT3 has also been associated with transcription regulation and represents a novel component of the core self-renewal and pluripotency circuitry that is conserved in mouse and human embryonic stem cells with CNOT1 and CNOT2 17 . The chromosomal gene set enrichment analysis further confirmed that genes located at 19q13 were expressed at higher levels in uterine carcinosarcoma and contributed to its poor prognosis because of its epithelial-to-mesenchymal transition characteristics 18 . In our microarray, many genes involved in the retinoic acid signaling pathway (such as RARRES2, CRABP2, MYC, TGIF, SUPV3L1, SET, ILF2, PTPRCAP, GPRK6, and SEC4L) were up-regulated in MDR HL-60[R] cells. In previous studies, RARRES2 (retinoic acid receptor responder 2 or chemerin) was over-expressed in malignant mesothelioma cells 19 and adrenocortical tumors 20 and was a biomarker for insulin resistance in type 2 diabetes mellitus patients 21 . TGIF functions as a transcriptional co-repressor and regulates developmental signaling by retinoic acid; it may also repress other RXRdependent transcriptional responses 22 . Moreover, TGIF's levels are inversely correlated with survival in patients with acute myelogenous leukemia, and its knockdown inhibits the differentiation of myeloid cell lines and increases apoptosis 23,24 . On the other hand, retinoic acid-resistant neuroblastoma cell lines also showed over-expression of MYC with ATRA in culture medium 25 , and SET (SET nuclear oncogene), a potential marker for HNSCC that is associated with drug resistance, is up-regulated in 97% of tumor tissue samples and HNSCC cell lineages 26 .
Next, we found that the genes involved in oxidative phosphorylation and metabolism (ACAT1, ATP5G3, ARF4L, HEAB, PDHA1, and GNPI) were differentially expressed in HL-60[R] cells; other genes involved in protein synthesis metabolism, such as eukaryotic translation initiation factors (EIF2S1 and EIF3S9), transcription and elongation factors (GTF3A, CNOT3, RUNX1, and TCEA1), and splicing factors (SFRS3 and SFRS10), were also over-expressed in MDR HL-60[R] cells. ACAT1 expression may be a prognostic marker in prostate cancer: it was expressed at a significantly higher level in cancerous cores than in adjacent benign cores and was specifically effective at differentiating between indolent and aggressive forms of cancer 27 . Moreover, targeting ACAT1 with avasimibe (ACAT inhibitor) could be an efficient treatment for glioblastoma because it can inhibit ACAT1 expression and induce cell apoptosis 28 . The expression of ARF4L (adenosine diphosphate-ribosylation factor 4-like), a glioma-associated antigen, is controlled by the activated Akt/mTOR pathway, which is a downstream effect of the loss of PTEN 29 . HEAB (cleavage and polyadenylation factor I subunit 1) is a protein-coding gene that contains an adenosine triphosphate/guanosine triphosphate-binding motif that is homologous to the adenosine triphosphate-binding transporter superfamily or guanosine triphosphatebinding proteins 30 . RUNX1 may play a critical role in chemotherapy response in acute megakaryocytic leukemia by regulating the phosphoinositide 3-kinase/Akt pathway 31 , and ectopic RUNX reduces intracellular long-chain ceramides in NIH3T3 fibroblasts and elevated extracellular sphingosine 1 phosphate. RUNX expression also opposed the activation of c-Jun-NH(2)-kinase and p38 (mitogenactivated protein kinase) and suppressed the onset of apoptosis in response to exogenous tumor necrosis factor a 32 . SFRS3 (SRp20) is a serine-and arginine-rich splicing factor and proto-oncogene that is critical for cell proliferation and tumor induction and maintenance. Increased expression of SFRS3 in rodent fibroblasts promoted immortal cell growth and transformation 33 . However, depletion of SFRS10 resulted in apoptosis of the neural progenitor cells as well as disorganization of the cortical plate 34 .
It is interesting that many of the genes related to DNA repair and cell survival (NME1, DDX1, YWHAQ, PSMC1, PSMB6, E2-EPF, FEN1, RFC5, EXO1, UNG, XRCC5, ADPRT, POLR2F, Rpo1-2, and GTF2H3) are over-expressed in MDR HL-60[R] cells. Of note, POLR2F, XRCC5, and NME1 were reported to be repressed in retinoid-induced cell differentiation 4 . Most of these genes are related to mismatch repair in eukaryotes (FEN1, RFC5, and EXO1), DNA double-strand break repair (XRCC5 and PARP1), and the nucleotide excision repair pathway (POLR2F and GTF2H3) on IPA analysis, and they are commonly up-regulated in ATRA-resistant HL-60[R] cells. Among these DNA repair genes, FEN1 (Flap endonuclease 1) is www.nature.com/scientificreports SCIENTIFIC REPORTS | 4 : 5577 | DOI: 10.1038/srep05577 a highly conserved structure-specific nuclease that catalyzes a specific incision to remove 59 flaps in double-stranded DNA substrates; it plays an essential role in key cellular processes, such as DNA replication, repair, and mutation 35 . ADPRT (PARP1) is a critical DNA repair enzyme that is involved in DNA single-strand break repair via the base excision repair pathway. PARP inhibitors have been shown to sensitize tumors to DNA-damaging agents and selectively kill BRCA-deficient cancers 36,37 . XRCC5 (X-Ray repair complementing defective repair in Chinese hamster cells 5) was over-expressed in cisplatin-resistant ovarian cancer cell lines 38 ; it also affected chemosensitivity 39 and was associated with the MDR phenotype 40 . An immunohistochemical analysis verified significant co-expression of MDR1 and NME1 in human epithelial ovarian carcinoma 41 . Increased NME1 mRNA levels were associated with resistance to initial chemotherapy in acute monocytic leukemia 42 and helped cells become resistant to oxidative stress 43 . On the other hand, higher levels of UNG (uracil-DNA glycosylase) were associated with pemetrexed resistance, and induction of UNG protein confirmed that up-regulation of the base excision repair enzyme is a feature of acquired pemetrexed resistance 44 . DDX1 is a member of the DEAD box protein family, which has RNase activity, plays an RNA clearance role at DNA double-strand break sites, and facilitates the templateguided repair of transcriptionally active regions of the genome 45 .
We also demonstrated that about 69 genes, including MYO5A, PAK1, FGR, PET112, GLS, ALB, MYD88, APOC2, and KLKB1, were up-regulated in HL-60 cells treated with ATRA and were repressed in HL-60[R] cells. The genes in cells with ATRA-induced differentiation were involved in cellular compromise, developmental and hereditary disorders, small molecule biochemistry, cell death and survival, cellular movement, hematological system development, and functional networks (Table 2). Next, we found that ATRA induced cell differentiation via pathways, including FcY receptormediated phagocytosis in macrophages and monocytes (MYO5A, PAK1, and FGR), L-glutamine biosynthesis II (PET112), glutamine degradation I (GLS), interleukin-12 signaling and production in macrophages (ALB, MYD88, and APOC2), and acute-phase response signaling pathways (KLKB1, ALB, and MYD88) ( Table 3). These functional networks and pathways will help us understand the dis-ease response to retinoic acid's biological effects in promyelocytic leukemia chemotherapy.
In summary, in this study, we characterized the genes that are involved in DNA repair and cell survival, oxidative phosphorylation and metabolism, and the retinoic acid signaling pathways in ATRAresistant HL-60[R] cells; we further revealed novel coordinated changes that occurred in resistant cells that allowed them to survive the cell differentiation and apoptosis elicited by ATRA chemotherapy. Thus, it is possible that blocking DNA repair and cell survival signaling pathways will not only enhance ATRA chemotherapy and improve the outcomes of patients with APL but also reduce the risk of second primary tumors.
Cytotoxicity assay. We performed an MTT (3-(4,5 dimethylthiazol-2-yl)-2,5 tetrazolium bromide) assay, which is based on the enzymatic (mitochondrial dehydrogenase) reduction of the tetrazolium salt, MTT, to a colored formazan product by viable cells 46 . Cells were plated in 96-well microassay culture plates at a cell density of 10 4 cells/well and grown overnight at 37uC in a 5% CO 2 incubator. Test compounds were then added to the wells to achieve a final concentration of 10 26 to 10 24 M. Control wells were prepared by adding 100 ml of culture medium with no cells. The plates were incubated at 37uC in a 5% CO 2 incubator for 72 h, and 20 ml of the stock MTT dye solution (5 mg/ml) was added to each well. After 4 h of incubation, 100 ml of DMSO was added to solubilize the MTT formazan. The optical density of each well was then measured with a microplate spectrophotometer at a wavelength of 570 nm. The IC 50 was determined from the plots of the percentage viability vs. the dose of compound added.
Cell cycle analysis. To determine the cell cycle distribution, we plated 0.2 3 10 6 HL-60 and HL-60[R] cells in 60-mm dishes and treated them separately with 1.0 mM ATRA for 2 days. Cells were then collected and fixed in 95% ethanol, washed in 1% bovine serum albumin and phosphate-buffered saline, resuspended in 1.0 g/ml RNase and 50 mg/ml propidium iodide, incubated for 30 min in the dark at 37uC, and analyzed by flow cytometry using FACSCalibur. The data were analyzed using the ModFit DNA analysis program.
Comparative genomic hybridization and digital image analysis. Comparative genomic hybridization (CGH) was performed essentially as described previously 47 CGH was performed using a digital image analysis system that contained a Zeiss Axioplan 2 microscope equipped with a Sensys cooled-charged device camera (Photometrics, Ltd., Tucson, AZ, USA), as previously described 48 . High copy number amplification of a chromosomal region was defined as a tumor5reference or HL-60[R]5HL-60 ratio of 1550.
RNA extraction and gene expression profiling. RNA was isolated with Trizol LS (Invitrogen, Carlsbad, CA, USA) and purified using the RNeasy Mini Kit (Qiagen, Valencia, CA, USA). RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, DE, USA). Microarray experiments were carried out using whole human genome oligonucleotide arrays with 20 mg of total RNA starting material, according to the manufacturer's protocol (National Cancer Institute [NCI] array). Total RNA from each sample was synthesized into double-stranded cDNA with SIII reverse transcriptase (Invitrogen) using an oligonucleotide d(T) primer. The double-stranded cDNA from HL-60 cells was labeled with Cy3 monofunctional reactive dye (Amersham Biosciences), and that from HL-60 [R] cells or ATRA-treated HL-60 cells was labeled with Cy5 monofunctional reactive dye (Amersham Biosciences). The probe was hybridized to an NCI oligonucleotide array (Hs-Operon V3-v1p24.gal) containing 36 K human transcripts (NCI Microarray Facility, Advanced Technology Center, Gaithersburg, MD) overnight at 42uC. For each treatment, the arrays were also queried with probes produced via reverse labeling, and the data were consistent with those obtained with initial standard labeling. Microarray slides were scanned with a GenePix 4000B microarray scanner (Axon Instruments, Union City, CA).
Microarray data and pathway analysis. The microarray images were analyzed with GenePix 5.1 software, and the subsequent gene lists and associated expression values were loaded into mAdb (NCI Microarray Facility). Fluorescence ratios were normalized for each array using a single scaling factor so that the median fluorescence ratio of well-measured spots on each array was 1.0. After flagging the bad spots, the mean log 2 -transformed ratio of resistant versus sensitive cells was calculated from triplicate experiments. The mean data were calculated using the antilog as the ratio of gene expression measures of resistant cells to parental cells and analyzed using Partek Genomics Solution software, as reported in our previous studies 49 . To determine the specific pathways on the basis of changes in gene expression, we used the Ingenuity Pathway Analysis (IPA) (http://www.ingenuity.com) commercial gene pathway analysis web tool.
Statistical analysis. All statistical values are presented as means 6 S.D. Data were analyzed using Student's t-test. The results were considered significant at p , 0.05.