Identification and validation of a histone modification-related gene signature to predict the prognosis of multiple myeloma

Background: Multiple myeloma (MM) is an incurable plasma cell malignancy with high heterogeneity. Current staging systems, including the International Staging System (ISS) and Revised ISS (R-ISS), have limited prognostic accuracy. Given the role of histone modifications in MM progression, we developed a histone modification-related (HMR) prognostic model to improve MM risk stratification.

Methods: Gene expression and mutation data were downloaded from the Gene Expression Omnibus database and the Cancer Genome Atlas. Prognostic HMR-related genes were identified through a combination of univariate Cox regression, least absolute shrinkage and selection operator Cox regression, and random survival forest analysis. The genes identified were then used to construct the HMR risk score model using multivariate Cox regression. The model was validated using Kaplan-Meier survival, time-dependent receiver operating characteristic curves analysis. A nomogram combining the HMR score with clinical features was developed. Functional enrichment, immune infiltration, somatic mutation, and drug sensitivity analysis were conducted to explore the biological relevance of the model.

Results: Seven HMR genes with prognostic significance were identified. The HMR risk score stratified patients into high-risk and low-risk groups, with significant survival differences. The model demonstrated favorable predictive performance, and was shown to be an independent prognostic factor. The nomogram showed good calibration and discriminative ability, offering a practical tool for individual patient risk assessment. Functional analysis revealed that the HMR risk score is associated with dysregulated cell cycle progression, proliferation, and immunosuppression in MM, which may contribute to disease progression and drug resistance. Moreover, drug sensitivity analysis indicated potential associations between the HMR score and response to specific therapeutic agents, highlighting its potential role in guiding personalized treatment.

Conclusion: We developed an HMR gene signature that has potential for prognostic prediction and may help guide personalized treatment strategies in MM.

Introduction

Multiple myeloma (MM) is a malignant plasma cell disease that accounts for about 10% of all hematological cancers, and is the second most common hematologic cancer worldwide (Siegel et al., 2022; Zhuge et al., 2025). Over the last 30 years, MM’s worldwide incidence and mortality rates have more than doubled (Zhuge et al., 2025). The clinical manifestations include renal impairment, hypercalcemia, lytic bone lesions, and anemia (Kyle and Rajkumar, 2008). Despite the continuous emergence of novel therapeutic agents and due to the heterogeneity in pathogenesis, the prognosis for MM patients remains poor (van de Donk et al., 2021; Attal et al., 2017; Karimi et al., 2025). The widely used International Staging System (ISS) is based on albumin (ALB) and β2-microglobulin (B2M) (Girija Sivasankaran et al., 2025). However, it lacks genetic and molecular markers, limiting its accuracy in guiding individualized treatment and prognosis assessment (Hagen et al., 2022). To further improve the prognosis, it is crucial to identify and integrate additional prognostic biomarkers to refine risk stratification, improve outcome prediction, and optimize treatment selection.

Epigenetic modifications are reversible heritable changes in gene expression that do not alter the DNA sequence. Epigenetic variations are achieved through covalent chemical modifications of chromatin, including DNA and histone modifications, chromatin remodeling, and non-coding RNAs (Nagano and Fraser, 2009; Chervona and Costa, 2012; Lee and Kim, 2022). Among these, histone modifications, such as methylation, acetylation, phosphorylation, adenylylation, ubiquitination, and ADP-ribosylation occur at histone tails and play a crucial role in regulating gene expression (Kouzarides, 2007). These modifications regulate gene expression by altering chromatin structure or influencing transcription factor binding to gene promoters (Li et al., 2007). Histone modifications are closely associated with various biological processes, including gene expression regulation, cell differentiation and development, the cell cycle, tumorigenesis, immune cell infiltration, and cancer prognosis (Chervona and Costa, 2012; Ji et al., 2025; Sharma et al., 2010; Bi et al., 2020).

Histone modifications play a crucial role in regulating key gene expression and influencing disease progression in MM (Ismail et al., 2022). Abnormal expression of histone-modifying enzymes, such as histone deacetylases (HDACs) and histone methyltransferases (HMTs), disrupts transcriptional balance in MM cells, thereby affecting cell proliferation, apoptosis, and drug resistance. Targeting histone modifications has emerged as a promising therapeutic strategy, with HDAC inhibitors and EZH2 (Enhancer of Zeste Homolog 2) inhibitors showing significant potential in MM treatment (Nylund et al., 2021; Pu et al., 2024). Among the genetic abnormalities in MM, the t(4;14) translocation occurs in 15%–20% of patients, leading to FGFR3 (fibroblast growth factor receptor 3) and MMSET (multiple myeloma SET domain) overexpression. MMSET, a member of the NSD (Nuclear receptor-binding SET domain) family, possesses a SET domain responsible for histone methylation. Overexpression of MMSET in MM cells leads to globally elevated H3K36 dimethylation, accompanied by a significant reduction in H3K27 methylation. Additionally, MMSET binds to the IRF4 (interferon regulatory factor 4) promoter and regulates its expression, affecting cell growth, adhesion, and chromatin accessibility (Xie et al., 2015; Martinez-Garcia et al., 2011). Polycomb repressive complex 2 (PRC2) maintains the silent state of target genes such as HOX genes through trimethylation of H3K27, contributing to development and differentiation as well as tumorigenesis (Margueron and Reinberg, 2011). EZH2 is a core component of PRC2 that confers the HMT activity (Chase and Cross, 2011). Overexpression of EZH2 is linked to MM progression and poor prognosis (Pawlyn et al., 2017). EZH2 inhibitors reduce global H3K27me3 level and trigger apoptosis in MM cells (Hernando et al., 2016). HDAC inhibitors also impair MM cell growth and survival. Preclinical studies have demonstrated that HDAC inhibitors trigger apoptosis, as well as induce cell cycle arrest, in MM cells (Maiso et al., 2006).

In this study, we developed a histone modification-related (HMR) prognostic model comprising of seven genes to predict clinical outcomes in MM. Our findings not only elucidate the prognostic value of HMR genes but also establish a foundation for future research on targeted therapies strategies in MM.

Materials and methods

Data collection

Gene expression and clinical data of GSE24080, GSE136337, GSE136324, and GSE2658 were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The data were preprocessed with normalization using the R package “limma”. The RNA sequence data and related somatic mutation data were acquired from the MMRF-CoMMpass project on the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). Patients with complete survival data and an overall survival time greater than 1 month were included in this study. Patients who did not meet these criteria or had incomplete clinical data were excluded from the analysis.

Dataset GSE24080 was used as the training cohort for model construction, while GSE136337, GSE2658, and MMRF-CoMMpass were used as validation cohorts. The MMRF-CoMMpass dataset provided somatic mutation data. GSE136324, which contains RNA extracted from the whole bone marrow (WBM) of MM patients, was used for immune infiltration analysis. Detailed information on all datasets used in this study is provided in Supplementary Table S1.

This study strictly adheres to data usage regulations and conducts secondary analysis based on publicly available data, thus requiring no additional ethical approval.

Construction and validation of a HMR prognostic risk model

HMR genes were extracted from the “GOMF_HISTONE_MODIFYING_ACTIVITY” gene set in the Gene Set Enrichment Analysis (GSEA) database (http://www.gsea-msigdb.org/gsea/msigdb). After intersecting these genes with those detected in the GSE24080, GSE136337, GSE136324, GSE2658, and MMRF-CoMMpass, a total of 173 genes were selected for further analysis.

The GSE24080 dataset was used as the training cohort to create the HMR prognostic risk model. Univariate Cox regression was performed to identify genes related to prognosis (p value <0.01 and false discovery rates (FDR) < 0.05). The candidate genes identified were then subjected to two complementary feature selection methods to enhance robustness. First, Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression was performed using the “glmnet” R package (v 4.1-8) to minimize overfitting and select features with non-zero coefficients. Ten-fold cross-validation was conducted to determine the optimal lambda value, resulting in the selection of 13 genes. Second, Random Survival Forest (RSF) analysis was carried out using the randomForestSRC (v 3.3.3) R package with 10 trees (ntree = 10) and one split per node (nsplit = 1). Variable importance was evaluated using the minimal depth method, and the top-ranked genes were retained. The intersection of genes identified by both LASSO and RSF methods yield seven genes: SUZ12, KAT2A, AURKA, BUB1, UTY, SUV39H2, and PCGF5. Pairwise Spearman correlation analysis was performed among the seven candidate genes. To further evaluate multicollinearity, the variance inflation factor (VIF) was calculated for each gene.

Subsequently, these genes were included in a multivariate Cox proportional hazards regression model to construct the final prognostic signature. A prognostic risk score was calculated for each patient as a linear combination of expression levels weighted by multivariate Cox coefficients: HMR score $={\sum }_{i=1}^n{\beta }_i\times {\mathit{Exp}}_i$, where ${\beta }_i$ is the coefficient of gene i, and ${\mathit{Exp}}_i$ denotes its normalized expression level. The formula was applied using the “predict” function in the “survival” package (v 3.8.3), with coefficients derived from the training cohort.

The optimal cutoff value for stratifying patients into high-risk and low-risk groups was determined using survminer R package (v 0.5.0), which identifies the threshold that best separates survival outcomes based on the log-rank test. To ensure comparability across datasets, the quantile corresponding to the optimal cutoff in the training cohort was calculated and applied to the validation cohorts to define high-risk and low-risk groups.

Finally, the prognostic value of the HMR risk score was evaluated using Kaplan–Meier (KM) survival analysis and time-dependent receiver operating characteristic (ROC) curves, using “survminer”, “survival”, and “timeROC” R packages.

Establishment of predictive nomogram

Univariate and multivariate Cox regression analyses were performed on the HMR risk score and clinical characteristics to identify independent prognostic predictors. Variables with a p-value <0.05 were selected for further analysis. A combined model was constructed using the “rms” and “regplot” R packages to generate a nomogram. The predictive performance of the model was evaluated using calibration curves, the concordance index (C-index), and time-dependent ROC curves (using the “timeROC” package, v 0.4). The clinical benefit of the model was further assessed by decision curve analysis (DCA) using the “ggDCA” package (v 1.2).

Furthermore, the Wilcoxon rank-sum test was used to evaluate differences in HMR risk score across various clinical subgroups of MM patients.

Differential expression genes (DEGs) analysis and pathway enrichment analysis

For GSE24080, DEGs between high-risk and low-risk groups were identified using the “limma” package (v3.62.2), with thresholds of |logFC| > 0.58 and adjusted p-value <0.05. For the MMRF-CoMMpass cohort, DEGs was conducted using the “DESeq2” package (v1.46.0), applying a threshold of |log2FoldChange | > 1 and adjusted p-value <0.05. The identified DEGs were then subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses using the “clusterProfiler” package (v4.14.4). To further validate the biological relevance, Gene Set Enrichment Analysis (GSEA) was performed using the “fgsea” package (v1.32.2), based on the preranked gene list ordered by the signal-to-noise ratio derived from differential expression between high-risk and low-risk groups.

Evaluation of immune infiltration populations

We used a modified method described by Danziger et al. (2020) to calculate the immune cell infiltration levels in the bone marrow. Briefly, the whole bone marrow expression matrix was deconvoluted to calculate the fraction of 27 different cell types in WBM, including four types of plasma B cells. Then, the tumor immune microenvironment specific gene expression matrix was calculated according to the method described by Li et al. (2024).

Evaluation of somatic mutations

Samples from the MMRF-CoMMpass project were divided into high-risk and low-risk groups as described above. The distribution of gene mutations across different groups was then visualized using the “maftools” R package (v 2.22.0).

Drug sensitivity prediction

Drug sensitivity differences between the low-risk and high-risk groups were predicted using the “pRRophetic” R package (v 0.5) (Geeleher et al., 2014a). The dataset within the “pRRophetic” package is derived from the “cgp 2016” initiative, encompassing gene expression matrices and drug treatment information (Geeleher et al., 2014b). Wilcoxon test was used to identify drugs exhibiting different sensitivity between groups (p < 0.05).

Statistical analysis

Statistical analysis was performed with R software 4.4.2. Continuous variables were compared using the Student’s t-test or Wilcoxon rank-sum test, depending on the normality test. Categorical variables were analyzed using the Chi-square test or Fisher’s exact test. Prognostic gene selection and model construction were performed using univariate and multivariate Cox regression, LASSO Cox regression, and RSF analysis. Survival differences were assessed by KM survival analysis and log-rank tests. Model performance was evaluated using time-dependent ROC curves and C-index. Clinical usefulness was assessed by DCA. Correlation analyses were performed using Spearman’s rank correlation test. A two-sided p-value <0.05 was considered statistically significant.

Results

Patient selection and statistics

This study included a total of 2790 MM patients from five datasets, with statistical details provided in Supplementary Table S1. The GSE24080 dataset (n = 556) was used as the training cohort to construct a HMR prognostic risk model. The GSE136337 (n = 424), GSE2658 (n = 555), and MMRF-CoMMpass project (n = 853) datasets were selected as validation cohorts. Somatic mutation data were extracted from the MMRF-CoMMpass project dataset (n = 804). The overall design and workflow of this study are illustrated in (Figure 1). Additionally, the GSE136324 dataset (n = 402) was used for immune infiltration analysis.

Figure 1

Figure 1. Study design for HMR risk model development.

Development of an HMR prognostic risk model

To construct a HMR risk score, we identified seven HMR genes in the GSE24080 training cohort through univariate Cox regression, LASSO Cox regression analysis, and random forest analysis, as shown in (Figures 2A–C). Univariate Cox regression showed that SUZ12, KAT2A, AURKA, BUB1, and SUV39H2 were associated with poor prognosis (HR > 1, p < 0.05), indicating their potential role as high-risk genes. In contrast, UTY and PCGF5 were associated with favorable prognosis (HR < 1, p < 0.05), suggesting a protective effect (Figure 2D). Consistently, KM survival analysis based on median expression level further confirmed that high expression of SUZ12, KAT2A, AURKA, BUB1, and SUV39H2 correlated with worse prognosis, while high expression of UTY and PCGF5 suggested better prognosis (Supplementary Figure S1A–G).

Figure 2

Figure 2. Construction of HRM prognostic risk model. (A,B) LASSO Cox regression analysis for variable selection. (C) Venn diagram shows the intersection of genes identified by LASSO regression and Random Forest analysis. (D) Forest plot of univariate Cox regression analysis in the GSE24080 training cohort. (E) Forest plot of multivariate Cox regression analysis in the training cohort. (F) Kaplan–Meier survival curve of high-risk and low-risk in the training cohort. (G) The distribution of risk scores and survival status of MM in the training cohort. (H) Heatmap of the HRM related genes in the training cohort. (I) Time-dependent ROC analysis in the training cohort.

Prior to model construction, pairwise correlations among the seven candidate genes were assessed using Spearman’s rank correlation analysis. The correlation coefficients ranged from −0.17 to 0.51 (Supplementary Figure S2A), indicating no strong collinearity. To further evaluate multicollinearity, the variance inflation factor (VIF) was calculated for each gene, all of which were <5 (Supplementary Figure S2B), confirming the absence of significant multicollinearity.

Multivariate Cox regression revealed that AURKA and UTY are independent prognostic risk factors in MM patients (Figure 2E). Based on the regression coefficient calculated by Multivariate Cox regression (Supplementary Figure S2C), the HMR risk score was calculated using the following formula: HMR risk score = (0.2052429×SUZ12 expression) +(0.2399438×KAT2A expression) +(0.3476679×AURKA expression) +(0.0905674×BUB1 expression) −(0.2257278×UTY expression) +(0.2372008×SUV39H2 expression) −(0.2765753×PCGF5 expression). Using this formula, we calculated the HMR risk score for each individual in both the training and validation cohorts. Subsequently, patients in the training dataset were stratified into high-risk and low-risk groups based on the optimal cutoff value. For the validation dataset, the corresponding quantile derived from the training dataset cutoff was used to define the risk groups.

Validation of the HMR prognostic model

In the training dataset, KM survival analysis revealed that patients in the high-risk group had worse overall survival (Figure 2F). Risk distribution plot showed a higher proportion of deceased patients and shorter survival time in the high-risk group, further indicating a poorer prognosis (Figure 2G). The gene expression heatmap showed that high-risk genes (SUZ12, KAT2A, AURKA, BUB1, and SUV39H2) were upregulated in the high-risk group. Conversely, low-risk genes (UTY and PCGF5) exhibited higher expression in the low-risk group (Figure 2H). To evaluate the predictive performance of the HMR risk score model, we conducted a time-dependent ROC analysis. The area under the curve (AUC) at 1-year, 2-year, and 3-year intervals was 0.708 (95% CI: 62.5%–78.98%), 0.717 (95% CI: 65.56%–77.88%), and 0.739 (95% CI: 68.53%–79.19%), respectively (Figure 2I).

These findings were further validated in three independent validation cohorts (GSE2658, GSE136337, and MMRF-CoMMpass). Consistent with the training set, high-risk patients in these datasets had significantly worse survival outcomes, as shown by KM curves (Figure 3A), and exhibited similar patterns in risk distribution (Figure 3B), gene expression (Figure 3C), and time-dependent ROC performance (Figure 3D). Collectively, these results support the robustness of the HMR prognostic model across multiple datasets.

Figure 3

Figure 3. Validation of HRM prognostic risk model. (A) Kaplan–Meier survival curve of high-risk and low-risk in the GSE2658, GSE136337 and MMRF validation cohort. (B) The distribution of risk scores and survival status of MM in the GSE2658, GSE136337 and MMRF validation cohort. (C) Heatmap of the HRM related genes in the GSE2658, GSE136337 and MMRF validation cohort. (D) Time-dependent ROC analysis in the GSE2658, GSE136337 and MMRF cohort.

HMR risk score is an independent prognostic factor in MM

To explore the relationship between the HMR risk score and clinical characteristics, we divided the patients in the GSE24080 training cohort into different subgroups based on clinical variables. The results of the Wilcoxon rank sum test showed that the HMR risk score was significantly elevated in female patients (P < 0.001), in those with higher levels of lactate dehydrogenase (LDH) (P < 0.001), beta-2-microglobulin (B2M) (P < 0.001), ISS stage (P < 0.001), those with cytogenetic abnormalities (Cyto Abn) (P < 0.001), and those with lower levels of albumin (ALB) (P < 0.05) (Supplementary Figure S1H–N). These results indicated that the HMR risk score correlated with adverse clinical parameters in MM.

Consistently, when comparing clinical characteristics between the high-risk and low-risk groups stratified by the HMR score, the high-risk group showed a significantly higher frequency of these unfavorable features (Supplementary Table S2), further confirming the close association between the HMR risk score and aggressive disease phenotypes in MM.

Next, Univariate and multivariate Cox regression analyses were performed to identify survival associated variables, including clinical parameters such as age, gender, ISS stage, LDH, Cyto Abn and HMR risk score in both training and validation cohort (Supplementary Figure S3). Notably, the HMR score demonstrated a p-value of <0.001 in both univariate and multivariate cox analyses across all cohorts, highlighting its potential as a robust and independent prognostic factor for MM patients.

Nomogram construction for the clinical application of the HMR risk model

We further constructed a nomogram model that integrates clinical features and the HMR risk score to provide a quantitative tool for survival prediction (Figure 4A). Notably, the HMR risk score played a major role in the nomogram model, reflecting its strong prognostic value.

Figure 4

Figure 4. Establish and validation of predictive nomogram to evaluate clinical application. (A) Combined nomogram based on HRM risk score and clinical features in training (GSE24080) and validation (GSE136337 and MMRF) cohorts. (B) Calibration plot of nomogram plot in training (GSE24080) and validation (GSE136337 and MMRF) cohorts. (C) Time-dependent ROC curve in training (GSE24080) and validation (GSE136337 and MMRF) cohorts. (D) DCA curve in training (GSE24080) and validation (GSE136337 and MMRF) cohorts.

The nomogram demonstrated strong predictive performance in both the training and validation cohorts, with C-index value of 0.728 (95% CI: 0.708-0.748) in the training cohort, 0.674 (95% CI: 0.653-0.695) in the GSE136337 validation cohort and 0.756 (95% CI: 0.740-0.772) in MMRF-CoMMpass validation cohort, indicating good discriminative ability (Figure 4B). Calibration curves showed that the predicted survival probabilities at 1-year, 2-year, and 3-year survival were highly consistent with the actual survival rates in both cohorts, further confirming the reliability of the model (Figure 4B).

In the 3-year AUC analysis, the nomogram achieved an AUC of 0.772 in the trainig cohort. The HMR risk score achieved an AUC of 0.740, outperforming LDH (0.669) and ISS stage (0.621), emphasizing its superior prognostic relevance. Similarly, the nomogram achieved an AUC of 0.724 and 0.759 in GSE136337 and MMRF-CoMMpas validation cohort respectively (Figure 4C). Furthermore, Decision curve analysis (DCA) demonstrated that the nomogram provided the highest net benefit across different decision thresholds, with the HMR risk score ranking second and consistently outperforming other clinical variables in both training and validation datasets (Figure 4D).

Together, these findings indicate that the nomogram is a robust and reliable tool for individualized survival prediction in MM patients across both the training and validation cohorts. Notably, the HMR risk score showed superior prognostic performance compared to conventional clinical variables and contributed the most to the nomogram model, underscoring its pivotal role in risk stratification.

GSEA enrichment analysis based on HMR risk score

Differential expression analysis identified 269 DEGs (|logFC| > 0.58, adjusted p-value <0.05), including 216 upregulated and 53 downregulated genes in GSE24080 dataset, and 663 DEGs (|log2FoldChange| > 1, adjusted p-value <0.05), including 359 upregulated and 304 downregulated genes in MMRF-CoMMpass dataset. Volcano plot visualization showed a distinct separation between upregulated and downregulated genes (Figures 5A,B). Gene Ontology (GO) analysis revealed that the significantly enriched biological processes (BP) were mainly involved in chromosome segregation and nuclear division. In the cellular component (CC) category, spindle and chromosomal region were predominantly enriched. The molecular function (MF) category was mainly associated with tubulin binding and microtubule binding. KEGG pathway analysis identified significant enrichment in cell cycle (Figure 5C,D). Consistently, GSEA demonstrated that gene sets related to cell cycle and proliferation were significantly enriched in the high-risk group (Figures 5E–H).

Figure 5

Figure 5. Differential expression and functional enrichment analysis based on HMR risk groups. (A,B) Volcano plot analysis of DEGs associated with HMR risk score in GSE24080 and MMRF datasets. (C,D) Gene Ontology (GO) and KEGG pathway enrichment analysis of DEGs. GSEA enrichment of DEGs associated with HMR risk score revealed a significant enrichment of gene set related to cell cycle (E,F) and proliferation (G,H) in the high-risk group.

Association between HMR risk score and tumor mutation burden (TMB)

To investigate the relationship between TMB and prognosis in MM, we first performed KM survival analysis. The results showed that patients with high TMB had a worse overall survival compared to those with low TMB (P = 0.032) (Figure 6A), suggesting that TMB is an adverse prognostic factor.

Figure 6

Figure 6. Characteristics of somatic mutation between HMR high-risk and low-risk groups. (A) Kaplan-Meier survival curves comparing overall survival between high and low TMB groups. (B) Boxplot comparing TMB levels between HMR high-risk and low-risk groups (left). Scatter plot and correlation analysis demonstrating a positive correlation between TMB and risk score (right). (C,D) Oncoplot visualizing mutation landscape of prognostic genes (KRAS, NRAS, TP53, and MYC) across risk groups. (E–H) Lollipop plots illustrating the distribution and frequency of mutations in KRAS, NRAS, TP53, and MYC in HMR high-risk and low-risk groups.

We then assessed TMB distribution between HMR risk groups. The high-risk group exhibited a significantly higher TMB than the low-risk group (P = 7.6e-05), and a positive correlation was observed between TMB and HMR risk score (R = 0.13, P = 0.00021) (Figure 6B), indicating that genomic instability increases with higher HMR risk scores.

To further characterize genetic alterations in high-risk patients, we analyzed somatic mutations in key MM prognostic genes KRAS, NRAS, TP53, and MYC. The oncoplot and lollipop plot showed higher mutation frequencies of KRAS, NRAS, and TP53 in the high-risk group, particularly for TP53, indicating its key role in MM progression (Figures 6C–H).

Correlation between HMR risk score and immune infiltration in bone marrow

The immune infiltration analysis revealed that the abundance of neutrophils was significantly decreased in the high-risk group than in low-risk group (Figure 7A). Consistently, gene expression analysis showed that neutrophil-associated genes, including LTF, LCN2, CTSG, ITGAM, CRISP3, CA2 and PADI2 were upregulated in the low-risk group, while immune response and inflammation-related genes such as CXCL9, CXCL10, LBP, MMP1, SOCS1, TNFAIP1, and IL10RB were upregulated in high-risk group (Figure 7B).

Figure 7

Figure 7. Immune infiltration and immune-related gene expression analysis between high-risk and low-risk groups. (A) Box plot showing the relative abundance of immune cell types in HMR high-risk and low-risk groups. (B) Box plots show immune-related gene expression levels in HMR high-risk and low-risk groups.

Prediction of chemotherapy sensitivity based on the HMR risk model

To evaluate the sensitivity of high-risk and low-risk patients to chemotherapeutic agents, we predicted the IC50 values of various chemotherapy drugs and pathway inhibitors using the pRRophetic algorithm. The analysis revealed that high-risk patients were more sensitive to agents like vorinostat (a HDAC inhibitor), cytarabine (a chemotherapy drug commonly used in the treatment of various types of leukemia), RO-3306 (a small molecule CDK1 (Cyclin-Dependent Kinase 1) inhibitor), bortezomib and lenalidomide (targeted therapies used in the treatment MM). Conversely, bleomycin (a chemotherapy drug primarily used in the treatment of certain types of cancer), OSI-930, and sorafenib (VEGFR2 inhibitors) showed higher IC50 values in high-risk patients, suggesting that these patients may resistant to these drugs (Figures 8A–F). These findings suggest that the HMR may could serve as a valuable tool in guiding personalized treatment decisions, helping to optimize chemotherapy drug selection for patients.

Figure 8

Figure 8. Prediction of drug sensitivity in HMR high-risk and low-risk groups. Boxplot comparing the drug IC50 levels between HMR high-risk and low-risk groups (upper). Scatter plot and correlation analysis demonstrating a correlation between drug sensitivity and risk score (lower). (A–H) Vorinostat, Bleomycin, Cytarabine, RO-3306, OSI-930, Sorafenib, Bortezomib and Lenalidomide.

Discussion

In this study, we developed and validated an HMR risk prognostic model based on seven epigenetic regulators using five independent datasets including 2790 MM patients. The HMR risk score was significantly correlated with ISS stage, LDH levels, B2M, ALB, and cytogenetic abnormalities. Multivariate Cox regression analysis confirmed it as an independent prognostic factor for MM. The model effectively stratified patients into high-risk and low-risk groups, with the high-risk group exhibiting significantly worse overall survival. Integration HMR score with clinical parameters improved predictive accuracy and clinical decision benefit. These findings underscore the critical role of epigenetic regulation in MM progression and patient survival.

The seven genes comprising the HMR model (SUZ12, KAT2A, AURKA, BUB1, SUV39H2, UTY, and PCGF5) are critically involved in epigenetic modification, tumor progression, supporting their prognostic value in MM. SUZ12, a core component of the Polycomb repressive complex 2 (PRC2 complex), catalyzes H3K27 trimethylation (H3K27me3) to silence gene expression. SUZ12-mediated H3K27me3 modification inhibits HDAC1 expression, and modulate docetaxel resistance in lung adenocarcinoma to (Jiang et al., 2022a). In T-cell acute lymphoblastic leukemia (T-ALL), SUZ12 mutations activate oncogenic signaling pathways, increasing sensitivity to HDAC6 inhibitors (Broux et al., 2019). High SUZ12 and H3K27me3 expression correlates with poor survival in soft tissue sarcomas (Cho et al., 2018). KAT2A, a histone acetyltransferase, promotes cell proliferation and invasion by acetylating H3K9 and activating E2F1 in colorectal cancer (Han and Chen, 2022). It also enhances resistance to hormone therapy in prostate cancer, further supporting its role in oncogenesis (Lu et al., 2021). Aurora kinase A (AURKA) phosphorylates NSD2 at S56, enhancing its methyltransferase activity and promoting chemotherapy resistance in t(4;14) MM patients (Jiang et al., 2022b). It interacts with BRD4, and targeting both BRD4 and AURKA reduces tumor growth (Xu et al., 2020). AURKA promotes the degradation of tumor suppressor p53, further boosting cell proliferation and oncogenic activity (Park et al., 2018). Furthermore, the AURKA inhibitor AT9283 induced apoptosis and inhibited growth in MM (Santo et al., 2011). BUB1 plays a key role in mitosis by enhancing its kinase activity through TIP60-mediated acetylation, leading to H2A T120 phosphorylation, ensuring accurate chromosome segregation. Non-acetylated BUB1 mutations impair this process, leading to genomic instability (Sun et al., 2024). Additionally, BUB1 promotes GC cell proliferation and metastasis through METTL3-mediated m6A methylation, highlighting its critical role in both mitosis and cancer progression (Wang et al., 2024). SUV39H2 mediates H3K9me3 modification, promoting tumor progression (You et al., 2024), poor survival outcomes (Zheng et al., 2018), and chemoresistance (Vougiouklakis et al., 2018; Wang et al., 2019). UTY (Ubiquitous Transcribed Y), the Y-chromosome paralog of UTX (KDM6A) encode lysine demethylases (KDMs), that catalyzes H3K27me3 demethylation in histone H3 (Black et al., 2012). KDM6A and UTY have significant tumor suppressor activity in several types of cancer (Gozdecka et al., 2018; Shi et al., 2021). PCGF5, a PRC1 component, facilitates Xist-mediated gene silencing through H2A ubiquitination and H3K27 methylation (Almeida et al., 2017). PCGF5 is crucial for neural differentiation in mouse ESCs. Its absence activates the SMAD2/TGF-β pathway, impairing differentiation and weakening H2AK119ub1 and H3K27me3, maintaining gene repression (Yao et al., 2018). Collectively, these genes modulate key histone modifications such as H3K27me3, H3K9me3, H2AK119ub1, and H3K9ac, which play important role in transcriptional regulation, DNA damage response, and chromatin accessibility (Alzrigat et al., 2018). Their dysregulation may disrupt normal gene expression, contributing to increased cell proliferation, impaired apoptosis, and therapeutic resistance, which are features associated with poor prognosis in MM.

Functional enrichment analysis revealed significant upregulation of chromosomal segregation, cell cycle, and proliferation pathways in the high-risk group. These processes promote chromosomal instability (CIN), leading to intratumoral heterogeneity and MM progression (Chen et al., 2025). Furthermore, the high-risk group exhibited a higher TMB, particularly with increased mutations in key oncogenic genes such as KRAS, NRAS, TP53, and MYC. Mutations of the RAS/mitogen-activated protein kinase (MAPK) pathway, including NRAS, KRAS, or BRAF, are found in up to 50% of newly diagnosed MM patients (Walker et al., 2015), and contribute to proteasome inhibitor (PI) resistance (Shirazi et al., 2020). TP53 mutation lead to decreased expression of p53 and associated with poor outcomes and drug resistance in MM (Jovanović et al., 2018). These findings suggest that epigenetic dysregulation may promote genomic instability and disease progression in MM.

Immune infiltration analysis showed reduced neutrophils in the high-risk group, potentially impairing anti-tumor immunity (Chan et al., 2023; Zhu et al., 2024). In parallel, gene expression profiling showed that neutrophil-related genes, including LTF, LCN2, CTSG, ITGAM, and PADI2, were consistently upregulated in the low-risk group, suggesting an enrichment of neutrophil-mediated innate immune activity in these patients. In contrast, the high-risk group demonstrated elevated expression of inflammatory and immunomodulatory genes, such as CXCL9, CXCL10, SOCS1, and IL10RB. These immune signatures may help explain the differential prognosis between groups and may guide personalized immunotherapy approches based on HMR score.

Furthermore, drug sensitivity analysis highlighted the potential of the HMR risk score to guide personalized treatment strategies. High-risk patients exhibited increased sensitivity to HDAC inhibitors (e.g., vorinostat), CDK1 inhibitors, proteasome inhibitors (e.g., bortezomib), and immunomodulatory agents (e.g., lenalidomide), while showing resistance to certain chemotherapy drugs and kinase inhibitors. This highlights the potential for personalized treatment strategies based on HMR risk scores. A recent study showed that AZD1775, a Wee1 inhibitor targeting the cell cycle checkpoint, synergizes with the HDAC inhibitor Vorinostat to induce DNA damage and apoptosis in leukemia cells, including p53-wild type and deficient AML (Zhou et al., 2015). Targeting both epigenetic modifiers (e.g., HDAC inhibitors) and cell cycle regulators (e.g., CDK1 inhibitors) may have synergistic effects in high-risk MM patients.

However, several limitations in this study should be acknowledged. First, although the model was validated by three independent datasets, further validation in prospective clinical cohorts is needed. Second, while this study relied on transcriptomic data, the lack of matched epigenomic datasets limited our ability to link gene expression with chromatin regulatory mechanisms. Integration of multi-omics data in future studies will provide a more comprehensive understanding of the model’s biological basis. Third, functional characterization of the seven genes remains to be fully elucidated. Future studies employing in vitro and in vivo experiments are needed to elucidate the mechanistic role of these genes in MM pathogenesis and disease progression.

Conclusion

In conclusion, our study developed and validated a novel HMR risk signature for MM, which exhibited independent prognostic value across multiple cohorts. Integration the HMR score with clinical parameters improved predictive accuracy and clinical decision benefit in both training and validation cohorts. Additionally, the HMR score showed associations with genomic instability, the tumor immune microenvironment, and drug sensitivity, suggesting its potential relevance to MM biology. These findings may contribute to a better understanding of MM progression and offer insights that could inform personalized therapeutic strategies.

Data availability statement

The datasets used and analysed during the current study were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) and the MMRF-CoMMpass project on the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/).

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

JL: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Visualization, Writing – original draft, Software. SL: Data curation, Validation, Writing – review and editing. YQ: Funding acquisition, Writing – review and editing, Methodology. YF: Writing – review and editing, Resources, Validation. ZZ: Resources, Writing – review and editing, Software, Supervision. LZ: Resources, Supervision, Writing – review and editing, Conceptualization, Data curation.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. The study was supported by the Zhejiang Provincial Natural Science Foundation of China under Grand No. LQ24H200001, Medical Science and Technology Project of Zhejiang Provincial Health Commission Grand 2023KY353, 2024KY462, and by Shaoxing City Health Science and Technology Plan Grand No. 2022SY020.

Acknowledgments

All authors would like to appreciate the TCGA database and the GEO database for sharing data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2025.1613631/full#supplementary-material

SUPPLEMENTARY FIGURE S1 | (A–G) Kaplan-Meier (KM) curves compare survival rates between the high-expression and low-expression groups, which were defined based on the median expression level of AURKA, BUB1, KAT2A, PCGF5, SUV39H2, SUZ12 and UTY. (H-N) Violin plots illustrating the relationship between HMR risk score and age, gender, lactate dehydrogenase (LDH), β2-microglobulin (B2M), albumin (ALB), cytogenetic abnormalities (Cyto Abn), and ISS stage.

SUPPLEMENTARY FIGURE S2 | (A) Pairwise correlation analysis among the seven genes. (B) Bar plot showing the variance inflation factor (VIF) values for each gene. (C) Bar plot showing the coefficients for each gene obtained from multivariate Cox regression.

SUPPLEMENTARY FIGURE S3 | (A) Univariate and (B) multivariate Cox regression analyses identifying prognostic variables in both the training (GSE24080) and validation cohorts (GSE136337 and MMRF-CoMMpass).

References

Almeida, M., Pintacuda, G., Masui, O., Koseki, Y., Gdula, M., Cerase, A., et al. (2017). PCGF3/5-PRC1 initiates polycomb recruitment in X chromosome inactivation. Science. 356 (6342), 1081–1084. doi:10.1126/science.aal2512

PubMed Abstract | CrossRef Full Text | Google Scholar

Alzrigat, M., Párraga, A. A., and Jernberg-Wiklund, H. (2018). Epigenetics in multiple myeloma: from mechanisms to therapy. Seminars Cancer Biol. 51, 101–115. doi:10.1016/j.semcancer.2017.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Attal, M., Lauwers-Cances, V., Hulin, C., Leleu, X., Caillot, D., Escoffre, M., et al. (2017). Lenalidomide, bortezomib, and dexamethasone with transplantation for myeloma. N. Engl. J. Med. 376 (14), 1311–1320. doi:10.1056/NEJMoa1611750

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, L., Wang, X., Li, J., Li, W., and Wang, Z. (2020)2025). Epigenetic modifications in early stage lung cancer: pathogenesis, biomarkers, and early diagnosis. MedComm 6 (3), e70080. doi:10.1002/mco2.70080

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, J. C., Van Rechem, C., and Whetstine, J. R. (2012). Histone lysine methylation dynamics: establishment, regulation, and biological impact. Mol. Cell. 48 (4), 491–507. doi:10.1016/j.molcel.2012.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Broux, M., Prieto, C., Demeyer, S., Vanden Bempt, M., Alberti-Servera, L., Lodewijckx, I., et al. (2019). Suz12 inactivation cooperates with JAK3 mutant signaling in the development of T-cell acute lymphoblastic leukemia. Blood 134 (16), 1323–1336. doi:10.1182/blood.2019000015

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, Y. T., yue, T. H., Lu, Y., Zhang, C., Cheng, C. S., Wu, J., et al. (2023). Pancreatic melatonin enhances anti-tumor immunity in pancreatic adenocarcinoma through regulating tumor-associated neutrophils infiltration and NETosis. Acta Pharm. Sin. B 13 (4), 1554–1567. doi:10.1016/j.apsb.2023.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Chase, A., and Cross, N. C. P. (2011). Aberrations of EZH2 in cancer. Clin. Cancer Res. 17 (9), 2613–2618. doi:10.1158/1078-0432.CCR-10-2156

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Agustinus, A. S., Li, J., DiBona, M., and Bakhoum, S. F. (2025). Chromosomal instability as a driver of cancer progression. Nat. Rev. Genet. 26 (1), 31–46. doi:10.1038/s41576-024-00761-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Chervona, Y., and Costa, M. (2012). Histone modifications and cancer: biomarkers of prognosis? Am. J. Cancer Res. 2 (5), 589–597.

PubMed Abstract | Google Scholar

Cho, Y. J., Kim, S. H., Kim, E. K., Han, J. W., Shin, K. H., Hu, H., et al. (2018). Prognostic implications of polycomb proteins ezH2, suz12, and eed1 and histone modification by H3K27me3 in sarcoma. BMC Cancer 18 (1), 158. doi:10.1186/s12885-018-4066-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Danziger, S. A., McConnell, M., Gockley, J., Young, M. H., Rosenthal, A., Schmitz, F., et al. (2020). Bone marrow microenvironments that contribute to patient outcomes in newly diagnosed multiple myeloma: a cohort study of patients in the total therapy clinical trials. PLoS Med. 17 (11), e1003323. doi:10.1371/journal.pmed.1003323

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N. J., and Huang, R. S. (2014a). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. 15 (3), R47. doi:10.1186/gb-2014-15-3-r47

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014b). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Girija Sivasankaran, R., Edachalil Veeraraghavan, M., and Sidharthan, N. (2025). Documentation of international staging system (ISS) and revised international staging system (R-ISS) for multiple myeloma at a tertiary care center: a clinical audit. Ind. J. Clin. Biochem. 40 (1), 151–154. doi:10.1007/s12291-023-01152-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gozdecka, M., Meduri, E., Mazan, M., Tzelepis, K., Dudek, M., Knights, A. J., et al. (2018). UTX-Mediated enhancer and chromatin remodeling suppresses myeloid leukemogenesis through noncatalytic inverse regulation of ETS and GATA programs. Nat. Genet. 50 (6), 883–894. doi:10.1038/s41588-018-0114-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagen, P., Zhang, J., and Barton, K. (2022). High-risk disease in newly diagnosed multiple myeloma: beyond the R-ISS and IMWG definitions. Blood Cancer J. 12 (5), 83. doi:10.1038/s41408-022-00679-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, X., and Chen, J. (2022). KAT2A affects tumor metabolic reprogramming in Colon cancer progression through epigenetic activation of E2F1. Hum. Cell. 35 (4), 1140–1158. doi:10.1007/s13577-022-00707-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernando, H., Gelato, K. A., Lesche, R., Beckmann, G., Koehr, S., Otto, S., et al. (2016). EZH2 inhibition blocks multiple myeloma cell growth through upregulation of epithelial tumor suppressor genes. Mol. Cancer Ther. 15 (2), 287–298. doi:10.1158/1535-7163.MCT-15-0486

PubMed Abstract | CrossRef Full Text | Google Scholar

Ismail, N. H., Mussa, A., Zakaria, N. A., Al-Khreisat, M. J., Zahidin, M. A., Ramli, N. N., et al. (2022). The role of epigenetics in the development and progression of multiple myeloma. Biomedicines 10 (11), 2767. doi:10.3390/biomedicines10112767

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, Y., Xiao, C., Fan, T., Deng, Z., Wang, D., Cai, W., et al. (2025). The epigenetic hallmarks of immune cells in cancer. Mol. Cancer 24 (1), 66. doi:10.1186/s12943-025-02255-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, M., Qi, F., Zhang, K., Zhang, X., Ma, J., Xia, S., et al. (2022a). MARCKSL1-2 reverses docetaxel-resistance of lung adenocarcinoma cells by recruiting SUZ12 to suppress HDAC1 and elevate miR-200b. Mol. Cancer 21 (1), 150. doi:10.1186/s12943-022-01605-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, H., Wang, Y., Wang, J., Wang, Y., Wang, S., He, E., et al. (2022b). Posttranslational modification of Aurora a-NSD2 loop contributes to drug resistance in t(4;14) multiple myeloma. Clin. Transl. Med. 12 (4), e744. doi:10.1002/ctm2.744

PubMed Abstract | CrossRef Full Text | Google Scholar

Jovanović, K. K., Escure, G., Demonchy, J., Willaume, A., Van de Wyngaert, Z., Farhat, M., et al. (2018). Deregulation and targeting of TP53 pathway in multiple myeloma. Front. Oncol. 8, 665. doi:10.3389/fonc.2018.00665

PubMed Abstract | CrossRef Full Text | Google Scholar

Karimi, F., Aghaei, M., and Saki, N. (2025). Impact of genetic polymorphisms on treatment outcomes of proteasome inhibitors and immunomodulatory drugs in multiple myeloma. Curr. Treat. Options Oncol. 26, 197–212. doi:10.1007/s11864-025-01295-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kouzarides, T. (2007). Chromatin modifications and their function. Cell. 128 (4), 693–705. doi:10.1016/j.cell.2007.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kyle, R. A., and Rajkumar, S. V. (2008). Multiple myeloma. Blood 111 (6), 2962–2972. doi:10.1182/blood-2007-10-078022

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. E., and Kim, M. Y. (2022). Cancer epigenetics: past, present and future. Seminars Cancer Biol. 83, 4–14. doi:10.1016/j.semcancer.2021.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., Carey, M., and Workman, J. L. (2007). The role of chromatin during transcription. Cell. 128 (4), 707–719. doi:10.1016/j.cell.2007.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. R., Wang, C., and Cheng, C. (2024). Identifying high-risk multiple myeloma patients: a novel approach using a clonal gene signature. Int. J. Cancer 155 (9), 1684–1695. doi:10.1002/ijc.35057

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, D., Song, Y., Yu, Y., Wang, D., Liu, B., Chen, L., et al. (2021). KAT2A-mediated AR translocation into nucleus promotes abiraterone-resistance in castration-resistant prostate cancer. Cell. Death Dis. 12 (8), 787. doi:10.1038/s41419-021-04077-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Maiso, P., Carvajal-Vergara, X., Ocio, E. M., López-Pérez, R., Mateo, G., Gutiérrez, N., et al. (2006). The histone deacetylase inhibitor LBH589 is a potent antimyeloma agent that overcomes drug resistance. Cancer Res. 66 (11), 5781–5789. doi:10.1158/0008-5472.CAN-05-4186

PubMed Abstract | CrossRef Full Text | Google Scholar

Margueron, R., and Reinberg, D. (2011). The polycomb complex PRC2 and its mark in life. Nature 469 (7330), 343–349. doi:10.1038/nature09784

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Garcia, E., Popovic, R., Min, D. J., Sweet, S. M. M., Thomas, P. M., Zamdborg, L., et al. (2011). The MMSET histone methyl transferase switches global histone methylation and alters gene expression in t(4;14) multiple myeloma cells. Blood 117 (1), 211–220. doi:10.1182/blood-2010-07-298349

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagano, T., and Fraser, P. (2009). Emerging similarities in epigenetic gene silencing by long noncoding RNAs. Mamm. Genome 20 (9), 557–562. doi:10.1007/s00335-009-9218-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Nylund, P., Atienza Párraga, A., Haglöf, J., De Bruyne, E., Menu, E., Garrido-Zabala, B., et al. (2021). A distinct metabolic response characterizes sensitivity to EZH2 inhibition in multiple myeloma. Cell. Death Dis. 12 (2), 167. doi:10.1038/s41419-021-03447-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, J. W., Chae, Y. C., Kim, J. Y., Oh, H., and Seo, S. B. (2018). Methylation of aurora kinase a by MMSET reduces p53 stability and regulates cell proliferation and apoptosis. Oncogene 37 (48), 6212–6224. doi:10.1038/s41388-018-0393-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pawlyn, C., Bright, M. D., Buros, A. F., Stein, C. K., Walters, Z., Aronson, L. I., et al. (2017). Overexpression of EZH2 in multiple myeloma is associated with poor prognosis and dysregulation of cell cycle control. Blood Cancer J. 7 (3), e549. doi:10.1038/bcj.2017.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Pu, J., Liu, T., Wang, X., Sharma, A., Schmidt-Wolf, I. G. H., Jiang, L., et al. (2024). Exploring the role of histone deacetylase and histone deacetylase inhibitors in the context of multiple myeloma: mechanisms, therapeutic implications, and future perspectives. Exp. Hematol. Oncol. 13 (1), 45. doi:10.1186/s40164-024-00507-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Santo, L., Hideshima, T., Cirstea, D., Bandi, M., Nelson, E. A., Gorgun, G., et al. (2011). Antimyeloma activity of a multitargeted kinase inhibitor, AT9283, via potent aurora kinase and STAT3 inhibition either alone or in combination with lenalidomide. Clin. Cancer Res. 17 (10), 3259–3271. doi:10.1158/1078-0432.CCR-10-3012

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, S., Kelly, T. K., and Jones, P. A. (2010). Epigenetics in cancer. Carcinogenesis 31 (1), 27–36. doi:10.1093/carcin/bgp220

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, B., Li, W., Song, Y., Wang, Z., Ju, R., Ulman, A., et al. (2021). UTX condensation underlies its tumour-suppressive activity. Nat. 2021 Sept. 597 (7878), 726–731. doi:10.1038/s41586-021-03903-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirazi, F., Jones, R. J., Singh, R. K., Zou, J., Kuiatse, I., Berkova, Z., et al. (2020). Activating KRAS, NRAS, and BRAF mutants enhance proteasome capacity and reduce endoplasmic reticulum stress in multiple myeloma. Proc. Natl. Acad. Sci. U. S. A. 117 (33), 20004–20014. doi:10.1073/pnas.2005052117

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2022). Cancer statistics, 2022. CA Cancer J. Clin. 72 (1), 7–33. doi:10.3322/caac.21708

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, M., Yang, B., Xin, G., Wang, Y., Luo, J., Jiang, Q., et al. (2024). TIP60 acetylation of Bub1 regulates centromeric H2AT120 phosphorylation for faithful chromosome segregation. Sci. China Life Sci. 67 (9), 1957–1969. doi:10.1007/s11427-023-2604-8

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Donk, NWCJ, Pawlyn, C., and Yong, K. L. (2021). Multiple myeloma. Lancet 397 (10272), 410–427. doi:10.1016/S0140-6736(21)00135-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Vougiouklakis, T., Saloura, V., Park, J. H., Takamatsu, N., Miyamoto, T., Nakamura, Y., et al. (2018). Development of novel SUV39H2 inhibitors that exhibit growth suppressive effects in mouse xenograft models and regulate the phosphorylation of H2AX. Oncotarget 9 (61), 31820–31831. doi:10.18632/oncotarget.25806

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, B. A., Boyle, E. M., Wardell, C. P., Murison, A., Begum, D. B., Dahir, N. M., et al. (2015). Mutational spectrum, copy number changes, and outcome: results of a sequencing study of patients with newly diagnosed myeloma. JCO 33 (33), 3911–3920. doi:10.1200/JCO.2014.59.1503

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Cheng, L., Yang, X., Chen, X., Miao, Y., Qiu, Y., et al. (2019). Histone methyltransferase SUV39H2 regulates cell growth and chemosensitivity in glioma via regulation of hedgehog signaling. Cancer Cell. Int. 19, 269. doi:10.1186/s12935-019-0982-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Shen, K., Wang, J., Yang, K., Zhu, J., Chen, Y., et al. (2024). BUB1 potentiates gastric cancer proliferation and metastasis by activating TRAF6/NF-κB/FGF18 through m6A modification. Life Sci. 353, 122916. doi:10.1016/j.lfs.2024.122916

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Z., Bi, C., Chooi, J. Y., Chan, Z. L., Mustafa, N., and Chng, W. J. (2015). MMSET regulates expression of IRF4 in t(4;14) myeloma and its silencing potentiates the effect of bortezomib. Leukemia 29 (12), 2347–2354. doi:10.1038/leu.2015.169

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J. L., Yuan, Y. J., Lv, J., Qi, D., Wu, M. D., Lan, J., et al. (2020). Inhibition of BRD4 triggers cellular senescence through suppressing aurora kinases in oesophageal cancer cells. J. Cell. Mol. Med. 24 (22), 13036–13045. doi:10.1111/jcmm.15901

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, M., Zhou, X., Zhou, J., Gong, S., Hu, G., Li, J., et al. (2018). PCGF5 is required for neural differentiation of embryonic stem cells. Nat. Commun. 9 (1), 1463. doi:10.1038/s41467-018-03781-0

PubMed Abstract | CrossRef Full Text | Google Scholar

You, J., Xue, H., Chao, C., Zhang, Z., Tan, X., Wang, X., et al. (2024). Histone methyltransferase SUV39H2 supports nasopharyngeal carcinoma cell metastasis by regulation of SIRT1. Environ. Toxicol. 39 (11), 4974–4983. doi:10.1002/tox.24370

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Y., Li, B., Wang, J., Xiong, Y., Wang, K., Qi, Y., et al. (2018). Identification of SUV39H2 as a potential oncogene in lung adenocarcinoma. Clin. Epigenetics 10 (1), 129. doi:10.1186/s13148-018-0562-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Zhang, Y., Chen, S., Kmieciak, M., Leng, Y., Lin, H., et al. (2015). A regimen combining the Wee1 inhibitor AZD1775 with HDAC inhibitors targets human acute myeloid leukemia cells harboring various genetic mutations. Leukemia 29 (4), 807–818. doi:10.1038/leu.2014.296

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, M., Wang, S., Qu, K., Lu, F., Kou, M., Yao, Y., et al. (2024). The trogocytosis of neutrophils on initial transplanted tumor in mice. iScience 27 (5), 109661. doi:10.1016/j.isci.2024.109661

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuge, L., Lin, X., Fan, Z., Jia, M., Lin, C., Zhu, M., et al. (2025). Global, regional and national epidemiological trends of multiple myeloma from 1990 to 2021: a systematic analysis of the global burden of disease study 2021. Front. Public Health 13, 1527198. doi:10.3389/fpubh.2025.1527198

PubMed Abstract | CrossRef Full Text | Google Scholar

Glossary

MM Multiple myeloma

HMR histone modification-related

GEO Gene Expression Omnibus

TCGA The Cancer Genome Atlas

LASSO least absolute shrinkage and selection operator

ISS International Staging System

HDACs histone deacetylases

HMTs histone methyltransferases

FGFR3 fibroblast growth factor receptor 3

MMSET multiple myeloma SET domain

NSD Nuclear receptor-binding SET domain

PRC2 Polycomb repressive complex 2

WBM whole bone marrow

KM Kaplan–Meier

ROC receiver operating characteristic

DCA decision curve analysis

H3K27me3 H3K27 trimethylation

T-ALL T-cell acute lymphoblastic leukemia

DEGs Differential expression genes

GO Gene Ontology

BP biological processes

CC cellular component

MF molecular function

TMB tumor mutation burden

KEGG Kyoto Encyclopedia of Genes and Genomes

GSEA Gene Set Enrichment Analysis

B2M β2-microglobulin

LDH lactate dehydrogenase

ALB albumin

HGB Hemoglobin

Cyto Abn cytogenetic abnormalities

CDK1 Cyclin-Dependent Kinase 1

AURKA Aurora kinase A

γ-H2AX phosphorylated H2AX

UTY Ubiquitous Transcribed Y

CIN chromosomal instability