A Gene Expression–Based Predictor for Myeloma Patients at High Risk of Developing Bone Disease on Bisphosphonate Treatment

Multiple myeloma is a plasma cell malignancy characterized by the frequent development of osteolytic bone lesions, which can lead to subsequent skeletal-related events (SRE). SREs, defined as pathologic fracture, spinal cord compression, radiation to bone and skeletal surgery, are common in myeloma with up to 51% of patients developing an SRE within 21 months (1). In addition to the impact on quality of life, myeloma bone disease is also associated with impaired survival. A few studies have investigated correlations of a range of parameters with bone disease development and have found associations with disease stages, as well as molecular subgroups of myeloma (2–4). However, none of these studies have been able to predict patients at risk of developing bone disease while on treatment.

The basic mechanism underlying myeloma bone disease is an uncoupling of bone resorption from bone formation as a consequence of increased osteoclastic activity and inhibition of osteoblast function (5). Several osteoclast activating factors, including macrophage inflammatory protein 1 alpha and receptor activator of NF-κB ligand (RANKL), have been found to be important in regulating bone resorption. RANKL is modulated by a decoy receptor, osteoprotegerin (OPG). Myeloma cells induce RANKL expression in bone marrow stromal cells and, meanwhile, inhibit stromal derived OPG, favoring bone resorption (6). More recently, global gene expression profiling has pointed to a central role for another important pathway, the WNT pathway, in pathogenesis of bone disease. The WNT signaling antagonist DKK1 is an inhibitor of osteoblast differentiation (7) and also increases bone resorption via an increase in RANKL/OPG ratio (8). The sFRP family proteins are also soluble WNT inhibitors, sFRP-2 and sFRP-3/FRZB have been investigated as possible mediators of osteoblast inhibition in myeloma (9–10). The dysregulated microenvironment as a consequence of bone disruption leads to the upregulation of growth factors for myeloma clone such as insulin-like growth factor 1 (IGF-1), creating a vicious cycle of bone destruction and increased myeloma growth (6, 11, 12).

Bisphosphonates are the current standard drugs for the prevention and treatment of myeloma bone disease (13). These drugs inhibit osteolysis by suppressing osteoclastic activity (14). But even with the most effective treatment, zoledronic acid, 30% of patients still go on to develop an SRE by 40 months (15). Thus, there is an unmet medical need in these patients and the challenge is to effectively identify these patients at presentation, so that their treatment can be modified. In this work, we carried out global gene expression profiling analyses on presenting samples from a randomized clinical trial of multiple myeloma to investigate the molecular basis of bone disease and to develop a predictor able to identify patients at high risk of developing SREs while receiving bisphosphonate treatment.

Bone marrow aspirates were obtained from newly diagnosed myeloma patients in MRC Myeloma IX study during standard diagnostic procedures following informed consent. The study was approved by the MRC Leukaemia Data Monitoring and Ethics committee (MREC 02/8/95, ISRCTN68454111). The design, patient evaluation, and endpoints of this trial have been reported previously (15), but in summary, the trial recruited 1,960 patients and prospectively recorded bone disease at presentation as well as SRE occurrence until relapse. Presenting bone disease was defined as presence of any lytic lesion and/or any pathologic fracture by X-ray. The following events were defined as SREs in this study: pathologic fracture, spinal cord compression, radiation to bone, and skeletal surgery. Patients were randomized into 2 bisphosphonate arms, either sodium clodronate or zoledronic acid, to compare their efficacy on reducing SRE occurrence.

The patients included in the gene expression profiling analysis (n = 261) were representative of patients entered into the overall trial (Supplementary Table SA1). The samples were split into a training set (n = 205) and a test set (n = 56). The training set and test set are comparable in terms of baseline bone disease (70% vs. 65%) and bisphosphonate treatment received (clodronate 45% vs. 52%).

Plasma cells were selected to a purity of more than 90% as assessed by morphology using CD138 microbeads as previously described (16). RNA from the purified cells was extracted, amplified, and hybridized to Affymetrix GeneChip Human Genome U133 Plus 2.0 arrays using methods previously described (16). The microarray data have been deposited in NIH Gene Expression Omnibus database under the accession number GSE15695.

A total of 261 samples were included in this analysis after quality control of arrays was done using the Bioconductor package affyQCReport. Gene expression signals were quantified using Robust Multi-array Average (RMA) normalization (Bioconductor package justRMA). A threshold of 7.2 (log₂ transformed from 150) was assigned to small values, and we retained the 21,887 probesets with at least 1% of samples above this threshold. The detailed process defining the expression threshold has been described previously (17).

All analyses were done in R 2.10.1 and Bioconductor. Differentially expressed genes between patient groups were selected using significance analysis of microarray (SAM; Bioconductor package samr; ref. 18), with a 1,000-permutation adjustment and 5% false discovery rate, which produces a t-statistic score.

The correlations between bone disease groups and various clinical parameters were investigated using Fisher's exact test for categorical parameters and Mann–Whitney U test for continuous variables. Any parameters statistically associated with SRE occurrence were tested in a multivariate logistic regression model, together with the gene signature to develop an on-treatment SRE predictor. Stepwise model selection (both directions) by Akaike's information criterion was used to optimize the gene signature. A χ² test (R function ANOVA) was used to test the statistical difference between logistic models. Once an optimal model was obtained, the probability of developing an on treatment SRE can be calculated as below:

The predictive power (sensitivity/specificity) of a model was tested using a Receiver Operating Characteristic (ROC) method and corresponding area under curves (AUC) were calculated. The predictive power was also tested by Harrell's concordance index (c-index; ref. 19) in a Cox model when taking time to events into account. AUC and Harrell's c-index can be interpreted as the chance of getting the prediction correct. The base value 0.5 means that the chance of choosing who is going to get an SRE is 50–50, which implies no predictive ability, whereas 1.0 denotes perfect performance.

Cumulative incidence analysis of SREs was carried out using relapse/death as a competing risk event. Gray's test (20)was used to compare the cumulative incidence curves. The difference between Kaplan–Meier curves was tested using Log rank test. The independence of bone disease from other prognostic factors was tested using multivariate Cox regression.

Using the complete trial dataset (n = 1,960), we investigated the correlation of baseline clinical and FISH parameters with presenting bone disease. At diagnosis, 72.1% of the patients presented with bone disease and were significantly associated with younger age (median 65 vs. 67 years, P = 0.003), higher hemoglobin (median 10.7 vs. 10.3 g/dL, P < 0.001), higher calcium (median 2.4 vs. 2.38 mmol/L, P = 0.01), lower β₂-microglobulin (median 4.2 vs. 4.5 mg/L, P = 0.04), and lower creatinine (median 98 vs. 106 μmol/L, P < 0.001). There was also a significantly higher proportion of hyperdiploidy, a lower percentage of t(4;14) and MAF translocations [either t(14;16) or t(14;20)] in the bone disease group (Table 1). The samples included in the expression training set (n = 205) followed the same trends as those in the whole trial, although only age (median 64.5 vs. 69 years, P = 0.02), creatinine (median 99 vs. 111 μmol/L, P = 0.01), hyperdiploidy and t(4;14; Supplementary Table SA2) reached the statistical significance threshold due to smaller sample size.

Patients with presenting bone disease had a significantly shorter overall survival (OS) compared with those without bone disease in the complete trial dataset (P = 0.009, median 45.5 vs. 51.6 months, Fig. 1A); however, there was no statistical difference between the 2 groups for progression-free survival (PFS; data not shown). Further analysis showed that survival from relapse was significantly shorter for those with presenting bone disease (P < 0.001, median 12.2 vs. 23.4 months, Fig. 1B).

Age, hemoglobin, β₂-microglobulin, albumin, creatinine, t(4;14), MAF translocations, 13q deletion, 17p deletion, and 1q gain were shown to significantly affect OS (P < 0.05) in univariate analyses; thus we looked at the impact of bone disease on survival when taking these factors into account. Multivariate analysis confirmed that bone disease was an independent prognostic factor for OS (P = 0.04, HR = 1.3, 95% CI: 1.01–1.68).

By comparing patients with presenting bone disease versus those without in the training set, 50 genes were identified as significantly differentially expressed by SAM at 5% false discovery rate (Supplementary Table SA3). The most significantly differentially expressed gene was a WNT inhibitor DKK1. By using Gene Ontology biological process terms and KEGG pathways as a tool to investigate differentially affected pathways, other associated genes were found mainly involved in growth factor signaling, apoptosis, and regulation of gene transcription.

In the complete trial dataset 28.3% of patients (545/1,924) developed 1 or more SREs before relapse. The shape of the cumulative incidence curve is unusual with a very rapid onset of events within the first 2 months (Fig. 2A). The nature of the first SREs varied over time. Within the first 2 months from presentation, 82% of events were “radiation to bone” or “bone surgery”; beyond this point the SREs were mainly “pathologic fracture” and “spinal cord compression” (Fig. 3), suggesting that the SREs within 2 months are closely associated with bone disease before presentation and, therefore, may not reflect a subgroup more likely to develop SREs after treatment. Consequently, we excluded the events that occurred within the first 2 months, with 324 cases remaining (Fig. 2B). Defining an optimal time point within which to examine SRE occurrence is also important, and in this context, our data shows that 93.2% (302/324) of first SREs occurred within 2 years. Therefore, we used this time point as a cut-off to define a comparator group of patients not developing SREs before relapse (n = 1,304). Seventy-three patients who died within the first 2 years from causes other than myeloma were excluded from analysis. Consistent with the overall trial data, among the 204 patients with SRE data available in the training set, 60 patients (29.4%) developed 1 or more SREs before relapse.

We analyzed the whole trial dataset to identify clinical and FISH abnormalities associated with SRE occurrence during treatment (n = 302 vs. n = 1,304). It was found that higher calcium (median 2.45 vs. 2.4 mmol/L, P = 0.002), lower creatinine (94 vs. 101 μmol/L, P = 0.003) and the presence of bone disease at presentation (86% vs. 65.4%, P <0.001) were significantly associated with SRE occurrence. A higher proportion of hyperdiploidy and lower percentage of t(4;14) were also associated with SRE development (Table 1, Fig. 2C and D). Patients in the training set (n = 204) followed the same trends as those in the overall trial, although only higher calcium (median 2.5 vs. 2.4 mmol/L, P = 0.002), the presence of baseline bone disease (84.2% vs. 66.4%, P = 0.04) and hyperdiploidy (Supplementary Table SA2) reached a statistical significance. Therefore these factors were incorporated into a gene signature to develop an SRE predictor.

We compared the gene expression data of 38 patients who developed an SRE between 2 months and 2 years with those of 134 cases without any SREs before relapse in the training set. Ten patients, who died within the first 2 years from causes other than myeloma, were excluded. SAM identified 14 genes overexpressed in the SRE group (false discovery rate < 0.05; Table 2). Some of these were IFN-induced genes, and the others mainly involved in cell signaling and mitosis.

Each expression values of gene were dichotomized, using the 75th percentile as a threshold between high and low expression. A 7-gene signature was optimized from the 14 genes by stepwise selection in a multivariate logistic model. This approach reduces the number of genes used in the predictor, yet at the same time maintains the predictive power of the model. The derived 7-gene model showed good predictive performance in a ROC analysis (AUC = 0.804). Cumulative incidence analyses confirmed each of these genes could significantly discriminate the risk of developing an on-treatment SRE (P ≤ 0.001; Supplementary Fig. SA1).

We tested the impact of other SRE-associated factors (high calcium, baseline bone disease, and hyperdiploidy) on this model to examine whether an improvement in the predictive power of the 7-gene model could be made. When the 3 factors were added individually, only high calcium (the upper normal limit 2.6 mmol/L used as cut-off) statistically improved its predictive power (P = 0.036). This suggests that hyperdiploidy and presenting bone disease may not be independent of the 7-gene signature. Indeed, higher expression of NDC80, PRKAA1, EPSTI1, and GJB6 were shown to be positively correlated with hyperdiploidy, and LCP2 as well as GJB6 are positively correlated with baseline bone disease (P < 0.05).

Subsequently, we generated an “on-treatment SRE” predictor comprising 7 genes plus calcium, which predicts early SRE occurrence on bisphosphonates before relapse: z = −3.261+ 1.111*calcium (H) + 0.67*NDC80 (H) + 0.613*LCP2 (H) + 0.744*USP10 (H) + 0.433*PRKAA1 (H) + 1.08*PPP4R2 (H) + 1.122*EPSTI1 (H) + 1.112*GJB6 (H).

ROC method showed that the 7-gene plus calcium model had an AUC of 0.813 with an optimal sensitivity and specificity of 87% and 72%, respectively (Fig. 4A). Its predictive capability has also been confirmed when taking time to event occurrence into account, with a Harrell's c-index of 0.805.

To validate the optimal model, it was applied to a further test set of cases (n = 56). The definitions used to define groups for comparison were consistent with that of the training set, and the derived ROC curve had an AUC of 0.955, which effectively validates the model (Fig. 4B). The predictor was also tested in subsets of patients with or without bone disease at presentation, treated with either zoledronic acid or clodronate, and carried out well in each of the settings with AUCs greater than 0.8 (Fig. 4C–F).

We have developed a gene expression–based predictor that is clinically applicable to identify patients at high risk of developing an early SRE before relapse despite being on bisphosphonate treatment. According to this approach, a patient with 4 or more of the 8 risk factors has a risk of developing an SRE of 62%, which comprises 19% of patients; in contrast, a patient having 3 factors or less has a probability of developing an SRE of only 11%, comprising 81% of patients (Supplementary Table SA4). The general utility of the predictor has been shown irrespective of treatment used or the presence of baseline bone disease. Although this study comprises 2 arms either using zoledronic acid or sodium clodronate, these 2 drugs share the same mechanism in reducing SREs by inhibiting osteoclast activity. Therefore, the group of patients identified by this approach reflects those not responding well to bisphosphonate treatment. These patients may need additional or changing to other types of treatment for bone disease. Recently a DKK1 neutralizing antibody (BHQ880, NCT00741377) and an activin A inhibitor (ACE-011) are being or has been tested in phase I/II trials for multiple myeloma patients (21, 22). They have been shown potentially to both inhibit osteoclast function and increase osteoblast differentiation (22, 23). BHQ880 has also been shown to inhibit multiple myeloma cell growth in a number of in vivo multiple myeloma models (23, 24). Therefore patients identified by the SRE predictor would particularly benefit from these novel therapies, which are expected to significantly improve the care of these patients.

Fourteen genes were upregulated in patients having SREs within 2 years before relapse (Table 2). Three of them were IFN-induced genes; the others were mainly involved in cell signaling and mitosis. The high expression of IFN-induced genes has been reported previously to be linked to subgroups of myeloma patients with poor prognosis (25). The identified genes could give insight into the biology underlying the SRE development for further investigation.

When the first 2-month SREs were included, in addition to the 14 genes outlined, a myeloid gene signature was also identified. The biological basis for this myeloid signature is uncertain; however, it has been reported previously in 2 other large gene expression datasets (3, 26). There is evidence that certain subsets of myeloma cells can express myelomonocytic antigen (27–29), and data from a number of studies have shown that malignant plasma cells can transdifferentiate into osteoclast-like cells and directly participate in bone resorption (30–34). This is consistent with other observations that osteoclasts may not be solely responsible for myeloma bone disease as there is a discrepancy between the severity of bone disease and the number of OCs in the osteolytic lacunae (6). Thus our observation may provide another piece of support to this hypothesis and suggest that a subset of myeloma cells, with a myeloid expression pattern, could play a direct role in bone resorption process. However, the direct role of malignant plasma cells in myeloma bone disease needs to be confirmed in additional in vivo studies before this hypothesis can be generally accepted.

Interestingly, some favorable prognostic factors such as younger age, higher hemoglobin, lower β₂-microglobulin, and lower creatinine are positively associated with presenting bone disease. The explanation for this association is uncertain but could be due to a mechanism whereby patients with bone disease presented earlier as a consequence of the symptoms caused by bone disease. There is also a higher proportion of hyperdiploidy, a lower percentage of t(4;14) and MAF translocations in the bone disease group. Although the associations between hyperdiploidy, t(4;14), MAF, and bone disease have been reported previously (35, 36), the mechanisms underlying these associations are unclear. Potentially, it could be that the natural history of myeloma in the two groups (hyperdiploidy and the adverse chr14 translocations) is different; therefore, the time taken to develop clinical disease may have a significant impact on the likelihood of having bone disease at presentation. Interestingly, multiple myeloma cells with hyperdiploidy lesions have been found to be uniquely dependent on their interaction with the bone marrow environment to trigger increased expression of cyclin D1 (37). In contrast, the hyperdiploid subgroup and cases with cyclin D1 overexpression are underrepresented in plasma cell leukemia (38), in which microenvironmental interactions are clearly less important. In addition, we found that DKK1 and FRZB were expressed at higher levels in the hyperdiploid group and lower levels in the t(4;14) group (data not shown), which may contribute to the bone disease. It is interesting that although presenting bone disease is associated with favorable prognostic features, patients with bone disease have shorter OS compared with those without bone disease. Indeed, multivariate analysis confirmed the independent prognostic effect of presenting bone disease. This observation suggests that the presence of bone disease significantly contributes to the impaired outcome of these patients or, alternatively, the biology of myeloma in this group is distinctly different.

DKK1 is the most significantly differentially expressed gene between patients presenting with bone disease and those without. It has been consistently reported to be produced by myeloma cells, and its expression level is correlated with the extent of presenting bone disease, using various detection methods (microarray, ELISA, RT-PCR, Western blot, and immunostaining), in different datasets (7, 39, 40). Evidence derived from several studies shows that DKK1 not only inhibits osteoblast differentiation and function but also increases osteoclast activity via an increase in RANKL/OPG ratio (7, 8, 41, 42). Three other secreted WNT inhibitors sFRP-1-3 were available for analysis in our dataset; however, only sFRP-3/FRZB was expressed at a higher level in patients with bone disease (P < 0.05, although not identified by SAM). This finding is consistent with the molecular classification of a low bone disease group, characterized by downregulation of both DKK1 and FRZB (3, 26). Both of these genes have also been reported to be upregulated in multiple myeloma samples in comparison with MGUS or normal plasma cells (7, 43–45). In keeping with prior reports, our data suggests that both FRZB and DKK1 contribute to myeloma bone disease. We also looked at the association of another myeloma cell-derived osteoclast activating factor macrophage inflammatory protein 1 alpha with bone disease; however, we could not identify a significant correlation.

Pathway analysis of the genes significantly associated with bone disease showed that the majority were involved in growth factor signaling, apoptosis, or regulation of transcription. These observations suggest that the biology of the myeloma cells from patients with bone disease is significantly different from those lacking bone disease. Genes shown to be downregulated include WNK1 and SOCS3, which are negative regulators of growth factors released from bone destruction process such as IGF-1 (6). IGF-1 is the major growth factor for myeloma cells; its serum level and the expression of its receptor IGF-1R in myeloma cells have been shown to be strong indicators of adverse prognosis in myeloma patients (11, 12). WNK1 encodes a cytoplasmic serine–threonine kinase, which negatively regulates IGF-1 signaling via the phosphoinositide 3-kinase pathway (46). SOCS3 is a cytokine-inducible negative regulator of cytokine signaling, which has been shown to inhibit interleukin-6 (IL-6) and the IGF-1–induced signaling pathway (47). The downregulation of these 2 genes in patients with bone disease leads to enhanced IGF-1 and IL-6 signaling, which may explain the adverse outcome of these patients.

In conclusion, in this article, we have developed an SRE predictor in patients on treatment for myeloma who are also receiving bisphosphonates. The signature is biologically relevant and takes account of both clinical and molecular parameters associated with on-treatment SRE development. This predictor can identify individuals, which comprises up to 20% of the myeloma patients who are not responding well to bisphosphonates and have a 62% chance of developing an SRE within 2 years of presentation. Patients identified using such a predictor would be suitable for additional therapeutic interventions such as bone anabolic agents to reduce SRE incidence. Such an approach could also be useful to define an entry criterion for patients inclusion in clinical trials aimed at reducing myeloma bone disease.

A.J. Child has received unrestricted educational research grants for Myeloma IX from Celgene, Pharmion, Novartis, and Ortho Biotech. G.H. Jackson has received honoraria from speakers bureau for Celgene and Ortho Biotech. G.J. Morgan has received honoraria from speakers bureau for Celgene, Novartis, and Johnson and Johnson and is also a consultant and is on advisory board of Celgene, Novartis, and Johnson and Johnson. The other authors disclosed no potential conflicts of interest.

Data and samples were collected as part of the MRC Myeloma IX phase III trial—clinical trial registration (ISRCTN68454111). Overall trial statistics were carried out by the Clinical Trials Research Unit, University of Leeds, Leeds, United Kingdom. Fluorescence in situ hybridization and cytogenetics were carried out by Leukaemia Research Fund UK Myeloma Forum Cytogenetics Group, Wessex Regional Genetics Laboratory, Salisbury, United Kingdom.

The authors thank the colleagues who also contributed to either discussion or data management: Kevin D. Boyd, Christopher Wardell, David C. Johnson, and Nicholas J. Dickens.

The work was supported by Myeloma UK, Cancer Research UK, the Bud Flanagan Leukaemia Fund, and the Biological Research Centre of the National Institute for Health Research at the Royal Marsden Hospital.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.