We used expression data of 229 MM patients and 20 healthy adults donors from GSE6477 and GSE13591 datasets. First, we used the robust multiarray standard (RMA) method to perform background correction, quantile normalization, and log transition for raw data (Supplementary Figure 1). As shown in Figure 1A, volcano plots were used to show the distribution of differential expressed genes between MM patients and healthy donors from GSE6477 to GSE13591. Red dots represent significantly upregulated genes (p < 0.05, —log₂FC— > 1), and blue dots represent downregulated genes (p < 0.05, —log₂FC— < −1). Heat maps showed the top 40 differential expressed genes in these two GEO datasets (Figure 1B). Each column denoted a MM or donor sample, and each row represented a significant expressed gene. We totally identified 1,017 downregulated genes and 578 upregulated genes from GSE6477 dataset. Again, we obtained 137 downregulated genes and 84 upregulated genes from GSE13591 dataset. Among them, we obtained common DEGs from the 2 datasets including 92 downregulated and 31 upregulated genes (Figure 1C).

To further explore the biological function of DEGs identified above, we used GO enrichment analysis using DAVID bioinformatics resources (Huang et al., 2009). The up- or downregulated genes-associated functional categories are shown in Table 1. As shown in Figure 2A, the top 5 biological processes (BPs) contains “immune response,” “phagocytosis, engulfment,” “positive regulation of B cell activation,” “B cell receptor signaling pathway,” and “phagocytosis, recognition.” The top 5 molecule functions (MFs) includes “immunoglobulin receptor binding,” “MHC class II receptor activity,” “antigen binding,” “peptide antigen binding,” and “protein heterodimerization activity.” In the cellular components (CCs) group, the top 5 CCs contains “external side of plasma membrane,” “blood microparticle,” “MHC class II protein complex,” “immunoglobulin complex, circulating,” and “transport vesicle membrane.” Additionally, the bubble map of the entire DEGs was drawn using the ggplot2.R package for the significant biological processes based on p-value, and the results indicated that the DEGs were closely associated with immune-related biological processes, such as “positive regulation of B cell activation,” “MHC class II receptor activity,” “immunoglobulin receptor binding,” and “immunoglobulin complex, circulating” (Figure 2B).

We first analyzed the correlation between DEGs and KEGG processes and found that they were significantly enriched in “systemic lupus erythematosus,” “influenza A,” “antigen processing and presentation,” “Th cell differentiation,” “hematopoietic cell lineage,” “intestinal immune network for IgA production,” and “graft-vs.-host disease” (Figure 3A). Then, we conducted the Reactome pathway enrichment analyses, and the results showed that the enrichment of downregulated DEGs were mainly in the immune system, adaptive immune system, antigen processing and presentation, major histocompatibility complex (MHC) class II antigen presentation, translocation of ZAP-70 to immunological synapse, phosphorylation of CD3 and T-cell receptor (TCR) zeta chains and graft-vs.-host disease, while upregulated DEGs only contained a small amount of genes and were mostly enriched in senescence-associated biological events (Figure 3B and Table 2). Taken together, the downregulated DEGs from the GSE6477 to GSE13591 GEO datasets play an important role in MM progression.

The above data indicated that the downregulated DEGs were significantly related with the immune-associated progression. Thus, we analyzed these downregulated DEGs by evaluating the immune system process in ClueGo. As shown in Figure 4A, downregulated DEGs were mainly enriched in antigen processing and presentation of exogenous peptide antigen via MHC class II, B-cell receptor signaling pathway, positive regulation of lymphocyte activation, peptide antigen assembly with MHC class II protein complex, and regulation of lymphocyte activation. Furthermore, we also conducted an interrelation analysis between pathways in the BPs of downregulated DEGs. As expected, we observed an enrichment of the downregulated DEGs largely in antigen processing and presentation of exogenous peptide antigen via MHC class II, peptide antigen assembly with MHC class II protein complex, regulation of B-cell receptor signaling pathway, positive regulation of immune system process, and B-cell differentiation (Figure 4B).

STRING is an online analysis tool, which is often used to study and integrate interaction between proteins (Szklarczyk et al., 2017). Based on the DEGs identified above, we constructed the protein–protein interaction (PPI) network of DEGs expression in MM. After removing the separated nodes, we described a complex PPI network of entire DEGs, and the final PPI network contained 68 nodes and 222 edges (Figure 5A). Subsequently, we utilized cluster analysis of PPI network in Cytotype MCODE; we identified two significant modules: module 1 and module 2. Module 1 contained 14 nodes and 72 edges (Figure 5B); the DEGs in module 1 were mainly enriched in adaptive immune system, immune system, and MHC class II antigen presentation (Table 3). Inversely, module 2 includes 8 nodes and 13 edges (Figure 5C); the DEGs in module 2 were not relevant for immune system (Table 4).

GSEA is a strong analytical method to identify the differential genes enrichment in several biological events (Saxena et al., 2006). Specifically, in GSE6477 cohort, the low expression of CD79A, CD19, and CD81 were in direct proportion to the B-cell receptor signaling pathway. The downregulated HLA-DRA, HLA-DRB4, HLA-DRB1, and CD19 genes were enriched in response to hematopoietic cell lineage (Figure 5D). In another GSE13591 dataset, we found that 12 downregulated genes were positively correlated with antigen processing and presentation, including HLA-DRB4, CTSB, HLA-DQB1, HLA-DRB1, HLA-DRA, HLA-DPB1, LGMN, HSPA1B, HLA-DQA1, HSPA1A, HLA-DMA, and HLA-DPA1 (Figure 5E).

The DEGs mentioned in the PPI network analysis were further scored in top 10 by three methods in CytoHubba, which were selected as key hub genes of MM (Figure 6A). The common genes included CD19, CTSB, FBL, GAPDH, RRM2, and TXN, which may play a vital role in MM progression (Figure 6B). Subsequently, our Kaplan–Meier analysis showed that CD19, CTSB, and FBL were not relevant to the OS of MM patients (Figures 6C–E, p > 0.05). On the contrary, GAPDH, RRM2, and TXN were closely associated with MM survival (Figures 6F–H, p < 0.05). Among them, the housekeeping gene GAPDH did not participate in MM progression. As a major player in the antioxidant pathway, TXN probably played a critical role in reducing oxidized cysteine residues and cleaving disulfide bonds (Cao et al., 2020; Yang et al., 2020). Surprisingly, it has been reported that silencing RRM2 inhibits multiple myeloma by targeting the Wnt/β−catenin signaling pathway (Liu et al., 2019), implying that RRM2 may be a new therapeutic target for MM. Therefore, we conducted various experiments to evaluate the biological functions by commercial specific RRM2 inhibitor osalmid in MM cell lines.

We further evaluated the mRNA expression level of RRM2 between MM patients and normal donors using the Oncomine database; the results showed a significant elevation of mRNA level of RRM2 in MM patients (Figure 7A). As a novel identified RRM2 inhibitor, osalmid has been approved for treating cholecystitis, inflammation, and postcholecystectomy syndrome, which inhibited RRM2 activity by binding to hydrogen bond of RRM2 (Liu et al., 2016). Herein, we used CCK8 cell viability assay to explore the effect of osalmid on MM cells; the results showed that osalmid inhibited RPMI-8226 and U266 cell viabilities in a concentration-dependent manner (Figure 7B). RRM2 is a member of ribonucleotide reductases, which is required for DNA synthesis. Thus, we further evaluated the impact of osalmid on DNA damage. Our data confirmed that RRM2 inhibitor osalmid could activate the expression of γH2Ax, which is a hallmark of DNA damage (Figure 7C). Subsequently, flow cytometry assay indicated that RRM2 inhibitor osalmid induced MM cell cycle S phase arrest in a dose-dependent manner (Figures 7D–F). Although osalmid could also inhibit normal cells, it was not relevant to the OS of acute myeloid leukemia (AML), diffuse large B-cell lymphoma (DLBCL), and follicular lymphoma patients (Supplementary Figures 2, 3). Totally, RRM2 could be an important experimental new agent for treatment of MM.

GSE6477 and GSE13591 were obtained from GEO¹ database, which were measured in the GPL96 platform (Affymetrix Human Genome U133A Array). We further analyzed the DEGs using online GEO2R software, and p < 0.05 and —log₂FC— > 1 were used as cutoff criteria.

GO term was used to analyze the biological events including biological process, cellular component, and molecular function by DAVID online biological information database (Huang et al., 2007). Signaling pathway analysis was carried out using another online database, REACTOME. p < 0.05 was considered as statistically significant and to achieve significant enrichment.

The STRING database contained known and predicated protein–protein association data for many organisms, including Homo sapiens (Szklarczyk et al., 2017). The Cytoscape software V3.5.1 was used to show the PPI network and the most significant module based on the MCODE score and node number (Shannon et al., 2003).

GSEA is a computational method to evaluate whether a set of genes show the significant difference in various biological events (Subramanian et al., 2007). The analysis was conducted with default settings. p < 0.05 and false discovery rate (FDR) q < 0.25 were considered as statistically significant.

Indicated cells were plated at 1 × 10⁴ cells per 96-well and treated with different concentrations of osalmid. The cell survival rate was detected using CCK8 kit according to the protocol (KeyGEN, Nanjing, China).

Indicated cells with different dose of osalmid were lysed in the radioimmunoprecipitation assay (RIPA) buffer, and the detailed procedure was conducted according to the previous report (Yao et al., 2017). The following antibodies were used: β-actin and γH2Ax (Proteintech, Wuhan, China).

The distribution of indicated cells were analyzed by flow cytometry. The experiments were performed according to the procedure described previously (Yao et al., 2018).

The results are presented as mean ± SD of three independent experiments. The statistical significance was assessed by Student’s t-test with #, ^∗p < 0.05 and ##, ^∗∗p < 0.01.