Identification of Potential Prognostic Biomarker for Predicting Survival in Multiple Myeloma Using Bioinformatics Analysis and Experiments

Multiple myeloma (MM) is a malignant disease of plasma cells, which remains incurable because of its unclear mechanism and drug resistance. Herein, we aimed to explore new biomarkers and therapeutic targets in MM. After screening differentially expressed genes (DEGs) in GSE6477 and GSE13591 dataset, we performed Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses of DEGs using DAVID online database. The results indicated that the downregulated DEGs were mainly enriched in the immune-associated biological process. The protein–protein interaction network was constructed by STRING database, on which we performed module analysis and identified key genes. Gene set enrichment analysis (GSEA) and Kaplan–Meier analysis showed that RRM2 could be a novel biomarker in MM diagnosis. We further confirmed that novel RRM2 inhibitor osalmid inhibited MM cell proliferation and triggered cell cycle S phase arrest. Targeting RRM2 was expected to develop new therapeutic strategies for malignant MM.

Following the improvement of microarray and RNAsequencing technology, many gene expression studies about MM have been reported in the past few years. In 2016, the Gene Expression Omnibus (GEO) database is built and maintained by the National Center for Biotechnology Information (NCBI) (Clough and Barrett, 2016), which collects lots of highthroughput microarray, chips submitted by researchers from all over the world. The differentially expressed genes (DEGs) Abbreviations: MM, multiple myeloma; GEO, gene expression omnibus; DEGs, differentially expressed genes; GO, gene ontology; BPs, biological processes; MFs, molecule functions; CCs, cellular components; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes; MCODE, Molecular Complex Detection; GSEA, gene set enrichment analysis; IMiDs, immunomodulatory drugs; RRM2, ribonucleotide reductase M2; OS, overall survival. between MM patients and healthy donors are screened using the GEO database; then, the biological process and gene regulatory networks are further analyzed (Zhou et al., 2015;Yan et al., 2019). Eventually, to verify the clinical significance of these critical genes, the overall survival (OS) correlation between critical genes and MM patients is explored in the website of Kaplan-Meier plotter.
In the present study, an integrated bioinformatics analysis was performed to pursue the DEGs in MM patients compared with healthy adult donors. First, we downloaded two messenger RNA (mRNA) gene expression profiles [GSE6477 (Chng et al., 2007) and GSE13591 (Agnelli et al., 2009)] from the GEO database and identified the common DEGs from the two microarrays. We further performed the gene ontology (GO) function, KEGG pathway enrichment, and protein-protein interactions (PPI)

Identification of DEGs in Two GEO Datasets
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 2 FC| > 1), and blue dots represent downregulated genes (p < 0.05, |log 2 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).

GO Functional Analysis of DEGs
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).

Signaling Pathway Enrichment Analysis of DEGs
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 Interrelation analysis between pathways in the BPs of downregulated DEGs. Each node represents a GO term, and lines reflect correlations between these terms. The node size is positively associated with term/pathway significance. The node color reflects the functional enrichment groups, with similar function clustering in same color. (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 Immune-Associated Functional Enrichment Analysis of Downregulated DEGs 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). Frontiers in Genetics | www.frontiersin.org

PPI and GSEA Analysis of Immune-Related Function of Hub Genes
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).

Identification of Highly Expressed Hub Genes by CytoHubba Methods
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 , 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.

RRM2 Inhibitor Osalmid Showed Antitumor Activity
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.

DISCUSSION
Multiple myeloma incidence has been increasing over the past few decades, making it a disease of considerable clinical and economic impact (Ferlay et al., 2015). In this study, we screened out DEGs based on GSE6477 and GSE13591 datasets. The upregulated DEGs were strikingly enriched in senescence-associated biological events, and the downregulated DEGs in immune system, adaptive immune system, antigen processing and presentation, MHC class II antigen presentation, translocation of ZAP-70 to immunological synapse, phosphorylation of CD3 and TCR zeta chains, and graft-vs.-host disease. It is well known that immunomodulation plays an important role in MM progression, implying that the downregulated DEGs could be new targets in developing novel immunomodulatory drugs.
In recent years, immunomodulatory drugs (IMiDs) are broadly applied in the treatment of multiple myeloma. At present, three IMiDs, namely, thalidomide, lenalidomide, and pomalidomide, have been approved for the treatment of MM (Singhal et al., 1999). Cereblon (CRBN), a common target for IMiDs, was part of E3 ubiquitin ligase complex. IMiD drugs bound to CRBN and triggered a change in CRBN targets, ultimately initiating therapeutic activity (Zhu et al., 2011;Chamberlain et al., 2014). The resistant mechanism to IMiDs was due to downregulation of CRBN expression or CRBN mutations (Kortum et al., 2016). In our study, we identified several hub genes using gene set enrichment analysis. Meanwhile, we conducted Kaplan-Meier analysis; lower expression levels of HLA-DPB1, HSPA1A, HSPA1B, and TNFRSF4 showed worse overall survival (OS) in MM patients.
Bioinformatics is a data-driven branch of science that has been commonly used for high-through data analysis and involves a large number of powerful analysis software packages and database (Xia, 2017). We further identified top 10 genes using Betweenness, Bottleneck, and Stress methods in Cytohubba. Kaplan-Meier analysis indicated that GAPDH, RRM2, and TXN were closely associated with MM survival (p < 0.05). Among them, RRM2 has a crucial role in mediating cancer cell survival and progression (Liu et al., 2016;Tang et al., 2020). Therefore, we used a new identified RRM2 inhibitor osalmid to silence RRM2 expression and evaluate its role in MM progression. As expected, silencing RRM2 by osalmid treatment led to cell viability inhibition and cell cycle arrest. Taken together, targeting RRM2 may provide a novel therapeutic strategy for the treatment of multiple myeloma.

CONCLUSION
In summary, integrated bioinformatics analysis of datasets of MM patients and healthy adult donors was performed. We identified several genes associated with immune-related biological processes, such as B-cell receptor signaling pathway, cytokine-cytokine receptor interaction, hematopoietic cell lineage, and antigen processing and presentation. Among them, GAPDH, RRM2, and TXN genes were all upregulated, and Kaplan-Meier analysis showed that these genes were closely associated with worse overall survival (OS) of MM. Gene enrichment analysis and biological experiments indicated that RRM2 could be a potential biomarker in MM progression.

GEO Datasets Selection and Data Processing
GSE6477 and GSE13591 were obtained from GEO 1 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 2 FC| > 1 were used as cutoff criteria.

GO and KEGG Pathway Enrichment Analysis of DEGs
GO term was used to analyze the biological events including biological process, cellular component, and molecular function by DAVID online biological information database 1 https://www.ncbi.nlm.nih.gov/geo/ (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.

PPI Network Analysis of DEGs
The STRING database contained known and predicated proteinprotein 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).

Gene Set Enrichment Analysis
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.

Cell Viability Assay
Indicated cells were plated at 1 × 10 4 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).

Western Blot
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).

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

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

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.