Systematic exploration of prognostic alternative splicing events related to immune microenvironment of Clear Cell Renal Cell Carcinoma

Clear cell renal cell carcinoma (ccRCC) accounts for almost 75% of all renal cancers, with high heterogeneity and poor prognosis. There is increasing evidence that alternative splicing (AS) is associated with tumorigenesis and immune microenvironment. However, studies on the relationship between AS events and immunity in ccRCC are still few but needed. The transcriptional data and clinicopathological information of ccRCC patients were downloaded from The Cancer Genome Atlas (TCGA) database. We grouped patients according to the ESTIMATE algorithm and identied for differentially expressed AS events (DEASs). The functional enrichment analysis and unsupervised consensus analysis were performed to investigate the association between AS events and immune features. The regulatory network of survival-related AS events and intersection SFs was used to explore potential mechanisms. Correlation


Introduction
Renal cell carcinoma (RCC) is one of the most common urological neoplasms, whose annual morbidity and mortality are increasing worldwide [1]. It is estimated that more than 403,000 people worldwide are initially diagnosed with RCC and 175,000 patients will die of the disease every year [2]. Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype of RCC, accounting for almost 75% of all renal cancers, and is the cause of most cancer-related deaths [3,4]. Due to the occult onset of RCC, many patients have been diagnosed in advanced stage, missed the opportunity of operation [4]. With the exception of the widely approved drug -sunitinib, limited treatment options for advanced RCC results in a poor prognosis, with a 5-year survival rate of less than 10% and a median survival of about one year [4][5][6][7]. New research showed that the emergence of immunotherapy, such as anti-CTLA-4 and anti-PD-1 antibodies, held great promise for the treatment of advanced kidney cancer, changing the traditional treatment pattern [8,9]. However, only a small percentage of patients show sustained responses[8, 10,11], leading us to wonder whether there may be some inhibitory mechanism between tumor microenvironments that reverses the effects of immunotherapy and helps tumors evade immune surveillance.
Tumor microenvironment is a complex environment composed of tumor cells, immune cells, stromal cells and other components, as a hotbed for tumorigenesis, which is closely related to tumor growth and evolution [12]. Previous studies had shown that cytotoxic CD8+ T cells and CD4+ helper T cells can target antigenic tumor cells to prevent tumor growth and prolong the outcome of many tumors, including RCC [13][14][15]. On the contrary, some studies had shown that the survival of RCC patients with high intracellular CD8+ T cells was worse than that of RCC patients with low intracellular CD8+ T cells [16]. Consistent with this, the recruitment of CD4+ T cells in tumor microenvironment (TME) could promote RCC proliferation by regulating TGFβ1/YBX1/-HIF2α signals [17]. Besides, Li et al. commented that TME can exert metabolic stress on immune cell in ltration, resulting in local immunosuppression and limited anti-tumor immune recovery [18]. Considering the inherent heterogeneity and high antigenicity of the tumor microenvironment of ccRCC, we believe that further exploring the interaction between genes and cells of TME is very necessary to help break through the bottleneck of immunotherapy of renal cell carcinoma.
The alternative splicing that occurs during RNA splicing procedure, regulated by splicing factors (SFs), takes many different patterns, such as Exon Skip (ES), Mutually Exclusive Exons (ME), Retained Intron (RI), Alternate Promoter (AP), Alternate Terminator (AT), Alternate Donor site (AD), and Alternate Acceptor site (AA), leading to biodiversity at both the RNA and protein levels [19][20][21]. Remarkably, about 95% of exon genes in human undergo alternative splicing, and uncontrolled splicing, such as abnormal expression and / or mutation of SFs, often leads to a series of diseases and even cancer [22][23][24][25]. Previous studies had shown that AS events were associated with tumorgenesis, metastasis, and poor prognosis of ccRCC, and its splicing regulatory network had hidden the potential genetic mechanism of ccRCC, which proved the great value of AS in developing novel anti-cancer strategies for RCC [26][27][28][29]. In addition, some research had proved that tumor-speci c splice may produce a new class of tumor-speci c neoantigens with strong immunogenicity and induce the production of anti-tumor immunity [30]. while, it may also cause changes in the biological behavior of tumor cells due to functional dysfunction to evade the immune surveillance of the host [31]. However, studies on the relationship between aberrant splicing events and immunity in ccRCC are still few but needed.
In the present study, the transcriptional data and relevant clinicopathological information of patients with ccRCC from The Cancer Genome Atlas (TCGA) were integrated. We applied the ESTIMATE algorithm to calculate the immune and stromal score of ccRCC patients and identi ed the differentially expressed AS events and differentially expressed SFs. Then the association between AS events and survival status, immune features, biological process, clinicopathological variables and SFs were further investigated.
Collectively, our study enriched the research content of the relationship between AS events and immune in ltration spectrum of ccRCC, which provided a theoretical basis for the development of new targets for immunotherapy and clinical diagnosis and treatment of ccRCC.

Data collection
The transcriptional data of 515 ccRCC patients were downloaded from TCGA data portal (https://tcgadata.nci.nih.gov/tcga/). Meanwhile, we obtained the data of AS events from the TCGA SpliceSeq database. There has been a wide consensus that the goal of using Percent Spliced In (PSI) ranging from 0-1 is to quantify events [32]. Then a strict set of screening conditions (sample percentage with PSI value of 75 and average PSI value of 0.05) was set to screen out suitable AS events included in consequent analyses. Besides, a list of 52 human SFs was extracted from the SpliceAid 2 database [33]. What's more, the relevant clinicopathological information were also collected from TCGA database, in which the OS of all patients were greater than 30 days.
Calculation and survival difference of immune and stromal scores ESTIMATE algorithm was carried to calculate immune and stromal scores for the mRNA expression data [34]. According to the results of ESTIMATE algorithm, 515 ccRCC patients were divided into high/low immune-score groups and high/low stromal-score groups taken the median as the cutoff value. Then, Kaplan-Meier curves were plotted to show the survival difference between each group.
Identi cation of differentially expressed AS events (DEASs) and functional enrichment analysis To screen out the critical AS events that affect survival, we performed differential analysis with a restricted condition of /logFC/>0 and FDR<0.01. Heatmaps and volcano plot were generated using the pheatmap package and ggplot2 in R software, respectively. Consequently, intersection AS events that were upregulated or downregulated in both immune-and stromal-score groups, namely DEASs, were selected out for further analysis, which were exhibited using Venn plot. Furthermore, functional enrichment analysis of DEASs was employed to evaluate the potential function and biological process from the Gene Ontology (GO) terms including the cellular component (CC), biological process (BP), and molecular function (MF) categories, as well as the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Construction of prognostic signatures models
In order to illuminate the prognostic role of the DEASs in ccRCC outcome, we employed the univariate Cox regression analysis to preliminarily screen the overall survival (OS)-related and progression-free survival (PFS)-related DEASs. The least absolute shrinkage and selection operator (LASSO) regression was then applied to identify the nal elimination of potential predictors with non-zero coe cients to make the nal models simplest [35]. Additionally, we further used the multivariate Cox regression analysis to con rm the DEASs involved in nal prognostic signatures, namely OS-related signature and PFS-related signature. According to the results of the multivariate Cox regression analysis, we calculated the risk scores of ccRCC patients. The risk score formula was as follow:riskscore = ∑ n i = 0 PSI × β i , where β is the regression coe cient. We then calculated the risk scores of ccRCC and split patients into high-/low-risk group in basis of the median risk score. Kaplan-Meier survival curves were plotted to display the survival difference and receiver operating characteristic (ROC) curves of 1-, 3-and 5-years were generated to exhibit the discrimination of predictive signatures.

Validation of the independence of prognostic signatures models
We integrated corresponding clinical information included Sex, Age, Grade, AJCC Stage, T, N, M with risk score to validate the independence of prognostic signatures models by univariate and multivariate Cox regression analyses. The analyses were based on two aspects of OS and PFS.

Exhibition of microscopic pro le of immune in ltration
To explore the in uence of DEASs on the formation and layout of the immune microenvironment of ccRCC patients, we further performed the hierarchical consensus clustering to classify the TCGA ccRCC cohort in an unbiased and unsupervised manner with the support of ConsensusClusterPlus package [36]. In our study, the clustering settings used were as follows: maxK= 3; reps=50; pItem=0.8; pFeature=1; title=workDir; clusterAlg="km"; distance="euclidean"; seed=888. The optimal number of clusters was selected according to the Elbow method and the Gap statistic. The association between clusters and clinical survival was analyzed. Additionally, Gene Set Variation Analysis (GSVA) enrichment analysis were employed to shed light on the difference in biological process between clusters by using "GSVA" R packages [37]. Furthermore, single-sample Gene-Set Enrichment Analysis (ssGSEA) algorithm was adopted to quantify the relative abundance of each cell in ltration in the TME of ccRCC, like activated B cell, activated CD4 T cell, activated CD8 T cell, activated dendritic cell and so on, with the aim of exploring the immunological explanation of the survival difference.

Identi cation of differentially expressed SFs (DESFs) and intersection SFs
It has gradually become a consensus that the dysregulation of AS events are orchestrated by a limited number of SFs, whose in uence on tumorigenesis cannot be ignored [38,39]. Differential analysis with a restricted condition of /logFC/>0 and FDR<0.01 was performed to pick out differentially expressed SFs (DESFs) for subsequent research. Heatmaps and volcano plots were generated using the pheatmap package and ggplot2 in R software, respectively. Collectively, intersection SFs that were upregulated or downregulated in both immune-and stromal-score groups, namely DESFs, were selected out, which were exhibited using Venn plot.

Correlation analysis between survival-related AS events and intersection SFs
Based on the result of multivariate Cox regression analysis of OS and PFS, we gained the AS events associated with survival. Next, Spearman correlation analysis between survival-related AS events and intersection SFs was carried out. What's more, the potential SF-AS regulatory network was generated among the signi cant correlation pairs by Cytoscape (version 3.6.1).

Cell line and transfection
Two renal clear cell cancer cell lines including 786O and A498 were purchased from ATCC. The ORF of MBNL1 was inserted into the pcDNA3.1 (+) vector after PCR ampli cation, while an untreated pcDNA3.1 (+) vector was used as a control. Lipofectamine 2000 (Invitrogen) was used for plasmid transfection.
Stable transfected cell lines were established by G418 screening. The expression of MBNL1 in transfected cells was detected by RT-qPCR and Western blot.

qRT-PCR and Western blot
The concentration and purity of RNA and protein were determined after extraction. Reverse Transcription of RNA into cDNA was performed by using Reverse Transcription System (Promega). The diluted cDNA was detected by using the primer of MBNL1 and GAPDH, and the primer sequence of MBNL1 was: 5'-GCTGTTAGTGTCACACCAATTCG-3', 5'-AGGCGATTACTCGTCCATTTTC-3'. Western blot of the protein was performed as follows: after electrophoresis, membrane transfer and sealing, the protein bands were incubated overnight with primary antibodies at 4℃. The primary antibodies used were MBNL1 (66837-1ig, Proteintech), GAPDH (SC-47724, Santa). The membrane was washed with washing buffer solution (TBST) for 3 times. After washing, a secondary antibody was added and incubated for 1 hour. The secondary antibody used was goat anti-mouse (BL001A, Biosharp). Then the membrane was washed by TBST for 3 times, and the electrochemical luminescence was used for detection. Finally, the membrane was scanned and analyzed on the Gel Imager System.

Cell proliferation, Transwell migration and invasion assay, cell cycle and apoptosis
Stably transfected Cell lines of 786O and A498 were cultured in 96-well plates, about 2000 cells were added to each well. Cell Counting Kit-8 (CCK-8; Beyotime) was used for dectection. The absorbance of each well at 450nm was measured to assess the cell proliferation capability. Transwell chambers (3422, Corning) were placed in the wells of 24-well plate which were lled with culture medium containing serum. The lower chamber was immersed in the culture medium, and cell suspension with the serum-free medium was added to the upper chamber. After 48 hours, 4% paraformaldehyde (Biosharp) was used to x the cells on the lower chamber, stained with crystal violet, and counted by a microscope. Cell cycle and apoptosis were analyzed by ow cytometry (FACSCalibur Instrument and CELLQUEST Software, Becton Dickinson). In cell cycle experiment, cells were xed with 70% ice ethanol, treated with RNase A (Sigma), and stained with Propidium iodide (PI; BD Pharmingen). In cell apoptosis experiment, cells were stained with annexin V-uorescein isothiocyanate (BD Pharmingen) and PI. The proportions of cells distributed in each cell cycle and the apoptotic cells were calculated.

Statistical analysis
All statistical analyses were conducted in the R software 3.6.1. The UpSet plot was generated by UpSetR (version 1.3.3) to visualize the intersections between the seven types of AS events [40]. All statistical tests with p < 0.05 (Two side) were statistically signi cant.

Result
Relationship of immune/stromal scores with ccRCC prognosis, identi cation and functional enrichment analysis of DEASs Details of this study design are illustrated in Figure 1 as a owchart. A total of 515 ccRCC patients with complete clinical data and transcriptome data from the TCGA database were included into the follow-up study. On the basis of ESTIMATE analysis, immune scores were distributed between -687.33 and 3339.36, stromal scores ranged from -1428.66 to 1976.04. The ccRCC patients were classi ed into high-and lowimmune/stromal score groups on the boundary of median score. Kaplan-Meier curves displayed that the low immune-score group had longer OS than the high low immune-score group (Figure 2A). It was worth mentioning that the PFS of the low immune-score group seemed to be longer than PFS of the high immune-score group, the OS and PFS of the low stromal-score group seemed to be longer than those of the high stromal-score group, although there was no statistical signi cance ( Figure 2B-2D). These results proved to some extent that immune/stromal score were related to the survival of ccRCC, which was worthy of further data mining.
Considering that a gene may have multiple types of AS events, UpSet plot was generated to visualize the intersecting sets of each AS type, which was illustrated in Supplementary Figure 1. We then conducted differential analysis of AS events expression in the high and low groups based on immune and stromal scores with the purpose of nding other explanations for survival difference. Speci cally, the up-regulated and down-regulated AS events in high/low immune-score group and stromal-score group were showed in volcano plots ( Figure 2E-2F). Considering the signi cant impact of intersection AS events on the prognosis of ccRCC, 400 AS events that were up-regulated in both immune and stromal groups and 461 AS events that were down-regulated in both immune and stromal groups were framed, as shown in Venn diagram, and selected as DEASs for subsequent study ( Figure 2G-2H), which referred to Supplement 1 for details.
In order to assess the potential function and biological process of DEAS events, GO analysis and KEGG pathway were performed. We exhibited the top 30 results of GO analysis incorporated with the top 10 results in BP, the top 10 results in CC and the top 10 results in MF ( Figure 2I). Interestingly, we found that in addition to the enrichment of immune-related response, such as " immune response−regulating cell surface receptor signaling pathway " and "Fc receptor signaling pathway", there were many metabolism-related pathways, particularly GTPase-related, such as "GTPase binding", "GTPase regulator activity" and "GTPase activator activity". The KEGG pathways demonstrated that these DEASs were associated with immune-related pathways, like "Natural killer cell mediated cytotoxicity", "T cell receptor signaling pathway", "Chemokine signaling pathway", "NOD−like receptor signaling pathway", "Fc gamma R−mediated phagocytosis", "TNF signaling pathway", "Fc epsilon RI signaling pathway", "B cell receptor signaling pathway". Besides, some famous cancer-related pathways, such as "Ras signaling pathway", "Rap1 signaling pathway", "Apoptosis", "ErbB signaling pathway" and "Ferroptosis", were enriched ( Figure  2J).

Establishment of two prognosis-related signatures
In order to identi ed survival-related DEASs, we next conducted the univariate Cox regression analysis to screen target AS events from the perspective of OS and PFS (Supplement 2, Supplement 3).In total, 455 OS-related AS events and 464 PFS-related AS events were put into LASSO analysis to obtain the the most prognostic DEASs, respectively ( Figure 3A-D). Subsequently, two prognosis-related signatures related to nal models were constructed by multivariate Cox regression, namely OS-related signature and PFSrelated signature (Table1, Table 2). According to the results of multivariate Cox regression, the formula of OS-related signature was listed as follow: risk score= QKI*-4.617+ SAA2*-0.930+ EPB41L5*2.152+ METTL21A*-1.589+ EPC2*-3.423+ ARHGAP24*1.985+ CTNNB1*-2.365+ FOXN3*1.762+ FAHD2A*1.648+ SHC1*-1.683. Similarly, the formula of PFS-related signature was expanded as follow: risk score=QKI*-4.970+ SAA2*-0.611+ PDCD5*1.622+ C18orf32*-6.596+ SLK*-1.964+ RGS3*1.151+ METTL21A*-0.968+GLIPR1*-1.832+IDH1*-2.601+ARHGAP24*1.930+SF1*3.374. The risk scores of 507 ccRCC patients with completed OS information and 505 ccRCC patients with completed PFS information were calculated. Bounded by the medium risk score, all patients were split into high-risk group and low risk group. Kaplan-Meier curves re ected the fact that high risk group had worse OS than the low-risk group, which was the same as the result of PFS ( Figure 3E, F). Additionally, ROC curves were then drawn with the aim of validating the predictive capacity of these models. The AUC value of OS-related signature for predicting 1-, 3-, and 5-year OS ranged 0.787 from 0.812, and the AUC value of PFS-related signature for predicting 1-, 3-, and 5-year OS were all greater than 0.775, which presented satisfactory performances ( Figure 3G, H).
Inspired by the strong predictive power of the prediction models, we wanted to know whether the risk score had the ability to predict independently. The univariate and multivariate Cox regression analysis were performed to validate the independence of risk scores by combining relevant clinical information, such as Sex, Age, Grade, AJCC Stage, T, N, M, with risk scores (Table 3, Table 4). Surprisingly, risk score was an independent prognostic factor in both OS and PFS, independent of other clinical factors, which was consistent with the results of age.
AS-based clusters signi cantly corelated with prognosis and immune features Driven by the outcomes of enrichment analysis, the unsupervised consensus analysis was conducted to assess the impacts of AS events on immune microenvironment of ccRCC. We then classi ed the TCGA ccRCC patients into three clusters in the help of the consensus matrix heatmap, speci cally, cluster A (n = 146, 27.8%), cluster B (n = 182, 35.9%) and cluster C (n = 179, 35.3%) ( Figure 4A). What's more, Kaplan-Meier analysis of the relationship between clusters and prognosis revealed that different clusters kept different survival patterns. Obviously, cluster A was associated with the best prognosis, while cluster C tended to carry the poorest prognosis, either from the perspective of OS or PFS ( Figure 4B-4C). The rare conclusion that the low immune-score group had a better prognosis further led us to wonder about the speci c effect of AS events on the immune in ltration microenvironment and prognosis of renal cell carcinoma. Through GSVA analysis between clusters, we found that compared with cluster A, cluster C and cluster B were enriched more signi cantly in immune-related pathways, such as "Intestinal immune network for IgA production", "Fc gamma R−mediated phagocytosis", "T cell receptor signaling pathway", "Chemokine signaling pathway", "Cytokine-cytokine receptor interaction", "Natural killer cell mediated cytotoxicity" (Figure 4D-4E). Additionally, cluster C were enriched more signi cantly in immune-related pathways than cluster B, such as "Intestinal immune network for IgA production" and "Natural killer cell mediated cytotoxicity"( Figure 4F).What's more, the subsequent microscopic display of immune cell in ltration also showed that the immune in ltration level of cluster C was the highest, but that of cluster A was the lowest, whether it was immune activated cells or immunosuppressive cells, which was consistent with the previous conclusion ( Figure 4G).

Identi cation of intersection SFs and the regulatory network of survival-related AS events and intersection SFs
A total of 52 SFs were investigated in the present study. To screen out SFs related to tumor immune microenvironment of ccRCC, differential analysis of SFs expression in the high and low groups based on immune and stromal scores were performed, whose analysis method was the same as the differential analysis of AS events. Theses heatmaps exhibited the SF expression between the high immune-score group and low immune-score group, the high stromal-score group and low stromal-score group ( Figure  5A, B). Speci cally, the up-regulated and down-regulated SF in high/low immune-score group and stromal-score group were showed in volcano plots ( Figure 5C, D). Taking intersection SFs as the core of research direction, 3 SFs that were up-regulated in both immune and stromal groups and 1 SFs that was down-regulated in both immune and stromal groups were chose for subsequent study, as shown in Venn diagram ( Figure 5E, F). Furthermore, we found the 4 intersection SFs, including YBX1, MBNL1, ELAVL4 and PCBP2, were positively correlated with immune cell in ltration, including the B cell, CD4+ T cell, CD8+ T cell, macrophage, neutrophil and dendritic cell (Supplementary Figure 2).
SFs are RNA-binding proteins that in uence exon selection and splicing site selection by recognizing cisregulatory elements in precursor mRNAs, suggesting that these survival-related AS events may be regulated by some key SFs [41]. We collected 17 survival-related AS events according to the result of multivariate Cox regression analysis of OS and PFS. Subsequent correlation analysis between survivalrelated AS events and intersection SFs allowed us to nalize 16 SF-related AS events, as shown in Figure  5H. A total of 16 AS events including 9 AS events with low risk (favorable prognosis, green dots) and 7 AS events with high risk (inferior prognosis, red dots) were signi cantly correlated with expression levels of these 4 SFs (bule dots). The negative correlation between AS with good prognosis and SF was the majority (red line), and the positive correlation between AS with poor prognosis and SF was the majority (green line). Motivated by the regulatory network, correlation matrix was further constructed to gure out the relation between intersection AS events and intersection SFs ( Figure 5I). Based on the results of the above regulatory network, correlation analysis and literature review, we selected MBNL1 from 4 hub SFs for subsequent validation in vitro.
Tumor suppressor phenotype of MBNL1 in 786O and A498 cells Stable expression cell lines 786O and A498 were established, and the expression of MBNL1 were veri ed through western blot and RT-qPCR ( Figure 6A). CCK8 assay showed that the proliferation capability of MBNL1 overexpressed cells was declined ( Figure 6B). Meanwhile, the migration and invasion of MBNL1 overexpressed cells were found to be reduced in Transwell assay ( Figure 6C). In the cell cycle and apoptosis experiments, it was found that the overexpressed MBNL1 cells were arrested in the G2/M phase, and the proportion of apoptosis was increased in the cells of transient transfection with MBNL1 expression ( Figure 6D-6E).

Discussion
With the rapid development of high throughput sequencing, it is possible to screen the differences and changes at the micro level through big data analysis to dissect tumor development [42,43]. In recent years, the effects of AS events on various malignant tumors have been demonstrated from different biological perspectives, such as immunology [44,45]. However, the interpretation of the immune perspective of AS events in ccRCC has been emphasized only in a few cases. In the present study, we downloaded the transcriptional data and the relevant clinicopathological information of patients with ccRCC from TCGA and calculated the immune and stromal score of patients by applying the ESTIMATE algorithm. We then identi ed the DEASs based on the immune/stromal group and further investigated the immunobiological response through enrichment analysis and unsupervised cluster analysis, which had certain guiding signi cance for clinical immunotherapy of ccRCC. Furthermore, Cox regression analysis were performed to generate two prognosis-related signatures according to DEASs, namely OS-related signature and PFS-related signature. In the assessment of modeling, the two signatures showed good prediction ability, whose AUC were both greater than 0.75. Additionally, we also focused on the correlation between AS events and SFs and screened out the hub SFs, including YBX1, MBNL1, ELAVL4 and PCBP2.
MBNL1 was selected for cell phenotypic validation and was found to have a signi cant tumor suppressor phenotype.
Intriguingly, we found that there were four overlapped AS events in the multivariate Cox regression of OA and PFS, which were ARHGAP24, METTL21A, QKI and SAA2, respectively. ARHGAP24, also known as FilGAP, is commonly expressed in most tissues, with the highest level of expression in the kidney [46]. It is also a Rac speci c member of the Rho GTPase activating protein family, and is a functional target for cancer cell migration and invasion [46,47]. Previous study had proved that overexpression of ARHGAP24 could inhibit cell proliferation, promote tumor apoptosis and inhibit metastasis of renal cell carcinoma [48,49], which is inconsistent with our conclusion. It may be related to different cell types, and further research is needed.For QKI (quaking), a member of the signal transduction and activation of RNA (STAR) protein family, is related to cancer regulatory mechanism [50].Consistent with our result, the overexpression of QKI can inhibit tumor progression and prolong survival, which is considered to be a tumor suppressor [51][52][53][54] and has a promising application in the development of new targets for tumor targeted therapy. Regarding SAA2, Serum Amyloid A 2, is the most prominent members of the acute phase response (APR) [55]. Jae W. Lee et al had shown that SAA2 levels rised in the circulatory system of patients with liver metastasis, and the high level of circulating SAA2 was associated with poor prognosis, which could promote the formation of metastatic niche [56]. As for METTL21A, no direct studies have been conducted on its relationship with tumor prognosis.
Noteworthily, we found that enrichment analysis results of DEASs were signi cantly enriched in immunerelated pathways, such as "Natural killer cell mediated cytotoxicity", "T cell receptor signaling pathway", "Chemokine signaling pathway", "NOD−like receptor signaling pathway", "Fc gamma R−mediated phagocytosis", "TNF signaling pathway", "Fc epsilon RI signaling pathway", "B cell receptor signaling pathway", which indicated that AS events played a nonnegligible role in shaping the immune microenvironment of ccRCC. In addition, we also realized that several metabolism-related pathways are also involved, like "Rap1 signaling pathway". Jae-Jun Kim et al speculated that cAMP/EPAC/Rap1 pathway might be involved in ETV2-related pathological angiogenesis, which provides convenient metastasis channels and metabolic supply for tumor metastasis [57,58]. Encouraged by the results of enrichment analysis, the unsupervised consensus analysis was further adopted to evaluate the shaping effect of AS events on the immune microenvironment among clusters of ccRCC. According to the ndings of the unsupervised consensus analysis, there was an interesting phenomenon that the higher the level of immune cell in ltration in the immune microenvironment of ccRCC, the more immune-related pathways were enriched, but the worse the survival. Speci cally, cluster A had the best prognosis but the most impoverished immune microenvironment, which seems to contradict the conventional thinking. Nakano et al. reported that due to the high heterogeneity of RCC, so with tumor grade was higher and higher, the lymphocyte in ltration was more and more, but tumor in ltrating lymphocytes itself in tumor tissue didn't mean the antitumor immune effect, did not represent that they were mature lymphocytes[59, 60], and their function may be swallowed by more and more aggressive tumors, resulting in poor prognosis [16]. In other words, the proper performance of anti-tumor immunity did not depend on the number of cells, the key lay in the cell quality, and was signi cantly related to the tumor proliferation activity, which may explain why the more advanced tumor, the worse the effect of immunotherapy. Collectively, the unusual complexity of the immune microenvironment of ccRCC mixed and interwoven by multiple forces provided some thinking for the rational application of clinical immunotherapy.
There is no doubt about the role of SFs in regulating the splicing process [39,61,62]. In the present study, the underlying correlation between AS events and SFs in ccRCC patients was elucidated. Four hub SFs were identi ed in regulatory network, among which MBNL1 has the function of inhibiting tumor cell proliferation and metastasis in 786O and A498 cells, suggesting that the hub SFs may play an important role in the development of ccRCC.
Although the ndings of this study enriched our understanding of the relationship between AS events and the tumor microenvironment of ccRCC, there were some limitations in this study.First, the contents of this study were only derived from the TCGA database, and the prognostic model lacked external validation.Second, the analysis of immune microenvironment stopped at the cellular level, lacking further protein level analysis of immune targets or metabolic levels.

Conclusion
Firstly, we applied the ESTIMATE algorithm to calculate immune and stromal scores, which were highly correlated with clinical outcomes of ccRCC. Secondly, we established ccRCC risk prediction models based on AS events to accurately predict the prognosis of patients, which was helpful for clinical personalized decision making. In addition, a comprehensive analysis of DEAS events provides an immunological perspective to elucidate the possible mechanisms that determine the clinical outcome of ccRCC. More importantly, some hub SFs were screened out which may play an important role in the occurrence and progression of tumors by regulating the corresponding AS events. To sum up, this study not only enriched the research content of the relationship between AS events and immune in ltration of TME, but also elaborated the potential mechanism between AS and SF to a certain extent, providing new insights for the application of clinical immunotherapy and the development of new targets for ccRCC.

Declarations
The dataset generated and analyzed during the current study is available in The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/).

Ethics statement and consent to participate
This study is based on public data and do not refer to the ethical approval and the informed consent.

Consent for publication
Not applicable. 55. Urieli-Shoval S, Linke R, Matzner Y. Expression and function of serum amyloid A, a major acute-phase protein, in normal and disease states. Current opinion in hematology. 2000;7(1):64-9.  Tables   Table 1 Identi cation of differentially expressed alternative splicing (DEAS) events involved in nal prognostic signature by multivariate analysis from the aspect of OS.  Figure 1 The work ow of the study.