Pan-Cancer Profiling of Intron Retention and Its Clinical Significance in Diagnosis and Prognosis

1. Introduction

Alternative splicing (AS) is a way for eukaryotes to produce multiple transcripts from a single gene by changing the composition of exons during RNA processing[1]. Intron retention (IR) is a type of AS and affects about 80% of human genes[2]. IR transcripts tend to be degraded by nonsense-mediated decay (NMD) or exosome pathway. Therefore, IR coupled with RNA surveillance can regulate gene expression and affect various biological processes [3,4]. IR transcripts can be detained in the nucleus and wait for splicing signals, such as in the T cell activation process responding to external stimulus[5]. IR can also produce novel proteins with different functions[6,7,8] or subcellular locations[9], and play essential roles in many key biological conditions[3,10,11]. IR is under sophisticated regulation that can be affected by sequence variations, splicing factors, epigenetic and transcriptional regulatory mechanisms, etc[12].

Tumors generally have 30% more aberrant alternative splicing events than normal tissues [13] and can give rise to many tumor-specific transcripts associated with oncogenic functions and drug resistance[14,15,16,17]. There are numerous relevant studies on exon skipping in cancer, while only a handful of them have identified features and functions of IR. Dvinge and Bradley reported that compared to adjacent normal tissues, more introns were retained in tumor tissues for 15 cancer types, and many of them could be detected in the cytoplasm [18]. Evidence showed that IR can promote cancer development. Somatic mutations in tumors can trigger IR and these mutations are enriched in tumor suppressor genes (TSGs)[19]. Inactivation of histone H3K36 methyltransferase SETD2 reduced DVL2 IR, and as a result activated Wnt signaling to promote colon cancer predisposition[20]. Similarly, ZRSR2 loss increased retention of a minor intron of LZTR1 and resulted in enhanced RAS signaling, potentially driving leukemia[21]. IR produces epitopes presented on MHC I, making it a potential source of tumor-specific antigens (neoantigens)[22]. However, these researches have mainly focused on individual or patient-specific IR events, a systematic survey of recurrent IR alterations at a pan-cancer level, especially from the prognostic perspective, is still lacking.

In the present study, we profiled the IR landscape of 33 cancer types in the Cancer Genome Atlas (TCGA), identified IR events that were differentially regulated between normal and tumor samples, and discovered IR events associated with survival. Some informative introns not only show potential in accurate diagnosis and prognosis with machine learning methods but also are involved in cancer pathology validated in our experiments (Figure 1). Many of these introns were recurrently retained in multiple patients across cancers, serving as a resource of potential diagnostic and prognostic biomarkers, even as promising targets for new therapy.

2. Materials and methods

2.1. RNA-seq bam download and IR quantification

We downloaded RNA-seq bam files using the GDC data transfer tool, including 33 cancer types and over 10,000 tumor and adjacent normal samples in the TCGA project.

Stringtie (version 2.1.3)[23] was used to quantify gene expression. IRFinder (version 1.3.1)[2] was used to identify and quantify IR. And before applying IRFinder, samtools (version 1.9) was used to sort bam files by read pairs. The genome version was hg38 and the gene annotation version was gencode v35.

We performed quality control on both intron and sample levels. Filtering out introns that overlap with known exons (flagged by IRFinder as “known-exon”) resulted in 243,151 introns covering 21,520 genes genome-wide. Before quantifying IR, at least 4 junction reads (split reads) spanning flanking exons were required to make sure that the isoform was expressed. IRratio was used to measure the IR level, and it was calculated by dividing the median read depth of an intron by itself plus the number of reads spanning flanking exons. If the coverage of an intron was less than 20%, it was likely to be completely spliced and its IRratio was assigned to 0. If the coverage of an intron was above 70% and the median read depth was above 3, the intron was likely to be retained and its IRratio was kept, which must exceed 0. Otherwise, the retention state could not be accurately decided for an intron and its IRratio was set as missing. IRFinder automatically decided if an RNA-seq sample was suitable for IR detection based on the ratio of the number of reads that map to intergenic regions to the number of reads that map to coding regions. When the ratio exceeded 10%, it gave a warning message, and such samples were not used for further analysis. Over 90% of TCGA RNA-seq samples passed this QC.

The majority of TCGA RNA-seq was unstranded, so the read orientation could not be determined, which may lead to false positive IR detection[2,24]. Nevertheless, we made use of a small number of strand-specific RNA-seq samples in TCGA which passed sample quality control (n = 34) to generate a “whitelist” of retained introns, meaning they were more likely to be genuine and reliable. Specifically, we obtained a maximum IRratio of each intron in these samples, and introns with a maximum IRratio above 0.08 comprised the “whitelist” (n = 47,026). That is, these introns were retained in at least one stranded RNA-seq sample, which made them more reliable. Thus, differentially retained introns and survival-associated introns were restricted to only “whitelist” introns for higher confidence.

2.2. Differential IR and differential gene expression

Because IR levels did not follow a normal distribution across individuals, we used paired Wilcoxon rank-sum test to detect differential IR events (DIRs) in paired tumor and normal samples (n > 15) for 14 cancer types. A differentially retained intron should have a P value less than 0.05 and the difference of median IRratios between tumor and normal samples should be larger than 0.1.

In terms of differentially expressed genes (DEG) detection, DESeq2[25] was used, and we selected DEGs with Padjusted < 0.05 and |log2FC| > 1.

2.3. Dimensionality reduction and visualization

If the quality control mentioned above was not passed, the IRratio was set to a missing value. We generated IRratio matrices, kept IR events with missing rates less than 30%, and imputed missing values with a mean value from the remaining samples. Principle component analysis was then performed. We extracted the first 100 principal components and further analyzed them with package Rtsne[26] to draw tSNE plots, with perplexity set to 30.

2.4. Functional enrichment

We used the package clusterprofiler[27] to enrich gene ontology (GO) biological process terms and KEGG pathways for target gene sets. GO terms related to cancer hallmarks were retrieved from two previous studies[28,29].

2.5. Sequence features analysis

We classified introns into four types to compare their sequence features. Constitutive introns (n = 48,344) were introns with median IRratio of 0 in tumor and normal samples across all cancers; not regulated IR were introns retained in tumor or normal samples in any cancer but showed no significant difference between tumor and normal samples (n = 10,887); down DIRs were introns with reduced retention level in tumor samples in any cancer (n = 401); up DIRs were introns with increased retention level in tumor samples in any cancer (n = 669).

Conservation analysis: hg38 version of phastCons30way.bw was downloaded from UCSC, and bwtool[30] was used to calculate a mean phastCons[31] score around boundaries of the above 4 types of introns separately.

Splice score of 5′ and 3′ sites were calculated using the maximum entropy modeling method[32].

Intron GC content, length and relative gene position were all analyzed based on hg38 and gencode v35. To analyze the distribution of introns in genes, the relative gene position of an intron was calculated by dividing the rank of the intron (5′ to 3′ orientation) by the total intron count of this transcript.

We adjusted a nonsense-mediated decay (NMD) prediction rule from a previous study[33]. Specifically, when an intron was retained in a protein coding gene, it was very possible to introduce a premature stop codon (PTC) and elicits NMD unless under following conditions. If the intron resides in a 5′ or 3′ untranslated region (UTR), or the intron was close to start codon (< 200 nt), or it was the last intron, NMD would not be elicited. If no PTC was produced, or the PTC was located within 55 nt upstream of the last exon-exon junction, NMD would not be elicited either. Otherwise, the IR transcript was prone to be targeted by NMD.

2.6. Random Forests model

We merged DIRs from different cancer types, and filtered out the ones with pan-cancer missing rates over 30% and the ones that were inconsistently up or down-regulated in different cancers. This resulted in 273 DIRs to be used in Random Forests models for pan-cancer modeling. R package randomForest[34] was applied to the train model, with mtree = 500, mtry = 3 and proximity = TRUE. A hundred times four-fold cross validation was used to calculate a pooled area under the curve (AUC) for the training set (paired tumor and normal samples from 14 cancers), and the receiver operating curve (ROC) in randomly selected one run was drawn using a R package pROC[35]. We used the Rfcv function to predict model performances with a sequentially reduced number of DIRs (ranked by importance), with cv.fold = 5 and step = 0.9.

2.7. Survival analysis

For introns that had valid IRratios in at least 50% of patients and at least 5% of valid IRratios were above 0.1, we restricted our analysis to patients with IRratios over 0 and classified them into IR-high or low groups based on median IRratio. Then we performed Univariate Cox regression and selected intron related to overall survival (OS) or disease-free survival (DFS), and an unadjusted p value less than 0.05 was considered significant. Similarly, gene expression associated with survival was identified by selecting genes with a median expression level over 1 TPM, dividing patients into high and low expression groups based on the median cutoff, and performing univariate Cox regression.

2.8. LASSO regression to build a prognostic model

To build an IR-based prognostic model for each cancer type, we selected introns with missing rates less than 20% (missing values were later filled with mean) and perform LASSO regression with R package glmnet[36,37]. Candidate introns and corresponding coefficients were derived with the λ parameter associated with the minimum mean error or with one standard error. Intron retention risk (IRR) score was calculated as the sum of IRratios multiplied by corresponding coefficients of candidate introns. And then patients were divided into IRR-low and high groups based on median value.

2.9. Cell culture and lentiviral transfection

HEK293 (human embryonic kidney 293 cells), H1299 cells, A549 cells were purchased from National Collection of Authenticated Cell Cultures and cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) (Invitrogen, 11960044) supplemented with 10% FBS (Gibco), streptomycin (100 μg/ml), and penicillin (100 U/ml) at 37℃/5% CO2.

We applied lentivirus transfection-mediated gene-silencing strategy to stably knockdown target gene STN1. HEK293T cells were transfected as described[38]. Specifically, shRNAs (shRNA sequences were listed in Supplementary Table S7) were obtained from Sigma-Aldrich, then annealed and cloned into pLKO.1 vector. HEK293 cells grown in 6-well plates were transfected with 1 μg constructed vectors or empty control vectors pLKO.1 with VSVG and gag/pol encoding plasmids by using Lipofectamine 2000 (Invitrogen). Then, the virus supernatant was harvested in 24 hours to infect A549 cells in six-well plates seeded one day ahead. After incubation for one more day, A549 cells were screened by 2.5 μg/ml puromycin for 24 hours. The surviving cells were cultured for two more days and then harvested for following RNA extraction and cell proliferation assays.

2.10. RNA preparation, RT-PCR and qRT-PCR

Total RNA was extracted from cells with TRIzol (Invitrogen) according to the manufacturer’s instructions (Invitrogen) and reversely transcribed into cDNA with oligo-(dT) primer (Supplementary Table S7) using FastQuant RT Kit (Tiangen).

For RT-PCR, 100 ng cDNA was used for each PCR reaction, and PCR products were detected by agarose gel electrophoresis. Gene expression at the RNA level was quantified by qRT-PCR using a 2× SYBR mix (Vazyme). GAPDH served as an internal control. Then, the reaction was run on the Bio-Rad CFX manager machine. The IR and spliced transcripts of STN1 were quantified by qRT-PCR using primers specifically targeting the retained intron and exon junction, respectively. All primer sequences were listed in Supplementary Table S7.

2.11. RNA stability assay, and isolation of nuclear and cytoplasmic fractions

A549s cells were treated with 10µg/ml actinomycin-D (Act D; Sigma-Aldrich, Inc., St. Louis, MO, USA, A4262) for 0, 2, 4, and 6 hours, respectively. The total RNAs were extracted, and RNA levels were quantified by qRT-PCR as described above.

For nuclear and cytoplasmic fractionation, A549 cells were cultured in a T25 flask until 90% confluence. Then, cells were trypsinized, washed twice with cold PBS, and then cells were fractioned by Nuclear/Cytosol Fractionation Kit (Phygene, PH1466). Following RNA extraction and qRT-PCR were carried out as described above.

2.12. Psi-CHECK2 constructs and dual luciferase assay

Luciferase reporter gene expression vectors were bought from Promega and were prepared according to the manufacturer’s protocol. Briefly, the 590bp intron retained 5′ UTR and 146 spliced 5′ UTR was PCR amplified using STN1_5′ UTR_Forward (Nhe I) 5’-TACGACTCACTATAGgctagcggggtcgtcgccgccag -3’ and STN1_5′ UTR_Reverse (Nhe I) 5’-TTGGAAGCCATGGTGgctagccaggctgcatcaagaggca-3’. PCR fragments were cloned into the NheI-restricted site of the psi-CHECK2 vector. The STN1 5′ UTR was inserted directly upstream of the Renilla gene. The psi-CHECK2 construct containing the mutated 5′ UTR IR fragment was synthesized by Tsingke company, where the three start codons within the intron (181ATG, 283ATG, and 392ATG) were all mutated to AGC. DH5α cells were transformed with the three distinct constructs and cultured on ampicillin-containing media.

One day before transfection, 1 × 10⁵ HEK293T cells were seeded into each well of a 24-well culture plate in 500μl DMEM supplemented with 10% FBS. Cells at 70% confluency were transfected with three psi-CHECK2 constructs with Lipofectamine 2000 reagent. 24 hours after transfection, growth media were removed and cells were washed gently with PBS. Passive lysis buffer (Promega) (100ul/well) was added and gently shaken for 15 min at room temperature, then cell lysates were harvested for dual luciferase assay. The activities of firefly and Renilla luciferase were measured using the Dual-Luciferase® Reporter 1000 Assay System (Promega) according to the manufacturer’s instructions. A total of 25μl cell lysate was transferred into a white opaque 96-well plate. The luminescence obtained for the mutated and wild-type constructs was normalized with the internal control firefly luciferase signal. Each experiment was performed in triplicate, and three independent experiments were performed. Quantitation of the reporter gene assay was calculated as mean ± SEM. Student’s t test was used to determine significant differences between each mutated construct compared with the wild-type construct.

2.13. Colony formation assay, transwell migration assay and cell proliferation assay

For the colony formation assay, cells were cultured in the six-well plate at a density of 800 cells per well. Cells were cultured under normal culture conditions for 15 days. For fixation of the cell, after the medium supernatant was removed, the cells were treated with 4% paraformaldehyde and were stained with 1% crystal violet (SigmaAldrich) for 15 min. Then, the plates were washed with phosphate-buffered saline (PBS) and photographed.

For the transwell migration assay, we used a 24-well transwell chamber (Corning). Cells were suspended in non-serum DMEM and then seeded in the top chamber of the transwell with a density of 1 × 104 per chamber, and 300 μl fresh complete DMEM (10% FBS) was added to the bottom chamber. After incubating for 48 hours, use PBS to wash the cells in the top chamber twice, and then fix the cells with 4% paraformaldehyde for 15 min, followed by staining with 1% crystal violet (Sigma-Aldrich) for 30 min. After washing and wiping off the cells in the inner side of the top chamber, the migratory cells adhering to the bottom surface of the membrane were observed, photographed and then counted by ImageJ software.

Cell proliferation assay was performed on cultured cells at four time points (24, 48, 72, and 96 h, respectively). A total of 100 μL cells were seeded in a 96-well plate with at least 1000 cells per well and four replicates for each time point. Cell Counting Kit-8 (CCK-8) reagent (Dojindo) was added according to the manufacturer’s protocol. Then, cells were incubated at 37°C for 2 hours and the absorbance at 450 nm was measured with a microplate reader (TECAN).

3. Results

3.1. Landscape of intron retention in 33 cancer types

RNA-seq data (in bam format) of the primary tumor and adjacent normal tissues from 33 TCGA cancer types, were downloaded and processed with IRFinder[2]. A total of 10,189 samples passed quality control and were used in this study (Figure 2A, Supplementary Table S1). IRratio was used to measure the retention level of each intron. Median IRratios of each intron among tumor samples were used to represent the IR level in each cancer type, and the correlation coefficients between different cancers were generally higher than 0.6. While acute myeloid leukemia (LAML), esophageal carcinoma (ESCA), and stomach adenocarcinoma (STAD) had lower similarity IR patterns with the other cancer types (Figure 2B).

Next, we computed the average number of retained introns (IRratio > 0.1) for each sample (Figure 2C, Supplementary Table S1). Notably, LAML tumor samples had an average of 27,148 retained introns, which was an-order-of-magnitude more than many other cancers. The sequencing depth and read length of LAML were comparable to that of other cancers (Supplementary Table S1), and no significant degradation of RNA was observed (Supplementary Figure S1A and B). Manual inspection showed that often multiple introns were retained within a gene and these signals did not appear to come from genomic contamination (Supplementary Figure S1C-E). In line with this, Dvinge and Bradley have also observed a much stronger IR increase in LAML than in other cancers [18]. In total, 13 cancer types exhibited significantly increased numbers of retained introns comparing to normal samples, while breast invasive carcinoma (BRCA) showed the opposite trend. Our results validate the findings of Dvinge and Bradley, and again suggest that inefficient intron removal is a common phenomenon in many cancers except breast cancer.

3.2. Differentially retained introns between tumor and normal tissues

To accurately compare IR between normal and tumor samples, we analyzed 14 cancer types each has at least 15 paired normal and tumor samples (Supplementary Figure S2A). Principle component analysis (PCA) showed varying degrees of difference in global IR patterns for normal and tumor samples across cancer types (Supplementary Figure S2B). We used paired Wilcoxon rank-sum test to detect introns significantly up- or downregulated in tumor tissues compared to matched adjacent normal tissues for each cancer type (P < 0.05) and required a difference of median values between normal and tumor samples greater than 10%.

A total of 988 differential introns were detected across 14 cancers, most of which were retained in both normal and tumor samples but with altered retention levels (Supplementary Figure S3A), while a handful of them were mainly retained in tumor samples and spliced in normal samples (cancer associated IR), or the other way around (i.e. normal associated IR) (Figure 3A). Colon Adenocarcinoma (COAD) and BRCA were detected with the most differential IR events (DIRs), but most DIRs were upregulated in COAD and downregulated in BRCA. DIRs were often specific to one type of cancer (Supplementary Figure S3B), and DIR genes usually had only one intron that was differentially retained (Supplementary Figure S3C). IR alteration between tumor and normal samples demonstrated a weak negative correlation with mRNA expression change (R = -0.17, P = 2.7e-11, Spearman correlation), and on average less than 15% of DIR genes were differentially expressed according to IR change (Figure 3B). Their mild correlation reflects the known function for IR to finetune gene expression through RNA degradation[3,5,10,11,39,40,41].

Mutation-induced IR has been widely reported to inactivate tumor suppressor genes (TSGs)[19], but such mutations were found in only a minority of patients. However, DIRs affected multiple patients. Although not enriched in COSMIC TSGs or oncogenes[42], DIRs were found in some cancer-related genes. For example, retention of intron 15 of LZTR1 was upregulated in COAD, STAD, and uterine corpus endometrial carcinoma (UCEC) tumor samples (Supplementary Figure S4A). LZTR1 is a TSG that negatively regulates RAS signaling through ubiquitination[43,44,45], and its minor intron retention is associated with tumorigenesis in leukemias[21]. Since the expression of LZTR1 did not significantly change between tumor and normal samples, elevated IR, implying a decreased level of normally spliced functional products. Another TSG that was affected by DIR was ERCC4 (Supplementary Figure S4B), which functions in nucleotide excision repair[46,47]. DIR was also found in oncogenes. Tumor samples of BRCA exhibited lower retention of intron 3 in CSF3R (Supplementary Figure S5), which is a highly mutated oncogene in chronic myeloid leukemia[48,49]. DIRs in different cancers could affect various biological pathways, such as DNA damage and cell cycle checkpoint in lung squamous cell carcinoma (LUSC) (Supplementary Figure S3C). These results indicate that differentially retained introns detected in multiple patients were found in genes with various biological functions including those relevant to carcinogenesis, though their functional consequences and detailed mechanisms require further study.

3.3. Sequence features of differentially retained introns

Many sequence features that are known to predispose introns to retention were also discovered in our retained introns, including higher GC content, shorter intron lengths, and biased distribution at the 3′ end of the gene, and lower likelihood of inducing NMD (Figure 3C,D; Supplementary Figure S6A). Notably, introns differentially retained in cancers were shorter compared to not regulated ones (the median intron lengths of the down-regulated, up-regulated, and unregulated IR were 456, 452, and 934 bp, respectively) (Figure 3C). On the other hand, no significant difference was observed in splice signal strength between different groups of introns (Figure 3D). This may be puzzling because many previous studies have found that retained introns usually have weaker splice sites[3,18,50]. However, Zhang and colleagues have also reported that introns differentially retained between prostate cancer and normal tissues have splice signals comparable with control introns[51]. Retained and constitutive introns demonstrated almost the same conservation within the intronic boundaries sequence, again supporting their comparable splice signals strength. However, constitutive introns had the most conserved flanking exonic boundaries, followed by down and up DIRs, indicating their higher selective pressure than unregulated retained introns (Supplementary Figure S6B).

For some cancers, up- and down-DIRs had different sequence features which may underlie different behaviors of these introns (increased or decreased retention levels) (Supplementary Figure S6). For instance, COAD and liver hepatocellular carcinoma (LIHC) had weaker 3′ splice signals in up-DIRs than down-DIRs. Differences in GC content or intron length were also observed in several cancers, such as lung adenocarcinoma (LUAD) and thyroid carcinoma (THCA).

3.4. Differential IR events have diagnostic potential

We used tSNE[26] to visualize IR patterns in over 10,000 samples from 33 TCGA cancer types. Genome-wide introns (n = 37,845) had limited tissue origin specificity and tumor vs. normal specificity (Figure 3E,F), whereas DIRs (n = 375) exhibited stronger specificity (Figure 3G,H). Similar results were seen when we restricted to 14 cancer types where DIRs were detected (Supplementary Figure S8). This suggests that DIRs may be able to distinguish tumor and normal samples, and then Random Forests[34] were applied to test IR’s potential as cancer biomarkers. After removing DIRs that were inconsistently up and downregulated in different cancer types and DIRs that had a missing rate exceeding 30% in all samples, 273 DIRs were left for model training. All TCGA samples were divided into 3 groups: 1) a training set including 1,210 paired tumor and normal samples from 14 cancers, which were the same samples used for the DIRs detection; 2) the first validation set including 5,738 unpaired tumor and normal samples from the same 14 cancers; 3) the second validation set including 3,190 samples from the rest 19 cancer types (Figure 4A). A hundred times 4-fold cross validation with the training set yielded a pooled area under curve (AUC) of 0.958, and the result of one randomly selected run was shown in Figure 4B. When the whole training set was used to train the model, the AUCs for the first and second validation sets were 0.936 and 0.838 respectively (Figure 4C,D). In addition, comparable performance can be achieved with the top 30 DIRs (Figure 4E,F). Two representative DIRs were displayed here, which were down-regulated (ZDHHC8) and up-regulated (ZNF800) in 5 cancer types respectively (Figure 4G,H). We also test the ability of DIRs as biomarkers for each cancer type, and the AUC of Random Forests models was all above 0.9 (Supplementary Figure S9). Taken together, DIRs demonstrated powerful diagnosis potential in a pan-cancer level and for a specific cancer type.

3.5. Identify prognostic IR events across cancers

We expanded our analysis to all primary tumor samples from 33 cancers in TCGA, to explore the prognosis potential of introns with varying degrees of retention in patients. Specifically, we investigated introns that had a valid IRratio for over 50% of patients within a cancer type, and the IRratio should be larger than 0.1 in at least 5% of these patients. And then univariate Cox regression was applied to patients with IRratio > 0. The results showed that 29 cancers had over 100 IR events associated with overall survival or disease-free survival (P < 0.05), of which LAML has the highest number of prognostic IR events (Figure 5A). Interestingly, prognostic introns of some cancers were predominantly associated with favorable or adverse outcomes, such as LAML and prostate adenocarcinoma (PRAD). A total of 11,522 prognostic IR events were detected among 33 cancers, of which 67.4% were specific to one cancer type, and 32.6% were shared by different cancers (Figure 5B, Supplementary Table S4).

Noteworthy, the number of prognostic IR events was not proportional to the number of genes whose expression was associated with survival (Supplementary Table S3). Moreover, only ~19% of prognostic IR genes had expression levels also associated with prognosis (Supplementary Figure S9). The discordance suggests that IR has unique prognostic features independent of gene expression. On average 42.8% of prognostic introns were negatively correlated with RNA expression across all cancer types, and 6.5% strong (r ≤ -0.5) negative association (Supplementary Figure S9).

3.6. Prognostic introns affect genes involved in tumorigenesis

Gene Ontology (GO) terms enriched for prognostic IR genes from 33 cancer types include chromatin organization, histone modification, cell cycle arrest, and GTPase activity regulation (Figure 5C). From the perspective of individual cancer types, there were 25 cancers with prognostic IR genes enriched for GO terms associated with cancer hallmarks (Figure 5D, Supplementary Table S5). Sustained angiogenesis, tissue invasion and metastasis, and genome instability and mutation were the three most frequently altered hallmarks. Prognostic IR genes in many cancer types were also enriched in KEGG pathways that are cancer related, such as the MAPK signaling pathway and PD-1 checkpoint pathway (Supplementary Figure S11). Compared with non-prognostic IR genes, prognostic IR genes in lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), kidney renal papillary cell carcinoma (KIRP) and LAML were significantly enriched for COSMIC cancer genes (Supplementary Figure S12). The increased last intron retention of MYC, a famous oncogene, was associated with a favorable prognosis in BLCA and negatively correlated with expression (Supplementary Figure S13). In a less studied gene, STN1, while the mRNA level did not associate with survival in LUAD (Supplementary Figure S14A), intron retention in its 5′ UTR demonstrated a significant connection with improved survival (Figure 6A). Moreover, this intron was retained at a lower level in tumor samples (Supplementary Figure S14C). To investigate the function of this 5′ UTR intron retention, we carried out a series of experiments as below.

A high IR level in the STN1 5′ UTR was also observed in the A549 lung adenocarcinoma cell line, with a spliced transcript to IR transcript ratio of approximately 2:1 (Figure 6B,C). The IR transcripts had a similar degradation rate as the spliced transcripts, which was consistent with our observation in TCGA LUAD samples that IR levels did not correlate with mRNA levels (Figure 6D and Supplementary Figure S14B). Although IR transcripts were more frequently distributed in the nucleus than spliced transcripts, more than 50% of IR transcripts were also localized in the cytoplasm (Figure 6E). By closely examining the intronic sequence of STN1 5′ UTR, we found three potential upstream open reading frames (uORFs) as predicted by ORFinder of NCBI (Supplementary Figure S15). Because uORFs are known to reduce protein expression[52], we speculated that this IR might diminish the translation efficiency of STN1, thereby reducing its protein levels. This hypothesis was confirmed by dual luciferase assay, which demonstrated that the IR in STN1 5′ UTR significantly impeded the translation of the main ORF, and mutating all three start codons within the intron to AGC nearly restored the translation efficiency (Figure 6F and G).

Given that STN1 is a telomere replication-related gene and that cells rely on telomerase for maintaining telomere levels[53], we postulated that the cell viability may be affected by STN1 levels, which is readily influenced by its 5′ UTR IR. To validate this, we knocked down STN1 to simulate a high level of 5′ UTR IR. STN1 knockdown significantly reduced the level of proliferation markers, such as CDK1 and MKI67 in A549 cells (Figure 6H), ultimately leading to impaired colony formation, deficient cell migration, and cell proliferation (Figure 6I and J). These results provide evidence for the regulatory role of intron retention of the STN1 5′ UTR in lung adenocarcinoma cells. Through the production of uORFs, as well as more detailed transcripts in the nucleus, the 5′ UTR IR reduced translated protein levels which may ultimately affect tumor growth and progression.

3.7. IR enables accurate risk stratification in multiple cancers

The current standard of care assessment of many cancers relies on clinical biomarkers and the detection of specific genetic mutations or chromosomal structural variation. Take acute myeloid leukemia, for example, cytogenetic screening and target gene sequencing-based stratification is effective for half of LAML patients. The rest of the patients are difficult to assign to a defined risk group since they are cytogenetically normal and do not carry mutations in known risk genes[54]. As a result, revising prognostic models using additional molecular features to improve stratification is an ongoing effort for LAML as well as other cancers.

Because LAML has the highest number of prognostic introns, we first tried to build an IR prognostic model with the least absolute shrinkage and selection operator (LASSO) regression[36] for LAML. The resultant model was made up of 4 IR signatures and was designated as IR risk (IRR) (Figure 7A, Supplementary Figure S16). Based on the IRR score, we can divide LAML patients from TCGA into two risk groups with significantly diverged overall survival outcomes (HR = 5.24, P = 1.36e-10, Figure 7B). We used univariate Cox regression to compare the contribution of IRR and other reported clinical and molecular predictors[54,55], and IRR turned out to be the most powerful predictor (Figure 7C, Supplementary Table S6). Although age over 60 at diagnosis is a known risk for LAML, IRR still enabled effective stratification in patients under or over 60 years old (Figure 7D). Similarly, we were able to construct a prognostic model based on 4~30 IR events with LASSO regression for 32 cancer types, except for testicular germ cell tumors (TGCT), achieving good risk stratification (Figure 8). This demonstrates the potential of IR as powerful prognostic markers for many cancers.

4. Discussion

Aberrant AS is one of the hallmarks of cancer, and tumor tissues commonly have about 30% more AS events than normal tissues[13]. By producing tumor-specific proteins or by altering the production of normal proteins, aberrant AS can lead to the activation of proto-oncogenes or inactivation of TSGs, ultimately affecting cell growth and differentiation, angiogenesis, tissue invasion, and metastasis[56]. Therefore, the study of aberrant splicing not only helps to understand the mechanisms of cancer initiation and development but also has potential clinical implications[57]. There is an intriguing imbalanced IR pattern in multiple cancers where tumors exhibit a significant increase in IR compared to normal samples, which is not observed in other AS types[18]. IR may contribute to cancer development by inactivating TSGs in cancer patients[19]. Moreover, IR can encode novel proteins and may be an important source of tumor-specific antigens[22]. In this study, we systematically quantified and profiled IRs in 33 cancers from TCGA and tentatively explored their clinical relevance. So far as we know, this is a comprehensive analysis of IR regarding the largest number of cancer types.

We identified differential intron retention events (DIRs) by comparing paired tumor and normal tissues in multiple cancer types. We found that the splicing signals of differentially retained introns were almost as strong as constitutive introns, which was consistent with Zhang et al.’s finding[51]. One possible explanation is that DIRs were recurrent in tumor and/or normal samples, and they may have similar biological importance as constitutive introns. Interestingly, DIRs were shorter than both constitutive introns and unregulated retained introns. Zhang et al. have reported that shorter exons are more likely to be excluded in cancers, and are possibly regulated by elevated transcription and dysregulation of some RBPs in cancer cells[58]. In addition, most genes are spliced co-transcriptionally[59], and increased RNA Pol II accumulation in retained introns has also been reported[3,60]. Whether shorter introns are more sensitive to transcription and other splicing alterations in cancer deserves further study.

Compared with adjacent normal tissues, dozens to hundreds of introns per cancer showed up or down-regulation in tumor tissues, resulting in a total of 988 DIRs across 14 cancer types. Some of them were cancer-type specific (such as CSF3R), while some others (such as LZTR1 and ERCC4) showed consistent alterations across multiple cancers. We further identified 30 DIRs that stratified tumor samples and normal samples, demonstrating their diagnostic potential. A recent study exerted intron splicing events generated by SF3B1 mutations that are specific to tumor patients, to design synthetic introns and achieve targeted clearance of tumor cells [61]. Similarly, DIRs in our study may also serve as promising candidates for therapeutic synthetic introns, since they were widespread in multiple patients and even in multiple cancers (e.g., IR of ZDHHC8 in Figure 4G). On the other hand, because retained introns can also encode peptides located on cancer cell surfaces[22], commonly retained introns may be potential targets for “off-the-shelf” immunotherapy.

We discovered several introns with the potential to be survival indicators. For example, we experimentally validated a functional prognostic intron in LUAD. Specifically, the 5′ UTR intron regulates the translation efficiency of STN1, which is essential for cancer cell proliferation through maintaining telomere replication and genome stability (Supplementary Figure S17). These insightful results also indicate potential novel therapeutic strategies to combat LUAD.

There have been ongoing efforts in exploring new biomarkers for risk stratification in cancers, such as gene expression[55,62,63,64,65]. However, splicing has been reported to outperform gene expression analysis in predicting survival in multiple tumor types[66,67]. The potential of alternative splicing (AS) in prognosis has also been demonstrated in various cancers, including ovarian cancer[68], colorectal cancer[69], lung cancer[70], esophageal cancer[71], liver cancer[72] and adrenocortical carcinoma[73]. We explored the prognostic power of IR in 33 cancer types, and most of the cancers (n=32) can be stratified with less than 30 introns. The performance of the IR-based prognostic model in TCGA LAML cohorts outperformed clinical and other molecular predictors. IR has been reported to indicate prostate cancer aggressiveness[51] and pancreatic cancer clinical outcomes[66,74]. Our results further suggest that IR can serve as an accurate and powerful biomarker in multiple cancers. Moreover, IR (including other AS types) is usually quantified in a relative rather than an absolute method, which may be less influenced by different processing procedures.

DIRs and prognostic IR events were two types of informative IR characterized in this study. Differentially retained introns across many cancer types were heterogenous and influenced genes of various functional categories. Prognostic introns affected genes involved in cancer-related pathways, including DNA damage and cell cycle regulation, angiogenesis, cancer cell invasion and metastasis. In DLBC, KIRP and LAML, prognostic IR genes were significantly enriched for COSMIC cancer genes. IR has been recognized as a widespread mechanism of TSG inactivation[19], and we found that IR may also regulate the activity of oncogenes, such as MYC (Supplementary Figure S13). Of note, one of the limitations of our study is that we still lack experimental assessments of the IR functions, since IR transcripts have diverse fates[12]. The second limitation is that only a few patients in TCGA dataset have matched normal tissues. More normal samples as background or negative controls will add to finding more informative IR events in future analysis.

6. Conclusions

Overall, we systematically characterized IR at a pan-cancer level, explored its clinical relevance and provided experimental evidence of IR in cancer pathology. Our results revealed a rich resource for IR that had potential clinical applications including cancer biomarkers and even therapeutic targets.