Investigation of a gene signature to predict response to immunomodulatory derivatives for patients with multiple myeloma: an exploratory, retrospective study using microarray datasets from prospective clinical trials

Summary

Background Immunomodulatory derivatives (IMiDs), along with proteasome inhibitors, are key components of treatment regimens for multiple myeloma. Nonetheless, outcomes vary among treated individuals. Drug-specific gene-expression profile (GEP) signatures that aid the prediction of favourable and unfavourable outcomes can provide patients with the most effective therapy for their individual disease. We aimed to develop and validate a gene expression signature to suggest which patients would benefit most from IMiD-based therapies.

Methods For this exploratory retrospective study, we selected a cohort of patients with newly diagnosed or relapsed or refractory multiple myeloma who were treated in clinical trials with IMiD-containing regimens. Cohorts were eligible if they had publicly available GEP data from patients’ bone marrow plasma cells, with long-term follow-up and clinicopathological data. In the development stage of the model, we identified 176 IMiD response genes that were differentially expressed before and after IMiD exposure using pharmacogenomic GEP data from patients who had bone marrow samples taken before and 48 h after a test dose exposure with thalidomide (n=42), lenalidomide (n=18), or pomalidomide (n=18). 14 of these genes had p values less than 0·05 for associations with progression-free survival in patients who received thalidomide in induction and maintenance therapy in the Total Therapy (TT) 2 trial (ie, the training cohort). We combined the 14 genes to create a continuous IMiD-14 score and an optimal cutoff. The subgroup with an IMiD-14 score higher than the cutoff was deemed to be IMiD-resistant. We obtained validation cohorts from four studies of IMiD combination regimens: the TT3a trial (thalidomide in induction and maintenance), the TT3b trial (thalidomide in induction and lenalidomide in maintenance), the TT6 trial (thalidomide in induction and lenalidomide in maintenance), and the vincristine, doxorubicin, and dexamethasone (VAD) group of the HOVON65/GMMG-HD4 trial (thalidomide in maintenance). The primary endpoint was to show the prognostic value of the IMiD-14 gene signature for progression-free survival.

Findings In the training cohort, progression-free survival was significantly shorter in the 83 patients with IMiD-14 high scores than in the 92 patients with IMiD-14 low scores; 3 year progression-free survival was 52% (95% CI 42–64) for the IMiD-14 high group versus 85% (78–92) for the IMiD-14 low group, with a hazard ratio (HR) of 2·51 (95% CI 1·72–3·66; p<0·0001). These findings were supported by the results in the validation cohorts, TT3a (115 patients with IMiD-14 high vs 160 patients with IMiD-14 low; 3 year progression-free survival 63% [95% CI 55–73] vs 87% [82–92]; HR 1·54 [1·11–2·15], p=0·010), TT3b (77 patients with IMiD-14 high vs 89 patients with IMiD-14 low; 62% [52–74] vs 80% [72–89]; HR 2·07 [1·28–3·34], p=0·0024), TT6 (20 patients with IMiD-14 high vs 36 patients with IMiD-14 low; 39% [22–68] vs 74% [61–90]; HR 2·40 [1·09–5·30], p=0·026), and the VAD group of HOVON65/GMMG-HD4 (65 patients with IMiD-14 high vs 77 patients with IMiD-14 low; 16% [9–28] vs 54% [44–67]; HR 2·29 [1·52–3·45], p<0·0001).

Interpretation Our results suggest that the IMiD-14 model has prognostic value in patients with multiple myeloma who are treated with IMiDs. Some genes in the model could provide novel targets for IMiD resistance and therapeutic intervention. The IMiD-14 model warrants evaluation in prospective studies.

Introduction

Immunomodulatory derivatives (IMiDs) have become an integral component of treatment for multiple myeloma. Initially approved in the relapsed or refractory setting, and later introduced into the upfront setting with and without autologous stem-cell transplantation, IMiD-based treatment
combinations have produced unprecedented disease response and progression-free survival outcomes. Despite these therapeutic advances, success is limited by disease relapse and progression that occur in almost all patients because of innate or acquired drug resistance. The precise mechanisms of IMiD resistance are not entirely clear.

Lenalidomide and other IMiDs have been shown to mediate antimyeloma effects through the E3 ubiquitin ligase cereblon, which increases degradation of the transcription factors Ikaros (IKZF1) and Aiolos (IKZF3).1–3 However, attempts to accurately predict resistance based on the RNA or protein concentrations of individual predictive markers such as cereblon, IKZF1, or IKZF3 have had contradictory results.4–G Several gene-expression profile (GEP) signatures have also been developed to stratify patients with multiple myeloma into different risk groups on the basis of GEP samples from different clinical trials.7–10 These prognostic signatures are either focused on bortezomib-response genes or are more general to the multiagent chemotherapies used in the corresponding trials. None have been specific to IMiDs.

The purpose of this study was to develop a GEP signature focused on IMiD response genes and to predict which patients would benefit most from IMiD-based therapies. The primary endpoint was to show a correlation between the IMiD-based signature with progression-free survival.

Methods

Study design and participants

In this retrospective exploratory study, we selected a subset of patients with multiple myeloma (newly diagnosed or relapsed or refractory) from prospective treatment trials that used combination regimens that included one of the three IMiDs (thalidomide, lenalidomide, or pomalidomide) and had GEP microarray data from purified bone marrow plasma cells.11–17 Patients selected for analysis needed to have data for clinicopathological variables and progression-free survival. We used no other inclusion or exclusion criteria.

To select the IMiD response genes for model development, we used GEP-based pharmacogenomics of paired CD138-purified plasma cells, which we obtained from a cohort of patients with multiple myeloma11,12 who had bone marrow samples taken before and 48 h after the administration of test doses of thalidomide (400 mg per day for 2 days; 42 newly diagnosed cases),11 lenalidomide (25 mg or 50 mg per day for 2 days; 18 relapsed or refractory cases),11 or pomalidomide (4 mg per day for 2 days; 18 relapsed or refractory cases).12 We compared the GEP- based pharmacogenomic data from before and after IMiD administration using significance analysis of microarrays. We deemed genes that had significantly changed expression levels after 48 h with all three IMiDs (p<0·05 by paired significance analysis of microarrays and with a change in the same direction) to be IMiD response genes. For the training cohort, we used data from patients in the thalidomide group of the Total Therapy (TT) 2 trial,13 in which thalidomide was given in induction and consolidation with autologous transplants, and maintenance (appendix pp 3, 4 for details of the trial). For validation cohorts, we used four independent multiple myeloma datasets14–17 for which both baseline GEP of purified plasma cells and survival data were available. All were from clinical trials in which patients received IMiD-containing regimens: the TT3a trial (thalidomide in induction, consolidation after tandem autologous transplants, and maintenance),14 the TT3b trial (thalidomide in induction and consolidation with tandem autologous transplants, and lenalidomide in maintenance),15 the TTG trial (thalidomide in induction, lenalidomide in consolidation with tandem autologous transplant, lenalidomide in maintenance),1G and the vincristine, doxorubicin, and dexamethasone (VAD) group of the HOVONG5/GMMG-HD4 trial (who received maintenance with thalidomide).17 More details about the training and validation cohorts are provided in the appendix (pp 3, 4). This study was exempt from institutional review board approval as it involved collection of publicly available data or de-identified information recorded by the investigators.

Procedures

The plasma cell isolation methods used in each cohort are described in the appendix (pp 3, 4). Specimens were obtained at the time of inclusion in the study, before the start of the treatment specified in the respective study protocol. We downloaded the gene expression datasets used in this study from NIH Gene Expression Omnibus (accession number GSE854G, GSE2G58, GSE57317, and GSE19784) and EBI ArrayExpress (under accession number E-MTAB-2441 and E-TABM-1138). We used the MAS5 algorithm to preprocess and normalise the Affymetrix expression data.

We combined baseline GEPs and progression-free survival data from the training cohort to build an IMiD-resistance gene signature. We used IMiD response genes that had a p values less than 0·05 on univariate Cox regression analysis for progression-free survival in the training set to calculate the IMiD-resistance score, which was the difference between the average log2-scale expression of prognosis-unfavourable (hazard ratio [HR] >1) and prognosis-favourable (HR<1) genes. We then established an optimal cutoff for the IMiD-resistance score with the running log-rank test, so that patients with scores higher than the cutoff were deemed to be IMiD resistant.

Outcomes

The primary objective of the study was to determine the discriminatory ability of IMiD GEP model to predict treatment outcome. The primary endpoint was to assess this model’s prognostic ability with respect to progression-free survival and the secondary endpoint was to determine its prognostic ability for overall survival. The secondary objectives were to establish whether the IMiD GEP model contributes additional information about likely clinical outcome beyond the information provided by standard clinical and pathological factors and known IMiD resistance genes.
Since IMiDs and proteasome inhibitors are essential components of many anti-myeloma regimens, we investigated whether combining IMiD GEP model with previously published 80-gene GEP model (GEP80; derived from bortezomib response genes)8 provided better prognostication for patients receiving both classes of drugs.

Statistical analysis

We used the MAS5 algorithm to evaluate differences between risk groups. We estimated progression-free survival and overall survival using Kaplan-Meier techniques and used the log-rank test to evaluate differences between subgroups. We estimated HRs comparing risk subgroups with Cox proportional hazards models. We did univariate and multivariate analyses of covariates and time-to-event outcomes by Cox regression. We evaluated the prognostic performances of the IMiD-14 model, and other IMiD-related bimomarkes, cereblon, Ikaros, and Aiolos using time- dependent receiver operating characteristics (ROC) curves for censored data. We used time-dependent ROC analysis for censored data and the area under the curve (AUC) as an indication of model performance, with an AUC of 0·5 indicating no predictive ability and an AUC of 1 representing perfect predictive ability. We did the statistical analyses with R version 3.1.0 and used the Ingenuity Pathways Analysis version 3GG01845 online tool for pathway analysis.

Role of the funding source

The funders of the study had no role in study design, data collection, data analysis, data interpretation, or writing of the report. QZ, MB, and SZU had access to the data. The corresponding author had full access to all data and had final responsibility for the decision to submit for publication.

Results

By applying paired significance analysis of microarrays and a p value cutoff of 0·05, we identified 17G genes with expression levels that changed significantly 48 h after a test dose of an IMiD relative to baseline (appendix pp 5–9). Among these 17G genes, 14 had p values less than 0·05 in univariate Cox regression analysis of progression-free survival in the thalidomide group of the TT2 trial (table 1). We combined these 14 genes to create an IMiD-based risk score (IMiD-14), which we defined as the average log2-scale differential expression of the four prognosis-unfavourable genes (HR >1) and the 10 prognosis-favourable genes (HR <1) genes. We then established an optimal cutoff of –1·075 for the IMiD-14 score using the running log-rank test, so that a patient would be classified as having IMiD-resistant disease if the score was higher than the cutoff and IMiD-sensitive disease otherwise (appendix p 2).

The training set included 175 patients, of whom 83 had an IMiD-14 high score and 92 had an IMiD-14 low score. We found a significant survival differences between the IMiD-14 high-score and IMiD-14 low-score groups in the training set (the thalidomide group from TT2; figure 1). Both progression-free survival and overall survival were significantly higher in the IMiD-14 low group than in the IMiD-14 high group. More patients in the IMiD-14 low group were progression-free and alive at 3 years than those in the IMiD-14 high group (figure 1).

Next, we examined whether the IMiD-14 model was able to discriminate outcomes in four independent validation cohorts in which patients with multiple myeloma were treated with IMiDs (figure 1). As in the training set, both 3 year progression-free survival and 3 year overall survival were significantly higher in the IMiD-14 low group than in the IMiD-14 high group in patients from TT3a, TT3b, and TTG, and in the VAD group from HOVONG5/GMMG-HD4 (figure 1).

To identify pathways that were affected after IMiD intake, we did pathway analysis on the 17G genes that responded to IMiD administration (appendix pp 5–9). The top five most affected canonical pathways were the nuclear factor of activated T cell (NFAT) pathway in immune response regulation, phospholipase C signalling, interferon signalling, phosphoinositide 3-kinase (PI3K) signalling in B lymphocytes, and integrin signalling (appendix pp 10, 11). Our analysis showed significant enrichment for genes of the NFAT signalling pathway, which are involved in the regulation of the immune response and drug resistance.

We combined the gene signatures from the IMiD-14 and GEP80 models and used them to divide the patients in the TT3a, TT3b, and TTG trials, who all received regimens including both an IMiD and a proteasome inhibitor, into different risk subgroups. We divided patients into four subgroups: IMiD-14 low and GEP80 low, IMiD-14 high and GEP80 low, IMiD-14 low and GEP80 high, and IMiD-14 high and GEP80 high. Patients with the IMiD-14 high and GEP80 high signature had the worst outcomes (figure 2).

No genes were common to both the IMiD-14 model and GEP80 model. Given that three of the four prognosis-unfavourable genes in IMiD-14 model are on chromosome 1q, we explored whether IMiD-14 risk score is associated with 1q amplification. We used PSMD4 expression as a surrogate for 1q amplification on the basis of our previous work showing a high correlation between fluorescent in-situ hybridisation-derived 1q copy number and PSMD4 expression.8,18 We found that GEP80-defined risk, but not IMiD-14-defined risk, was strongly correlated with PSMD4 expression level (appendix p 12). The prognostic value of the IMiD-14 model did not seem to be associated with an increase in 1q copy number-linked PSMD4 expression.

When we combined the training (TT2 thalidomide group) and test (TT3a, and TT3b) datasets with sufficient data to analyse, univariate analysis showed that several standard clinicopathological variables, including the patient’s age and albumin, β2-microglobulin, creatinine, haemoglobin, and lactate dehydrogenase concentrations were associated with progression-free survival and overall survival (table 2). High-risk designation on the IMiD-14 or GEP80 model, GEP-derived TP53 deletion, high PSMD4 expression, International Staging System (ISS) stage III, and the presence of cytogenetic abnormalities also represented adverse features for progression-free survival and overall survival in univariate analysis (table 2). In the multivariate Cox regression analysis, IMiD-14-defined high risk, ISS stage III, high lactate dehydrogenase, cytogenetic abnormalities, and GEP80-defined high risk predicted inferior progression-free survival and overall survival (table 2).

We then compared the prognostic performance of the IMiD-14 model with the IMiD-related biomarkers cereblon, Ikaros, and Aiolos, which have previously been reported to predict IMiD resistance.1–3 The IMiD-14 model showed consistently good prognostic performance over the follow- up duration in both the training set (TT2 thalidomide group) and the four independent test sets (TT3a, TT3b, HOVONG5/GMMG-HD4 VAD group, and TTG; appendix p 13), with an average AUC of 0.G8 (95% CI 0·G7–0·G8) in all cohorts. The IMiD-14 model outperformed stratification based on the expression of cereblon, Ikaros, or Aiolos, which yielded lower AUCs. The exception was Aiolos expression during the first several years of the TT3b follow- up. Additionally, the cereblon and IMiD-14 model had very similar AUC in HOVONG5/GMMG-HD4.

To determine whether the IMiD-14 model has some risk discriminatory power in certain non-IMiD combination therapy datasets, we analysed patients in the TT2 study who did not receive thalidomide. Patients who had a high IMiD-14 score had worse progression- free survival than those with a low score (appendix p 14).

Discussion

We identified a signature of 14 genes from IMiD-specific GEP pharmacogenomic data, which was associated with clinical outcomes after IMiD-based treatment for multiple myeloma. The IMiD-14 model was cross-validated in four independent GEP datasets from clinical trials of treatment regimens that included IMiDs. In each of these studies, patients with high IMiD-14 scores had poor progression-free survival and overall survival compared with patients who had low IMiD-14 scores. Other commonly used molecular signatures such as GEP707 and EMC929 have been used to stratify patients with multiple myeloma into different risk categories, and can identify 10–30% of patients tested as being at high risk of relapse. However, the focus of these signatures differs from the IMiD-14 model in that they are based on outcomes from more than one chemotherapeutic regimen rather than outcomes specific to one class of drugs. The IMiD-14 model, however, was derived from a set of 14 survival-discriminatory genes that were altered specifically after short-term exposure to single agent thalidomide, lenalidomide, or pomalidomide, which suggested that these genes represent specific tumour or tumour microenvironment responses to IMiDs. In both training and testing datasets, patients with high IMiD-14 scores, who were theoretically resistant to IMiDs, had worse progression-free survival and overall survival than did patients with low IMiD-14 scores after IMiD- containing therapies.
An unresolved issue is whether the IMiD-14 gene profile is specific for IMiDs or applicable to different chemotherapy regimens. The model consists of 14 IMiD-responsive genes and the prognostic value of the signature was successfully validated in four independent publicly available multiple myeloma GEP datasets with large cohorts of patients who received IMiD-containing therapies, which suggest that the signature is an IMiD-related prognostic gene expression signature. However, patients in our training cohort received an IMiD in combination with other potent anti-myeloma agents such as dexamethasone, doxorubicin, and high-dose melphalan. Therefore, it is difficult to isolate the individual drug effect with precision. We found that the IMiD-14 model has some risk discriminatory power in certain non- IMiD combination therapy datasets. It is possible that GEP profiles are heterogeneous and reflect the biological complexity of chemotherapy responses, since multiple resistance mechanisms could be common to several drugs. The best strategy to address this possibility would be to comprehensively investigate the predictive value of the IMiD-14 signature in the setting of prospective single- agent IMiD maintenance trials.

In real-world scenarios, patients with multiple myeloma are often treated with therapy consisting of an IMiD or a proteasome inhibitor, or combinations of these drugs when appropriate. We combined the IMiD-14 model with a previously described proteasome inhibitor-related prognostic model derived from bortezomib response genes, the GEP80 signature.8 By combining these two drug-specific gene signatures, we stratified patients into four risk subgroups based on IMiD and proteasome inhibitor resistance and sensitivity. For patients treated with both an IMiD and a proteasome inhibitor, the subgroup who had high scores for both IMiD-14 and GEP80 had the worst progression-free survival and overall survival, whereas the other three subgroups could not be clearly separated across all tested datasets. In the multivariate Cox regression analysis, IMiD-14-defined high risk was an independent adverse predictor for progression-free survival and overall survival, in addition to GEP80-defined high risk and ISS stage III. Hence, the IMiD-14 model would seem to complement the GEP80 model as a valuable prognostic tool in patients with multiple myeloma who are treated with IMiDs and proteasome inhibitors. For patients predicted to be resistant to both drugs, the addition of new classes of drugs or later generation IMiDs or proteasome inhibitors to treatment regimens might be worthwhile.

Furthermore, the IMiD-14 model outperformed the other previously described individual IMiD-resistance biomarkers, cereblon, Ikaros, and Aiolos, in predicting survival outcomes. The canonical pathway analysis of 17G IMiD response genes showed enrichment for the NFAT-regulated immune response pathway and other pathways involved in phosphate C, interferon, PI3K, and integrin signalling. NFAT proteins are transcriptional activators of interleukin 2, a key regulator of T-cell immune response.19 IMiDs cause nuclear translocation of NFAT2 and AP-1 via the activation of PI3K signalling, resulting in interleukin-2 secretion and T-cell proliferation.20 Data from other studies suggest that overexpression of the exportin XPO1/CRM1 can mediate nuclear export of NFAT, resulting in termination of its action.21 In addition to NFAT, XPO1 can mediate the nuclear export of other key cargo proteins, including cyclin B, MAPK, and AP-1.22 Notably, in the IMiD-14 model, XPO1 was the top- ranking gene that was associated with poor prognosis. High expression of XPO1 has also been associated with poor survival in other cancers.23 KPT-330 is an XPO1 antagonist that is in clinical development and is being investigated in multiple clinical trials in patients with relapsed or refractory multiple myeloma. Emerging data suggest that KPT-330 is a promising agent for patients with penta-refractory disease (resistant to daratumumab, carfilzomib, bortezomib, pomalidomide, and lenalidomide)24 and several trials are investigating KPT-330 in combination with existing therapies across the broader population in multiple myeloma (NCT02343042, NCT031105G2, NCT02831G8G, and NCT02780G09).

Other genes in the IMiD-14 model have been implicated in functions involving cell growth, migration, differentiation, and metabolism. One gene worthy of note in the prognosis-favourable group is TNFRSF7. The protein encoded by this gene CD27 is a co- stimulatory receptor that is expressed on the surface of T, B and, natural killer (NK) cells, providing a target to enhance NK-mediated tumour clearance, while generating adaptive immune response by IMiDs. Other genes of interest include LAMA5 and ITGAG, which have roles in cell-matrix interactions and mediate cell adhesion to the endothelium. Treatment with IMiDs has shown to decrease the expression of integrin subunits, integrin receptors, or both.25 Functional studies are needed to investigate the role of specific genes in the IMiD-14 model to improve understanding of the biological basis of IMiD resistance and responsiveness that could be targeted to develop more effective therapeutic combinations.

Our data show that different risk groups can be identified with the IMiD-14 model. These risk groups were associated with significantly different outcomes in patients with multiple myeloma who received polychemotherapy that included IMiDs. The IMiD-14 model adds new knowledge to the field, and we are further validating this gene signature in prospective clinical studies. The addition of the IMiD-14 model to traditionally used clinical and pathological factors could provide valuable information to establish which patients might benefit from IMiDs, so that patients predicted to derive less clinical benefit from IMiDs could be offered alternative treatment regimens including new classes such as monoclonal antibodies, other immunotherapies, and epigenetic therapies. Unlike traditional multiple myeloma gene signatures, which predict the overall effect of a certain combination therapy, drug-specific gene signatures derived from pharmacogenomics studies, such as the IMiD-14 and GEP80 models, are useful complementary tools in the era of personalised medicine.