ARTICLE https://doi.org/10.1038/s41467-020-20715-x OPEN Targeting aberrant DNA methylation in mesenchymal stromal cells as a treatment for myeloma bone disease Antonio Garcia-Gomez1,2,11 ✉, Tianlu Li1,2,11, Carlos de la Calle-Fabregat1,2, Javier Rodríguez-Ubreva 1,2, Laura Ciudad1,2, Francesc Català-Moll1,2, Gerard Godoy-Tena1, Montserrat Martín-Sánchez3, Laura San-Segundo3, Sandra Muntión3, Xabier Morales 4, Carlos Ortiz-de-Solórzano 4, Julen Oyarzabal5, 1234567890():,; Edurne San José-Enériz6, Manel Esteller 7,8,9,10, Xabier Agirre 6, Felipe Prosper 6, Mercedes Garayoa 3 & Esteban Ballestar 1,2 ✉ Multiple myeloma (MM) progression and myeloma-associated bone disease (MBD) are highly dependent on bone marrow mesenchymal stromal cells (MSCs). MM-MSCs exhibit abnormal transcriptomes, suggesting the involvement of epigenetic mechanisms governing their tumor-promoting functions and prolonged osteoblast suppression. Here, we identify widespread DNA methylation alterations of bone marrow-isolated MSCs from distinct MM stages, particularly in Homeobox genes involved in osteogenic differentiation that associate with their aberrant expression. Moreover, these DNA methylation changes are recapitulated in vitro by exposing MSCs from healthy individuals to MM cells. Pharmacological targeting of DNMTs and G9a with dual inhibitor CM-272 reverts the expression of hypermethylated osteogenic regulators and promotes osteoblast differentiation of myeloma MSCs. Most importantly, CM-272 treatment prevents tumor-associated bone loss and reduces tumor burden in a murine myeloma model. Our results demonstrate that epigenetic aberrancies mediate the impairment of bone formation in MM, and its targeting by CM-272 is able to reverse MBD. 1 Epigenetics and Immune Disease Group, Josep Carreras Leukaemia Research Institute (IJC), Badalona, 08916 Badalona, Barcelona, Spain. 2 Chromatin and Disease Group, Cancer Epigenetics and Biology Programme (PEBC), Bellvitge Biomedical Research Institute (IDIBELL), 08908 L’Hospitalet de Llobregat, Barcelona, Spain. 3 Centro de Investigación del Cáncer, IBMCC (Universidad de Salamanca-CSIC) and Hospital Universitario de Salamanca-IBSAL, 37007 Salamanca, Spain. 4 Imaging Platform, Center for Applied Medical Research (CIMA), University of Navarra, IDISNA, Ciberonc, 31008 Pamplona, Spain. 5 Small Molecule Discovery Platform, Molecular Therapeutics Program, Center for Applied Medical Research (CIMA), University of Navarra, 31008 Pamplona, Spain. 6 Division of Hemato-Oncology, Center for Applied Medical Research (CIMA), University of Navarra, IDISNA, Ciberonc 31008 Pamplona, Spain. 7 Josep Carreras Leukaemia Research Institute (IJC), Badalona, Barcelona, Catalonia, Spain. 8 Centro de Investigacion Biomedica en Red Cancer (CIBERONC), Madrid, Spain. 9 Institucio Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Catalonia, Spain. 10 Physiological Sciences Department, School of Medicine and Health Sciences, University of Barcelona (UB), Barcelona, Catalonia, Spain. 11These authors contributed equally: Antonio Garcia- Gomez, Tianlu Li. ✉email: nonogarcia85@gmail.com; eballestar@carrerasresearch.org NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x M ultiple myeloma (MM) is an incurable hematological Results malignancy characterized by clonal expansion of BM-derived MSCs of distinct MM stages exhibit altered DNA plasma cells in the bone marrow (BM) that accounts methylation profiles. We first obtained genome-wide DNA for 1% of all cancers1,2. Nearly 90% of myeloma patients suffer methylation profiles of BM-derived MSCs isolated at different from skeletal-related events during the course of the disease, stages of MM (newly diagnosed MGUS, high-risk SMM, and including severe bone pain, hypercalcemia, pathological frac- MM) and healthy controls. DNA methylation changes were tures, and spinal cord compression3, that not only affect the identified using two different statistical approaches (Fig. 1A): (i) quality of life but also their overall survival4. Myeloma- detection of differentially methylated CpG positions (DMPs) associated bone disease (MBD) is characterized by an increase based on differences in DNA methylation means between the in bone-resorptive activity and number of osteoclasts (OCs), as patient (MGUS, SMM, and MM) and healthy MSCs (Δβ ≥ 0.1 well as impairment of bone-forming activity and differentiation and **p < 0.01) (Supplementary Data 1); and (ii) detection of of osteoblasts (OBs), which ultimately lead to the development differentially variable CpG positions (DVPs) based on differences of osteolytic lesions5. in variance of DNA methylation levels (qval < 0.05 and *p < 0.05) In most cases, symptomatic myeloma is preceded by between the sample groups (Supplementary Data 2). In regards to sequential asymptomatic stages of monoclonal gammopathy of DMPs, the largest number of altered CpGs was found in MSCs undetermined significance (MGUS) and smoldering myeloma from patients of the SMM stage compared to healthy donors (SMM), with increasing BM plasmocytosis and monoclonal (Supplementary Fig. 1A, B). On the other hand, we observed the component as well as augmented risk of progression to active highest number of DVPs in comparison to healthy donors in MM6,7. The biological behavior and clinical outcome of MM MSCs isolated from MGUS followed by SMM and MM patients are partly dependent on genetic and epigenetic abnormalities of (Supplementary Fig. 1C, D), supporting the notion that these tumor subclones that arise from MGUS and SMM stages8. stochastic and heterogeneous DNA methylation patterns are However, the clinical stability of MGUS cases, despite dis- associated with early stages of carcinogenesis, as previously playing shared genetic lesions with MM cells, suggests that the reported24,25. We also observed that the majority of identified BM microenvironment may critically modulate disease DMPs and DVPs are disease stage-specific, although the progression6,9,10. In this regard, it has been widely shown that a asymptomatic stages showed a moderate proportion of overlap complex and bidirectional relationship exists between MM cells (Fig. 1A). and the BM niche, which results in oncogenesis support, ane- Given that myeloma is a multi-stage disease, we then analyzed mia, immunosuppression, and uncoupling of the bone remo- the accumulative changes of DNA methylation associated deling process11. with MM progression by selecting DMPs (Supplementary Data 3) Mesenchymal stromal cells (MSCs) are an essential cell type and DVPs (Supplementary Data 4) that were found either only in in the formation and function of the BM microenvironment, the MM stage, shared by SMM and MM and in all three being the progenitors of bone-forming OBs, adipocytes, and stages (Fig. 1A). With this analysis, we identified 872 hyper- chondroblasts, as well as the hematopoietic-supporting stroma and 951 hypomethylated DMPs, and 260 hyper- and 318 components of the BM12. It is well-documented that BM- hypomethylated DVPs. derived MSCs from MM patients contribute to MM progression Analyzing the distribution of MM progression-associated (reviewed in ref. 11) and show an impaired ability to differ- CpGs in relation to CpG islands (CGI), we observed a significant entiate into OBs13,14. Moreover, MM-MSCs are considered enrichment of CpGs in open sea regions in the hypomethylated inherently abnormal, as their dysfunctionality remains even DMP data set (Fig. 1B) and in CpG islands in the hypermethy- following ex vivo culture in the absence of MM cells15. Fur- lated DVP data set (Supplementary Fig. 1E). Utilizing publicly thermore, bone lesions persist in many MM patients even after available chromatin state maps of BM-derived MSCs from therapeutic remission, suggesting a long-term defect in MSCs healthy individuals26, we found a significant enrichment of both that inhibit their ability to properly differentiate into functional hyper- and hypomethylated DMPs sites that correspond to OBs16. enhancers (Fig. 1C). In addition, we observed an enrichment in Previous studies described that MSCs from MM patients are flanking transcription start sites (TSS) in the hypomethylated set, cytogenetically normal17,18, but show alterations in their and regions repressed by Polycomb Group (PcG) in hypermethy- transcriptional13,19 and proteomic11 profiles even in the absence lated CpGs (Fig. 1C). On the other hand, hypermethylated DVPs of myeloma cell interaction. This suggests that epigenetic were enriched in TSSs, bivalent regions, and regions repressed by mechanisms could be governing the tumor-promoting functions Polycomb (Supplementary Fig. 1F). of MSCs and their prolonged OB suppression in MM. In fact, To determine whether these MM progression-associated loci Adamik and colleagues20 reported abnormal recruitment of shared any common DNA elements, we performed a search for chromatin remodelers in MSCs from myeloma patients, con- enriched transcription factor (TF)-binding sites in these regions tributing to the transcriptional repression of Runx2, a master using the HOMER algorithm27. We observed a significant over- regulator of OB differentiation. Yet, there is a lack of information representation of binding sites for the Runt and Tead family in about DNA methylation-related mechanisms that may contribute differentially hypermethylated DMPs associated with MM to MM progression and subsequent bone defects. DNA methy- progression (**p < 0.01; Fig. 1D). These results suggest that lation is an essential epigenetic modification involving the addi- key transcription factors involved in the upregulation of tion of a methyl group to the 5-carbon of the cytosine ring by a osteogenic genes, such as RUNX228 or TEAD229, may family of DNA methyltransferase (DNMT) enzymes21, which has participate in aberrant DNA hypermethylation. Since DNA been described to play a critical role in MSC lineage determina- methylation has been originally linked to transcriptional tion22, as well as in tumor progression and immunosuppression repression, these results suggested that the hypermethylation in other cancer types23. of these regions could compromise the ability of MSCs to In this study, we identify DNA methylation alterations in undergo proper OB differentiation. On the other hand, DMP MSCs of MM patients mediated by MM cells resulting in the sites that experienced aberrant DNA hypomethylation were dysregulation of osteogenesis, and this is reversed by the treat- highly enriched in binding motifs of the bZip and Homeobox ment with CM-272, a dual inhibitor of DNMTs and the histone families (**p < 0.01; Fig. 1D). In this respect, the loss of DNA methyltransferase G9a. methylation could be selectively driving the occupancy of TF 2 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE A DMPs MGUS SMM DNA methylation MGUS 572 514 9618 SMM MM 445 303 315 1063 CpG island context se HD y lth MM ea ea MGUS EPIC Chromatin states is H D BM-MSCs SMM array Progression Motif enrichment MM DVPs MGUS SMM DNA methylation Gene expression 840 253 521 53 86 28 497 se y lth MM ea ea is H B C D log (odd ratio) FDR Significant hyper- 1 1 NO 1e−06 0 Relation to CGI 1e−12 YES hypo- −1 1e−18 −2 Island N_Shelf 1e−24 EnhG ZNF/Rpts Enh Quies TssA ReprPC Tx Het TxWk TxFlnk TssAFlnk ReprPCWk BivFlnk EnhBiv TssBiv N_Shore OpenSea S_Shelf S_Shore 100 D TF family NF1:FOXA1 bHLH RUNX1 hyper- Nrf2 hypo- hypo- bZIP 1.5 2.00 CTF 75 NFE2L2 EBF Percentage Bach1 HOXA2 Forkhead Gata1 HOXA1 HMG Fold of enrichment RUNX-AML Homeobox 50 1.4 HOXC9 1.75 HTH RUNX2 NF−E2 Jun-AP1 Fra1 RUNX Fosl2 MADS JunB Hoxb4 NR PRDM1 Six2 AP−1 Runt Mef2d Tcf7 1.3 Gata2 Atf3 25 TEA CEBP 1.50 ERE Fra2 BATF TEAD Gata6 Pdx1 Zf Hand2 Foxa2 PBX2 Gata4 Bach2 TEAD4 Foxo3 Mef2b 0 1.2 TEAD3 Gata3 Foxa3 Pit1+1bp Foxl2 FOXK1 FOXA1 p value 1.25 bg hyper- hypo- FOXP1 FOXM1 Pdx1 FOXA1 0.05 to 1e−03 FOXK2 TEAD Bcl6 Foxo1 EBF1 1e−03 to 1e−05 Bcl6 TEAD4 Fox:Ebo SCL 1.1 Hoxa13 Foxf1 AMYB TEAD1 BMYB NF1−halfsite 1e−05 to 1e−08 TEAD3 0 25 50 75 100 0 25 50 75 100 < 1e−8 % of windows containing motif % of windows containing motif E OSM SFRP2 IL7 CHSY1 COL4A1 GFI1 cg25666403 cg00705808 cg13896190 cg03199330 cg27657439 cg11954479 1.0 1.0 1.0 1.0 1.0 0.5 ** * ** 0.9 ** * * * * 0.8 0.9 * 0.4 ** 0.8 0.8 * * 0.8 * β value 0.6 0.8 0.3 0.7 0.6 0.6 0.4 0.7 0.2 0.6 0.4 0.4 0.5 0.2 0.6 0.1 0.4 0.0 0.2 0.2 0.5 0.0 D S M M D S M M D S M M D S M M D S M M D S M M U H U U H H SM M SM M SM M U U H U H H SM M SM M SM M G G G G G G M M M M M M that have been reported as negative regulators of OB (Supplementary Fig. 1G). In all, these findings suggested that differentiation such as HOXA230 and ATF331. In addition, we MM progression-associated DNA methylation changes in observed transcriptional deregulation of some members of MSCs might be mediated by the sequential activity of specific these TF families using expression array data from BM-derived TF families, which are also functionally deregulated in MM32. MSCs of healthy controls, MGUS, SMM, and MM patients. Furthermore, other genes that play important roles in the Some of these TFs were specifically downregulated in MSCs of pathophysiology of MM (such as the cytokines IL6 and OSM) active myeloma (RUNX2 and TEAD2), others were already and associated MBD (secreted factors such as RANKL, SFRP2, downregulated in precursor myeloma stages (HOXC9 and IL7, CHSY1, COL4A1, and the transcriptional repressor GFI1) CEBPD), whereas other TFs, including HOXA2 and ATF3, did were also found to alter their DNA methylation levels (Fig. 1E not change their expression in any myeloma stages and Supplementary Fig. 1H). NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x Fig. 1 High-throughput stepwise DNA methylation changes in BM-derived MSCs associated with MM progression. A Workflow depicting the methodological approach for selecting DNA methylation changes in bone marrow-derived mesenchymal stromal cells (BM-MSCs) from monoclonal gammopathy of undetermined significance (MGUS; n = 10), smoldering myeloma (SMM; n = 8), and multiple myeloma (MM; n = 9) patients versus healthy controls (HD; n = 8). An example of a CpG site experiencing increased mean (differentially methylated position, DMP) or variance (differentially variable position, DVP) in the disease versus the control condition is shown. Venn diagrams show the number of DMPs or DVPs resulting from each comparison. B Distribution of DNA methylation changes in relation to CpG islands (CGI), including shores (south, S; north, N), shelves (south, S; north, N), and open sea regions for differentially hyper- or hypomethylated CpG sites. C Enrichment analysis of differentially hyper- and hypomethylated CpG sites located in different genomic regions, annotated by 15 chromHMM states. Color scale refers to log odd ratio and circle size refers to p-value significance. D Bubble plot representation of HOMER transcription factor (TF) motif enrichment analysis of differentially hyper- and hypomethylated CpGs in MSCs during MM progression (left and right panel, respectively). Color range depicts different transcription factor families and circle size refers to p-value significance. E Box plots showing β-values from DMPs obtained from the EPIC array in MSCs from healthy donors and MGUS, SMM, and MM patients of relevant genes involved in the pathogenesis of MM and associated bone disease. HD is represented in dark blue, MGUS in light blue, SMM in orange, and MM in red. eBayes-moderated ANOVA t-test was performed to calculate statistical significance (*p < 0.05, **p < 0.01, and ***p < 0.005). Aberrant DNA methylation is associated with differential association between DNA methylation at gene promoter and gene Homeobox gene expression in MSCs at different MM stages. expression (Fig. 2D and Supplementary Fig. 2A, B). To further investigate the relationship between differential DNA We then validated the aforementioned DNA methylation and methylation and gene expression, we mapped the DMPs to the gene expression changes in an independent cohort of BM-derived most proximal gene. Using expression array data from BM- MSCs from different MM disease stages by pyrosequencing and derived MSCs of healthy controls, MGUS, SMM, and MM real-time quantitative PCR. Among the differentially methylated patients (Supplementary Data 5), differential expression of DMP- genes of the Homeobox family, we selected HOXA2, -A4, -A6, associated genes was identified using a cutoff of *p < 0.05 com- -A9, -A10, -B7, -C9, -C10, and PITX1 on the basis of their paring MGUS/SMM/MM to healthy controls for both DNA reported role in MSC pluripotency (Fig. 2E, F and Supplementary methylation and gene expression (Fig. 2A and Supplementary Fig. 2C, D)40. Furthermore, we validated differentially methylated Data 6). Gene ontology (GO) analysis revealed that the genes genes with osteogenic roles in the myeloma context, including displaying both differential methylation and expression were RUNX2 and IBSP (Fig. 2E, F and Supplementary Fig. 2C, D). In enriched in functional categories important in cell fate commit- most cases, we observed that DNA methylation negatively ment and bone phenotype (Fig. 2B). The most enriched func- correlated with gene expression. tional category corresponded to genes from the Homeobox family. Within the Homeobox family, we found the subset of Hox genes that encodes a large family of highly conserved TFs Healthy MSCs change their DNA methylation profile to one responsible for driving the correct differentiation of MSCs33, partially resembling that of MSCs from MM patients upon namely genes belonging to the HOXA-to-D clusters. Further- interaction with MM plasma cells. To address the potential more, we observed a significant enrichment in genes reported to contribution of MM cells in mediating aberrant DNA methyla- be downregulated in MM-MSCs (Fig. 2B)13. Integration of tion changes in MSCs, we evaluated whether the epigenetic methylation and gene expression data corresponding to the changes observed in MM-MSCs could be mimicked in vitro by Homeobox and bone formation-related genes revealed that both direct contact of healthy MSCs with MM cells. Thus, we co- DNA hyper- and hypomethylation events were associated with cultured BM-derived MSCs from healthy donors with the human both downregulation and upregulation of gene expression in MM cell line MM.1S for 2 weeks Subsequently, MSCs were sorted different genomic locations (Fig. 2C). Specifically, hypermethy- by CD13+ expression and subjected to DNA methylation analysis lated genes that showed a reduced expression in patient MSCs (Fig. 3A). include positive regulators of OB differentiation such as RUNX2 Under these conditions, MM.1S cells were able to induce the or NRP234 (Fig. 2C and Supplementary Table 1) In contrast, expression of genes known to be upregulated in MM-MSCs negative regulators of osteogenesis such as SFRP235 or NFATC236 (IL1B, IL6, and HGF) in HD-MSCs compared with mono- were hypomethylated and consequently upregulated in patient cultured HD-MSCs (Fig. 3B). Additionally, we validated the MSCs (Fig. 2C). In all, these factors could potentially contribute inhibitory effect of MM cells in MSC-to-OB differentiation and to impaired osteoblastogenesis associated with bone disease in observed a decrease in both ALP activity and OB mineralization MM and this is summarized in Supplementary Table 1. in OBs differentiated in the presence of conditioned media from Upon a closer inspection of several Homeobox-associated the MM.1S cell line as compared to OBs differentiated alone genomic regions, we observed a negative correlation between (Fig. 3C). DNA methylation of promoters and gene expression. Specifically, We then investigated the DNA methylation profiling of MSCs the HOXA gene cluster showed aberrant hypomethylation at the from healthy donors generated upon interaction with MM cells. HOXA4 promoter, and its gene expression was upregulated at We observed that 142 CpGs that change their methylation levels different disease stages. Conversely, gene promoters of HOX-A6, upon co-culture with MM.1S cells were shared with aberrant -A7, -A9, -A10, and -A11 displayed hypermethylation and these DNA profiles found in MSCs isolated from MGUS/SMM/MM genes were downregulated in MGUS/SMM/MM (Fig. 2D and patients (Fig. 3D and Supplementary Data 7). Although this Supplementary Fig. 2A). A similar pattern of an inverse accounted for a small percentage of DMPs identified in MM association between methylation and expression was observed patient MSCs, GO analyses revealed enrichment in Homeobox in the HOXB and HOXC gene cluster, where HOXB5, -C5, and genes and categories related to bone formation, similar to that -C8 were aberrantly hypomethylated and upregulated, whereas observed in primary patient MSCs (Fig. 3E). Specifically, we HOXC9, -C10, and -C11 were hypermethylated and down- found that healthy MSCs exposed to MM cells underwent gains regulated in patients (Fig. 2D and Supplementary Fig. 2A). Other (HOXA9, ACVR2A, EBF2) and losses (HOXA2, HOXA3, Homeobox genes such as TBX5, PITX1, or EMX2 were also HOXC5) of DNA methylation in the direction of those observed reported as regulators of bone formation37–39 and showed an for MSCs from myeloma patients (Fig. 3F). 4 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE hypo- MGUS SMM MM HD A C MAB21L2 hyper- MAB21L1 MGUS EBF2 SMM HOXB5 HDAC9 MM HD promoter−TSS SFRP2 BMPR1B HOXC10 PITX1 MGUS SMM HOXC10 PITX1 HOXA11 HOXA9 TWIST1 HOXC9 HOXC11 RUNX3 95 44 119 EMX2 DOWN- HOXC11 Promoter-TSS HOXA10 MAB21L2 MAB21L1 53 TBX15 HOXC8 11 35 FGFR2 exon PITX1 HOXA7 EMX2 EN1 ISL1 MAB21L1 74 HOXD4 NFATC1 PITX1 TBX15 5' HDAC9 HOXA11 RARA GLI3 HOXA9 MM HOXC9 ISL1 EMX2 HOXC11 exon HOXC10 RUNX1 IRX5 SOX9 FGFR2 HOXA6 MGUS SMM TBX5 intron CADM1 RUNX3 EMX2 MSX2 HOXC11 MAB21L1 HOXA10 96 64 133 TBX15 MAB21L2 EN1 HOXA9 UP- EBF2 PITX1 58 HOXC9 HOXB5 IL15 18 35 TWIST1 JAG1 NRP2 HOXA4 COL12A1 −1 0 1 intron RUNX1 IGF1R BMPR1B 66 Gene expression SOX9 HOXA6 EMX2 ALX3 HOXC11 RUNX2 SOX6 3' PITX1 HOXC8 HOXD4 MM HDAC9 MAB21L1 ZEB2 EMX2 TBX5 HOXC11 B EBF2 MAB21L1 PITX1 IRX5 TTS ISL1 HOXC5 HOXA10 MAB21L2 TTS HOXA9 HOXC9 HOXA4 HOXA7 JAG1 TBX5 homeobox domain EN1 HOXA4 FGFRL1 GO annotation EBF2 OSR1 skeletal system development TBX5 MSX2 EN1 NFATC2 OSR1 SFRP2 HOXD4 abnormal skeleton morphology JAG1 intergenic SNAI2 HOXC8 intergenic GLI2 HOXC8 EMX2 pattern specification process NRP2 HOXC11 HOXC10 EPHB3 TWIST1 IRX5 genes down-regulated in MM-MSCs HOXC9 SEMA3A CADM1 RUNX1 SOX6 ISL1 TCF3 0 50 100 150 200 HOXA9 HOXA6 ISL1 PITX1 -log (p value) non−coding HOXC9 HOXA10 non−coding HOXC11 EMX2 TBX5 HOXA7 MAB21L2 D hyper- hyper- DMP DVP hypo- hypo- 50 kb 50 kb DMP MGUS DVP MGUS SMM SMM MM MM HOXA cluster A1 A2 A3 A4 A5 A6 A7 A9 A10 A11 A13 HOXC cluster C13 C12 C11 C10 C9 C8 C6C5 C4 100 kb 10 kb MGUS MGUS SMM SMM MM MM HOXB cluster B1 B2 B3 B4 B5 B7 B9 B13 PITX1 B6 B8 E HOXA6 HOXA9 HOXA10 HOXB7 HOXC9 PITX1 RUNX2 cg01210554 cg05139500 cg15864691 cg07547765 cg09045007 cg069335074 cg26366735 15 100 60 80 100 * * * ** * * ** ** * * DNA methylation ** ** 80 * 100 100 60 * 10 80 40 60 60 40 50 40 50 5 40 20 20 20 20 0 0 0 0 0 0 D S M M D S M M D S M M D S M M D S M M D S M M D S M M U U U H H H U U SM M U U SM M H SM M H H H SM M SM M SM M SM M G G G G G G G M M M F M M M M HOXA6 HOXA9 HOXA10 HOXB7 HOXC9 PITX1 RUNX2 2.5 8 15 15 10 Relative gene expression ** * * 30 * 8 * * ** 2.0 * ** ** 6 ** 6 10 10 1.5 20 5 4 4 1.0 5 5 10 2 0.5 2 0 0 0.0 0 0 0 0 D S M M D S M M D S M M D S M M D S M M D S M M D S M M U H U U U U U U H H H H H H SM M SM M SM M SM M SM M SM M SM M G G G G G G G M M M M M M M Moreover, we validated the DNA methylation and gene gene expression coupled with hypermethylation of these loci expression changes in healthy MSCs driven by co-culture with (Fig. 3G, H). MM.1S with another MM cell line, RPMI-8226. Here, we Altogether, these results support the notion that MM cells not observed a similar effect of co-culture with RPMI-8226, in only are capable of inducing changes in the global methylome of which there was a clear inhibition of Homeobox and osteogenic MSCs but also have a significant impact at specific osteogenic loci. NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x Fig. 2 DNA methylation changes are associated with differential gene expression of Homeobox genes in MSCs from MGUS, SMM and MM patients. A Venn diagrams showing differentially methylated or downregulated (upper) and upregulated (lower panel) genes when comparing MGUS (blue), SMM (orange), and MM (red) samples with healthy individuals. B Gene ontology (GO) enrichment analysis of CpG sites undergoing DNA methylation and gene expression changes in MSCs of patients compared to controls using the GREAT online tool. A binomial test was performed to calculate statistical significance. C Heatmaps showing gene expression of Homeobox and other OB-related genes associated with differentially hyper- (left) or hypomethylated (right) CpG sites. Heatmaps are grouped according to the genomic location (promoter, TSS, 5′, exon, intron, intergenic, 3′ or non-coding region) of analyzed CpG sites. Color scale ranging from light yellow to dark blue represents low to high expression levels. D Scheme depicting differentially methylated and variable CpG sites located in the Homeobox genes (HOX-A, -B, and -C clusters) and PITX1. Dark blue lines indicate hypomethylated DMPs, light blue lines indicate hypomethylated DVPs, dark red lines indicate hypermethylated DMPs and light red lines indicate hypermethylated DVPs associated with MGUS, SMM, and MM condition. E DNA methylation analysis by pyrosequencing of selected CpGs located at the promoter regions and F gene expression of HOXA6, -A9, -A10, -B7, -C9, PITX1, and RUNX2 in MSCs from healthy controls (dark blue; n = 17), MGUS (light blue; n = 8), SMM (orange; n = 8), and MM (n = 16) patients. Gene expression was normalized against RPL38. Box plots represent median ± IQR and whiskers represent maximum and minimum. Statistical significance was calculated using unpaired two-tailed Student t-tests (*p-value < 0.05, **p-value < 0.01, and ***p-value < 0.005). Dual targeting of DNMTs and G9a restores Homeobox gene utilizing a transwell system avoiding contact between both cell expression in vitro and promotes osteogenic differentiation of types, we observed that soluble factor secreted by MM.1S and mesenchymal precursors. Gene expression analysis of DNMTs in RPMI-8226 cell lines were sufficient to change the expression of MSCs from HD and MM patients co-cultured with MM cells several OB-relevant genes in healthy MSCs, including RUNX2, obtained from a previous study41 showed an aberrant upregula- SPP1, IBSP, and HOXB7 (Fig. 5A), concordantly to direct co- tion of the DNA methyltransferase DNMT1 (Fig. 4A). DNMT1 culture. Furthermore, treatment with the dual inhibitor CM-272 interacts with the methyltransferase G9a to coordinate DNA and was able to partially reverse those changes in gene expression H3K9 methylation during cell replication42 promoting tran- mediated by soluble factors secreted by MM cells (Fig. 5A). scriptional silencing of target genes. Moreover, G9a can suppress Changes in gene expression were accompanied by inverse transcription by inducing DNA methylation in addition to its changes in DNA methylation in some of the genes (IBSP, activity as a chromatin remodeler43. In this regard, we hypothe- HOXB7) (Fig. 5B), which were also observed to at least be sized that the dual inhibition of DNMT1 and G9a could reactivate partially reversed by CM-272. However, for some other genes hypermethylated and silenced genes of MSCs from MM patients such as RUNX2 and SPP1, transwell co-culture with MM cells preserving their osteogenic potential and therefore preventing induced minimal effect on DNA methylation, suggesting that myeloma-associated bone loss. Thus, we utilized a dual inhibitor direct cell–cell contact may be required. Moreover, CM-272 was of DNMTs and G9a, termed CM-272, which has been previously able to partially reverse the MM cell lines-mediated inhibitory described to have a potent therapeutic response, both in vitro and effect on OB mineralization (Fig. 5C). Altogether, these results in vivo, in other neoplasias44–47. suggest that MM cells at least partially exert its effects on MSCs We first checked the effect of CM-272 on the cell viability of through secretory mechanisms, and treatment with CM-272 was mesenchymal progenitors and we selected a dose (50 nM) with no able to reserve these effects through the inhibition of DNA significant toxicity in order to perform further experiments methylation. (Fig. 4B). CM-272 treatment was able to restore the expression of Homeobox genes (HOX-A6, -A9, -A10, -C9, PITX1, and RUNX2) CM-272 not only controls tumor burden but also prevents the that were epigenetically repressed in MSCs from MM patients myeloma-associated bone loss. To test the effect of CM-272 in (Fig. 4C). Mechanistically, we observed a loss of DNA methyla- the context of MBD, we used an established murine model of tion in the promoter region of the majority of the aforementioned bone marrow-disseminated myeloma. After equivalent engraft- genes after CM-272 treatment in MM-MSCs (Fig. 4D). We then ment of myeloma cells (RPMI-8226-luc) was verified by biolu- checked the levels of the inactive chromatin mark H3K9me2, a minescence measurement, mice were treated for 4 weeks with hallmark of methyltransferase G9a activity, at these gene CM-272 as described in Methods. Compared with the vehicle promoters upon CM-272 treatment. The chIP-qPCR analysis control group, CM-272 controlled tumor progression as mea- showed a decrease in H3K9me2 levels at the promoter regions of sured by bioluminescence (Fig. 6A) or by serum levels of hIgλ Homeobox genes after CM-272 treatment (Fig. 4E). Taken secreted by RPMI-8226 cells (Fig. 6B). Representative microCT together, our results suggest that CM-272 acts in vitro by images at the metaphyses of distal femurs showed a tumor- inhibition of both DNMT and G9a methyltransferase activity. associated bone loss in vehicle-treated mice, in contrast with Next, we addressed whether targeting DNMT and G9a may trabecular structures observed in CM-272-treated animals have a role in regulating osteogenic differentiation. For this (Fig. 6C). In the vehicle control group, 3D reconstruction images purpose, we cultured MSCs from myeloma patients in osteogenic of distal femurs revealed a marked bone loss evidenced by a thin media to obtain differentiated OBs in the presence or absence of trabecular network (in red) but also by loss of cortical bone (in CM-272. As observed in Fig. 4F, G, CM-272 was able to increase gray) in vehicle-treated mice (Fig. 6D). By contrast, CM-272- ALP activity in early-stage OBs. Furthermore, CM-272 treatment treated mice presented a gain in both trabecular and cortical bone was able to upregulate the relative expression of several late bone (Fig. 6D). This was also reflected by bone morphometric para- formation markers (namely, bone siaploprotein, osteopontin, and meters that resulted in increased trabecular bone volume, occu- osteocalcin) in MSCs from myeloma patients (Fig. 4H). pancy, and connectivity and reduced trabecular separation in Previous research has described that MM cells exert their effect CM-272-treated animals, as compared with vehicle control on MSCs through both direct cell–cell contact and soluble factor (Fig. 6E). Finally, these findings correlated with a significant mechanisms11,48. Our main data show that direct myeloma-MSC increase in serum levels of the bone formation marker P1NP co-culture conditions are able to induce changes in the MSC analyzed after CM-272 treatment compared to untreated control methylome; however, it would be of particular interest to (Fig. 6F). In summary, these data demonstrate that CM-272 investigate whether MM cells may also mediate the same changes exerts in vivo anti-myeloma activity along with bone-anabolic only through soluble factor mechanisms. For this purpose, effects in human MM-bearing mice. 6 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE A 3 10 10 3 10 2 10 2 MSC 14 days 10 1 10 1 HD-MSCs MM.1S 0 0 10 10 SSC-A SSC-A CD13+ MSC sorting 0 200 400 600 800 1000 0 1 2 3 10 10 10 10 FSC-A CD13-PE hyper- hypo- B C D HD-MSC HD-MSC HD-MSC + MM.1S HD-MSC + MM.1S HD-MSC + MM.1S 4 1.5 -M + M B M Relative gene expression C -O B -O D D H 3 H *** nM Alizarin Red 1.0 Normalized ALP activity DNA methylation 2 −2 0 2 Mineralization 0.5 1 HD-MSC 0 0 HD1 HD2 HD3 HD1MM HD2MM HD3MM HD1 HD2 HD3 HD1MM HD2MM HD3MM -M + IL1B IL6 HGF B M B -O M C -O D D H H E F hyper- hypo- 1.0 1.0 * * HD-MSC 0.8 0.8 * HD-MSC + MM.1S skeletal system development GO annotation β values 0.6 0.6 homeobox domain abnormal skeleton development 0.4 0.4 * * * pattern specification process genes down-regulated in MM compared to MGUS 0.2 0.2 0 2 4 6 8 0.0 0.0 ACVR2A HOXA9 EBF2 HOXA3 HOXA2 HOXC5 -log (p value) G HOXA6 HOXA9 HOXA10 HOXB7 HOXC9 PITX1 RUNX2 IBSP cg01210554 cg05139500 cg15864691 cg07547765 cg09045007 cg069335074 cg26366735 cg24986876 1.5 1.5 1.5 1.5 1.5 1.5 1.5 2.0 * DNA methylation * * 1.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 MSC 1.0 MSC + MM.1S MSC + RPMI-8226 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4 1.5 3 1.5 1.5 2.0 2.0 2.0 * * DNA methylation 3 ** * 1.5 1.5 1.5 1.0 2 1.0 1.0 2 1.0 1.0 1.0 0.5 1 0.5 0.5 1 0.5 0.5 0.5 0 0.0 0 0.0 0.0 0.0 0.0 0.0 H HOXA6 HOXA9 HOXA10 HOXB7 HOXC9 PITX1 RUNX2 IBSP Relative gene expression Relative gene expression 3 2.0 2.0 2.0 2.0 2.0 2.0 10 * 8 1.5 1.5 * 1.5 * 1.5 1.5 1.5 MSC 2 6 1.0 1.0 1.0 1.0 1.0 MSC + MM.1S 1.0 4 MSC + RPMI-8226 1 0.5 0.5 0.5 0.5 0.5 0.5 2 0 0.0 0.0 0.0 0.0 0.0 0.0 0 3 2.0 2.0 2.0 2.0 2.0 3 10 * * 8 1.5 * 1.5 ** 1.5 1.5 1.5 * 2 * * 2 * 6 1.0 1.0 1.0 1.0 1.0 4 1 1 0.5 0.5 0.5 0.5 0.5 2 0 0.0 0.0 0.0 0.0 0.0 0 0 To further examine the in vivo effect of CM-272 on DNA Myeloma-bearing mice that were treated with CM-272 displayed methylation of myeloma-associated MSCs, we performed reduced a partial reversion of aberrant hypermethylation of MSCs caused representation bisulfite sequencing (RRBS) analysis of MSCs by the presence of myeloma cells (Fig. 6G). These DNA isolated from vehicle- and CM-272-treated myeloma-bearing methylation changes occurred at genomic loci enriched for genes mice using healthy mice as controls. First, we observed significant involved in cell commitment and differentiation, such as alterations in the DNA methylome of MSCs from myeloma- Homeobox genes (Fig. 6H). Specifically, we were able to identify bearing mice compared to healthy mice (Supplementary Data 8). CpGs that experienced a gain in DNA methylation in vehicle- NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x Fig. 3 MSCs from healthy donors recapitulate DNA methylation changes observed in MSCs from MM patients upon interacting with MM plasma cells. A Scheme depicting workflow (left panel) and sorting strategy (right panel) for selecting CD13 + MSCs after 14 days of co-culture with the MM.1S cell line. B Gene expression analysis of IL1B, IL6, and HGF, normalized against RPL38, of the hMSC-TERT cell line (HD-MSC) co-cultured with the MM.1S cell line for 14 days. Expression levels were normalized against MSCs in monoculture. Data are represented as the mean ± SEM from three independent experiments. C ALP activity and matrix mineralization were assessed in differentiated osteoblasts from HD-MSCs in the presence (HD-OB + CM-MM) or absence (HD- OB) of conditioned media from the MM.1S cell line (CM-MM). Representative images of each experimental condition are shown. Barplot represents the mean ± SEM from six independent experiments, and a paired two-tailed t-test was performed to evaluate statistical significance (***p < 0.005). D Heatmap showing differentially methylated CpG sites (eBayes-moderated paired t-test *p < 0.05) in sorted HD-MSCs (three independent donors) in monoculture (HD1-3) or co-cultured with the MM.1S cell line (HD1-3MM) for 14 days that overlaps with previously identified DMPs. The color scale from blue to red represents low to high methylation levels. E GO enrichment analysis of DMPs in HD-MSCs co-cultured with the MM.1S cell line overlapping with myeloma- associated DMPs analyzed using the GREAT online tool. p-values were calculated using a binomial test. F Bar plots showing β-values obtained from the DNA methylation array presented in D, representing mean and ±SEM of three independent experiments. G DNA methylation and H gene expression levels of DMPs validated in HD-MSCs monoculture (blue) or co-cultured with MM.1S (orange) or RPMI-8226 (red) cell lines as indicated. Gene expression data were normalized against RPL38 and all data were normalized against HD-MSC monoculture. Statistical significance was calculated using paired one-tailed Student t-tests and bar plots represent mean ± SEM of 3–5 independent experiments (*p-value < 0.05, **p-value < 0.01, and ***p-value < 0.005). treated MSCs compared to healthy controls at the same genomic MM arises from early transcriptional deregulation of loci previously identified in human MM-MSCs, including Homeobox genes, and altered DNA methylation may be the HOXA7, -B4, -B7, -B9, and -D10 (Fig. 6I and Supplementary primary mediator in this process. Nevertheless, we cannot Data 9). Importantly, CM-272 treatment was able to restore the overlook the limitations of our in vitro studies, as other cell DNA methylation levels at these loci to resemble that of healthy types of the BM microenvironment may also play important mice, which was concomitant with the reduced tumor burden as roles in perpetuating the methylome alterations observed well as bone loss recovery observed in these mice (Fig. 6I). in MSCs. Although the biology of MBD is relatively well described, there Discussion is still a lack of pharmacological treatments to improve bone loss. The pathogenic transition from premalignant stages to active Clinically approved bone-modifying agents for the treatment of MM is complex and not well understood. One example of this MBD include bisphosphonates56, which inhibit bone resorption complexity is that although all MM cases emerge from the pre- by suppressing OC activity, and denosumab57, a monoclonal existing asymptomatic MGUS/SMM stage, not all MGUS antibody against the osteoclastogenic cytokine RANKL. However, progress into MM and may exist as a stable and independent these drugs only target the OC compartment, and bone disease disease. Nevertheless, despite being an asymptomatic stage, persists due to the absence of bone formation. Thus, therapeutic transformed plasma cells in MGUS present cytogenetic agents targeting OBs are needed. In this study, we demonstrated a alterations similar to that of myeloma plasma cells, as well as strategy for treating MBD by targeting aberrant DNA methyla- significant abnormalities in bone remodeling49,50. This indi- tion in MSCs. Firstly, we showed that co-culture of healthy MSCs cates that both genetic and microenvironmental alterations with MM cell lines yielded epigenetic and transcriptional changes exist from the early stages of the disease. In our study, we similar to that observed for MSCs from myeloma patients, and show that epigenetic alterations in MSCs already occur in the treatment with CM-272 was able to at least partially reverse these early asymptomatic stages of both MGUS and SMM, and changes. Additionally, this agent promoted the ability of MSCs to although many alterations are shared between all stages, the differentiate into OBs. These in vitro effects on bone were mir- majority of DNA methylation changes are specific to each rored in a mouse model of disseminated MM. Of note, CM-272 stage. These results are in accordance with previous studies treatment not only prevented bone loss by bone-anabolic effects that indicate the existence of stage-specific epigenetic altera- but also showed anti-myeloma activity. This is in line with pre- tions during MM progression in malignant plasma cells51,52. vious reports showing that DNMTs are targets for the treatment This phenomenon could be explained by the expansion of sub- of MM58–60 and also for improving the osteogenic differentiation populations of MSCs during MM disease progression, which ability of MSCs61. Additionally, we cannot discard the possibility may favor tumor development and drug resistance, similar to that the observed effects on tumor growth inhibition may be a what was observed to occur in MM cells53,54. consequence of the impairment of the cross-talk between MSCs Deregulation of methylome in MM-MSCs mediates tran- and MM cells. Moreover, the dual targeting effects of CM-272 scriptional and phenotypical alterations. Interestingly, many also inhibit the dimethylation of H3K9, which has been described genes of the Homeobox family displayed both epigenetic and to be crucial in the establishment of DNA methylation42,43. It is transcriptional dysregulation in patient MSCs, and these therefore rational to envision that the bone-anabolic effects changes were observed in earlier stages of the disease. In this mediated by CM-272, both in vitro and in vivo, involves the regard, members of the HOX family have been recently reversion of aberrant hypermethylation at Homeobox loci and described to be key drivers of OB differentiation, in which other OB-related genes in the MSC population. Nevertheless, it is their expression is fine-tuned by demethylation of their pro- possible that reduced tumor burden could be partially responsible moters during the osteogenic process55. Furthermore, we for restoring the bone-forming capacities of MM-MSCs. observed that healthy MSCs exposed to MM cells, similarly to In summary, our findings highlight the existence of aberrant that observed in patient MSCs, not only displayed an altered DNA methylation patterns in the BM-derived MSC population methylome but also showed impaired MSC-to-OB differ- which may impact myeloma progression and the development of entiation, as previously described20. We also observed that MBD. Moreover, our preclinical results support the idea that some of these methylome changes in MSCs occur in the therapeutic targeting of aberrant DNA methylation would result absence of direct cell–cell contact with MM cells, suggesting in an anti-myeloma effect and preservation of the appropriate the contribution of secretory mechanisms. Hence, our results osteogenic differentiation of MSCs to combat myeloma bone suggest that the impairment of osteogenesis in all stages of disease. 8 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE A B C CM-272 50 nM 11 150 5 HD-MSC Relative gene expression ** % MTT incorporation 10 * HD-MSC + MM cell line 4 MM-MSC Arbitrary units MM-MSC + MM cell line 100 9 3 8 2 * 8 * 50 ** * 7 1 vehicle 6 5 0 0 DNMT1 DNMT3A DNMT3B HOXA9 HOXC9 PITX1 RUNX2 HOXA10 HOXA6 0 25 50 100 250 nM CM-272 D HOXA6 HOXA9 HOXA10 HOXC9 PITX1 RUNX2 cg01210554 cg05139500 cg15864691 cg09045007 cg069335074 cg26366735 1.5 1.2 1.2 1.2 1.2 1.2 Relative methylation * 1.0 * 1.0 * 1.0 * 1.0 1.0 * 1.0 0.8 0.8 0.8 0.8 0.8 0.5 0.6 0.6 0.6 0.6 0.6 0.0 veh CM-272 veh CM-272 veh CM-272 veh CM-272 veh CM-272 veh CM-272 E HOXA6 HOXA9 HOXA10 HOXC9 PITX1 RUNX2 H3K9me2 relative enrichment 10 5 6 2.5 15 5 *** 2.0 * 4 H3K9me2 8 4 * 4 10 IgG 6 3 1.5 * 3 4 2 1.0 2 2 5 2 1 0.5 1 0 0 0 0.0 0 0 veh CM-272 veh CM-272 veh CM-272 veh CM-272 veh CM-272 veh CM-272 F G H CM-272 50 nM * % ALP activity (nM/min/μg) 150 10 Relative gene expression * vehicle * 100 * CM-272 25nM 5 * 50 CM-272 50 nM vehicle 0 0 vehicle CM-272 25nM CM-272 50 nM SPP1 IBSP BGLAP Methods 02156), the Clinical Ethics Committee for drug research in the Salamanca Health Participants. BM samples were obtained from the iliac crest of patients with newly Area (CEIC 73/07/2015), the Clinical Research Committee of the Bellvitge Uni- diagnosed MGUS (n = 10), SMM (n = 8), and MM (n = 9), according to the versity Hospital (ref. PR076/15) and the Research Ethics Committee of the Uni- International Myeloma Working Group criteria. BM samples from healthy controls versity of Navarra (ref. 2017.218). Clinical characteristics of MGUS, SMM, and (n = 8) were obtained from participants undergoing orthopedic surgery not related MM patients are listed in Supplementary Table 2. to oncology disease. Each sample was obtained after receiving the informed written consent of all participating subjects and following approval from the committees listed below for obtaining them and for the study protocol using them. The study Inhibitor. CM-272 (dual DNMTs and G9a inhibitor) was synthesized at the Center was approved by the Cancer Research Center–IBMCC Review Board (CICIC 2015/ for Applied Medical Research (University of Navarra)44,47,62. NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x Fig. 4 CM-272 treatment reactivates Homeobox gene expression and promotes the osteogenic differentiation of MSCs from MM patients. A Expression of DNA methyltransferases 1, 3A, and B (DNMT1, -3A, and -3B) in HD and MM-MSCs co-cultured for 24 h with the MM.1S cell line comparing with monoculture as assessed by the GeneChip Human Gene 1.0 ST Array. Box plots present mean ± SEM for 8 healthy donors and 9 MM patients. p-values were calculated by paired two-tailed Student t-test (*p < 0.05, **p < 0.01). B MSCs from MM patients were treated with the indicated doses of CM-272 for 72 h and subjected to MTT assay for viability. Mean and SEM are indicated on the line chart from 3 independent experiments. C Real-time RT- PCR was performed to determine the expression of hypermethylated Homeobox genes (HOX-A6, -A9, -A10, -C9, PITX1, RUNX2) in MM-MSCs treated with vehicle or 50 nM of CM-272 for 7 days. Box plots represent median ±IQR, with whiskers representing minimum and maximum, of 9 independent experiments. A paired two-tailed Student’s t-test was performed to calculate statistical significance (*p < 0.05, **p < 0.01). D DNA methylation analysis by pyrosequencing of selected CpGs located at the promoter regions of Homeobox genes in MM-MSCs treated with vehicle (gray) or CM-272 (blue) for 7 days. Bar plots present mean ± SEM for 3 independent experiments and paired two-tailed Student’s t-tests were performed (*p < 0.05). E ChIP assays showing the H3K9me2 (blue) enrichment at the promoter regions of Homeobox genes in MM-MSCs treated with vehicle or CM-272 for 7 days. IgG (gray) was used as a negative control. Data are shown as a relative enrichment of the bound fraction with respect to the input DNA. Bar plots present mean ± SEM for 3 independent experiments and paired two-tailed Student’s t-tests were performed (*p < 0.05, **p < 0.01, and ***p < 0.005). ALP activity was assessed in MM-MSCs (n = 3) cultured in osteogenic media in the presence of 25 nM (light blue) and 50 nM (dark blue) of CM-272, compared to vehicle (gray), by F p-NPP hydrolysis and G NBT-BCIP. H Expression of osteoblastogenesis markers IBSP (bone sialoprotein), BGLAP (osteocalcin) and SPP1 (osteopontin) was checked by qRT-PCR in MM-MSCs cultured in osteogenic media in the absence (vehicle) or presence of CM-272. For F and H data are shown as mean values ± SEM from three different experiments. Statistically significant tests (paired two-tailed Student’s t-tests) are represented as *p < 0.05, **p < 0.01, and ***p < 0.005 between vehicle and CM-272 condition. Cell lines. The human multiple myeloma cell line MM.1S was provided by Dr. Probes were annotated using IlluminaHumanMethylationEPICmanifest v0.3.067. Steven Rosen (Northwestern University, Chicago, IL), whereas RPMI-8226 cells Data quality was assessed using the standard pipeline from the minfi package. The were purchased from the American Type Culture Collection. The human data were quantile-normalized and chromosomes X and Y were removed to avoid mesenchymal stem cell (MSC) line immortalized by expression of the telomerase technical and biological bias. Furthermore, we discarded the DNA methylation reverse transcriptase gene (hMSC-TERT) was a generous gift from Dr. D Campana changes associated with the long-term culture of BM-MSCs based on the previous (Department of Pediatrics, Yong Loo Lin School of Medicine, National University studies68. M values (log2-transformed β-values) were utilized to obtain a p-value of Singapore, Singapore). Both cell lines were cultured in RPMI 1640 medium between sample groups by an eBayes-moderated paired t-test using the limma supplemented with 10% heat-inactivated fetal bovine serum (FBS), 100 U/ml package, in which age and sex were added in the interaction matrix. For the penicillin and 100 mg/ml streptomycin, and 1% L-glutamine. All the cell culture analysis of MSCs isolated from MM patients and healthy controls, we considered a media and reagents were purchased from Invitrogen (Paisley, UK). All cell types probe to be differentially methylated when the difference between the mean of β at were cultured at 37 °C in a humidified atmosphere in the presence of 5% disease versus control was over 10% (Δβ ≥ 0.1) and the statistical test was sig- CO2–95% air. nificant (**p < 0.01). In addition, we used the iEVORA algorithm69, provided by matrixTests v0.1.9 (https://CRAN.R-project.org/package=matrixTests), to desig- nate a probe as differentially variable. This algorithm detects the homogeneity of Bone marrow-derived MSC isolation and culture. MSCs were isolated and variances using the Bartlett’s test (FDR < 0.05) and then selects those probes whose characterized as described by Garayoa et al.17. Briefly, bone marrow aspirates were t-test is significant (*p < 0.05) in order to regularize the variability test which is obtained from the iliac crest and subjected to centrifugation on Ficoll-Paque (GE overly sensitive to single outliers. Healthcare, Uppsala, Sweden) to obtain mononuclear cells (BMMCs). BMMCs To evaluate the contribution of various covariates, including age and sex, we were plated and plastic-adherent cells were expanded until passage 3 (P3) in low- performed either a Pearson correlation or Wilcoxon signed-rank test depending on glucose DMEM supplemented with 10% heat-inactivated fetal bovine serum, 100 whether the covariate of interest was continuous or categorical. This is represented U/ml penicillin, 100 mg/ml streptomycin, and 1% L-glutamine. Selected MSCs from in Supplementary Fig. 4, in which a covariate with a p-value < 0.05 was considered both MM patients (n = 4) and healthy donors (n = 4) at P3 were tested to meet to significantly contribute to DNA methylation. minimal criteria as defined by the International Society for Cellular Therapy for For the direct co-culture of healthy MSCs with MM cell line, samples were multipotent mesenchymal stromal cells63. Specifically, MSCs were evaluated by normalized utilizing Noob and Quantile normalizations provided by minfi. The FACS for positive expression of CD73, CD90, CD105, CD44, and CD166 and paired analysis was performed and a probe was considered differentially negative staining for HLA-DR and hematopoietic markers (CD19, CD34, and methylated if Δβ was more than 10% and p-value was <0.01. CD45) (Supplementary Fig. 3A). In addition, the capability to differentiate into osteoblast, adipocyte, and chondrocyte was assessed (Supplementary Fig. 3B–D). Analyses and experiments were performed with MSCs at P3, with a maximum of 3 weeks at each passage. Gene ontology, motif, and chromatin state analysis. Functional annotation enrichment analysis was performed using GREAT tool v4.0.4 (http://great.stanford. edu/public/html)70 by mapping differentially methylated CpG site to the single DNA and RNA isolation and quantification. Genomic DNA was isolated by the nearest gene. CpGs annotated in the EPIC 850K array were used as background. proteinase K method or using the Maxwell® RSC Cell DNA Purification Kit GO categories with p-value of <0.01 were considered significantly enriched. (Promega) for samples containing low cell number. RNA was isolated using For TF binding motif analysis, HOMER motif discovery software v4.5 was Maxwell® RSC simplyRNA Cells Kit (Promega) according to the manufacturer’s used27, where a 500 bp-window upstream and downstream of the differentially instructions. DNA and RNA were quantified using Qubit® DNA Assay Kit methylated CpG sites was applied. CpGs annotated in the EPIC array were used as (Invitrogen) or NanoDrop ND-1000, respectively. background. To analyze chromatin states associated with DMPs, ChromHMM26 data sets DNA methylation and gene expression profiling using arrays. DNA samples from healthy donor MSCs were downloaded from the UCSC Genome Browser were bisulfite-converted using an EZ DNA methylation kit (Zymo Research, (https://genome.ucsc.edu/). Overlap was performed in R using the GenomicRanges Orange, CA) and hybridized onto an Infinium® MethylationEPIC BeadChip array package v1.42.071, where CpGs annotated in the EPIC array were used as (Illumina, Inc.). The array platform allows the assessment of DNA methylation background. status at >850,000 CpG sites at single-nucleotide resolution and covers 99% of RefSeq genes and 95% of CpG islands with an average number of six probes per island. Gene expression array normalization and analysis. Data processing and nor- RNA samples were obtained from healthy donors (n = 8), MGUS (n = 10), malization were carried out using the R statistical language. Background correction SMM (n = 10), and MM patients (n = 24) at diagnosis and 100 ng of excellent was performed using Robust Microarray Analysis (RMA) normalization provided quality RNA (RIN > 9) was hybridized onto a GeneChip Human Gene 1.0 ST by oligo package v1.54.172 and probes were annotated utilizing the hugene10st- (Affymetrix). transcriptcluster.db R package v8.7.073. Average expression was calculated for probes mapping to the same gene. For comparisons between groups, eBayes- Quality control, data normalization, and detection of differentially methylated moderated t-test provided by the limma R package65 was applied, where a p-value and variable CpGs. Methylation array data were processed in the statistical lan- < 0.05 was considered statistically significant. DMPs were mapped to the nearest guage R v4.0 in RStudio 1.3 (https://rstudio.com) using methods from the Bio- gene utilizing the GREAT online tool, and overlap with differentially expressed conductor libraries minfi (v1.36.0), lumi (v2.42.0), and limma (v3.46.0)64–66. genes were performed by overlapping gene names. 10 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE A RUNX2 SPP1 IBSP HOXB7 2.5 3 4 2.0 C Relative gene expression MSC 2.0 3 1.5 MSC + MM.1S MSC 2 1.5 ** * ** MSC + RPMI-8226 OB * * 2 1.0 OB + MM.1S 1.0 * 1 OB + RPMI-8226 1 0.5 0.5 2.0 0.0 0 0 0.0 CM-272 - - + - - + - - + - - + 1.5 mM Alizarin Red 2.5 3 5 2.5 Relative gene expression 1.0 2.0 4 2.0 * 2 1.5 * 3 1.5 0.5 1.0 ** * 2 1.0 1 0.0 0.5 1 0.5 CM-272 - - - + - + 0.0 0 0 0.0 CM-272 - - + - - + - - + - - + B RUNX2 SPP1 IBSP HOXB7 cg26366735 cg02549628 cg24986876 cg07547765 10 4 1.5 ** MSC 8 1.0 3 * DNA methylation MSC + MM.1S 1.0 6 MSC + RPMI-8226 2 0.5 4 0.5 1 2 0.0 0.0 0 0 CM-272 - - + - - + - - + - - + 8 1.5 1.5 1.0 6 * DNA methylation 1.0 1.0 4 0.5 0.5 0.5 2 0.0 0.0 0 0.0 CM-272 - - + - - + - - + - - + Fig. 5 CM-272 restores the gene expression levels and suppression of mineralization in MSCs from healthy donors exposed to MM cells. HD-MSCs (MSC; blue) were co-cultured with MM.1S (orange) or RPMI-8226 (red) cell lines separated by a transwell system in the presence (darker shade) or absence (lighter shade) of CM-272. A Relative expression and B DNA methylation of genes RUNX2, SPP1, IBSP, and HOXB7 was assessed. Gene expression was by normalization against RPL38. Bar plots represent mean ± SEM of 3–5 independent experiments and statistical significance was calculated by paired one-tailed Student’s t-test (*p < 0.05, **p < 0.01). C Mineralization was assessed by alizarin red staining in differentiated OBs from HD-MSCs co-cultured with MM.1S or RPMI-8226 cell lines in the presence or absence of CM-272. Representative micrographs show matrix mineralization by alizarin red staining of corresponding differentiated OBs. Data are represented as the mean ± SEM from three independent experiments and paired two-tailed Student’s t-tests were performed (*p < 0.05, **p < 0.01). Bisulfite pyrosequencing. For total DNA extraction, cells were lysed using lysis protein L38 (RPL38) used for normalization. Relative quantification of the target buffer (50 mM Tris pH 8.8, 10 mM EDTA pH 8.3, 100 mM NaCl, 1% SDS) in the gene expression was calculated by the comparative threshold cycle (Ct) method. presence of Proteinase K (Roche). Repeated centrifugation was performed to separate nucleic acids from lipids, in which DNA was subsequently precipitated using isopropanol and washed with 75 % ethanol. 100-300 ng of isolated DNA were Co-culture system and MSC sorting. MSCs from healthy donors at passage 3 (8 × bisulfite (BS)-converted using EZ DNA Methylation-Gold™ Kit (Zymo Research, 103 cells/cm2) or the hMSC-TERT cell line (10 × 103 cells/cm2) were first cultured CA, USA) according to manufacturers’ instructions. BS-converted DNA (~10 ng) in 100 mm culture dishes until they reached ~85% confluency, and then MM.1S was used as a template for amplification by conventional PCR using IMMOLA- cells (1:3 MSC:MM.1S ratio) were added in RPMI 1640 medium supplemented SETM DNA Polymerase kit (Bioline, London, UK). PCR primers were designed with 10% FBS and antibiotics. MM cells were changed twice a week until day 14 with the PyroMark Assay Design v2.0.2 software (Qiagen). PCR products were when MSCs were recovered by trypsinization and flow cytometry-based sorting of pyrosequenced with the PyromarkTM Q24 system (Qiagen), according to the CD13+ cells (BD Biosciences). For transwell experiments, MSCs from healthy manufacturer’s protocol. donors were seeded on bottom chambers and MM.1S and RPMI-8226 cells were seeded on PET membrane inserts containing 1 μm size pores to allow an exchange of soluble molecules. Transwell experiments were performed as in direct co- Quantitative real-time PCR (qRT-PCR). Reverse transcription was performed cultures. For OB differentiation studies, MSCs from healthy donors were cultured using the Transcriptor First Strand cDNA Synthesis Kit (Roche) according to the in an osteogenic medium supplemented with 20% of conditioned media from the manufacturer’s instructions. qRT-PCR was performed using LightCycler® 480 II MM.1S cell line. This medium was changed twice a week until day 10 (ALP System with LightCycler® 480 SYBR Green Mix and data were analyzed with activity) or day 20 (OB mineralization). For the isolation of mouse MSCs, cells were LightCycler® 480 II Software, version 1.5, all provided by Roche. Reactions were stained with a combination of CD45-PE, Ter-119-PE, Sca-1-FITC, and PDFGRα- performed in triplicate for both the target and the housekeeping gene ribosomal APC (BD Biosciences) as previously reported74 (Supplementary Fig. 5). Cell sorting NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x A B 2.5 Total flux (x 10 photons/sec) vehicle 2000000 2.0 CM-272 * Human Igλ (ng/ml) veh 1500000 1.5 10 1000000 1.0 * CM-272 500000 * 0.5 * ** 0 72 e 0.0 cl -2 hi 0 3 6 9 15 18 21 25 M ve Radiance 20,000 40,000 60,000 C (photons/sec/cm2/sr) days of treatment C D transversal plane sagittal plane vehicle CM-272 transversal plane vehicle sagittal plane CM-272 E F Trabecular Trabecular Trabecular Trabecular Serum P1NP volume occupancy connectivity separation 0.5 12 200 8 8 ** * * * Tb connectivity (mm) 0.4 150 Tb occupancy (%) Tb spacing (mm) Tb volume (mm ) 6 6 3 P1NP (ng/ml) 8 0.3 100 4 4 0.2 50 4 2 2 0.1 0 0.0 0 -50 0 0 72 e 72 e e 72 72 72 e e cl cl cl cl cl -2 -2 hi -2 -2 hi -2 hi hi hi M ve M ve ve M M M ve ve C C C C C 72 G H I y e th cl -2 al hi M he ve C **** Normalized DNA methylation HOXA7 **** **** HOXB4 GO annotation 1 cell fate commitment HOXB7 cell fate specification mesoderm development HOXB9 0 −1 0 1 anterior/posterior pattern specification Normalized 0 2 4 6 8 10 DNA methylation -1 HOXD10 -log (p value) 72 y e th cl -2 al hi M he ve C experiments were performed by the Flow Cytometry Core Facility at Germans ascorbic acid) and assayed as in Garcia-Gomez et al.75. Briefly, primary MSCs Trias i Pujol Research Institute utilizing FACSAria II cell sorter and analyzed using (P2–3) were cultured in osteogenic medium for the analysis of alkaline phosphatase BD FACSDiva version 6.1.1 (BD Biosciences, San Jose CA). (ALP) activity, expression of osteogenic markers (day 10), and formation of mineralized-nodules formation (day 20). ALP activity was determined by hydro- lysis of p-nitrophenylphosphate (Sigma-Aldrich) into p-nitrophenol and NBT/ OB differentiation assays. OBs were generated from mesenchymal precursors by BCIP substrates (Roche), whereas mineralization was assessed by quantitative culture in osteogenic medium (containing 5 mM β-glycerophosphate and 50 mg/ml measurement of Alizarin Red (Sigma-Aldrich) staining and absorbance was 12 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE Fig. 6 CM-272 prevents tumor-associated bone loss besides reducing multiple myeloma tumor burden. RPMI-8226-luc cells (8 × 106) were intravenously injected into NSG mice. After 4 weeks, mice were randomized into 2 groups [receiving vehicle (gray) and CM-272 (red); n = 6/group] and treated for additional 4 weeks with dosing and regimen schedules as specified in Supplementary Methods. Tumor dissemination was checked by A bioluminescence measurement and B serum levels of human Igλ secreted by RPMI-8226-luc cells at specified time points. Line plots represent mean and SEM. Box plots represent median ±IQR, with whiskers representing the minimum and maximum. Statistical significance was determined utilizing paired two-tailed Student’s t-test (*p < 0.05, **p < 0.01). C Representative microCT cross-sections at the metaphyses of distal femurs in a vehicle and CM-272- treated mice in transversal (upper) and sagittal (down) planes. D Transversal (left) and sagittal (right) planes of corresponding 3D renderings from microCT images at distal femurs (trabecular bone in red, cortical bone in gray). E Trabecular bone morphometric parameters from microCT images were quantitated for trabecular bone volume, occupancy, connectivity, and separation. F Serum levels of the bone formation marker P1NP were quantified by ELISA. Graphs represent mean values ± SEM with whiskers representing minimum and maximum values. CM-272-treated (red) mice were compared to the vehicle group (gray), where *p < 0.05; **p < 0.01 versus the vehicle control group. G Box plots showing DNA methylation levels of pooled MSCs obtained from healthy, vehicle- and CM-272-treated animals corresponding to hypermethylated CpGs between healthy and tumor-bearing animals. H GO enrichment analysis of CpG sites undergoing DNA hypermethylation changes in vehicle-treated MSCs versus MSCs from healthy mice. I Heatmap showing normalized DNA methylation levels of individual CpGs at selected Homeobox loci among animal groups. Data pooled from mice (n = 6) for each group with sufficient RRBS coverage (≥5 valid sequencing reads per CpG). <0.01 by paired two-tailed Student’s t-test. G Box plots showing mean ± SEM, with whiskers representing minimum and maximum values, of DNA methylation levels of pooled MSCs obtained from healthy (blue), vehicle- (gray), and CM- 272-treated (red) animals corresponding to hypermethylated CpGs between healthy and tumor-bearing animals. A paired two-tailed Student’s t-test was performed to calculate statistical significance (****p < 0.001). H GO enrichment analysis of CpG sites undergoing DNA hypermethylation changes in vehicle-treated MSCs versus MSCs from healthy mice. p-values were calculated utilizing a binomial test. I Heatmap showing normalized DNA methylation levels of individual CpGs at selected Homeobox loci among animal groups. The color scale ranges from white to red, representing low to high levels of DNA methylation. Data pooled from mice (n = 6) for each group with sufficient RRBS coverage (≥5 valid sequencing reads per CpG). measured by Multiskan Sky Microplate Spectrophotometer via SkanIt PC software images were acquired using a Quantum-GX microCT (Perkin Elmer) with the (ThermoFisher). following parameters: 80 kVp X-ray source voltage, 120 μA current, and the high- resolution scan protocol for a total acquisition time of 14 min and a gantry rotation of 360 degrees. The tomographic three-dimensional images containing the entire MTT assay. MSCs were seeded in 96 well plates and treated with increasing bone yielded a total of 512 slices, with isotropic 50 microns voxel size and a concentrations of CM-272. MTT was added at a final concentration of 0.5 mg/ml resolution of 512 × 512 pixels per slice. To perform the bone histomorphometry and incubated at 37 °C for 1 h. Cells were then washed with PBS and incubated in analysis a (10 × 10 × 10 mm) ROI containing the bone metaphysis was defined and the dark for 10 min in the presence of dimethyl sulphoxide. Absorbance at 570 nm subsequently reconstructed from the original scan at a resolution of 20 microns per was measured utilizing the Multiskan Sky Microplate Spectrophotometer. voxel using the Quantum 3.0 software. Analysis of trabecular microarchitecture in the distal femur was carried out Flow cytometry antibodies. Prior to analysis and sorting by flow cytometry, using ImageJ v1.8.076. First of all, cortical and trabecular bones were segmented distinct amounts of MSCs were stained with fluorochrome-conjugated antibodies. from the CT volume. To this end, the following steps were followed: (i) Antibody specifications and concentrations used were the following: CD13-PE (BD segmentation of the entire bone volume by thresholding the original volume to Biosciences, 347406), CD45-PE (eBioscience, 12-0451-81), Ter-119-PE obtain a 3D binary mask; (ii) segmentation of empty volumes inside the cortical (eBioscience, 12-5921-81), Sca-1-FITC (eBioscience 11-5981-81), PDFGRα-APC volume (trabecular-free zones) using logical operators over filled vs. unfilled (BD Biosciences, 562777), CD44-FITC (BD Biosciences, 347943), CD19-PerCP versions of the result of step i; (iii) segmentation of the interior volume of the (BD Biosciences, 332780), CD90-FITC (BD Biosciences, 555595), HLA-DR-PerCP cortical bone by applying 10 morphological dilations followed by 10 morphological (BD Biosciences, 347402), CD14-FITC (BD Biosciences, 345784), CD166-PE (BD erosions to the 3D mask obtained in step ii; (iv) segmentation of the cortical bone Biosciences, 559263), Cd45-PercCPcy5.5 (BD Biosciences, 332784), CD34-FITC by performing an XOR logical operation between the masks obtained in steps i and (Invitrogen, 11-0349-42), CD73-PE (BD Biosciences, 550257), CD105-APC (R&D iii; and finally, (v) segmentation of the trabecular bone by performing an AND System, FAB10971A). logical operation between the masks obtained in steps i and iii. The final cortical and trabecular bone segmentations were further refined by applying a median filter to remove noise in the respective 3D masks. Cortical and trabecular bone volumes Chromatin immunoprecipitation (ChIP)-quantitative PCR. MSCs (15 × 103 cells were then calculated by applying the segmentation masks on the original volume. per IP) were cross-linked with 1% formaldehyde for 15 min and subjected to From the obtained trabecular masks, histomorphometry parameters were chromatin immunoprecipitation after sonication. ChIP-qPCR assays were per- calculated using the BoneJ plugin (version 1.4.2)77. Finally, bone 3D reconstruction formed using LowCell ChIP kit™ protein A (Diagenode) and the antibody (5 μg) and visualization were performed using Amira 5.2 software (ThermoFisher against H3K9me2 (H3K9me2 Abcam ChIP-grade, clone:mAbcam 1220, Ref: Scientific). ab1220, Lot:GR45436-1). Data are represented as the ratio of the bound fraction over the input for each histone modification or factor. IgG was used as a negative control. Primer sequences were designed as close as possible from the CpG ELISA. Serum levels of human Igλ (indicating tumor burden) and N-terminal undergoing methylation changes. Primer sequences are shown in Supplementary propeptide of type I procollagen (P1NP) (indicating bone formation) were mea- Table 3. These experiments were performed with three biological replicates of each sured in mice sera using the Human Lambda ELISA kit (Bethyl Laboratories, origin. Texas, USA) and Rat/Mouse PINP EIA kit (Immunodiagnostic Systems, East Boldon, UK), respectively, following manufacturers’ instructions. Absorbance was In vivo model. Animal experiments were conducted according to relevant ethical measured using the Multiskan Sky Microplate Spectrophotometer. regulations for the use of laboratory animals and after acquired permission from the University of Salamanca Committee for animal experimentation (ref # 0000061). BALB/c-Rag2null IL2rγnull (BRG) mice // or NOD-scid IL2rγnull (NSG) Reduced representation bisulfite sequencing (RRBS). Sorted MSCs from three mice were bred and maintained in the SPF area of the University of Salamanca groups (healthy, vehicle-treated, and CM-272-treated mice) were pooled in order Animal Facility with controlled environment conditions (20–23 °C, 12:12 light/ to obtain a significant number of cells for performing the RRBS-seq. Isolated DNA dark cycles, 30–70% relative humidity) and fed ad libitum. CM-272 was solubilized from each pool was subjected to the RRBS pipeline as previously described78. In in 0.9% saline solution. RPMI-8226-luc cells (8 ×106) were injected intravenously brief, purified DNA was digested with MspI and subjected to bisulfite conversion. into 8-week-old NOD-SCID-IL-2Rγ−/− (NSG) mice (Charles River Laboratories) Following PCR amplification, RRBS libraries were generated from sequenced DNA and tumor development was monitored by noninvasive bioluminescence imaging following previously published procedures (http://code.google.com/p/bsmap/ (BLI) with a Xenogen IVIS 50 system (Caliper Life Sciences). After 4 weeks, ani- downloads/). Downstream normalization and analyses were performed following a mals were randomized into two groups (n = 6/group) receiving vehicle (0.9% saline previously published pipeline (http://rrbs-techdev.computational-epigenetics.org/). solution) or CM-272 (5 mg/kg, 5 times/week by intraperitoneal injection). CpG annotation and GO enrichment analysis was performed utilizing the GREAT online tool. Microcomputed tomography analysis. One femur of each animal was fixed in 10% formalin in order to preserve bone microarchitecture. 3D X-ray tomographic Primers. All primers used are listed in Supplementary Table 3. NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x Statistical analysis. Data are expressed as mean ± SEM and the n value for each 21. Goll, M. G. & Bestor, T. H. Eukaryotic cytosine methyltransferases. Annu. Rev. in vitro assay is specified in the corresponding figure legend. Statistical analyses Biochem. 74, 481–514 (2005). were carried out with Prism version 6.0 (GraphPad) and were performed using a 22. Cakouros, D. et al. Specific functions of TET1 and TET2 in regulating two-tailed Mann–Whitney U test or Student’s t-test. mesenchymal cell lineage determination. Epigenet. Chromatin 12, 3 (2019). 23. Garcia-Gomez, A., Rodríguez-Ubreva, J. & Ballestar, E. Epigenetic interplay Reporting summary. Further information on research design is available in the Nature between immune, stromal and cancer cells in the tumor microenvironment. Research Reporting Summary linked to this article. Clin. Immunol. 196, 64–71 (2018). 24. Teschendorff, A. E. et al. Epigenetic variability in cells of normal cytology is associated with the risk of future morphological transformation. Genome Med. Data availability 4, 24 (2012). The DNA methylation and expression data supporting the findings of this study have 25. Teschendorff, A. E. et al. The dynamics of DNA methylation covariation been deposited in the NCBI’s Gene Expression Omnibus database under the accession patterns in carcinogenesis. PLoS Comput. Biol. 10, e1003709 (2014). code GSE137419. This SuperSeries (GSE137419) is composed of the following SubSeries: 26. Ernst, J. & Kellis, M. ChromHMM: automating chromatin-state discovery and GSE137360 (methylation), GSE137369 (expression), and GSE137416 (methylation II). characterization. Nat. Methods 9, 215–216 (2012). All the other data supporting the findings of this study are available within the article and 27. Heinz, S. et al. Simple combinations of lineage-determining transcription its supplementary information files and from the corresponding author upon reasonable factors prime cis-regulatory elements required for macrophage and B cell request. Source data are provided with this paper. identities. Mol. Cell 38, 576–589 (2010). 28. Ducy, P., Zhang, R., Geoffroy, V., Ridall, A. L. & Karsenty, G. Osf2/Cbfa1: a transcriptional activator of osteoblast differentiation. Cell 89, 747–754 (1997). Received: 9 October 2019; Accepted: 14 December 2020; 29. Håkelien, A.-M. et al. The regulatory landscape of osteogenic differentiation. Stem Cells 32, 2780–2793 (2014). 30. Kanzler, B., Kuschert, S. J., Liu, Y. H. & Mallo, M. Hoxa-2 restricts the chondrogenic domain and inhibits bone formation during development of the branchial area. Development 125, 2587–2597 (1998). 31. Jeong, B.-C. ATF3 mediates the inhibitory action of TNF-α on osteoblast References differentiation through the JNK signaling pathway. Biochem. Biophys. Res. 1. San-Miguel, J. F. & Kantarjian, H. M. Multiple myeloma and chronic Commun. 499, 696–701 (2018). leukaemias in 2014: Improved understanding of disease biology and 32. Giuliani, N. et al. Myeloma cells block RUNX2/CBFA1 activity in human bone treatment. Nat. Rev. Clin. Oncol. 12, 71–72 (2014). marrow osteoblast progenitors and inhibit osteoblast formation and 2. Kumar, S. K. et al. Multiple myeloma. Nat. Rev. Dis. Primers 3, 17046 (2017). differentiation. Blood 106, 2472–2483 (2005). 3. Raje, N. & Roodman, G. D. Advances in the biology and treatment of bone 33. Seifert, A., Werheid, D. F., Knapp, S. M. & Tobiasch, E. Role of Hox genes in disease in multiple myeloma. Clin. Cancer Res. 17, 1278–1286 (2011). stem cell differentiation. World J. Stem Cells 7, 583 (2015). 4. Saad, F. et al. Pathologic fractures correlate with reduced survival in patients 34. Verlinden, L. et al. Nrp2 deficiency leads to trabecular bone loss and is with malignant bone disease. Cancer 110, 1860–1867 (2007). accompanied by enhanced osteoclast and reduced osteoblast numbers. Bone 5. Panaroni, C., Yee, A. J. & Raje, N. S. Myeloma and bone disease. Curr. 55, 465–475 (2013). Osteoporos. Rep. 15, 483–498 (2017). 35. Oshima, T. et al. Myeloma cells suppress bone formation by secreting a 6. Dhodapkar, M. V. MGUS to myeloma: a mysterious gammopathy of soluble Wnt inhibitor, sFRP-2. Blood 106, 3160–3165 (2005). underexplored significance. Blood 128, 2599–2606 (2016). 36. Kaneshiro, S. et al. Bruton tyrosine kinase (Btk) suppresses osteoblastic 7. Sundararajan, S., Kumar, A., Korde, N. & Agarwal, A. Smoldering multiple differentiation. J. Bone Miner. Metab. 33, 486–495 (2015). myeloma: emerging concepts and therapeutics. Curr. Hematol. Malig. Rep. 11, 37. Pizard, A. et al. Connexin 40, a target of transcription factor Tbx5, patterns 102–110 (2016). wrist, digits, and sternum. Mol. Cell. Biol. 25, 5073–5083 (2005). 8. Amodio, N., D’Aquila, P., Passarino, G., Tassone, P. & Bellizzi, D. Epigenetic 38. Lanctôt, C., Moreau, A., Chamberland, M., Tremblay, M. L. & Droui, J. modifications in Multiple Myeloma: recent advances on the role of DNA and Hindlimb patterning and mandible development require the Ptx1 gene. histone methylation. Expert Opin. Ther. Targets 21, 91–101 (2016). Development 126, 1805–1810 (1999). 9. Das, R. et al. Microenvironment-dependent growth of preneoplastic and 39. Wei, X.-F., Chen, Q.-L., Fu, Y. & Zhang, Q.-K. Wnt and BMP signaling malignant plasma cells in humanized mice. Nat. Med. 22, 1351–1357 (2016). pathways co-operatively induce the differentiation of multiple myeloma 10. Ghobrial, I. M., Detappe, A., Anderson, K. C. & Steensma, D. P. The bone- mesenchymal stem cells into osteoblasts by upregulating EMX2. J. Cell. marrow niche in MDS and MGUS: implications for AML and MM. Nat. Rev. Biochem. 120, 6515–6527 (2018). Clin. Oncol. https://doi.org/10.1038/nrclinonc.2017.197 (2018). 40. Iyyanar, P. P. R. & Nazarali, A. J. Hoxa2 inhibits bone morphogenetic protein 11. Garcia-Gomez, A., Sanchez-Guijo, F., Del Cañizo, M. C., San Miguel, J. F. & signaling during osteogenic differentiation of the palatal mesenchyme. Front. Garayoa, M. Multiple myeloma mesenchymal stromal cells: contribution to Physiol. 8, 929 (2017). myeloma bone disease and therapeutics. World J. Stem Cells 6, 322–343 41. Garcia-gomez, A. et al. Transcriptomic profile induced in bone marrow (2014). mesenchymal stromal cells after interaction with multiple myeloma cells: 12. Friedenstein, A. J., Petrakova, K. V., Kurolesova, A. I. & Frolova, G. P. implications in myeloma progression and myeloma bone disease. Oncotarget Heterotopic of bone marrow. Analysis of precursor cells for osteogenic and 5, 8284–8305 (2014). hematopoietic tissues. Transplantation 6, 230–247 (1968). 42. Estève, P.-O. et al. Direct interaction between DNMT1 and G9a coordinates 13. Corre, J. et al. Bone marrow mesenchymal stem cells are abnormal in multiple DNA and histone methylation during replication. Genes Dev. 20, 3089–3103 myeloma. Leukemia 21, 1079–1088 (2007). (2006). 14. Xu, S. et al. Impaired osteogenic differentiation of mesenchymal stem cells 43. Tachibana, M., Matsumura, Y., Fukuda, M., Kimura, H. & Shinkai, Y. G9a/ derived from multiple myeloma patients is associated with a blockade in the GLP complexes independently mediate H3K9 and DNA methylation to deactivation of the Notch signaling pathway. Leukemia 26, 2546–2549 (2012). silence transcription. EMBO J. 27, 2681–2690 (2008). 15. Reagan, M. R. & Ghobrial, I. M. Multiple myeloma mesenchymal stem cells: 44. San José-Enériz, E. et al. Discovery of first-in-class reversible dual small characterization, origin, and tumor-promoting effects. Clin. Cancer Res. 18, molecule inhibitors against G9a and DNMTs in hematological malignancies. 342–349 (2012). Nat. Commun. 8, 15424 (2017). 16. Rajkumar, S. V. Multiple myeloma: 2016 update on diagnosis, risk- 45. Bárcena-Varela, M. et al. Dual targeting of histone methyltransferase G9a and stratification, and management. Am. J. Hematol. 91, 719–734 (2016). DNA-methyltransferase 1 for the treatment of experimental hepatocellular 17. Garayoa, M. et al. Mesenchymal stem cells from multiple myeloma patients carcinoma. Hepatology 69, 587–603 (2019). display distinct genomic profile as compared with those from normal donors. 46. Segovia, C. et al. Inhibition of a G9a/DNMT network triggers immune- Leukemia 23, 1515–1527 (2009). mediated bladder cancer regression. Nat. Med. 25, 1073–1081 (2019). 18. Giuliani, N. et al. Bone osteoblastic and mesenchymal stromal cells lack 47. Rabal, O. et al. Discovery of reversible DNA methyltransferase and lysine primarily tumoral features in multiple myeloma patients. Leukemia 24, methyltransferase G9a inhibitors with antitumoral in vivo efficacy. J. Med. 1368–1370 (2010). Chem. 61, 6518–6545 (2018). 19. Todoerti, K. et al. Distinct transcriptional profiles characterize bone 48. Xu, S., De Veirman, K., De Becker, A., Vanderkerken, K. & Van Riet, I. microenvironment mesenchymal cells rather than osteoblasts in relationship Mesenchymal stem cells in multiple myeloma: a therapeutical tool or target? with multiple myeloma bone disease. Exp. Hematol. 38, 141–153 (2010). Leukemia 32, 1500–1514 (2018). 20. Adamik, J. et al. EZH2 or HDAC1 inhibition reverses multiple myeloma- 49. Kristinsson, S. Y. et al. Monoclonal gammopathy of undetermined significance induced epigenetic suppression of osteoblast differentiation. Mol. Cancer Res. and risk of skeletal fractures: a population-based study. Blood 116, 2651–2655 15, 405–417 (2017). (2010). 14 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20715-x ARTICLE 50. Capp, J.-P. & Bataille, R. Multiple myeloma exemplifies a model of cancer Acknowledgements based on tissue disruption as the initiator event. Front. Oncol. 8, 355 (2018). We thank CERCA Program/Generalitat de Catalunya and the Josep Carreras Foundation 51. Agirre, X. et al. Whole-epigenome analysis in multiple myeloma reveals DNA for institutional support. E.B. was funded by the Spanish Ministry of Science and hypermethylation of B cell-specific enhancers. Genome Res. 25, 478–487 Innovation (grant numbers SAF2014-55942-R and SAF2017-88086-R), co-funded by (2015). FEDER funds/European Regional Development Fund (ERDF)—a way to build Europe, 52. Heuck, C. J. et al. Myeloma is characterized by stage-specific alterations in and a Senior Research Award from the Multiple Myeloma Research Foundation DNA methylation that occur early during myelomagenesis. J. Immunol. 190, (MMRF). C.O.-d.-S. was funded by the Spanish Ministry of Science, Innovation and 2966–2975 (2013). Universities, under grant RTI2018-094494-B-C22 (MCIU/AEI/FEDER, UE). M.G. 53. Ledergor, G. et al. Single cell dissection of plasma cell heterogeneity in received financial support from the Spanish FIS-ISCIII (PI15/02156 and PI19/01384) and symptomatic and asymptomatic myeloma. Nat. Med. 24, 1867–1876 (2018). FEDER. A.G.G is funded by a postdoctoral contract of the Asociación Española Contra el 54. van Nieuwenhuijzen, N., Spaan, I., Raymakers, R. & Peperzak, V. From MGUS Cáncer (AECC). F.P. was funded by grants from Instituto de Salud Carlos III (ISCIII), to multiple myeloma, a paradigm for clonal evolution of premalignant cells. PI17/00701 and PI19/01352, TRASCAN (EPICA and Immunocell), Fundació La Marató Cancer Res. 78, 2449–2456 (2018). de TV3, the Accelerator award CRUK/AIRC/AECC joint funder-partnership, CIBER- 55. da Silva, R. A. et al. HOXA cluster gene expression during osteoblast ONC (CB16/12/00489) and co-financed with FEDER funds and Fundación Ramón differentiation involves epigenetic control. Bone 125, 74–86 (2019). Areces (PREMAMM). 56. Terpos, E., Ntanasis-Stathopoulos, I. & Dimopoulos, M. A. Myeloma bone disease: from biology findings to treatment approaches. Blood 133, 1534–1539 (2019). Ethical declaration We have complied with all relevant ethical regulations for work with human participants 57. Raje, N. et al. Denosumab versus zoledronic acid in bone disease treatment of in which informed written consent was obtained and the collection of samples and the newly diagnosed multiple myeloma: an international, double-blind, double- study protocol for the use of those samples was approved by the Cancer Research dummy, randomised, controlled, phase 3 study. Lancet Oncol. 19, 370–381 Center–IBMCC Review Board (CICIC 2015/02156), the Clinical Ethics Committee for (2018). drug research in the Salamanca Health Area (CEIC 73/07/2015), the Clinical Research 58. Zhou, W., Chen, H., Hong, X., Niu, X. & Lu, Q. Knockdown of DNA Committee of the Bellvitge University Hospital (ref. PR076/15) and the Research Ethics methyltransferase-1 inhibits proliferation and derepresses tumor suppressor Committee at Universidad de Navarra (code 2017.218). Bone marrow samples were genes in myeloma cells. Oncol. Lett. 8, 2130–2134 (2014). obtained conforming to ethical standards and good clinical practice for research work 59. Das, D. S. et al. A novel hypoxia-selective epigenetic agent RRx-001 triggers with principles expressed in the Declaration of Helsinki. All animal work was conducted apoptosis and overcomes drug resistance in multiple myeloma cells. Leukemia according to relevant national and international guidelines for animal research and 30, 2187–2197 (2016). approved by the Bioethics Committee of the University of Salamanca (ref # 0000061). 60. Harada, T. et al. HDAC3 regulates DNMT1 expression in multiple myeloma: therapeutic implications. Leukemia 31, 2670–2677 (2017). 61. Yan, X. et al. 5-Azacytidine improves the osteogenic differentiation potential Author contributions of aged human adipose-derived mesenchymal stem cells by DNA A.G.-G., T.L., J.R.-U., and E.B. conceived and designed experiments; A.G.-G., T.L., C.C.-F., J. demethylation. PLoS ONE 9, e90846 (2014). R.-U., L.C., G.G.-T., and L.S.-S. performed experiments; A.G.-G., T.L., and F.C.-M. performed 62. Rabal, O. et al. Detailed exploration around 4-aminoquinolines chemical space biocomputing analysis; M.M., L.S.-S., S.M., X.M., J.O., C.O.-d.-S., E.S.J.-E., M.E., X.A., F.P., and to navigate the lysine methyltransferase G9a and DNA methyltransferase M.G. participated in the data acquisition (performed patient selection, provided drugs, biological spaces. J. Med. Chem. 61, 6546–6573 (2018). samples, animals, and facilities); A.G.-G., T.L., J.R.-U., and E.B. analyzed and interpreted the 63. Dominici, M. et al. Minimal criteria for defining multipotent mesenchymal data; A.G.-G., T.L., and E.B. wrote the paper. All authors read and approved the final stromal cells. The International Society for Cellular Therapy position manuscript. statement. Cytotherapy 8, 315–317 (2006). 64. Aryee, M. J. et al. Minfi: A flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, Competing interests 1363–1369 (2014). The authors declare no competing interests. 65. Ritchie, M. E. et al. Limma powers differential expression analyses for RNA- sequencing and microarray studies. Nucleic Acids Res. 43, e47 (2015). Additional information 66. Du, P., Kibbe, W. A. & Lin, S. M. lumi: A pipeline for processing Illumina Supplementary information is available for this paper at https://doi.org/10.1038/s41467- microarray. Bioinformatics 24, 1547–1548 (2008). 020-20715-x. 67. Hansen, K. D. IlluminaHumanMethylationEPICmanifest: manifest for Illumina’s EPIC methylation arrays (2016). Correspondence and requests for materials should be addressed to A.G.-G. or E.B. 68. Koch, C. M. et al. Pluripotent stem cells escape from senescenceassociated DNA methylation changes. Genome Res. 23, 248–259 (2013). Peer review information Nature Communications thanks Stuart Rushworth and the 69. Teschendorff, A. E. et al. DNA methylation outliers in normal breast tissue other, anonymous, reviewer(s) for their contribution to the peer review of this work. identify field defects that are enriched in cancer. Nat. Commun. 7, 10478 (2016). Reprints and permission information is available at http://www.nature.com/reprints 70. McLean, C. Y. et al. GREAT improves functional interpretation of cis- regulatory regions. Nat. Biotechnol. 28, 495–501 (2010). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 71. Lawrence, M. et al. Software for computing and annotating genomic ranges. published maps and institutional affiliations. PLoS Comput. Biol. 9, e1003118 (2013). 72. Carvalho, B. S. & Irizarry, R. A. A framework for oligonucleotide microarray preprocessing. Bioinformatics 26, 2363–2367 (2010). 73. MacDonald, J. hugene10sttranscriptcluster.db: Affymetrix hugene10 Open Access This article is licensed under a Creative Commons annotation data (chip hugene10sttranscriptcluster) (2017). Attribution 4.0 International License, which permits use, sharing, 74. Houlihan, D. D. et al. Isolation of mouse mesenchymal stem cells on adaptation, distribution and reproduction in any medium or format, as long as you give the basis of expression of Sca-1 and PDGFR-α. Nat. Protoc. 7, 2103–2111 appropriate credit to the original author(s) and the source, provide a link to the Creative (2012). Commons license, and indicate if changes were made. The images or other third party 75. Garcia-Gomez, A. et al. Dasatinib as a bone-modifying agent: anabolic and material in this article are included in the article’s Creative Commons license, unless anti-resorptive effects. PLoS ONE 7, e34914 (2012). indicated otherwise in a credit line to the material. If material is not included in the 76. Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 article’s Creative Commons license and your intended use is not permitted by statutory years of image analysis. Nat. Methods 9, 671–675 (2012). regulation or exceeds the permitted use, you will need to obtain permission directly from 77. Doube, M. et al. BoneJ: free and extensible bone image analysis in ImageJ. the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Bone 47, 1076–1079 (2010). licenses/by/4.0/. 78. Gu, H. et al. Preparation of reduced representation bisulfite sequencing libraries for genome-scale DNA methylation profiling. Nat. Protoc. 6, 468–481 (2011). © The Author(s) 2021 NATURE COMMUNICATIONS | (2021)12:421 | https://doi.org/10.1038/s41467-020-20715-x | www.nature.com/naturecommunications 15