Comprehensive characterization of the epigenetic landscape in Multiple Myeloma

Background: Human multiple myeloma (MM) cell lines (HMCLs) have been widely used to understand the molecular processes that drive MM biology. Epigenetic modifications are involved in MM development, progression, and drug resistance. A comprehensive characterization of the epigenetic landscape of MM would advance our understanding of MM pathophysiology and may attempt to identify new therapeutic targets. Methods: We performed chromatin immunoprecipitation sequencing to analyze histone mark changes (H3K4me1, H3K4me3, H3K9me3, H3K27ac, H3K27me3 and H3K36me3) on 16 HMCLs. Results: Differential analysis of histone modification profiles highlighted links between histone modifications and cytogenetic abnormalities or recurrent mutations. Using histone modifications associated to enhancer regions, we identified super-enhancers (SE) associated with genes involved in MM biology. We also identified promoters of genes enriched in H3K9me3 and H3K27me3 repressive marks associated to potential tumor suppressor functions. The prognostic value of genes associated with repressive domains and SE was used to build two distinct scores identifying high-risk MM patients in two independent cohorts (CoMMpass cohort; n = 674 and Montpellier cohort; n = 69). Finally, we explored H3K4me3 marks comparing drug-resistant and -sensitive HMCLs to identify regions involved in drug resistance. From these data, we developed epigenetic biomarkers based on the H3K4me3 modification predicting MM cell response to lenalidomide and histone deacetylase inhibitors (HDACi). Conclusions: The epigenetic landscape of MM cells represents a unique resource for future biological studies. Furthermore, risk-scores based on SE and repressive regions together with epigenetic biomarkers of drug response could represent new tools for precision medicine in MM.


Introduction
Multiple myeloma (MM) is a B cell neoplasia characterized by the accumulation of clonal plasma cells in the bone marrow. Despite the survival improvement provided by current treatments, the Ivyspring International Publisher majority of patients relapse and eventually become resistant to all treatments. The median overall survival of MM patients is 6 years [1]. MM is characterized by a high degree of biological, genetic and intra-clonal heterogeneity [2,3]. The development of MM is usually accompanied by a series of genetic alterations, such as cytogenetic abnormalities, primary and secondary chromosomal translocations and oncogenic activation. Identifying these alterations is highly valuable in understanding the pathogenesis of MM and predicting outcome and therapeutic response. Chromosomal alterations, including t(4; 14) translocation, 17p13 deletion and 1q21 gain correlate with poor disease prognosis [4]. Several mutations also negatively impact on survival, including mutations in CCND1 and genes involved in DNA repair pathways (TP53, ATM, ATR and ZNFHX4) [5].
The epigenetic control of gene expression plays an essential role in the regulation of cell fate and cell identity. Aberrant changes in key regulatory chromatin features, such as DNA methylation and histone post-translational modifications (PTMs), are involved in MM pathophysiology and drug resistance [6]. The histone PTM landscape is particularly dynamic and constantly evolving [7]. Structural changes in active euchromatin or silenced heterochromatin are controlled by epigenetic enzyme complexes of chromatin writers, readers and erasers. Histone deacetylases (HDAC) are dysregulated in MM with an aberrant overexpression of class I HDACs in a subset of patients, which is associated with a shorter overall survival [8]. HDAC inhibitors (HDACi) are now used in the treatment of several hematologic malignancies, including MM [9]. Moreover, epigenetic biomarkers for MM can be used as predictive and prognostic indicators to guide diagnosis and treatment [10][11][12][13]. These data suggest that individualizing therapy based on epigenetic biomarkers has the potential to increase therapeutic efficacy for MM patients.
Human MM cell lines (HMCLs) have been widely used for the understanding of MM biology, identifying target genes, screening of anti-myeloma drugs and, more recently, characterizing the MM mutational landscape [14,15]. In the past few years, we have derived a large cohort of patient-derived HMCLs that remain dependent on the addition of exogenous MM growth factors, reflecting primary tumor conditions [16]. Using these myeloma cell lines as well as commercial myeloma cell lines, we described their molecular diversity by analyzing their gene expression profile and mutational landscape and showed that HMCL molecular diversity reflects part of the molecular heterogeneity of primary MM cells [14,16]. However, the global epigenetic landscape of HMCLs has never been described.
A comprehensive characterization of the epigenetic landscape of HMCLs would advance our understanding of MM pathophysiology and its analysis might lead to the identification of new therapeutic targets for future application in MM. In this study, we present the first epigenetic landscape of HMCLs. We generated chromatin immuno precipitation sequencing (ChIP-seq) maps for a set of active (H3K4me1, H3K4me3, H3K27ac, H3K36me3) and inactive (H3K9me3, H3K27me3) histone PTMs on 16 HMCLs, representative of the molecular heterogeneity of MM [14,16]. Differential analysis of PTM profiles on HMCLs highlighted links between histone modifications and cytogenetic abnormalities. We have identified super-enhancers (SE) associated with genes involved in the biology of MM, such as MAF, MYC, CCND1, CCND2, TRAF3 or NSD2, and repressive domains, characterized by the co-localization of the H3K9me3 and H3K27me3 modifications on promoters of genes with potential tumor suppressor function. We used genes belonging to these two subsets and associated with prognostic value to build two distinct scores that can predict MM patient survival. Finally, we have identified H3K4me3 active promoters associated with the sensitivity of HMCL cells to lenalidomide and the HDACi romidepsin. From these data, we have developed epigenetic biomarkers based on this H3K4me3 PTM predicting the response to lenalidomide and romidepsin HDACi. Together, these data constitute a rich source of information that might guide analyses towards the identification of novel biomarkers and therapeutic targets for MM.
Bone marrow samples were collected after patient's written informed consent in accordance with the Declaration of Helsinki and institutional research board approval from Montpellier University Hospital. Bone marrow was collected from 69 patients at diagnosis and 28 patients at relapse, this cohort was called Montpellier cohort. MM cells of patients were purified using anti-CD138 MACS microbeads (Miltenyi Biotec, Bergisch Gladbach, Germany). We also used RNA sequencing data of 674 newly diagnosed MM patients with longitudinal follow-up from the Multiple Myeloma Research Foundation CoMMpass trial (NCT01454297; version IA11a), termed in the following CoMMpass cohort.
Normal plasma cells were generated using a 3-step in vitro model starting from purified memory B cells from 3 different healthy donors as reported [17,18].

Gene expression profiling and statistical analysis
HMCLs, patients MM cells and normal plasma cells RNA sequencing was done as previously described [12]. The RNA sequencing (RNA-seq) library preparation was done with 150 ng of input RNA using the Illumina TruSeq Stranded mRNA Library Prep Kit. Paired-end RNA-seq were performed with Illumina NextSeq sequencing instrument (Helixio, Clermont-Ferrand, France). RNA-seq read pairs were mapped to the reference human GRCh37 genome using the STAR aligner [20]. All statistical analyses were performed with the statistics software R (version 3.6.2; https://www.rproject.org), and R packages developed by Bioconductor project (https://www.bioconductor. org) [21]. The expression level of each gene was summarized and normalized using the DESseq R/Bioconductor package [22].

Whole Exome sequencing (WES) and variant calling
WES of human MM cell line was performed as described above [14]. WES public data were used for KMS-12-BM [28]. Reads were mapped to the reference human hg19 using the Bowtie 2 aligner version 2.3.2 [29]. SAMtools version 1.5 was used to convert SAM files to sorted BAM files [30]. Indel realignment, base quality recalibration and variant calling steps were completed with GATK 3.8-1 [31]. ANNOVAR was used to annotate variants [32]. Then, filters described in Vikova et al. were applicated [14].
ChIP-seq for H3K4me1, H3K4me3, H3K9me3, H3K27ac, H3K27me3 and H3K36me3 histone marks of 4 MM patients were performed using purified MM cells (MMC) of 4 MM patients as previously described [33]. LiftOver function from the rtracklayer R package (version 1.52.1) was used to change the genomic coordinate of MMC peaks, from hg38 to hg19. Peaks identified in MMC were compared to HMCL peaks for each histone marks using the findOverlaps function from the regioneR R package (version 1.24.0).
Reads were mapped to the human reference genome (hg19) using Bowtie2 (version 2.3.2) [29]. Peak-calling was performed using MACS2 (version 2.1.2) [34]. Peak annotation and differential binding analyses of ChIP-seq peak data were performed using ChIPseeker [35] and DiffBind [36] R/bioconductor packages, respectively. Heatmaps and average profiles were generated by deeptools (version 3.5.0) [37]. Super-enhancers (SE) were identified using the ROSE (Rank Ordering of Super-Enhancers) algorithm [38,39] based on H3K27ac and H3K4me1 ChIP-seq signals. TSS exclusion zone size was adjusted to 2500 bp to exclude promoter regions. Transcriptionally active genes were assigned to super-enhancers using a simple proximity rule (50 kb window). Scores were built using our previously published methodology [10]. Genes composing super-enhancer and repressive domain scores were selected for their prognostic significance using Maxstat R function and the Benjamini Hochberg multiple testing correction. Score values correspond to the sum of the Cox beta-coefficients of each gene, weighted by ± 1 if patient MM cells signal for the gene of interest is above or below the Maxstat reference value of this gene. The statistical significance of differences in overall survival between patients' groups was calculated using the log-rank test and survival curves were plotted using the Kaplan-Meier method. CoMMpass and Montpellier cohorts were used as training and validation cohorts, respectively. Statistical differentially bound regions between drug-sensitive and -resistant HMCLs were identified with DiffBind R/bioconductor package (FDR ≤ 0.05).

Epigenetic landscape of HMCLs
We performed ChIP-seq analysis using 16 HMCLs characterized by molecular and genetic variability, which reflect part of the heterogeneity of the disease (Table 1) [14,16]. To evaluate the epigenetic landscape of these HMCLs, two repressive marks (H3K9me3 and H3K27me3) and four active marks (H3K4me1, H3K4me3, H3K27ac and H3K36me3) were selected. The six histone marks used in this study are those proposed by the International Human Epigenome Consortium as the most informative ones and constitute an essential part of the reference epigenome (http://ihec-epigenomes. org/research/reference-epigenome-standards/). Specifically, H3K36me3 is related to transcriptional elongation in the gene body, H3K4me3 to poised/active promoters, H3K4me1 to poised/active enhancers, H3K27ac to determine active regulatory elements (both enhancers and promoters), H3K9me3 to detect heterochromatin and H3K27me3 to detect Polycomb-repressed chromatin. The annotation of the peaks reveals that the enhancer-related H3K4me1 modification is mainly found in distal intergenic (45.4% ±2.6%) and intron (35.4% ±2.4%) regions whereas H3K27ac, a modification known to be associated with both active enhancers and promoters, is mainly found in distal intergenic (34.8% ±4.0%) and promoter (32.9% ±4.6%) regions ( Figure 1A). H3K4me3 is mainly located at promoters, as expected for a mark of actively transcribed gene promoters, but also in intergenic regions (promoters: 42.1% ±6.8%; intergenic: 36.2% ±5.0%). The H3K4me3 and H3K27ac marks exhibit the same profile surrounding the transcription start site (TSS) with a major peak located downstream the TSS ( Figure 1B-D). The H3K4me1 mark is located both upstream and downstream to the TSS ( Figure 1B-D). On enhancer regions, the majority of H3K27ac peaks overlapped with H3K4me1 peaks, revealing active enhancers, whereas non-overlapping H3K4me1 peaks were associated with poised enhancers ( Figure 1D) [40]. The H3K36me3 modification is identified downstream to the TSS, on the body of active genes with a distribution of intron and exon regions representing 69.8% of total regions. The repressive marks, H3K9me3 and H3K27me3, respectively detecting HP1-and Polycomb-repressed chromatin, were essentially found in intergenic regions (61.6% ±6.2% and 54.1% ±4.3%, respectively), with distinct localization, but also on gene promoters ( Figure 1A-D). A small fraction of these promoters was enriched by H3K9me3 and H3K27me3 co-localization. Moreover, we compared our results with the ChIP-seq data of 4 MM patients [33]. We identified that more than 80% of MM patients peaks were commonly found in HMCL data for H3K4me1, H3K4me3, H3K9me3, H3K27me3 and H3K36me3 histones marks. Concerning H3K27ac, 70% of overlap was identified comparing HMCLs and primary MMC from patients (Supplementary Table 1). Figure 2A shows the global profile of histone modifications on the 16 HMCLs. The clustering based on peak localizations distinguishes each mark from all others and repressive (H3K9me3 and H3K27me3) from active (H3K4me1, H3K4me3, H3K27ac and H3K36me3) marks. TSSs of actively transcribed genes are marked by H3K4me3 and H3K27ac modifications, which explains the high similarity between H3K4me3 and H3K27ac profiles ( Figure 1C and 1D) [40]. When each mark is observed independently, the clustering of HMCLs reveals subgroups associated to different cytogenetic abnormalities or HMCLs transcriptional classification [16]. Differential analysis on H3K9me3 reveals a highly similar cluster, composed of KSM-12-BM, OPM2, RPMI8226, XG19 and SKMM2, characterized by the association of del1p, del13q and del17p cytogenetic abnormalities ( Figure 2B and Supplementary Table 2  (r > 0.5). The t(4;14) translocation results in IgH/MMSET hybrid transcripts inducing overexpression of the MMSET gene [41]. There is a correlation between the number of H3K27me3 peaks and MMSET mRNA expression in the groups of HMCLs without t(4;14) translocation (R 2 = 0.33; P-value < 0.02) (Supplementary Figure 1A) whereas no correlation is observed in t(4;14) cell lines (data not shown). However, the percentage of reads in H3K27me3 peaks is correlated to the MMSET mRNA expression in all HMCLs (R 2 = 0.48; P-value < 0.01) (Supplementary Figure 1B). These results suggest that H3K27me3 localizations are very similar in t(4; 14) HMCLs and that the global enrichment of H3K27me3 increases with the deregulation of MMSET without increasing the total number of H3K27me3 sites on the genome. Popovic et al. showed that MMSET overexpression led to a global loss of H3K27me3 and an enrichment of this modification on specific loci [42]. Overexpression of MMSET is associated with a shift in the genomic localization of the H3K27 methyltransferase EZH2, imposing EZH2 and H3K27me3 accumulation at specific loci [43]. Differentially bound site analysis revealed around 46000 sites significantly differentially associated to the  Table 4). Global analyses of H3K4me1, H3K4me3, H3K27ac and H3K36me3 modifications do not provide a clear link between cytogenetic abnormalities and histone modifications, while H3K27ac and H3K36me3 clustering tends to group HMCLs by HMCL molecular subgroups ( Figures 2E-G). We also studied the link between histone modification profiles and frequently mutated genes in HMCLs (Supplementary Figure 3). Interestingly, HMCLs with TP53 bi-allelic events are characterized by a specific H3K9me3 profile (Supplementary Figure 4).  Patients of the CoMMpass cohort (n = 674) were ranked according to the increased SE-risk score and a maximum difference in OS was obtained with SE-risk score of -4.0 splitting patients into high-risk (n = 283; red curve) and low-risk (n = 391; green curve) groups. SE-risk score also had a prognostic value in an independent cohort of 69 patients (Montpellier cohort). To compute the SE-risk score in Montpellier cohort, the proportion of patients of the CoMMpass cohort for each gene was applied on the Montpellier cohort to determine gene cut-point. The same gene beta-coefficient and score cut-point was applied to distinguish high-risk (n = 34: red curve) from low-risk (n = 35; green curve) groups. (C) SE-associated gene score correlates with the response to GSK525762 in HMCLs. Linear regression analysis of the SE-associated gene score in function of the IC50 of GSK525762 in 8 HMCLs. Coefficient of determination R 2 represents the square of the Pearson correlation coefficient (r) (Pearson correlation test).

SE-risk score predicting survival in MM patients
H3K4me1 and H3K27ac modifications are known to be associated with enhancer regions and, more specifically, active enhancers can be identified by co-occupancy of H3K4me1 and H3K27ac [40]. Super-enhancers consist of a set of large enhancer domains displaying physical proximity (+/-12.5kb) and are associated with genes that control and define cell identity [38,39]. Using the MM1S cell line, Lovén et al. identified super-enhancers associated with key MM genes, including MYC, IRF4, PRDM1, XBP1, CCND2, PIM1, BCL-xL and MCL1 [39].We identified a number ranging from 607 to 2510 predicted super-enhancers in each of the HMCLs (Supplementary Table 5). These super-enhancers differed from typical enhancers in both size and H3K4me1 and H3K27ac levels ( Figure 3A). Enhancers tend to loop to and associate with adjacent genes in order to activate their transcription [44]. Using a simple proximity rule, we assigned all TSSs to super-enhancers within a 50 kb window [45]. Among the SE-associated genes identified, we found genes identified in the study of Lovén et al. (MYC, IRF4, PRDM1, XBP1 and CCND2) but also other key MM genes, such as ACTG1, MAF, CCND1, TRAF3 and NSD2 (MMSET). Among them, we found SE-associated genes targeted by MM translocations (Table 1 and Figure 3B).
We used these SE-associated genes to build a score predicting survival in MM patients. First, we selected SE-associated genes with strong expression in HMCLs, overexpressed in plasma cells of MM patients compared to normal plasma cells and associated with shorter overall survival in MM patients using the Maxstat algorithm [46] Figures 5 and 6). Then, the prognostic information of these 28 SE-associated genes was combined to build a SE-risk score. This score is defined by the sum of the beta coefficients of the Cox model for each prognostic gene, weighted by ±1 if the patient plasma cells signal for a given gene is above or below the Maxstat reference value of this gene as previously described [47]. Using two independent cohorts, the SE-risk score had a prognostic value when used to split patients into two groups using the Maxtstat R function. The score splits patients into a high-risk group (score > -4.0) and a low-risk group (score < -4.0) in the CoMMpass (P-value < 0.0001; high-risk group: 58.0% of patients; low-risk group: 42.0% of patients) and Montpellier (P-value < 0.005; high-risk group: 50.7% of patients; low-risk group: 49.3% of patients) cohorts ( Figure 4B). Investigating the link between the SE-risk score and cytogenetic abnormalities, SE-risk score values were significantly higher in MM patients with del(13q), del(17p), del(1p), 1q gain and t (12;14) in the CoMMpass cohort (Supplementary Figure 7). We evaluated the evolution of the score from diagnosis to relapse in 14 MM patients with paired samples. No significant difference between the two groups was identified (Supplementary Figure 8). However, 2 MM patients presented a striking increase in SE-risk score from diagnosis to relapse. We next investigated the SE-risk score value distribution according to Affymetrix GEP-based risk scores previously reported in MM. The SE-risk score values were significantly higher in high-risk patients defined by RS score [25], and IFM score [27]. The SE-risk score values were also significantly higher in intermediate and high-risk patients compared to low-risk patients defined by the International Staging System (ISS) [48]. Furthermore, high-risk patients defined with the SE-risk score demonstrated a significant increase in the percentage of proliferating MM cells determined by BrdU incorporation [49] (median = 1.55%, range: 0 -17.3%) compared to the low group (median = 0.7%; range: 0 -7.3%) (Supplementary Figure 9). Bromodomain and extra-terminal (BET) inhibitors have been shown to repress super-enhancer-associated transcription. To evaluate the link between the score and the response of MM cells to bromodomain and extra-terminal (BET) inhibitors, we calculated the SE-associated gene score in HMCLs. Interestingly, we found a significant anti-correlation between the score and the response to a BET inhibitor in clinical development in hematological cancers (GSK525762) in HMCLs (R 2 = 0.4081; P-value < 0.05) ( Figure 4C). The same tendency was obtained using JQ1 and OTX015 BET inhibitors (Data not shown). Using other scores based on gene expression profiling (EZ [12], HA [11] and HRS [26] scores) calculated in 8 HMCLs, we did not find a significant correlation between these scores and the BET inhibitor sensitivity (Supplementary Figure 10). Nevertheless, a correlation was found between HRS score and the response to GSK525762 (R 2 = 0.489; P-value < 0.05). High SE-associated gene score allows the identification of high-risk MM patients that could benefit from BET inhibitors. Conversely, high HRS score is associated with resistance to GSK525762 BET inhibitor. Patients of the CoMMpass cohort (n = 674) were ranked according to the increased H3K9me3/H3K27me3 score and a maximum difference in OS was obtained with increased H3K9me3/H3K27me3 score of -3.45 splitting patients into high-risk (n = 390; red curve) and low-risk (n = 284; green curve) groups. High-risk group for individual gene composing H3K9me3/H3K27me3 score is associated with low gene expression (negative beta-coefficient) whereas high-risk group for H3K9me3/H3K27me3 score is associated with high score because of the multiplication of the beta-coefficient (negative value) and the weight (+1 or -1 if the signal of MM patient gene is above or below the Maxstat reference value, respectively) of genes. H3K9me3/H3K27me3 score also had a prognostic value in an independent cohort of 69 patients (Montpellier cohort). To compute the H3K9me3/H3K27me3 score in Montpellier cohort, the same method described before was used.  Figure  5), suggesting that genes may play a role in disease progression and independence of MM cells from micro-environment. Moreover, DEMETER2 [50] and CERES [51] dependency scores calculated from public datasets of RNAi [52,53] and CRISPR/Cas9-based [51] screening (Dependency Map data, Broad Institute, www.depmap.org) indicated that YWHAQ, IL10, HK2 and THY1 are significant essential MM genes (Supplementary Table 6 and Supplementary Figure  11). Indeed, THY1 knockdown presented the larger DEMETER2 dependency score difference between HMCLs and all other cell types investigated. Altogether, these data show that a signature of genes associated with super-enhancers in MM could be defined allowing the identification of high-risk MM patients that could benefit from treatment with BET inhibitors.

H3K9me3 and H3K27me3 co-localization in promoters can be associated with prognostic value in MM
Several reports have demonstrated that H3K9me3 and H3K27me3 modifications are mutually exclusive in constitutive heterochromatin. It was previously suggested that H3K9me3 may crosstalk with the Polycomb mediated H3K27me3 modification to cooperate in gene repression ( Figure 1D) [54][55][56]. We decided to find gene promoters enriched by both H3K9me3 and H3K27me3 modifications to build a score based on genes associated with epigenetic transcriptional repression and identify new potential tumor suppressor genes. First, H3K9me3 and H3K27me3 overlapping regions on gene promoters were identified for all HMCLs. To identify potential tumor suppressor genes, we selected genes significantly under-expressed in plasma cells of MM patients compared to normal plasma cells, and associated them with poor outcome when their expression is low in MM cells of patients. We identified eighteen genes : ARHGEF5, BIVM, DEF8,  GRID2IP, HDAC9, HSPA1L, KDM4C, NLRP2, P4HA3,  PAG1, PM20D1, RMND5A, SEMA6A, SFMBT2, THEMIS2, TPRKB, ZFP2 and ZNF5188B (Supplementary Figures 12 and 13). The promoters of SEMA6A and ARHGEF5 genes were most frequently enriched by both H3K9me3 and H3K27me3 modifications in HMCLs ( Figure 5). We combined the prognostic information from these eighteen genes into a tumor suppressor-based risk score as described above. This score splits patients into two groups in the training cohort (CoMMpass cohort, n=674, P-value < 0.0001): a low-risk group (score < -3.45) of 57.9% MM patients associated with global high expression of these genes and a high-risk group (score > -3.45) of 42.1% MM patients associated with global low gene expression ( Figure 5B). This score also splits MM patients into low and high-risk groups (62.3% and 37.7% of patients, respectively) in the validation cohort (Montpellier cohort, n=69, P-value < 0.01). In the CoMMpass cohort, the score was significantly higher in non-hyperdiploid, del13q, del17p, 1q gain, t(4;14), t(12;14) and t(14;16) subgroups of patients (Supplementary Figure 14). Of the eighteen genes, thirteen genes (ARHGEF5, DEF8, GRID2IP, HDAC9, HSPA1L, KDM4C, P4HA3, PAG1, SEMA6A, SFMBT2, THEMIS2, ZFP2 and ZNF518B) were down-regulated in HMCLs compared to primary MM cells (Supplementary Figure 12). We evaluated the link between the score and the response of MM cells to EZH2 and SUV39H1/2 inhibitors, EPZ-6438 [57] and chaetocine [58], respectively. Interestingly, we found a significant anti-correlation between the score and the response to EPZ-6438 (R 2 = 0.4348; P-value < 0.05) ( Figure 5C) and chaetocine (R 2 = 0.4128; P-value < 0.05) ( Figure 5D) in HMCLs.

H3K4me3 bound sites linked to lenalidomide and romidepsin resistance in HMCLs
Epigenetic modifications, including histone modifications, participate to MM pathogenesis and chemotherapy resistance. To better understand the epigenetic mechanisms involved in drug resistance, we explored histone marks in drug-resistant and -sensitive HMCLs. H3K4me3, a modification associated with promoters of actively transcribed genes, was chosen to identify epigenetic changes involved in drug resistance. We analyzed the relationship between the H3K4me3 landscape of HMCLs and their response to conventional drugs including bortezomib, melphalan, dexamethasone, lenalidomide and HDACi [14]. Of major interest, we identified 5903 significantly differentially H3K4me3-enriched sites between lenalidomide-resistant (AMO1, JJN3, KMS-12-BM, SKMM2, XG1, XG5, XG7, XG12, XG20 and XG21 HMCLs) and lenalidomide-sensitive cell lines (OPM2, RPMI8226, XG2, XG6, XG13 and XG19 HMCLs) ( Figure 6A) (FDR ≤ 0.05). Regarding gene promoters, 203 H3K4me3 differentially bound sites were significantly enriched in lenalidomide-sensitive HMCLs compared to lenalidomide-resistant HMCLs (FDR ≤ 0.05) whereas 136 H3K4me3 differentially bound sites were significantly enriched in the lenalidomide-resistant group compared to the lenalidomide-sensitive group ( Figure 6B). Examples of significantly differentially enriched regions are presented in Figure 6C. Gene set enrichment analysis revealed that H3K4me3 differentially bound sites enriched in lenalidomide-resistant HMCLs were associated to interferon signaling and cytokine signaling in the immune system ( Figure 6D). Differential sites enriched in lenalidomide-sensitive HMCLs were significantly associated with an NFKBIA target gene signature.
Lenalidomide targets the Cereblon complex, formed by CUL4, RBX1, DDB1 and CRBN proteins [59]. Interestingly, we observed a H3K4me3 enrichment on the promoter of CUL4B in lenalidomide-sensitive compared to resistant HMCLs (Supplementary Figure 15). The presence of CUL4B isoform 1 appeared to be associated with lenalidomide sensitivity. To validate this observation, we generated XG1 cell line derivatives stably expressing GFP tagged-Cul4B isoform 1 (Cul4B-GFP cell line) or the GFP alone. The overexpression of Cul4B isoform1 expression was validated both by western blot and RT-qPCR ( Figure 7A and 7B). Cells were sorted by flow cytometry based on their GFP levels in order to keep those with high expression of Cul4B-GFP ( Figure 7C). Nuclear localization of the ectopic Cul4B-GFP was confirmed by microscopy ( Figure 7D). Sensitivity to lenalidomide of the XG1 Cul4B-GFP and the control XG1-GFP cell lines were studied by proliferation assays. MM cells expressing Cul4B-GFP were significantly more sensitive to lenalidomide compared to the control ( Figure 7E). Moreover, Cul4B-GFP expression also increased cell sensitivity to pomalidome, a third generation IMiD that also targets CRBN [59] ( Figure 7F). Altogether, these data indicate that the expression of Cul4B isoform 1 is related to the response to IMiDs in MM cells and support the potential interest of ChIP-seq analyses to identify biomarkers that could predict response to these drugs. To confirm the potential of ChIP sequencing data to identify epigenetic mechanisms involved in drug resistance, we investigate the bound sites differentially enriched in H3K4me3 between romidepsin-resistant (AMO1, XG7, XG13, XG19 and XG20) and romidepsin-sensitive (OPM2, RPMI8226, SKMM2, XG1, XG5, XG6 and XG12) cell lines (Supplementary Figure 16A) (FDR ≤ 0.05). Among the 2098 H3K4me3 significantly differentially bound sites, 250 bound sites were localized on gene promoters (Supplementary Figure 16B). Romidepsin is a class I histone deacetylase inhibitor currently evaluated in combination with lenalidomide in phase I/II (NCT017755975) in patients with relapse or refractory multiple myeloma. Among the 250 H3K4me3 bound sites located on gene promoters, 180 were significantly enriched in romidepsin-sensitive HMCLs associated with an HDAC3 and PRDM6 target gene signature (Supplementary Figure 16D). We found 70 differentially bound sites enriched in the romidepsin-resistant group compared to the romidepsin-sensitive group. Gene set enrichment analysis revealed a significant association with H3K27me3 and polycomb complex subunit (SUZ12 and EED) target genes. That indicates that ChIP-seq strategies could be an interesting tool to identify biomarkers associated with drug response or resistance in MM. Further investigations are needed to establish potential future clinical applications.
Moreover, among the 6 histones modifications, only the differentially bound sites identified with H3K4me1, H3K4me3 and H3K9me3 allowed the distinction between resistant and sensitive cell lines for lenalidomide and romidepsin drugs (Supplementary Figure 17 and 18). Using H3K9me3, only few bound sites were significantly differentially enriched between resistant and sensitive groups. It could be interesting to carry out more investigations on the H3K4me1 differentially bound sites to identify super-enhancers associated to drug resistance. Altogether, these data underline the interest of ChIP sequencing to identify epigenetic biomarkers related to drug response in MM.

Discussion
In this study, we have characterized the landscape of histone modifications in a large collection of HMCLs representative of MM molecular heterogeneity [16]. Some histone modification profiles, such as H3K27me3, were clearly related to cytogenetic abnormalities. In t(4;14) HMCLs, MMSET overexpression results in the genome-wide redistribution of the SET domain protein EZH2 and its associated H3K27me3 mark [60]. Gene expression profiling revealed deregulated genes involved in cell cycle (e.g. CCNE2), apoptosis (e.g. BAX and BCL2) and DNA repair (e.g. ATM and GADD45A) [61]. Moreover, we showed that promoter of genes involved in the negative regulation of cell cycle (e.g. CEND1 and CDKN1C) and the induction of apoptosis (e.g. TP73 and BNIP3) were enriched in the H3K27me3 mark in t(4;14) HMCLs. Methylation of the BNIP3 promoter was also described and associated with poor overall survival in MM patients [62]. These results confirmed the interplay between MMSET, EZH2 and H3K27me3 [61,63] and revealed new potential therapeutic targets for t(4;14) MM patients associated with poor prognosis [64]. We also observed that TP53 bi-allelic events are associated with a H3K9me3 signature. SUV39H1, the histone methyltransferase responsible for the H3K9me3 mark, is a target of p53 repression and p53 target promoters are enriched in the H3K9me3 repressive mark [65], suggesting that TP53 mutation could modulate H3K9me3 level on these promoters and thus the p53 apoptotic response.
It is widely believed that H3K9me3 and H3K27me3 do not co-occur at the same loci, but some ChIP-seq data indicate that these two marks can be found together at a subset of developmentally regulated gene in mouse embryonic stem cells, extra-embryonic lineages and human differentiated cells [54][55][56]. This dual repression at specialized regulators may point to the importance of maintaining their silencing to confer a selective advantage in tumoral cells. Here, we identified 18 potential suppressor tumor genes enriched in H3K9me3 and H3K27me3 repressive marks on their promoter. To our knowledge, the tumor suppressor function of these genes has not yet been described in MM. However, some of them are involved in cell proliferation inhibition. The NLRP2 gene, coding for a member of the nucleotide binding and leucine-rich repeat receptor (NLR), is an inhibitor of NF-κB pathway inhibitor which plays a key role in survival and proliferation of myeloma cells [94,95]. PAG1 encodes for a type III transmembrane adaptor protein and is an inhibitor of Src tyrosine kinases known to promote proliferation in MM [96,97]. From these 18 genes we built a score identifying a high-risk group of MM patients with global low expression of putative suppressor genes associated with H3K9me3 and H3K27me3 marks. We showed that this score was correlated with the response to EPZ-3864 and chaetocin, a SUV39H1/2 protein inhibitor, in HMCLs. High-risk MM patients identified with this score could benefit from combined therapeutic targeting of EZH2 and SUV39H1 histone methyltransferases. Moreover, the therapeutic interest of EZH2 inhibitors has already been demonstrated in HMCLs and primary myeloma samples [98]. Concerning SUV39H1, high expression level is associated with a poor prognosis in MM patients and SUV39H1 inhibitor also exhibited anti-MM effects both in HMCLs and primary samples [99].
Despite the improvement of MM patient survival through the development of novel agents, including new generation of IMiDs or proteasome inhibitors, monoclonal antibodies and HDAC inhibitors, the acquisition of drug resistance is the major limitation of MM therapy. The great majority of MM patients ultimately relapse and the treatment of relapse/refractory MM remains a major challenge. In this study, we analyzed the relationship between differential enrichment of H3K4me3 on gene promoters in HMCLs and response to treatment. We developed an epigenetic biomarker based on this histone modification, which predicts lenalidomide and romidepsin responses in HMCLs. IMiDs are known to target the cereblon complex. Cereblon is composed of CUL4, RBX1, DDB1 and CRBN proteins and induces the ubiquitination of B-cell transcription factors IKZF1 and IKZF3 in presence of IMiDs [59]. Among the genes encoding these proteins, we observed a significant difference in H3K4me3 enrichment on the promoter of CUL4B in lenalidomide-sensitive HMCLs compared to resistant HMCLs (Supplementary Figure 15). In other studies, loss of CUL4B has been described as conferring resistance to lenalidomide in lymphoma and myeloma cell lines [100,101]. Here, we have shown that the presence of the CUL4B splicing variant 1 is indeed associated with sensitivity to lenalidomide and pomalidomide.
This study provides a comprehensive characterization of the MM epigenetic landscape, representing a unique resource for future biological studies which could help in identifying novel critical epigenetic modifications involved in MM progression and drug resistance. Furthermore, risk-scores based on super enhancers and repressive regions together with epigenetic biomarkers of drug response could represent new tools for precision oncology in MM.

Author Contributions
EA performed research, the data analyses and participated in the writing of the paper. LH, HDB, SB, GP, NR, GR, SO, MK and AB participated in the research or in data analyses. GCartron and LV participated in clinical data analysis and participated in the writing of the paper. AMM, JIMS and GCavalli participated in the research and in the writing of the paper. JM supervised the research and the writing of the paper.