Identification of potential glucocorticoid receptor therapeutic targets in multiple myeloma.

Glucocorticoids (GC) are a cornerstone of combination therapies for multiple myeloma. However, patients ultimately develop resistance to GCs frequently based on decreased glucocorticoid receptor (GR) expression. An understanding of the direct targets of GC actions, which induce cell death, is expected to culminate in potential therapeutic strategies for inducing cell death by regulating downstream targets in the absence of a functional GR. The specific goal of our research is to identify primary GR targets that contribute to GC-induced cell death, with the ultimate goal of developing novel therapeutics around these targets that can be used to overcome resistance to GCs in the absence of GR. Using the MM.1S glucocorticoid-sensitive human myeloma cell line, we began with the broad platform of gene expression profiling to identify glucocorticoid-regulated genes further refined by combination treatment with phosphatidylinositol-3'-kinase inhibition (PI3Ki). To further refine the search to distinguish direct and indirect targets of GR that respond to the combination GC and PI3Ki treatment of MM.1S cells, we integrated 1) gene expression profiles of combination GC treatment with PI3Ki, which induces synergistic cell death; 2) negative correlation between genes inhibited by combination treatment in MM.1S cells and genes over-expressed in myeloma patients to establish clinical relevance and 3) GR chromatin immunoprecipitation with massively parallel sequencing (ChIP-Seq) in myeloma cells to identify global chromatin binding for the glucocorticoid receptor (GR). Using established bioinformatics platforms, we have integrated these data sets to identify a subset of candidate genes that may form the basis for a comprehensive picture of glucocorticoid actions in multiple myeloma. As a proof of principle, we have verified two targets, namely RRM2 and BCL2L1, as primary functional targets of GR involved in GC-induced cell death.

Introduction

For over 40 years, GCs have been used to successfully treat myeloma patients, and they have continued to be included as part of combination therapies for the last 20 years, highlighting the important contribution of this drug to the treatment of myeloma [Cavo et al., 2011; Rajkumar et al., 2002]. In addition, expression of the glucocorticoid receptor (GR) in myeloma patients correlates with better overall patient survival, underscoring the importance of this signaling pathway in myeloma clinical outcome [Heuck et al., 2012]. Despite the long-term use of this drug, the mechanism of action is largely unknown. With the growing trend of combining therapeutics that target different pathways, it becomes imperative to understand the primary downstream targets for GCs to enable efficacious use of this drug. In addition, understanding the GC targets that are required for cell death will provide the basis for the development of novel therapeutics that induce cell death in the face of GC resistance caused by loss of the GR.

GCs mediate their biological effects through the GR. We and others [Moalli et al., 1992] have demonstrated that the intact receptor is required for apoptosis, as hematologic cell lines with absent or mutant receptors evade the cytotoxic effects of GCs [Greenstein et al., 2002; Kofler, 2000]. In canonical GC signaling, ligand binding induces a conformational change in GR, which releases the receptor from its chaperones, and allows translocation to the nucleus, where it then binds to a consensus DNA sequence termed the glucocorticoid response element (GRE). Interaction of GR with the GRE stimulates binding of nuclear receptor co-regulators, which modulates transcription of target genes. In addition to induction of gene expression, GR also acts to repress transcription. The trans-repressive functions of the GR have been linked to the anti-inflammatory and pro-apoptotic actions of GCs [Zanchi et al., 2010]. Proposed mechanisms for repression include GR tethering to and inhibiting growth induction transcription factors through protein-protein interactions or by interactions with growth induction transcription factors at composite GREs [Yamamoto et al., 1998]. However, despite extensive gene expression array analyses, specific genes that mediate cell death have not been well defined, perhaps in part because recent data support GCs acting through the GR to impact a complex network of activators [Beato et al., 1995; Hollenberg et al., 1985] [Martin, 2003; Wilson et al., 2013].

Initial identification of the glucocorticoid receptor as a hormone-activated transcription factor [Beato et al., 1995; Hollenberg et al., 1985] was the basis of the hypothesis that glucocorticoid-induced cell death involved the activation of specific “death” genes. Several laboratories, including ours, have conducted GC-regulated gene expression analyses in a variety of cell lines of hematologic origin. Over 900 different genes have been reported to be GC-regulated, but only about 70 have appeared in more than one publication [Schmidt et al., 2004]. Due to the difficulties in identifying a single primary death-inducing gene, it has been suggested that multiple cell type-dependent mechanisms may exist. Rather than a single conserved canonical pathway causing GC-induced cell death, a network of GC-induced gene action interacting with other signaling pathways may be triggering cell death [Schmidt et al., 2004]. As an example, we have noted that GILZ (glucocorticoid induced leucine zipper), a protein involved in GC-mediated cell death, is synergistically induced by a combination of GC and PI3 kinase pathway inhibition [Grugan et al., 2008]. More recently, systems biology approaches indicate that cellular signaling pathways, such as glucocorticoid signaling, are complex, integrated pathways with primary nodes of interaction leading to complex downstream signaling events [Martin, 2003; Wilson et al., 2013].

To add further complexity, gene expression array data captures regulation of mRNAs that may be primary targets of GR but also may be targets of effectors downstream of GCs. The location and sequence of specific GR binding regions (GBR) in the genome have been published for other cell types [John et al., 2011; So et al., 2007; Yu et al., 2010a]. To distinguish primary GC targets that are directly regulated by GR from the indirect targets that are impacted downstream in our multiple myeloma model, we employed chromatin immunoprecipitation (ChIP) combined with massively parallel sequencing (Seq). The location of GBRs, by proximity to gene coding regions, will identify potential GC-regulated genes. Combining the location data of GC-regulated genes and mRNA expression data sets will identify primary target candidates of GC activity [John et al., 2011; So et al., 2007; Yu et al., 2010a]. The ChIP-Seq studies reported here identify potential direct gene targets of GR that could be targeted to overcome GC-resistance.

ChIP-Seq studies of GR have recently been published for lung cancer cells [So et al., 2007], mouse breast cancer cells [John et al., 2011], mouse pituitary cells [John et al., 2011], and adipocytes [Yu et al., 2010a], with some startling results. GR ChIP-Seq data indicates that in each cell type examined, approximately 90% of the GBRs were located more than 5 KB from the transcription start site (TSS) and not proximal to the TSS, as had been assumed in the classical model of hormone action. These data imply that the GBR is not acting as a promoter, as was initially thought, but rather in the capacity of an enhancer which regulates the activities of a distant promoter through looping [Biddie et al., 2010] to interact with basal transcription machinery. When the GRE mapping data is combined with DNase I hypersensitive sites (DHS) mapping of active chromatin conducive to initiation of transcription, the vast majority of GBRs occur at a DHS, as expected [John et al., 2011]. Interestingly, DHS mapping and the associated GBRs are radically different between cell types, indicating that active GBRs are cell type-specific, which is dictated by the pre-existing chromatin landscape [John et al., 2011]. These data indicate that ChIP-Seq is critical in the identification of active GBRs in myeloma cells.

In this report we have enriched the population of GR-regulated genes that contribute to cell death by building on our previous observations that combine a low dose of PI3 kinase inhibitor with glucocorticoid treatment to enhance cell death in myeloma cell lines [Grugan et al., 2008]. To assess clinical relevance, we interrogated gene expression array data from myeloma patients to identify those genes that are over-expressed in myeloma and down-regulated by our combined drug treatment. To further identify those GC-regulated genes that may be primary targets for GR binding in a myeloma setting, we performed ChIP-Seq using a myeloma cell line. We used bioinformatics platforms to combine these data sets to identify candidate genes involved in GC-induced cell death. Gene set enrichment analyses identified signaling pathway enrichment for GC-regulated genes as candidates for GC-induced cell death. Ultimately, selected candidate genes were examined for participation in the cell death mechanism through genetic manipulation and assessing the impact on glucocorticoid sensitivity as proof of principle. In the future, this validated gene set can provide the framework to identify and develop critical components of GC-induced cell death as therapeutic targets.

Materials and Methods

Cell culture: MM.1S and MM.1RL cell lines were previously established from a myeloma patient who had been treated with glucocorticoids and become resistant to that treatment. MM.1S cells were sub-cloned from the initial patient isolate through dilution cloning and are sensitive to glucocorticoids, as measured by the induction of cell death [Goldman-Leikin et al., 1980]. MM.1RL cells were isolated through step-wise incubation with increasing concentrations of glucocorticoids and are resistant to the glucocorticoid-induced cell death [Moalli et al., 1992]. MM.1RL cells do not express the glucocorticoid receptor, while MM.1S cells express the glucocorticoid receptor at levels comparable to normal B cells. Cells are cultured as previously described [Tessel et al., 2011] in RPMI-1640 supplemented with 10% fetal bovine serum, 2 mm glutamine, 100 U/mL penicillin, 100 mg/mL streptomycin, 2.5 μg/mL fungizone and 0.5 μg/mL plasmocin (InvivoGen) and maintained in a 370 C incubator with 5% CO2.

Cell death assays: Annexin/DAPI staining was performed according to manufacturer directions (BD Biosciences,). Stained cells were analyzed with the LSR Fortessa flow cytometer (BD Biosciences). A total of 20,000 events were collected for experiments. The data were analyzed using FlowJo software (Tree Star).

Immunoblotting: Immunoblotting was performed as described previously [Tessel et al., 2011]. The antibodies used were BCL-XL (Abcam, ab2568), RRM2 (Santa Cruz, E-16 sc-10846), β-actin (Cell Signaling, 8H10D10) and α-tubulin (Cell Signaling, DM1A #3873S).

Quantitative real-time PCR (qRT-PCR): Reverse transcription polymerase chain reaction (RT-PCR) was performed as previously described [Rovedo et al., 2011] with the following modifications. Total RNA was extracted from cells using RNeasy Kit (Qiagen) according to manufacturer instructions. Complementary DNA was synthesized from 1μg of RNA using SuperScript Vilo cDNA Synthesis Kit (Invitrogen). qRT-PCR was performed using Taqman gene expression assays (BCL2L1, RRM2, and GILZ) according to manufacturer protocols (Applied Biosystems, Grand Island, NY) in an ABI 7900 Fast instrument. RNA expression was normalized to three independent control genes (YWHAZ, EIF4A2 and RPL13A; Primer Design, Southampton, UK) using the geNorm normalization algorithm (http://medgen.ugent.be/~jvdesomp/genorm/).Eachassaywasrunintriplicatefrombiologictriplicates.

Gene expression microarray: Total RNA was extracted using Trizol (Life Technologies) and purified with the RNeasy Mini Kit (Qiagen) following manufacturer instructions. RNA was reverse transcribed and the microarray hybridization was performed using the Illumina Gene Expression Sentrix BeadChip HumanHT-12_V4 (Illumina) at the Genomics Core, Northwestern University. Microarray scanned images were imported to Illumina® GenomeStudio for data quality control and the transcriptome profile data was quantile-normalized by the Bioconductor lumi package [Du et al., 2008]. Significant gene expression changes were determined using the t-test (p<0.05) and fold change exceeding 1.25. For hierarchical clustering, gene expression was transformed to z-scores, with respect to the control values, using the R statistical system. The data set has been deposited in GEO with the following link: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=upyjqwwmdlurdqh&acc=GSE64704

We performed Gene Set Enrichment Analysis (GSEA) of primary transcriptome responses, as well as of gene signatures derived via integrative analysis, using the software resources available at the Molecular Signature Database (MSigDB, http://www.broadinstitute.org/gsea/msigdb/index.jsp). A Multiple Myeloma patient gene signature was generated from a publicly-available gene set, GSE6691 [Gutierrez et al., 2007]. In that study, the myeloma patients are newly-diagnosed and untreated. Significant gene expression changes were determined using the t-test (p<0.05) and fold change exceeding 1.25. Comparison with transcriptomic response of MM cells to treatment with GC, PI3K inhibitor, and combination treatment was performed using the R statistical system, as previously described by our group [Bliss-Moreau et al., 2015]; similar to the analysis strategy used in our prior publication, we did not apply multiple hypothesis testing corrections for the primary gene expression analysis, thereby enabling the integrative analysis to yield a rich set of biologically-interesting targets.

ChIP-Seq: Chromatin immunoprecipitation was performed as described previously [Yu et al., 2010b], with the following modifications. MM.1S and MM.1RL cells were treated with or without 1 μM dexamethasone (Dex) for 2 h prior to cell harvest. Protein crosslinking to chromatin was achieved through incubation with 1% formaldehyde and chromatin isolated as previously described. Chromatin was sheared by ten 12 second cycles of sonication and lysate cleared by incubation with protein-A linked agarose beads. Antibody against GR (Santa Cruz SC-8992X) or non-specific IgG was added at a concentration of 3 to 5 μg per sample and incubated overnight at 4o C with rotation. Antibody-bound chromatin was concentrated by incubation with Protein A-linked agarose beads. Crosslinking was reversed by incubation with RNase A and protein digested with Proteinase K. DNA was purified by ethanol precipitation and used for library construction.

The libraries were sequenced on a HiSeq 2000 Sequencing System with a single-end 40 cycle run, yielding 17-31 million reads per sample. Sequenced reads were mapped to the human genome (UCSC build hg19/NCBI build 37) using the Burrows-Wheeler Alignment tool (BWA) [Li and Durbin, 2009]. High-resolution genome-wide maps were derived and visualized in the UCSC Genome Browser (http://genome.ucsc.edu/). The model-based analysis of ChIP-Seq data macs2 [Zhang et al., 2008] peak calling package was used to identify enriched regions, using default parameters. Motifs were inferred using the homer [Heinz et al., 2010] and the meme-chip [Machanick and Bailey, 2011] tools. Genes targets of GR were determined by computing the overlap between the ChIP-Seq peaks and the genes using bedtools [Quinlan and Hall, 2010] and the RefSeq gene definition considering 100,000 base-pair windows around each gene. The data set has been deposited in GEO with the following link: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=gbcvksqedhudrux&acc=GSE63209

Parallel sequencing: Chromatin immunoprecipitated DNA was prepared for sequencing following a modified version of the Illumina Genomic DNA protocol and as described previously [Yu et al., 2010b](above). In brief, sheared DNA fragments undergo end repair by T4 DNA polymerase to fill in 5′ overhangs, Klenow polymerase to remove 3′ overhangs, and T4PNK to phosphorylate the 5′-OH. For Illumina adapter ligation, a single adenine nucleotide overhang must be added to the polished DNA. A 1:30 dilution of the Adaptor Oligo Mix was used in the ligation step. The resulting adapted constructs are amplified by 18 cycles of PCR enriching the DNA fragments with adapters flanking each end. The amplified product was quantified using a spectrophotometer at 260 nm (Nanodrop), and then a portion of final sample is size separated on a 6% TBE (Invitrogen) to verify correct target range amplification. The validated DNA library is submitted at a 10 nmol/L concentration for sequencing. All protocols for Illumina library sequence preparation, sequencing, and quality control are provided by Illumina.

Statistics: Significance was calculated by Student’s T-tests. Paired T-tests were used for determining difference in gene expression (Figure 4), and unpaired T-tests were used for calculating difference in cell viability (Figure 5). Significance: * = p≤0.05, ** = p≤0.01, *** = p≤0.0001. Each experiment consists of three independent replicates and is presented as the mean ± the standard deviation.

Figure 4

Confirmation of gene and protein regulation by GC/PI3K inhibitor treatment.

(A) GSEA was performed on the gene signature of all Dex/LY294002 down-regulated genes and the narrower signature of down-regulated genes with potential GR binding sites. Gene sets that were significantly enriched in both signatures are shown. Direct GR binding and regulation of genes within these gene sets may affect the down-regulation of additional genes in the set. Black arrows indicate gene sets that include BCL2L1 and blue arrows indicate gene sets containing RRM2. (B and C) MM.1S and MM.1RL were treated with 1 μM Dex, 25 μM LY294002, or in combination for 6 h and BCL2L1 (B) and RRM2 (C) expression was measured by qRT-PCR. (D) BCL-XL and RRM2 protein expression was quantified by immunoblot and normalized to β-actin and α-tubulin, respectively.

Figure 5

Knock-down of BCL-XL and RRM2 sensitizes cells to GC-induced apoptosis and is required for cell survival.

(A and D) Confirmation of BCL-XL (A) and RRM2 (D) knock-down (KD) by siRNA at 24 and 96 h post-transfection. (B and E) Cell death measured by Annexin V staining of KD and nonspecific control (NSC) MM.1S and MM.1RL cells 96 h post-transfection. (C and F) 24 h post-transfection, KD and NSC MM.1S cells were treated with increasing concentrations of Dex for 72 h. Cell death was measured by Annexin V staining.

Results

Identification of genes modulated by GC and further modified by PI3K inhibitor treatment

We had previously observed that the combination treatment of glucocorticoids with a sub-lethal concentration of the PI3K inhibitor LY294002 causes synergistic cell death in the MM.1S human myeloma cell line [Grugan et al., 2008]. In addition, we demonstrated that this combination caused a synergistic increase in gene expression of GILZ (glucocorticoid induced leucine zipper), which contributes to glucocorticoid-induced cell death. Based on these observations, we hypothesized that the subset of glucocorticoid-modulated genes that are further modulated by combination with the PI3K inhibitor would enrich for those genes involved in GC-induced cell death. To identify the cell death gene set, we performed a gene expression array with total RNA isolated from MM.1S cells treated either with glucocorticoids, PI3K inhibitor, a combination of GC + PI3K inhibitor or vehicle control (Figure 1). The concentration of PI3K inhibitor is based on previously-published dose-response in MM.1S cells [Grugan et al., 2008] where induction of GILZ expression is observed, but there is no increased cell death. Similarly, there is no observed cell death in the GC-resistant myeloma cell line MM.1RL. We chose to examine gene expression at a time point of 3 h, where GILZ expression is elevated, but well before initiation of cell death, in order to identify early events in cell death initiation. As shown in the heat map of gene expression (Figure 1 A), we observed a greater number of down-regulated genes in all treatments as opposed to up-regulated genes similar to published data of GC regulation of gene expression at early time points in myeloma cell lines [Chauhan et al., 2002b]. Using the relaxed cut-off of 1.25-fold change, we found 834 genes down-regulated by Dex, 2235 genes down-regulated by inhibition of PI3K, and 3180 genes down-regulated by the combination drug treatment. In contrast, we observed 388 genes up-regulated by Dex, 832 genes up-regulated by PI3K inhibition, and 1407 genes up-regulated by the combination drug treatment (Figure 1 B). While genes that are up-regulated by combination treatment may be involved in the induction of cell death, we chose to focus on those genes that are down-regulated by combination treatment, due to current therapeutic approaches that favor inhibition of over-expressed proteins through application of small molecule inhibitors. Interestingly, many of the combined treatment down-regulated genes – those genes down-regulated by GC treatment and further down-regulated by combined GC and PI3K inhibition – are known to be involved in cell cycle progression, transcription, translation, metabolism and repression of apoptosis, as indicated by gene set enrichment analysis [Subramanian et al., 2005] (GSEA, q<0.05) (Figure 2). These are all pathways that are under active investigation for the development of pharmaceutical inhibitors.

Figure 1

Identification of genes regulated by combined GC and PI3K inhibition in MM.1S cells.

(A) Heat map of gene expression following 3-h treatment of MM.1S cells with either vehicle control, 1 μM Dex, 25 μM LY294002, or a combination of Dex and LY294002. (B) Venn Diagram illustrating overlap of genes up- and down-regulated by GCs and PI3K inhibition alone and in combination.

Figure 2

Identification of potential death-associated pathways.

GSEA was performed on a GC/PI3K inhibitor down-regulated gene set composed of genes down-regulated by GC treatment and further down-regulated by combined GC and PI3K inhibitor treatment. The repression of these pathways/processes may be important for GC-induced myeloma cell death.

Clinical relevance of drug-induced down-regulated genes is supported by myeloma patient gene expression pattern

To assess the identified gene sets for potential clinical relevance, we used a publicly-available gene expression data set (GSE6691), and we compared the gene expression from the normal B lymphocytes (8 samples) and normal plasma cells (5 samples) to gene expression in myeloma patients (12 samples) [Gutierrez et al., 2007]. The myeloma patients were newly-diagnosed and untreated. Here, we found an almost equivalent number of deregulated genes expressed in MM patients; namely, 2850 genes were down-regulated and 2895 genes were up-regulated compared to control cells. Not surprisingly, when GSEA was performed on these gene sets, we found enrichment for B cell and immune signaling pathways, as well as enrichment of apoptosis, cell death, cancer biology and cell metabolism pathways (GSEA, q<0.05) (Figure 3 A and B, and see Supplementary Material Figure 1). Using bioinformatics platforms, we next integrated the MM patient gene expression profile (GEP) with our myeloma cell line combination treatment GEP. We identified a gene set that is up-regulated in myeloma patients and down-regulated by our combination drug treatment in MM.1S myeloma cell line (Figure 3 C), with the hypothesis that this gene set would be amenable to therapeutic targeting in myeloma patients through small molecule inhibitors. We observed 266 genes to be down-regulated by the combination treatment. Of those, 97 genes were up-regulated in the myeloma patient GEP data set. In contrast, of 125 genes up-regulated by combination treatment in myeloma cells, 25 of those genes were down-regulated in the MM GEP.

Figure 3

Clinical relevance of drug-induced combination treatment down-regulated genes is supported by multiple myeloma patient gene expression pattern.

(A) Gene Expression Profile of myeloma patients in comparison to normal B cells and plasma cells. (B) GSEA of GO Biological Processes gene sets was performed on the myeloma patient gene signature determined in (A). Processes involving metabolism, cell death, and immune response were enriched. Top 30 most significant gene sets shown. (C) Number of genes with correlation between MM patients and combination treatment regulation with GC/PI3K inhibitor combination treatment. Importantly, 97 genes were found to be up-regulated in MM patients, but down-regulated following Dex/LY294002 treatment.

Identification of direct GR gene targets in the MM.1S genome

To refine the list of potential therapeutic targets, we sought to identify combination treatment-regulated genes potentially under the direct transcriptional control of GR. These candidate direct targets are likely to be nodes of control that impact a network of downstream targets to induce cell death, making them therapeutically-relevant targets. To identify the potential direct targets of GR transcriptional regulation, ChIP-Seq was performed. We treated MM.1S and MM.1RL myeloma cell lines with vehicle control or 1μM Dex for 2 h and isolated GR-bound chromatin associated DNA as described in Materials and Methods. The MM.1RL cells do not express GR and were included as a negative control to confirm antibody specificity. In addition, the no hormone group of MM.1S cells indicates the baseline of GR-chromatin interaction in the absence of hormone. The isolated DNA was used to generate tagged libraries of fragments for next generation sequencing.

We identified 8689 high-confidence GR binding sites in the Dex-treated MM.1S group. In comparison, the MM.1S control group had 177 peaks, the MM.1RL control group had 26 peaks, and the Dex-treated MM.1RL group had 9 peaks, indicating the high level of specificity in the MM.1S Dex-treated group (see Supplementary Material Figure 2). A consensus binding site was identified which closely parallels the consensus GR binding site, with specific binding predominantly within 2500 bp of the transcription start site (1708 peaks out of a total of 8688 exhibit the canonical GRE motif, p=10-1165). Other motifs highly-enriched at GR binding sites include the ETS motifs Fli1 (p<10-507), ETS1 (p< 10-447) and ETV1 (p<10-432), and IRF motifs IRF4 (p<10-455) and IRF1 (p<10-432). To pair the GR binding sites to genes, the binding information was overlaid on the UCSC genome browser, and genes with binding were defined as containing a GR binding site at a distance up to 100kB from the transcriptional start site of each gene. This parameter was used because GR has been described to bind to enhancer regions up to 100kB from the transcription start site [John et al., 2011]. While this cut-off may not identify all GC-regulated genes with GR binding sites, it provides a starting point for analysis and identification of potential candidate genes. Using this data, we were able to select those Dex/LY294002 combination treatment-regulated genes that may be direct transcriptional targets of GR and clinically-relevant based on GEP of MM patients. We identified 57 genes down-regulated by GC/PI3K inhibitor treatment that are candidates for GR binding in myeloma (Table 1) and 22 that are up-regulated (Table 2). Therefore, we have identified the potential direct GR transcriptional targets whose repression correlates with Dex-induced myeloma cell death. Of note, motifs enriched at the GR peaks within 100kB from the genes presented in Table 1 and Table 2 include GR (p<10-17 for Table 1, p<10-22 for Table 2) and Fli1/ETS (p<10-14 for Table 1, p<10-9 for Table 2) motifs.

Table 1

Genes with GR binding sites and negative correlation between Myeloma patients and combination drug treatment groups.

SYMBOL	DEFINITION	FC MM up	FC Treated down
XBP1	X-box protein 1	8.93	0.53
RRM2	ribonucleotide reductase M2	4.47	0.50
RBM47	RNA binding motif protein 47	4.17	0.53
ARPP19	cAMP-regulated phosphoprotein, 19kDa	3.91	0.61
IFNAR2	Interferon receptor 2	3.53	0.77
ST6GALNAC4	ST6-N-acetylgalactosaminide alpha-2,6 sialytransferase	3.25	0.61
PREB	Prolactin regulatory element binding	2.97	0.42
PPIF	peptidylprolyl isomerase F	2.79	0.40
CCND2	Cyclin D2	2.70	0.61
TLE3	transducin-like enhancer of split 3	2.48	0.53
PTS	6-pyruvolyltetraheydropterin synthase	2.38	0.51
SRP19	signal recognition particle 19kDa	2.14	0.70
BCL3	BCL3	2.09	0.36
GNL3	guanine nucleotide binding protein-like 3	2.02	0.77
CAV2	caveolin 2	1.97	0.35
TFB2M	transcription factor B2, mitochondrial	1.96	0.57
TRIAP1	TP53 regulated inhibitor of apoptosis 1	1.92	0.42
CARM1	coactivator-associated arginine methyltransferase 1	1.90	0.54
MYC	Myc	1.82	0.08
CISD1	CDGSH iron sulfur domain 1	1.81	0.62
CSNK1E	casein kinase 1, epsilon	1.80	0.53
BCL2L1	BCL-XL	1.78	0.20
LEPROTL1	Leptin receptor overlapping transcript-like 1	1.75	0.27
COQ10B	coenzyme Q10 homolog B	1.72	0.67
ZNF22	zinc finger protein 22 (KOX 15)	1.68	0.40
PRKD2	protein kinase D2	1.63	0.38
AKAP8	A kinase (PRKA) anchor protein 8	1.61	0.37
DDX21	DEAD box polypeptide 21	1.61	0.43
WNT5B	Wnt 5B	1.60	0.45
FOXK2	forkhead box K2	1.60	0.53
RAB8A	RAB8A, member RAS oncogene family	1.59	0.57
DDX3X	DEAD box polypeptide 3, X-linked	1.58	0.71
FAM82A2	Family with sequence similarity 82, member A2	1.57	0.53
RRP1	ribosomal RNA processing 1 homolog	1.57	0.34
TNPO2	transportin 2	1.55	0.42
METTL1	methyltransferase like 1	1.53	0.38
FAM168B	family with sequence similarity 168, member B	1.53	0.64
MECR	mitochondrial trans-2-enoyl-CoA reductase	1.52	0.59
NFKBIE	NFkappaB inhibitor epsilon	1.50	0.27
RARA	Retinoic acid receptor alpha	1.49	0.16
BIRC5	survivin	1.46	0.58
PLA2G4A	Phospholipase A2, group IVA	1.46	0.17
YRDC	yrdC domain containing, mitochondrial protein	1.46	0.62
NOP16	NOP19 nucleolar protein homolog	1.40	0.20
HRAS	HRAS	1.40	0.54
PNO1	partner of NOB1 homolog	1.39	0.06
SEPHS1	selenophosphate synthetase 1	1.36	0.47
UBIAD1	UbiA prenyltransferase domain containing 1	1.34	0.38
ADO	2-aminoethanethiol (cysteamine) dioxygenase	1.34	0.25
DKC1	dyskerin	1.34	0.48
AKAP1	A kinase (PRKA) anchor protein 1	1.32	0.42
APOL2	apolipoprotein L, 2	1.31	0.53
RPP25	ribonuclease P/MRP 25 kDa subunit	1.30	0.15
HS3ST1	heparan sulfate (glucosamine) 3-O-sulfotransferase 1	1.30	0.11
NOLC1	nucleolar and coiled-body phosphoprotein 1	1.27	0.59
PBK	PDZ binding kinase	1.26	0.38
IL12A	interleukin 12A	1.25	0.31

Genes which are up-regulated in Myeloma patients (fold change (FC) > 1.25) and down-regulated with combination treatment in MM.1S myeloma cells (FC< 0.80) are listed.

Table 2

Genes with GR binding sites and negative correlation between Myeloma patients and combination drug treatment groups.

Gene Symbol	Name	FC MM down	FC treated up
FRAT1	frequently rearranged in advanced T-cell lymphomas	0.44	29.14
TLE1	transducin-like enhancer of split 1	0.21	23.72
TXNIP	thioredoxin interacting protein	0.27	13.79
PIK3IP1	PI3K interacting protein 1	0.47	10.91
CXCR4	chemokine (CXC motif) receptor 4	0.45	5.01
KLHL24	Kelch-like family member 24	0.55	4.43
TMEM2	transmembrane protein 2	0.64	3.73
SOS1	son of sevenless homolog 1	0.75	3.38
SOX30	SRY(sex determining region Y)-box 30	2.91	3.37
RNF122	ring finger protein 122	0.77	3.27
OGFRL1	opioid growth factor receptor-like 1	0.76	2.43
C1orf56	Methylated in normal T-lymphocytes (MENT)	0.52	2.23
TAX1BP3	Tax1(human T-cell leukemia virus type 1)	0.52	2.05
AIDA	axin interactor, dorsalization-associated	0.64	1.93
STK38	serine/threonine kinase 38	0.32	1.90
CDC42SE1	CDC42 small effector 1	0.57	1.88
LAPTM5	lysosomal protein transmembrane 5	0.14	1.78
PPOX	protoporphyrinogen oxidase	0.79	1.71
SLC2A3	GLUT3	0.76	1.69
WIPI2	WD repeat domain, phosphoinositide interacting 2	0.65	1.59
OSBPL2	oxysterol binding protein-like 2	0.65	1.49
AIM2	absent in melanoma 2	0.57	1.38

Genes which are down-regulated in Myeloma patients (fold change (FC) < 0.80) and up-regulated in MM.1S myeloma cells (FC> 1.25) are listed.

GR modulates anti-apoptotic and nucleotide metabolism pathways

We hypothesize that these candidate transcriptional targets of GR may be acting as nodes to regulate the downstream events of glucocorticoid-induced cell death. We performed GSEA of GR-regulated genes using the set of genes that are down-regulated by combination treatment and have potential GR binding sites to gain a global understanding of pathways potentially involved in cell death that may be directly modulated by glucocorticoids (Figure 4 A). Gene sets involving the cell cycle, nucleotide metabolism, transcription and major growth-promoting signaling pathways were enriched for down-regulated genes, indicating these mechanisms are repressed by GR actions following glucocorticoid treatment.

From our list of candidate genes (Table 1), we chose two genes to validate their role in glucocorticoid-induced cell death, as proof of principle for our identification approach. We chose two distinct genes within down-regulated pathways identified in the GSEA analysis, BCL2L1 and RRM2, to validate as GR-targeted, death-inducing genes. B-cell lymphoma-extra-large (BCL2L1/BCL-XL) is a member of the BCL2 family of anti-apoptotic proteins, and inhibitors of the BCL2 family and inhibition of this protein specifically has been demonstrated to cause death in myeloma cell lines [Chauhan et al., 2007]. In our GEP of combination drug treatment, BLC2L1 was one of the most repressed genes, with a 0.20-fold change of expression compared to the control. Additionally, BCL2L1 expression is increased 1.78-fold in myeloma patients compared to healthy controls (Table 1). Strong GR binding was identified within 10kB of the promoter (see Supplementary Material Figure 2). To confirm regulation of BCL-XL by combined Dex and LY294002 treatment, we measured mRNA by qRT-PCR and protein expression by immunoblot (Figure 4 B and D). We observed a decrease in both mRNA and protein expression following glucocorticoid treatment that was further decreased by the addition of PI3K inhibitor LY294002.The validation of BCL2L1 as a primary target for glucocorticoid regulation involved in glucocorticoid-induced cell death validates our target identification process.

Next, we tested our identification process with a gene in a complementary pathway that had not yet been identified as a primary GR target. Ribonucleotide reductase M2 (RRM2) catalyzes the formation of deoxyribonucleotides from ribonucleotides, which is essential for nucleic acid production for DNA replication and repair [Zhou et al., 2013]. Small molecule inhibitors are in development targeting RRM2 as a potential therapeutic target in breast cancer [Shah et al., 2014]. In our screening, we observed RRM2 expression is significantly overexpressed in MM patients (4.47-fold change, Table 1); RRM2 expression was reduced upon combination treatment compared to the control in our gene expression array (0.50-fold change), and contains a potential GR binding site within 100kB of the transcription start site (see Supplementary Material Figure 2). The potential GR binding site for RRM2 is in the promoter region for the KLF11 gene. None of our treatments resulted in a change in KLF11 expression (data not shown); therefore, we believe it is reasonable to assign this GR binding site to regulation of RRM2. The reduction of RRM2 mRNA and protein levels following combined Dex and LY294002 treatment was confirmed by qRT-PCR and immunoblot (Figure 4 C and D). BCL-XL and RRM2 are key players in distinct pathways essential for cell survival, and therapeutic targeting of both proteins may produce a more robust response in myeloma patients than targeting each pathway individually. The dramatic decrease in transcript and protein expression of both BCL-XL and RRM2 following Dex treatment, as well as their potential direct transcriptional control by GR, strongly suggests their involvement in the Dex-induced cell death mechanism.

BCL-XL and RRM2 are required for Dex-sensitive and -resistant multiple myeloma cell survival

We next sought to determine if the GC-induced repression of BCL-XL and RRM2 contributes to Dex-induced cell death of myeloma cells. MM.1S cells and MM.1RL cells, which are Dex-resistant due to loss of GR, were transfected with siRNA targeting either BCL-XL or RRM2 to deplete protein expression. Knock-down efficiency was confirmed by qRT-PCR (data not shown) and immunoblotting (Figure 5 A and D). Using densitometry of immunoblots, knock-down efficiency was determined after 96 h, the time point when cell viability was measured. BCL-XL expression was decreased 80% compared to the non-specific control (p<0.05) in MM.1S cells and decreased 75% in MM.1RL cells (p<0.05). RRM2 expression was decreased approximately 40% compared to the non-specific control in both MM.1S and MM.1RL cells (p<0.05). BCL-XL knock-down caused an average 1.5-fold increase in death of MM.1S cells and an average 1.9-fold increase in MM.1RL death (Figure 5 B). It is interesting to note that Dex-resistant cells were slightly more sensitive to BCL-XL loss. RRM2 knock-down resulted in a 1.5-fold increase in MM.1S cell death and a 1.4-fold increase in MM.1RL cells (Figure 5 E). Myeloma cells are sensitive to BCL-XL and RRM2 loss, indicating both proteins are required for cell survival. Moreover, loss of either protein induces cell death in MM.1RL cells, suggesting their usefulness as therapeutic targets for Dex-resistant myeloma.

To investigate whether regulation of BCL-XL and RRM2 can modulate the Dex-induced cell death mechanism, each protein was knocked-down by siRNA in MM.1S cells, and cells were treated with a range of Dex concentrations. Knock-down of either protein sensitized MM.1S cells to Dex-induced cell death at concentrations approximately 50% less than the amount of Dex required for maximal apoptosis without gene knock-down (Figure 5 C and F). No sensitization was seen in Dex-resistant MM.RL cells (data not shown). Altogether, these results suggest that down-regulation of BCL-XL and RRM2 by activated GR contribute to the Dex-induced myeloma cell death response, and targeting these proteins may provide a novel replacement therapy for Dex-resistant myeloma.

Discussion

Although glucocorticoids can be an effective therapeutic for the treatment of myeloma, efficacy of this approach is limited by loss of GR expression. Therefore, delineation of the GR targets which direct myeloma cell death and therapeutic intervention at that level may obviate the need for an active GR in advanced myeloma. Here, we have taken a systematic approach to the genomic dissection of GR-induced cell death specific to myeloma through the integration of multiple data sets. To focus on the GR actions that induce myeloma cell death, we co-treated with a PI3K inhibitor which we have previously shown to induce synergistic cell death in myeloma cells. Utilizing publicly- available data sets from myeloma patients, we assessed clinical relevance of our myeloma cell line GEP and finally examined genomic GR binding sites in myeloma cells to identify potential primary GR targets that could be clinically-relevant for treatment of myeloma. This global genomic dissection of gene expression changes following treatment with glucocorticoids in combination with PI3K inhibitors has identified druggable pathways that can be exploited for therapeutic advantage in advanced myeloma.

Microarray analyses of glucocorticoid-induced changes in gene expression have been performed on a number of tissues and cell lines in an attempt to identify key genes regulated for development of new therapeutic approaches. We observed approximately twice as many down-regulated genes as up-regulated at the time point of 3 h of incubation. Chauhan et al. [Chauhan et al., 2002a] also observed a preponderance of GC-induced down-regulated genes in myeloma cell lines at 3 h of incubation. In mice injected with synthetic GCs 3 h prior to sacrifice, examination of GC-regulated genes in the liver indicate 70% of the genes were down-regulated compared to 30% up-regulated [Phuc Le et al., 2005]. In contrast, in breast cells [Wu et al., 2004] or in human T-cell lymphoma cell lines or mouse thymus [Wang et al., 2003], the majority of GC-regulated genes were up-regulated rather than down-regulated. Because glucocorticoids impact many tissue types and regulate a wide spectrum of biological responses, it is not unexpected that there should be differences in gene expression between tissue types.

Inhibition of BCL2L1 has been previously observed to cause cell death in myeloma, and therapeutics are in development to target this protein in myeloma [Chauhan et al., 2007]. Our identification of BCL2L1 as a GR target in myeloma and subsequent confirmation of its participation in GC-induced cell death validates our genomics approach. In addition to myeloma, overexpression of BCL-XL has been shown to play a role in additional cancers. For example, BCL-XL may be involved in resistance to common breast cancer therapeutics, such as 5-fluoruracil and doxorubicin [Wang et al., 2015]. Another study demonstrated that reduction of BCL-XL in colorectal cancer improves KRAS-dependent apoptosis resistance [Okamoto et al., 2015]. Inhibition or reduction of BCL-XL expression also induces colorectal cancer cell apoptosis and reduces invasion [Liu et al., 2014; Ye et al., 2015] . Inhibition of BCL-XL in small-cell lung cancer sensitized cells to radiation therapy [Loriot et al., 2014]. Overall, overexpression of BCL-XL is critical to cell survival in various cancer types. Our study confirms these findings and provides further evidence for the use of BCL-XL inhibitors to treat myeloma and other cancers.

Our identification of RRM2 as a GR target and possible therapeutic target in myeloma is novel. Targeting of RRM2 has not been previously explored in myeloma, and this approach could complement current therapeutics. RRM2 is overexpressed in nearly all cancer types that are publicly-available in the Oncomine database. This includes bladder, brain, breast, cervical, colorectal, esophageal, gastric, head/neck, kidney, liver, lung, ovarian, pancreatic and prostate cancers, as well as lymphoma, melanoma and sarcoma. Recent publications implicate the overexpression of RRM2 with higher grade tumors and poor patient survival. For example, RRM2 may act as a prognostic marker for breast cancer and for tamoxifen resistance [Putluri et al., 2014]. RRM2 may also serve as marker for cervical cancer, with its overexpression correlating with higher grade cancers and predicting poor survival [Su et al., 2014]. Another study implicated RRM2 in promoting angiogenesis in cervical tumors [Wang et al., 2014]. RRM2 was found to be overexpressed in an array of gastric cancer cell lines and in primary gastric tumor cells and correlated with tumor grade and poor survival time. Knock down of RRM2 with siRNA (similar to our study) in gastric cancer cell lines decreased proliferation, colony formation and invasion and induced apoptosis [Kang et al., 2014]. Another study found RRM2 expression may act as a predictor of prostate cancer recurrence [Huang et al., 2014]. The current literature validates our findings regarding the importance of RRM2 in cell survival. Furthermore, our study supports the growing body of evidence that RRM2 may be an important bio-marker and target in cancers.

In summary, we have developed a biology-guided approach to screening for clinically-relevant GR targets and used this approach to identify 57 candidate genes as potential direct targets of GC actions which are over-expressed in myeloma patients and inhibited by combined treatment with GCs and inhibition of PI3K. As proof of principle for this screening process, we have validated two of the targets as being required for myeloma cell survival and whose inhibition sensitizes cells to GC-induced cell death. Overall, our results identify nearly 400 candidate genes under the transcriptional control of GR, the majority being down-regulated upon Dex treatment. We confirmed two candidate transcriptional targets of GR, namely BCL-XL and RRM2, that are up-regulated in MM patients, are required for cell survival, and inhibition of either protein causes death in myeloma cells sensitive and resistant to Dex, which validated our screening pipeline. Our data suggest targeting these proteins as mono-therapies or in combination may provide novel therapeutic approaches for multiple myeloma and provide the means to identify additional GR transcriptional targets whose inhibition contributes to Dex-induced cell death. Future studies will focus on investigating the therapeutic potential of those genes identified in this screening process.