key: cord-0020750-4sp3quua
authors: Dong, Chuanpeng; Cesarano, Annamaria; Bombaci, Giuseppe; Reiter, Jill L.; Yu, Christina Y.; Wang, Yue; Jiang, Zhaoyang; Zaid, Mohammad Abu; Huang, Kun; Lu, Xiongbin; Walker, Brian A.; Perna, Fabiana; Liu, Yunlong
title: Intron retention-induced neoantigen load correlates with unfavorable prognosis in multiple myeloma
date: 2021-09-09
journal: Oncogene
DOI: 10.1038/s41388-021-02005-y
sha: 245326124831fe15e21b19bed2956e5728bab8f5
doc_id: 20750
cord_uid: 4sp3quua

Neoantigen peptides arising from genetic alterations may serve as targets for personalized cancer vaccines and as positive predictors of response to immune checkpoint therapy. Mutations in genes regulating RNA splicing are common in hematological malignancies leading to dysregulated splicing and intron retention (IR). In this study, we investigated IR as a potential source of tumor neoantigens in multiple myeloma (MM) patients and the relationship of IR-induced neoantigens (IR-neoAg) with clinical outcomes. MM-specific IR events were identified in RNA-sequencing data from the Multiple Myeloma Research Foundation CoMMpass study after removing IR events that also occurred in normal plasma cells. We quantified the IR-neoAg load by assessing IR-induced novel peptides that were predicted to bind to major histocompatibility complex (MHC) molecules. We found that high IR-neoAg load was associated with poor overall survival in both newly diagnosed and relapsed MM patients. Further analyses revealed that poor outcome in MM patients with high IR-neoAg load was associated with high expression levels of T-cell co-inhibitory molecules and elevated interferon signaling activity. We also found that MM cells exhibiting high IR levels had lower MHC-II protein abundance and treatment of MM cells with a spliceosome inhibitor resulted in increased MHC-I protein abundance. Our findings suggest that IR-neoAg may represent a novel biomarker of MM patient clinical outcome and further that targeting RNA splicing may serve as a potential therapeutic strategy to prevent MM immune escape and promote response to checkpoint blockade.

Multiple Myeloma (MM) is characterized by the clonal expansion of malignant plasma cells in the bone marrow [1] . Recent therapeutic advances have extended overall survival, but most MM patients ultimately relapse [2] . Immune checkpoint blockade (ICB) therapy has revolutionized the treatment of many solid tumors by harnessing the immune system for effective anticancer treatment [3] . In these diseases, clinical response to ICB therapy is associated with the presence of tumor-specific antigenic peptides, or neoantigens [4] , a source of potential neoepitopes that can be loaded onto major histocompatibility complex (MHC) class I molecules to generate an antitumor immune response [5] .

Cytotoxic T-cells recognize tumor neoantigens as foreign and kill the presenting tumor cells, which initiates an antitumor immunological memory that hinders tumor recurrence. An important source of cancer neoantigens is somatic DNA mutations in the genome's coding regions [6] and the mutation-derived neoantigen load in several types of solid tumors corresponds with a better prognosis [7] [8] [9] [10] . However, MM has a relatively low mutation frequency. In contrast to solid tumors, mutation-derived neoantigen load in MM has been associated with unfavorable outcome [11, 12] .

Another potential source of tumor neoantigens is aberrant RNA splicing [13] [14] [15] [16] . Alternative splicing is a regulatory mechanism that generates multiple mRNA transcripts from a single gene and significantly expands proteome diversity [17] . Consequently, disruption of splicing mechanisms has a large impact on the transcriptome and is a significant driver of disease [18] . Intron retention (IR) occurs when the spliceosome fails to remove specific introns from pre-mRNA molecules, and they remain in the mature polyadenylated mRNA. In normal cells, IR functions to further reduce the levels of relatively low abundance transcripts that are not needed in specific cell types, such as the expression of developmentally regulated genes [19, 20] . This type of regulation has been termed transcriptome-tuning and is brought about through both nuclear RNA degradation and nonsense-mediated mRNA decay [21] .

IR occurs more frequently in nearly all cancer types compared with normal control tissues, even in the absence of DNA mutations in genes encoding proteins involved in splicing.

Additionally, in cancer cells, transcripts with IR are present at relatively high levels in cytoplasmic mRNA [22] . These transcripts are translated and degraded by the nonsense-mediated decay (NMD) pathway, a translation-coupled mechanism that eliminates mRNAs containing premature translationtermination codons [23] . Although most IR transcripts are subject to NMD-induced degradation, this process does not occur until after the pioneer round of translation, which can result in the production of neopeptides that bind to MHC molecules [24, 25] . Therefore, we hypothesized that IR-neoAgs in MM might impact immune response.

Herein, we used RNA-seq data from the MM Research Foundation (MMRF) CoMMpass study to identify IR events and predict IR-neoAgs. We found cells in bone marrow aspirates from MM patients exhibited high levels of IR events. However, consistent with the findings that high mutation-neoantigen load predicts unfavorable prognosis, high IR-neoAg load was correlated with shorter overall survival (OS) in MM. To investigate why high IR-neoAg load was not correlated with better MM patient survival, we performed gene set enrichment analysis on MM samples with high versus low IR-neoAg load. This analysis revealed that high IR-neoAg load was positively associated with the expression of T-cell inhibitory molecules, such as those involved in interferon (IFN) and tumor necrosis factor (TNF) alpha signaling activity. In addition, flow-cytometric analyses of four MM cell lines showed an inverse correlation between IR levels and MHC-II abundance, while treatment with a splicing inhibitor increased MHC-I protein abundance, especially in MM cells bearing high IR levels.

Genes involved in spliceosome activities are differentially expressed between MM and normal plasma cells To investigate whether the expression of genes involved in RNA splicing was altered in MM compared to normal plasma cells, we analyzed differentially expressed genes using RNA-seq data of plasma cells from five newly diagnosed MM patients (NDMM) and five healthy controls (GSE110486). These results showed that the spliceosome pathway was among the top upregulated pathways in MM (Fig. 1A) , where 67 out of 126 genes in the spliceosome pathway were upregulated significantly in MM samples. Gene set enrichment analysis also demonstrated that the spliceosome pathway was enriched in MM samples with a normalized enrichment score of 1.46 (p value < 0.001, FDR = 0.016, Fig. 1B) . We also found that the 230 upregulated differentially expressed genes identified in GSE110486 were also highly enriched in the NDMM samples from the MMRF cohort, as compared to the normal plasma cells (Fig. S1A ). In addition, the expression levels of 11 out of 12 serine and arginine-rich splicing factor (SRSF) protein genes, a conserved family of proteins involved in RNA splicing, were upregulated in MM (Fig. S1B ). Additional analysis of the MMRF data suggested that the increased expression of each of these 12 SRSF family genes was associated with decreased overall survival time (Fig. 1C) .

IR events are more common in MM compared to control plasma cells Accumulating studies provide strong evidence that IR is an important source of tumor neoantigens [15] . We sought to (Fig. 1D ).

Emerging evidence suggested that IR events in the cancer genome can be a source for immunogenic peptides [24] . Therefore, we investigated the potential for IR events to produce neoantigens in MM. To begin to address this question, we filtered the IR events that also occurred in normal plasma cells from the events identified in MMRF RNA-seq data. IR events occurring in normal plasma cells were removed because they were not expected to produce immunogenic peptides due to host immune tolerance. To identify the IR events in the healthy plasma cells, we analyzed RNA-seq data from the 13 plasma cell samples in GSE110486 and GSE114816. We detected a total of 9715 IR events that appeared in at least one healthy control sample; these IR events were eliminated from the list of events identified in the MM samples (Table S1 ). After filtering the normal IR events, the average number of MM-specific IR events per sample was 1009 and ranged from 21 to 4138 (Table S1 ). Interestingly, gene ontology analysis of 450 genes harboring MM-specific IR events that occurred in more than half of the NDMM samples showed that these genes were enriched in pathways involving RNA processing and RNA transport ( Fig. 2A) .

To computationally predict IR-neoAgs, we first determined the HLA-I genotype of each MMRF patient using the RNA-seq data. A total of 178 unique HLA-A/B/C alleles were identified from 767 individual patients of the MMRF cohort (HLA alleles and their frequencies are listed in Supplementary Table S2 ). Next, the retained intron sequences were translated into protein sequences, which were then segmented into 8-11 amino acid peptides, where at least one amino acid was translated from the intronic region. Any peptide that could also be generated from normal proteins was further removed. The remaining IR-derived neopeptides were then evaluated for their predicted binding affinity with the set of patient-specific HLA alleles using NetMHCpan (v4.1). Peptides with a NetMHCpan predicted rank score less than two (the default cutoff from NetMHCpan) were selected as IR-neoAgs. IR-neoAgs were called for 893 RNA-seq samples from the MMRF cohort (including both newly diagnosed and relapsed samples). Not surprisingly, the number of IR events and the IR-neoAg load were highly correlated (Spearman correlation rho = 0.97, p value < 0.0001, Fig. 2B ). We further evaluated whether any HLA allele presented more IR-neoAgs than other alleles at the population level ( was detected in 10.8% of samples. The most common allele, HLA-A02:01 which was detected in 19.4% of the samples, presenting 31,161 IR-derived peptides. Notably, 24,680 out of 479,685 of the IR-neoAgs were shared across more than 5% of the multiple MM samples. This observation would suggest that there might be potential for developing cancer vaccines in the future based on IR-neoAgs. We also found that 20 neoantigens occurring in more than 80% of NDMM samples were preferentially presented by HLA-C alleles (Fig. 2D ), suggesting neoantigens presented by HLA-C alleles could be prioritized for cancer vaccine development.

IR-neoAg load correlates with unfavorable clinical outcome We next asked whether IR-neoAg load was associated with overall survival (OS) in the MMRF cohort. Kaplan-Meier survival analysis revealed that NDMM patients with higher than the median IR-neoAg load had significantly shorter OS (log-rank test, P = 0.027, Fig. 3A ). When considering the expression levels of IR-neoAgs, we observed an even more substantial prognostic effect, with a p value reaching 0.006 (Fig. 3B) . Similarly, higher than the median IR-neoAg load predicted shorter OS for MM patients at the time of first relapse (log-rank test, P = 0.002, Fig. 3C , n = 60). Notably, relapsed MM samples with lower IR-neoAg load had a higher 2-year OS rate compared to patients with higher IR-neoAg load (OS 0.85 vs. 0.57).

To determine whether IR-neoAg load was associated with clinical features of MM, we asked whether IR-neoAg load correlated with the International Staging System (ISS) [26] , which is a reproducible predictor of MM outcome. We did not find that IR-neoAg load was associated with the ISS disease stage in NDMM from the MMRF cohort (one-way ANOVA P = 0.724, Fig. S2A ). To determine whether the addition of IR-neoAg load to the ISS stage improved the prediction of OS, we performed survival analysis on patients stratified by disease stage and IR-neoAg load. This analysis showed that stage II MM patients with higher than the median IR-neoAg load had significantly shorter OS than stage II patients with low IR-neoAg load (log-rank test, P = 0.002, Fig. 3D) . A similar trend was observed with stage I patients, although the association did not reach statistical significance (P = 0.14). IR-neoAg load had no apparent prognostic value for OS in stage III MM patients (P = 0.486). In addition to ISS stages, chromosomal hyperdiploidy (HRD) is widely used in defining genetic subtypes of MM patients, and HRD-myeloma is associated with better survival compared to nonhyperdiploid (nHRD) MM [27] . Although we observed that higher than the median IR-neoAg load was apparently associated with shorter OS in both HRD and nHRD MM patients, these associations did not reach statistical significance (Fig. S2B) .

We further examined whether the prognostic performance of IR-neoAg load was independent of other clinical factors. We used multivariate Cox analysis to test the performance of IR-neoAg load after adjusting for other clinical factors, including age, sex, P53 status, ISS stage, as well as the revised ISS stage after adjusting for lactate dehydrogenase (LDH) level, chromosomal aberrations, and other factors. In the multivariate analysis, the hazard ratio of high versus low IR-neoAg load for OS in NDMM was 1.491 (p value = 0.027; 95% CI 1.056 to 2.492) ( Table 1 ), indicating that the IR-neoAg load offers prognostic power that is independent of other clinical factors.

Higher T cell inhibitory signals associate with IR-neoAg and poor prognosis in MM Our observation that higher IR-neoAg load was associated with shorter OS is consistent with previous reports of mutation-derived neoantigen load in MM [11, 12] . However, this finding is the reverse of previously reported observations that high mutationderived and IR-neoAg loads are associated with longer OS in patients with solid tumors, including melanoma [7] , lung cancer [8] , breast cancer [10] , and pancreatic cancer [6] . In addition, there is increasing evidence that T cells present in the MM microenvironment show an exhausted and suppressed phenotype [28] . This would suggest that additional changes in MM plasma cells may affect the anti-MM immune response.

To test this hypothesis, we conducted differential expression and gene set enrichment analysis on RNA-seq data from the MMRF cohort, comparing samples from NDMM patients with either higher or lower than the median IR-neoAg load. Notably, we observed a significant enrichment of the pathways related to T-cell suppression. We found that IFN gamma signaling and TNFα signaling via NF-κB pathways were upregulated in patients with high IR-neoAg load (Fig. 4A) . These pathways are involved in the recruitment of T-regulatory (Treg) cells that control cytotoxic T-cell killing. In addition, the B-cell receptor (BCR) signaling pathway was significantly enriched in patients with high IR-neoAg loads. Previous studies demonstrated that sustained activation of BCR signaling plays critical roles in B-cell malignancies [29] . This result suggests that molecular features in cells with higher IR-neoAg load might contribute to T-cell and B-cell dysfunction in MM.

Based on this finding, we postulated that increased expression of T-cell co-inhibitory molecules in MM cells exhibiting high IR-neoAg load might be a partial explanation for the reduced antitumor immunity and thereby facilitate cancer immune evasion [30] . These co-inhibitory molecules function as brakes to inhibit T-cell activation. Higher expression levels of co-inhibitory ligands on the cancer cell surface can negatively impact T-cell function. To begin to address this question, we first analyzed the expression levels of 20 co-inhibitory genes identified by Dufva and colleagues [31] , which include genes for eight B7 ligands, six enzymes impacting T-cell activity, and six other genes from the butyrophilins and CD226 family (Table S3) . We found that these coinhibitory genes tend to have higher expression levels in NDMM samples with higher IR-neoAg load (Fig. 4B) . We found that CD274 (PD-L1) expression was 1.3-fold higher in patients with high IR-neoAg load (adjusted p value < 0.0001), suggesting there could be a stronger immune suppression in patients with higher IR-neoAg load. Next, we analyzed co-inhibitory gene expression from 29 MM cell lines compared to other cancer cell lines in the Cancer Cell Line Encyclopedia (CCLE). Surprisingly, we found that the average expression level of PD-L1 in the MM cell lines was lower than most other types of cancer cell lines (Fig. S3) , which might partially explain why anti-PD1 therapy has had a limited response rate in MM. Other B7 co-inhibitory ligands, such as CD86, CD80, and HHLA2, showed high expression levels in myeloma cell lines relative to the other cancer cell lines, implying that these B7 B T-cell signaling co-inhibitory genes were upregulated in patients with high IR-neoAg load. The genes with adjusted p value < 0.05 are labeled in purple; the blue bar indicates downregulated genes in NDMM patient samples with high IR-neoAg load and the red bar indicates upregulated genes. C Kaplan-Meier survival curves showing overall survival in MMRF cohort patients with high (>= 75%) and low (< 75%) expression of HHLA2 and high (above the median) or low (below the median) IR-neoAg load.

ligands might serve as potential targets for immune checkpoint therapy. Kaplan-Meier survival analysis revealed that the patients with higher HHLA2 and IR-neoAg load had the worst outcome (Fig. 4C) , which provides further support that HHLA2 may be a druggable target for treating MM in the future [32] .

MHC molecules encoded by the HLA-I and HLA-II genes are essential components in IR-neoAg presentation on the cell surface. Therefore, we investigated the relationship between IR events (IR levels) and MHC protein abundance in four MM cell lines, namely JJN3, U266, KMS11, and AMO1 cells. These cell lines were selected because KMS11 and U266 had the highest levels of IR, while JJN3 and AMO1 had the lowest levels of IR based on RNA-seq data from the CCLE consortium (Fig. 5A ). We measured MHC-I and MHC-II cell surface abundance in these MM cells by flow cytometry before and after treatment with the splicing inhibitor pladienolide-B for 96 h. As demonstrated in Fig. 5B , the basal cell surface level of MHC-II was lower in the KMS11 and U266 cell lines bearing higher IR levels, compared to the JJN3 and AMO1 cell lines bearing lower IR levels. Low MHC-II abundance in MM cells with high IR levels is consistent with our observation in MM patients where higher IR levels and low HLA-II gene mRNA expression was associated with worse clinical outcomes (data not shown).

Next, we investigated whether splicing activities, as measured by the mRNA expression levels of genes encoding key splicing factors and regulators, correlated with the mRNA expression levels of the HLA genes encoding MHC-I and MHC-II complexes in the MMRF RNA-seq data. We observed a negative correlation between the expression levels of MHC-I genes and spliceosome pathway activities, as measured by ssGSEA enrichment scores (Fig. 5C , Spearman correlation rho = −0.245, p < 0.001). No correlation was observed between the expression of genes encoding MHC-II molecules and spliceosome pathway activities (Fig. S4A) .

To investigate whether low MHC-I expression might be a result of increased splicing activity, we treated MM cell lines with the splicing inhibitor pladienolide-B, which targets SF3B1, a gene encoding subunit 1 of the splicing factor 3b protein complex, and measured MHC-I cell surface expression by flow cytometry. We found that MHC-I expression levels were significantly increased in three of the 4 MM cell lines, including both cell lines with higher IR levels (Fig. 5D ). As shown in Fig. 5E , MHC-I cell surface abundance in KMS11 MM cells exhibiting high IR increased following pladienolide-B treatment. This finding strongly suggests that modulation of splicing activity may regulate the abundance of MHC-I class proteins along with the antigen presentation potential in MM cells. Consistent with the lack of correlation between HLA-II gene expression and spliceosome pathway activity (Fig. S4A) , no significant changes in MHC-II protein abundance were observed in the MM cell lines after splicing inhibition (Fig. S4B ).

In this study, we demonstrate that intron retention is an important source of neoantigens in multiple myeloma, which impacts patient clinical outcome. We showed that newly diagnosed MM samples exhibited more intron retention events than normal plasma cells and that higher IR-neoAg load was significantly associated with unfavorable survival in both newly diagnosed and relapsed MM. Our findings indicate that bioinformatic predictions of immune recognition of neoantigens arising from genomic or transcriptomic alterations in MM might not be useful in selecting patients for immune checkpoint therapy. Further, our analyses revealed that poor outcome in MM patients with high IR-neoAg load is associated with higher expression levels of checkpoint genes and elevated IFN signaling activity, which implies strong T-cell suppression. Therefore, our results suggest a potential mechanism for MM cell immune evasion despite having an increased neoantigen load compared to normal plasma cells. Whereas high neoantigen load generally predicts favorable survival and higher likelihood of response to checkpoint blockade in many solid tumors such as breast cancer [10] , lung cancer [8] , glioblastomas [9] and melanoma [7] , we found that a high neoantigen load in MM patients was associated with poor prognosis [11] . In addition, the immune context of the bone marrow microenvironment is more complex compared with solid tumors, where cytokines and immune cell components in the bone marrow provide a unique seedbed for myeloma cell growth [33] . Therefore, the underlying mechanisms that allow for MM cell immune escape are apparently different from other tumors.

In addition to somatic DNA mutations, RNA alternative splicing, including intron retention, has been reported to be a novel source of neoantigens [24] . Numerous studies have reported that the splicing machinery is dysregulated in multiple cancer types, including bladder cancer [34] , breast cancer [35] , melanoma [36] , prostate cancer [37] and hematological cancers [38, 39] . In addition, intron retention events have been observed frequently in prostate cancer [40] and pancreatic cancer [41] . Yang et al. reported that blood cells have a high level of splicing diversity compared to other tissues, next to testis, brain, and muscleskeletal tissue, in the GTEx transcriptional data [42] . IR events represent a large proportion of alternative splicing events in blood tissue. These findings prompted us to investigate whether IR-neoAg could contribute to immune responses in MM, in particular to antigen presentation and T-cell mediate responses. We found that higher IR-neoAg load was significantly associated with shorter survival time, both in newly diagnosed and relapsed MM. This finding was further strengthened when the expression levels of IR-neoAg were considered (p-value reached 0.006).

Over the past decade, immune checkpoint blockade (ICB) therapy has revolutionized cancer therapy in several tumor types [43] . However, response to the immune checkpoint inhibitor pembrolizumab (anti-PD1) has been limited in MM [44] . Clinical response to ICB has been closely linked with the abundance of tumor-specific neoantigens, the presence of cytotoxic T-cell infiltration, and distinct tumor microenvironment profiles. Previous reports have demonstrated an increase in the mutationderived neoantigen load in MM and have also confirmed a neoantigen T-cell response in relapsed patients with MM [12] . These results implied that a T-cell mediated immune response might be suppressed or impaired in MM. Zelle-Rieser et al. reported that CD8 + T-cells expressed several molecules associated with T-cell exhaustion (PD-1, CTLA-4, CD160) as well as the T-cell senescence marker CD57 at the MM tumor site [28, 45] . Our results showed that higher IR-neoAg load was positively correlated with higher expression levels of T-cell inhibitory molecules and genes belonging to the Tregs activating pathway. Dufva et al. reported that decreased HLA-II gene expression might be a potential immune evasion mechanism in hematological cancers [31] . We also found that gene expression levels of HLA-II genes were significantly lower in newly diagnosed MM compared with healthy control cells.

Despite these RNA-seq-based observations, direct evidence of IR-neoAg presentation on MM cells using immune-peptidomics technology could strengthen our conclusion. We hypothesize that standard MM treatment options do not generate an effective immune response that leverages the neoantigen immunotherapeutic potential. Indeed, we observed lower levels of MHC-II activity in MM cell lines with higher intron retention. In addition, we observed that MHC-I activity appeared to be inhibited in cells with elevated expression levels of splicing factors, a hallmark of MM, and that inhibition of spliceosome activity resulted in increased MHC-I activity. Collectively, these two mechanisms may partially explain why higher IR-induced neoantigen load in MM samples was not associated with better prognostic outcome. This result also suggests that splicing inhibitors could possibly boost the efficacy of immune checkpoint blockade therapy in MM by activating MHC-I presentation [46, 47] . Further analysis with integrated multi-omics data from different aspects of the immune landscape is needed to further understand the potential determinants of responsiveness to cancer immunotherapies in MM.

In conclusion, while neoantigen load has been associated with favorable survival in many solid cancers, our study strongly suggests that IR-neoAg load may serve as a clinically relevant risk factor that negatively impacts myeloma patient survival. Our analysis provides evidence that MM cells bearing high levels of IR-neoAgs also present T-cell inhibitory gene signatures, which may offset the neoantigen load in eliciting a cytotoxic T cell response. Moreover, we found that aberrant RNA splicing may also regulate MHC abundance and thus, contribute to MM immune escape. Our findings highlight the need to integrate multi-omics data to uncover the immune context and understand the factors that determine responsiveness of MM to immunotherapies. Also, this works suggests that targeting splicing may represent an additional therapeutic strategy to promote anti-MM immune response.

The raw data from the MMRF study was obtained through an authorized data access request for dbGaP study accession: phs000748.v7. p4. RNA-seq data from 893 samples, including both newly diagnosed and relapsed subjects, were downloaded and converted to fastq format using SRA-tools (v2.10.0). Curated survival and clinical data were downloaded from the UCSC Xena cancer browser (http://xena.ucsc.edu) [48] . The revised International Staging System (R-ISS) was calculated as defined by the International Myeloma Working Group [49] , by considering the presence of del(17p), t(4;14), and t(14;16) and information on serum β2-microglobulin, albumin, and lactate dehydrogenase levels. B2M mutations and the status of TP53 in baseline samples were obtained from Dr. Brian Walker, as described previously [50] .

Two other RNA-seq studies with normal plasma cells were retrieved from the Gene Expression Omnibus (GEO). Data of bone marrow-derived plasma cells from five healthy individuals and five newly diagnosed MM patients were obtained from GSE110486 [51] . Data of plasma cells from bone marrow or tonsil of another eight normal subjects were acquired from GSE114816 [52] .

For all RNA-seq data, an initial sequence-level quality assessment was performed using FastQC (v0.11.5). The alignment-free quantification tool Salmon (v1.2.1) [53] was used to quantify the expression of gene transcripts from RNA-seq data using the reference transcriptome built from Gencode (GRCh38, v32) gtf annotation as the index. The gene-level transcript abundance was calculated using the tximport package in R.

The normalized gene expression data of 887 cancer cell lines (dated 2018.09.29) and their annotations (dated 2018.12.26) were downloaded from the Cancer Cell Line Encyclopedia (CCLE) [54] data portal (https:// portals.broadinstitute.org/ccle/data).

To quantify the IR events for MM samples, RNA-seq reads were aligned to the GRCh38 reference genome using STAR (v2.7.2) [55] . Uniquely mapped RNA-seq reads were used to quantify the expression levels of retained introns using HTseq [56] package. Additional criteria were applied to filter the identified IR events: (1) read counts for both the intron region and its flanking exon regions were > 10; (2) read coverage of the intron was comparable to its flanking exons, such that the transcripts per million (TPM) ratio of introns to flanking exons was > 0.05 and < 0.5. MM-specific IR events were further selected by removing the events that were observed in normal plasma cells using the same filtering criteria.

We used arcasHLA (v1.1) to infer HLA class-I genotypes from the RNA-seq data [57] . Sequences from the retained introns were translated into peptides by extending the open reading frame of the upstream exon using the standard codon table. The translated peptides were segmented into 8 to 11 amino acid lengths that contained at least one intron-encoded amino acid. For each patient, NetMHCpan4.1 was used to estimate the binding affinity of the IR-derived neopeptides with the patient's HLA alleles [58] . A binding affinity rank score less than 2 (default parameter of NetMHCpan4.1) was regarded as a neoantigen candidate. Expression levels for each neoantigen were determined by the abundance of IR events that generated the specific neoantigen.

Differentially expressed genes in plasma cells between MM and healthy bone marrow samples were identified using the limma [59] package in R. A total of 1790 gene sets from MsigDB version 7.2, including KEGG, REACTOME, and HALLMARK gene sets, were used for enrichment analysis. Fisher's exact test was used to test for pathway enrichment significance, and the p-value was adjusted for multiple hypothesis correction using the Bonferroni method [60] . ClusterProfiler was used to visualize the pathway enrichment results [61] . Single-sample gene set enrichment analysis (ssGSEA) was used to assess the pathway activity in each individual using GSVA [62] package in R, using default parameters.

Cell culture of MM cells and spliceosome inhibition KMS11, U266, JJN3, AMO1 MM cell lines were kindly provided by Dr. David Roodman and cultured at 37°C in a humidified atmosphere containing 5% CO 2 and maintained in RPMI media supplemented with 10% fetal bovine serum (FBS). Cells were tested for mycoplasma infection monthly as a regular lab routine. Pladienolide-B (Cayman Chemical Company, Ann Arbor, MI; cat# 16538) was dissolved in dimethyl sulfoxide and used at the following concentrations: 0, 0.1, 1, 5, 10, and 100 nM. Flow-cytometric analyses were performed at 96 h post-treatment, two independent biological replicates were analyzed for each treatment.

The following flow-cytometry antibodies were used: HLA-DR, DQ, DP-APC (Biolegend, San Diego, CA cat# 361714), Isotype control-APC (Biolegend cat# 400222), MHC class I-PE (LSBio, Seattle, WA cat# LS-C751033-0.1). Isotype control-PE (Abcam, Waltham, MA cat# ab91357). Flow cytometric data were acquired using the LSR II flow cytometer (BD Biosciences, San Jose, CA) and analyzed with FlowJo software.

Survival analysis and Cox-proportional hazard comparison were performed using the R package Survival with log-rank test and hazard ratio statistical tests [63] . Significant differences in the value of the two given groups were assessed using the Mann-Whitney-Wilcoxon test [64] . Statistical analyses were performed in R (v4.0.2).

Global burden of multiple myeloma: a systematic analysis for the global burden of disease study

Concurrent lipidomics and proteomics on malignant plasma cells from multiple myeloma patients: Probing the lipid metabolome

Fundamental mechanisms of immune checkpoint blockade therapy

Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade

Neoantigens in cancer immunotherapy

Exploiting the neoantigen landscape for immunotherapy of pancreatic ductal adenocarcinoma

Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma

Predictive biomarkers of response to PD-1/PD-L1 immune checkpoint inhibitors in non-small cell lung cancer

The combination of neoantigen quality and T lymphocyte infiltrates identifies glioblastomas with the longest survival

HLA class-I and class-II restricted neoantigen loads predict overall survival in breast cancer

High somatic mutation and neoantigen burden are correlated with decreased progression-free survival in multiple myeloma

Mutation-derived Neoantigen-specific T-cell Responses in Multiple Myeloma

Cancer-specific splicing changes and the potential for splicing-derived neoantigens

RNA transcription and splicing errors as a source of cancer frameshift neoantigens for vaccines

Alternative mRNA splicing in cancer immunotherapy

IRIS: Big datainformed discovery of cancer immunotherapy targets arising from pre-mRNA alternative splicing

Expansion of the eukaryotic proteome by alternative splicing

Defective control of pre-messenger RNA splicing in human disease

Coordinated regulation of neuronal mRNA steady-state levels through developmentally controlled intron retention

Orchestrated intron retention regulates normal granulocyte differentiation

Widespread intron retention in mammals functionally tunes transcriptomes

Intron retention as a component of regulated gene expression programs

Nonsense-mediated mRNA decay (NMD) mechanisms

Intron retention is a source of neoepitopes in cancer

Major source of antigenic peptides for the MHC class I pathway is produced during the pioneer round of mRNA translation

International staging system for multiple myeloma

Understanding the role of hyperdiploidy in myeloma prognosis: which trisomies really matter?

T cells in multiple myeloma display features of exhaustion and senescence at the tumor site

Targeting the B-cell receptor signaling pathway in B lymphoid malignancies

Co-inhibitory molecules of the B7-CD28 family in the control of T-cell immunity

Immunogenomic Landscape of Hematological Malignancies

HHLA2 is a member of the B7 family and inhibits human CD4 and CD8 T-cell function

The Immune Microenvironment in Multiple Myeloma: Friend or Foe? Cancers (Basel)

Polypyrimidine tract binding protein 1 promotes lymphatic metastasis and proliferation of bladder cancer via alternative splicing of MEIS2 and PKM

The splicing factor SRSF1 regulates apoptosis and proliferation to promote mammary epithelial cell transformation

Splicing factor SRSF6 promotes hyperplasia of sensitized skin

Dysregulation of the splicing machinery is directly associated to aggressiveness of prostate cancer. EBioMedicine

Aberrant expression of splicing factors in newly diagnosed acute myeloid leukemia

SnapShot: splicing alterations in cancer

Intron retention is a hallmark and spliceosome represents a therapeutic vulnerability in aggressive prostate cancer

Intron retention is a robust marker of intertumoral heterogeneity in pancreatic ductal adenocarcinoma

A systematic survey of human tissue-specific gene expression and splicing reveals new opportunities for therapeutic target identification and evaluation. bioRxiv

Checkpoint blockade in Hodgkin and non-Hodgkin lymphoma

Lessons learned from checkpoint blockade targeting PD-1 in multiple myeloma

Blocking IFNAR1 inhibits multiple myeloma-driven Treg expansion and immunosuppression

Therapeutically increasing MHC-I expression potentiates immune checkpoint blockade

Pharmacologic modulation of RNA splicing enhances anti-tumor immunity

Visualizing and interpreting cancer genomics data via the Xena platform

Revised International Staging System for Multiple Myeloma: A Report From International Myeloma Working Group

A highrisk, Double-Hit, group of newly diagnosed myeloma identified by genomic analysis

Aberrant hyperediting of the myeloma transcriptome by ADAR1 confers oncogenicity and is a marker of poor prognosis

Long non-coding RNAs discriminate the stages and gene regulatory states of human humoral immune response

Salmon provides fast and bias-aware quantification of transcript expression

The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity

STAR: ultrafast universal RNA-seq aligner

HTSeq-a Python framework to work with highthroughput sequencing data

arcasHLA: highresolution HLA typing from RNAseq

NetMHCpan-4.1 and NetMHCIIpan-4.0: improved predictions of MHC antigen presentation by concurrent motif deconvolution and integration of MS MHC eluted ligand data

limma powers differential expression analyses for RNA-sequencing and microarray studies

STATISTICAL QUESTION Multiple significance tests: the Bonferroni correction

clusterProfiler: an R package for comparing biological themes among gene clusters

GSVA: gene set variation analysis for microarray and RNA-seq data

Cox Analysis of Survival Data with Non-Proportional Hazard Functions

A class of Mann-Whitney-Wilcoxon type statistics

The computational algorithm and source code allowing for reproduction of the intron retention-induced neoantigen quantification in this manuscript are available at https://github.com/cpdong/IntronNeoantigen.

The authors declare no competing interests.

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41388-021-02005-y.Correspondence and requests for materials should be addressed to Fabiana Perna or Yunlong Liu.Reprints and permission information is available at http://www.nature.com/ reprintsPublisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.