Comprehensive Analysis of Prognostic Alternative Splicing Signature Reveals Recurrence Predictor for Papillary Thyroid Cancer

Background: Alternative splicing (AS) plays a key role in the diversity of proteins and is closely associated with tumorigenicity. The aim of this study was to systemically analyze RNA alternative splicing (AS) and identify its prognostic value for papillary thyroid cancer (PTC).

Methods: AS percent-splice-in (PSI) data of 430 patients with PTC were downloaded from the TCGA SpliceSeq database. We successfully identified recurrence-free survival (RFS)-associated AS events through univariate Cox regression, LASSO regression and multivariate regression and then constructed different types of prognostic prediction models. Gene function enrichment analysis revealed the relevant signaling pathways involved in RFS-related AS events. Simultaneously, a regulatory network diagram of AS and splicing factors (SFs) was established.

Results: We identified 1397 RFS-related AS events which could be used as the potential prognostic biomarkers for PTC. Based on these RFS-related AS events, we constructed a ten-AS event prognostic prediction signature that could distinguish high-and low-risk patients and was highly capable of predicting PTC patient prognosis. ROC curve analysis revealed the excellent predictive ability of the ten-AS events model, with an area under the curve (AUC) value of 0.889; the highest prediction intensity for one-year RFS was 0.923, indicating that the model could be used as a prognostic biomarker for PTC. In addition, the nomogram constructed by the risk score of the ten-AS model also showed high predictive efficiency for the prognosis of PTC patients. Finally, the constructed SF-AS network diagram revealed the regulatory role of SFs in PTC.

Conclusion: Through the limited analysis, AS events could be regarded as reliable prognostic biomarkers for PTC. The splicing correlation network also provided new insight into the potential molecular mechanisms of PTC.

Introduction

Thyroid cancer is the most rapidly increasing malignancy worldwide in both men and women (1). Papillary thyroid cancer (PTC), the most common thyroid carcinoma type, comprise 80% of all cases (1). PTC has several subtypes, including classical papillary cancers, less aggressive variants, such as follicular, oxyphilic, and cribriform-morular variants, and more aggressive variants, such as diffuse-sclerosing, tall-cell, columnar-cell and solid variants (2). In general, PTC has an excellent prognosis, with a 5-year survival rate over 97%, and PTC tumors measuring less than 1 cm have a 10-year disease-specific survival rate over 99% (3). Although PTC is associated with low mortality, the incidence of disease recurrence or metastasis is 20-30%, and is even higher in patients with the more aggressive variants (4, 5). It is important to assess the PTC recurrence risk accurately for ensuring patients to receive the most appropriate treatment strategy. Over the past few decades, great efforts have been made in exploring prognostic biomarkers for PTC, in particular gene markers, such as mutations in the BRAF, RAS, PIK3CA, P53, PTEN, P53 and ALK genes (6, 7). Although these studies showed promising results, they only focused on the factors driven by mutation and the level of transcription while ignoring the diversity of RNA isoform resulting from posttranscriptional modifications. Currently, there is no consensus on assessing the prognosis of PTC patients.

Alternative splicing (AS) is a significant molecular posttranscriptional modification mechanism that converts mRNA into different RNA transcripts which are then translated into different protein products, thus greatly increasing the diversity of protein species (8, 9). Recent studies have shown that more than 90% of human genes have AS modifications and these modifications play an important role in biological processes (10, 11). The dysregulation of AS is involved in a variety of physiological and pathological processes, including tumorigenesis. Since tumor cells tend to generate sub-isoform changes, which lead to the functional loss of tumor suppressor genes and the activation of oncogenic genes (12). These multiple AS events are conducive to tumor cell proliferation, invasion and metastasis, drug resistance and immune escape (13). For example, exon 13 skipping in CD46 and the exon 13-containing CD46 isoform play opposite roles in bladder cancer development, and exon 13 skipping remarkably accelerated DNA synthesis, cancer cell proliferation, migration and invasion (14). KLF6-SV1, an oncogenic alternatively-spliced isoform of KLF6 produced by alternative 5′ splice sites, is often highly expressed in various human malignancies including non-small cell lung cancer and hepatocellular carcinoma (15–17). BRCA1/2 germline mutations are most commonly seen in breast and ovarian cancer patients who benefit from treatment with PARP inhibitors (PARPis) or platinum compounds, but BRCA1-Δ11q splice variants lacking the majority of exon 11 contribute to therapeutic resistance (18, 19). We found that prognostic models constructed from AS events data had good efficiency for evaluating the survival time of adrenocortical carcinoma, cervical cancer and prostate cancer patients (20–22).

In addition, some studies have shown that AS events could be intricately regulated by key splicing factors (SFs) (23). The abnormal expression of SFs cause subversive alteration in tumor-specific AS events, which affects the initiation and progression of carcinoma (24). In recent years, the development of genome-wide sequencing technology has provided new opportunities to explore and identify tumor-specific molecules and prognostic markers (25, 26). A comprehensive analysis of AS events and underlying SF-AS regulatory networks can provide new insight into the molecular mechanism of PTC and prognosis-related biomarkers for PTC patients. Preliminary studies in AS have provided evidence of prognosis value, while the function and mechanism of AS in PTC remains unknown.

In our study, we revealed a large number of RFS-related AS events in PTC through a systemic analysis of the AS events of all genes in the PTC cohort from the TCGA SpliceSeq dataset. We constructed a prognostic prediction model based on the identification of RFS-associated AS events, and presented the clinicopathological characteristics and a nomogram of AS prognostic predictors, which could predict the recurrence-free survival rate of PTC patients. Finally, development of an SF-AS relationship network diagram showed the potential regulatory mechanisms involved in PTC recurrence and patient prognosis.

Materials And Methods

The flowchart of the study is shown in Figure 1.

FIGURE 1

Figure 1 The flowchart of the study.

Acquisition of AS Data

The percent splice-in (PSI) data of AS events in PTC were downloaded from the TCGA SpliceSeq (https://bioinformatics.mdanderson.org/TCGASpliceSeq/) database, which provides the overview of AS events across 33 types of tumors based on TCGA RNA-seq data. AS events for 7 types have been identified so far, namely 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) (27). The annotation of the AS event consists of the parent gene symbol, the unique ID number and the splicing type. PSI values are used to quantify AS events, and its range is from zero to one. In order to acquire a credible AS events data set, we set a strict screening condition that the proportion of samples contain PSI values over 75%. We filtered out the data with AS events missing rate over than 20%, then replaced the missing value with median values.

Analysis of RFS-Related AS Events, Function and Pathway Enrichment Analysis and Gene Interaction Network

The clinical information of the PTC cohort was obtained from the cBioPortal (http://www.cbioportal.org) database, including recurrence-free survival status and time. Patients were divided into high and low PSI subgroups based on the median value of PSI, and then univariate Cox regression analysis was used to explore RFS-associated seven types AS events respectively, p values less than 0.05 considered as statistically significant. UpSetR mapping was used to analyze the interaction between the RFS-related AS events for each splicing type and corresponding parent genes. Target genes network were constructed via the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org) and Cytoscape (version 3.7.1). Database for Annotation, Visualization and Integrated Discovery (DAVID) online functional annotation tool (https://david.ncifcrf.gov/tools.jsp) was used to complete genetic function and pathway enrichment analysis, and use RStudio drawing.

Construction of the AS Model for Predicting Recurrence of PTC Patients

First, LASSO regression analysis was performed on the RFS-associated AS events obtained of 7 types by univariate Cox regression analysis. In order to avoid overfitting of the model, multivariate Cox regression analysis was used to further screen the candidate AS events and identify independent prognostic predictors. We calculated the risk score for each patient based on each predictor and the calculation formula is as follows: Risk score = PSI_{AS event1} × coefficient _{AS event1}+ PSI_{AS event2} × coefficient _{AS event2}+· · · + PSI_{AS eventn} × coefficient _{AS eventn.} According to the median risk score, PTC patients were divided into high and low risk subgroups, and Kaplan-Meier analysis was used to evaluate the accuracy of each prognostic prediction signature. In addition, the receiver operating characteristic (ROC) curve by the survival ROC package was used to calculate the corresponding area under the curve (AUC) value. Furthermore, the cBioPortal online database was used to analyze mutations and expression changes in corresponding parental genes.

The Verification of Prognostic Value of AS Predictor

The modeling dataset was random divided into two validation datasets (50 percent vs 50 percent, n=215), Kaplan-Meier survival curve and ROC curve were used to evaluate the performance of the model. Besides, we also performed a pan-cancer survival analysis based on data from TCGA.

To further analyze the independent risk factor associated with recurrence of PTC, AS prognostic predictor signature along with all clinicopathological variable mentioned above were performed by univariate Cox regression analysis. The candidate variables were subjected to multivariate regression analysis to screen out independent prognostic predictors.

In addition, we analyzed the clinicopathological characteristics of the high- and low-risk subgroups. Judge and verify the prognostic performance of the final AS prediction model in the stratified survival analysis, such as age, sex, histologic subtype, tumor grade, lymph node grade, and pathological stage.

Analysis of RFS-Related SFs and Construction of SF-AS Relationship Network

Splicing factors (SFs) were obtained from SpliceAid 2 (www.introni.it/spliceaid.html) database. The normalized mRNA expression data of the SFs were obtained from UCSC Xena (https://xena.ucsc.edu) database. The Protein expression level of SFs was obtained from The Human Protein Atlas (https://www.proteinatlas.org/) database. Univariate Cox regression analysis was used to screen out RFS-related SFs. Spearman correlation analysis was used to detect the relationship between RFS-related AS events and SFs, P value less than 0.05 and the correlation coefficient greater than 0.4 as cutoff value. Finally, Cytoscape is used to construct a potential SF-AS relationship network diagram.

Results

A Complete Overview of AS Events in the TCGA PTC Cohort

Through integrating all AS events of PTC patients from the TCGA SpliceSeq database, we discovered 37833 AS events involving 18231 genes, including 10219 ESs in 3904 genes, 9127 APs in 3653 genes, 8597 ATs in 3753 genes, 3683 AAs in 2592 genes, 3190 ADs in 2240 genes, 2787 RIs in 1865 genes, and 232 MEs in 224 genes (see Supplementary Figure 1A). The figures showed that one gene can produce multiple types of AS events in PTC patients. Among these 7 types of AS events, the most frequent splicing type was ES, while the least type was ME (see Supplemental Figure 1A).

Detection of RFS-Related AS Events and Analysis of Function and Pathway Enrichment

The survival and clinical information for PTC was obtained from the cBioPortal database (Supplementary Table 1). There was a total of 430 PTC patients with available recurrence-free survival time data and complete clinical information in our analysis. In PTC cohort, univariate Cox analysis of all AS events revealed that 1396 AS events were significantly related to the RFS (P<0.05, Supplementary Tables 2). In order to better visualize the intersection of different types of AS events and corresponding parent genes, an UpSet plot was constructed, as shown in Supplementary Figure 1B. Interestingly, we found that one gene can produce 3 different types of AS events in this study. The different types of prognoses associated AS events, except ME, in the top 20 genes are clearly exhibited in Figure 2. Next, we performed functional and pathway enrichment of 989 parent genes of RFS-associated AS events. The results showed that a total of 130 GO terms and 3 KEGG terms were significantly involved in prognosis (p<0.05), and Figures 3A–D showed the top 10 GO functional enrichment and KEGG pathways. To further explore the biological association between the corresponding paternal genes in PTC, we used STRING and Cytoscape to create a gene interaction network. Figure 3E shows a network diagram of the parental genes. The larger the node, the greater the degree of association with other genes, and the top 3 genes identified were UBA52, UBB and RPL31, they may be closely related to the occurrence and progression of PTC.

FIGURE 2

Figure 2 The top 20 RFS-associated AS events. (A) Volcano map of AS event, red dots represent RFS- related AS. The top 20 AS events related to recurrence outcomes in different splice types in PTC, including (B) AA, alternate acceptor site. (C) AD, alternate donor site. (D) AP, alternate promoter. (E) AT, alternate terminator. (F) ES, exon skip. (G) ME, mutually exclusive exons. (H) RI, retained intron.

FIGURE 3

Figure 3 The gene interaction network of RFS-associated AS events, functional and pathway analysis. (A) GO biological processes (BP) enrichment. (B) GO cellular component (CC) enrichment. (C) GO molecular function (MF) enrichment. (D) KEGG pathway enrichment analysis. (E) Parent genes interaction network.

Establishment of AS Recurrence Prediction Model for PTC Patients

We performed LASSO regression analysis for the significant RFS-associated AS events in each AS type (Supplementary Figure 2A–G). In order to avoid model overfitting, the above results of each AS type were further analyzed by multivariate Cox regression analysis to screen out the most suitable predictor for AS recurrence models. Seven types of AS models (AA, AT, ME, RI, AD, AP and ES) were constructed, and the formula corresponding to each model was shown in Table 1. Based on the formula, we calculated the risk score of each patient and divided into high and low risk groups. The Kaplan-Meier survival analysis showed that the recurrence model of each AS type had good predictive power to distinguish between good and poor survival results (Figures 4A–G). To further evaluate and compare the efficiency of the model, ROC curves were used to calculate the AUC value predicting the 1-year, 3-year, 5-year and 10-year recurrence-free survival rate (Figures 4a–g). The AUC values ​​of the seven types of models at different times did not exceed 0.04. The largest AUC value of the ROC for the 1-year, 5-year and 10-year RFS rate was obtained with the AA prognostic predictor (0.860, 0.824 and 0.827 respectively), and the largest AUC value for the 3-year survival rate was obtained with the AP model (0.825). Importantly, a ten-AS event predictor was obtained by the overall analysis of prognostic-related AS events using LASSO regression and multivariate Cox regression analysis (Supplementary Figure 2H). The calculation used for the risk score is shown in Tables 1 and  2, the high-risk group showed a worse significant survival outcome than low-risk (Figure 5A). The 1-year, 3-year, 5-year and 10-year AUC values of the ROC curve for the combined model were calculated as 0.923, 0.916, 0.900 and 0.889 respectively. These values were higher than those obtained by the seven separate models individually, suggesting that the mixed AS model had the highest-level performance among all prognostic models (Figure 5B). The distribution of patient survival status and survival time, risk score for the prognostic predictors and the PSI of the ten AS events for final recurrence model, as illustrated in Figure 5C, the results showed that the shorter the patient’s survival time and the more recurrent cases, the higher the risk score of the model was significantly higher (P <0.05, Figure 5C).

TABLE 1

Table 1 Formula of each prognostic signature for PTC.

FIGURE 4

Figure 4 Construction of Kaplan-Meier survival curve and ROC curve and calculation of AUC values ​​for recurrence prognostic predictors. (A–G) Kaplan-Meier survival curve for AA, AD, AP, AT, ES, ME and RI prediction models. (a–g) ROC curve for AA, AD, AP, AT, ES, ME and RI prediction models.

FIGURE 5

Figure 5 The ten-AS events signature as the best prognostic predictive model in PTC cohort. (A) Kaplan-Meier curve of the ten-AS prognostic predictor. (B) ROC curve of the ten-AS prognostic predictor. (C) The upper part showed the heatmap of PSI score of AS events involved in the prognostic predictor. The middle part was risk score of each individual. The bottom was the recurrence status and RFS time of each PTC patients.

TABLE 2

Table 2 Prognostic predictors for PTC.

In addition, parental genetic alteration of the ten-AS event model is shown in Figure 6A. The mutation of these ten genes rarely appeared in PTC patients from the TCGA dataset, but the mRNA expression level of most of the genes were altered; for example, NUD16 expression was decreased in 71% of PTC samples compared to normal tissue (Figure 6A). We detected the relationship between the expression level of parental genes and the RFS rate of PTC patients. There was a statistically significant relationship between PTC patient’s prognosis and the expression of SPHK2, SLC22A17, NUDT16, FXN, ADIRF, MARK3 and MTURN. The representative survival curves showed that NUDT16, MTURN and FXN had the most changes, and high expression was a favorable prognostic factor (Figures 6B–D, Supplementary Figure 3). The changes in mRNA may be caused by AS events, but AS events are not limited to changes in mRNA levels as they are also involved in the specific functions of protein regions.

FIGURE 6

Figure 6 The parent genetic alteration in PTC cases. (A) In the PTC cohort, waterfall plot of parent genes variation and expression from the ten-AS models. (B–D) Kaplan-Meier survival curves of NUDT16, MTURN and FXN gene expression.

The Efficiency of AS Prognostic Predictor in Stratified Clinicopathologic Subgroups in PTC Patients

According to the ten-AS prognostic model, PTC patients were divided into two risk levels (high or low). The clinicopathologic characteristics of the two groups as shown in Table 3, PTC patients with high risk tended to be over the age of 55 and had tumors with a higher grade (p<0.05). Moreover, we analyzed the predictive performance of ten-AS prognostic model in stratified PTC patients (Table 3). The prognostic model identified high-risk patients with worse RFS rates in each subgroup except T1 and Stage II, which may result from the small number of endpoint events in these two subgroups (Table 3). In order to further explore the potential factors related to recurrence in PTC patients, the clinicopathologic variables along with the risk score of the ten-AS prognostic predictor were subjected to univariate Cox regression analysis. Age <55, female sex, histologic type-classical PTC, pathologic T1, pathologic N1 and pathologic stage I were set as references. The results showed that age, tall cell type, T3, T4, N1b, stage III, stage IV and risk score were significantly related to PTC recurrence (Table 4). Furthermore, the meaningful factors were analyzed by multivariate Cox regression analysis, and the results uncovered that the ten-AS prognostic model was the only independent recurrence prognostic factor (Table 4). The risk score of the ten-AS model, the independent predictive factor, was used to establish nomogram (Supplementary Figure 4). Our results suggested that the ten-AS prognostic predictor had a better efficiency in predicting PTC recurrence than clinicopathological characteristics, and predictive value in stratified subgroups well.

TABLE 3

Table 3 Clinicopathology feature of the final AS signature and prognostic analysis in stratified PTC cohorts.

TABLE 4

Table 4 Univariate and multivariate Cox regression analysis for clinicopathology variables.

The Verification of Prognostic Value of the Ten-AS Signature

Internal validation in two datasets shown good performance of the ten-AS signature, PTC patients with high glycolysis scores exhibited worse prognosis (see Supplementary Figure 5). Besides, we also evaluated the prognostic values of the ten-AS prognostic model in various cancers. In our results, the ten-AS prognostic model also applied to prostate adenocarcinoma (PRAD) and lung adenocarcinoma (LUAD), PRAD or LUAD patients with high-risk scores exhibited worse prognosis (see Supplementary Figure 6).

Detection of RFS-Associated SFs and Construction of SF-AS Relationship Network

In order to explore the upstream regulatory factors of dysregulated AS, the expression of 71 SFs was extracted from level 3 RNA-seq data of TCGA PTC. The results for univariate Cox regression analysis exhibited that 5 SFs (KHSRP, NOVA2, PTBP2, SRSF3 and RBM9) were significantly correlated to the RFS rate of PTC patients (Supplementary Table 3). The recurrence-free survival time curve with high and low expression of these 5 SFs shown as Figures 7A–E, among them, KHSRP was oncogenic factor, while NOVA2, PTBP2, SRSF3 and RBM9 were tumor inhibitor. We further searched The Protein Atlas database to detect the protein level of the 5 SFs in PTC. The immunohistochemistry (IHC) results showed that KHSRP, SRSF3 and RBM9 were located in nucleus, and the expression of SRSF3 and RBM9 were significantly lower in cancer than normal thyroid, while there was no significantly different between KHSRP, NOVA2 and PTBP2 in carcinoma and normal tissues (Figure 7F). The mRNA expression level of PTBP2 in normal tissues was higher than PTC, and NOVA2 expression was significantly involved in tumor stage, high expression in I-II stages and low expression in III-IV stage for PTC (Figures 7G, H). In addition, we used Spearman’s test to detect the relationship between the expression of these 5 SFs and PSI values ​​of RFS-associated AS events. The relationship network diagram showed that RFS-related 5 SFs (blue rectangles) were significantly associated with 117 AS events, with P value less than 0.05 and Spearman coefficient greater than or 0.4 as the cutoff value (p<0.05, Spearman≥0.4, Figure 8B). Interestingly, we found that the expression of KHSRP was positively correlated (red lines) with most of adverse survival prognostic AS events (red dots) but negatively correlated (green lines) with most of favorable AS events (green dots), however, the tumor suppressor SFs were inversely related to AS events. For example, Figure 8A exhibited the representative scatter plots of SFs and AS events correlation. Based on our preliminary exploration, we proposed the hypothesis that antineoplastic SFs play a key role in dysregulated AS, which may lead to tumor progression in PTC.

FIGURE 7

Figure 7 Validation the prognostic correlation of the targeting SFs and the expression in PTC tissues. Kaplan-Meier curves of (A) KHSRP (B) SRSF3 (C) RBM9 (D) NOVA2 (E) PTBP2 for PTC. (F) Immunohistochemistry (IHC) staining shown the expression of KHSRP, SRSF3, RBM9, NOVA2 and PTBP2 in normal thyroid tissues and PTC, data obtained from the HUMAN PROTEIN ATLAS database (HPA, https://www.proteinatlas.org/). (G) Comparison of the expression level of KHSRP, SRSF3, RBM9, NOVA2 and PTBP2 in normal thyroid tissues and PTC, data obtained from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) Student's t-test, ***p < 0.001; ns mean no significance. (H) Comparison of the expression level of KHSRP, SRSF3, RBM9, NOVA2 and PTBP2 in different stage (I, II, III and IV stage) of PTC from GEPAI database (http://gepia.cancer-pku.cn/).

FIGURE 8

Figure 8 The correlation between the targeting SFs with AS events in PTC. (A) Two representative scatter plots showed the relevance between the expression of SFs and AS events, (left) positive correlation between NOVA2 and CLSTN1.575.ES PSI value, (right) negative correlation between PTBP2 and EPB41L5.55145.AT PSI value. (B) Interaction network of SF-AS (P<0.05, Spearman>0.4). Blue triangles represent RFS-associated SFs. Red lines represent SFs positively correlated with AS events, and green lines represent SFs negatively correlated with AS events. Red dots represent adverse AS events, and green dots represent favorable AS events.

Discussion

In this study, we first recognized diversified AS events with prognostic power using PSI of PTC AS data obtained from the TCGA. By using TCGA data, reported studies involved to PTC prognostic have shown that long noncoding RNA (lncRNAs), microRNA, and methylation data can work as prognostic factors. Chen et al., using binding motif data from Ensembl Biomart, predicted transcription factors (TFs) for affected genes to construct a TF/lncRNA/mRNA network, which predicted PTC prognosis with an AUC of 0.794 (28). Wang et al. established an N6-methyladenosine (m6A) RNA methylation-related risk signature of disease-free survival for a total PTC cohort with an AUC of 0.817. These models have also shown favorable prognostic predictions (29). We explored all events and established a ten-AS event signature to evaluate the RFS rate of PTC patients. The AUC values of the ROC curve for the 1-year, 3-year, 5-year and 10-year recurrence-free survival rate of PTC patients were 0.923, 0.916, 0.900 and 0.889 respectively, and we obtained more efficient prediction values from this model than with others. Importantly, the present ten-AS event signature has been confirmed with universality in predicting the prognosis of a wide range of tumors, including PRAD and LUAD patients.

AS is a key regulatory factor in the diversity of protein translation and gene phenotype, which is not only involved in normal physiological process but also plays an important role in the occurrence and development of human diseases, including PTC. For example, the alternatively spliced variant of thyroid stimulating hormoneβ (TSHβ), TSHβv (exon 2 deleted, exon 3 retained) has been associated with autoimmune thyroiditis in humans, which is also a high-risk factor for thyroid carcinoma (30). Circadian clock-independent AS events that play an important role in the homeostasis of the endocrine system, such as alternatively spliced Clock and Bmal1, are regulated by thyroid hormone receptor-associated protein 3 (THRAP3) and are closely associated with endocrine diseases including PTC (31, 32). Therefore, we could design different primers to evaluate the presence of AS events and types by PCR experiments and verified them by sanger sequence, which is relatively simple and effective.

In the ten-AS event prognostic prediction signature, some parental genes have been reported to play a key role in oncologic progression. Two spliced variants of MARk3 (exon 16 included and exon 16 skipped) are differentially expressed by neural progenitors and neuronal cells and contribute to the important molecular regulation of cortical development (33). TNFSF13, a tumor necrosis factor, plays a significant role in tumor development and autoimmune diseases, and hypoxia promotes the retention of the intron of TNFSF13 and suppresses the spliced isoform in MCF7 cells, which may contribute to a tumor suppressor effect (34, 35). SEC14L1 with 3 alternatively spliced exons spanning exon 11 was specifically expressed in human peripheral blood leukocytes, and different protein isoforms may show differential expression in breast and ovarian cancer development (36). Nonetheless, few studies have reported the functional characteristic and of other parental genes in this prognostic signature. Moreover, we found that changes in the mRNA levels of the parental genes in most PTC samples were associated with patient prognosis. AS events can affect the level of transcription and proteome expression. Whether the change in mRNA levels is caused by the corresponding AS events needs to be verified with further experiments. However, there were no statistically significant associations between some genes and prognosis, and the loss or gain of regions resulting from AS events might produce meaningful biological behaviors. Therefore, the underlying molecular mechanisms of these AS events in the final model is unclear, and further functional experimental research is necessary.

In addition, we also explored the correlation between clinicopathological characteristics and the RFS rate of PTC patients, and the results of the univariate analysis demonstrated that age greater than 55 years, tall cell variant PTC, T3 and T4, lateral neck lymph node metastasis and pathological stage III and IV are indicative of poor prognosis. However, further multivariate analysis showed that the risk score of the ten-AS model was the only independent prognostic factor of PTC. In addition, we found that the subgroup of high-risk AS signatures was associated with age and tumor stage, which also showed that tumor-specific AS events play a role in cancer progression and metastasis. Moreover, biological function enrichment and pathway analysis of RFS-related AS events showed that cell-cell adhesion and the transforming growth factor beta receptor signaling pathway promote PTC tumor cell growth, invasion and metastasis (37, 38). KEGG enrichment revealed ribosome and transcriptional mis-regulation in cancer that was associated with the tumorigenesis and prognosis of PTC (39). Therefore, we hypothesized that carcinoma-related outcomes due to changes in AS may involve these pathways.

SFs are the main regulators of AS events and influence splicing sites by recognizing and binding precursor mRNAs. In our study, we identified five SFs related to the prognosis of PTC patients. KHSRP was reported to be oncogenic in non-small lung cancer, colorectal cancer and PTC (40–42). Overexpression of KHSRP activated IFN-αJAK-STAT1 signaling pathway and induced lung cancer cell invasion and metastasis (40). KHSRP might be a target mRNA regulated by the STAU1-mediated mRNA decay (SMD) pathway in PTC; however, the detailed mechanism is unclear (41). NOVA2, a key AS regulator of vascular morphogenesis, was overexpressed in lung carcinoma but its expression was negatively correlated with the prognosis of PTC patients in our analysis (43, 44). Finally, an obvious trend for SF-AS correlation network that the most of favorable prognostic AS events were positively correlated with the tumor suppressor SFs, while negatively correlated with oncogenic SFs expression, however, there was opposite relationship between adverse AS events and SFs. The role of 5 SFs in PTC and the regulation of alternative splice events remains to be verified by more experiments. This study provides a deeper understanding of the mechanism of SFs in the regulation and associated splicing patterns, which will help us to further explore the potential mechanism of AS events in the development and progression of PTC.

Although well performance of the present model has been watched, some limitations are inevitably existed in this study. First, this research lacks repeatability data that could be obtained from assessing the established prognostic predictors in other independent cohorts of PTC patients. Second, the prognostic significance of these potential therapeutic targets and diagnostic biomarkers for PTC still needs to be validated with further biological function experiments, mouse model and clinical trial. Nevertheless, our comprehensive analysis of recurrence-related SFs and AS events provides new knowledge and a new perspective for studying intrinsic molecular mechanisms and identifying potential therapeutic strategies for PTC.

Conclusion

We performed a systematic analysis of AS events in PTC and constructed prognostic signatures that can be used to predict the recurrence-free survival rate for PTC patients. The ten-AS event signature involves genes including SPHK2, SEC14L1, SLC22A17, CCL14, NUDT16, FXN, ADIRF, MARK3, TNFSF13 and MTURN, which can affect the prognosis and biological progression of PTC. The identification of prognosis-related AS events and SF regulatory network increases the understanding of the underlying mechanisms of PTC development and provides a new avenue for developing treatment strategies.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author Contributions

Conceptualization, ML, PH, SC, and RK. Resources, RK, ML, PC, HYH, NT, D-jO and BW. Data Curation, ML, RK, PC, HH, D-jO, Y-xZ, and NT. Writing, ML, D-jO and PH. Supervision, D-jO, PH, SC, and RK. Funding Acquisition, PH and SC. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the National Natural Science Foundation of China (81974423, 81902729), the Key Research and Development Programme of Hunan Province (2019SK2031), the Natural Science Foundation of Hunan Province of China (2020JJ5904) and China Postdoctoral Science Foundation (2020M672517).

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.

The handling editor declared a shared affiliation with the authors at time of review.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

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

References

1. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2020. CA Cancer J Clin (2020) 70(1):7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Schneider D, Chen H. New Developments in the Diagnosis and Treatment of Thyroid Cancer. CA Cancer J Clin (2013) 63(6):374–94. doi: 10.3322/caac.21195

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yu X, Wan Y, Sippel RS, Chen H. Should All Papillary Thyroid Microcarcinomas be Aggressively Treated? An Analysis of 18,445 Cases. Ann Surg (2011) 254(4):653–60. doi: 10.1097/SLA.0b013e318230036d

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Brown R, de Souza J, Cohen E. Thyroid Cancer: Burden of Illness and Management of Disease. J Cancer (2011) 2:193–9. doi: 10.7150/jca.2.193

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Mazzaferri E, Jhiang S. Long-Term Impact of Initial Surgical and Medical Therapy on Papillary and Follicular Thyroid Cancer. Am J Med (1994) 97(5):418–28. doi: 10.1016/0002-9343(94)90321-2

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Cell J. Integrated Genomic Characterization of Papillary Thyroid Carcinoma. Cell (2014) 159(3):676–90. doi: 10.1016/j.cell.2014.09.050

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu M, Chen P, Hu HY, Ou-Yang DJ, Khushbu RA, Tan HL, et al. Kinase Gene Fusions: Roles and Therapeutic Value in Progressive and Refractory Papillary Thyroid Cancer. J Cancer Res Clin Oncol (2021) 147(2):323–37. doi: 10.1007/s00432-020-03491-5

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Keren H, Lev-Maor G, Ast G. Alternative Splicing and Evolution: Diversification, Exon Definition and Function. Nat Rev Genet (2010) 11(5):345–55. doi: 10.1038/nrg2776

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Tress M, Abascal F, Valencia A. Alternative Splicing May Not Be the Key to Proteome Complexity. Trends Biochem Sci (2017) 42(2):98–110. doi: 10.1016/j.tibs.2016.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Matera A, Wang Z. A Day in the Life of the Spliceosome. Nat Rev Mol Cell Biol (2014) 15(2):108–21. doi: 10.1038/nrm3742

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kahles A, Lehmann K, Toussaint N, Hüser M, Stark S, Sachsenberg T, et al. Comprehensive Analysis of Alternative Splicing Across Tumors From 8,705 Patients. Cancer Cell (2018) 34(2):211–24.e6. doi: 10.1016/j.ccell.2018.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Dvinge H, Guenthoer J, Porter P, Bradley R. RNA Components of the Spliceosome Regulate Tissue- and Cancer-Specific Alternative Splicing. Genome Res (2019) 29(10):1591–604. doi: 10.1101/gr.246678.118

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bonnal S, López-Oreja I, Valcárcel J. Roles and Mechanisms of Alternative Splicing in Cancer - Implications for Care. Nat Rev Clin Oncol (2020) 17(8):457–74. doi: 10.1038/s41571-020-0350-x

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zeng J, Xu H, Huang C, Sun Y, Xiao H, Yu G, et al. CD46 Splice Variant Enhances Translation of Specific mRNAs Linked to an Aggressive Tumor Cell Phenotype in Bladder Cancer. Mol Ther Nucleic Acids (2021) 24:140–53. doi: 10.1016/j.omtn.2021.02.019

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhang N, Yan Q, Lu L, Shao J, Sun Z. The KLF6 Splice Variant KLF6-SV1 Promotes Proliferation and Invasion of non-Small Cell Lung Cancer by Up-Regultating PI3K-AKT Signaling Pathway. J Cancer (2019) 10(22):5324–31. doi: 10.7150/jca.34212

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Hu K, Zheng QK, Ma RJ, Ma C, Sun ZG, Zhang N. Kruppel-Like Factor 6 Splice Variant 1: An Oncogenic Transcription Factor Involved in the Progression of Multiple Malignant Tumors. Front Cell Dev Biol (2021) 9:661731. doi: 10.3389/fcell.2021.661731

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Lopez-Canovas JL, Del Rio-Moreno M, Garcia-Fernandez H, Jimenez-Vacas JM, Moreno-Montilla MT, Sanchez-Frias ME, et al. Splicing Factor SF3B1 Is Overexpressed and Implicated in the Aggressiveness and Survival of Hepatocellular Carcinoma. Cancer Lett (2021) 496:72–83. doi: 10.1016/j.canlet.2020.10.010

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wang Y, Bernhardy AJ, Cruz C, Krais JJ, Nacson J, Nicolas E, et al. The BRCA1-Delta11q Alternative Splice Isoform Bypasses Germline Mutations and Promotes Therapeutic Resistance to PARP Inhibition and Cisplatin. Cancer Res (2016) 76(9):2778–90. doi: 10.1158/0008-5472.CAN-16-0186

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Gelli E, Colombo M, Pinto AM, De Vecchi G, Foglia C, Amitrano S, et al. Usefulness and Limitations of Comprehensive Characterization of mRNA Splicing Profiles in the Definition of the Clinical Relevance of BRCA1/2 Variants of Uncertain Significance. Cancers (Basel) (2019) 11(3):295. doi: 10.3390/cancers11030295

CrossRef Full Text | Google Scholar

20. Liang W, Sun F. Prognostic Alternative mRNA Splicing in Adrenocortical Carcinoma. Front Endocrinol (Lausanne) (2021) 12:538364. doi: 10.3389/fendo.2021.538364

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Zhao J, Chang L, Gu X, Liu J, Sun B, Wei X. Systematic Profiling of Alternative Splicing Signature Reveals Prognostic Predictor for Prostate Cancer. Cancer Sci (2020) 111(8):3020–31. doi: 10.1111/cas.14525

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ouyang D, Yang P, Cai J, Sun S, Wang Z. Comprehensive Analysis of Prognostic Alternative Splicing Signature in Cervical Cancer. Cancer Cell Int (2020) 20:221. doi: 10.1186/s12935-020-01299-4

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Tripathi V, Ellis J, Shen Z, Song D, Pan Q, Watt A, et al. The Nuclear-Retained Noncoding RNA MALAT1 Regulates Alternative Splicing by Modulating SR Splicing Factor Phosphorylation. Mol Cell (2010) 39(6):925–38. doi: 10.1016/j.molcel.2010.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Ladomery M. Aberrant Alternative Splicing Is Another Hallmark of Cancer. Int J Cell Biol (2013) 2013:463786. doi: 10.1155/2013/463786

PubMed Abstract | CrossRef Full Text | Google Scholar

25. de Klerk E, t Hoen P. Alternative mRNA Transcription, Processing, and Translation: Insights From RNA Sequencing. Trends Genet (2015) 31(3):128–39. doi: 10.1016/j.tig.2015.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Katz Y, Wang E, Airoldi E, Burge C. Analysis and Design of RNA Sequencing Experiments for Identifying Isoform Regulation. Nat Methods (2010) 7(12):1009–15. doi: 10.1038/nmeth.1528

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ryan M, Wong W, Brown R, Akbani R, Su X, Broom B, et al. TCGASpliceSeq a Compendium of Alternative mRNA Splicing in Cancer. Nucleic Acids Res (2016) 44:D1018–22. doi: 10.1093/nar/gkv1288

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Chen Y, Jiang B, Wang W, Su D, Xia F, Li X. Identifying the Transcriptional Regulatory Network Associated With Extrathyroidal Extension in Papillary Thyroid Carcinoma by Comprehensive Bioinformatics Analysis. Front Genet (2020) 11:453. doi: 10.3389/fgene.2020.00453

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Wang X, Fu X, Zhang J, Xiong C, Zhang S, Lv Y. Identification and Validation of M(6)A RNA Methylation Regulators With Clinical Prognostic Value in Papillary Thyroid Cancer. Cancer Cell Int (2020) 20:203. doi: 10.1186/s12935-020-01283-y

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Klein J. Novel Splicing of Immune System Thyroid Stimulating Hormone β-Subunit-Genetic Regulation and Biological Importance. Front Endocrinol (Lausanne) (2019) 10:44. doi: 10.3389/fendo.2019.00044

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang D, Jones R, James P, Kitahara C, Xiao Q. Associations Between Artificial Light at Night and Risk for Thyroid Cancer: A Large US Cohort Study. Cancer (2021) 127(9):1448–58. doi: 10.1002/cncr.33392

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Marcheva B, Perelis M, Weidemann BJ, Taguchi A, Lin H, Omura C, et al. A Role for Alternative Splicing in Circadian Control of Exocytosis and Glucose Homeostasis. Genes Dev (2020) 34(15-16):1089–105. doi: 10.1101/gad.338178.120

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Liu J, Geng A, Wu X, Lin R, Lu Q. Alternative RNA Splicing Associated With Mammalian Neuronal Differentiation. Cerebral Cortex (2018) 28(8):2810–6. doi: 10.1093/cercor/bhx160

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Xiao Y, Motomura S, Podack E. APRIL (TNFSF13) Regulates Collagen-Induced Arthritis, IL-17 Production and Th2 Response. Eur J Immunol (2008) 38(12):3450–8. doi: 10.1002/eji.200838640

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Han J, Li J, Ho JC, Chia GS, Kato H, Jha S, et al. Hypoxia Is a Key Driver of Alternative Splicing in Human Breast Cancer Cells. Sci Rep (2017) 7(1):4108. doi: 10.1038/s41598-017-04333-0

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kalikin LM, Bugeaud EM, Palmbos PL, Lyons RH Jr, Petty EM. Genomic Characterization of Human SEC14L1 Splice Variants Within a 17q25 Candidate Tumor Suppressor Gene Region and Identification of an Unrelated Embedded Expressed Sequence Tag. Mamm Genome (2001) 12(12):925–9. doi: 10.1007/s00335-001-2073-3

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Vasko V, Espinosa A, Scouten W, He H, Auer H, Liyanarachchi S, et al. Gene Expression and Functional Evidence of Epithelial-to-Mesenchymal Transition in Papillary Thyroid Carcinoma Invasion. Proc Natl Acad Sci U S A (2007) 104(8):2803–8. doi: 10.1073/pnas.0610733104

PubMed Abstract | CrossRef Full Text | Google Scholar

38. He J, Jin Y, Zhou M, Li X, Chen W, Wang Y, et al. Solute Carrier Family 35 Member F2 Is Indispensable for Papillary Thyroid Carcinoma Progression Through Activation of Transforming Growth Factor-β Type I Receptor/Apoptosis Signal-Regulating Kinase 1/Mitogen-Activated Protein Kinase Signaling Axis. Cancer Sci (2018) 109(3):642–55. doi: 10.1111/cas.13478

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Jeong S, Kim I, Kim H, Choi M, Lee J, Jo YJE, et al. Liver X Receptor β Related to Tumor Progression and Ribosome Gene Expression in Papillary Thyroid Cancer. Endocrinol Metab (Seoul) (2020) 35(3):656–68. doi: 10.3803/EnM.2020.667

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Yan M, Sun L, Li J, Yu H, Lin H, Yu T, et al. RNA-Binding Protein KHSRP Promotes Tumor Growth and Metastasis in non-Small Cell Lung Cancer. J Exp Clin Cancer Res (2019) 38(1):478. doi: 10.1186/s13046-019-1479-2

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Gou Q, Gao L, Nie X, Pu W, Zhu J, Wang Y, et al. Long Noncoding RNA AB074169 Inhibits Cell Proliferation via Modulation of KHSRP-Mediated CDKN1a Expression in Papillary Thyroid Carcinoma. Cancer Res (2018) 78(15):4163–74. doi: 10.1158/0008-5472.CAN-17-3766

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Caiazza F, Oficjalska K, Tosetto M, Phelan J, Noonan S, Martin P, et al. KH-Type Splicing Regulatory Protein Controls Colorectal Cancer Cell Growth and Modulates the Tumor Microenvironment. Am J Pathol (2019) 189(10):1916–32. doi: 10.1016/j.ajpath.2019.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Angiolini F, Belloni E, Giordano M, Campioni M, Forneris F, Paronetto M, et al. A Novel L1CAM Isoform With Angiogenic Activity Generated by NOVA2-Mediated Alternative Splicing. eLife (2019) 8:e44305. doi: 10.7554/eLife.44305

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Li C, Liu H, Niu Q, Gao JJC. Circ_0000376, a Novel circRNA, Promotes the Progression of Non-Small Cell Lung Cancer Through Regulating the miR-1182/NOVA2 Network. Cancer Manag Res (2020) 12:7635–47. doi: 10.2147/cmar.S258340

PubMed Abstract | CrossRef Full Text | Google Scholar