Comprehensive Analysis of m6A RNA Methylation Regulators in the Prognosis and Immune Microenvironment of Multiple Myeloma

Background N6-methyladenosine is the most abundant RNA modification, which plays a prominent role in various biology processes, including tumorigenesis and immune regulation. Multiple myeloma (MM) is the second most frequent hematological malignancy. Materials and Methods Twenty-two m6A RNA methylation regulators were analyzed between MM patients and normal samples. Kaplan–Meier survival analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis were employed to construct the risk signature model. Receiver operation characteristic (ROC) curves were used to verify the prognostic and diagnostic efficiency. Immune infiltration level was evaluated by ESTIMATE algorithm and immune-related single-sample gene set enrichment analysis (ssGSEA). Results High expression of HNRNPC, HNRNPA2B1, and YTHDF2 and low expression of ZC3H13 were associated with poor survival. Based on these four genes, a prognostic risk signature model was established. Multivariate Cox regression analysis demonstrated that the risk score was an independent prognostic factor of MM. Enrichment analysis showed that cell cycle, immune response, MYC, proteasome, and unfold protein reaction were enriched in high-risk MM patients. Furthermore, patients with higher risk score exhibited lower immune scores and lower immune infiltration level. Conclusion The m6A-based prognostic risk score accurately and robustly predicts the survival of MM patients and is associated with the immune infiltration level, which complements current prediction models and enhances our cognition of immune infiltration.


INTRODUCTION
Multiple myeloma (MM) is a neoplastic hematological malignancy characterized by clonal proliferation of malignant plasm cells in bone marrow (BM). MM has been the second most frequent hematological malignancy and accounts for 1.79% of all new cancer cases and 2.11% of all cancer deaths worldwide (1). The treatment arsenal in the battle of MM began with chemotherapy, and was rediscovered with proteasome inhibitor (PI)-based triple-drug combination, and extended with the splendid results achieved by immunotherapy, such as monoclonal antibodies (mAbs) and chimeric antigen receptorengineered T cells (CAR-T). The survival of MM patients has dramatically prolonged in the last decade. Nonetheless, the medical needs of MM patients remain unmet, especially highrisk patients (2). Two leading obstacles for MM treatment are high recurrence rate and highly heterogeneous cell population whose subclone content evolves over time. Not only are genotypes and phenotypes of MM responsible for the high tumor heterogeneity but also the tumor microenvironment is implicated in it. Similar to what is observed in chemotherapy, MM presents capability to escape from immunotherapy, which is associated with high genetic instability and the protection of bone marrow microenvironment (BMME) corrupted by MM cells (3). Considerable research emphasizing the role of BMME in hematological cancers explicitly demonstrated that cellular components of bone marrow, including immune cell and stromal cell, can support the proliferation of MM cells and shield MM cells from the attack of chemotherapy and immune system (4). The alteration of immune infiltration level in BMME involves the development and recruitment of multiple immunosuppressive cells, such as regulatory T cells (Treg), myeloid-derived suppressor cells (MDSCs) and tumor-associated macrophages (TAMs) (5). In addition, the dysregulation of immune checkpoint genes, such as programmed cell death 1 ligand1/L2 and TIGIT, and the loss or downregulation of human leukocyte antigen (HLA) are also involved in the immune-suppressive status of BMME (6)(7)(8).
In the present study, we used the gene expression data and clinical information from the Gene Expression Omnibus (GEO) database to explore differentially expressed and prognosis-related m6A regulators. Then, we established a four-gene prognostic risk signature model using least absolute shrinkage and selection operator (LASSO) Cox regression analysis to evaluate the survival outcomes of MM patients. Finally, enrichment analysis and the evaluation of immune infiltration level were performed. Our study systematically dissected the abnormal expression status of m6A-related genes and revealed the role of m6A RNA methylation modification in the prognosis and immune infiltration of MM, which complements current prediction models and enhances our cognition of immune infiltration.

Data Preparation
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.

Protein-Protein Interaction Network and Spearman Correlation Analysis
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.

Construction of the Prognostic Risk Model
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.

Establishment and Assessment of the Nomogram
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.

Differential Expression Analysis
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).

Evaluation of Immune infiltration Level
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 immuneassociated 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
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 posthoc 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.

PPI Network and Spearman Correlation Analysis Among m6A RNA Methylation Regulators
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).

Construction and Validation of the Prognostic Risk Model Based on Four m6A Methylation Regulators
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-  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 The PPI network of the differentially expressed m6A-related genes. Yellow, green, and blue nodes represent "writer," "eraser," and "reader," respectively. Node size represents the expression value. Edge size represents the combined score (low score to orange and high score to green). (D) The Spearman correlation analysis of the differentially expressed m6A-related genes. The size of each circle denotes the Spearman correlation coefficient (negativity to red and positivity to blue). *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. . 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%   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.

Subgroup Analysis
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.

Correlation Between the Prognostic Risk Score and Clinical Characteristics
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.

Association of the Risk Score With the Immune Microenvironment
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 highrisk group ( Figure 8A). Immune checkpoint markers, including CD274 and VTCN1, were significantly downregulated in patients with high risk score ( Figure 8B).

DISCUSSION
N6-methyladenosine RNA modification, the most remarkable epitranscriptomic modification, has been widely considered to be involved in tumorigenesis and immune regulation, while its role in MM progression and bone marrow immune microenvironment is little known yet. In this study, we conducted the comprehensive bioinformatics analysis to explore the significance of m6A RNA methylation regulators in MM. We identified the aberrant expression status of m6A RNA methylation regulators in MM: RBMX and HNRNPC were significantly upregulated, while METTL3, METTL14, METTL16, ZC3H13, KIAA1429, FTO, and IGF2BP3 were significantly downregulated in MM patients compared to normal plasm samples. Furthermore, our results elucidated that the overexpression of HNRNPA2B1, HNRNPC, and YTHDF2, and the downregulation of ZC3H13 were negatively associated with the overall survival. Based on these four genes, a prognostic risk signature model was established via LASSO Cox regression analysis. Kaplan-Meier survival analysis showed that MM patients in the high-risk group had significantly shorter OS than patients in low-risk group. Two validation cohorts confirmed the predictive efficiency of the m6A risk score by ROC curve analysis. In addition, univariate and multivariate Cox regression analyses manifested that the risk score was an independent prognostic factor. To determine the clinical feasibility of the prognostic risk signature in MM, we assessed the correlation between the risk score and significant clinical characteristics. Subgroup survival analysis upon age, gender, ISS stage, cytogenetic abnormality, and different therapeutic strategies showed that the risk score was precise and robust for the individualized medical demands of MM patients. Additionally, we found that advanced ISS stage and the existence of cytogenetic abnormalities were positively associated with the risk score.
Consistent with published articles, HNRNPC can serve as an oncogene, which is overexpressed in various cancers. HNRNPC knockdown arrested the proliferation of breast cancer cell lines via activating the antitumor cascade of interferon response initiated by the cytoplasmic RNA sensors retinoic acid-inducible gene I (RIG-I) (38). ZC3H13 is involved in the assembly of m6A methyltransferase complex. In colorectal cancer, ZC3H13 acting as a tumor suppressor could suppress cancel cells proliferation and invasion via inactivating Ras-ERK signaling pathway (39). YTHDF2, as a member of m6A "reader," can selectively recognize and bind to m6A-containing RNA to initiate its degradation. The overexpression of YTHDF2 in hepatocellular carcinoma (HCC) promoted stemness and metastasis via increasing OCT4 translation (40). In glioblastoma (GBM) stem cells, YTHDF2 was preferentially expressed, and it enhanced GBM growth via YTHDF2-MYC-IGFBP3 axis (41). HNRNPA2B1 was highly expressed in multiple types of cancer. HNRNPA2B1-deficient breast cancer cells showed inhibited growth and increased apoptosis rate via dampening the phosphorylation of STAT3 (42  Then, to figure out the underlying mechanisms of the risk score, we conducted the differential expression analysis between high-and low-risk groups. Enrichment analysis including overrepresentation analysis (ORA) and GSEA was subsequently conducted. The ORA results showed that oncogenesis-related (J) Twenty-eight immune-related gene sets were performed for ssGSEA. Seventeen immune cell gene sets were significantly different between patients with high and low risk score. Red star represents higher enrichment score; blue star represents lower enrichment score compared to high-risk MM patients. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
terms, including mitotic cell cycle phase transition, transcriptional misregulation in cancer, and cell cycle were significantly enriched in upregulated DEG, while immune-related terms, such as response to chemokine, regulation of leukocyte activation, regulation of inflammatory response, regulation of immune effector process, leukocyte proliferation, and migration, were substantially enriched in downregulated DEGs. GSEA results indicated that pathways correlated with tumorigenesis and proteasome inhibitor resistance (MYC targets, unfold protein response, and proteasome) were positively correlated with the high risk score. Since considerable immune-related processes were enriched in downregulated DEGs, we assumed that the m6A risk signature was correlated with the immunosuppressive microenvironment. Thus, we evaluated the immune infiltration level by ESTIMATE algorithm and ssGSEA. As expected, patients in the high-risk group had lower immune score. Sixteen of 28 immune cell subsets, including antitumor immune cells (effector memory CD8 T cell, natural killer cell, natural killer T cell, type 1 T helper cell, and type 17 helper cell) and protumor immune cells (MDSC, regulatory T cell and immature dendritic cell), were significantly reduced in the high-risk group. Moreover, the expression of HLA genes, as cell surface markers for the recognition of immune system, were generally decreased in MM patients with high risk score. Inhibitory immune checkpoint molecules, involved in the immune evasion, PD-L1, and VTCN1, were also significantly downregulated in patients with high risk score. Taken together, we found that the risk score is correlated with the immunosuppressive bone marrow microenvironment from cellular level to molecular level.
The heterogeneity of the tumor immune microenvironment can result in diverse dimensions, such as prognosis and therapeutic response to immunotherapies. However, the underlying mechanism by which tumor cells "reprogrammed" immune microenvironment into the culprit for oncogenesis and drug resistance remains unclear. It was not until 2017 that m6A RNA modification was first elucidated to be involved in the regulation of mammal immune system. METTL3/METTL14deficient T cells failed to maintain homeostasis and remained in the naive state through IL-7/STAT5/COCS pathway (43). Su et al. reported that forced expression of FTO contributed to the elevated expression of immune checkpoint genes such as PD-L1, PD-L2, and leukocyte immunoglobulin-like receptor subfamily B4 (LILRB4) via inhibiting YTHDF2-dependent degradation in AML cell lines (44). Deletion of the demethylase ALKBH5 inhibited the recruitment of immunosuppressive cells (Treg and MDSC) through targeting the key enzyme of lactate transport MCT4/SLC16A3 and decreasing the extracellular lactate concentration, which subsequently resulted in the sensitization of anti-PD-1 therapy in melanoma cells (45). Zhang et al. comprehensively evaluated the m6A modification in gastric cancer and established the m6Ascore to quantify m6A modification patterns for individual patients. They observed that low m6Ascore was characterized by the immune activation with increased efficacy of immune checkpoint blockade treatment, while high m6Ascore was characterized by stroma activation (46). Jin et al. constructed an eight-gene m6A risk signature model in adrenocortical carcinoma to evaluate the prognosis and tumor immune microenvironment (47). Our study also revealed that the m6A-based risk score can accurately predict survival outcomes and evaluate the immune infiltration level.
However, our study does have some limitations. First, the majority of patients included in our study were white; the prognostic value of m6A-related genes remains to be further verified in Asian and African cohorts. Second, given that our results are based on the mRNA expression data, the expression at protein level should be further clarified. Third, some important clinical information, including R-ISS stage, mSMART risk stratification, somatic mutation, and copy number variation, was not available in the GEO database, by which in-depth analysis was inaccessible. Fourth, since the lack of significant clinicopathological information in two validation cohorts, the extrapolation of our nomogram result was not further validated. Thus, experiment verification and clinical investigation with innovative design and orchestrated implement are needed to translate these descriptive results into clinical benefits.
In conclusion, the four m6A regulators (ZC3H13, HNRNPA2B1, HNRNPC, and ZC3H13)-based prognostic risk score can accurately and stably predict the survival of MM patients; and the risk score is tightly associated with the impaired immune infiltration level, which complements current prediction models. The m6A-based prognostic risk score enhances our cognition of immune infiltration and addresses the indispensable role of RNA methylation to cause the heterogeneous tumor microenvironment. Further research is required to understand the overwhelmingly complex immune regulation network in bone marrow microenvironment.

DATA AVAILABILITY STATEMENT
The public datasets used in this article can be downloaded in the GEO database. The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
RL: conceptualization, data curation, methodology, software, formal analysis, writing-original draft, and visualization. YS: conceptualization, writing-review and editing, supervision, and funding acquisition. JH: conceptualization, formal analysis, methodology, supervision, validation, and writing-review and editing. XW: validation and writing-review and editing. DW: methodology, software, visualization, and writing-review and editing. MZ: conceptualization, methodology, software, validation, visualization, and writing-review and editing. AH: conceptualization, project administration, validation, writingreview and editing, and supervision. JB: conceptualization, project administration, validation, writing-review and editing, and supervision. All authors contributed to the article and approved the submitted version.