In vitro and ex vivo gene expression profiling reveals differential kinetic response of HSPs and UPR genes is associated with PI resistance in multiple myeloma

Extensive inter-individual variation in response to chemotherapy (sensitive vs resistant tumors) is a serious cause of concern in the treatment of multiple myeloma (MM). In this study, we used human myeloma cell lines (HMCLs), and patient-derived CD138+ cells to compare kinetic changes in gene expression patterns between innate proteasome inhibitor (PI)-sensitive and PI-resistant HMCLs following test dosing with the second-generation PI Ixazomib. We found 1553 genes that changed significantly post treatment in PI-sensitive HMCLs compared with only seven in PI-resistant HMCLs (p < 0.05). Genes that were uniquely regulated in PI-resistant lines were RICTOR (activated), HNF4A, miR-16-5p (activated), MYCN (inhibited), and MYC (inhibited). Ingenuity pathway analysis (IPA) using top kinetic response genes identified the proteasome ubiquitination pathway (PUP), and nuclear factor erythroid 2-related factor 2 (NRF2)-mediated oxidative stress response as top canonical pathways in Ix-sensitive cell lines and patient-derived cells, whereas EIF2 signaling and mTOR signaling pathways were unique to PI resistance. Further, 10 genes were common between our in vitro and ex vivo post-treatment kinetic PI response profiles and Shaughnessy’s GEP80-postBz gene expression signature, including the high-risk PUP gene PSMD4. Notably, we found that heat shock proteins and PUP pathway genes showed significant higher upregulation in Ix-sensitive lines compared with the fold-change in Ix-resistant myelomas.


Introduction
Multiple myeloma (MM) remains an incurable disease with 5-year survival rate of 48.5% (NCI-SEER Cancer statistics) 1 . Extensive inter-individual variation in response to chemotherapy is a serious cause of concern in the treatment of myeloma 2,3 . Not all patients respond equally well to treatment and those who do often develop resistance over the course of treatment. Drug resistance may therefore be categorized into: (1) innate resistance already present in drug-naive patients who never respond to treatment, or (2) emerging/acquired resistance where a patient's tumor ultimately undergoes relapse or "acquires" the ability to resist therapy in course of treatment despite good response to initial treatment 3,4 . Proteasome inhibitors (PIs) are standard-of-care chemotherapeutic agents for myeloma that impede tumor metastasis and angiogenesis by accelerating unfolded protein response (UPR) and by interfering with the nuclear factor-κB-enabled regulation of cell adhesion-mediated drug resistance 3,5,6 .
Alterations in gene expression levels have been shown to be associated with response to cancer drugs, including PIs in myeloma 7,8 . In our previous study, we used a panel of 50 human myeloma cell lines (HMCLs) representing the broad spectrum of biological and genetic heterogeneity of myeloma, and generated in vitro drug response profile for the four PIs: bortezomib (velcade/Btz), carfilzomib (kyprolis/Cfz), oprozomib (Opz), and ixazomib (ninlaro/Ix) as single agents. As PI sensitivity in these HMCLs were highly correlated (given a common target), this suggested that any of these four PIs could be used as surrogates in models to understand common PI resistance. Further, we used machine learning-based computational approaches to derive a gene expression profiling (GEP) signature predictive of baseline PI response in myeloma 9 .
Here, we were interested in understanding how gene expression patterns may be differentially influenced in sensitive vs resistant myelomas following exposure to PIs (kinetic PI response). Therefore, to characterize differences in kinetic PI response, we compared transcriptome profiles of HMCLs in vitro as well as patient primary tumor cells ex vivo representative of highly PI-sensitive and PI-resistant myelomas between baseline and post treatment following test dosing with the model PI drug ixazomib (Ix) 3 . Our results provide significant insights into the wide heterogeneity by identifying key differences in PI-specific kinetic responses.

Drugs
Ix and Bz was procured from Takeda (Takeda Pharmaceuticals Inc., Deerfield, IL, USA). Cfz and Opz were obtained from Amgen (Thousand Oaks, CA, USA). Drugs were dissolved in dimethyl sulfoxide and stored at −20°C.

Cell lines
Twelve (12) myeloma cell lines-established, characterized, and authenticated by our collaborators at The Translational Genomics Research Institute (TGen) and NIH representing variabilities in PI-response in myeloma patients were selected from our HMCL panel described previously 10 .

Patient samples
CD138-selected plasma cells were derived from myeloma patients enrolled in a phase-2 Ix clinical trial at Mayo Clinic, Minnesota (Mayo-Ix; NCT01415882). These were relapsed myeloma patients had less than six cycles of prior treatment with a PI-based regimen and were not refractory to the PI drug bortezomib (Bz) 9 .
In vitro chemosensitivity assays Cells (HMCLs and patient-derived primary cells) were counted using countess automated cell counter (Invitrogen, Carlsbad, CA, USA) and seeded in 96-well plates at a concentration of 4 × 10 5 cells per ml. Twenty-four (24 h) hours later, Ix, representative PI drug, was added as a single agent in increasing concentrations. Forty-eight (48 h) hours post Ix treatment, CellTiter-Glo luminescent cell viability assay (Promega, Madison, WI, USA) was used to perform cell viability assays according to manufacturer's instructions using Synergy 2 Microplate Reader (BioTek, Winooski, VT, USA). Survival curves were generated, percent survival values were calculated at each drug concentration and normalized to untreated controls. Half-maximal inhibitory concentration (IC 50 ) values were computed by nonlinear regression using sigmoidal dose-response equation (variable slope). Area under survival curve (AUSC) was calculated using trapezoidal rule 9 .
Cell plating for kinetic Ix treatment and RNA-sequencing HMCLs and patient-derived primary myeloma tumor cells were plated at a density of 4 × 10 5 cells per mL. In all, 15 nM (IC 50 of the most sensitive HMCL) and 40 nM (median IC 50 of patient cells) of Ix were added to HMCLs and PMCs, respectively, after 24 hours of incubation. Cells were collected 24 hours post-treatment as baseline (0 nM/ untreated) and post-treatment (treated). High-quality RNA was extracted from untreated and treated myeloma cells using QIAshredder and RNeasy kit (Qiagen). RNA concentration and integrity were assessed using the Nanodrop-8000 and Agilent 2100 Bioanalyzer and stored at −80°C. RNA integrity number/RIN threshold of eight was used for RNA-seq analysis. RNA-seq libraries were constructed using Illumina TruSeq RNA sample Preparation kit v2. Libraries were then size selected to generate inserts of~200 bp. RNA sequencing was performed on llumina's HiSeq 2000 next-generation high-throughput sequencing system using 50 bp paired-end protocol with depth of >20million reads per sample. Average quality scores were thoroughly above Q30 for all libraries in both R1 and R2.

RNA-Seq data processing
RNA-seq data were normalized and fragments per kilobase million values were used in further analysis using Partek Genomics Suite and Galaxy data analysis software, an open source, web-based platform that provides tools necessary to create and execute RNA-seq analysis. In brief, RNA-seq data analysis pipeline was developed using Galaxy software workflow. Quality control (QC) check on the RNA-seq raw reads was performed using FastQC tool followed by read trimming to remove base positions that have a low median (or bottom quartile) score. Tophat2 tool mapped processed RNA-seq reads to the hg19 human genome build. Picard's CollectInsertSizeMetrics tool was applied on the initial tophat2 run to obtain estimated insert sizes, which was then used to calculate mean inner distance between mate pairs (mean = esti-mated_insert-size−2×read_length). Tophat2 was re-run using corrected mean value and Cufflinks tool was to assemble the reads into transcripts.

Bioinformatics and statistical analysis
Gene expression data were pre-processed using Galaxy platform and filtered using the following criteria: genes with mean FPKM < 1 were removed. GEP data were analyzed further using a combination of Galaxy and Partek Genomics Suite v7.0 to perform differential expression testing to identify GEP signatures of kinetic PIresponse. Mean fold-change >|1| and p < 0.05 was considered as threshold for reporting significant differential gene expression.
All statistical analyses were performed using R for statistical computing and graphics, version 3.4.2, and GraphPad Prism version 7.0. All tests were two-sided and p < 0.05 was considered statistically significant.
Differentially expressed (DE) gene analysis was performed between two groups of gene expression data sets (e.g., treated vs untreated). For HMCL data, we first detected DE genes for: (1) all the 12 HMCLs combined between treated (15 nM) and untreated groups; (2) six sensitive HMCLs and six resistant HMCLs between treated (15 nM) and untreated groups; separately.
Similarly, for patient data, we also performed DE analysis for: (1) eight patients and (2), separately, among four sensitive patients and four resistant patients between treated (40 nM) and untreated groups. Owing to the small sample size, we use limma, an empirical Bayesian method, to detect the DE genes. The advantage of limma compared with traditional t test is that limma provides a moderated t test statistic by a shrinking the variance statistics, therefore, improves the statistical power. Thus, we obtained the p value for each DE gene and further provided FDR based on the p value using Benjamini-Hochberg procedure. Genes with P < 0.05 are considered as DE genes.

Ingenuity pathway analysis
Ingenuity pathway analysis (IPA) software (QIAGEN) was be used to identify the most significantly affected molecular pathways predicted to be activated or inhibited in response to secondary drugs and in stem-cell based PI resistance modeled on list of significantly differentially regulated genes. IPA is a web-based software application that integrates and interprets the data derived from differential mRNA expression analysis to: (i) identify the most significantly affected canonical pathways predicted to be activated or inhibited; (ii) predict upstream regulator molecules like microRNA and transcription factors, which may be causing the observed gene expression changes; (iii) analyze downstream effects and biological processes that are increased or decreased; (iv) predict causal networks; and (v) perform predictive toxicology analysis using toxicogenomics approaches (IPA-Tox) 11 .

Western blotting
Myeloma cells were cultured with (treated) or without (untreated) Ix (15 nm), harvested, washed, and lysed using radioimmunoprecipitation assay lysis buffer containing 50 mM Tris-HCl, pH 7.5, 150 mM NaCl, 1% NP40, 5 mM EDTA, 1 mM DTT, phosphatase, and protease inhibitors cocktail (Sigma) and incubated on ice for 15 min. Samples were then centrifuged at 14,000 rpm at 4°C for 30 mins. The supernatant was then aspirated and quantified using Pierce BCA Protein Assay Kit (Thermo Scientific). Samples were solubilized in sodium dodecyl sulfate polyacrylamide gel electrophoresis sample buffer, and equal amounts of protein were loaded per lane of 10% SDS-PAGE and transferred onto polyvinylidene difluoride membranes (Millipore; Billerica, MA). Membranes were blocked in tris-buffered saline (TBS) with SuperBlock blocking buffer (Thermo Fisher). Membranes were incubated with primary antibodies for HSP40 (DNAJB1), HSPA1b (HSP70), and HSP90alpla and secondary antibodies in TBS with 0.2% Tween 20 and 2.5% bovine serum albumin. Immunoreactivity was detected by Chemiluminescent HRP substrate (Bio-Rad) and the exposed image was captured using a ChemiDoc MP Imaging System (Bio-Rad). Densitometry analysis was performed (in triplicates) using Image J software.

In vitro and ex vivo cell survival assays and patient outcome
Response to Ix in HMCLs and patient-derived human myeloma cells are provided in Fig. 1a, b. As some of the test cells (PMCs and HMCLs) did not achieve IC 50 , we used AUSC values to assess relative sensitivity to Ix and performed binary grouping into Ix-sensitive and Ixresistant subgroups, as highlighted in the figure. In vitro response to Ix in eight primary myeloma patient PCs is shown in Fig. 1c, d. Patient outcome data in terms of response to PI-containing chemotherapy regimen showed all the patients in the low Ix AUC category (in vitro PIsensitive) presented with some form of response (CR/PR/ VGPR), whereas only one patient in the PI-resistant (high in vitro cytotoxicity) category showed complete response (Fig. S1a). Median progression-free survival values were 371 days for the PI-resistant patients and 858.5 for the Ixsensitive patients (Hazards ratio >1.5) (Fig. S1b).

Differential kinetic GEP of HMCLs
In order to identify kinetic changes in gene expression patterns common to all the HMCLs we first performed a combined analysis of all 12 HMCLs (Ix-sensitive + Ix-resistant) considered together between baseline and 24 h post Ix treatment of 15 nM. A total of 583 DE genes displayed p value > 0.05 (|fold-change | >1; p < 0.05) and 304 genes at |fold-change | >1.5 (p < 0.05). Figure 2a shows a heatmap of top genes in 12 myeloma HMCLs (six sensitive, six resistant) that showed significant change in expression of in response to Ix treatment, 24 hours following drug exposure (Ix-15 nM). IPA analysis of the 583 genes identified the proteasome ubiquitination pathway showing differential gene expression of top genes in 12 myeloma cell lines (six sensitive, six resistant) that showed significant de-regulation (|fold-change| > 1.5; p < 0.05) in response to Ixazomib treatment, 24 h following drug exposure (Ixazomib-15nM). Columns represent cell lines and rows represent genes. Prior to hierarchical clustering, gene expression values were filtered (samples with mean FPKM < 1 or max FPKM < 1 were removed) and z score normalized. Columns are ordered by Ix-response followed by Ix-treatment. b Heat map of Top 50 genes that changed significantly (p < 0.05) between baseline and post-treatment Ix-sensitive myeloma cell lines (n = 6). Columns are ordered by Ix-treatment. c Heat map of top genes (p < 0.05) that changed significantly between baseline and post-treatment Ix-resistant myeloma cell lines (n = 6). Columns are ordered by Ix-treatment.
However, we observed that the GEP data were highly influenced by the DE genes in sensitive HMCLs. Therefore, to understand major differences in differential kinetic gene expression patterns between Ix-sensitive and Ix-resistant lines, analysis of untreated baseline and post 15 nM Ix treatment was performed on the six most Ixsensitive and six most Ix-resistant HMCLs, separately (Fig. S2).
When Ix-sensitive lines were considered alone, 1553 genes changed significantly post-treatment with any change (|fold-change| > 1; p < 0.05). In total, 1167 genes at |fold-change| >1.5 (p < 0.05). The top 50 most highly regulated genes for each condition are listed in Table S1 a-d. A heatmap for the top 50 genes in Ix-sensitive lines is provided in Fig. 2b.
On the other hand, kinetic response profiling in Ixresistant HMCLs revealed at p < 0.05, only seven genes showed any significant change (|fold-change| >1; p < 0.05) (Fig. 2c). Comparison of the distinct kinetic response profiles of the Ix-sensitive and Ix-resistant HMCLs revealed three common DE genes (PIR, TNKS, PSMD2) that changed significantly in both the lists at p < 0.05 as shown in the venn diagram in Fig. 3a. Four genes were unique to the Ix-resistant HMCLs only-ANO10, ERCC6, LRRC28, PPIP5K1.
IPA analysis of the most DE genes in Ix-sensitive and Ixresistant HMCLs in response to Ix-test dosing showed that the top canonical pathway was the PUP (p = 2.64E-23 and 1.96E-08 for sensitive and resistant HMCLs, respectively) (Fig. S3). IPA-Tox analysis in Ix-sensitive HMCLs revealed 33 genes overlapped with the top pathway NRF2mediated oxidative stress response (P = 6.26E-05), whereas nine genes overlapped with oxidative stress (P = 9.51E-03).
A comparison of the major PUP and stress response pathway genes between the Ix-sensitive and Ix-resistant lists revealed that the three of the top four genes that showed significantly different kinetic response were DNAJB1 (HSP40), HSPA1B (HSP70), and HSP90AA1 with sensitive to Resistant fold-change ratios of 31.6, 6.3, and 4.4, respectively. Furthermore, all the seven PUP genes (BIRC6, PSMA1, PSMB3, PSMB5, PSMC4, PSMC6, PSMD2) that were common between the two lists that showed >1.5-fold higher upregulation in Ix-sensitive lines (Table 1).
On the other hand, kinetic gene expression analysis of Ix-resistant PMCs with lower four AUSCs demonstrated 683 genes with significant |fold-change| (>1) between baseline vs post-treatment with significant p value (p < 0.05). Heatmap for the top 50 most differentially regulated genes is represented in Fig. 4c. The gene lists are provided in Table S1c (Ix-sensitive) and S1d (Ix-resistant). When compared, 107 kinetic genes were common between the Ix-sensitive and Ix-resistant patients, whereas 576 genes were unique in Ix-resistant PMCs only, as demonstrated in the venn diagram (Fig. 3b).
Supplementary figure S4 shows venn diagrams obtained from the kinetic response profiling of all the subgroups (HMCLs and patients; PI-sensitive and resistant) and depicting the comparison of all the significant gene lists (p < 0.05).

Ix induced higher kinetic upregulation of HSPs in sensitive cells compared to resistant cells
As our in silico study showed that Ix-sensitive cells have significantly elevated of HSP70b levels post treatment compared with Ix-resistant cells, we evaluated heat shock protein (HSP) expression in a panel of PI-sensitive and -resistant cells at protein level to validate further the genomic data. Following treatment with Ix for 24 h at 15 nM concentration, Ix-sensitive cells showed higher expression of HSP40 (DNAJB1), HSPA1b (HSP70), HSP90alpla compared with the increase in expression in Ix-resistant cells (Fig. 5). Thus, western blotting confirmed fold-change differences in kinetic response between sensitive vs resistant lines.

Discussion
Extensive inter-individual variation in response to chemotherapy is a serious cause of concern in the treatment of MM. Using a panel of HMCLs, we have earlier presented a GEP signature that could not only distinguish good and poor PI-response in the HMCL panel, but was also successful in stratifying exceptional response to PIbased chemotherapy in myeloma clinical data sets 9 .
Although a recent study used CRISPR-based screening approaches to identify genes associated with PI response in myeloma, no study so far has combined in vitro and ex vivo pre-and post-treatment gene expression profiles to compare PI treatment-induced gene signatures between sensitive and resistant myelomas. In this study, we used six most PI-sensitive and six most innate PIresistant myeloma lines from our HMCL panel above and identified in vitro the kinetic changes in gene expression patterns following Ix test dosing. Subsequently, we also performed an ex vivo analysis of kinetic gene expression changes between untreated and treated myeloma tumor cells derived from ongoing Ix-containing clinical trials. Interestingly, our ex vivo Ix-sensitivity data corroborated with patient PI response from follow-up data. Our analysis identified sets of genes that changed significantly post-treatment in Ix-sensitive and comparatively Ixresistant myelomas with significant similarities and some notable differences. These therefore represent distinct signatures of kinetic response to Ix in sensitive and innate resistant myeloma and highlight major differences in pathway regulation associated with differential treatment responses. Notably, we find a number of genes and pathways that have been previously associated with response and resistance are, in fact, differentially expressed in response to the drug in sensitive and resistant myeloma PCs.
We observed strikingly lower number of genes that changed significantly post-treatment in Ix-resistant HMCLs when compared with the Ix-sensitive HMCLs. A large number of these under-perturbed genes belonged to the PUP, suggesting a possible mechanism of PI resistance where this pathway is not sufficiently affected following PI treatment in relatively resistant HMCLs. For example, the gene PSMC6, the lack of which has been associated with cells having reduced chymotrypsin-like proteasome activity by PIs 12 . In our study, the upregulation of PSMC6 was more than twofold higher in the sensitive lines compared with Ix-resistant HMCLs. In addition to the PUP, the kinetic genes that were differentially responsive to PI treatment in Ix-sensitive HMCLs and patient-derived primary cells included NRF2mediated oxidative stress response and oxidative stress response genes, that also corroborated with our current finding on kinetic GEP response in HMCLs reported here 13 . The NRF2-mediated oxidative stress response pathway has been shown earlier to control ER-stress induced apoptosis in myeloma in UPR 14 . Furthermore, among Ix-sensitive models, XBP1 inhibition and NFE2L2 activation were predicted as top upstream regulators based on the observed significantly regulated gene expression changes. Earlier studies have shown that XBP1 regulation is associated with significant inhibition of myeloma tumors which in turn is associated with chronic endocrine reticulum stress from UPR 15 . Thus, it was particularly interesting that XBP1 expression is differentially regulated in response to Ix. As described below, we found other genes whose expression has been associated with response and resistance are, in fact, differentially regulated by treatment with PIs.
The HERC1 gene belongs to the HERC family that has ubiquitin ligase activity dependent on the protein substrates, HERC proteins may act as a tumor suppressor or oncoprotein in specific cancer types. For example, HERC4 is believed to contribute to carcinogenesis of solid tumors such as lung cancer, but it suppresses the proliferation of myeloma cells 14 .
RICTOR is a scaffold protein that activates mTORC2 complex constituting mammalian Target of Rapamycin or mTOR, an intracellular serine/threonine kinase involved in intracellular cell survival, tumor growth, and drug resistance. Regulation of mTOR2 has been shown to influence myeloma angiogenesis 16 .
Major microRNA genes expressed uniquely in Ixresistant patient cells included miR-16-5p (activated). An earlier study that used microarray profiling of tumorderived exosomal RNA showed downegulation of miR-16-5p in circulating exosomes of patients resistant to PIs 17 . A number of studies have suggested that tumor-derived exosomal miRNA have a crucial role in the cross-talk between the tumor and its microenvironment and is also predictive of drug resistance in myeloma patients 17,18 . One of the most notable findings of our study was that eIF2 signaling was the most significantly affected canonical pathway unique to the Ix-resistant cells. eIF2 signaling is involved in regulating ternary complex formation and protein synthesis 19 . Earlier studies have shown that inhibition of the 26 S proteasome results in a rapid decrease in the rate of protein synthesis owing to phosphorylating alpha subunit of the eukaryotic translation initiation factor 2 (eIF2α) by the heme-regulated inhibitor kinase (HRI) 20 . Further, dysregulation of the PERK-eIF2α pathway was involved in the efficiency of and sensitivity to target PIs [21][22][23] .
MYC genes are highly deregulated proto-oncogenes in human cancer. Aberrant MYC expression is frequently associated with cancer cell vulnerability and poor prognosis 24 . Several studies have therefore proposed the therapeutic targeting of aberrant MYC expression in hematological cancers including myeloma 24 . In our study, we found the MYC family genes were predicted to be inhibited based on the gene expression profile of the Ixresistant cells. This may not be surprising, as we have previously demonstrated that the impact of myc expression on survival is very context dependent. High levels of myc expression in our previously reported transgenic mouse model in fact was associated with high cell death in the plasma cells expressing it 25 . Thus, the decrease in expression of myc pathway genes may well be associated with resistance to cell death.
IL-17A has been shown to influence myeloma cell proliferation and tumor growth through IL-17A receptor (IL-17RA) 4 expressed on tumor cells 26 . Moreover, anti-IL-17A mAb could significantly inhibit myeloma cell-growth and survival both in vitro and in vivo 27 . Interestingly, in our IPA analysis, IL-17A signaling was the top canonical pathway for the 149 genes uniquely expressed in PI-resistant patient cells following treatment with PI.
IPA's upstream analysis based on the genes found consistently downregulated with increase in dose (0 vs 20 nM vs 40 nM) predicted inhibition of the gene MITF as the top upstream regulator (Z score = −4.191; P = 1.08E-04). Earlier studies have shown that MITF is a direct target of miR-137, and miR-137-targeting of MITF regulates drug sensitivity in myeloma cells by reducing c-MET expression and suppressing AKT phosphorylation, accompanied by an increase in p53 expression 28 . Another study has found miR-137 directly targets and negatively regulates the protein expression of Enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2)-a histone methyltransferase that we have earlier shown to influence B-cell development, disease progression and treatment outcomes in myeloma [29][30][31] .
When the DE gene sets representing signatures of Ix were uploaded as separate data sets into the IPA/ IPA software, the following revelations came forward. Although both the gene lists (Ix-sensitive and Ix-resistant myelomas) identified the PUP as the top canonical pathways, the overlap was around two times higher for the subgroups that respond to Ix treatment compared with the Ix-resistant subgroups. Furthermore, the UPR pathway was ranked higher among Ix-sensitive subgroup compared with the top canonical pathways observed in the subgroups that were relatively more resistant to Ix/PIs treatment. We also found, both in HMCLs and patient cells, that the treatment-induced increase in the involvement of the heat shock proteins and PUP was PI-dose dependent (not shown). HSPs play a crucial role in a number of biological processes involved in the maintenance of cellular protein homeostasis in a normal cells as wells as in cancer cells, including protein folding, cellular proliferation, differentiation, survival, metastasis, invasion, and angiogenesis [32][33][34][35][36][37][38][39] . In myeloma cells, HSPs play important role by activating the cytoprotective heat shock response during PI treatment with HSP27, HSP40, and HSP70 being highly unregulated following Btz treatment 39,40 . Btz has also been shown to induce time-dependent increased expression of HSP27, HSP70, and HSP90, whereas inducing apoptosis through proteasome inhibition in MM1S cells 39,40 . Upregulation of HSP90 have been observed in hematological malignancies including myeloma which is required for the stability and function of oncoproteins thereby supporting tumor development [41][42][43] . Our western blotting data also confirmed our transcriptomic analysis findings where Ix increased HSP40, HSP70b1, and HSP90α protein expression in Ix-sensitive HMCLs with little or no effect in Ix-resistant HMCLs. HMCLs representing intermediate PI response (IM-9, ARP-1, MOLP-8, MC CAR, U266, and RPMI-8226) showed variable kinetic PI-response levels.
Shaughnessy et al. performed Btz test dosing on a training set of 142 Total Therapy 3 A (MM-TT3A) patients and a validation set of 128 TT3B patients and identified a GEP80-postBz kinetic expression signature that included 80 highly survival-discriminatory genes 44 . When our kinetic Ix-response profiles in HMCLs and patients were compared with Shaughnessy's GEP80-postBz, 10 genes (primarily involved in PUP) were found similar between all the subgroups (Fig. S5a), including PSMD4-the novel high-risk gene in myeloma patients 44 . Further, comparison of GEP80-postBz with our Ixsensitive-only in vitro and ex vivo kinetic GEP signatures revealed 16 and 12 common (PUP) genes, respectively (Fig. S5b).
At last, for the most Ix-resistant HMCLs (LP1 and UTMC2), we added an additional Ix dose of 150 nM (10 times the test-dose of 15 nM) and performed gene expression analysis. We were particularly interested in DE genes that were consistently upregulated (fold-change > 1) or downregulated (fold-change < 1) among the three groups (no treatment vs 15 nM vs 150 nM).
Our results showed 331 genes were consistently downregulated (p < 0.05) between 0 vs 15 nM and 15 nM vs 150 nM treatments. In all, 493 genes that were consistently upregulated between 0 vs 15 nM (upregulated) and 15 nM vs 150 nM (upregulated) treatments. Thirty-two of these genes belong to the PUP (P = 1.12E-15).
Our results thus provide key insights into the difference between sensitive and resistant myeloma in terms of kinetic PI response. At present, we are using CRISPRbased gene knockdown of Hsp genes to explore novel mechanisms of PI resistance. This will eventually lead to future studies on the understanding of the biology governing the pharmacogenomics of differential PI response and the development of novel therapeutic approaches to circumvent the challenge.