Alternative Splicing in Multiple Myeloma is Driven by Translocations and Biallelic Events in DIS3.
Previously we determined that mutations in SF3B1 resulted in increased AS (Fig. 1A) in MM patient samples(12). However, only 1.6% of patients have this mutation, indicating that AS may not be as important in MM than in other hematological malignancies, or that there are other mechanisms specific to MM that are at play. To determine the full extent of alternative splicing in MM, we performed UMAP clustering of 598 NDMM samples using the top 1% most-variate events, Fig. 1B-D and Supplementary Fig. 1. A clear segregation could be observed between samples with or without IgH translocations (Fig. 1B). AS events also clustered within the translocation groups with clear clusters for the t(4;14), t(11;14), and t(MAF) (consisting of t(14;16), t(14;20), and t(8;14)) groups. This is in line with other genomic studies where the primary translocation events drive clustering of gene expression and DNA methylation data(30) (31, 32). Using our previous annotations for this dataset we also noticed that biallelic DIS3 samples were separate from mono-allelic and non-mutated samples (Supplementary Fig. 1B), and were associated with the t(4;14) and t(14;16) groups. In contrast, despite being the most frequently mutated gene KRAS mutated samples were evenly distributed throughout the whole population, Fig. 1D, indicating no role of this common mutation in driving AS. No significant difference (Two-sided Mann-Whitney U test) between the number of events with positive and negative dPSI were observed in any comparison (Supplementary Fig. 2).
To determine if other genomic factors could be involved in generating AS in MM we performed 16 comparisons that cover all major cytogenetic subgroups (IgH translocations and hyperdiploidy) and the most common mutations found in MM (BRAF, KRAS, NRAS, DIS3, TENT5C, and TP53), SF3B1 mutation and SF3B1 mutation hotspots, as well as those with biallelic abnormalities including DIS3, TENT5C, and TP53, plus one combination event where there is an association of DIS3 mutations in the t(4;14) (t(4;14) + DIS3 mutation). Sample size and comparison groups are described in Supplementary Table 1.
Differentially spliced events were identified in each comparison and further compared against a set of normal plasma cell RNA-seq data (N = 5) to remove any cell-type specific events. Significant AS events in each comparison identified 4 045 AS events (|dPSI|>10%, TPM > 1, P < 0.05), Fig. 1E and Supplementary Tables 2 and 3. Samples with SF3B1 mutation hotspots contained the greatest number of AS events (N = 862), and samples with any SF3B1 mutation had approximately half as many. Samples with primary IGH translocations had an equivalent number of AS events to those with SF3B1 mutations, with t(14;16) having the most (N = 587) followed by t(11;14) (N = 366), and t(4;14) (N = 256). The proportion of AS event types was similar among all comparisons (Fig. 1F), with alterative first (AF) exons being the dominant event type, followed by skipped exons (SE), and alternative last exons (AL). The least frequent event types were mutually exclusive (MX) exons. The different AS events for the main translocation groups are shown along with representative examples from each group, Fig. 2A-G.
DIS3 at 13q and TENT5C at 1p are frequently mutated in MM and play a part in RNA processing. Both genes are on chromosomal regions that are frequently deleted in MM patient samples and so we also looked at the effect of monoallelic and biallelic alterations at the loci. We identified increased AS in the biallelic state compared to the monoallelic state for both DIS3 (P = 0.02; one-sided Mann-Whitney U test) and TENT5C (P = 0.01; one-sided Mann-Whitney U test), Fig. 1B. The same was true for TP53, albeit at lower levels. As DIS3 mutations were also enriched in the t(4;14) subgroup, we compared t(4;14) samples with and without DIS3 mutations to determine the effect on AS. This analysis indicated an increased number of AS events as a result of DIS3 mutations on a t(4;14) background. As del(13q) is common in MM (~ 50% of patients) and is enriched in t(4;14) MM (95% with del(13q)), we compared samples with or without del(13q) which resulted in very few (N = 32) AS events, indicating that del(13q) is not a confounder in this comparison. As expected, the frequently mutated genes KRAS, NRAS, and BRAF, which encode RAS pathway proteins, had negligible effects on AS (Fig. 1E).
Alternative Splicing Events Affect Key Genes and Pathways in Multiple Myeloma
We further examined the AS events among all comparisons (Fig. 2). In the t(4;14) comparison, 256 AS events were identified (Fig. 2A). These included YIPF1, a gene that encodes a protein and plays a role in the Golgi reassembly and glycan synthesis. An alternative transcript, YIPF1-203, was identified with an extra exon near the 3’ end compared to YIPF1-201 (canonical transcript) and YIPF1-202 with the exon skipped (Fig. 2B). As a result, YIPF1-203 was predicted to undergo nonsense mediated decay (Supplementary Fig. 3B), while YIPF1-201 and YIPF1-202 were protein coding transcripts(33). Although expression of YIPF1 at the gene level was not different, the YIPF1-203 transcript had significantly higher expression in the t(4;14) while the YIPF1-201 and YIPF1-202 transcripts had significantly higher expression in the non-translocated samples (Fig. 2C). Therefore, even though the overall gene expression level of YIPF1 were similar, this gene may lose its function in the t(4;14) due to the high expression of the nonsense mediated decay YIPF1-203 transcript.
In the t(14;16) comparison, SLAMF7, which encodes a cell surface immunotherapeutic antigen that can be targeted in MM(34) and immune cells(35), underwent a mutually exclusive exon event generating transcripts SLAMF7-204 and SLAMF7-201 (Fig. 2D-E). SLAMF7-204 was expressed at significantly lower levels in t(14;16) samples (log2FC=-0.52, P = 3x10− 8) while SLAMF7-201 showed no significant expression difference between t(14;16) and non-translocated samples (Fig. 2F). Compared to the canonical SLAMF7-203 transcript, SLAMF7-201 contained different amino acid sequences at the intracellular domain regions (258–335) (Supplementary Fig. 3A). However, its extracellular and transmembrane domain were intact and thus predicted to encode a transmembrane protein (Supplementary Fig. 3A)(36). SLAMF7-204 encoded an isoform without a transmembrane domain, and therefore was predicted to encode a secreted protein (Supplementary Fig. 3A)(36). Cell surface bound SLAMF7 is targeted by Elotuzumab therapy while secreted SLAMF7 is used in diagnosis, monitoring and assessment of response of such therapy(37). Significantly less SLAMF7-204 would result in less secreted SLAMF7, and as a result, the response of Elotuzumab therapy may not be correctly assessed for t(14;16) patients.
Alternative Splicing Events are Highly Unique Among Comparisons
When comparing the AS events identified in the different genomic comparisons, we found that most of them are unique with little overlap, Fig. 2G and Supplementary Fig. 4, except similarities within the translocation comparisons, and within the DIS3 comparisons. In addition, there were a higher number of consistent events between the t(4;14) or t(14;16) and biallelic DIS3 comparisons, perhaps reflecting the enrichment of this abnormality within these translocation groups.
Within the DIS3-related comparisons, mutation alone identified the fewest AS events (Fig. 2H), and other alterations built on top of these. For example, 88% of AS events in the DIS3 mutation comparison were consistently identified in the biallelic DIS3 comparison and 56% in the t(4;14) + DIS3 comparison. In addition, the biallelic DIS3 or t(4;14) + DIS3 comparisons identified significantly more additional unique AS events. This may indicate that complete loss of wild type DIS3 has a profound effect on RNA splicing and that the co-occurrence of DIS3 mutations in the t(4;14) may have a biological role in AS, perhaps through epigenetic coordination with NSD2.
Event-level heterogeneity also led to functional heterogeneity. Pathway analysis of AS genes indicated functional differences within the subgroup analysis (Fig. 2I and Supplementary Fig. 5). Common pathways can be found among translocation groups including metabolic process, immune process, and cell mobility, even though most of the categories are upregulated in t(14;16) and t(4;14) but inconsistently dysregulated in t(11;14). t(14;16) contributes the most number of dysregulated pathways with cell differentiation and metabolic process being the most upregulated.
DIS3 abnormalities had a profound negative effect on many pathways with down regulation of genes. In particular, RNA processing pathways were down regulated across all DIS3 categories. The combination of t(4;14) + DIS3 dysregulated the most pathways, all of which were downregulated. Biallelic DIS3 contributed the second most, which included immune and metabolic processes, RNA processing and protein processing.
Alternative Splicing Frequency is Determined by the Activity of the Non-Homologous End Joining Pathway
We speculated that spliceosome expression regulates the splicing frequency of transcripts. For each genetic comparison, we calculated the log2(Fold-Change) of all spliceosome genes and performed unsupervised clustering. Three heat shock proteins, HSPA1A, HSPA1B and HSPA6 had the highest correlations of their expression with splicing frequency (\(\rho >0.6\), Spearman correlation, Fig. 3A). All three genes encode proteins that belong to the Hsp70 protein family that contains various protein homologs, including Hsp73, Hsp70B’ and Hsp72(38). Hsp73 is the constitutive form and is encoded by the spliceosome gene HSPA8(38). Hsp72 (encoded by HSPA1A and HSPA1B) and Hsp70B’ (encoded by HSPA6) are induced when the cell is under stress(38). The Hsp70 protein family is a key part in maintaining protein homeostasis and cell survival(39). Moreover, in the spliceosome, the three homologs act as part of the Prp19 complex (Prp19C)(40, 41) which is involved in multiple cellular processes, including splicing, transcription elongation, DNA repair and protein degradation. We speculated that these processes regulate the activity of Prp19C or Hsp70, which further regulates the spliceosome, resulting in a splicing frequency change.
To verify the hypothesis, we empirically defined high and low expression groups based on curated gene lists(41, 42) for each pathway (Supplementary Fig. 7) and conducted differential splicing analysis between the two groups. Among the five pathways, DNA repair reflected significantly more (P < 0.05, one-sided Mann-Whitney U test) AS events (N = 1 197) than other pathways (Fig. 3B), and also more than any of the cytogenetic groups that were previously defined. We further determined the importance of the DNA repair pathway by investigating the three main mechanisms responsible: homologous recombination (HR), non-homologous end joining (NHEJ) and microhomology-mediated end joining (MMEJ)(43). For each pathway, samples were divided into groups based on unsupervised hierarchical clustering of pathway gene expression (Supplementary Fig. 7) and differential AS events were identified. The NHEJ pathway comparison identified significantly more (P < 0.05, one-sided Mann-Whitney U test) AS events (N = 2 999) than the other two pathway comparisons (Fig. 3C). Moreover, samples in the low NHEJ expression group were clustered together (Fig. 3D).
Among the 12 genes in the NHEJ pathway, 11 showed a significant expression difference between the high and low groups with up to a 4-fold change in expression (Fig. 3G and Supplementary Fig. 8B). XRCC4 and LIG4, which together form the DNA ligase IV complex that is essential for NHEJ(44), had the highest fold-change in expression and correlation with splicing frequencies (\(\rho =0.56\) and \(\rho =0.63\) Spearman correlation, Fig. 3E). The expression of the ligase complex showed a consistently high correlation (\(\rho =0.63\)) with splicing frequencies (Fig. 3E) and strengthens the association between NHEJ activity and AS frequency in NDMM.
Samples in the NHEJ-high expression group were enriched for mutations in TP53 (15% vs. 6%, P = 2x10− 5, Fisher’s exact test), biallelic TENT5C (21% vs. 10%, P = 1x10− 9) and biallelic TP53 (10% vs. 4%, P = 3x10− 5) events (Fig. 3G). The combination of these abnormalities with high expression of the NHEJ pathway resulted in significantly more differential AS events than the abnormality alone, for example, biallelic TP53 + high NHEJ compared against samples without biallelic TP53 abnormalities results in more differential AS events than biallelic TP53 alone (Fig. 3F). By comparing the NHEJ-low, -medium and -high groups, a significantly elevated expression of 66 of the 69 documented splicing factors was observed (Supplementary Fig. 8A).
To investigate the association between NHEJ activity and DNA instability, structural variations (SVs) that occurred in samples within the -high and -low group were compared and there was a significant difference between them (mean 87 vs. 19, respectively, one-sided Mann-Whitney U-test P = 6x10− 4). Moreover, the NHEJ-high group were associated with a worse outcome compared to the others (Logrank test statistic = 7.3, P = 0.007), Fig. 3H which remained significant in a multivariate analysis with other well-known risk factors in MM (HR = 1.8, CI=(1.3,2.2), P = 1.6x10− 5), including 1q gain/amp, 17p deletion, high-risk translocations, TP53 mutation, and the previously defined Double-hit, Fig. 3I(5).
Here we bring together the association of DNA instability, an increased number of SVs, enrichment of abnormalities in TP53 and TENT5C, increased activity in the NHEJ pathway, which in turn affects spliceosome component expression and AS of RNA transcripts.
Alternative Splicing Events Are Associated with Outcome Within Cytogenetic Subgroups
We subsequently conducted a survival analysis and evaluated the effect on outcome associated with each AS event both within their genomic subgroup and among all NDMM samples, Supplementary Table 4. Within subgroup comparisons, 26 and 62 events showed an association with either PFS or OS, respectively (P < 0.05 univariate Cox regression). Among them, AS events in the t(4;14) were most frequent accounting for approximately half of those associated with outcome. When well-known high-risk factors were incorporated along with AS events in the analysis, 96 and 62 AS events showed a significant association with PFS and OS, respectively (P < 0.05 multivariate Cox regression). AS events associated with outcome were enriched in the t(14;16) (58% and 77%) and biallelic-TP53 (62% and 78%) groups.
Examples of these events included an AF event in RHOQ within the t(11;14) comparison. RHOQ encodes a GTPase that binds to a variety of effector proteins to regulate cellular responses. The mean PSI of this event was used to divide all t(11;14) samples into two subgroups, and although there was no gene level expression difference between the two t(11;14) groups (Fig. 4A-B), the ‘PSI < mean’ group was associated with a significantly worse OS compared to the ‘PSI > mean’ group (Fig. 4C). Furthermore, ‘PSI < mean’ maintained its association with high-risk in a multivariate analysis (Fig. 4D), suggesting that this AF event can independently define a high-risk subgroup in t(11;14) patients.
Alternative Splicing Events Are Associated with High-Risk Multiple Myeloma
We further investigated the independent prognostic effect of these events in a multivariate model with existing high-risk markers across all newly diagnosed patient samples. The top 200 expressed AS events with the highest variation across the entire population were selected and the mean PSI was used to divide the population into ‘<mean’ and ‘>mean’ groups. A backward elimination approach was used to select the predictors of poor PFS using Cox regression and as a result the final model contained 16 predictors in which 12 of them were AS markers (Fig. 4E) and reflected a high overall prediction performance (C-index = 0.69).
Among the AS markers, eight were associated with unfavorable prognosis while five were indicators of favorable prognosis. We have previously shown that accumulation of poor prognosis markers results in a reduced outcome and so for the ten markers associated with unfavorable prognosis (7 AS markers, NHEJ-high, Double-Hit and gain 1q), we tested whether patients with an increasing number of markers would result in decreased outcome. The samples were divided into three subgroups based on the number of poor prognosis markers they carried: low risk (0–1 markers), medium risk (2–5 markers) and high risk (6–7 markers) groups (no samples contained more than 7 markers). As the samples contained more markers, an additive risk was observed (Fig. 4F) with the high-risk group (N = 27) reflecting a significantly worse median PFS compared to the low-risk group (N = 520) (median PFS = 502 vs. 1 378 days, respectively, P = 1x10− 13, Logrank test). The high-risk group constituted ~ 5% of the entire population and were significantly enriched for gain 1q, Double-Hit, NHEJ-high and five AS markers of unfavorable prognosis (Fig. 4G). Details of the AS events in the model can be found in the Supplementary Data and Supplementary Fig. 6, but included several genes implicated in MM biology, including TANK and CFLAR.
AS Levels are Highly Correlated with Expression of Splicing Factors
To understand the potential roles of splicing factors in regulating AS in MM we systematically studied the association between the expression of 69 splicing factors and the splicing levels of 25 AS events identified in all three translocation subgroups (Fig. 5B). An empirical threshold (\(\left|\rho \right|>0.3\), Spearman correlation) was set to preserve high correlations, from which three correlation networks were constructed: t(11;14) Fig. 5A, t(14;16) (Supplementary Fig. 9A) and t(4;14) (Supplementary Fig. 9B). Edge and node values of overlapped correlation networks can be found in Supplementary Table 5. Among the three networks, most events were regulated by the combinational effects from the up- or down-regulation of splicing factors, while very few events were only associated with a single splicing factor. A large number of pairwise-consistent correlations existed (Fig. 5C) with 55% of correlations in t(11;14) and 35% in t(4;14) samples being consistent with other comparisons. Ten correlations were found to consistently exist among all three networks, eight of which are with an AL event in ZFAS1, whose splicing level was consistently lower in the three translocation comparisons compared to non-translocated samples. For example, the correlation between AL_ZFAS1 splicing and SRSF6 expression remained consistently high in all three translocations (Fig. 5D) as well as in hyperdiploid comparisons (Fig. 5E).
However, other correlations existed either uniquely or in a pairwise fashion. For instance, the correlation between AL_ZFAS1 and SRSF7 was consistently top-ranked among t(4;14) (Fig. 5F) and t(14;16) samples (Fig. 5G), while the correlation between AL_ZFAS1 and SF1 was uniquely found in t(11;14) (Fig. 5H). Nonetheless, different transcripts are their roles of this gene has not been resolved. This AL event occurred between two transcripts: ZFAS1-203 and ZFAS1-210 (Supplementary Fig. 9C), which have different last exons. ZFAS1-210 was annotated as one of the representative transcripts in GENCODE(45) and was significantly higher in t(14;16) while ZFAS1-203 was significantly lower in t(4;14) and t(11;14).