Transcriptional profiles define drug refractory disease in myeloma.

Identifying biomarkers associated with disease progression and drug resistance are important for personalized care. We investigated the expression of 121 curated genes, related to immunomodulatory drugs (IMiDs) and proteasome inhibitors (PIs) responsiveness. We analyzed 28 human multiple myeloma (MM) cell lines with known drug sensitivities and 130 primary MM patient samples collected at different disease stages, including newly diagnosed (ND), on therapy (OT), and relapsed and refractory (RR, collected within 12 months before the patients' death) timepoints. Our findings led to the identification of a subset of genes linked to clinical drug resistance, poor survival, and disease progression following combination treatment containing IMIDs and/or PIs. Finally, we built a seven-gene model (MM-IMiD and PI sensitivity-7 genes [IP-7]) using digital gene expression profiling data that significantly separates ND patients from IMiD- and PI-refractory RR patients. Using this model, we retrospectively analyzed RNA sequcencing (RNAseq) data from the Mulltiple Myeloma Research Foundation (MMRF) CoMMpass (n = 578) and Mayo Clinic MM patient registry (n = 487) to divide patients into probabilities of responder and nonresponder, which subsequently correlated with overall survival, disease stage, and number of prior treatments. Our findings suggest that this model may be useful in predicting acquired resistance to treatments containing IMiDs and/or PIs.

INTRODUCTION

The introduction of immunomodulatory drugs (IMiDs) and proteasome inhibitors (PIs) changed the therapeutic paradigm for treatment of multiple myeloma (MM) due to their specific antimyeloma mechanisms. While the majority of patients receiving combination chemotherapy including IMiD and/or PI initially respond, most of them eventually develop resistance. Understanding the underlying mechanisms of nonresponsiveness and identifying biomarkers associated with drug resistance and disease progression have become critical for personalized medicine and development of novel therapeutic strategies.

IMiDs mediate anti‐MM effects by binding to the E3 ubiquitin ligase cereblon (CRBN) [1, 2, 3], which subsequently increases degradation of the transcription factors Ikaros (IKZF1) and Aiolos (IKZF3), culminating in downregulation of IRF4 and MYC expression leading to inhibition of MM cell growth [4, 5]. IMiD resistance in MM has been linked to deletion, functional mutation, or dysregulation of CRBN and the proteins directly and indirectly associated with CRBN or IMiD‐mediated signaling [2, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Resistance to PIs in MM has also been extensively studied [15, 16, 17, 18, 19, 20, 21] and attributed to mutation and dysregulation of proteasome subunits [15] .

In clinical practice, IMiDs and PIs are usually given in combination with each other as well other classes of drugs (such as dexamethasone) as standard of care. Developing a method to monitor both IMiD and PI drug sensitivity and detect disease progression during treatment is important for precision medicine. In the present study, we sought to identify transcriptional changes associated with treatment nonresponse by studying the candidate genes previously associated with IMID and/or PI sensitivity to determine whether measuring expression levels of such genes could serve as a biomarker for disease progression.

NanoString nCounter technology, a direct multiplexed measurement of gene expression based on digital color‐barcoding technology [22], is a flexible, reproducible, and robust method when used for molecular subtyping of diffuse large B‐cell lymphoma [23, 24, 25]. In this study, we employed this technology to investigate the transcriptional expression of 121 genes potentially associated with IMiD or PI response, in both human myeloma cell lines (HMCLs) and primary MM patient samples collected at different times of the disease.

MATERIAL AND METHODS

Study design

We investigated the transcriptional expression of 121 genes selected from published literature [2, 4, 7, 14, 26–33] and an in‐house database linked to IMiD and PI response and resistance (Figure 1A). We screened 28 HMCLs with known drug sensitivities to IMiDs and PIs and 156 primary MM patient samples, including 41 patients with paired samples, collected at various stages of disease evolution, including newly diagnosed and untreated (ND), collected on therapy (OT), and late relapsed and refractory (RR, bone marrow samples from treated patients taken within the 12 months preceding their death) (Figure 1B, Supporting Information Data 1). Differential gene expression between sample groups with expected distinct drug response profiles and disease stages (sensitive vs. resistant HMCLs and ND vs. during treatment or vs. RR patient samples) were analyzed.

FIGURE 1

Collated gene list and primary multiple myeloma (MM) samples or human multiple myeloma (MM) cell lines (HMCLs) for NanoString profiling. (A) Genes comprising the CodeSet were selected based on previous studies. (B) Patient materials were selected and grouped based on the stage of disease activity when samples were collected. Numbers in brackets indicate the number of probes for each gene or number of patients in each group

MM cell lines and human MM cells

All HMCLs used in this study were provided by the Bergsagel laboratory and fingerprinted to confirm their identity [34]. The cells were cultured in RPMI1640 medium with 5% fetal calf serum. Isogenic IMID and PI‐sensitive and ‐resistant cell lines were previously generated in our laboratory [7, 33]. The generation of OCIMY5/Vec and OCIMY5/CRBN by coculture with drug was also previously described [28].

Primary human MM cells were recovered from bone marrow aspirates collected from all Mayo Clinic sites. Informed consent was given in writing for collection and research use under Institutional Review Board approval (IRBs 919‐04, 15–009436, 18–003198, 2207‐02) in accordance with the Declaration of Helsinki. After collection, CD138+ cells were isolated by immunomagnetic bead selection (RoboSep; Stemcell Technologies).

Sample preparation and RNA extraction

Total RNA from HMCLs and CD138‐selected plasma cells from MM patient bone marrows were isolated using RNeasy Mini kit and the AllPrep DNA/RNA Kit (Qiagen), respectively. After spectrophotometric quantification (NanoDrop2000, Thermo‐Fisher Scientific), RNA was stored at −80°C until use.

Collated gene list

A unique list including 121 gene candidates was selected for NanoString profiling (Figure 1A). We first selected 26 probes targeting CRBN, genes with altered transcription in cell lines with low versus normal CRBN expression, genes encoding proteins associated with CRBN activity [26], or additional genes linked to IMiDs activity and sensitivity (such as IKZF1, IKZF3, IRF4, and MYC). Since CRBN isoforms, including the isoform lacking exon 10, have been associated with IMiD sensitivity [35], four probes targeting different exon junctions of CRBN were included, with one (CRBN 3) designed to span the exon 10/11 junction. Twenty‐one genes were selected by analyzing baseline gene expression levels associated with drug response in a cohort of 44 refractory MM patients before initiation of pomalidomide and dexamethasone therapy on a phase 2 clinical trial [28, 36, 37] (Supporting Information Data 2) and from the isogenic lenalidomide‐sensitive/resistant HMCL XG1 pair (XG1/XG1LenRes) with normal CRBN levels [7] (Geo123506). Twenty of these gene targets were identified between the responders and nonresponders in both data sets.

An additional 10 genes included in the exploratory panel were selected based on data from 59 MM patients that exhibited differential responses to a first line of treatment containing IMiDs in the CoMMpass data (Supporting Information Data 3, generated as part of the Multiple Myeloma Research Foundation Personalized Medicine Initiatives). These 10 genes were differentially expressed in samples that showed complete response, partial response or stable disease; five of which were also identified in the XG1/XG1LenRes analysis. Additional previously reported 24 genes noted as predictive markers of response to IMiDs and/or PIs in MM were included [30, 31, 32]. Since a recent study in our laboratory identified an upregulation of proteasome subunit genes in PI‐resistant cell lines when compared with their isogenic‐sensitive cell lines [33], 40 probes for genes encoding the multiple proteasome subunits were included.

NanoString CodeSet design and expression quantification

The 121informative genes potentially associated with IMiD and PI response, along with 11 housekeeping genes (Supporting Information Data 4), were combined to generate the exploratory CodeSet for this study. The experiments were performed with nCounter Elements XT reagents in accordance with the manufacture's recommendations. An input of 100 ng total RNA was used for all lines and samples. The collected data were first evaluated for quality control, followed by technical normalization using synthetic controls and biological normalization via housekeeping genes. Data were initially analyzed using nSolver 4.0 software and an advanced analysis software plugin (version 2.0, R‐based statistical tool) to detect and visualize differentially expressed genes.

Ranking predictive probes and development of models for predicting drug sensitivity and disease progression

Differentially expressed genes between ND (n = 51) and RR samples (n = 57) were selected using edgeR [38]. The expression correlations between the differentially expressed genes were shown by Pvclust [39]. Each gene was then individually analyzed using a generalized linear model [40] to filter the top q‐associated probes with different response outcomes. To build a multivariate ordinal model for prediction, with the 121 gene probes and annotated ND or RR for each patient sample, a linear logistic regression model was built by using R package bhGLM [40], followed by step Akaike information criterion (AIC) [41] for optimization. The resulting prediction model was based on gene expression from seven genes using following calculation:

Z = −32.1129 + 0.3773DIRAS1 + 2.2558CRBN_2 + 0.8810CD53 + 1.6225PSMD14 – 0.6304CEP55 + 0.9362SK2 ‐ 1.7078PSMA7

Based on the value Z, the probability of “responder” is calculated by sigmoid function

The predicted probability from the model ranges from 0 to 1, where a higher value meant a higher chance of being a responder.

Assessment of the performance of established model

We first evaluated the performance of seven‐gene model (named MM‐IP‐7, stands for multiple myeloma‐IMiD and PI sensitivity ‐7 genes) by five‐fold cross‐validation as described (https://rdrr.io/cran/cvAUC/man/ci.cvAUC.html). Using this model, we retrospectively analyzed existing RNA sequencing (RNAseq) data from the primary samples of the Multiple Myeloma Research Foundation (MMRF) coMMpass (n = 578, ND patients) and Mayo Clinic MM registry (n = 487, collected from patients at different disease stages) datasets to correlate the MM‐IP‐7 results with other clinical data such as survival, disease stage, and number of treatment protocols. Briefly, we submitted RNAseq data to MM‐IP‐7 to calculate probabilities (by ranking scores) and then compared estimated results with other clinical data in each dataset. Since the RNAseq data are highly correlated but have different scales when compared to NanoString, the probability of estimate from this analysis is based on rank order rather than actual cut point criteria.

RESULTS

Validation and quality control of NanoString expression profiling in MM

We first confirmed performance by the nCounter instrument by profiling of MM cells using testing a CodeSet of 43 of the 121 selected genes (Supporting Information Data 5). We demonstrated that nCounter technology was able to generate reproducible results from two biological repeats (MM1.S, Figure 2A). As expected, the NanoString platform also detected CRBN downregulation and IL6 upregulation in two established lenalidomide‐resistant HMCLs, when compared with their isogenic‐sensitive cell lines (Figure 2B,C), consistent with previous observations [7]. A gene expression heatmap of normalized data from four pairs of lenalidomide isogenic HMCLs showed that each isogenic cell line pair clustered together as expected. Further analysis of the expression data using the nSolver 4.0 software (NanoString, Seattle, WA) identified downregulation of CRBN as a significant change in three resistant cell lines (Figure 2D, Figure S1), consistent with previously published data [7].

FIGURE 2

Demonstrating NanoString technology as a sensitive, reliable and reproducible method to quantitate gene expression changes in myeloma cells. (A) Correlation of two biological repeats generated from the NanoString profiling of multiple myeloma1 (MM1).S cell lines. (B and C) NanoString profiling detected the downregulation of CRBN mRNA and upregulation of IL6 mRNA in two different lenalidomide isogenic‐resistant cell lines, consistent with the previous RNA sequencing (RNAseq) data. (D) Heatmap view of the normalized data from four pairs of isogenic introduction of immunomodulatory drugs (IMiDs)‐sensitive/resistant cell lines. (E) Detection of lenalidomide‐mediated transcriptional response in lenalidomide‐sensitive cell line, OCIMY5/Cereblon (CRBN)

Then, using the complete CodeSet, we measured the relative expression levels of the 121 selected genes in untreated and lenolidomide‐treated OCIMY5/vec, which expresses a very small amount of CRBN and is resistant to IMiDs and OCIMY5/CRBN, which has forced CRBN expression and is sensitive to IMiDs. We were able to accurately detect known lenalidomide‐regulated gene expression changes (Figure 2E) in OCIMY5/CRBN. Therefore, NanoString technology was confirmed as a sensitive and reproducible method to quantitate gene expression changes in MM cells.

Identification of differentially expressed genes in MM patients with samples collected at different time

We next measured the differential expression of the 121 genes in all primary MM samples and HMLCs, grouped by known or likely drug sensitivity and resistance profiles. Forty‐three genes were identified that had significantly different expression (p ≤ 0.05) between 52 ND and 57 RR samples (Supporting Information Data 6). In addition to the expected CRBN, we identified six genes (TMEM107, DIRAS1, CD53, TNFRSF13C, LTBP1, and FOS) as most significantly downregulated in RR samples; meanwhile another seven genes (PRR11, CEP55, BIRC5, KPNA2, DEPDC1, PSMB4, and ETV4) were shown as most significantly upregulated (Figure 3A, Supporting Information Data 6). We noticed that all probes against the different CRBN isoforms detected CRBN downregulation in RR samples, but no upregulation of the isoform with exon 10 deletion was identified.

FIGURE 3

Detection of the differentially expressed genes between newly diagnosed and late stage, relapsed/refractory samples. Volcano plot displaying each gene's ‐log10 (p‐value) and log2 fold change with the selected covariate. Highly statistically significant genes fall at the top of the plot above the horizontal lines, and highly differentially expressed genes fall to either side. Horizontal lines indicate various p‐value thresholds. The 20 most statistically significant genes are labeled in the plot. Top 16 differentially expressed genes are shown in the table beside each plot. (A) Fifty late/ relapsed and refractory (RR) samples (bone marrow samples taken from treated patients within the 12 months preceding their death) were compared with 52 newly diagnosed samples. (B) Eleven paired late/RR samples and newly diagnosed (ND) samples were compared

We then analyzed 22 paired samples from 11 patients, which were collected at both ND and RR stages, identifying 45 genes as differentially expressed between ND and RR stages (p ≤ 0.05, Supporting Information Data 7). In addition to confirming most changes above, the RR samples also showed downregulation of IFITM1 and PSMC4 and upregulation of ITPRIPL2, PBK, PSMD4, and CTAG1B as their most differentially expressed genes (Figure 3B, Supporting Information Data 7). These changes were not detected, or detected at a lower significance, when comparing paired samples at ND with a secondary sample collected during treatment but before disease progression; this dataset included samples from patients treated with IMiDs (eight pairs, Figure 4A and Supporting Information Data 8), IMiDs + PI (14 pairs, Figure 4B and Supporting Information Data 9), and solely PIs (three pairs, Figure S2). When comparing ND samples with paired “on active treatment” samples, downregulation of CRBN, CD53, IFITM1 and upregulation of PRR11, CEP55, and BIRC5 was demonstrated (Figure 4B). A similar trend of upregulation of PRR11, ETV4, and BIRC5 was also identified in later relapsed samples compared with early samples collected during treatment from five patients (Figure S3).

FIGURE 4

Detection of differentially expressed genes between newly diagnosed multiple myeloma (MM) and samples harvested during active treatment volcano plot displaying each gene's ‐log10 (p‐value) and log2 fold change with the selected covariate. Highly statistically significant genes fall at the top of the plot above the horizontal lines, and highly differentially expressed genes fall to either side. Horizontal lines indicate various p‐value thresholds. The 20 most statistically significant genes are labeled in the plot. Top 16 differentially expressed genes are shown in the table. (A) Eight paired samples harvested at the time of diagnosis and during or after treatment with introduction of immunomodulatory drugs (IMiDs)‐based therapy (no proteasome inhibitors [PIs] were used) were compared. (B) Fourteen paired samples harvested at the time of diagnosis and during or after treatment with IMiDs and PIs were compared

We next explored HMCLs with known response to IMiDs and PIs [42]. By comparing gene expression of six IMiD‐sensitive and nine IMiD‐resistant HMCLs, 22 genes were identified as differentially expressed (Supporting Information Data 10, Figure S4A). Six changes in IMiD‐resistant HMCLs were consistent with those identified in RR samples from MM patients, including upregulation of PRR11, HN1, RFC3, PSMB2, PSMD14, and downregulation of SKA2. When five paired PI‐sensitive and ‐resistant isogenic cells lines were compared, we identified changes in the expression of seven proteasome subunit genes, including upregulation of PSMB5 in resistant cell lines (Figure S4B), consistent with our previous observations [33].

Identifying predictive probes and establishing a predictive model

Using the NanoString profiling data obtained from all ND (n = 52) and RR patients (n = 57), we further evaluated the predictive value of each differentially expressed gene for correlation with drug resistance and disease progression. Using edgeR software, we identified 45/121 differentially expressed genes between ND and RR samples (p ≤ 0.01) (Supporting Information Data 11). The correlation between those differentially expressed genes was identified by clustering analysis, for example, the expression of PRR11 was found to cluster together with the expression of BIRC5, CEP55, PBK, and DEPDC1 (Figure 5A). We identified 31/45 genes as significant predictors for separating ND from RR patients (p < 0.05, Figure 5B), with CRBN, PRR11, CD53, BIRC5, DIRAS1, DEPDC1, and CEP55 being the most differentially expressed. Finally, using R‐package BhGLM, we built a multivariable ordinal model (MM‐IP‐7) that contained seven‐associated predictors, CRBN, CEP55, DIRAS1, SKA2, CD53, PSMA7, and PSMD14 (Figure 6A). The performance of MM‐IP‐7 was evaluated by five‐fold cross‐validation resulting in an area under curve (AUC) = 0.91 (Figure 6B). Using MM‐IP‐7, we also analyzed RNAseq data from MMRF CoMMpass (n = 578) and Mayo Clinic MM patient registry (n = 487), found that model prediction correlated with OS (CoMMpass data, Figure 6C), disease stage, and treatment (Mayo Clinic MM patient data, Figure 6D,E). As expected, patient samples that classified as “responders” through the MM‐IP‐7 were enriched for longer OS, ND samples, samples without treatment, and fewer prior therapies.

FIGURE 5

Hierarchical clustering of 45 differentially expressed genes between newly diagnosed (ND) and late/ relapsed and refractory (RR) samples and identification of predictive probes. (A) The expression pattern of 45 differentially expressed genes between the ND and late/RR samples (p ≤ 0.01) were analyzed by pvclust. Values at branches are approximately unbiased (AU) p‐values (red) and bootstrap probability (BP) values (green). Clusters with AU ≥ 90 are indicated by the rectangles. (B) Predictive genes were identified by analysis of 45 differentially expressed genes between the ND and late/RR samples (p ≤ 0.01) using single gene GLM model regression with coefficient p‐value ≤ 0.05

FIGURE 6

Establishing the predictive model based on the differentiated expressed genes between newly diagnosed (ND) and Late/relapsed and refractory (RR) samples. (A) A seven‐gene predictive model (multiple myeloma [MM]‐IMiD and PI sensitivity‐7 genes [IP]‐7) was built based on a linear logistic regression with R package BhGLM. (B) Area under curve (AUC) plot with 95% confidence interval resulted from five‐fold cross‐validation of established model. (C) The established model was employed on RNA sequencing (RNAseq) data from CoMMpass dataset for responder/nonresponder prediction. The scores based on this seven‐gene expression in each sample were calculated and ranked. The survival data from 20% samples that ranked at each side of probability of response were compared; the samples on the “nonresponder” probability side have a shorter survival compared with the samples on the “responder” probability side. (D and E) The established model was also employed on mRNAseq data from the Mayo Clinic MM primary patient dataset. The scores were calculated in the samples that grouped by different stage and treatment protocols. (D) Analysis demonstrated that newly diagnosed (ND) patients’ samples more frequently have “responder” probabilities as compared to samples taken during therapy (other) or at refractory and end stages (ES). (E) Compared the patients with treatments (1 or 2 or >3 prior treatment protocols), the patients with no treatment or less treatment have more “responders” probabilities

DISCUSSION

We investigated the expression of 121 selected candidate genes that were previously demonstrated to associate with IMID and/or PI response to determine whether such genes could be combined into a biomarker assay of disease progression following combination chemotherapy regimens. We identified a subset of genes with changes in transcriptional expression during treatment and disease progression. The top downregulated gene in refractory and late stage disease was CRBN. We not only identified downregulation of CRBN in samples harvested during or after treatment with an IMiD‐based therapy, but also observed that downregulation became increasingly prominent in the relapsed and end stages samples. These results are consistent with both our own and other's prior reports, which demonstrated that with progressive IMiDs exposure, CRBN copy number and expression was reduced in MM cell lines and patient samples [7, 43]. Although several CRBN probes were designed to detect various isoforms of CRBN including the one with exon 10 deletion (associated with IMID resistance), no increased expression of this CRBN isoform was identified in resistant samples.

In addition to CRBN, in RR samples we also detected non‐CRBN‐related transcriptional changes identified in patient samples classified as nonresponders to pomalidomide (Supplemental Data 2) and in the XG1LenRes IMiD‐resistant isogenic HMCL, including downregulation of IFITM1 and upregulation of PRR11, BIRC5, PBK, DEPDC1, RNFT2, and CEP55. IFITM1 is one of the interferon‐stimulated genes upregulated by IMiD treatment and suggested to be involved in IMID‐mediated antimyeloma activity [44, 45]. Upregulation of BIRC5, PRR11, PBK, DEPDC1, RNFT2, and CEP55 appeared subsequent to drug exposure in our study. BIRC5 and CEP55 are genes incorporated within the gene expression‐based proliferation indices [46], which is associated with poor prognosis. In the MMRF CoMMpass dataset, we found that high expression of PRR11, PBK, and DEPDC1, like BIRC5, correlated with short survival and poor drug response (Figure S5), implying upregulation of this group of genes may associate with disease progression and reduced drug response. Indeed, we demonstrated that inhibition of PBK reduced myeloma cell growth and enhanced lenalidomide and bortezomib activity in the IMiD‐resistant cell lines JJN3 and XG1LenRes (Figure S6). This result agrees with a recent study, demonstrating that inhibition of PBK overcame lenalidomide resistance [47]. Other upregulated genes identified in RR samples include ETV4. Upregulation of ETV4 has been demonstrated in MM patient samples at the time of acquire IMiD resistance [14].

Based on the NanoString profiling data from ND and RR samples, we identified a subset of genes (Figure 5B) whose transcriptional expression may serve as predictive biomarkers for treatment nonresponse and disease progression, including several known genes involved in IMiD and PI sensitivity, including CRBN [11], CD53 [32], and PSMD4 [48]. We also built a multivariate ordinal model, MM‐IP‐7, for predicting drug resistance and disease progression. We demonstrated that classification into responder or nonresponder from this model correlated well with treatment and disease stage from a clinical dataset encompassing over 400 patients, suggesting that this model may have potential to monitor the emergence of drug resistance and disease progression. Interestingly, we also found that the MM‐IP‐7 prediction, in ND patients, is associated with OS, suggesting it may be used to identify patients that will be less responsive with faster progression. We recognize that a potential limitation of our study is that end stage patients were treated with combinations of multiple drugs, rather than a uniformly treated cohort; therefore some of our observed gene expression changes may be also linked to multiple‐drug resistance during disease progression. It will be important to extend our study to test more samples harvested at early relapse, especially from uniform therapy.

In the current study, we also demonstrated that NanoString nCounter technology is a reliable method for accurately detecting gene expression in MM samples. As described in non‐Hodgkin's lymphoma [23–25, 49], this method requires a small amount of RNA, is reasonably fast, and adaptable to a clinical diagnostic laboratory suggesting its potential application in MM bio‐marker development.

In summary, we developed a NanoString gene expression profiling‐based prediction model, which may serve as a useful tool for clinical investigation and therapy selection. This is the first study to build a multivariate ordinal model to predict clinical drug resistance, poor survival, and disease progression following treatment with IMIDs and PI, by using a comprehensive gene list generated from multiple other publications and the robust NanoString technology, which is currently in clinical practice.

CONFLICT OF INTEREST

AKS has a consulting role for Amgen, Bristol‐Myers Squibb, Glaxo Smith Klein, Janssen, Ono, Regeneron, Skyline Diagnostics, and Tempus Inc. and has received research funding from Amgen, Celgene, and Janssen. The remaining authors declare no competing financial interests.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.