Global methylation analysis identifies prognostically important epigenetically inactivated tumor suppressor genes in multiple myeloma

It is becoming clearer that multiple myeloma (MM) is not a single disease but rather comprises multiple diseases with differences in outcome that are defined by their genetic makeup. The identification of factors that mediate these differences in disease biology is of central importance if we are to improve treatment outcome further. To date there has been a focus on DNA changes such as translocations, copy number abnormalities, and mutational spectrum, but more recently the importance of epigenetic changes that affect myeloma cell biology has been realized.^1-3 

Focal epigenetic changes characterize interindividual differences in many cancer types, including hematologic malignancies.^4-6  DNA methylation changes at gene promoters and transcription start sites can be associated with heterochromatin formation and long-term gene silencing, a mechanism that is involved in the inactivation of tumor suppressor genes.^6,7  Recently, further complexity in this system has been noted with important methylation changes of nonpromoter and intergenic regions, but the functional implications of these modifications are unknown.^6,8,9 

MM is characterized by pronounced variation in interindividual DNA methylation patterns that exceeds the methylation variability seen in several solid cancers.⁴  But so far, differential methylation associated with patient prognosis has been noted only at a few genes, such as CDKN2A, TGFBR2, or SPARC.^10-14  In previous work, we investigated DNA methylation patterns associated with disease progression and with molecular subtypes of the disease. We found specific, focal hypermethylation changes in clinically aggressive subtypes, such as plasma cell leukemia (PCL) and the prognostically unfavorable group with translocation t(4;14), suggesting that methylation changes can affect disease biology and that further genes remain to be identified.² 

In the present work, we analyzed the association of genomewide differential DNA methylation with prognosis in 159 patients with myeloma that had been treated in the Medical Research Council (MRC) Myeloma IX trial to gain insights into DNA methylation changes relevant for disease biology and progression. The identified epigenetically modified tumor suppressor genes may contribute to the molecular characterization of a tumor in an individualized treatment approach.

The MRC Myeloma IX trial recruited 1970 newly diagnosed patients with symptomatic myeloma. The median follow-up of the trial at the time point of analysis was 5.9 years. The design, patient evaluation, and end points of this trial have been published elsewhere.^15,16  The study was approved by the MRC Leukaemia Data Monitoring and Ethics committee (MREC 02/8/95, ISRCTN68454111). Informed consent was obtained in accordance with the Declaration of Helsinki.

CD138-positive cells from bone marrow aspirates were obtained from newly diagnosed myeloma patients in the MRC Myeloma IX study (n = 161), from individuals with monoclonal gammopathy of undetermined significance (n = 5), and from patients with PCL (n = 31) following standard diagnostic procedures after informed consent was obtained. Plasma cell selection, purification controls and DNA and RNA extraction were performed as previously described.²  Briefly, cells were selected to a purity of >90% by the use of CD138 microbeads (Miltenyi Biotech). DNA and RNA were extracted using commercially available kits (RNA/DNA mini kit or Allprep kit; QIAGEN) according to manufacturer’s instructions. Interphase fluorescence in situ hybridization analysis was performed on purified myeloma cells as previously described.¹⁷ 

DNA methylation was assessed using the Illumina Infinium HumanMethylation27 BeadArray platform as previously described.²  Briefly, 500 ng of DNA was bisulfite-converted using the Zymo EZ DNA methylation kit (Zymo Research) and hybridized to the Infinium arrays following the manufacturer’s protocol. The fluorescence signals generated for unmethylated and methylated cytosine nucleotides by single-nucleotide extension of locus-specific methylation probes was transformed into a β value ratio ranging from 0 to 1.0, equivalent to 0% to 100% methylation, for the 27 578 interrogated CpG residues. Principle component analysis of the array data identified 2 outlier cases that were excluded from further analyses. Methylation data have been deposited in the National Center for Biotechnology Information’s Gene Expression Omnibus (GEO) and are accessible through GEO series accession number GSE GSE21304 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21304).

Gene expression data were generated on the Affymetrix HG-U133 Plus 2.0 platform as previously described.¹⁸  Matched gene expression and DNA methylation data were available for 115 cases. The expression values were Robust Multi-array Average normalized and log2-transformed using the statistical software package R. The associated microarray datasets have been deposited into GEO under accession number GSE15695 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE15695).

DNA from the following human myeloma cell lines (HMCL; n = 11) was used for methylation profiling on the Illumina array platform: NCI-H929, KMS11, KMS12BM, KMS26, LP1, MM1S, MM1R, RPMI-8226, U-266, JIM-3, and OPM2. Cell line identities were confirmed for all cell lines by performing short tandem repeat profiling.

The cell line KMS11 was used for demethylation experiments and treated with 5-aza-2′-deoxycytidine (DAC; Sigma-Aldrich) at a final concentration of 200 nM for 4 days with DAC freshly added to the cell cultures every 24 hours. DNA and RNA were harvested at baseline, after 96 hours, and 17 days after the removal of DAC and extracted using the Allprep kit (QIAGEN). Gene-specific methylation was analyzed by bisulfite pyrosequencing and gene expression by reverse transcription-polymerase chain reaction (RT-PCR).

For RT-PCR, 1 µg of RNA was reverse transcribed using Superscript III (Life Technologies) followed by thermal cycling for 40 cycles (30 second each: 95°C, 58°C, 72°C) using the Platinum Taq DNA polymerase (Life Technologies). A total of 10 µL of PCR products were used for gel electrophoresis and were visualized on a 2% agarose gel. RT-PCR primer sequences are provided in the supplemental Materials and methods.

For bisulfite pyrosequencing, 500 ng of DNA were bisulfite converted using the EZ DNA methylation kit (Zymo Research) per manufacturer’s protocol. Pyrosequencing assays were designed using the PSQ Assay Design software (QIAGEN) for the interrogated sequence to match the exact genomic locations of the Illumina Infinium methylation array probes. Details of pyrosequencing and primer sequences are provided in the supplemental Materials and methods.

Methylation values were processed using GenomeStudio (Illumina). K-means clustering and statistical analyses were performed in R 2.12.2 and Bioconductor, Graphpad Prism Version 5.01, and IBM SPSS Statistics Version 20. The association with overall survival (OS) for low and high methylation groups was estimated for each methylation probe using the Kaplan-Meier methods (log-rank test). Results were corrected for multiple testing using the Benjamini-Hochberg method with a false discovery rate of 5%. The independence of the methylation groups from other established risk factors was tested by multivariate Cox regression analysis.

We analyzed DNA methylation data of 159 myeloma presentation samples using the Infinium HumanMethylation27 array, interrogating 27 578 individual CpG sites across the genome.²  Of these, 16 268 sites were differentially methylated among the 159 MM patients and used for downstream analyses. The majority of methylation probes on the array is located in CpG-dense areas, mapping to promoters, transcription start sites, or the first exon of genic regions. DNA methylation in these genomic regions follows a binary distribution across a population, with the majority of probes showing low methylation and a subset of probes having medium to high methylation values; indeed, we observed a similar binary distribution of methylation values in our dataset. We hypothesized that this binary distribution represents 2 underlying biologic states and aimed to categorize methylation data for each probe into a “low” or “high” methylation group. To take differences in the distribution of methylation values between probes into account, we applied a k-means clustering algorithm to all 16 268 methylation probes, creating an individual cutoff value separating low from high methylation.

We next investigated associations with OS for low and high methylation groups for all 16 268 categorized probes using log-rank testing with the Kaplan-Meier method. The results of our analysis were corrected for multiple testing using the Benjamini-Hochberg algorithm with a false discovery rate of 5%. We focused the analysis on associations with OS because progression-free survival may vary between different treatment regimens in myeloma, and several treatment arms were tested in the Myeloma IX trial. However, we confirmed that significant associations between methylation and OS were consistently reflected by progression-free survival. A total of 207 methylation probes were significantly associated with OS, termed prognostically significant, differentially methylated regions (p-DMRs). These probes mapped to 195 genes, with a maximum of 4 probes mapping to 1 gene. The full list of the identified p-DMRs can be found in supplemental Table 1. Genomic preferential overlap analysis revealed a more-than-expected overlap by chance of the p-DMRs with genes that are targeted by the polycomb repressive group protein EZH2 and marked by the repressive histone H3K27me3 mark in lymphatic B cells as well as in human embryonic stem cells (supplemental Table 2). There was also a significant overlap of p-DMRs with “poised promoters,” a chromatin state defined by the Encyclopedia of DNA Elements project that characterizes promoters that can rapidly and dynamically change their accessibility and, subsequently, gene expression of the associated gene.¹⁹  We also investigated overlap of p-DMR probe locations with transcription factor binding sites and found SP1-, NFKB-, and EGR1-binding sites in the B-lymphocyte cell line GM12878 to overlap more than expected by chance (supplemental Table 3). These transcription factors are known to play important roles in myeloma biology.^20-23 

With the aim of identifying p-DMRs that have potentially important disease-modulating functions (termed pd-DMRs), we combined methylation with Affymetrix gene expression array data to identify genes in which methylation changed gene expression. The majority of p-DMR genes were not expressed or not differentially expressed across the 115 samples in which we had data available. Interestingly, genes that were differentially expressed and methylated showed a nonlinear correlation between methylation and expression that followed a binary “on/off” gene expression pattern. More precisely, low methylation was associated with a “permissive” state in which the gene was either expressed or not expressed. On the contrary, high methylation was associated with silenced gene expression or a “locked” expression state.

We further filtered p-DMR genes with differential methylation and expression for genes with a known function in cancer biology. This way, we identified pd-DMRs mapping to the following genes (in alphabetical order): CD38 (Infinium array probe ID cg10169812); neural cell adhesion molecule 1 (NCAM1, cg13692433); glutathione peroxidase 3 (GPX3, cg07699362); pyruvate dehydrogenase kinase 4 (PDK4, cg04498739); Ras, dexamethasone-induced 1 (RASD1, cg06990298); retinol-binding protein 1 (RBP1, cg23363832); secreted protein, acidic, cysteine-rich (SPARC, cg25913233); and transforming growth factor, β induced (TGFBI, cg21034676).

Kaplan-Meier survival curves for pd-DMRs with median OS and median progression-free survival for all cases as well as methylation/expression correlation scatter plots for 115 cases are shown in Figures 1 and 2 and supplemental Figure 1. In multivariate analyses, each of the identified pd-DMR genes was tested against established myeloma risk markers with respect to OS. The pd-DMRs at GPX3 (P = .021; hazard ratio [HR] 3.1; confidence interval [CI] 1.2-7.9], RBP1 (P = .007; HR 2.2; CI 1.2-4.0), SPARC (P = .001; HR 2.7-5.1), and TGFBI (P = .003; HR 2.6; CI 1.4-4.8) retained independent statistical significance for OS.

We next tested which of these 4 genes were independent from each other in a multivariate model and found that SPARC and TGFBI were independently associated with OS (Table 1). A comparison of clinical, laboratory, and cytogenetic parameters between the low and high methylation groups for these 2 pd-DMR genes identified an association of high TGFBI methylation with t(4;14) samples but no other significant associations (Table 2).

Subgroup analyses demonstrated significant associations with OS for GPX3, RBP1, and SPARC methylation, both for the intensive and the nonintensive treatment arm of the MRC Myeloma IX trial. TGFBI methylation was prognostic in the intensive treatment arm with a trend toward significance in the nonintensive arm (supplemental Figure 2). Survival analyses of pd-DMRs in cytogenetic subgroups were hampered by small group sizes, but a significant association between SPARC methylation and OS within the group of hyperdiploid samples could be demonstrated (supplemental Figure 3). In 35% of all cases (55 patients), at least 1 of the 4 independent pd-DMRs was methylated. Of these, 55% (30 patients), 25% (14 patients), 11% (6 patients), and 9% (5 patients) were hypermethylated for 1, 2, 3, or all 4 pd-DMRs, respectively. Median OS was 57.1, 39.1, and 13.2 months (P < .0001) for cases with no, 1 and 2 or 3 and 4 hypermethylated pd-DMRs, respectively, indicating an additive effect of the pd-DMRs on tumor biology (supplemental Figure 3). Thus, 4 prognostically independent pd-DMRs with a potential direct impact on tumor biology were identified.

We next analyzed whether methylation of pd-DMRs changed with the progression of MM to PCL or HMCL. We found a statistically significant increase in methylation from MM to PCL and from MM to HMCL for all 4 pd-DMRs, namely GPX3 (mean β-value [MBV] MM 0.088; PCL 0.11; HMCL 0.82; P = .0001 for MM vs PCL and MM vs HMCL), RBP1 (MBV MM 0.19; PCL 0.26; HMCL 0.74; P < .005 for both comparisons), SPARC (MBV MM 0.13; PCL 0.18; HMCL 0.37; P < .005), and TGFBI (MBV MM 0.14; PCL 0.23; HMCL 0.68; P < .005) (Figure 3). Methylation for all 4 pd-DMRs was low in monoclonal gammopathy of underdetermined significance (MGUS). Thus, increasing methylation of pd-DMRs is associated with more aggressive phenotypes of malignant plasma cells.

In order to assess whether these epigenetic changes could be reversed therapeutically, we tested whether de-methylation at pd-DMR genes was associated with gene expression changes. The HMCL KMS11 was treated with low doses (200 nM) of DAC.²⁴  Gene expression analysis by RT-PCR showed that expression of GPX3, RBP1, SPARC, and TGFBI was induced after four days of DAC exposure. Interestingly, gene expression was silenced again when cells were cultured for additional 17 days after removal of DAC. This was accompanied by demethylation of the corresponding pd-DMR loci under DAC treatment and remethylation after drug removal (Figure 4).

This study investigated the association between DNA methylation and prognosis in MM using genome-wide DNA methylation profiling of 159 presentation myeloma cases that were treated in the MRC Myeloma IX trial. The analysis identified 207 p-DMRs mapping to 195 genes. Among these, we identified the epigenetically regulated genes (in alphabetical order) CD38, GPX3, NCAM1/CD56, PDK4, RASD1, RBP1, SPARC, and TGFBI with a potential impact on tumor biology and on pathways that are relevant for responses to therapies (named pd-DMRs). Multivariate analysis confirmed that the pd-DMR genes GPX3, RBP1, SPARC and TGFBI were associated with survival independent of established risk factors in myeloma.

The pd-DMR genes RASD1, TGFBI, SPARC, GPX3, and RBP1 are tumor modulating genes that are involved in the response pathways to therapeutic drugs or growth suppressing stimuli from the microenvironment, a fact that may explain the prognostic impact of their differential methylation in myeloma.

RASD1 encodes for a Ras GTPase with tumor suppressor functions that is directly induced by dexamethasone.²⁵  DNA hypermethylation has been linked to dexamethasone resistance in myeloma cell lines due to non-responsiveness of the hypermethylated RASD1 promoter to dexamethasone. Responsiveness and sensitivity to dexamethasone could be increased by pretreating cells with the demethylating agent DAC.²⁶  Thus, DNA methylation-associated formation of heterochromatin with subsequent decreased plasticity and responsiveness to growth-inhibitory stimuli might influence tumor biology and prognosis for RASD1 and other pd-DMRs.⁷ 

High expression of TGFBI increases the sensitivity toward chemotherapy and has been linked to a reduced metastatic potential in solid tumors.^27,28  Hypermethylation of TGFBI has also been associated with reduced response to chemotherapy in ovarian cancer.²⁹  Interestingly, TGFBI has been shown to be downregulated in the transformation from MGUS to MM.³⁰  Our data establish a prognostic role for TGFBI methylation in myeloma through its biology. Patients with hypermethylated TGFBI had a median OS of only 25.7 months compared with 50.9 months in cases with low methylation (HR 2.4). TGFBI hypermethylation was also associated with presence of the t(4;14) translocation. Of note, the prognostic significance of TGFBI methylation was independent of presence of t(4;14) in multivariate analysis, suggesting that TGFBI hypermethylation might contribute to the unfavorable prognosis in this translocation group.

The gene SPARC/Osteonectin encodes a matrix-associated protein that mediates interactions of cells with the extracellular matrix. SPARC function in cancer is context dependent with tumor suppressor and chemosensitizing functions as well as growth-promoting properties having been described.^31-34  Of interest for myeloma, SPARC expression has been linked to the response to lenalidomide in 5q-myelodysplastic syndrome.³⁵  In line with our finding, Heller et al reported DNA methylation of SPARC to be associated with an adverse prognosis in myeloma, a finding that not only underlines the prognostic significance of SPARC methylation but also validates our experimental approach.¹⁰  In contrast to their work, we tested the prognostic impact of SPARC in a multivariate model for OS and found it to be not only independent of established risk markers, but also of the other pd-DMRs (HR 2.6). Moreover, a prognostic impact of SPARC methylation in the subgroup of hyperdiploid patients could be demonstrated.

In prostate cancer, GPX3 has been shown to interact with p53-induced gene 3 and to induce cell death by sensitizing cells to reactive oxygen species, which are induced by several drugs used for therapy in myeloma.³⁶  Underlining its tumor-suppressive role, GPX3 is frequently deleted or silenced by hypermethylation in solid tumors.^37,38 

Our data also indicate a prognostic role of impaired retinoic acid (RA) signaling in myeloma. The gene RBP1/CRBP1 encodes for a protein involved in the intracellular storage and transfer of RA. Treatment with RA inhibits myeloma cell growth, induces apoptosis, and leads to reexpression of immune markers that are lost in the malignant transformation of plasma cells.^39-41  Interestingly, recent findings indicate that RA is produced by the microenvironment in solid tumors, acting as a tumor suppressor.⁴² 

Exposure of the myeloma cell line KMS11 to low doses of the demethylating agent DAC led to demethylation and reexpression of the pd-DMR genes RBP1, SPARC, and TGFBI. Thus, locking of gene expression by hypermethylation may be overcome by a DNA methyltransferase inhibitor treatment. Pretreatment with DAC in a clinical setting could resensitize tumor cells to therapies or microenvironmental stimuli in the context of pd-DMR hypermethylation.

Differential methylation of the genes CD38 and NCAM1/CD56 was associated with prognosis in our dataset in univariate analyses. For CD38, low methylation was associated with an unfavorable prognosis, similar to findings in chronic lymphocytic leukemia.⁴³  However, unlike in chronic lymphocytic leukemia, no association between CD38 surface expression and prognosis has been established in myeloma, making it unclear how CD38 methylation influences myeloma disease biology. NCAM1/CD56 mediates the interaction of tumor cells with their microenvironment. CD56 surface expression is frequently lost in PCL, and absent CD56 surface expression has been associated with shorter survival in myeloma, at least in the context of nonintensive chemotherapy.^44-46  However, no association of CD56 methylation and survival has yet been observed.

A significant proportion of the p-DMR genes mapped to genomic loci that are bound by the polycomb group repressive complex protein EZH2 in embryonic stem cells as well in B-cell lymphocytes. This is in line with a recent report in acute myeloid leukemia that found an association of prognostic DNA methylation patterns with stem cell chromatin marks,⁵  a fact that may be in part explained by the facilitated reprogramming capacity at polycomb marked genes.⁴⁷  We also observed an accumulation of genes that are known to play a role in early hematopoiesis in our set of p-DMRs, such as BAALC, SOX7, and CEBPA. It is unclear if epigenetic reprogramming of these genes may have already occurred at early differentiation stages of the cell that later develops into the myeloma clone, with the effects of the epigenetic memory becoming obvious only at the myeloma stage, as recently postulated for other epigenetic marks.^48,49  Alternatively, rare but highly relevant cell subsets of the tumor may express stem cell genes that are thought to be silenced in myeloma based on bulk expression data, and recent data indicate that investigation of gene expression in cell subpopulations may yield surprising results that cannot be anticipated by bulk expression analysis.⁵⁰ 

In conclusion, the integration of methylation and expression data led to the identification of prognostically relevant pathways that are disturbed by epigenetic reprogramming in myeloma. Hypermethylation at pd-DMR genes may likely affect the responses to antimyeloma therapy and perturb the interactions of the myeloma cell with the microenvironment. Assessment of DNA methylation of these pd-DMRs could provide a powerful tool for risk stratification as well as for individualized treatment approaches in myeloma.

This work as supported by a grant from the Deutsche Forschungsgemeinschaft (DFG KA3338/1-1) (M.F.K.) and a Myeloma UK program grant (D.C.J., B.A.W., C.P.W., and F.M.). F.E.D. is a Cancer Research United Kingdom Senior Clinical Fellow.

Contribution: M.F.K. and D.C.J. designed and carried out research, analyzed and interpreted data, and wrote the manuscript; P.W. contributed vital new analytical tools and analyzed data; B.A.W. designed and carried out research and analyzed data; C.P.W., A.B., F.M., L.M., and F.E.D. analyzed data and provided intellectual input; and G.J.M. obtained funding and wrote the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Gareth J. Morgan, Haemato-Oncology Research Unit, Division of Molecular Pathology, The Institute of Cancer Research, 15 Cotswold Rd, Sutton, Surrey, SM2 5NG, United Kingdom; e-mail: gareth.morgan@icr.ac.uk.