Landscape analysis of alternative splicing in kidney renal clear cell carcinoma and their clinical significance

Landscape of AS events profiles in TCGA-KIRC cohort

Comprehensive AS events profilings of 533 KIRC patients were obtained from TCGA-KIRC cohort. All AS events were divided into seven types, which are ES, AP, AT, AA, AD, RI, and ME. Separately, the seven types of special splicing patterns were shown in Figure 1A. According to integrated analysis of SpliceSeq tool, we identified 46276 AS events from 10577 genes, comprising of 18062 ESs in 6899 genes, 9472 APs in 3793 genes, 8611 ATs in 3762 genes, 3810 AAs in 2676 genes, 3265 ADs in 2295 genes, 2821 RIs in 1897 genes, and 235 MEs in 227 genes (Figure 1B). As a result, we have known that each gene may account for more than one AS event. And the results showed that ES ranks first in all AS events, then comes the AP and AT.

Figure 1. Landscapes of AS events profiles in TCGA-KIRC dataset. (A) Schematic diagram of seven types of AS events, that is AA (Alternate Acceptor site), AD (Alternate Donor site), AP (Alternate Promoter), AT (Alternate Terminator), ES (Exon Skip), ME (Mutually Exclusive Exons) and RI (Retained Intron). (B) Overview of AS events in TCGA-KIRC dataset.

Prognostic-associated AS events in TCGA-KIRC cohort

AS events were related to carcinogenesis, progression and prognosis in a variety of tumors. Therefore, we conducted Cox regression analysis to identify prognostic-associated AS events of KIRC. As a result, 3823 ES events in 2504 genes, 3587 AP events in 1853 genes, 3048 AT events in 1546 genes, 1002 AA events in 866 genes, 887 AD events in 755 genes, 1186 RI events in 868 genes, and 75 ME events in 72 genes were associated with KIRC prognosis. ES was the most frequent AS event related to OS, followed by AP and AT. Among them, more than half were unfavorable predictors, and the others were favorable. The result was shown as an upset plot (Figure 2A). According to the results, we found that one single gene had more than one AS type. The ES was also the most common prognostic AS type. The top 20 prognostic AS events related genes of seven AS patterns were visualized in forest plots (Figure 2B–2H).

Figure 2. Prognosis-related AS events profiles. (A) Intersection UpSet plot of seven types of prognosis-related AS events in KIRC. (B–H) Forest plots of HRs of the top 20 prognostic-associated seven types of AS events. The color scale on the right side indicates the P-value.

Pathway and functional enrichment analysis

We utilized Metascape to identify the enriched pathways of prognostic associated AS to conduct a deep investigation of the underlying mechanism of prognostic associated AS genes in KIRC. The pathway and functional enrichment analysis indicated that metabolism related pathways were the most frequently involved, including “Metabolism of amino acids and derivatives”, “nucleobase-containing compound catabolic process”, “peptide biosynthetic process”, etc. (Figure 3A). Protein-protein interaction network analysis revealed that these AS genes were concentrated in twenty MCODE components (Figure 3B).

Figure 3. Pathway and functional enrichment analysis. (A) Top 20 pathways and functional enrichment clusters. (B) The 20 MCODE of genes from prognostic-related AS events.

Development of AS-based predictive model

Because of the high mortality of KIRC, we aimed to develop an AS-associated prediction model to improve prognostic prediction of KIRC. To construct the prediction model precisely, we performed LASSO Cox regression analysis to identify the most optimal prognostic AS model of seven types of AS patterns. As shown in Figure 4, we figured out the optical lambda (Figure 4A–4H) from the seven types of AS and established clinical prediction models. Then, according to the median risk score of prediction models, patients were divided into either low or high-risk groups. And those in the high-risk group had poorer survival compared with the ones in the low-risk group from the Kaplan–Meier survival analysis (Figure 5A–5H). By using time-dependent ROC analysis, we evaluated the predictive accuracy of AS signatures. The time-dependent ROC analysis results showed that integrated AS signatures, AD predictive signatures, and AT predictive signatures had better performance than other AS signatures (Figure 5I). All these results showed that these AS signatures could be a reliable and robust predictor of prognosis for KIRC patients.

Figure 4. LASSO Cox model to construct prognostic-related AS signatures. (A–H) Indicated constructions of the most valuable prognostic-related AS signatures and the LASSO coefficients profiles of seven types of AS events. The vertical lines were drawn at the optimal values by the minimum criteria and the 1-SE criteria.

Figure 5. Kaplan-Meier survival analysis of prognostic-related AS signature. (A–H) Survival analysis of seven types of prognostic-related AS signature. (I) Time-dependent ROC for different survival prediction models.

Network of prognostic-associated alternative splicing events

Given all AS alteration patterns were broadly regulated by critical splicing factors (SFs), thus, we identified the prognostic-related SFs and explored the correlations between prognostic-related SFs and AS events. Univariate Cox regression analysis showed that 30 SFs were related to the survival of KIRC patients significantly (Supplementary Table 1). Most of the SFs predicted good outcomes of patients, such as FMR1 and HNRNPU (Figure 6B, 6E). Furthermore, analysis of correlation revealed that a total of 1369 prognostic-associated AS events were correlated to 23 prognostic-associated SFs significantly in KIRC patients (Supplementary Table 2). The results of correlation analyses were visualized with a splicing regulatory network (Figure 6A). Interestingly, we found that most of the prognostic-related SFs were positively (red lines) related to the adverse prognostic AS events (green dots) and negatively (green lines) related to favorable prognostic AS events (red dots), Representative correlations between SFs and specific AS events were presented in dot plots. For example, expression of FMR1 was positively correlated with AT of RBM15 (Figure 6C), and negatively correlated with AT of EIF4E2 (Figure 6D), while expression of HNRNPU was positively correlated with AP of NUMB (Figure 6F), and negatively correlated with AT of SAMD4B (Figure 6G).

Figure 6. AS correlation network in KIRC. (A) Splicing correlation network in KIRC. The expression of all the prognostic-related splicing factors (blue dots) was negatively (green line) or positively (red line) associated with PSI value of all the favorable prognosis and adverse AS events (red dots and green dots, respectively). (B, E) Using the GEPIA tool to analyze the prognosis of splicing factors FMR1 and HNRNPU, respectively. (C, D, F, G) Representative dot plots of correlations between expression of SFs FMR1 or HNRNPU and PSI value of AS events (P < 0.001).

Favorable effect of FMR1 on KIRC

To further validate the good outcomes of SFs in KIRC, we conducted in vitro experiments with KIRC tissues and cell lines. The results from TCGA and CPTAC databases indicated that FMR1 is downregulated in KIRC tissues compared to normal tissues both at protein and gene levels, and the expression of FMR1 is negatively related to stage and grade of KIRC (Figure 7A). Consistently, our clinical samples also verified that FMR1 expression is reduced in KIRC tissues compared to paracancerous tissues (Figure 7B). So, we chose FMR1 for the next experiments. FMR1 was overexpressed by transfecting FMR1 plasmid, and was knockdown by transfecting FMR1-target-siRNAs separately in ACHN cell. The overexpression and knockdown efficiency were confirmed by qRT-PCR (Figure 7C) and WB (Figure 7D). After overexpression of FMR1, the cell proliferation was significantly suppressed (Figure 7E). Furthermore, clonogenic survival ability of ACHN was also inhibited by FMR1 (Figure 7F). On the contrary, FMR1 knockdown promoted ACHN cells proliferation and clonogenic formation efficiency (Figure 7E, 7F). Cell migration contributed to tumor progression. Then, we conducted transwell migration assay, which showed that FMR1 overexpression could reduce ACHN cell migration rate, while FMR1 knockdown promoted the migration, compared to corresponding control group (Figure 7G). Statistical analysis of 3 independent experiments confirmed the results (Figure 7H). Afterall, our results revealed the favorable effect of SFs FMR1 on KIRC, but the underlying genetic mechanisms still need to be further explored.

Figure 7. Protective effect of FMR1 in KIRC. (A) The expressing level of FMR1 in KIRC and normal tissues. (B) The qRT-PCR detects FMR1 expression of paired paracancerous tissues and KIRC cancer tissues from Zhongnan Hospital. (C) Overexpression and knockdown efficiency in ACHN cells. (D) WB confirmed the overexpression and knockdown of FMR1 in ACHN cells. (E) MTT assay to investigate the proliferation of ACHN cells. (F) Clonogenic formation results from three independent experiments. (G) Transwell migration assay to investigate cell migration ability, scale bar = 100 μm. (H) Confirmed by statistical analysis. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Data acquisition and curation process

RNA-seq raw counts of kidney renal clear cell carcinoma (KIRC) and corresponding clinical data were achieved from The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga/). Then, SpliceSeq tool, a Java application, was used to analyze the mRNA splicing patterns of KIRC patients. Rating from zero to one, the Percent Spliced In (PSI) value was employed to evaluate AS events and calculate for seven types of alternative splicing events: Exon Skip (ES), Alternate Promoter (AP), Alternate Terminator (AT), Alternate Acceptor site (AA), Alternate Donor site (AD), Retained Intron (RI) and Mutually Exclusive Exons (ME). AS events with PSI value >75% were included for all analysis.

Prognostic associated AS events identification

To identify all prognostic AS events in KIRC, 537 KIRC patients in total with overall survival (OS) data were collected. The corresponding clinical information was extracted and summarized. Then, the univariate Cox regression was taken to identify prognostic factors of seven types of AS patterns. We used a novel visualization tool UpSet plot [33] to quantitatively analyze the interactive sets, to visualize intersections of seven types of prognostic associated AS events.

Gene functional enrichment analysis

We conducted pathway enrichment analysis of genes of survival-associated AS to detect the underlying mechanisms of survival associated alternative splicing genes in KIRC. Metascape (http://metascape.org/) is a web portal for gene analysis and annotation. We used the tools in this portal to search for deep insight into the biological functions of survival-associated AS genes [34]. The automated meta-analysis tool provided in this portal could be used to detect unique and common pathways within a group of orthogonal target-discovery studies. And the protein-protein interaction (PPI) analysis was also accessible based on BioGrid, enrichment heatmaps generation and interactive visualization of Gene Ontology (GO) networks.

Establishment of predictive models by LASSO Cox regression

By conducting Cox regression model with LASSO (least absolute shrinkage and selection operator) penalty, the LASSO Cox regression analysis could simultaneously achieve shrinkage and variable selection. We used LASSO Cox regression model to establish the most predictive models of KIRC based on AS events. To prove if the predictive models can make a difference between the long OS patients and the shorter OS ones, Kaplan–Meier curves were performed. And the receiver operator characteristic curves were performed to further evaluate the efficiencies of each predictive model by running survivalROC package in R.

Construction of splicing factor regulatory network

All the human splicing factors (SFs) were achieved from the SpliceAid database (http://www.introni.it/splicing.html) [35]. EdgeR package in R software was performed for normalized expression profiles of SFs in TCGA-KIRC dataset [36]. The expression values of SFs were log2(* + 1)-transformed for further analysis. We then conducted univariate Cox regression analysis to determine SFs associated with survival. Next, Pearson correlation analysis was conducted between the PSI value of AS events, which were obtained from previous prognostic signature construction, and the expression level of prognosis-associated SFs. Then, with a criterion of adjusted p < 0.05, we selected the significant correlation pairs and used Cytoscape (version 3.5.1) to construct the potential SF-AS regulatory network.

Ethical statement of KIRC tissues

Ten paired KIRC tumor and paracancerous tissues were obtained through surgery in Zhongnan Hospital of Wuhan University, and patients’ clinicopathological characteristics are detailed in Supplementary Table 3. The fresh tumors and paracancerous tissues were immediately immersed in liquid nitrogen for subsequent experiments. The pathological diagnosis of KIRC tumor tissues specimens were independently validated by two pathologists. All patients had signed informed consents before the study. Using of human KIRC tissues was approved by the Ethics Committee of Zhongnan Hospital of Wuhan University (Approval No. 2020102).

Cell culture and transfections

Human KIRC cell lines ACHN was provided by Chinese Academy of Sciences, China. MEM supplemented with 10% fetal bovine serum was used to cultivate the cells. We purchased siRNA and plasmid from GenePharma (China). The sense sequences were as follows: FMR1-si1, 5′-GTGTTAGTGGCTTCATCAGTT-3′; FMR1-si2, 5′-GCCTGATAGGCAGATTCCATT-3′; Control-siFMR1, 5′-UUCUCCGAACGUGUCACGUTT-3′. Cell transfection was mediated by Lipofectamine 2000 with either plasmid or siRNA.

RNA extraction and qRT-PCR

Human KIRC tumor and paracancerous tissues, and ACHN cells were taken to extract cell total RNA by following the instruction manual of Qiagen RNeasy Kit (Cat. #74101). The RNA was quality controlled and reverse transcribed to cDNA. qRT-PCR was then conducted using iQ™ SYBR^®-Green Supermix (Bio-Rad, USA). The primer sequences were as follows: FMR1 forward primer, 5′-CCAACAAACCTGCCACAAAAG-3′, reverse primer, 5′-GCACACATTTGCCGTAAGTCTT-3′; GAPDH forward primer, 5′-TGCACCACCAACTGCTTAG-3′, reverse primer, 5′-GATGCAGGGATGATGTTC-3′.

Protein extraction and WB

ACHN cells were lysed in RIPA buffer for 30 minutes on ice. After spinning for 15 minutes at 4°C, the total protein included in the supernatant was collected and boiled in water bath for denaturation. Then, the protein band was separated and detected as we previously described [37].

Cell phenotype experiment

After transfected with plasmid or siRNA, 3000/well cells were seeded into 96-well plates. Cell absorbance at 490 nm was measured every single day after treated with MTT and DMSO to evaluate cell viability. Similarly, 1000/well cells were seeded into 6-well plates for 10–14 days to perform clonogenic formation assay. The clone number was counted after the colony get fixed and stained with crystal violet.

For transwell migration assay, 3 × 10⁴ cells were seeded into the upper chamber with serum-free medium. Cells migrated to the lower chamber containing serum medium from the upper chamber after 24 hours incubating. Then, the cells were fixed, stained and counted for statistical analysis.

Statistical analysis

R software 3.5.0 was used to conduct all the statistical analyses. Two-tailed Student’s t-test was chosen to investigate if there is a statistical difference between two groups. The cutoff probability value was set at P < 0.05. χ2 test was used to analyze correlations between AS events and clinicopathological parameters. As to the survival differences, we performed Kaplan-Meier survival analysis between the low-/high-risk groups. By running survival package in R, we conducted a two-sided log-rank test. And by running survivalROC package in R, the time-dependent receiver operating characteristic (ROC) analysis was conducted to detect the prediction accuracy of the predictive model.

Availability of data and materials

Raw data were deposited in The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga/). The data are available from the corresponding author upon reasonable request.