RNA-Sequencing-Based Transcriptomic Score with Prognostic and Theranostic Values in Multiple Myeloma

Multiple myeloma (MM) is the second most frequent hematological cancer and is characterized by the clonal proliferation of malignant plasma cells. Genome-wide expression profiling (GEP) analysis with DNA microarrays has emerged as a powerful tool for biomedical research, generating a huge amount of data. Microarray analyses have improved our understanding of MM disease and have led to important clinical applications. In MM, GEP has been used to stratify patients, define risk, identify therapeutic targets, predict treatment response, and understand drug resistance. In this study, we built a gene risk score for 267 genes using RNA-seq data that demonstrated a prognostic value in two independent cohorts (n = 674 and n = 76) of newly diagnosed MM patients treated with high-dose Melphalan and autologous stem cell transplantation. High-risk patients were associated with the expression of genes involved in several major pathways implicated in MM pathophysiology, including interferon response, cell proliferation, hypoxia, IL-6 signaling pathway, stem cell genes, MYC, and epigenetic deregulation. The RNA-seq-based risk score was correlated with specific MM somatic mutation profiles and responses to targeted treatment including EZH2, MELK, TOPK/PBK, and Aurora kinase inhibitors, outlining potential utility for precision medicine strategies in MM.


Introduction
Multiple myeloma (MM) is an incurable malignant plasma cell disorder characterized by strong genetic heterogeneity impacting treatment response and clinical outcomes [1][2][3]. Pre-clinical studies have identified a compendium of mechanisms associated with MM cell treatment resistance. Given this heterogeneity, one of the current challenges is to precisely predict survival and treatment response according to patient molecular characteristics in order to develop personalized medicine. Gene expression profiling (GEP) has created major insights in both the understanding and clinical management of this disease. Indeed, analysis of patient transcriptomic data at diagnosis has highlighted substantial molecular heterogeneity between patients that is characterized by distinct gene signatures and

Statistical Analyses
All of the statistical analyses were performed with the R statistics software (version 3.2.3; available from https://www.r-project.org) (28 September 2021) and with R packages developed by the BioConductor project (available from https://www.bioconductor.org/) (28 September 2021) [24]. The expression level of each gene was summarized and normalized using the DESeq2 R/Bioconductor package [25]. Differential expression analysis was performed using the DESeq2 pipeline [25]. p-values were adjusted to control the global FDR across all comparisons with the default option of the DESeq2 package. Genes were considered differentially expressed with an adjusted p-value < 0.05 and a fold change >1.5. For the Montpellier cohort, Affymetrix U133P chips were also used, as previously described [22,26], to analyze GEP and to calculate previously published risk scores, including the RS score [8], UAMS HRS score [5], and IFM score [7].

EZH2 Inhibition in Primary MM Cells from Patients
Mononuclear cells from tumor samples of seven patients with previously untreated MM from the Montpellier cohort were cultured for 8 days in the presence of IL-6 (2 ng/mL) with and without 1 µM EPZ-6438. At day 8 of the culture period, the count of viable MMC was determined using CD138 staining by flow cytometry.

RNA-Seq-Based Gene Risk Score in Multiple Myeloma
We used the gene expression profiling (GEP) data of 674 newly diagnosed MM patients from the Multiple Myeloma Research Foundation's (MMRF) CoMMpass study. We also performed the RNA sequencing of purified MMC from 76 newly diagnosed MM patients treated by high dose therapy (HDT) and autologous stem cell transplantation. Using the Maxstat R function and the Benjamini-Hochberg multiple testing correction, 267 genes were found to have a prognostic value for overall survival (OS) (adjusted p value < 0.05) in the CoMMpass cohort of patients with previously untreated MM. These 267 prognostic genes comprised 142 genes associated with a poor outcome in MM and 97 genes associated with high expression related to a good prognostic value. The 267 genes were used to build an RNA-seq-based risk score. The RNA-seq-based risk score is defined by the sum of the beta coefficient derived from the Cox model for each prognostic gene weighted by −1 or +1 according to the MMC gene expression above or below the Maxstat defined cutpoint [29]. Patients from the CoMMpass cohort were ranked according to increased RNA-seq-based risk scores, and the Maxstat algorithm was used to find the cutoff associated with the maximum difference in the OS ( Figure 1A). The RNA-seq-based risk score split patients into a high-risk group (22.8%) and a low-risk group (77.2%) in the CoMMpass cohort (p-value = 1.7 × 10 −46 ) ( Figure 1B). The prognostic value of the RNA-seq-based risk score was validated in the Montpellier cohort (p-value = 2.8 × 10 −11 ) ( Figure 1C). When applied to a non-disease dataset, i.e., the GEP profiles from a normal B to plasma cell differentiation model, significantly higher score values were identified in pre-plasmablast, plasmablast, and plasma cell stages compared to in the memory B cell stage ( Figure 1D). These data demonstrate that the RNA-seq-based risk score is not a hallmark of proliferation since no significant difference was found between the proliferating pre-plasmablasts and the nonproliferating plasma cells ( Figure 1D). Furthermore, the RNA-seq-based risk score was not correlated with plasma cell labeling index (PCLI-percentage of MM cells in S phase [30]) at diagnosis in the Montpellier cohort (Supplementary Figure S1). Significantly higher score values were also observed in the HMCLs compared to in the MM cells of patients, indicating that a high score is associated with MM progression ( Figure 1D). Altogether, these data highlight that the RNA-seq-based risk score can identify newly diagnosed high-risk MM patients in independent cohorts.

High-Risk MM Patients
Identified with the RNA-Seq-Based Risk Score Are Characterized by Enrichment of Genes Related to Cell Proliferation, Growth Factor Signaling, MYC Pathway and Epigenetic Deregulation GSEA [31,32] analyses revealed that the genes associated with a poor prognostic value were significantly enriched in the genes related to interferon response, cell proliferation, hypoxia, IL-6 signaling pathway, interferon (IFN) response, stem cell genes, MYC, and epigenetic deregulation (Figure 2A). Among the epigenetics-related genes enriched in high-risk MM patients, EZH2 targets, HDAC targets, and DNA methylation target genes were identified. We next investigated the RNA-seq-based risk score value distribution according to the Affymetrix GEP-based risk scores described in MM. The RNA-seq-based risk score values were significantly higher in patients who were defined as high-risk based on their RS score [8], UAMS HRS score [5], and IFM score [7]. Furthermore, high-risk RNA-seq-based score patients demonstrated a significant increase in the percentage of proliferating MM cells (median = 0.7%; range: 0-7.3%) compared to the low-risk group (median = 1.55%, range: 0-17.3%) ( Figure 2B).

Association between RNA-Seq-Based Risk Score and Mutations in MM
We then analyzed the relationship between the RNA-seq-based risk score distribution, cytogenetic abnormalities, and mutational status in the CoMMpass cohort. The RNA-seq-based risk score values were significantly higher in MM patients with del(13q), del(17p), del(1p), 1qgain, t(4;14), t (12;14), and t(14,16) groups ( Figure 3A). The RNA-seq-based risk score values of the MM patients were also evaluated according to the status of the most frequent mutations identified in HMCLs [28]. The RNA-seq-based risk scores were significantly higher in the MM patients and were characterized by more ASXL1, ATM, BRAF, DIS3, EP300, FGFR3, KMT2B, LRP1B, MAP3K1, MAX, NOTCH2, NUP214, PRDM1, PTPRD, RB1, ROS1, SETD2, TP53, TRRAP, and ZFHX3 mutations compared to the unmutated MM patients ( Figure 3B). Interestingly, the RNA-seq-based risk scores were significantly higher in MM patients harboring double-hit TP53 mut /del(17p) or FGFR3 mut /t(4;14) compared to MM patients without double-hit ( Figure 3C-D). The prognostic value of the RNA-seq-based risk score was then compared with the molecular subgroups and mutations mentioned above using Cox analyses. In univariate Cox analyses, the RNA-seq-based risk score, del(17p), 1qgain, and t(12;14) subgroups as well as the ATR, CREBBP, MAP3K1, PMS1, and TP53 mutations showed a prognostic value in the CoMMpass cohort (Table 1). When the RNA-seq-based risk score and molecular subgroups were tested together (multivariate analysis), only the RNA-seq-based risk score retained prognostic value ( Table 2). When the RNA-seq-based risk score and mutation were tested together, the score and PMS1 gene mutation remained independent prognostic factors in the CoMMpass cohort (Table 3). CREBBP, MAP3K1, PMS1, and TP53 mutations showed a prognostic value in the CoMMpass cohort (Table 1). When the RNA-seq-based risk score and molecular subgroups were tested together (multivariate analysis), only the RNA-seq-based risk score retained prognostic value (Table 2). When the RNA-seq-based risk score and mutation were tested together, the score and PMS1 gene mutation remained independent prognostic factors in the CoMMpass cohort (Table 3).  . Wilcoxon test. ns: nonsignificant, * p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001, **** p-value < 0.0001.

RNA-Seq Risk Score Revealed New Genes Significantly Associated with MM Pathophysiology
We used public datasets of RNAi [33,34] and CRISPR-Cas9 [35] viability assay screening (Dependency Map data, Broad Institute, www.depmap.org) (28 September 2021) to identify the essential genes in myeloma cell lines compared to other cancer cell lines. Among the 142 genes associated with poor survival composing the RNA-seq-based risk score, we found seven genes (ATP8B1, FGFR4, FOXD4, MX1, NPTXR, TMEM171, and TNFRSF10B) with a significantly lower DEMETER2 score in the myeloma cell lines (n = 16) compared to other cancer cell lines (n = 695) ( Figure 4A and Table 4). The CERES score of four genes (ISG20, NDC1, SF3B3, UMPS) was significantly lower in the myeloma cell lines (n = 20) compared to in the other cancer cell lines (n = 769) ( Figure 4B and Table 4). These essentiality scores were calculated using the RNAi [36] and CRISPR-Cas9 [35] methods, respectively. Thus, these analyses identified essential MM genes that had not previously been considered.

RNA-Seq Risk Score Revealed New Genes Significantly Associated with MM Pathophysiology
We used public datasets of RNAi [33,34] and CRISPR-Cas9 [35] viability assay screening (Dependency Map data, Broad Institute, www.depmap.org) (28 September 2021) to identify the essential genes in myeloma cell lines compared to other cancer cell lines. Among the 142 genes associated with poor survival composing the RNA-seq-based risk score, we found seven genes (ATP8B1, FGFR4, FOXD4, MX1, NPTXR, TMEM171, and TNFRSF10B) with a significantly lower DEMETER2 score in the myeloma cell lines (n = 16) compared to other cancer cell lines (n = 695) ( Figure 4A and Table 4). The CERES score of four genes (ISG20, NDC1, SF3B3, UMPS) was significantly lower in the myeloma cell lines (n = 20) compared to in the other cancer cell lines (n = 769) ( Figure 4B and Table 4). These essentiality scores were calculated using the RNAi [36] and CRISPR-Cas9 [35] methods, respectively. Thus, these analyses identified essential MM genes that had not previously been considered. . DEMETER2 and CERES scores were calculated from RNAi-based and CRISPR-Cas9 viability screens, respectively. A lower score means that a gene is more likely to be dependent in a given cell line. A score of 0 is equivalent to a gene that is not essential, whereas a score of -1 corresponds to the median of all common essential genes. Boxplots in orange represent multiple myeloma cell lines, and boxplots in grey correspond to other cancer cell lines. ** p-value < 0.01, *** pvalue < 0.001, **** p-value < 0.0001. . DEMETER2 and CERES scores were calculated from RNAi-based and CRISPR-Cas9 viability screens, respectively. A lower score means that a gene is more likely to be dependent in a given cell line. A score of 0 is equivalent to a gene that is not essential, whereas a score of -1 corresponds to the median of all common essential genes. Boxplots in orange represent multiple myeloma cell lines, and boxplots in grey correspond to other cancer cell lines. ** p-value < 0.01, *** p-value < 0.001, **** p-value < 0.0001.

Association between RNA-Seq-Based Risk Score and Response to Treatment
We analyzed the relationship between the RNA-seq-based risk score and HMCL response to drugs. The Spearman correlation was assessed between the RNA-seq-based risk score and IC50 for 22 different drugs [19][20][21][22]28,37]. Among them, the MM cell responses to four drugs were found to be significantly associated with the RNA-seq-based risk score ( Figure 5). A significant negative correlation was observed between the risk score and the response to EZH2 (r = −0.66; p-value = 0.037), Aurora kinase (r = −0.61; p-value = 0.028), MELK (r = −0.6; p-value = 0.04), and TOPK/PBK (r = −0.83; p-value = 0.0015) inhibitors in HMCLs ( Figure 5A). Moreover, as identified in the HMCLs, a significant correlation between the RNA risk score and the EZH2 inhibitor activity on the primary cells of MM patients was observed (r = −0.91; p-value = 0.0041) ( Figure 5B). A high RNA risk score value is associated with a higher EZH2 inhibitor efficacy in patient myeloma cells. These results highlight that MM patients associated with higher RNA-seq-based risk score mays benefit from these drugs and from EZH2 inhibitor in particular.

Discussion
Here, we defined an RNA-sequencing-based transcriptomic signature with prognostic value in two independent cohorts of patients with MM. Despite a significant accumulation of knowledge related to MM drug resistance, there is a need to routinely integrate these data into clinical decision making. However, several profiling methods have been developed to provide information related to molecular classification and risk prediction. Different groups have combined GEP analysis with cytogenetics to delineate 10 different molecular subgroups with distinct prognostic values and clinical features [4,38,39]. These prognostic classifications have been associated with the clinical data and incorporated

Discussion
Here, we defined an RNA-sequencing-based transcriptomic signature with prognostic value in two independent cohorts of patients with MM. Despite a significant accumulation of knowledge related to MM drug resistance, there is a need to routinely integrate these data into clinical decision making. However, several profiling methods have been developed to provide information related to molecular classification and risk prediction. Different groups have combined GEP analysis with cytogenetics to delineate 10 different molecular subgroups with distinct prognostic values and clinical features [4,38,39]. These prognostic classifications have been associated with the clinical data and incorporated into a consensus statement by the International Myeloma Working Group (IMWG) [40]. Several groups including IFM [7], UAMS [5], our group [8], and HOVON [6] have developed microarraybased GEP prognostic signatures. More recently, large-scale clinical studies have leveraged NGS including WES and have targeted NGS panels in MM [41][42][43][44][45][46]. Microarrays are being phased out, and our study sought to validate the use of RNAseq as an alternative. There is currently a need to integrate clinically useful biomarkers to predict treatment response in association with the growing palette of anti-MM therapies that are available.
In this study, we built a 267-genes risk score using RNA-seq data. The RNA-seqbased risk score demonstrated a prognostic value in the two independent cohorts of newly diagnosed MM patients. MM patients with high-risk UAMS-or RS-microarraydefined GEP signatures demonstrated significantly higher RNA-seq-based risk score values ( Figure 2B). Targeted sequencing could be used for clinical applications of the RNA-seq score where cost is important. Recently, Jang JS et al. reported a single cell RNA-seq of 597 cells derived from 15 MM patients. They developed a gene expression signature associated with MM progression that demonstrated significant prognostic value using the microarray GEP data of the APEX trial dataset [47]. We investigated the overlap between the two reported gene signatures and six common genes that were identified, including BST2, FAM214A, HBP1, ISG20, KLF6, and RAB30. The low degree of overlap could be explained by the fact that the signatures were defined using bulk RNA-seq or single-cell RNA-seq. The number of expressed genes detected from single cell RNA-seq is typically lower compared to bulk RNA-seq. The high-risk patients identified with our score were characterized by the genes involved in several major pathways implicated in MM pathophysiology, including interferon response, cell proliferation, hypoxia, IL-6 signaling pathway, stem cell genes, MYC, and epigenetic deregulation (Figure 2A). c-MYC is a key regulator in MM with deregulations related to translocations, gains and amplification, mutations in RAS genes, and MYC transcription or translation activation [48]. Hypoxia is a specific feature of MM with a significant increase of hypoxia-inducible factor-1 (HIF-1) in the bone marrow of MM tumor-bearing mice [49]. This suggests that the inhibition of HIF-1-mediated transcription may represent an interesting target in MM. Recently, we reported that chetomin, an inhibitor of HIF-1/p300 interaction, exhibits specific antitumor activity in human myeloma cell lines and in the primary MM cells from patients [50]. This approach could present therapeutic interest for high-risk patients identified with the RNAseq risk score. IL-6 is one of the major MM growth factors [51]. Blocking IL-6 signaling was thus developed into a therapeutic approach for MM. Even if the first clinical trials did not demonstrate clear benefits, the development of IL-6 antagonism is still ongoing with clinical trials [52]. The clinical and biological role of IFN is controversial in MM. Several groups have reported that IFN-α inhibits MM cell growth [53], whereas other groups have shown that it is an MM growth factor [54]. IFN-α was used in the treatment of MM, and this was stopped due to the absence of reproducible clinical efficacy [55,56]. More recently, we reported that DNA methyltransferase inhibitors induce the overexpression of IFN-regulated genes [20]. Since high-risk MM patients identified by RNA-seq are characterized by gene signature enrichment related to DNA methylation, the identified IFN response genes could be related to epigenetic deregulations. Interestingly, high-risk RNAseq score-defined MM patients are characterized by a significant enrichment in the genes related to stem cell genes. Genes that unrelated to the cell cycle and that are overexpressed in pluripotent, hematopoietic, and mesenchymal stem cells have been reported to be significantly overexpressed in MM in association with a poor outcome [57]. Furthermore, RNA-seq-defined high-risk patients present significant enrichment in terms of repressive epigenetic modifications, including the overexpression of the Polycomb repressive complex PRC1 and PRC2 target genes, DNMT target genes, and HDAC target genes ( Figure 2A). These data underline major epigenetic remodeling in high-risk MM patients. Interestingly, a significant overlap between the PRC2 and DNA methylation target genes has been reported in MM, suggesting an overlap between these repressive chromatin marks to inactivate important MM tumor suppressor genes [21]. Transcriptional programs mediated by DNA methylation and HDAC are also associated with poor outcome and key biological deregulations [19,20]. Furthermore, a combination of epidrugs has demonstrated anti-MM cell cytotoxicity in preclinical studies [21,58]. Of particular interest, using our large cohort of MM cell lines, we found that our RNA-seq-based risk score was significantly correlated with the response to the EZH2 inhibitor. Moreover, these data were validated using primary samples from MM patients. These data suggest that MM patients with a high-risk RNA-seq score may respond to EZH2i in combination with conventional treatment. Additionally, we also found that HMCLs with high RNA-seq score values presented sensitivity to three oncogene kinase inhibitors: MELKi, TOPK/PBKi, and Aurora kinase inhibitor. Aurora A kinase inhibitor in combination with bortezomib is currently in clinical development to treat relapsed or refractory MM patients (NCT01034553) [59]. MELK and PBK were identified as being of therapeutic interest in MM [22,60]. The RNA-seq score may be of interest for patient stratification in clinical trials. The validation of these results using the primary MM cells of patients will be of interest.
Moreover, among the genes composing the RNA-seq-based risk score and associated with bad prognosis, 11 genes were identified as significant essential MM genes, including SF3B3, FGFR4, TNFRSF10B, UMPS, and NPTXR. SF3B3 encodes a subunit of the splicing factor 3b protein complex. SF3B3 regulates EZH2 alternative splicing, and its expression is associated with poor outcome in renal cell carcinoma [70]. Moreover, the alternative splicing of EZH2 seems to have an important role in the tumorigenesis of human renal cancer. Since EZH2 is overexpressed in myeloma cells in association with a poor outcome [21,71], it would be interesting to explore the potential role of alternative EZH2 splicing regulated by SF3B3 in plasma cell tumorigenesis and disease progression. Furthermore, the potential role of SF3B3 in splicing regulation in MM remains unknown. The protein encoded by the FGFR4 (fibroblast growth factor receptor 4) gene is a tyrosine kinase and cell surface receptor for fibroblast growth factors (FGF). Acting as proangiogenic and mitogenic cytokines, FGF/FGFR protects myeloma cells for oxidative stress-induced apoptosis, leading to myeloma cell survival and progression [72]. In Rhabdomyosarcoma, FGFR4-specific single-domain antibodies demonstrated a good specificity and affinity for targeting FGFR4-expressing cells and for blocking the FGF19-FGFR4-MAPK signaling axis [73]. These FGFR4 antibodies could be of therapeutic interest for high-risk MM patients who have been identified with the RNA-seq-based risk score. TNFRSF10B, also called TRAIL-R2 or DR5, encodes a protein belonging to the TNF-receptor superfamily. The tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) induces apoptosis through the activation of DR5 [74]. It was already reported that the anti-DR5 monoclonal antibody, lexatumumab, induces myeloma cell death [75]. Moreover, because DR5 is transcriptionally regulated by p53, the efficiency of lexatumumab is increased by p53, inducing stress in myeloma cells [76]. UMPS (uridine monophosphate synthetase) is the last enzyme in the novo pyrimidine biosynthetic pathway. The inhibition of UMPS by 5-aminoimidazole-4-carboxamide-1-beta-riboside treatment decreases UMP levels and leads to pyrimidine starvation and myeloma cell death [77]. These data identify pyrimidine biosynthesis as a potential molecular target for future therapeutic targeting in MM. The neuronal pentraxin receptor is a type II transmembrane protein that functions as a trans-synaptic organizer and anchors neuronal pentraxin complexes to plasma membranes [78]. The role of NPTXR remains unclear in cancers. NPTXR is associated with cancer progression and shorter survival in gastric and colorectal cancers [79,80]. Silencing NPTXR using a monoclonal Ab against NPTXR inhibits gastric cancer cell proliferation and leads to cell apoptosis [79]. It could be interesting to investigate the function of NPTXR in myeloma cells.
Despite improvement in MM patient survival due to the clinical use of novel agents, the acquisition of drug resistance remains a major limitation of MM therapy. Indeed, the great majority of MM patients relapse and eventually become resistant to all treatments. We developed an RNA-seq-based risk score that allows the identification of high-risk MM patients that may benefit from EZH2, MELK, TOPK/PBK, and Aurora kinase inhibitors. The RNA-seq-based risk score may be used to implement precision medicine strategies in MM.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/jpm11100988/s1, Figure S1: Correlation between PCLI and RNA-seq based risk score, Table S1: List of the 267 prognostic genes constituting the RNA-seq based risk score.
Author Contributions: E.A. and V.V. designed and performed the research, analyzed the data, and wrote the paper; A.B., C.B. and A.K. performed the research and analyzed the data; N.R., L.V., C.H., G.R. and G.C. analyzed the data; O.E. analyzed the data and wrote the paper; J.M. designed and supervised the research and wrote the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.