The gene expression matrix and clinical information of all samples were downloaded from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). The raw count data was preprocessed with normalization and log2 transformation. The overall design and workflow of this study is given in Supplementary Figure S1 . Briefly, The GSE47552 was applied for the identification of differentially expressed m6A-related genes. The GSE47552, GSE6477, and GSE13591 were used for the evaluation of risk score for MM diagnosis. Three independent series GSE24080, GSE9782, and GSE57317 were identified using the following selection criteria: (1) >40 subjects, (2) available survival data, and (3) confirmed MM patients. Three hundred thirteen MM patients from the GSE24080 was used as the training cohort for model construction, and the GSE9782 and GSE57319 were used as validation cohorts. Detailed information of all datasets used in this study are shown in Supplementary Table S1 .

We systematically analyzed the expression of 22 m6A-related genes according to previously published articles (27), including 7 writers (METTL3, METTL14, WTAP, KIAA1429, RBM15, RBM15B, and ZC3H13), 3 erasers (FTO, ALKBH5, and ALKBH3), and 12 readers (IGF2BP1/2/3, YTHDF1/2/3, YTHDC1/2, HNRNPA2B1, HNRNPC, and RBMX) in GSE47552. The heatmap was generated via the PHEATMAP package (version 1.0.12; https://cran.rstudio.com/web/packages/pheatmap/index.html).

The PPI network among the differentially expressed m6A regulatory genes was constructed using the STRING online database (www.string-db.org). Cytoscape software (version 3.8.2) was used for visualization (28). Spearman correlation analysis was performed to demonstrate the association among different m6A-related genes based on the Spearman’s correlation coefficient (r) value between −1 and +1. The CORRPLOT package (version 0.84; https://github.com/taiyun/corrplot) was used to visualize the correlation matrix.

First, all MM patients in the training cohort were divided into two groups according to the mean expression level of each gene. Kaplan–Meier survival analysis with log-rank test was used to evaluate the survival rate by the “SURVIVAL” package (version 3.2-7; https://github.com/therneau/survival). We selected four genes (ZC3H13, HNRNPA2B1, HNRNPC, and YTHDF2) showing significant correlation with the overall survival (OS) as prognostic-related genes. Second, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis to select the optimal weighting coefficient and build the risk model with the “glmnet” package (version 4.1-1; https://glmnet.stanford.edu/). Optimal values of penalty parameter lambda were determined by 1,000-fold cross-validation via the minimum criteria (29). The risk score formula was as follows: risk score = (−0.4820058*expression value of ZC3H13) + (0.4641464*expression value of HNRNPA2B1) + (0.3691959*expression value of HNRNPC) + (0.2454545*expression value of YTHDF2). Third, we divided MM patients of the training cohort into high- and low-risk groups upon the optimal cutoff of the risk score with the “Survminer” package (version 0.4.9; https://rpkgs.datanovia.com/survminer/index.html) and analyzed the survival rate by Kaplan–Meier survival analysis. Subsequently, the predictive performance of the m6A risk score was assessed using ROC curves and the value of area under the ROC curve (AUC). The reliability and stability of this prognostic risk model was further ensured in two validation cohorts.

Univariate and multivariate Cox regression analyses were employed to evaluate the prognostic efficiency of other clinical features and to identify whether the risk score was an independent prognostic factor. Briefly, variables with p < 0.10 in univariate analysis were eligible for multivariate Cox regression analysis with forward likelihood ration (LR) method, and p < 0.05 was considered statistically different. Then, the nomogram with independent prognostic factors was constructed using the “rms” package (version 6.2-0; https://hbiostat.org/R/rms/). The predictive performance of the nomogram was evaluated with calibration curves and ROC curves.

Three hundred thirteen MM patients from the GSE24080 were stratified into high- and low-risk groups according to the cutoff value in the survival analysis. Differential expression analysis was applied using the LIMMA package (30) (version 3.46.0; http://bioinf.wehi.edu.au/limma/). An absolute log2 fold change (FC) > 0.6 and adjusted p < 0.05 was used as the cutoff for differentially expressed genes (DEGs).

The online website WebGestalt (http://www.webgestalt.org/) (31) was applied to analyze enriched biological terms of DEGs, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. Gene set enrichment analysis (GSEA) was conducted with GSEA software (version 4.1.0). Official gene sets, including Hallmark (H), KEGG (C2), and Gene ontology (C5), were downloaded from GSEA website (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp) for enrichment. A permutation number of 1,000 was adopted. A p-value cutoff of 0.05 with a false discovery rate (FDR q-value) < 0.05 was considered statistically significant.

The ESTIMATE algorithm (32) was used to calculate the immune score, stromal score, and tumor purity of each sample for preliminary evaluation. Next, we obtained gene sets marking each tumor microenvironment infiltration immune cell type from previous studies (33, 34). We performed single sample GSEA (ssGSEA) to derive the enrichment score of each immune-associated gene set to represent the relative abundance of each infiltrating cell in each sample by using “GSVA” package (version 1.38.2; https://github.com/rcastelo/GSVA).

Statistical analysis was conducted using SPSS 26.0, R Studio (version 1.4.1103; https://rstudio.com/) and GraphPad Prism software (version 8.0.2). Continuous variables were described as the mean ± standard deviation (SD); Mann–Whitney test and Student’s t-test were used to compare the difference in subgroups as appropriate. One-way ANOVA followed by Bonferroni post-hoc comparison were employed to compare the difference of more than two subgroups. p < 0.05 was considered statistically significant if not specified. All statistical tests were two-sided.

The gene expression of 22 m6A-related genes in 5 normal plasm cell (NPC) samples and 41 MM patients from the GSE47552 was used for differentially expressed analysis. The heatmap showed that the expression of nine m6A methylation regulators (FTO, HNRNPC, IGF2BP3, KIAA1429, METTL16, METTL14, METTL3, RBMX, and ZC3H13) were differentially expressed in MM patients compared to normal volunteers ( Figure 1A ). We observed that RBMX and HNRNPC were significantly upregulated, while METTL3, METTL14, METTL16, ZC3H13, KIAA1429, FTO, and IGF2BP3 were significantly downregulated in MM patients ( Figure 1B ). Furthermore, we used the online database cBioPortal (http://www.cbioportal.org/) to determine the mutant frequency of m6A-related genes in MM samples; however, no frequent mutation was found ( Supplementary Figure S2 ) (35, 36).

To further investigate the relationship among nine m6A RNA methylation regulators, PPI network was established ( Figure 1C ), which demonstrated that KIAA1429 and METTL3 were hub genes. In addition, IGF2BP3 had no protein–protein association with others, which was hidden in the figure. Moreover, the correlation of each gene expression level was analyzed by the Spearman correlation analysis ( Figure 1D ). The correlation between METTL3 and METTL14 was most positively significant (r = 0.26), and the relationship between METTL3 and IGF2BP3, FTO, and RBMX were most significantly negative (r = −0.16).

Three hundred thirteen IgG MM patients from the GSE24080 were divided into low- and high-expression groups according to the mean expression level of each m6A-related gene. We found that the overexpression of YTHDF2 (p = 0.00248), HNRNPC (p = 0.02291), and HNRNPA2B1 (p = 0.02719) was associated with the poor OS, and survival outcomes of patients with high ZC3H13 expression were significantly prolonged than that of patients with low expression (p = 0.01062, Figures 2A–D ). However, the expression levels of other genes were not correlated with OS ( Supplementary Figure S3 ). Subsequently, based on the result of LASSO Cox regression analysis using the penalized maximum likelihood estimator and min criteria for the optimal lambda value, four m6A methylation regulators, including ZC3H13, HNRNPA2B1, HNRNPC, and YTHDF2, were identified to construct the risk model ( Figures 2E, F ). Then, we calculated the risk score for each patient of the training cohort and stratified them into low- and high-risk groups according to the optimal cutoff value (7.060, calculated by R program). Kaplan–Meier survival analysis showed that MM patients in high-risk group had significantly shorter OS than patients in low-risk group (p < 0.0001, Figure 2G ). The ROC curve analysis was used to investigate the predictive specificity and sensitivity of the risk score. The AUC values of the risk score for 1-, 2-, 3-, 4-, and 5-year survival were 0.5529 (95%CI, 0.4202–0.6856, p = 0.4286), 0.6020 (95%CI, 0.5089–0.6952, p = 0.0351), 0.5979 (95%CI, 0.5170–0.6789, p = 0.0185), 0.6390 (95%CI, 0.5646–0.7133, p = 0.0004), and 0.6573 (95%CI, 0.5785–0.7360, p = 0.0002), respectively ( Figure 2H ).The AUC values indicated that the risk score had a great discrimination ability for the prognosis of MM patients.

Next, Kaplan–Meier survival analysis and ROC curve for the risk score were conducted in two validation cohorts. The OS of low-risk patients was significantly longer than that of high-risk patients in the GSE9782 (p < 0.0001, Figure 3A ) and GSE57317 (p < 0.0001, Figure 3C ). The AUC values of the risk score for 1- and 2-year survival were 0.5566 (95%CI, 0.4827–0.6305, p = 0.1329) and 0.6388 (95%CI, 0.5558–0.7219, p = 0.0019) in the GSE9782 ( Figure 3B ). Parallel values in the GSE57137 were 0.8488 (95%CI, 0.7224–0.9753, p = 0.0061) and 0.8333 (95%CI, 0.6841–0.9825, p = 0.0014, Figure 3D ). Our results demonstrated the great efficiency of the prognostic risk score for MM patients.

The development of MM has been widely identified as a multistage process starting with the premalignant condition monoclonal gammopathy of unknown significance (MGUS), followed by the intermediate stage smolder MM (SMM), and some progressed into the refractory and/or relapse MM (RRMM) or plasma cell leukemia (PCL). We observed that the risk score was significantly increased with the development of MM ( Figures 3E–G ). The ROC curves showed that the risk score could discriminate between MM and normal individuals with high sensitivity and specificity; the AUC values were 0.9383 (95%CI, 0.8854–0.9913, p = 0.0009), 0.7005 (95%CI, 0.5703–0.8306, p = 0.0149), and 1.000 (95%CI, 1.000–1.000, p = 0.0003) in GSE13591, GSE6477, and GSE47552, respectively ( Figure 3H ).

The GSE24080 was used for univariate and multivariate Cox regression analyses to identify the independent prognostic factor of the risk score and other relevant clinicopathological factors, including age, gender, race, lactic dehydrogenase (LDH), albumin (ALB), hemoglobin (HGB), ISS stage, bone marrow plasm cell (BMPC), bone lesion, and cytogenetic abnormality. The results showed that LDH [hazard ratio (HR), 2.040; 95%CI, 1.327–3.137; p = 0.001], ISS stage (HR, 1.869; 95%CI, 1.217–2.869; p = 0.004), BMPC (HR, 1.601; 95%CI, 1.027–2.497; p = 0.038), cytogenetic abnormalities (HR, 1.819; 95%CI, 1.202–2.753, p = 0.005), bone lesion (HR, 2.323; 95%CI, 1.411–3.826; p = 0.001), and the risk score (HR, 1.928; 95%CI, 1.278–2.907; p = 0.002) were significantly associated with OS ( Table 1 and Figure 2I ), which indicated that the risk score based on four m6A regulators is an independent prognostic factor for MM patients. To quantify the OS predication, we integrated the risk sore and independent clinical features (LDH, ISS stage, BMPC, bone lesion, and cytogenetic abnormalities) to construct a nomogram ( Figure 4A ). The calibration curves illustrated a great agreement between predicted and actual survival probabilities ( Figures 4B–D ). The nomogram risk score system greatly improved the prognostic risk signature model. The AUC values of the nomogram risk score for 1-, 2-, 3-, 4-, and 5- year survival were 0.8087 (95%CI, 0.7327–0.8847; p < 0.0001), 0.7904 (95%CI, 0.7242–0.8566; p < 0.0001), 0.7627 (95%CI, 0.6991–0.8262; p < 0.0001), 0.7696 (95%CI, 0.7091–0.8301; p < 0.0001), and 0.8072 (95%CI, 0.7457–0.8688; p < 0.0001), respectively ( Figure 4E ).

Subgroup analysis upon age, gender, International Staging System (ISS), cytogenetic abnormality in the training cohort, and upon therapeutic strategies in both the training cohort and the validation cohort GSE9782 was performed to evaluate the accuracy and stability of the risk score for different patients. The risk score presented stable predictive efficiency in all age and gender (p < 0.05, Figures 5A–D ). Furthermore, higher risk score was correlated with poorer OS in patients with cytogenetic abnormalities (p = 7e−04, Figure 5E ), while no significant survival difference was observed in patients without cytogenetic abnormalities (p = 0.05931, Figure 5F ). ISS I- and ISS III-stage MM patients with high risk score had significantly shorted OS than those with low risk score (p = 0.00442 and 0.00652, respectively; Figures 5G, H ), while no significant difference was observed in ISS II stage (p = 0.1365; Figure 5I ). In addition, subgroup analysis for therapeutic regimens verified the robustness of the prognostic risk signature in different treatment choices, including total therapy 2 (TT2, a combination treatment with thalidomide, p = 0.00026, Figure 5J ), TT3 (a combination treatment with the addition of bortezomib, p = 0.00434, Figure 5K ), dexamethasone (DEX, p = 0.02434, Figure 5L ), and PS341 (a kind of proteasome inhibitor, p = 0.00113, Figure 5M ), which demonstrated that the prognostic risk score meets individualized therapeutic needs of MM patients.

We then analyzed the correlation between the risk score and the clinicopathological features of MM patients in the GSE24080 cohort ( Table 2 ). We observed significantly increased risk scores regarding advanced ISS stage (p = 0.013). Moreover, MM patients with cytogenetic abnormalities were observed to have higher risk scores than those cytogenetically normal (p < 0.0001). Subsequently, we analyzed the relationship between the expression of the four prognostic risk signature genes and clinical features ( Supplementary Table S2 ). We observed that high expression of YTHDF2 and HNRNPA2B1 were associated with advanced ISS stage. The higher expression of HNRNPA2B1 was correlated with higher BMPC ratio. Additionally, MM patients with cytogenetic abnormalities were observed to have higher HNRNPC expression and lower ZC3H13 expression.

In the training cohort, we conducted the differential expression analysis between high- and low-risk groups. A total of 156 upregulated and 90 downregulated genes were identified to be significantly associated with the risk score ( Supplementary Table S3 and Supplementary Figure S4 ). Then, the functional enrichment analysis observed that immune-related biological processes and malignant hallmarks were enriched ( Figures 5A–D ). The upregulated genes showed significantly enriched GO terms, including response to interleukin-6, mitotic cell cycle phase transitions, and KEGG terms, including transcriptional misregulation in cancer, p53 signaling pathway, and cell cycle ( Figures 6A, B ). GO terms, including response to chemokine, regulation of inflammatory response, and leukocyte proliferation, and KEGG terms, including nuclear factor (NF)-kappa B signaling pathway, hematopoietic cell lineage, and cytokine–cytokine receptor interaction, were significantly enriched in downregulated genes ( Figures 6C, D ). Furthermore, GSEA found significantly enriched terms related to tumorigenesis and drug resistance in high-risk group ( Supplementary Table S4 ), including spliceosome (NES = 2.354, p < 0.0001, FDR < 0.0001, Figure 6E ) and proteasome (NES = 1.881, p = 0.025, FDR = 0.039, Figure 6F ) in KEGG gene set; MYC targets (NES = 2.304, p < 0.0001, FDR < 0.0001, Figure 6G ), E2F targets (NES = 1.794, p = 0.019, FDR = 0.029, Figure 6H ), and unfold protein response (NES = 1.785, p = 0.024, FDR = 0.023, Figure 6I ) in HALLMARK gene set; ribonucleoprotein complex biogenesis (NES = 2.337, p < 0.0001, FDR = 0.001, Figure 6J ), spliceosomal complex assembly (NES = 2.324, p = 0.002, FDR = 0.006, Figure 6K ), and ribosome biogenesis (NES = 2.219, p < 0.0001, FDR = 0.010, Figure 6L ) in GOBP gene set.

Given that considerable immune-related terms were substantially enriched in DEGs based on the risk score, we speculated that the risk score was correlated with immune infiltration level. To verify our hypothesis, we assessed the immune infiltration level by the ESTIMATE algorithm. Then, MM patients were divided into high- and low-score groups according to the median value. Kaplan–Meier survival analysis revealed that MM patients with higher immune score and lower tumor purity had longer OS (p = 0.01673 and 0.00896, respectively, Figures 7A–C ). Patients with high risk score in the training cohort ( Figures 7D–F ) and validation cohort ( Figures 7G–I ) had lower immune score and higher tumor purity, which indicated poorer prognosis. Furthermore, the implementation of ssGSEA was used for identifying the alteration of immune cell population upon the risk score. We observed that 16 of 28 immune cell subsets were significantly decreased in the high-risk group ( Figure 7J ). Antitumor immune cells, such as effector memory CD8 T cell, natural killer cell, natural killer T cell, type 1 T helper cell, and type 17 helper cell, were substantially reduced. However, immunosuppressive cells, including MDSC, regulatory T cell, and immature dendritic cell were also decreased in patients with high risk score.

Next, we studied the expression of HLA system and immune checkpoints, including CD274 (PD-L1), PDCD1LG2, VTCN1 (v-set domain-containing T cell activation inhibitor), CD276, and IDO1 (indoleamine 2,3-dioxygenase 1) (37), which were found on the surfaces of cells to help the recognition of immune system and involved in the immune activation and inhibition. We observed that the HLA family, including HLA-DMA, HLA-DMB, HLA-DOA, HLA-DOB, HLA-DPA1, HLA-DPB1, HLA-DPB2, HLA-DQA1, HLA-DQB1, HLD-DRA, HLA-DRB1, HLA-DRB6, and HLA-E, showed significant decrease in high-risk group ( Figure 8A ). Immune checkpoint markers, including CD274 and VTCN1, were significantly downregulated in patients with high risk score ( Figure 8B ).