Alternative Splicing Events and Splicing Factors Are Prognostic in Adrenocortical Carcinoma

Alternative splicing is involved in the pathogenesis of human diseases, including cancer. Here, we investigated the potential application of alternative splicing events (ASEs) and splicing factors (SFs) in the prognosis of adrenocortical carcinoma (ACC). Transcriptome data from 79 ACC cases were downloaded from The Cancer Genome Atlas database, and percent spliced-in values of seven splicing types were downloaded from The Cancer Genome Atlas SpliceSeq database. By the univariate Cox regression analysis, 1,839 survival-related ASEs were identified. Prognostic indices based on seven types of survival-related ASEs were calculated by multivariate Cox regression analysis. Survival curves and receiver operating characteristic curves were used to assess the diagnostic value of the prognostic model. Independent prognosis analysis identified several ASEs (e.g., THNSL2| 54469| ME) that could be used as biomarkers to predict the prognosis of patients with ACC accurately. By analyzing the co-expression correlation between SFs and ASEs, 188 highly correlated interactions were established. From the protein interaction network, we finally screened six hub SFs, including YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4, whose expression levels were significantly related to the overall survival and prognosis of ACC. Our findings provide a reliable model for predicting the prognosis of ACC patients based on aberrant alternative splicing patterns.

Introduction

Alternative splicing of pre-messenger RNA (mRNA) produces transcript isoforms for 95% of human genes, increases protein diversity, and provides functional diversity at various regulation level (Pan et al., 2008). There are seven types of splicing patterns, including alternate acceptor site (AA), alternate donor site (AD), alternate promoter (AP), alternate terminator (AT), exon skip (ES), mutually exclusive exons (ME), and retained intron (RI), as listed in The Cancer Genome Atlas (TCGA) SpliceSeq database (Ryan et al., 2016). Splicing factors (SFs) are involved in the removal of introns to create mature mRNAs, a process catalyzed by a large complex termed spliceosome (Yan et al., 2019). Alterations in SF expression lead to missplicing of key cancer-associated genes (Anczukow and Krainer, 2016; Lee and Abdel-Wahab, 2016). Aberrant alternative splicing events (ASEs) have been frequently observed in cancers and is recognized as an important signature for tumorigenesis and related pathologies, such as initiation and development of cancer (Oltean and Bates, 2014; Chen and Weiss, 2015; Lee and Abdel-Wahab, 2016), cancer metabolism (Kozlovski et al., 2017), cancer immunotherapy (Frankiw et al., 2019), cancer drug resistance (Siegfried and Karni, 2018), and so on.

Adrenocortical carcinoma (ACC) is a rare aggressive tumor with poor prognosis and less than 40% survival rate in 5 years (Jouinot and Bertherat, 2018; Crona and Beuschlein, 2019). Recent studies highlighted that specific molecular signatures could predict the survival and prognosis of ACC patients, which came from genomic approaches, including transcriptome, exome or whole-genome sequencing, chromosome alterations, methylome, and miRnome (Barreau et al., 2013; Patel et al., 2013; Assie et al., 2014; Szabo et al., 2014; Jouinot and Bertherat, 2018). However, only limited studies have focused on the potential roles of alternative splicing patterns in the pathogenesis of ACC. The present study aims to investigate the relationship between seven types of ASEs and SFs with the prognosis of ACC. Our findings provide a new path to identify potential targets for the diagnosis and treatment of ACC.

Results

Overview of Alternative Splicing Events in Adrenocortical Carcinoma

Transcriptome data from the TCGA database provide identity information for up to 56,754 transcripts, which represents a key resource for exploring ASEs in tumors concurrently deposited in the database. Individual ASE is assigned a unique annotation in the TCGA SpliceSeq database; for instance, in the term CIRBP| 46443| ES, CIRBP is the official gene symbol, 46443 is the unique ID number of a specific ASE, and ES represents the type of the alternative splicing pattern. Full datasets of splicing events of 79 ACC cases have been deposited in TCGA SpliceSeq, which contains 34,419 ASEs corresponding to 8,994 parent genes. The splicing event data are presented as the percent-spliced-in (PSI), a score to evaluate the level of specific ASEs. The data were firstly filtered by deleting the unreliable ASEs with mean PSI of less than 0.05 and a standard deviation of PSI of less than 0.01. This generated a total of 22,521 ASEs from 8,040 parent genes. As summarized in the UpSet plot (Figure 1A), ES and AT were the top splicing types with high frequency in most genes, whereas ME was the least frequent splicing type.

FIGURE 1

Figure 1. Overview of ASEs and SR-ASEs profiling in ACC. (A) UpSet plot of top 50 ASEs with frequency in ACC. (B) Volcano plot of 22,521 ASEs. Red points represent the 1,839 SR-ASEs under P < 0.01. Transverse axis is the Z score of univariate Cox regression analysis. (C) UpSet plot of top 50 SR-ASEs with frequency. (D) Bubble plot of the top 20 most significant SR-ASEs. Size of point represents -log10(P-value), and the color of points represents P-value. Terms of ASEs contain three parts: the gene name, a unique splicing event ID number, and alternative splicing type. AA, alternate acceptor site; AD, alternate donor site; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.

Identification of Survival-Related Alternative Splicing Events

To explore the relationship between the alternative splicing pattern and the prognosis of ACC, we performed a univariate Cox regression analysis by comparing ASEs with the overall survival of ACC patients. A total of 1,839 ASEs were significantly associated with overall survival in ACC patients (P < 0.01), termed survival-related ASEs (SR-ASEs; Figure 1B). The upSet plot showed that ES, AT, and AP were the most frequent SR-ASEs, whereas only a small number of genes displayed a combination of multiple splicing forms (Figure 1C). Supplementary Figure 1 listed the top 20 most significant SR-ASEs in each splicing pattern according to their Z score and P-value. Taking all types together, splicing events within the parent genes CIRBP, BLOC1S1, TRAFD1, UNG, EIF6, METTL15, CMC2, HM13, KLHL7, TECPR2, DNAJC12, DUT, and MPND were highly related to the overall survival of ACC patients (Figure 1D). Survival curves with a cutoff at the median PSI value of the top four ASEs showed a remarkable difference (Figures 2A–D).

FIGURE 2

Figure 2. Kaplan–Meier survival curves and ROC curves for five alternative splicing type: CIRBP| 46443| ES (A), BLOC1S1| 22229| AP (B), TRAFD1| 24580| AD (C), CIRBP| 46445| ES (D), and THNSL2| 54469| ME (E). Seventy-nine ACC patients were divided into high- and low-risk groups based on the median of PSI. Red line indicates the high PSI score group, and the blue line indicates the low PSI score group.

Independent Prognostic Predictors of Survival-Related Alternative Splicing Events in Adrenocortical Carcinoma

To identify independent prognostic indices and to explore the relationship between aberrant types of ASEs and ACC survival outcomes, we performed the multivariate Cox regression analysis for each splicing type to build a prognostic model. Lasso regression analysis was done to select the most significant SR-ASEs by avoiding overfitting. The median value of risk scores was then used to stratify the 79 ACC patients into low- and high-risk groups. Kaplan–Meier method was used to analyze the efficacy of the prognostic indices to predict the overall survival. The plotted survival curves and receiver operating characteristic (ROC) curves were shown in Figure 3. Significant differences in survival curves were observed in individual splicing type as well as all together (ALL), indicating that each alternative splicing type could be recognized as an independent prognostic indicator (Figure 3). The greatest difference in overall survival curves was observed in ES type (P = 6.684e-11), which is the most frequent splicing type among ASEs and SR-ASEs (Figures 1A,C). The area under the curve (AUC) of each ROC curve was more than 0.7, indicating the predictive efficiency of the eight models. We found that the AUC for AA type is 0.978, which is higher than all others (Figure 3), indicating that the prognostic indices based on AA type demonstrated the greatest efficacy in stratifying patients.

FIGURE 3

Figure 3. Kaplan–Meier survival curves and ROC curves for each alternative splicing type: AA (A), AD (B), AP (C), AT (D), ES (E), ME (F), RI (G), and ALL (H) model. Seventy-nine ACC patients were divided into high- and low-risk groups based on the median of prognostic indices. Red line indicates the high-risk group, and blue line indicates the low-risk group. Prediction time of ROC curve is 1 year. AA, alternate acceptor site; AD, alternate donor site; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.

We also performed multivariate Cox regression analysis to evaluate the effect of other clinical parameters, including sex, tumor stage, T and N stages in the tumor–node–metastasis (TNM) classification in Table 1. The significant hazard ratios (HRs) for T stage and tumor stage were 3.349 (95% confidence interval [CI]: 2.075–5.403) and 2.886 (CI: 1.819–4.579), respectively (Table 1). HRs of AA, AD, AP, AT, ES, RI, and ALL models are also shown in Table 1. Among the seven types, the HR for ME type was 1.270 (CI: 1.093–1.476), the highest HR value (Table 1).

TABLE 1

Table 1. Univariate Cox regression analysis of AA, AD, AP, AT, ES, ME, RI, and ALL models.

Hazard ratios in each of AA, AD, AP, AT, ES, RI, and ALL models with clinical parameters are shown in Figure 4. The HRs for all the AS type were under a significant level (P-value < 0.05), and ME type has the highest HR value of 1.681 (CI: 1.344–2.103) (Figure 4F). Moreover, THNSL2| 54469| ME ranks that most significant event in the ME type (Supplementary Figure 1F). THNSL2| 54469| ME was the top significant SR-ASE in ME model to predict the prognostic status of ACC cases. Therefore, the most significant SR-ASEs in ME pattern is THNSL2| 54469| ME (Figure 2E). To test the accuracy of ME model of THNSL2| 54469| ME, the survival state of ACC patients could be significantly classified into high- and low-PSI groups according to the median value of PSI scores of THNSL2| 54469| ME (Figure 2E). These results indicate that THNSL2| 54469| ME could be the best independent prognostic indicator to predict the prognosis of ACC cases.

FIGURE 4

Figure 4. Multivariate Cox regression analysis for AA (A), AD (B), AP (C), AT (D), ES (E), ME (F), RI (G), and ALL (H) model. Hazard ratio is shown as hazard ratio (95% confidence interval). AA, alternate acceptor site; AD, alternate donor site; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.

Also, the AUCs for 1, 3, and 5 years of ME model are predicted to be 0.767 (Figure 3F), 0.792 (Figure 5A), and 0.829 (Figure 5B), respectively. The number of deaths increased in the high-risk group, with short survival time (Figures 5C,D). In summary, those results indicated that the constructed ME model had great efficacy in predicting the prognosis of ACC patients.

FIGURE 5

Figure 5. Overview of ME model. ROC curves of ME type with prediction time of 3 years (A) and 5 years (B). (C) Risk factors of ME model of 79 ACC cases. Patients were divided into low- and high-risk groups based on the median of prognostic indices. Red indicates high-risk group, and green indicates low-risk group. (D) Survival time and survival status of ME model of 79 ACC cases. Red indicates the deaths of the cases. Green indicates cases who are alive.

Enrichment Analysis of Survival-Related Alternative Splicing Events

To explore the biological processes and signal pathways related to alternative splicing in the progression of ACC, enrichment analysis of the SR-ASEs parent genes was performed by gene ontology and pathway analysis in Metascape. The most significantly enriched terms were regulation of mitotic cell cycle, cofactor metabolic process, covalent chromatin modification, cell cycle G2/M phase transition, and mitochondrion organization (Figure 6), pathways that are all involved in tumorigenesis.

FIGURE 6

Figure 6. Bar graph of enriched terms across parent genes of SR-ASEs by Metascape, colored by P-values.

Correlation Network of Survival-Related Alternative Splicing Events and Splicing Factors Expression

Alteration of ASEs is largely attributable to changes in SF expression. We then extracted the expression data of 404 SFs, summarized by a previous study (Seiler et al., 2018), from the transcriptome data of the 79 ACC patients. Principal component analysis (PCA) plots could discriminate the distribution pattern of SF expression levels according to survival state, tumor stage, and TNM classification (Supplementary Figure 2), suggesting an impact of altered SF expression on the ACC outcome.

As our current knowledge is incapable of dissecting the sequence specificity for each SF, we could not establish a direct network for SF-regulated ASEs. Thus, we analyzed the relationship between SFs and SR-ASEs based on their co-expression patterns using the Spearman method, which has been widely used in alternative splicing studies (He et al., 2018; Xiong et al., 2018; Zhang et al., 2020). A total of 188 highly correlated interactions between SFs and SR-ASEs were detected with a correlation coefficient larger than 0.65 (Figure 7A and Table 2). Hub SFs with no less than five SR-ASE connections were HSPA1B, YBX1, SRPK1, SART1, PRCC, ILF2, SNRPG, SNRPE, SF3B4, BUD13, INTS4, and CLK2 (Table 3). Among these SFs, interestingly, SNRPE was exclusively correlated with AT-type ASEs, suggesting a specific role in regulating terminal exon selection (Figure 7A). Moreover, HSPA1B was exclusively correlated with detrimental ASEs showing negative correlations with overall survival, indicating a potential causal role of HSPA1B in the progression of ACC (Figure 7A).

FIGURE 7

Figure 7. Correlation network of SR-ASEs and SFs. (A) Co-expression networks between parent genes of SR-ASEs and SFs. Yellow square indicates SFs, red circle indicates the poor prognostic SR-ASEs, and blue circle indicates the better prognostic SR-ASEs. Red line indicates that SFs were positively correlated with the SR-ASE, and blue line indicates SFs were negatively correlated with the SR-ASE. (B) Protein–protein interaction network of SFs.

TABLE 2

Table 2. Top 20 significant associations between SR-ASEs and SFs.

TABLE 3

Table 3. Two MCODE modules of protein–protein interaction network of parent genes of SR-ASEs and SFs.

The expression data showed that SF YBX1 was positively correlated with the PSI values of C16orf13| 32916| ES, COX4I1| 37906| RI, FHAD1| 747| AT, MYL6| 22381| AA, PI4K2A| 12728| AP, ASH2L| 83369| AP, PTAR1| 86546| AT, and IRF3| 51012| ES and was negatively correlated with the PSI values of FHAD1| 749| AT, AIG1| 77970| AT, PI4K2A| 12729| AP, STARD3NL| 79286| ES, and ASH2L| 83368| AP; SF SNRFE was positively correlated with the PSI values of VWA8| 25742| AT, MSI2| 42617| AT, and PELI3| 17034| AT and was negatively correlated with the PSI values of VWA8| 25741| AT, MSI2| 42616| AT, and PELI3| 17033| AT (Figure 7A). We then constructed the Kaplan–Meier survival curves for each ASEs related with YBX1 and SNRFE and found that high levels of PSI values for all the negatively correlated ASEs have a better prognostic except FHAD1| 747| AT, whereas the low levels of PSI values for all the positively correlated ASEs have a better prognostic (Figure 8). Interestingly, six pairs of ASEs from three parent genes were observed for these two SFs in the network. For example, ASH2L| 83368| AP and ASH2L| 83369| AP are two alternative splicing sites for ASH2L exon 1 selection. PI4K2A| 12728| AP and PI4K2A| 12729| AP are two alternative splicing sites for ASH2L exon 1 selection (Figures 8A,B). Our results indicate that promoter selections of ASH2L and PI4K2A are important for tumorigenesis of ACC. Similar results were also observed for the SF SNRFE (Figures 8C,D). These results suggest that terminator selections of MSI2 and PELI3 and VWA8 are important for the progression of ACC development (Figures 8C,D).

FIGURE 8

Figure 8. Kaplan–Meier survival curves for ASEs that correlated with splicing factors YBX1 and SNRFE. (A) ASEs that negatively correlated with YBX1. (B) ASEs that positively correlated with YBX1. (C) ASEs that negatively correlated with SNRFE. (D) ASEs that positively correlated with SNRFE. Seventy-nine ACC patients were divided into high- and low-risk groups based on the median of PSI. Red line indicates the high PSI score group, and blue line indicates the low PSI score group.

Considering that SFs would influence each other when regulating the work mode of the spliceosome, we constructed the protein–protein internecion network to illustrate the interactions among SFs using the Search Tool for the Retrieval of Interacting Genes database (Figure 7B). A total of 66 nodes and 485 edges were revealed, with the interaction score of 0.900 (Figure 7B). Module analysis was done to identify hub genes. Two modules were further identified by the app MCODE in Cytoscape. Module 1 has 30 genes with a score of 27, and module 2 has four genes with a score of 8 (Table 3). By combining the ASE correlation results and protein interaction networks, we finally identified six hub SFs, including YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4 (Table 4). Overall, these results indicate that these six hub genes may play an important role in the progression of ACC by regulating the pattern of SR-ASEs.

TABLE 4

Table 4. Counts of SFs correlated with SR-ASEs in correlation network.

Analysis of Hub Splicing Factors

To confirm whether these six genes, YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4, are high-risk factors, the overall survival and expression level of these genes were further investigated. The results showed that the survival rate of patients with high expression of hub SFs was significantly lower than that in the low expression group (P < 0.001) (Figure 9), and the results were consistent with those from the Gene Expression Profiling Interactive Analysis (GEPIA) database (Supplementary Figure 3).

FIGURE 9

Figure 9. Relationship between hub SFs and survival curve, tumor stage, T stage and lymph node state. (A) YBX1. (B) SART1. (C) PRCC. (D) SNRPG. (E) SNRPE. (F) SF3B4. *P < 0.05; **P < 0.01; ***P < 0.001 and n.s. indicates no significance.

The correlation between the expression levels of six SFs and clinical information, including clinical stages, tumor stage, and lymph nodes stage, was further analyzed. We found that the expression levels of SFs increased along with the clinical stages. The expression level of six hub SFs in tumor patients has an increasing trend with the tumor progress stage. Among six SFs, the level of SF3B4 in ACC cases with the tumor stage III and the level of PRCC and SF3B4 in ACC cases with stage IV were statistically different from those in stage I (Figure 9), which is consistent with the results from the GEPIA database (Supplementary Figure 3). The levels of the SFs in ACC cases with lymph node stage also increased, and the expression of YBX1 was statistically changed. The results indicate that the expression levels of these SFs are closely related to the survival time and prognosis of patients with ACC (Figure 9).

Considering the lack of normal control tissue for ACC in the TCGA database, two microarray datasets from Gene Expression Omnibus (GEO) database were used to compare the expression of SFs between tumor tissues and normal controls (Figure 10). The mRNA levels of YBX1 and SNRPE were increased in tumor samples in the two microarrays, although the expression levels only showed statistically significant differences in GSE19750 (Figure 10). Combining the mRNA level of YBX1 and SNRPE in different tissues and the effects of YBX1 and SNRPE on the survival curve of ACC cases (Figure 7), we could conclude that YBX1 and SNRPE could exert positive regulation in the progression of ACC.

FIGURE 10

Figure 10. Expression level of hub SFs of ACC or normal samples in the microarray: (A) GSE19750 and (B) GSE33371. *P < 0.05; ***P < 0.001 and n.s. indicates no significance.

We next performed a Cox regression analysis to evaluate the prognostic value of these six hub SFs and other clinical parameters. Results showed that the HRs of these six genes ranged from 1.003 to 1.156; although the HR is lower than the T and N stages, it statistically significant for all the six SFs (Supplementary Figure 4). The AUC of each ROC curve is higher than 70% (Supplementary Figure 4), indicating that all of these six SFs could be used to classify the stages of ACC.

Discussion

Adrenocortical carcinoma is a rare malignancy tumor with a poor prognosis. Currently, surgery is the only available curative treatment option (Crona and Beuschlein, 2019). Recent studies (Barreau et al., 2013; Patel et al., 2013; Assie et al., 2014; Szabo et al., 2014; Jouinot and Bertherat, 2018) highlighted that genomic approaches derived from TCGA database and GEO datasets could provide specific molecular signatures for the diagnosis and prognosis of ACC. At present, few studies focused on the different isoforms of alternative splicing in ACC (Bie et al., 2019). The present study is to explore the aberrant of ASEs and hub SFs to develop novel diagnostic and prognostic markers for ACC.

Seven types of ASEs were investigated in this study. ES type was the top splicing type with high-frequency ASEs and SR-ASEs (Figures 1A,C), indicating the ES is the dominant alternative splicing type in ACC.

In this study, we identified 1,839 SR-ASEs by univariate Cox regression analysis. Interestingly, different ASEs in the same gene could exert opposite functions in the overall survival of ACC, indicating that the parent genes of these ASEs may play an important role in ACC development (Figures 1D, 8). Because ME type has the highest HR value (Figure 4F) THNSL2| 54469| ME ranks the most significant event in the ME type (Supplementary Figure 1F). Therefore, THNSL2| 54469| ME could be used as an independent prognostic indicator to predict the prognosis of ACC cases (Figure 2E). THNSL2 encodes a threonine synthase-like protein and has multiple transcript variants. The function of THNSL2 and the alternative splicing of THNSL2| 54469| ME could be further investigated. Enrichment analysis revealed several important pathways, including regulation of the mitotic cell cycle and cell cycle G2/M phase transition, which could impact the occurrence and development of ACC (Figure 6).

Genes involved in the same biological process or signaling pathway are usually co-regulated in the disease context. Co-expression network analysis has been widely used to dissect the functional gene panels in large datasets, including alternative splicing studies (He et al., 2018; Xiong et al., 2018; Zhang et al., 2020). One hundred eighty-eight highly correlated interactions between SFs and SR-ASEs were identified (Figure 7). SR-ASEs that were positively correlated with SFs have a poor prognostic value, whereas SR-ASEs that were positively correlated with SFs have a poor prognostic value (Figures 7, 8). It provided a new insight for the molecular mechanism of alternative splicing in ACC.

The number of pairs of ASEs from one parent gene was observed for SFs in the network (Figure 7), and six pairs of ASEs related with YBX1 and SNRFE were illustrated in detail (Figure 8). We observed that the pairs of ASEs conferred the opposite function for the ACC progress shown in the overall survival curve (Figure 8), indicating that specific exons were important, such as the alternative promoters’ selection of ASH2L, and alternative terminator selection of MSI2. Moreover, among these six pairs of ASEs, only alternative promoters of ASH2L have been reported in embryonal carcinoma (Alagaratnam et al., 2013).

A previous study on ASEs in endometrial cancer also identified YBX1 as the hub SF. One pair of ASEs that correlated with YBX1 was identified in that study: DNAH9-AT-39292 and DNAH9-AT-39293 (Wang et al., 2019). We could conclude that, firstly, it seems YBX1 more specifically regulates the first exon, as the splicing type is AT type in both studies; secondly, YBX1, as the hub SF, might regulate specific genes in different cancer types.

Six hub SFs, including YBX1, SART1, PRCC, SNRPG, SNRPE, and SF3B4, were identified in this study. We found that the expression level of hub SFs was negatively correlated with the survival time and survival state of ACC patients (Figure 9). There was still some evidence that the expression level of hub SFs in tumor tissues was higher than that in normal tissues (Figure 10). The altered expression level of hub SFs has been reported in multiple types of cancer, such as YBX1 (Rossner et al., 2016; Chen et al., 2019), SART1 (Ishida et al., 2000; Sasatomi et al., 2000; Yutani et al., 2001), SNRPE (Tamura et al., 2007; Jia et al., 2011), and SF3B4 (Liu et al., 2018). Other studies have shown that SNRPE (Quidville et al., 2013) and SF3B4 (Shen and Nam, 2018) could develop a new therapeutic agent in cancer. Further study found that the inactivation of SF3B4 inhibited liver tumorigenesis in vitro and in vivo (Shen et al., 2018). In terms of molecular mechanism, SF3B4 triggers SF3b complex to splice tumor suppressor KLF4 transcript to non-functional skipped exon transcripts, downregulates the transcriptional activity of p27Kip1, and upregulates the transcriptional activation of Slug genes (Shen et al., 2018). However, the functions of hub SFs in ACC development need to be further studied.

Our present study performed a bioinformatic analysis of SR-ASEs and hub SFs and provided insight into the function of aberrant ASEs in ACC development and progression. SR-ASEs and hub SFs identified in this study could be potential targets for the diagnosis of ACC patients.

Materials and Methods

Data Collection

The transcriptome data of 79 ACC cases and the corresponding clinical information, including survival time, survival status, sex, tumor stage, T stage, and lymph node metastasis, were downloaded from TCGA¹. Seven types of ASEs of 79 ACC cases were downloaded from the TCGA SpliceSeq database² (Ryan et al., 2016). ASEs were quantified using the PSI, rating from 0 to 1. Splicing event data of 79 ACC cases included 34,419 ASEs, corresponding to 8,994 genes. To better track the AS events, the name of the ASEs contains three parts: the official gene symbol, a unique splicing event ID number, and splicing type. The number of ASEs for different types of ASEs are shown in the UpSet plot (Figure 1A).

Screening of Survival-Related Alternative Splicing Events

The missing value of PSI in seven AS event types was supplemented by impute.knn function using nearest neighbor averaging method with k = 10, rowmax = 0.5, colmax = 0.8, and other default setting in R. Then, the PSI data of 79 ACC cases were filtered by deleting the ASEs with mean PSI of less than 0.05 and standard deviation of PSI of less than 0.01. The filtered PSI data, including 22,521 ASEs from 8,040 genes in 79 ACC cases, were used for visualization and subsequent analysis. The survival time and survival status of 79 ACC cases were integrated with the filtered PSI data. Then, SR-ASEs were selected using univariate Cox regression analysis with a threshold set to a P-value < 0.01 and were visualized by volcano plot and bubble plot.

Construction of Prognosis Model of Survival-Related Alternative Splicing Events

Multivariate Cox regression analysis was performed to calculate the prognostic indices of each for each type of splicing pattern. To prevent over-fitting of the model, seven different types of SR-ASEs and ALL SR-ASEs were analyzed by lasso regression analysis, making the Lambda values at smaller level. The SR-ASEs with high correlation have been removed to guarantee the accuracy of the model. Risk factors were calculated using the following formula βSR-ASE1 × PSISR-ASE1 + βSR-ASE2 × PSISR-ASE2 + …… + βSR-ASEn × PSISR-ASEn, where β corresponded to the regression coefficient. The samples were stratified to high- and low-risk groups based on the median of risk scores. The survival curve and ROC curve of the 79 ACC cases were used to evaluate the accuracy of the model. Univariate and multivariate Cox regression analyses were used to determine whether the prognostic indices could be used as an independent prognostic factor for predicting the prognosis of ACC cases.

Enrichment Analysis of Survival-Related Alternative Splicing Events

To illustrate the biological functions and pathway associated with SR-ASEs of ACC cases, gene enrichment analysis for parent genes of SR-ASEs was performed in the online database Metascape³ (Zhou et al., 2019) with the default setting. Metascape is an online analytical tool for pooled gene annotation, enrichment analysis, and protein interaction analysis.

Construction of Splicing Correlation Network

Four hundred four SFs summarized from the previous study were used in this study for SF analysis (Seiler et al., 2018). They collected and prioritized the final list of 404 SF genes by compiling and filtering spliceosome and splicing related genes from three sources, which were all experimentally validated in the previous study (Barbosa-Morais et al., 2006; Hegele et al., 2012; Cvitkovic and Jurica, 2013). Source 1 (Hegele et al., 2012) was from a comprehensive yeast two-hybrid study using spliceosome components as bait. Source 2 (Barbosa-Morais et al., 2006) included 254 SFs and splicing related proteins. Source 3 was from SpliceosomeDB (Cvitkovic and Jurica, 2013). Of 404 SFs, 390 were extracted from the transcriptome data of 79 ACC cases. Co-expression network analysis calculated by spearman method was used in the construction of alternative splicing regulation networks to screen the important SFs. Cytoscape visualized a highly correlated interaction network with the threshold of correlation coefficient of 0.65.

Also, protein–protein interaction network of SFs was established to find the important SFs in the progression of ACC, using the online database Search Tool for the Retrieval of Interacting Genes⁴ (version 11.0). Then, the hub SFs were selected by the MCODE app of Cytoscape. Six prognostic-related hub SFs for patients with ACC were identified by combining the co-expression networks analysis and protein–protein internecion networks.

Analysis of Hub Splicing Factors

Correlations of the mRNA expression level of hub SFs with clinical data (survival curve, tumor stage, T stage, and lymph node state) of 79 ACC cases were analyzed. The online database GEPIA (Tang et al., 2017)⁵, which contains the TCGA database with the GTEx database, was further used to evaluate the SFs on clinical features. Also, three gene expression microarray data were downloaded from GEO⁶ for illustrating the expression levels of hub SFs between normal tissue and tumor tissues. GSE19750 (Demeure et al., 2013) includes 44 ACC samples and 4 control samples. GSE33371 (Heaton et al., 2012) includes 33 ACC samples and 10 control samples.

Data Availability Statement

The transcriptome data of 79 ACC cases and the corresponding clinical information can be downloaded from The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov/tcga/). Seven types of alternative splicing events (ASEs) of 79 ACC cases downloaded from TCGA SpliceSeq database (https://bioinformatics.mdanderson.org/TCGASpliceSeq/). Gene expression microarray data GSE19750 and GSE33371 were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo).

Author Contributions

ZW and LL conceived and designed the project. JL and YH downloaded and analyzed the data. JL and LL drafted the manuscript. All authors contributed to text revision and discussion.

Funding

This study was supported by the National Natural Science Foundation of China (no. 81722007) and the National Health Commission of China (no. 2017ZX1030440 2001-008).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We appreciate the generosity of TCGA, TCGA SpliceSeq, GEO, and GEPIA for sharing the huge amount of data. We also thank members of our laboratory for supporting this work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00918/full#supplementary-material

Footnotes

1. ^ https://tcga-data.nci.nih.gov/tcga/
2. ^ https://bioinformatics.mdanderson.org/TCGASpliceSeq/
3. ^ http://metascape.org
4. ^ http://string-db.org
5. ^ http://gepia.cancer-pku.cn/index.html
6. ^ http://www.ncbi.nlm.nih.gov/geo

References

Alagaratnam, S., Harrison, N., Bakken, A. C., Hoff, A. M., Jones, M., Sveen, A., et al. (2013). Transforming pluripotency: an exon-level study of malignancy-specific transcripts in human embryonal carcinoma and embryonic stem cells. Stem Cells Dev. 22, 1136–1146. doi: 10.1089/scd.2012.0369

PubMed Abstract | CrossRef Full Text | Google Scholar

Anczukow, O., and Krainer, A. R. (2016). Splicing-factor alterations in cancers. RNA 22, 1285–1301. doi: 10.1261/rna.057919.116

PubMed Abstract | CrossRef Full Text | Google Scholar

Assie, G., Letouze, E., Fassnacht, M., Jouinot, A., Luscap, W., Barreau, O., et al. (2014). Integrated genomic characterization of adrenocortical carcinoma. Nat. Genet. 46, 607–612. doi: 10.1038/ng.2953

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbosa-Morais, N. L., Carmo-Fonseca, M., and Aparicio, S. (2006). Systematic genome-wide annotation of spliceosomal proteins reveals differential gene family expansion. Genome Res. 16, 66–77. doi: 10.1101/gr.3936206

PubMed Abstract | CrossRef Full Text | Google Scholar

Barreau, O., Assie, G., Wilmot-Roussel, H., Ragazzon, B., Baudry, C., Perlemoine, K., et al. (2013). Identification of a CpG island methylator phenotype in adrenocortical carcinomas. J. Clin. Endocrinol. Metab. 98, E174–E184. doi: 10.1210/jc.2012-2993

PubMed Abstract | CrossRef Full Text | Google Scholar

Bie, J., Liu, K., Song, G., Hu, X., Xiong, R., Zhang, X., et al. (2019). ENST00000489707.5 is a preferred alternative splicing variant of PTK7 in adrenocortical cancer and shows potential prognostic value. Med. Sci. Monit. 25, 8326–8334. doi: 10.12659/MSM.919818

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., and Weiss, W. A. (2015). Alternative splicing in cancer: implications for biology and therapy. Oncogene 34, 1–14. doi: 10.1038/onc.2013.570

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Q. F., Li, W., Wu, P., Shen, L., and Huang, Z. L. (2019). Alternative splicing events are prognostic in hepatocellular carcinoma. Aging 11, 4720–4735. doi: 10.18632/aging.102085

PubMed Abstract | CrossRef Full Text | Google Scholar

Crona, J., and Beuschlein, F. (2019). Adrenocortical carcinoma - towards genomics guided clinical care. Nat. Rev. Endocrinol. 15, 548–560. doi: 10.1038/s41574-019-0221-227

CrossRef Full Text | Google Scholar

Cvitkovic, I., and Jurica, M. S. (2013). Spliceosome database: a tool for tracking components of the spliceosome. Nucleic Acids Res. 41, D132–D141. doi: 10.1093/nar/gks999

PubMed Abstract | CrossRef Full Text | Google Scholar

Demeure, M. J., Coan, K. E., Grant, C. S., Komorowski, R. A., Stephan, E., Sinari, S., et al. (2013). PTTG1 overexpression in adrenocortical cancer is associated with poor survival and represents a potential therapeutic target. Surgery 154, 1405–1416. doi: 10.1016/j.surg.2013.06.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Frankiw, L., Baltimore, D., and Li, G. (2019). Alternative mRNA splicing in cancer immunotherapy. Nat. Rev. Immunol. 19, 675–687. doi: 10.1038/s41577-019-0195-197

CrossRef Full Text | Google Scholar

He, R. Q., Zhou, X. G., Yi, Q. Y., Deng, C. W., Gao, J. M., Chen, G., et al. (2018). Prognostic signature of alternative splicing events in bladder urothelial carcinoma based on Spliceseq data from 317 cases. Cell Physiol. Biochem. 48, 1355–1368. doi: 10.1159/000492094

PubMed Abstract | CrossRef Full Text | Google Scholar

Heaton, J. H., Wood, M. A., Kim, A. C., Lima, L. O., Barlaskar, F. M., Almeida, M. Q., et al. (2012). Progression to adrenocortical tumorigenesis in mice and humans through insulin-like growth factor 2 and beta-catenin. Am. J. Pathol. 181, 1017–1033. doi: 10.1016/j.ajpath.2012.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegele, A., Kamburov, A., Grossmann, A., Sourlis, C., Wowro, S., Weimann, M., et al. (2012). Dynamic protein-protein interaction wiring of the human spliceosome. Mol. Cell 45, 567–580. doi: 10.1016/j.molcel.2011.12.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishida, H., Komiya, S., Inoue, Y., Yutani, S., Inoue, A., and Itoh, K. (2000). Expression of the SART1 tumor-rejection antigen in human osteosarcomas. Int. J. Oncol. 17, 29–32. doi: 10.3892/ijo.17.1.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, D., Wei, L., Guo, W., Zha, R., Bao, M., Chen, Z., et al. (2011). Genome-wide copy number analyses identified novel cancer genes in hepatocellular carcinoma. Hepatology 54, 1227–1236. doi: 10.1002/hep.24495

PubMed Abstract | CrossRef Full Text | Google Scholar

Jouinot, A., and Bertherat, J. (2018). MANAGEMENT OF ENDOCRINE DISEASE: adrenocortical carcinoma: differentiating the good from the poor prognosis tumors. Eur. J. Endocrinol. 178, R215–R230. doi: 10.1530/EJE-18-0027

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozlovski, I., Siegfried, Z., Amar-Schwartz, A., and Karni, R. (2017). The role of RNA alternative splicing in regulating cancer metabolism. Hum. Genet. 136, 1113–1127. doi: 10.1007/s00439-017-1803-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. C., and Abdel-Wahab, O. (2016). Therapeutic targeting of splicing in cancer. Nat. Med. 22, 976–986. doi: 10.1038/nm.4165

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Li, W., Pang, Y., Zhou, Z., Liu, S., Cheng, K., et al. (2018). SF3B4 is regulated by microRNA-133b and promotes cell proliferation and metastasis in hepatocellular carcinoma. eBio Med. 38, 57–68. doi: 10.1016/j.ebiom.2018.10.067

PubMed Abstract | CrossRef Full Text | Google Scholar

Oltean, S., and Bates, D. O. (2014). Hallmarks of alternative splicing in cancer. Oncogene 33, 5311–5318. doi: 10.1038/onc.2013.533

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, Q., Shai, O., Lee, L. J., Frey, B. J., and Blencowe, B. J. (2008). Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing. Nat. Genet. 40, 1413–1415. doi: 10.1038/ng.259

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, D., Boufraqech, M., Jain, M., Zhang, L., He, M., Gesuwan, K., et al. (2013). MiR-34a and miR-483-5p are candidate serum biomarkers for adrenocortical tumors. Surgery 154, 1224–1228. doi: 10.1016/j.surg.2013.06.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Quidville, V., Alsafadi, S., Goubar, A., Commo, F., Scott, V., Pioche-Durieu, C., et al. (2013). Targeting the deregulated spliceosome core machinery in cancer cells triggers mTOR blockade and autophagy. Cancer Res. 73, 2247–2258. doi: 10.1158/0008-5472.CAN-12-2501

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossner, F., Gieseler, C., Morkel, M., Royer, H. D., Rivera, M., Blaker, H., et al. (2016). Uncoupling of EGFR-RAS signaling and nuclear localization of YBX1 in colorectal cancer. Oncogenesis 5:e187. doi: 10.1038/oncsis.2015.51

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryan, M., Wong, W. C., Brown, R., Akbani, R., Su, X., Broom, B., et al. (2016). TCGASpliceSeq a compendium of alternative mRNA splicing in cancer. Nucleic Acids Res. 44, D1018–D1022. doi: 10.1093/nar/gkv1288

PubMed Abstract | CrossRef Full Text | Google Scholar

Sasatomi, T., Yamana, H., Shichijo, S., Tanaka, S., Okamura, T., Ogata, Y., et al. (2000). Expression of the SART1 tumor-rejection antigens in colorectal cancers. Dis. Colon Rectum. 43, 1754–1758. doi: 10.1007/bf02236863

PubMed Abstract | CrossRef Full Text | Google Scholar

Seiler, M., Peng, S., Agrawal, A. A., Palacino, J., Teng, T., Zhu, P., et al. (2018). Somatic mutational landscape of splicing factor genes and their functional consequences across 33 cancer Types. Cell Rep. 23, 282–296. doi: 10.1016/j.celrep.2018.01.088

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Q., Eun, J. W., Lee, K., Kim, H. S., Yang, H. D., Kim, S. Y., et al. (2018). Barrier to autointegration factor 1, procollagen-lysine, 2-oxoglutarate 5-dioxygenase 3, and splicing factor 3b subunit 4 as early-stage cancer decision markers and drivers of hepatocellular carcinoma. Hepatology 67, 1360–1377. doi: 10.1002/hep.29606

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Q., and Nam, S. W. (2018). SF3B4 as an early-stage diagnostic marker and driver of hepatocellular carcinoma. BMB Rep. 51, 57–58. doi: 10.5483/bmbrep.2018.51.2.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegfried, Z., and Karni, R. (2018). The role of alternative splicing in cancer drug resistance. Curr. Opin. Genet. Dev. 48, 16–21. doi: 10.1016/j.gde.2017.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Szabo, D. R., Luconi, M., Szabo, P. M., Toth, M., Szucs, N., Horanyi, J., et al. (2014). Analysis of circulating microRNAs in adrenocortical tumors. Lab. Invest. 94, 331–339. doi: 10.1038/labinvest.2013.148

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamura, K., Furihata, M., Tsunoda, T., Ashida, S., Takata, R., Obara, W., et al. (2007). Molecular features of hormone-refractory prostate cancer cells by genome-wide gene expression profiles. Cancer Res. 67, 5117–5125. doi: 10.1158/0008-5472.CAN-06-4040

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017). GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 45, W98–W102. doi: 10.1093/nar/gkx247

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Zheng, M., Wang, S., Nie, X., Guo, Q., Gao, L., et al. (2019). Whole genome analysis and prognostic model construction based on alternative splicing events in endometrial cancer. Biomed. Res. Int. 2019:2686875. doi: 10.1155/2019/2686875

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, Y., Deng, Y., Wang, K., Zhou, H., Zheng, X., Si, L., et al. (2018). Profiles of alternative splicing in colorectal cancer and their clinical significance: a study based on large-scale sequencing data. eBio Med. 36, 183–195. doi: 10.1016/j.ebiom.2018.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, C., Wan, R., and Shi, Y. (2019). Molecular mechanisms of pre-mRNA splicing through structural biology of the spliceosome. Cold Spring Harb. Perspect. Biol. 11:a032409. doi: 10.1101/cshperspect.a032409

PubMed Abstract | CrossRef Full Text | Google Scholar

Yutani, S., Shichijo, S., Inoue, Y., Kawagoe, N., Okuda, K., Kurohiji, T., et al. (2001). Expression of the SART1 tumor-rejection antigen in hepatocellular carcinomas. Oncol. Rep. 8, 369–372. doi: 10.3892/or.8.2.369

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Deng, Y., Zuo, Y., Wang, J., and Zhao, Y. (2020). Analysis of colorectal cancer-associated alternative splicing based on transcriptome. DNA Cell Biol. 39, 16–24. doi: 10.1089/dna.2019.5111

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10:1523. doi: 10.1038/s41467-019-09234-9236

CrossRef Full Text | Google Scholar