Whole-epigenome analysis in multiple myeloma reveals DNA hypermethylation of B cell-specific enhancers.

While analyzing the DNA methylome of multiple myeloma (MM), a plasma cell neoplasm, by whole-genome bisulfite sequencing and high-density arrays, we observed a highly heterogeneous pattern globally characterized by regional DNA hypermethylation embedded in extensive hypomethylation. In contrast to the widely reported DNA hypermethylation of promoter-associated CpG islands (CGIs) in cancer, hypermethylated sites in MM, as opposed to normal plasma cells, were located outside CpG islands and were unexpectedly associated with intronic enhancer regions defined in normal B cells and plasma cells. Both RNA-seq and in vitro reporter assays indicated that enhancer hypermethylation is globally associated with down-regulation of its host genes. ChIP-seq and DNase-seq further revealed that DNA hypermethylation in these regions is related to enhancer decommissioning. Hypermethylated enhancer regions overlapped with binding sites of B cell-specific transcription factors (TFs) and the degree of enhancer methylation inversely correlated with expression levels of these TFs in MM. Furthermore, hypermethylated regions in MM were methylated in stem cells and gradually became demethylated during normal B-cell differentiation, suggesting that MM cells either reacquire epigenetic features of undifferentiated cells or maintain an epigenetic signature of a putative myeloma stem cell progenitor. Overall, we have identified DNA hypermethylation of developmentally regulated enhancers as a new type of epigenetic modification associated with the pathogenesis of MM.

Multiple myeloma (MM) is an aggressive and incurable neoplasm characterized by clonal proliferation of plasma cells in the bone marrow and a marked clinico-biological heterogeneity (Morgan et al. 2012; Bergsagel et al. 2013). MM frequently arises from a premalignant state known as monoclonal gammopathy of undetermined significance (MGUS), although the precise molecular mechanisms involved in the progression from MGUS to MM are only partially understood. Various distinct genetic abnormalities have been observed in both MM and MGUS, including gene mutations, chromosomal rearrangements, or aneusomies (Bergsagel and Kuehl 2005; Chng et al. 2006; Chapman et al. 2011).

Besides genetic abnormalities, recent studies have shown that epigenetic changes such as DNA methylation play an important role in MM. Various reports indicate that DNA methylation patterns are capable of distinguishing normal plasma cells (NPCs) from MGUS and MM cells. The major differences between these entities can be explained by the occurrence of DNA hypomethylation in malignant plasma cells (Salhia et al. 2010; Walker et al. 2011; Heuck et al. 2013). In addition to this DNA hypomethylation, some studies have demonstrated aberrant DNA hypermethylation of promoter regions of different tumor suppressor genes in MM, although the classical CpG island methylator phenotype (CIMP) extensively observed in a wide variety of tumors (Issa 2004) has been rarely reported in MM (Martin et al. 2008). Hypermethylation of RASD1, for example, has been correlated with resistance of MM to dexamethasone (Nojima et al. 2009). Inappropriate DNA methylation of TNFRSF18 (also known as GITR) (Liu et al. 2013), MIR34B/C (Wong et al. 2011), or the combined inactivation of genes GPX3, RBP1, SPARC, and TGFBI (Kaiser et al. 2013) have been associated with poor prognosis, survival, and disease progression in patients with MM.

In spite of these significant findings, the high-throughput DNA methylation reports published so far in MM were restricted to the study of promoter regions (Nojima et al. 2009; Salhia et al. 2010; Walker et al. 2011; Heuck et al. 2013; Kaiser et al. 2013). Hence, the purpose of our study was to adopt a more extensive and unbiased analysis of the DNA methylome, including promoters, gene bodies, and intergenic regions in normal plasma cells (NPC) and plasma cells from MGUS and MM patient samples. Using this approach, we have identified that DNA methylation of B cell-specific enhancer regions is a new phenomenon associated with MM pathogenesis.

Results

The DNA methylome of MM is globally characterized by a large degree of heterogeneity

To analyze the DNA methylome and define general epigenetic signatures associated with plasma cell disorders, we applied the HumanMethylation450 BeadChip (Illumina) to purified plasma cells obtained from bone marrow samples of MM (n = 104) and MGUS (n = 16) patients as well as normal bone marrows (n = 3 pools from four donors each) and nontumoral tonsils (n = 8) (Supplemental Table 1). Unsupervised principal component analysis of the normalized DNA methylation data identified a clear distinction between NPC and MM samples, with a larger degree of heterogeneity in the myeloma plasma cells (Fig. 1A,B; Supplemental Fig. 1). Next, we further characterized this heterogeneity of MM by comparing it with other lymphoid malignancies. We calculated the median methylation value per case and measured the variability per entity. The coefficient of variation (CV) was significantly higher (P < 0.001) in MM (CV = 30.7) than in acute lymphoblastic leukemia (ALL; CV = 7.9), chronic lymphocytic leukemia (CLL; CV = 5.3), and diffuse large B cell lymphoma (DLBCL; CV = 10.4) (Fig. 1C). This analysis shows that the extreme heterogeneity of DNA methylation levels seems to be an epigenetic feature that is specific for MM. To further characterize this variable pattern, we sequenced the whole DNA methylome at a single base pair resolution (>51-fold coverage per sample) of NPCs from bone marrow and two MMs in which the HumanMethylation450 BeadChip indicated extreme hyper- or hypomethylation (white arrows in Fig. 1B; Supplemental Tables 1, 2; Supplemental Fig. 1). Indeed, these analyses validated the array data at the whole-genome level, with the genome of MM1 being hypermethylated and that of MM2 being hypomethylated as compared to NPCs (Fig. 1D–F). Hence, MM can show different methylomes, with a wide range of methylation levels that range from globally hypo- to hypermethylated.

Figure 1.

Whole-genome DNA methylation data of neoplastic and normal plasma cells. (A) Principal component analysis of HumanMethylation450 BeadChip data in MGUS, MM, and NPC from tonsils or bone marrow. (B) Unsupervised hierarchical clustering analysis of MM and MGUS patient samples as well as NPC samples from tonsil or bone marrow from healthy donors (the 20,000 CpG sites with most variable methylation values were used for this analysis). (HyperD) Hyperdiploid MM sample; [t(IGH)] MM sample with IGH translocation; [del(TP53)] MM sample with deletion of TP53. (C) Bar plot showing the coefficient of variation of median methylation values per case in ALL, CLL, DLBCL, MGUS, and MM. (D) Circular representation of DNA methylation levels of purified plasma cells of two MM patient samples (MM1 and MM2) compared to purified NPC obtained from bone marrow from healthy donors. Histograms represent CpG methylation levels averaged in 10-Mbp genomic windows. This representation with low resolution indicates that the shift of methylation levels in MM1 and MM2, as compared to NPCs, takes place throughout the genome and not only in particular chromosomal regions. (E,F) Scatter plots and density color codes for DNA methylation data of all autosomes. Pairwise comparisons of MM1 sample to NPC (E) and MM2 sample to NPC (F) are shown.

Detection of aberrant DNA methylation in MM and MGUS

To obtain a reliable whole-genome measure of the magnitude of epigenetic changes in MM, we compared the methylation levels of each MM sample with extreme methylation patterns and NPCs. Sample MM1 showed 1,504,700 differentially methylated regions (DMRs) covering a total of 2,158,824 CpG dinucleotides, whereas MM2 had 1,875,824 DMRs covering 4,308,958 CpGs. As previously suggested in Figure 1D, the majority of DMRs identified in MM2 were hypomethylated, whereas MM1 showed less hypomethylation and increased hypermethylation levels (Fig. 2A,B).

Figure 2.

Differential DNA methylation in MM and MGUS as compared to NPCs. (A,B) Differentially methylated CpGs using WGBS data in (A) MM1 versus NPC and (B) MM2 versus NPC. (C,D) Differentially methylated CpGs identified using the HumanMethylation450 BeadChip data in (C) MGUS versus NPC samples and (D) MM versus NPC samples. (E) Heatmap of significantly hypomethylated (left) or hypermethylated CpGs (right) in MM samples versus NPC samples. (F) Relative distribution of differentially hypo- or hypermethylated CpGs across different functional chromatin states of the genome using WGBS data (left) or HumanMethylation450 BeadChip data (right) as compared to their respective backgrounds (first column). The numbers inside each cell point to the percentage of CpGs belonging to a particular chromatin state. (WG) Whole-genome; (HyperM) hypermethylated CpGs; (HypoM) hypomethylated CpGs.

Next, we attempted to generate a consensus differential methylation signature of MGUS and MM in comparison with NPCs, using microarray data of a large series of samples (Fig. 2C–E). From 78,998 differentially methylated CpGs identified between MM and NPC, the great majority were hypomethylated (98.3%). In MGUS, we detected 9388 differentially methylated CpGs, of which 6874 (73.2%) were hypomethylated and 2514 (26.8%) hypermethylated (Fig. 2C–E; Supplemental Fig. 2). A direct comparison between MM and MGUS revealed that the progression from MGUS to MM seems to be associated with increased hypomethylation in the absence of additional hypermethylation events (Supplemental Fig. 3).

Applying an unsupervised analysis in order to compare the DNA methylation profiles of clinically relevant MM subtypes based on the presence or absence of cytogenetic alterations, i.e., hyperdiploidy, chromosomal translocations affecting the immunoglobulin heavy chain (IGH) locus and TP53 deletions, we did not find any apparent association (Fig. 1B; Supplemental Fig. 4). These results were confirmed using a supervised approach, in which we observed only a few differentially methylated CpGs between cytogenetic subgroups (Supplemental Fig. 5). However, these findings do not rule out the possibility that other genetic aberrations, not analyzed in our study, may be related to DNA methylation changes, as shown previously in MM as well as in other tumors (Weisenberger et al. 2006; Walker et al. 2011).

Hypo- and hypermethylated CpGs in MM are located in specific genomic regions and enriched in particular chromatin states

Hypomethylated sites in MM and MGUS mainly occurred in regions lacking CpG islands (CGI) and were predominantly located in gene bodies and intergenic regions. Although both whole-genome bisulfite sequencing (WGBS) and HumanMethylation450 BeadChip showed this effect, the enrichment, when compared with the background distribution, was particularly evident in the consensus signature (MM versus NPCs) derived from HumanMethylation450 BeadChip. Interestingly, hypermethylated sites also occurred in regions lacking CGIs (Supplemental Fig. 6).

In order to investigate the functional significance of aberrant DNA methylation in MM, we annotated hypo- and hypermethylated CpGs according to functional chromatin states recently described by ChIP-seq in immortalized mature B cells (IMBCs) (Ernst et al. 2011), which have a DNA methylation profile similar to normal memory B cells and plasma cells (Supplemental Fig. 7). As previously observed in colon cancer (Berman et al. 2011), the great majority (75%–80%) of hypomethylated CpGs in MM and MGUS, or those losing methylation in the progression from MGUS to MM, were located in heterochromatic and late-replicating regions such as lamina-associated domains (LADs) (Fig. 2F; Supplemental Figs. 3, 8). We analyzed the potential effect of this epigenetic change and we found neither a morphological alteration of the chromatin bound to the nuclear lamina by electron microscopy nor a detachment of hypomethylated regions from the nuclear lamina by 3D-FISH. We then studied the transcriptome of 11 MM cases and four NPCs by RNA-seq, and we observed that genes located in hypomethylated heterochromatic regions globally showed low expression levels both in NPCs and MMs regardless of their methylation status (Supplemental Fig. 9). Furthermore, focusing only on promoter regions undergoing demethylation, only 1.9% of the CpGs showed a significant correlation (P < 0.01) with gene expression. These data suggest that heterochromatin hypomethylation seems to have a low impact on gene expression.

Remarkably, a large proportion of hypermethylated CpGs in MM (n = 794; 59.9%; P < 1 × 10−12) (Supplemental Table 3) and MGUS (n = 1656; 65.9%; P < 1 × 10−12) were located within enhancer regions (Figs. 2F, 3A; Supplemental Fig. 10), and in particular, within intronic enhancers (Supplemental Fig. 11). To evaluate whether differential methylation of enhancers in MM may be functionally relevant, we analyzed the expression of their associated genes by RNA-seq. Overall, we observed a shift toward an inverse correlation between enhancer methylation and their associated gene expression as compared to nonenhancer intronic regions of the same genes (Fig. 3B–F; Supplemental Table 4), which mostly took place in the absence of DNA methylation changes in the corresponding promoter regions (Supplemental Fig. 12). In fact, of the 50 genes with the strongest inverse correlation between enhancer methylation and gene expression, only two (4%) showed concurrent promoter hypermethylation. We further analyzed the consequences of differential DNA methylation at enhancer regions and observed that treating a MM cell line with 5-azacytidine led to increased expression of SLC15A4 and PVT1 (Supplemental Fig. 13). Due to the unspecificity of this pharmacological treatment, we also performed reporter assays using a CpG-Free Luciferase Vector (Klug and Rehli 2006). We cloned three enhancer regions into this vector, which showed variable levels of correlation between DNA methylation and expression levels (Supplemental Table 5). Luciferase activity analyzed in transient transfection assays comparing unmethylated and in vitro methylated vectors revealed that upon enhancer methylation, the luciferase activity was reduced, further supporting our findings in primary tumors (Fig. 3D–F).

Figure 3.

Functional and transcriptional analysis of hypermethylated enhancer regions in MM. (A) Percentage of hypermethylated CpGs associated with enhancer regions in MGUS and MM patient samples using HumanMethylation450 BeadChip data. (B) Density plot of correlation coefficients between methylation levels of hypermethylated enhancers and the expression of their associated genes (cyan) using RNA-seq data. As control, intronic CpGs of the same genes were studied excluding the enhancer associated ones (black dotted line). For this analysis, we used 663 CpGs (out of 794) annotated to 574 genes with available gene expression, and 8956 CpGs in nonenhancer intronic regions of the same genes. (C) A snapshot of the UCSC Genome Browser showing the promoter (left) and the 3′ intronic enhancer region (right) of the SLC15A4 gene. Displayed tracks include the chromatin state characterization in IMBCs and ChIP-seq data for H3K27ac, H3K4me1, and H3K4me3. DNA methylation levels of NPCs and MM patient samples measured by WGBS and HumanMethylation450 BeadChip are also shown. (D–F) Correlation analysis between DNA methylation levels of the hypermethylated enhancers and expression of the associated gene (left). In addition, we display the luciferase reporter activity data of the analyzed enhancer region (right) located in the intron of (D) SLC15A4, (E) PVT1, and (F) NCOR2. (MP) Minimal promoter; (MET) methylated; (UNMET) unmethylated.

Hypermethylation of B cell-specific enhancers in MM is associated with enhancer decommissioning

After showing that DNA methylation is related to enhancer activity, we attempted to analyze the biological function of the genes with hypermethylated enhancers (631 genes associated with the 794 CpGs). A Gene Ontology analysis pointed to significant enrichment of various terms associated with B-cell differentiation (Supplemental Table 6). Furthermore, using transcription factor binding sites (TFBSs) identified by ChIP-seq in the ENCODE Project (Gerstein et al. 2012), we observed that hypermethylated enhancers overlapped with binding sites of transcription factors (TFs) associated with B-cell differentiation such as BCL11A, BATF, EBF1, and PAX5 (Supplemental Table 7). These data led us to investigate whether the DNA methylation levels of these enhancer regions are modulated during normal B-cell differentiation. We analyzed the DNA methylation profiles of seven B-cell differentiation stages, including hematopoietic multipotent progenitors, several immature B cells, naive B cells, germinal center B cells, and plasma cells. In addition, we included embryonic stem cells (ESCs) and IMBCs (Ernst et al. 2011) in the analysis, as well as other B-cell neoplasms such as CLL (n = 139) (Kulis et al. 2012), ALL (n = 46) (Busche et al. 2013), and DLBCL (n = 40). Surprisingly, the 794 enhancer-associated CpGs hypermethylated in MM were also methylated in ESCs and gradually lost DNA methylation, finally becoming completely unmethylated in germinal center B cells and terminally differentiated plasma cells (Fig. 4A). We found that ALLs appear to maintain a methylation pattern resembling their cellular origin (precursor B cell), whereas DLBCLs, similarly to MM, show hypermethylation of B cell-specific enhancers when compared to their cell of origin (i.e., germinal center B cells) (Supplemental Fig. 14).

Figure 4.

DNA methylation and chromatin features of hypermethylated enhancers in MM in the context of normal B-cell differentiation. (A) DNA methylation levels of 794 enhancer-associated CpGs in seven B-cell differentiation stages, MGUS and MM patient samples, as well as ESCs (H1), IMBCs (GM12878), and U266-MM cell lines. (B) ChIP-seq levels of H3K4me1 and H3K27ac, and DNase-seq data of 794 enhancer-associated CpGs in ESCs, IMBCs, NPCs, and the U266 MM cell line. (C) Density plot of H3K4me1, H3K27ac, and DNase levels in ESCs, IMBCs, NPCs, and the U266 MM cell line. Among the 794 enhancer-associated CpGs in MM, those hypermethylated in U266 are shown in the top panel, whereas those unmethylated in this cell line appear at the bottom. (MPP-CLP) Hematopoietic multipotent progenitors–common lymphoid progenitors; (ESCs) embryonic stem cells; (IMBCs) immortalized mature B cells; (NPCs) normal plasma cells; (U266-MM) multiple myeloma derived cell line U266.

Next, we studied whether enhancer hypermethylation in MM was accompanied by a change in the chromatin state. We performed ChIP-seq with antibodies against monomethylated lysine 4 (H3K4me1) and acetylated lysine 27 (H3K27ac) of histone H3 in NPCs and in the MM cell line U266, and also DNase-seq in the latter. Data were compared to ESCs (showing hypermethylated enhancers, similar to MM) and IMBCs (with unmethylated enhancers, similar to NPCs). As shown in Figure 4B, the genomic areas surrounding unmethylated enhancers in IMBCs show a clear enrichment of H3K4me1 and H3K27ac as well as DNase hypersensitive sites. In NPCs, the regions under study lost H3K27ac, but H3K4me1 showed a strong signal, indicating that they maintain their enhancer status. In contrast, methylated enhancers in ESCs and the U266 cell line actually represent decommissioned enhancers (Whyte et al. 2012; Rivera and Ren 2013) lacking all the enhancer-associated features (Fig. 4B,C). Interestingly, 125 of the 794 enhancers with significant hypermethylation in MM cases were unmethylated in the U266 MM cell line. Similarly to IMBCs, these unmethylated B-cell enhancers were accompanied by a chromatin structure linked to active enhancers in the MM cell line (Fig. 4C). This finding is further supported by the fact that enhancers (defined as such in IMBCs) hypomethylated in IMBCs, NPCs, and the MM cell line, but hypermethylated in ESCs, only showed active enhancer-related chromatin marks in the first three (Supplemental Fig. 15).

Down-regulation of B cell-specific transcription factors in MM is associated with enhancer hypermethylation

All the preceding results suggest that MMs acquire hypermethylation in enhancer regions that paradoxically undergo demethylation during normal B-cell differentiation (Fig. 4A). Since the process of B-cell differentiation is guided by expression of stage-specific TFs (Matthias and Rolink 2005; Kurosaki et al. 2010), the binding of TFs to their target enhancers may result in demethylation. Several lines of evidence support this model in which DNA demethylation during cell development is mediated by the presence of TFs, whose binding may block DNA methylation maintenance after cell division (Stadler et al. 2011; Lee et al. 2012; Feldmann et al. 2013; Hon et al. 2013; Kulis et al. 2013). Interestingly, here we have observed the reverse process in plasma cell disorders in which B cell-specific TF binding sites seem to reacquire a chromatin state that is typical of undifferentiated cells. Therefore, we postulate that down-regulation of B cell-specific TFs in MM cells could result in enhancer methylation and decommissioning. To test this hypothesis, we analyzed the expression levels of B cell-specific TFs in MMs and NPCs by RNA-seq and correlated them with the mean DNA methylation level of the enhancers with altered methylation in MM. Interestingly, for most of the B cell-specific TFs, we observed an inverse correlation between their expression and methylation levels of their target enhancers (Fig. 5A). Among these, the expression of STAT5, PAX5, NFATC1, and BATF showed a significant negative correlation with DNA methylation levels of enhancers in normal and neoplastic plasma cells (Pearson coefficient ranging from −0.79 to −0.62, P < 0.05) (Fig. 5A–D; Supplemental Table 8). Furthermore, variable expression of several of these TFs was coordinated in MM patients (Fig. 5E), suggesting that the expression level of a TF network, and not that of individual TFs, was associated with the methylation levels of enhancers. In line with this finding, the expression levels of genes regulated by enhancer methylation also showed a correlated expression among MM patients (Supplemental Fig. 16).

Figure 5.

B cell-specific TFs' expression correlates with DNA methylation levels of their binding sites in MM. (A) The upper part shows that hypermethylated enhancers are enriched for binding sites of B cell-specific TFs. The lower part displays the correlation coefficients between TF expression and mean DNA methylation level of their respective binding sites with hypermethylation in MM. (*) P < 0.05. (B–D) Scatter plots showing the association between TF expression and mean methylation level of their respective binding sites. Normal plasma cells are shown in blue and multiple myeloma samples are depicted in red. (E) Correlation matrix of expression levels of TFs among the 11 MM cases with available RNA-seq data. Only those TFs from A with at least 1.5-fold enrichment (in log2) were used for this analysis.

Discussion

In this study, we have characterized the DNA methylome of normal and neoplastic plasma cells using a double approach. On the one hand, we generated an epigenetic signature of common DNA methylation changes in MGUS and MM by comparing them to NPCs, using the HumanMethylation450 BeadChip in a large series of patients. On the other hand, we analyzed the whole DNA methylome at single base-pair resolution of two samples from the extremes of the DNA methylation spectrum of MM. This second approach allowed us to characterize the epigenetic heterogeneity in MM, to define the boundaries of DMRs, and to measure the magnitude of DNA methylation changes in MM at the whole-genome level.

As shown in other studies (Salhia et al. 2010; Walker et al. 2011; Heuck et al. 2013), we describe a clear distinction in DNA methylation patterns between malignant and normal plasma cells, mostly in the form of hypomethylation in the former. Interestingly, one of the major features of MM is the extreme heterogeneity of median DNA methylation levels, which ranges from globally hypo- to globally hypermethylated as compared to NPCs. This heterogeneity in MM is remarkably higher than in other malignancies of the B-cell lineage (Fig. 1C). Premalignant MGUS samples were less heterogeneous than MM but shared a similar hypermethylation signature. In contrast, hypomethylation in MM was much more extensive and heterogeneous than in MGUS, suggesting that it may be related to progression from MGUS to MM, as previously reported (Salhia et al. 2010; Walker et al. 2011; Heuck et al. 2013). In contrast to other studies (Walker et al. 2011), we did not observe the epigenetic variability in MM to be associated with particular genetic changes. However, we cannot rule out the possibility that some genetic features, not analyzed in our study, could be correlated with the heterogeneous DNA methylome. Although the causes of epigenetic heterogeneity in MM require further investigation, we postulate that this phenomenon could be the consequence of both the DNA methylation pattern of the founder plasma cell clone within the heterogeneous plasma cell pool and the proliferative history of the neoplastic cells in a particular case.

One of the most unexpected findings of our study is that DNA hypermethylation in MGUS and MM occurred in regions outside CGIs. This finding contrasts with the classical cancer-associated CIMP (Issa 2004; Fang et al. 2014) and recent whole-methylome reports, in which hypermethylation is invariably associated with CGIs and CGI shores close to promoter regions (Berman et al. 2011; Hansen et al. 2011; Kulis et al. 2012). Remarkably, a large proportion of hypermethylated CpGs in MM and MGUS were located within intronic regions with enhancer-related chromatin marks in nonmalignant mature B cells and plasma cells. Analyzing the methylation of these regions in the context of B-cell differentiation, we observed a gradual demethylation from stem cells to plasma cells. Hence, the observed methylation pattern of these regions in MM is similar to undifferentiated precursor and stem cells. This finding may reveal new potential insights in the biology of the disease. On the one hand, it may represent a de novo epigenetic reprogramming in MM, leading to the acquisition of a methylation pattern related to stemness. On the other hand, we cannot rule out that the observed pattern could reflect an epigenetic imprint of initial premalignant phases of the disease in progenitor cells.

Our analysis of the chromatin in B cells, plasma cells, and MM revealed that enhancer hypermethylation in MM was associated with enhancer decommissioning, as recently shown in prostate and breast cancer (Taberlay et al. 2014). Several studies have provided evidence indicating that the binding of specific TFs mediates the loss of DNA methylation at regulatory elements during the cell differentiation (Stadler et al. 2011; Lee et al. 2012; Feldmann et al. 2013; Hon et al. 2013; Kulis et al. 2013). Here, we have observed that the hypermethylation and decommissioning of specific B-cell enhancers in MM appear to be related to down-regulation of a network of B-cell TFs. Similarly to the globally heterogeneous methylation pattern, this hypermethylation of B-cell enhancers also varied among different cases and may be useful to detect disease subtypes. Interestingly, ∼20% of the MM patients express the B-cell marker CD20 (San Miguel et al. 1991). We speculate that this MM subtype may lack hypermethylation and decommissioning of B-cell enhancers, allowing part of the B-cell signature to be maintained. In conclusion, our results suggest that MMs show a widely heterogeneous DNA methylation pattern. Embedded into an extensive hypomethylated genome, we identified DNA hypermethylation of developmentally regulated enhancers as a new epigenetic feature of MM. The delineation of the precise role of enhancer hypermethylation in MM pathogenesis and its potential usefulness in the clinical setting as a new source of epigenetic biomarkers should be further investigated.

Methods

Patient samples and cell line

Purified plasma cells from bone marrow aspirations were obtained from newly diagnosed patients of MM (n = 104) and MGUS (n = 16). As negative controls, we purified plasma cells from healthy donors (eight tonsils and 12 bone marrows, the latter pooled into three samples of four donors each). The MM cell line U266 was kindly provided by Dr. Cigudosa (Spanish National Cancer Centre-CNIO, Madrid, Spain). Further details on the biological materials used in our study are described in the Supplemental Material.

DNA methylation analysis

Whole-genome bisulfite sequencing and HumanMethylation450 BeadChip (Illumina) were carried out in order to analyze the DNA methylation profile of MM, MGUS, normal plasma cells, and the U266 MM cell line. Details on the experimental procedures, annotation of CpG sites, detection of differentially methylated regions, and Gene Ontology analysis are described in the Supplemental Material.

Analysis of RNA-seq data

RNA samples from 11 MMs (also analyzed by Human Methylation450 BeadChip) and four NPCs from tonsils were sequenced using the standard Illumina protocol. Technical details are described in the Supplemental Material.

Luciferase reporter assay

pCpGL empty vector and pCpGL-CMV/EF1 were kindly provided by Dr. Rehli (University Hospital, Regensburg, Germany). To properly initiate the transcription, a CpG free minimal promoter (5′-AGAGGGTATATAATGGAAGCTTAACTTCCAG-3′) was also cloned in the pCpGL empty vector (see details in the Supplemental Material).

ChIP-seq and DNase-seq

ChIP-seq was performed in NPCs and in the U266 cell line, and DNase-seq was carried out only on the latter. Standard protocols generated within the Blueprint Consortium were followed. Protocol details can be found at http://www.blueprint-epigenome.eu/index.cfm?p=7BF8A4B6-F4FE-861A-2AD57A08D63D0B58.

Mapped reads of the ChIP-seq data for H3K4me1 and H3K27ac as well as DNase-seq data from the cell lines GM12878 and H1 were downloaded from the UCSC Genome Browser (GRCh37/hg19 Assembly; http://genome-euro.ucsc.edu).

Read density around the 794 CpGs hypermethylated in MM was calculated by seqMINER tool using a 5-kb window upstream and downstream from the CpGs of interest (Ye et al. 2011).

Analysis of transcription factor binding sites

We used ChIP-seq data from the ENCODE Project available at the UCSC Genome Browser (http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=wgEncodeHaibTfbs) to obtain information on transcription factor binding sites (TFBS) of a total of 163 TFs. The relative enrichment of each TFBS was calculated for the 794 hypermethylated CpGs located in enhancer elements in comparison with the background. A Fisher's exact test was used, and both an odds ratio and a P-value were assigned to each comparison.

Data access

WGBS, HumanMethylation450 BeadChip, and RNA-seq data have been submitted to the European Genome-phenome Archive (https://www.ebi.ac.uk/ega/) under accession numbers EGAD00001000672, EGAS00001000841, and EGAS00001001110, respectively. ChIP-seq and DNase-seq data from the U266 cell line have been submitted to the European Nucleotide Archive (https://www.ebi.ac.uk/ena/) under sample accession numbers SAMEA2165292 and SAMEA2165299, respectively.