Alternative splicing of ceramide synthase 2 alters levels of specific ceramides and modulates cancer cell proliferation and migration in Luminal B breast cancer subtype

Global dysregulation of RNA splicing and imbalanced sphingolipid metabolism has emerged as promoters of cancer cell transformation. Here, we present specific signature of alternative splicing (AS) events of sphingolipid genes for each breast cancer subtype from the TCGA-BRCA dataset. We show that ceramide synthase 2 (CERS2) undergoes a unique cassette exon event specifically in Luminal B subtype tumors. We validated this exon 8 skipping event in Luminal B cancer cells compared to normal epithelial cells, and in patient-derived tumor tissues compared to matched normal tissues. Differential AS-based survival analysis shows that this AS event of CERS2 is a poor prognostic factor for Luminal B patients. As Exon 8 corresponds to catalytic Lag1p domain, overexpression of AS transcript of CERS2 in Luminal B cancer cells leads to a reduction in the level of very-long-chain ceramides compared to overexpression of protein-coding (PC) transcript of CERS2. We further demonstrate that this AS event-mediated decrease of very-long-chain ceramides leads to enhanced cancer cell proliferation and migration. Therefore, our results show subtype-specific AS of sphingolipid genes as a regulatory mechanism that deregulates sphingolipids like ceramides in breast tumors, and can be explored further as a suitable therapeutic target.


Introduction
Sphingolipids help in maintaining the structural integrity of cell membranes, and also aid in numerous signaling processes in response to different stimuli [1][2][3] . Dysregulation of the sphingolipid pathway is one of the important contributing factors for breast cancer pathogenesis as it is implicated in various aspects of cancer initiation, progression, invasion, metastasis, and drug resistance 4,5 . Any alteration in the enzymes regulating the expression of sphingolipids play a vital role in cancer cell survival and apoptosis 6 . Therefore, dynamic metabolic interconversions among sphingolipids, attuning the cellular signaling mechanisms in cancer cells, make them a desirable yet challenging target for cancer therapy 7,8 .
Ceramides are key precursors of complex sphingolipids that are synthesized via de novo as well as through salvage pathway. De novo pathway involves conjugation of specific fatty acyl chains to sphinganine using ceramide synthases, whereas salvage pathway involves degradation of sphingomyelins and complex sphingolipids to ceramides 9 . Human breast cancer tissues usually have high expression of ceramide species as compared to normal tissue samples 10 . There is an increased expression of C16:0, C24:0, and C24:1 ceramides in malignant tumors as compared to benign and normal tissues, where C16:0 ceramides were associated with positive lymph node status 11 . Estrogen receptor (ER)-positive breast tumors showed high expression of C18:0 and C20:0 ceramides as compared to ER-negative tumor tissues 11 .
Mammalian cells express six ceramide synthase isoforms (CERS1-6) that show similar catalytic mechanism and intracellular localization in the endoplasmic reticulum but exhibit distinct fatty acyl-CoA specificity 12 . As CERSs are key enzymes of sphingolipid metabolism, many genetic, transcriptional, post-transcriptional, and posttranslational regulatory mechanisms regulating CERSs have been identified in various pathophysiological processes 13 . CERSs including CERS2, CERS4, and CERS6 are overexpressed in malignant breast tissues 14,15 Human bladder carcinoma patients with loss of CERS2 mRNA expression have poor prognosis, and is associated with tumor progression and invasion 16 . CERS2 mRNA levels are also found to be low among high-grade meningioma patients who suffer from tumor recurrence more frequently, and have increased risk of death 17 .
CERS2 is the only ceramide synthase that codes for very-long-chain ceramides. As altered expression of CERS2 is associated to different types of cancer 18 , deciphering the mechanisms responsible for dysregulation of CERS2 gene expression is critical for understanding its role in cancer progression. Exon 1 of CERS2 has been shown to regulate the CERS2 transcription as it allows the binding of Kruppel-like factor 6 (KLF6) and zinc finger transcription factor, Sp1 (ref. 19 ). This transcription factor-mediated upregulation of CERS2 is effective in suppression of metastasis of prostate cancer cells. Among different post-transcriptional regulatory mechanisms, miR-133a, miR-221, and miR-222 have been identified to bind with 3′-UTR regions of CERS2, and regulate its expression 20,21 . A truncated CERS2, tumor metastasis suppressor gene-1 (TMSG1), lacking the N-terminal 150 amino acid residues showed very low expression in several metastatic cancer cell lines and tissues compared to the non-metastatic ones 22 . Although other CERS2 splice variants have also been reported, CERS2 alternative transcript lacking exon 8 has not been studied in detail.
Alternative splicing (AS) has emerged as a crucial posttranscriptional regulatory mechanism in many disease conditions, especially cancer 23 . Reorganization and differential expression of specific transcript isoforms are beneficial for cancer pathogenesis, making them a valuable target for cancer therapy 24 . AS also modifies the balance between pro-survival and pro-apoptotic variants of many proteins like p53, caspase 2, FAS, BCL2L1, and caspase 9 (ref. 25,26 ). Recently, our group has shown that a combination chemotherapy induces the AS in glucocerebrosidase β (Gba1), and increases Gba1 expression in murine tumor tissues. This enhanced Gba1 expression catalyzes the degradation of glucosylceramides to ceramides responsible for tumor regression, and lowers the glucosylceramide levels critical for drug resistance 27 .
Recent studies have shown that subtype-specific qualitative and quantitative expression of different transcript isoforms via AS are integral to the molecular portrait of each subtype in breast cancer patients 28,29 . Differential expression of dysregulated transcript isoforms along with prominent isoform switching patterns have allowed stratification of subtypes in breast cancer patients 30 . We hypothesize that different breast cancer subtypes (Luminal A, Luminal B, HER2 + , and Basal) may possess unique RNA splicing signatures for genes of the sphingolipid pathway, and post-transcriptional regulation (alternative splicing) of these genes, especially ceramide synthases, may contribute to breast cancer development.
Therefore, we aimed to identify the AS events of sphingolipid genes in different breast cancer subtypes from TCGA-BRCA dataset 31 . We show that ceramide synthase 2 (CERS2) undergoes a unique AS event where Exon 8 is skipped in Luminal B subtype patients, and validate the differential expression of CERS2 isoforms in Luminal B representative cancer cells and tumor tissues. We further show that this AS event in CERS2 significantly affects survival in Luminal B patients of TCGA-BRCA cohort, and is a poor prognosis factor. Loss of Exon 8 contributes to the lack of catalytic activity and substrate specificity of CERS2 for very-long-chain ceramides. Finally, we demonstrate that loss of Exon 8 reduces the levels of very-long-chain ceramides, thereby affecting the cancer cell proliferation and migration.

Results
Breast cancer patients exhibit subtype-specific AS events RNA sequencing (RNA-Seq) data for TCGA-BRCA controlled dataset was downloaded from Genomics Data Commons (GDC portal, NIH) for Breast Invasive Carcinoma Project 31 . The data for 817 cases of Invasive Ductal Carcinoma (IDC) was classified into five different subtypes (Luminal A, Luminal B, HER2 + , Basal, and Normallike) based on their hormonal character and PAM50 gene expression (Fig. 1A, Data set 1). We analyzed the RNA-Seq data to identify transcriptome-wide splicing events in different breast cancer subtypes compared to Normal-like using a previously described bioinformatic pipeline 32 . We used 'Normal-like' as the reference based on the understanding that the "Normal-like" tumors probably had only a few tumor cells and a large number of normal breast epithelium 33 , and therefore, they have a profile measuring nearest to the normal breast epithelium. Our analysis specifically focussed on two AS events, cassette exon (CE) events where an exon is spliced-in or spliced-out, and intron retention (IR) events where an intron is retained in an isoform under a certain condition (Fig. 1A (Fig. 1C, D). AS in cancer-related genes can alter the desired protein expression and modulate the signal transduction pathways, thereby leading to transformation of cancer cells 34 . Gene enrichment analysis using DAVID showed that subtype-specific genes predicted to undergo AS are enriched for metabolic pathways, cell cycle arrest, transport and trafficking, DNA repair, protein homeostasis, and mRNA transport and splicing ( Supplementary Figures 1-4, Data sets 4-7). GO terms for the predicted AS genes showed biological processes such as mitosis, cell adhesion, GTPase and transcription factor-mediated signaling pathways, translation, and transport ( Supplementary Figures 1 35 . We used event-specific primer pairs for the predicted CE event, and found that levels of proteincoding (PC) and AS transcripts varied in different cell lines (Fig. 1E, F). The level of AS transcript is significantly lower in BT-474 (p < 0.05) and MDA-MB-453 (p < 0.05) cells as compared to MCF-10A cells, but higher in MDA-MB-231 cells (p < 0.001) (Fig. 1F, Data set 2, Supplementary Table 1).
We also observed a CE event targeting Exon 3 of a Ras homolog family member A (RHOA) in basal subtype (p < 0.0003) where Exon 3 is skipped. RHOA is a small GTPase protein known to regulate malignant transformation and motility of cancer cells 36 . Using event-specific primer pairs, we observed differential expression of PC and AS isoforms in all cell lines (Fig. 1G, H, Data set 2, Supplementary Table 1). The level of AS transcripts that include the predicted event (Exon 3 exclusion) is significantly higher in all subtype-specific cell lines with MDA-MB-231 cells showing the most abundant AS1 transcripts (p < 0.05), followed by BT-474 cells (p < 0.05), and MDA-MB-453 cells (p < 0.01) as compared to MCF-10A cells (Fig. 1H).

AS of sphingolipid pathway genes provide subtypespecific signature
Analysis of subtype-specific CE and IR events specifically for genes of sphingolipid pathway revealed 18 CE events in 10 genes and 5 IR events in 4 genes in Basal subtype, and 14 CE events in 10 genes and 3 IR events in 2 genes for Luminal A subtype ( Fig. 2A-D). We also observed 15 CE events in 10 genes and 3 IR events in 1 gene in Luminal B subtype, and 16 CE events in 10 genes and 4 IR events in 2 genes in HER2 + subtype ( Fig. 2A-D). Maximum number of sphingolipid genes that undergo AS were predicted to have one event while few selected genes (see figure on previous page) Fig. 1 Breast cancer patients exhibit subtype-specific alternative splicing (AS) events. A Schematic representation of the strategy for identifying alternatively spliced (AS) events in breast cancer subtypes (TCGA-BRCA dataset). B Number of AS events caused by differential usage of cassette exons (CE) or by intron retention (IR) for Luminal A, Luminal B, HER2 + , and Basal breast cancer subtypes, and total number of genes that they correspond to. Statistical test Irwin-Hall p-value summarization was used on individual splicing events. Details of samples for each subtype are in Data set 1 in supplementary information.  showed two or three events (Data set 8,9). Sphingolipid genes undergoing AS showed some common signatures for all the subtypes as well as a unique signature for each subtype (Fig. 2E, Data set 8,9). For example, common AS targets among subtypes include CE events in Fatty acid elongase 5 (ELOVL5), CERS5, and Hexosaminidase subunit alpha (HEXA) (Fig. 2F). In contrast, CE event in CERS4 is specific for Basal subtype (Fig. 2F), whereas CE event in CERS2 is unique for Luminal B subtype (Fig. 2F). Among different IR events, IR event in CERS5 is common between Luminal A and Basal subtypes, and there is an IR event in CERS4 that is unique for Basal subtype (Fig. 2G). As ceramide is the key precursor of sphingolipid pathway, AS events in CERS genes may regulate the subtypespecific synthesis of ceramides and other sphingolipids in breast cancer patients. Subtype-specific AS events in different CERS genes such as CERS2 in Luminal B subtype and CERS4 in Basal subtype suggest a strong subtypespecific post-transcriptional regulation of CERS genes that may play a role in breast cancer pathogenesis (Fig. 2F, G; Data set 8,9).  Figure 10), and is probably subjected to nonsense-mediated decay (NMD) through the NMD pathway 37 . Therefore, we focused on the rest of our study on the full-length PC and the alternatively spliced, AS1, transcript of CERS2.
We validated the AS1 event in Luminal B BT-474 cells, and compared it with HER2 + MDA-MB-453 and normal MCF-10A cells using specific primer pairs of similar efficiency for PC and AS1 transcripts (Supplementary Table 1). The PC transcript is highly abundant in MCF-10A cells in comparison to BT-474 and MDA-MB-453 cells (Fig. 3C), whereas the AS1 transcript is significantly more in BT-474 cells and not in MDA-MB-453 cells as compared to MCF-10A cells (Fig. 3D). As the abundance of CERS2 AS1 transcript is very low, it was not possible to use the canonical approaches of using one primer pair to identify both PC and AS1 transcript isoforms. We, therefore, used a strategy based on published methods, that has shown the use of boundary-spanning (exon-exon junction specific) primers to quantify the different isoforms of a transcript that differ in their abundance [38][39][40] . However, using RNA isolated from cells overexpressing PC and AS1 transcripts, both transcripts can be identified using the same primer pair. We used forward primer from exon 7 and reverse primer from exon 9, thereby giving a 129 bp reduction in band size for AS1 overexpressing BT-474 cells (Supplementary Figure 11).
CERS2 is a transmembrane protein of 380 amino acids with a 199 amino acids long Tram-Lag-CLN8 (TLC) domain (Fig. 3A). Lag1p motif, a conserved stretch of 52 amino acids (203-254), is a part of the 151 residues (151-301 amino acids) span of TLC domain that is predicted to be the minimal region required for acyl chain specificity 41,42 (Fig. 3A, Supplementary Figure 12). Lag1p motif, with two conserved histidine residues, is also suggested to be crucial for CERS catalytic activity and substrate binding (Supplementary Figure 13) 42 . Structurefunction relationships of CERS proteins have found that Lag1p motif confers catalytic activity, and assists in substrate binding to synthesize very-long-chain (C22:0, C24:0, C24:1) ceramides for CERS2 43 . As Exon 8 of CERS2 corresponds to a stretch of 43 amino acids (205-247) of Lag1p motif, AS1 transcript lacking Exon 8 would code for a protein devoid of almost the entire Lag1p domain (Supplementary Figure 12).
Immunoblot of MCF-10A, BT-474, MDA-MB-231, and MDA-MB-453 cell lysates with CERS2 antibody showed a higher expression of CERS2 protein (band closer to 43 kDa) arising from full-length PC transcript (Fig. 3E). Among different breast cancer cell lines, BT-474 cells showed the lowest CERS2 expression arising from fulllength PC transcript. Repeated attempts to visualize the protein product from AS1 transcript (5 kDa smaller than PC protein) on an immunoblot were unsuccessful due to low abundance of the AS1 transcript as compared to the PC transcript. However, we could detect both PC and AS1 protein products translated from overexpression of PC and AS1 transcripts in HEK 293T cells (Supplementary Figure 14), thereby validating that AS1 transcript gets translated. If AS1 transcript indeed gets translated in BT-474 cells, corresponding mRNAs will be associated with polysomes. Therefore, we analyzed the distribution of fulllength PC and AS1 mRNAs in the polysome fractions of BT-474 cells. As expected, full-length PC mRNA is distributed throughout the polysomes. Interestingly, we observed that AS1 mRNA is also enriched specifically in the polysome fractions (fractions 5-12) compared to free RNA or monosome fractions (fractions 1-4) (Fig. 3F).
To further validate the existence of stable proteins translated from PC and AS1 transcripts, total endogenous protein from BT-474 cells was immunoprecipitated using anti-CERS2 antibody-conjugated agarose beads. The immunoprecipitated complex was resolved on SDS-PAGE, and the gel region (between 45 and 35 kDa) was excised in form of three gel slices as the immunoblot probed with anti-CERS2 antibody showed the bands in this region (Fig. 3E, Supplementary Figure 15A).
The amino acid stretches corresponding to the junction of Exon 7-8 in the protein translated from PC transcript and corresponding to the junction of Exon 7-9 derived from AS1 protein sequence are rich in lysine (K) and arginine (R) residues. Therefore, gel slices corresponding to CERS2 proteins were in-gel digested with chymotrypsin yielding 9-14 amino acid peptides of PC and AS1 protein sequences. Fractions corresponding to three gel slices were analyzed by LC-MS/MS using Informationdependent acquisition (IDA) method and processed by Peptide Shaker software that confirmed~51% coverage for CERS2 PC/FL (full-length) and CERS2 AS1/SF protein (spliced form) (Fig. 4A Figure  15C). Therefore, these results confirm that both full-length PC and AS1 endogenous mRNAs are translated.

Luminal B tumor tissues show a higher expression of alternatively spliced CERS2
Next, we validated the predicted AS1 event of CERS2 in Luminal B tumor tissues excised from Indian patients to strengthen the pathological significance of this event across race and ethnicity. Tumor and adjoining normal breast tissues from 10 Luminal B cancer patients were collected from patients undergoing surgery (Supplementary Table 3). Expression of CERS2 PC and AS1 transcripts were determined by qRT-PCR in tumor and adjoining normal tissue pairs using event-specific primer pairs (Supplementary Table 1). We observed a~3-fold increase in the expression of AS1 transcripts in tumor tissues as compared to normal tissues (p = 0.09), whereas expression of PC transcript was not markedly different ( Fig. 5A). We also quantified the CERS2 expression using immunohistochemistry, and found a~1.7-fold increase in CERS2 protein expression in tumor tissues as compared to normal tissues (Fig. 5B). On staining the tissue sections with anti-ceramide antibody, we also observed an increase in the level of total ceramides (Supplementary Figure 16). To show that Luminal B tumor tissues express the alternate protein product coded by AS1 transcript, we performed western blot analysis using protein isolated from 4 Luminal B tumor tissue samples and adjoining normal tissue collected from patients (Fig. 5C, Supplementary Table 3). We observed the expression of protein from both PC (marked by blue asterisk) and AS1 (marked by red Asterix) coded transcripts in tumor tissue (Fig. 5C). However, adjoining normal tissue samples show less expression of both PC and AS1 protein products as apparent from a longer exposure of the western blot (  tissues of ER-positive patients 10,11 . Therefore, the increase in the expression of CERS2 protein in Luminal B tumor tissue, in all probability, arises from full-length PC and elevated AS1 transcripts.

Exon 8 skipping event in CERS2 forecasts a poor prognosis
Survival analysis, in general, involves the correlation of differential gene expression with clinicopathological data of patients 44 . However, CERS2 [p = 0.012, log fold change (FC) of 0.94] did not appear as a differentially expressed gene in Luminal B subtype tumors of TCGA-BRCA dataset (Data set 10). Recently, it has been reported that SURVIV (Survival analysis of mRNA Isoform Variation) method using alternative splicing-based survival predictors outperform gene expression-based survival predictors 45 . Therefore, we hypothesized that differential expression of AS1 transcripts of CERS2 can be a critical contributor to the patient survival over CERS2 gene expression. Accordingly, we analyzed the effect of AS1 (Exon 8 skipping) of CERS2 on survival. Exon 8 lies between 150967202 and 150967074 position of CERS2 on the negative strand of chromosome 1 in humans. Using an in-house script (see "Materials and methods" section), we found three flanking and two skipping junctions associated with Exon 8 (Supplementary Table 4). Upstream of Exon 8, there are two flanking junctions, JUNC51134 (connecting Exon 8 to 7) and JUNC511135 (connecting Exon 6 to 8). There is only JUNC511126 downstream to Exon 8 connecting Exons 8 and 9. There are two skipping junctions to Exon 8, JUNC511127 that connects Exon 7 to 9, and JUNC51128 that connects Exon 6 to 9 (Fig. 5E). So, we considered Exon 8 and its associated terms (metafeatures such as flanking junctions and skipping junctions), and performed survival analysis on TCGA-BRCA Luminal B cohort patients (Data set 1) using both univariate and multivariate Cox-PH model. Exon 8 and associated terms, individually, did not show any effect on survival in univariate Cox-PH analysis (p = 0.11-0.61) (Fig. 5F). However, multivariate Cox-PH analysis showed that all meta-features (Exon 8 (EX419712), JUNC51127, JUNC51128, JUNC51135, JUNC51134, and JUNC51126) combined together significantly affect patients' overall survival (Fig. 5F). The Cox regression p value for all these terms was found to be highly significant (p < 2e−16). Concordance statistics (C-index) gives probability of the prediction accuracy of a logistic regression model 46 . Our model utilizing Exon 8 and its meta-features is a good predictor of patient survival, as confirmed by its C-index (C = 1, s.e. = 0) (Fig. 5F). The result of multivariate Cox analysis shows a decrease in hazard ratio, and suggests that exclusion of Exon 8 in breast cancer patients is a poor prognosis factor (β-coefficient = −229) (Fig. 5F).

Overexpression of CERS2 (AS1) transcript attenuate levels of very-long-chain ceramides
To functionally validate the effect of AS1 transcripts on CERS2 catalytic activity, we overexpressed PC and AS1 transcript-encoding cDNAs in BT-474 cells by transfecting the cells with only vector or CERS2 OE (PC) or CERS2 OE (AS1) constructs. CERS2 OE (PC) construct allows the overexpression of protein-coding transcript whereas CERS2 OE (AS1) construct helps in overexpression of AS1 transcript. Overexpression of PC and AS1 transcripts led to the selective increase in the corresponding CERS2 transcript levels (Fig. 6A). We did not observe any change in expression of other CERSs on overexpression of CERS2 PC and AS1 transcripts (Fig. 6A). Immunoblot analysis confirmed the overexpression of protein products of CERS2 OE (PC) or CERS2 OE (AS1) constructs after 36 h of transfection (Fig. 6B), where we observed the protein band translated from AS1 transcripts in CERS2 OE (AS1) overexpressed BT-474 cells (Fig. 6B). The endogenous CERS2 PC gets also increased in the CERS2 OE (AS1) overexpressing cells, that may be due to a recalibration effort by the cells. CERS2 is known to interact and form dimer with other CERS enzymes like CERS5 and CERS6 (ref. 47 ). We, therefore, checked for any alterations in expression of CERS1, CERS5, and CERS6, and did not observe any appreciable change due to overexpression of CERS2 (Supplementary Figure 18).  Figure 19). Therefore, these results suggest that unlike CERS2 OE (PC) cells, CERS2 OE (AS1) BT-474 cells are unable to increase the levels of very-long-chain ceramides.
Overexpression of CERS2 (AS1) transcript fails to suppress proliferation and migration of BT-474 cells as effectively as PC transcript In addition to its fundamental function as a verylong-chain ceramide synthase, CERS2 has been shown to inhibit tumor invasion and metastasis [49][50][51] . Overexpression of CERS2 induces significant suppression in proliferation and migration of breast, prostate, and bladder cancer cells 52 (Fig. 7A). Therefore, these results suggest that AS1 transcript overexpressing cells have higher proliferation rates over PC transcript overexpressed cells, and overexpression of PC transcripts inhibit the proliferation of Luminal B subtype-specific cancer cells.
To elucidate the role of Lag1p domain on cell migration, we quantified the number of migrating cells using scratch wound assay for 36 h after overexpression of PC and AS1 transcripts, and compared with BT-474 and BT-474 (+Vector) cells. CERS2 OE (AS1) cells showed~1.77fold (p < 0.0001) increase in number of migrating cells compared to CERS2 OE (PC) cells (Fig. 7B, C). CERS2 OE (PC) cells showed~2.15-fold decrease (p < 0.00001) in number of migrating cells compared to BT-474 (+Vector) cells whereas CERS2 OE (AS1) cells did not show any significant change in number of migrated cells compared to BT-474 (+Vector) cells (Fig. 7B, C). There was >55% increase (p < 0.0001) in number of CERS2 OE (AS1) migrating cells as compared to CERS2 OE (PC) cells (Fig.  7B, C). Therefore, these results suggest that overexpression of PC transcript and increased levels of verylong-chain ceramides help in reduced migration of (see figure on previous page) Fig. 5 CERS2 differential isoform expression predict poor prognosis in Luminal B patients. A Comparison of differential isoform expression of CERS2 in tumor (n = 10) and matched control (n = 10) breast tissues from Luminal B patients using full-length PC and AS1 transcript specific primers show increased AS1 expression in tumor tissue. A pictorial representation of exon/intron position and primer sites are also shown. B Immunofluorescence images and quantification of CERS2 (mean ± SD, n = 7) show elevated levels in tumor tissue as compared to normal breast tissue sections. C, D Immunoblot showing the expression of proteins corresponding to PC and AS1 transcripts from 4 pairs of Luminal B tumor and adjoining normal tissue samples (C), and 12 patient tumor samples (D). Protein corresponding to PC transcript is marked by blue asterisk, and to AS1 transcript by red asterisk. E, F Survival Analysis on Luminal B patient data (TCGA-BRCA) using univariate and multivariate Cox proportional hazards (Cox-PH) regression. Inclusion of meta-features (exon 8, flanking junctions 51135, 51134, 51126 and skipping junctions 51127, 51128) (E), in multivariate analysis significantly affect patient survival (β coefficient = −229, p < 2e−16, and C-index = 1) (F). Data in A and B were analyzed by twotailed unpaired Student's t-test. To confirm that the difference in proliferation and migration between CERS2 OE (PC) and CERS2 OE (AS1) cells is not due to cell death, we performed cell death assay by Annexin-FITC/Propidium Iodide assay. Overexpression of CERS2 OE (PC) and CERS2 OE (AS1) transcripts induced some apoptosis in BT-474 cells. There was only a~5% increase in apoptosis in CERS2 OE (PC) overexpressed cells compared to CERS2 OE (AS1) cells (Supplementary Figure 20), whereas we observed 30% decrease in proliferation and >55% decrease in migration in CERS2 OE (PC) cells compared to CERS2 OE (AS1) cells. Therefore, these results reveal a significant contribution of skipping of Exon 8 in CERS2 towards cancer cell proliferation and migration in Luminal B subtypespecific cancer cells.

Discussion
Genome-wide approaches reveal that cancer development and progression involves massive changes in AS events affecting the tumor proteome. This is evident from the large repertoire of subtype-specific CE and IR events scored from our study on TCGA-BRCA RNA-Seq data. Events specific to sphingolipid genes predict that AS contributes distinctly to reprogramming of sphingolipids in breast tumors of different subtypes. This may be the possible link between high levels of sphingolipids like ceramides and sphingomyelins in tumor tissues, contrary to their antiapoptotic role established in in vitro and in vivo studies. Identification of unique and common AS events in ceramide synthase genes suggest that AS may be responsible for chain-length-specific accumulation of ceramides in subtype-specific breast tumor tissues. The increase in ceramides in response to AS events may therefore contribute to the lipid metabolic reprogramming required during tumor progression 53,54 .
Validation of one of the predicted AS transcripts for CERS2 in Luminal B-specific BT-474 cells showed the expression of AS1 transcript lacking Exon 8 along with full-length PC transcript. Both transcripts get translated endogenously as observed by the polysome assay and LC-MS/MS-based MRM identification of corresponding peptides spanning the junction between Exon 7-9 for AS1 and Exon 7-8 for PC-derived proteins. Lag1p domain is an essential part of TLC domain of CERSs that helps in selecting different acyl-CoAs to generate chain-lengthspecific ceramides. Recently, it has been shown that a 11 residue region (amino acid 299-309) in the last luminal loop between the putative transmembrane domains 5 and 6 confers acyl chain substrate specificity 55 . However, acyl-CoA specificity assay using two sets of chimeric proteins, where 11 residue region was swapped between CERS2 and CERS5 showed that CERS5 (299-309 CERS2) can utilize verylong-chain acyl-CoAs but CERS2 (299-309 CERS5) did not show any activity towards both short and long-chain acyl-CoAs. Therefore, apart from the 11 residue domain, there may be other determinants like Lag1p domain assigning the activity and substrate chain specificity to CERS enzymes. Our study shows that CERS2 OE (AS1) lacking almost the entire Lag1p domain loses its specificity in tethering the very-long-chain acyl-CoAs to the sphingosine backbone compared to CERS2 OE (PC) that have an intact Lag1p domain. Therefore, in BT-474 cells, CERS2 AS1 protein probably does not disrupt the activity of the PC-derived protein but disrupts the balance between very-long-chain ceramides and other short/long-chain ceramides. Thus, this post-transcriptional regulatory mechanism possibly disturbs the relative balance of different chain-length-specific "ceramide pools" in tumor tissues 56 .
A balance of short-chain (C16:0), long-chain (C18:0, C20:0), and very-long-chain (CC22:0, C24:0, C24:1) ceramides is critical for determining the fate of cancer cells depending on the cancer type 57 . Overexpression of CERS4 and CERS6 in breast and colon cancer cells, leading to an increase in short (C16:0) and long-chain ceramides (C18:0 and C20:0), inhibited cell proliferation and increased apoptosis 58 . At the same time, accumulation of long-chain C18:0 ceramides in human head and neck squamous cell carcinoma suppressed the proliferation of cancer cells, whereas short-chain C16:0 ceramides induced proliferation 4 . C16:0 ceramides have been shown to play a pivotal role in cell migration during epithelial to mesenchymal transition in different cancer cell lines 59 . In contrast, overexpression of CERS2, that synthesizes very-longchain ceramides, has been demonstrated to inhibit cancer cell proliferation and migration in prostate, breast, and bladder cancers 60,61 . Luminal B-specific BT-474 cells overexpressing PC transcript of CERS2 showed a significant decrease in cell proliferation and migration due to an increase in levels of very-long-chain ceramides. In contrast, overexpression of alternatively spliced CERS2 OE (AS1) did not alter the migration and proliferation of BT-474 cells significantly, reflecting the effect of a disbalance in the level of very-long-chain ceramides. Therefore, disequilibrium between short, long, and verylong-chain ceramides and the balance between AS transcripts of critical genes regulating the ceramide metabolism adds on to the myriad of factors that lead to cancer progression.
Luminal B subtype tumors are marked with increased expression of proliferation-related genes, more aggressive clinical behavior, and significantly unfavorable prognosis as compared to Luminal A subtype [62][63][64] . In spite of better prognosis after treatment, there is a higher proportion of local recurrence as well as a persisting risk of metastasis in Luminal B patients in comparison to non-Luminal cancers 65 . Most commonly-used survival methods for retrospective studies of clinical data undermine the contributions of differential AS of genes. SURVIV, a statistical method, identifies and associates mRNA isoform variation with survival time in both censored and uncensored large-scale RNA-Seq data. Exon inclusion events were found to be significantly correlated with survival time in TCGA ductal carcinoma patients when SURVIV was used. Combined effect of Exon 8 skipping event and associated meta-features of CERS2 in Luminal B patients from TCGA-BRCA cohort shows that exclusion of Exon 8 of CERS2 predicts poor survival. Therefore, contribution of post-transcriptional regulatory events needs to be considered for better interpretation of clinical outcome of cancer patients. It is prudent to mention that our data validated for Luminal B subtype will be further strengthened by future studies in Indian patient tissues of other subtypes.
In summary, our results show that AS of sphingolipid genes provide a unique signature for breast cancer subtypes. Luminal B-specific AS of CERS2 reduces the levels of very-long-chain ceramides that allows enhanced proliferation and migration of cancer cells, and a poor prognosis outcome for patient survival. Therefore, identification of post-transcriptional regulatory mechanisms responsible for sphingolipid metabolism can provide new therapeutic strategies for combating tumor progression.

Alternative splicing (AS) analysis
Differential alternative splicing (AS) in Luminal A, Luminal B, HER2 + , and Basal tumor tissues in comparison to normal-like tissues were identified by Exon Pointer (EP) and Intron Pointer (IP) algorithm as per published protocol using in-house pipeline 32 .

GO enrichment analysis
Enrichment analysis of all the genes predicted by AS analysis for CE and IR events was performed using the Database for Annotation, Visualization, and Integrated Discovery v6.8 (DAVID). All default parameters were used, and cut off was kept at p-value ≤ 0.05 (Fisher's exact test) and/or false discovery rate (FDR) of ≤15.0. The functional association was done using Biological processes (GO_BP_Direct), KEGG Pathways, and Protein domains. The raw data and values are provided in Data sets 4-7.

Survival analysis
We have designed and used a novel strategy to predict survival analysis using exon/intron level RNA expression counts combined with associated meta-features in a given splicing condition, i.e., flanking junction and skipping junction to the respective exons/introns. Since we are interested to find the cumulative and individual effect of two or more terms (exons and junctions) in a gene, we applied multivariate Cox regression model and univariate Cox regression model, respectively 66,67 . The resulted term with multivariate Cox-PH and univariate analysis with pvalue <0.05 were considered as significantly affecting the survival. Multivariate Cox-PH analysis gives individual pvalue for exon/intron and its associated terms (flanking and skipping junction). To find the combined effect of given exon/intron and associated terms, these p-values can be summarized using Irwin-Hall method as follows for a cassette event: where, (k ∈ H e ) denotes the set of terms (exon e, flanking junctions, and skipping junctions) where flanking junctions are expressed with exon e when included into a transcript for a given gene k, and skipping junctions are expressed when exon e is not included into transcript. For an intron retention event, (k ∈ H e ) will include intron e and skipping junction. For the event described in results (Fig. 4C), the Eq.
(2) will be as follows; In case the null hypothesis is true, the sum of uniform distribution will behave as uniform distribution according to Irwin-Hall method 68 . Given the X statistic is smaller than one, the corresponding probability distribution function (piecewise polynomial function) would be; where x is the summation of the p-values mentioned in Eq. (2). The probability to obtain a combined p-value is given by: where p-value combine is the summarized p-value given for exon or intron along with associated terms upon significantly affecting the survival.

Transfection protocol
BT-474 (or HEK 293T) cells (~2.5 × 10 5 /well) were seeded in a six-well plate in complete DMEM media. After 24 h, media was carefully removed from each well, and cells were washed with DPBS (Sigma, D5652). Cells were then incubated with a transfection mix comprising 1000 µL Opti-MEM (Gibco, 31985070), 5 µL Lipofectamine (Invitrogen, 11668019), and 2 μg of plasmid PBBL-FLAG vector (control vector) or CERS2 OE (PC) or CERS2 OE (AS1) plasmid in respective wells. After 6 h of transfection, media was replaced with antibiotic-free DMEM media containing 10% FBS, and cells were incubated for 24-48 h. Media was removed, and cells were rinsed carefully with 1X DPBS. To collect the cells for RNA isolation, 1 mL of RNAiso Plus (DSS Takara, 9109) was added to each well, and incubated for 5 min. The cells were then aspirated slowly using a pipette. To collect the cells for protein and lipid isolation, cells were scraped from the wells. The cell pellets were rinsed with DPBS by centrifugation at 2500 rpm for 5 min for three times. Pellets were then air-dried and stored at −80°C.

Quantitative real-time PCR
Frozen cell pellets in RNAiso Plus were homogenized in Trizol (Qiagen, 79306). RNA was purified by phenol: chloroform:isoamyl alcohol extraction followed by ethanol precipitation. RNA concentration was determined using NanoDrop 2000 (ThermoScientific, USA). The integrity and quality of ribosomal 28S and 18S were determined on the agarose gel. Any traces of genomic DNA were removed by DNase (Invitrogen, Am2238) treatment for 30 min at 37°C followed by heat inactivation using 50 mM EDTA (Sigma, E5134). cDNA synthesis and real-time PCR were done as described previously 27 . Relative quantitation of gene expression was done using β-actin as the endogenous reference gene for normalization. All primer sequences used for RT-PCR are listed in Supplementary Table 1.

Immunoprecipitation and mass spectrometry
BT-474 cells were lysed in RIPA lysis buffer, and total protein was isolated and quantified using BCA protein estimation kit. Protein equivalent to 1 mg (in 1000 μL lysate) was mixed with 20 μL (equivalent to 10 μg) of LASS2/CERS2 antibody-conjugated to agarose beads (SantaCruz, sc-390745 AC) at 4°C for overnight. Bound protein complexes were centrifuged, washed two times with 500 μL buffer, and the supernatant containing unbound protein was discarded. Anti-CERS2 antibodyimmunoprecipitated complexes were eluted in 30 μL SDS-PAGE Laemmli sample buffer. Protein complexes were resolved on 15% SDS-PAGE, and stained with 0.1% Coomassie stain (Sigma, 27816).
The MS/MS spectrum raw files (.wiff) were converted in .mgf format, and a search output file was created using the "SearchGUI 3.3.16" platform coupled with Andromeda search engine. Data were searched against the database that consisted of full-length protein-coding (PC/FL) and spliced form (AS1/SF1) of CERS2 protein sequence (Supplementary Figure 7). The search parameters for identification of interested peptide sequences of CERS2 protein from the digested samples were as follows: (a) chymotrypsin as a specific proteolytic enzyme (with up to two missed cleavages), (b) peptide mass error tolerance of 100 ppm, (c) fragment mass error tolerance of 0.  Table 2) generated from both PC/FL and AS1/SF-derived parent ions as mentioned above to construct a multiple reaction monitoring (MRM) method. The collision energy for fragmentation was calculated based on the straight line equation, CE = [Slope] × (m/z) + [intercept], where slope = 0.0625 and intercept = −3 were used for doubly charge peptides as per the suggestion given by the manufacturer. The peptides (Exon 7-8 junction in PC and Exon 7-9 junction in AS1) were identified by MRM scan on a triple quadrupole/linear ion trap mass spectrometer (6500 QTRAP, SCIEX, USA) coupled to a LC system. A Kinetex ® C18, 2.1 × 50 mm column (Phenomenex®, 00B-4601-AN) with a particle size of 5.0 μm was used with 5 μL injection volume, and 40°C oven temperature. The total run time optimized was 10 min. Solvents used were Solvent A (100% water with 0.1% formic acid) and solvent B (100% acetonitrile with 0.1% formic acid) at 0.3 mL/min flow rate. The gradient program used for separation was 5% B-50% B from 0 to 6 min, 50% B-80% B from 6 to 8 min, holding 80% B up to 8.5 min, 80% B-5% B from 8.5 to 9.0 min and keeping 5% B till 10 min. Data analysis and quantification were performed by Analyst ® (version 1.7; SCIEX, USA) and MultiQuant™ software (version 3.0.2; SCIEX, USA). Inclusion criteria include the women patients of all age (18-85 years) of any socioeconomic status who have given consent for tissue collection, patients with operable breast cancers (Stages I, II, IIIa) who will undergo adjuvant therapy, and the patients with estrogen receptor-positive (ER + ), progesterone (PR + ), and HER2 + status. Exclusion criteria include patients who are undergoing neo-adjuvant chemotherapy, and patients who are not fit to undergo surgery. Details of patients and tumor signature are mentioned in Supplementary Table 3.

RNA isolation and quantitative real-time PCR from tumor tissues
Patient tissue samples (~20 mg) in Allprotect Tissue Reagent (Qiagen, 76405) were homogenized in QIAzol Lysis reagent (Qiagen, 74804) using mortar pestle. After homogenization, the samples were incubated at room temperature for 5 min. Phase separation was done by addition of chloroform followed by centrifugation. The aqueous layer was collected, and ethanol was added. The samples were then applied to Qiagen RNeasy spin column (Qiagen, 74804), and centrifuged, washed, and eluted using 30 µL of RNase-free water. The RNA concentration was quantified using NanoDrop 2000 (ThermoScientific, USA) spectrophotometer at OD 260 , and purity was determined by OD 260/280 ratio. cDNA synthesis and realtime PCR were done as described above.

Immunohistochemistry
Patient tumor and normal tissues frozen in Allprotect tissue reagent were cut out and thawed. Tissues were soaked in 4% paraformaldehyde solution (Thomas Baker, 81847) for 24 h, put in cytomatrix (ThermoScientific, 6769006), and frozen at −80°C. Frozen tissue samples were then subjected to cryo-sectioning using cryotome to get 6-8 µm thin sections. The tissue sections were stained with CERS2 (1:400 dilution) (SantaCruz, sc-390745) or Ceramide antibody (1:10 dilution) (Sigma, C8104) using the protocol previously described 69 . Confocal imaging of the samples was performed with Leica TCS SP8. The sections were visualized at ×40 oil immersion using LAS AF software. Z-stacking was performed, and images were acquired for each section. The images were processed using LAS X software. Image Quantification data were analyzed using two-tailed unpaired Student's t-test.

LC-MS/MS method and data analysis
Cell pellets were washed in ice-cold PBS, harvested, and centrifuged. The cell pellets were resuspended in PBS (200 µL), and homogenised by sonication. An aliquot (20 µL) was taken for protein estimation by BCA protein estimation kit. One milligram protein equivalent of each sample was used for lipid isolation using protocol described previously 69 . Only organic phase (OP) samples were extracted for quantitation of ceramides 69 .
OP sphingolipids were resuspended and separated using solvent A (99:1, water/formic acid) and solvent B (99:1, methanol/formic acid), both with 5 mM ammonium formate. The resulting OP extracts for each sample were sonicated for 15 sec, vortexed, and centrifuged at 13,000 rpm for 5 min. The clear supernatant was transferred to the autoinjector vial for LC-MS/MS analysis. Ceramides were analyzed using high-pressure UHPLC liquid chromatography (ExionLC AC, SCIEX, USA) coupled to a hybrid triple quadrupole/linear ion trap mass spectrometer (4500 QTRAP, SCIEX, USA). A Kinetex ® C8, 2.1 × 50 mm column (Phenomenex®) with a particle size of 1.7 μm was used, and oven temperature was kept at 60°C. Total run time optimized was 10 min. Separation method used for ceramides was, Buffer A and B (20:80) for 0-1.2 min, 100% buffer B for 1.2-6.5 min, buffer A and B (20:80) for 6.5-10 min.
Ceramides were monitored using scheduled multiple reaction monitoring (MRM) with preidentified retention time. An MRM to enhanced product ion (EPI) scan was used. All parameters used were as described previously 69 . A standard curve was generated for absolute quantitation of ceramides. Ceramides in each sample were quantified using the same protocol as described previously 69 . Five technical replicates were run for each independent biological replicate (n = 4). Data were analyzed by ANOVA test.
CERS2 in vitro activity assay Preparation of membrane fractions HEK 293T cells were transfected with CERS2 OE (AS1) and CERS2 OE (PC) plasmids as mentioned above. After 36 h, the membrane fraction of HEK 293T cells was prepared 48 . Cells were suspended in buffer A [50 mM HEPES-NaOH, pH 7.4, 150 mM NaCl, 10% glycerol, 1 mM DTT, 1 mM PMSF, and a 1X complete protease inhibitor cocktail], and lysed by sonication (30 amp, three cycle of 20 s on and 30 s off cycle). Homogenate (20 μL) were taken for protein estimation. Equal quantity of all samples (1.0 mg of protein equivalent) were processed. Cell debris were removed by centrifugation at 300 × g, 4°C for 5 min. The supernatant was centrifuged at 100,000 × g (Optima XPN-100 Ultracentrifuge, Beckman Coulter, USA) at 4°C for 30 min. The pellets (total membrane fraction) were suspended in buffer A.
Estimation of Ceramide-d7 using LC-MS/MS LC-MS/MS was performed using the same conditions as mentioned above. Ceramides were monitored using Scheduled MRM. All parameters are described in Supplementary Table 5. Data analysis and quantification were performed using MultiQuant software (version 3.0.2; SCIEX, USA). Six technical replicates were run for each independent biological replicate (n = 4). Data were analyzed using ANOVA test.

Proliferation assay
BT-474 cells (~2.5 × 10 5 /well) were seeded in a six-well plate as mentioned above. After 24 h of incubation, cells were transfected (according to above mentioned procedure) with control vector, CERS2 OE (PC), and CERS2 OE (AS1) plasmids, and incubated for 24 h. Transfected BT-474 cells (5000 cells/well) were then seeded in a 96 well plate for 24, 48, and 72 h in complete DMEM media, and incubated at 37°C. Cell proliferation was quantified using MTT assay 70 . Absorbance reading was taken at 590 nm in iMark TM Microplate Absorbance Bio-RAD plate reader at different time points (24,48, and 72 h). The experiment was performed five independent times, and data were analyzed by two-way ANOVA.

Scratch wound migration assay 71
BT-474 cells (~2.5 × 10 5 /well) were seeded in a six-well plate as mentioned above. After 24 h of incubation, cells were transfected (as mentioned above) with control vector, CERS2 OE (PC), and CERS2 OE (AS1) plasmids, and incubated further at 37°C for 24 h. As cells attain around 90% confluency, the media was aspirated, and a scratch was made by sterile 200 μL pipette tip. Boundaries on each side of the scratch were demarcated by a borderline as shown in Fig. 7C. The wells were washed gently with 1X DPBS twice to get rid of the scraped floating cells, and complete DMEM media was added to the wells. The scratch was imaged immediately for the zero-time point at ×10 magnification in bright field using an inverted microscope (Nikon Inverted LED microscope, Eclipse Ts2, Minato, Tokyo, Japan), and cells were further incubated under similar conditions. The scratch was monitored, imaged, and documented at 6, 18, 24, and 36 h at the same position. A mark was made on the plastic plate as a reference point to ensure that the same area was imaged every time over the course of the experiment. The number of cells migrating into the scratch wound was monitored, and counted using Capture 2.0 image software. The experiment was performed five times independently, and data were analyzed by two-way ANOVA.
Annexin V-PI apoptosis assay BT-474 cells (1 × 10 5 ) were seeded in 24-well plate and incubated for 24 h. Cells were transfected with control vector or CERS2 OE (PC) or CERS2 OE (AS1) plasmids in respective wells as mentioned above. After 48 h, cells were trypsinized and washed with 1X PBS followed by a second wash with 1X binding buffer. Cells were further resuspended into 100 µL of 1X binding buffer. Mixture of Annexin V-FITC (5 µL) and propidium iodide (PI) (5 µL) (BD Biosciences, 556547) were added followed by incubation at room temperature for 15 min in the dark. The stained cells were acquired in BD FACSverse with appropriate gating.

Statistical analyses
All statistical analyses were performed using GraphPad Prism version 7.0 (GraphPad Software, La Jolla, CA, USA). Data were analyzed by ANOVA or two-tailed unpaired Student's t-test or two-way ANOVA.
TCGA-BRCA cohort data. We thank biorepository of Rajiv Gandhi Cancer Institute and Research Center (RGCIRC), Delhi. Graphical abstract Figure was drawn in BioRender.com.