Systematic Genome-Wide Profiles Reveal Alternative Splicing Landscape and Implications of Splicing Regulator DExD-Box Helicase 21 in Aggressive Progression of Adrenocortical Carcinoma

Alternative splicing (AS) in the tumor biological process has provided a novel perspective on carcinogenesis. However, the clinical significance of individual AS patterns of adrenocortical carcinoma (ACC) has been underestimated, and in-depth investigations are lacking. We selected 76 ACC samples from the Cancer Genome Atlas (TCGA) SpliceSeq and SpliceAid2 databases, and 39 ACC samples from Fudan University Shanghai Cancer Center (FUSCC). Prognosis-related AS events (PASEs) and survival analysis were evaluated based on prediction models constructed by machine-learning algorithm. In total, 23,984 AS events and 3,614 PASEs were detected in the patients with ACC. The predicted risk score of each patient suggested that eight PASEs groups were significantly correlated with the clinical outcomes of these patients (p < 0.001). Prognostic models produced AUC values of 0.907 in all PASEs’ groups. Eight splicing factors (SFs), including BAG2, CXorf56, DExD-Box Helicase 21 (DDX21), HSPB1, MBNL3, MSI1, RBMXL2, and SEC31B, were identified in regulatory networks of ACC. DDX21 was identified and validated as a novel clinical promoter and therapeutic target in 115 patients with ACC from TCGA and FUSCC cohorts. In conclusion, the strict standards used in this study ensured the systematic discovery of profiles of AS events using genome-wide cohorts. Our findings contribute to a comprehensive understanding of the landscape and underlying mechanism of AS, providing valuable insights into the potential usages of DDX21 for predicting prognosis for patients with ACC.


Introduction
Adrenocortical carcinoma (ACC) is a rare endocrine malignancy with an incidence of 0.7-2 new cases per million per year, and it accounts for 0.2% of all cancer mortalities in the United States (Kostiainen et al. 2019). Surgery is usually the first and most effective therapeutic strategy for treating patients with ACC. However, early diagnosis of ACC is relatively difficult, so most patients have metastasized tumors at initial diagnosis, resulting in few chances of surgical intervention and poor prognosis (Xu et al. 2019a,b,c;Jonasch et al. 2021). Currently, mitotane is the major medication approved by the United States Food and Drug Administration for ACC treatment for patients with metastasis (Schteingart et al. 2005). High-risk patients with ACC can be treated with mitotane alone or in combination with cytotoxic drugs, but these treatments have limited survival benefits (Jasim and Habra 2019). For example, the combination of mitotane and etoposide/doxorubicin/cisplatin (EDP) is now the standard treatment in advanced ACC, but its effects in improving survival time is still unsatisfactory (Libe 2015;Tierney et al. 2019).
Alternative splicing (AS) is a post-transcriptional process, in which different exonic regions of pre-mRNAs are spliced together and intronic regions are removed, thus producing different mRNA transcripts from the same gene (Maniatis and Tasic 2002;Urbanski et al. 2018;Xu et al. 2021). High-throughput sequencing technology and bioinformatics studies have found that about 35-60% of the entire human genome undergoes AS variants, and about 50% of human genes have multiple transcripts (Nilsen and Graveley 2010). There exist seven types of AS events: alternate acceptor site (AA), alternate donor site (AD), alternate promotor (AP), alternate terminator (AT), exon skip (ES), mutually exclusive exons (ME), and retained intron (RI) (Climente-Gonzalez et al. 2017). Recent studies showed that AS is a universal mechanism of functional regulation in the human genome, and that transcript variants of the same genes may play opposing roles in cancer progression (Lee andRio 2015, Climente-Gonzalez et al. 2017). Interestingly, the mRNAs of splicing factors (SFs) are frequently mutated and abnormally expressed in different cancers where they shift gene expression landscapes and contribute to oncogenesis (Rahman et al. 2020). However, the functions and targets of SFs during cancer progression and microenvironment remain elusive. Hence, it is critical to investigate specific SFs that directly regulate mRNA splicing-related pathways and lead to multiple oncogenic splicing changes, which may be promising candidate targets for cancer treatment.
Comprehensive management of AS alterations and SFs with clinical profiles has recently contributed to great breakthroughs in many cancers Meng et al. 2019). In this study, we aimed to identify prognosis-related AS events (PASEs) using the TCGA SpliceSeq dataset and clinicopathological data of patients with ACC by complex bioinformatics and real-world data. We identified specific SFs and performed functional enrichment and correlation analyses to determine their possible biological roles.

Acquisition, Normalization, and Processing of Raw Data Based on Multiple Databases
The RNA sequencing profiles and matched clinicopathological data of 76 ACC patients were obtained from TCGA database (Tomczak et al. 2015). University of North Carolina TCGA genome characterization center experimentally measured profiles of genes by Illumina HiSeq 2000 RNA Sequencing platform. In addition, the mRNA alternative splicing data of 76 ACC patients was downloaded from TCGA SpliceSeq dataset (https:// bioin forma tics. mdand erson. org/ TCGAS plice Seq/) (Ryan et al. 2012). The percent-spliced-in (PSI) value is a common, intuitive ratio for quantifying splicing events (Ryan et al. 2016). PSI value was defined as the ratio of inclusion/exclusion normalized read counts as a percentage of the total (both inclusion and exclusion) normalized read counts for that event, and it was calculated for different types of AS events . The PSI value can be estimated as the ratio of the density of inclusion reads (i.e., reads per position in regions supporting the inclusion isoform) to the sum of the densities of inclusion and exclusion reads (Wang et al. 2008). For each splicing event, a PSI value was used for quality control according to an intuitive ratio of the long form on total form presented, ranging from 0 to 1. In this study, PSI value was identified as follows: the average PSI value ≥ 0.05 and the minimum PSI standard deviation ≥ 0.01.

Identification and Matching of PASEs
Prognosis-related AS events were selected on the basis of PSI values and overall survival (OS) of ACC patients by univariate Cox regression analysis. The "UpSet" package in R software was used for quantitative intersective analysis of prognosis-related AA, AD, AP, AT, ES, ME, and RI events. Then, the top 20 significant events in seven different types of AS events were presented in bubble charts using R software. To explore critical PASEs and avoid overfitting and bias, "glmnet" package in R software was used for LASSO regression analysis among AA, AD, AP, AT, ES, ME, RI, and all AS events. A multivariate Cox regression analysis was used to calculate the risk score for each ACC patient on the basis of the corresponding PSI value of highly important prognostic AS ).

Survival Analysis and ROC Construction of PASEs
Clinical and pathological parameters, including age, gender, TNM stage, and follow-up data, of ACC patients were processed from TCGA database. A total of 76 ACC patients were divided into low-risk and high-risk groups based on the risk score. Log-rank test in separate curves and Kaplan-Meier method with 95% confidence intervals (95% CI) were utilized to perform the follow-up duration analysis using "Survival" package in R software. To further investigate significant independent clinicopathological factors, including age (ref. Low) of ACC, univariate and multivariate Cox regression analyses were performed using "Survival" package of eight groups. The receiveroperating characteristic curve (ROC) was constructed in AA, AD, AP, AT, ES, ME, RI, and all types of AS events using the "ROC" package in R software. The area under the curve (AUC) was calculated to assess the sensitivity and specificity value of all types of AS events for each prediction model.

Development of Interaction Networks of Splicing Factors
The splicing factors (SFs) that interact with sequence elements on RNA precursors were obtained from the SpliceAid2 database (www. intro ni. it/ splic eaid. html) (Giulietti et al. 2013). The correlation network between the transcriptional expression of SFs and PSI values of survival-related AS events was measured using Spearman's correlation coefficient, and it was visualized using Cytoscape software (Version 3.6.1). The threshold value of Spearman's R was set more than 0.8, and p value was adjusted less than 0.001. Protein-protein interaction (PPI) network of SFs and neighbor genes was constructed by GeneMANIA (http:// www. genem ania. org), which was utilized for generating hypotheses about gene functions (Zuberi et al. 2013).

Functional Annotations of Prognosis-Related AS and SFs
Function enrichments of SFs in Gene Ontology (GO) enrichment analysis, including biological process (BP), cellular components (CC), and molecular functions (MF), were annotated using bubble charts. Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathways enrichment analyses were also performed in selected SFs. Network of the most significant neighbor signaling pathways related to mRNA splicing was developed with correlation score more than 0.4.

Correlation Analysis of PASEs and SFs
The correlation among eight groups of PASEs and SFs was displayed in correlation heat map using "Corrplot" package in R software. Hierarchical clustering of identified SFs was constructed in the heatmap on the basis of the mRNA expression of SFs in 76 ACC patients.

Immunohistochemical (IHC) Staining and Evaluation
To further improve the clinical value of the study, a total of 39 ACC patients, who underwent surgery in Fudan University Shanghai Cancer Center (FUSCC; Shanghai, China) from April 2015 to May 2019, were enrolled in this study. Tissue samples were collected during surgery and available from FUSCC tissue bank. IHC staining of DExD-Box Helicase 21 (DDX21) was performed using a mouse monoclonal anti-DDX21 antibody [1:500, ab182156, Abcam, USA] in 39 ACC samples. Positive or negative staining of DDX21 protein in an FFPE slide was independently evaluated. The overall IHC score ranging from 0 to 12 was measured based on the multiply of the staining intensity and extent score, as previously described (Xu et al. 2019a,b,c). Low DDX21 expression group scores are from 0 to 3, and high DDX21 group scores are from 4 to 12.

Survival and Functional Enrichment Analysis of DDX21
Differential DDX21 expression levels and its clinical implications were performed based on TCGA cohort. Gene set enrichment analysis (GSEA) was used to explore the underlying involved signaling hallmarks of ACC microenvironment. Then, TCGA Splicing Variants database was used to identify splicing locations on exon and junction of DDX21 (Sun et al. 2018).

Cell Culture and Transfection
The human ACC cell lines (SW-13 and NCI-H295R) were obtained from the American Type Culture Collection (Rockville, MD). SW-13 and NCI-H295R cells were cultured in Dulbecco's modified Eagle's medium (DMEM) supplied with 10% fetal bovine serum (Gibco, US) and 100 U/mL penicillin in the atmosphere filled with 5% CO 2 at 37 Celsius. DDX21-specific shRNA (5′-CGC TCC TTG ATC AAC TCA AAT-3′) and negative control shRNA (5′-GCA CTA CCA GAG CTA ACT CAG ATA GTA CT-3′) were synthesized using Lipofectamine 8000 reagent (Invitrogen) according to the manufacturer's instructions. Cells were harvested for further analysis at 24 h after transfection.

Real-Time Quantitative PCR (RT-qPCR) Analysis
We performed RT-qPCR to measure total mRNA levels of splicing factor DDX21 and targeted genes, including CNOT4, KDM51, NFATc3, PCBP2, Tcf3, SPRY1 in human ACC cell lines, SW-13, and NCI-H295R. The primer sequences are shown in Supplementary Table S1. Primescript RT reagent kit (Takara, Otsu, Japan) and SYBR Premix Ex Taq Kit (Takara) were used (Xu et al. 2019a,b,c;Wang et al. 2020). GAPDH was examined as the internal control. The 2 −ΔΔCt method was applied to perform statistical analysis. Each experiment was repeated at least three times.

Statistical Analysis
In this study, R (Version 3.3.2) and RStudio (Version 1.2) were utilized to perform most of data analyses, including Upset plot, Cox regression analyses, LASSO regression analysis, Kaplan-Meier plots, ROC curves, risk plots, PPI network and functional annotations. All tests were two-sided, and p value less than 0.05 was taken as stastically significant.

Results
This study was performed in four stages with a schematic flowchart of systematic profiling shown in Fig. 1. First, 3,614 PASEs and 2,261 corresponding genes were identified after matching SpliceSeq sequences and clinical pathological data in TCGA-ACC patients. Second, survival analysis and ROCs were developed after LASSO regression analysis to avoid overfitting. Third, significant SF networks were constructed, and functional annotations were performed. Fourth, identification, validation, and prognostic implications of differential expressed DDX21 were performed.

Overview and Integration of AS Patterns in ACC
A schematic diagram of seven different types of AS patterns is displayed in Fig. 2A. In the ACC cohort, 23,984 AS events in 14,074 genes were found from 76 ACC samples after quality control, including 2,013 AA events in 1,560 genes, 1,719 AD events in 1,313 genes, 4,257 AP events in 2,427 genes, 4,909 AT events in 2,819 genes, 9,201 ES events in 4,624 genes, 74 ME events in 72 genes, and 1,811 RI events in 1,259 genes (Fig. 2B). These results indicated that ES pattern was the major type and consisted of approximately 58% of all AS events. In addition, the Upset intersection diagram suggested that there might be several AS event types in most genes (Fig. 2C).

Screening and Identification of PASEs
After univariate survival analysis, a total of 3,614 PASEs and 2,261 genes were collected from 23,984 AS patterns, as shown in an Upset intersection diagram (Fig. 2D). Next, the volcano plot indicated that PASEs and insignificant AS events were screened and identified (Fig. 2E). The top 20 significant PASEs in AA, AD, AP, AT, ES, ME, RI, and all types of AS events were displayed in bubble charts ( Fig. 2F-L). LASSO regression model was conducted in each 7 PASEs and all types of AS events, respectively (Supplementary Fig. 1). A total of 13 AA events, 13 AD events, 9 AP events, 8 AT events, 11 ES events, 11 ME events, 9 RI events, and 11 all AS events were screened and identified as the critical prognosis-related AS patterns.

Survival Analysis and ROC Construction of PASEs
To measure prognostic implications of each AS-related risk score, Kaplan-Meier analysis of AA, AD, AP, AT, ES, ME, RI, and all types of AS event groups suggested that eight groups were significantly correlated with prognosis of ACC patients (p < 0.001, Fig. 3A-H) with highrisk group in red and low-risk group in blue. ROC analysis indicated the AUC value of each group as follows: AUC for AA = 0.927, AD = 0.992, AP = 0.859, AT = 0.797, ES = 0.921, ME = 0.817, RI = 0.944, All = 0.907. All ROCs showed strong efficiency as prognostic models (Fig. 3I).

Survival Risk Assessment and Cox Regression Analysis
Survival risk assessment identified patients with high-risk score, and it suggested markedly relationship between the distribution of risk score and survival time in AA, AD, AP, AT, ES, ME, RI, and all of AS event groups. Meanwhile, BLOC1S1-AP, TRAFD1-AD, CIRBP-ES, and USMG5-ES were identified as significant differential AS events in all types of groups. Significant differential AS predictors of AA, AD, AP, AT, ES, ME, and RI are shown in Supplementary Fig. 2.
In addition, univariate and multivariate Cox analyses models of prognosis-related AS patterns of eight groups were constructed using forest plots (Fig. 4). Univariate Cox analyses showed that age (ref.

Interaction and PPI Networks of SFs
A total of eight SFs, including BAG2, CXorf56, DDX21, HSPB1, MBNL3, MSI1, RBMXL2 and SEC31B, were identified and matched as significant splicing factors in close relation to regulation of AS patterns in ACC. Interaction network between prognosis-related AS patterns and eight SFs is displayed in Fig. 5A. The significant nodes in PPI, activation, phosphorylation, or inhibition relationship with eight SFs were shown in circus plot (Fig. 5B). To detect whether the selected SFs could regulate predicted targeted genes, DDX21-specific shRNA was synthesized

Functional and Pathway Enrichment Analysis of SFs
GO enrichments analysis suggested that alterations of SFs in BP were mainly enriched in cargo loading into vesicle, cellular response to VEGF stimulus, mRNA processing and RNA splicing; alterations of SFs in CC were significantly enriched in chaperone complex, endoplasmic reticulum exit site, polysome and peptidase complex; alterations of SFs in MF were significantly enriched in protein binding involved in protein folding, regulation of snoRNA binding and rRNA binding  5D). Pathway analysis suggested that alterations in KEGG are mainly enriched in protein processing in endoplasmic reticulum, VEGF signaling pathway, mRNA surveillance pathway, amoebiasis and spliceosome; alterations in Reactome pathway enrichment analysis existed in AUF1 that bind and destabilize mRNA, regulation of HSF1-mediated heat shock response, regulation of mRNA stability by proteins that bind AU-rich elements, MAPK6/MAPK4 signaling, VEGFA-VEGFR2 pathway, positive epigenetic regulation of rRNA expression, and signaling by VEGF (Fig. 5E). The correlation specificity network between GO: RNA splicing pathway, a significant pathway in BP enrichment in blue, and other top 20 related enrichment biological processes in yellow were illustrated (Fig. 5F). Absolute correlation coefficient was greater than 0.04, and p value was less than 0.01. The gradation of red lines between these nodes represents different degrees of correlation.

Identification and Validation of Prognostic Value of DDX21 in ACC Patients
As displayed in Fig. 6A, unsupervised clustering correlation profiles of eight SFs and different types of PASEs based on PSI values were shown in a heat map. DDX21 and BAG2 showed the most significant associations with AS patterns with correlation value (Correlation value > 0.2). For example, DDX21 showed relatively strong correlation with PSI values for all types of AS events, especially in PSI_All group (r = 0.25). Next, elevated DDX21 expression was found significantly associated with shorter OS (p < 0.0001, HR 4.6, Fig. 6B) and DFS (p < 0.0001, HR 4.3, Fig. 6C). The expression of BAG2 was decreased in tumor compared with normal tissues, while no statistically significant association between BAG2 expression and prognosis was found in 76 ACC patients from TCGA cohort ( Supplementary Fig. 3A, B). In addition, DDX21 mRNA expression showed significantly positive relationship with advanced clinical stage (Spearman rho = 0.283, p = 0.0127), and peak DDX21 expression was found in stage IV (Fig. 6D). Immunohistochemical results suggested that DDX21 expression was found significantly highly expressed in ACC samples compared with normal adrenal gland samples (Fig. 6E). To investigate potential involved hallmarks, GSEA analysis showed a total of 100 up-and down-regulation genes (b), and it indicated that DDX21 was significantly involved in mitotic spindle (NES: p= 0.183, NOM: p = 0.018, Fig. 6G). Meanwhile, 100 genes with positive and negative correlation were plotted according to differential DDX21 expression groups (Fig. 6H).   BAG2, CXorf56, DDX21, HSPB1, MBNL3, MSI1, RBMXL2 and SEC31B, were identified and matched as significant splicing factors in close relation to regulation of AS patterns in ACC. B PPI network circus plot. C The expression levels of DDX21 and the targeted genes (CNOT4, KDM51, NFATc3, PCBP2, Tcf3 and SPRY1) in sh-DDX21 group compared with normal con-trol. D GO enrichments analysis of hub SFs. E KEGG and Reactome pathway analysis of hub SFs. F The correlation specificity network between GO:RNA splicing pathway, a significant pathway in BP enrichment, is in blue, and other top 20 related enrichment biological processes are in yellow. Absolute correlation coefficient was greater than 0.04, and p value was less than 0.01. The gradation of red lines between these nodes represents different degrees of correlation Furthermore, the association between the amplification of AS events frequency and DDX21 expression level according to different clinical stages is shown in Fig. 6I.
Interestingly, ACC could be divided into different subtypes due to its heterogeneous nature. To explore the description of the variable AS events in the molecular typing of ACC samples and the impact on the prognosis, this study classified the samples according to the cluster of clusters (CoCs) from four platforms (DNA copy number; mRNA expression; DNA methylation; miRNA expression) in TCGA-ACC cohort. It is suggested that the riskscore constructed using all significant AS event and DDX21 mRNA Fig. 6 Identification and validation of prognostic value of DDX21 in ACC patients. A Unsupervised clustering correlation profiles of eight SFs and different types of PASEs based on PSI values. B, C Prognostic value of DDX21 expression in 76 ACC patients from TCGA cohort. D Relationship between DDX21 mRNA expression and advanced clinical stage, and the peak of DDX21 expression was found in stage IV. E Immunohistochemical DDX21 expression in ACC samples compared with normal adrenal gland samples from FUSCC cohort. F GSEA analysis of 100 up-and down-regulation genes. G DDX21 was significantly involved in mitotic spindle. H 100 genes with positive and negative correlation were plotted according to differential DDX21 expression groups using GESA methods in the hallmarks of mitotic spindle. I Association between the amplification of AS events frequency and DDX21 expression level according to different clinical stages expression significantly increased with elevated CoCs level using ANOVA test (p < 0.05, Supplementary Fig. 4).

Prognostic Validation of DDX21 in ACC Patients from FUSCC Cohort
To further improve the clinical value of the study, we supplemented the IHC experiment and followed up, and both the survival and clinical data of 39 ACC patients matching the specimens were updated. The findings suggested that DDX21 expression was significantly correlated with aggressive progression and poor prognosis for ACC patients (OS: p = 0.005, HR 4.171; PFS: p = 0.020, HR 2.781; Fig. 7A, B). In ACC samples with identified pathological necrosis, DDX21 expression was markedly associated with aggressive progression and poor prognosis for ACC patients (OS: p = 0.011, HR = 8.556; PFS: p = 0.020, HR = 4.748; Fig. 7C, D). These real-world findings suggest that hub SFs, such as DDX21, could play vital roles in the development of ACC and serve as a potential clinical biomarker.

Discussion
Traditional modeling methods based on genomic characterization have limited clinical prediction abilities for patients with ACC, who have high mortality rates and complex treatment responses (Xu et al. 2019a,b,c). The role of AS in many biological processes has provided novel perspectives and a better understanding of diverse post-transcriptional regulation in diseases such as cancers (Rahman et al. 2020). However, the clinical implications of the AS pattern in individual patients with ACC are unknown, which requires an in-depth investigation. In this study, we comprehensively profiled the AS events in a powerful large-scale genome-wide cohort to elucidate the AS landscape in ACC.
Over the decades, the role of abnormal AS events has been considered as biomarkers and treatment targets in Fig. 7 Prognostic validation of DDX21 in ACC patients from FUSCC cohort. A, B Survival outcomes clarified by DDX21 expression detected using IHC analysis for 39 patients with ACC from FUSCC cohort. C, D Survival outcomes clarified by DDX21 expression detected using IHC analysis for 25 ACC patients with pathological necrosis from FUSCC cohorts diverse cancers, while the lack of current sequencing technologies and analysis processes hinders systematic analysis. Recently, increasing evidence suggested that RNA-seq is a reliable sequencing technique to be used in alternative splicing studies (Feng et al. 2013;Consortium 2014). Accumulating evidence also indicated that different types of AS events identified using RNA-Seq data play crucial roles in carcinogenesis of many neoplasms (Liu et al. 2018;Meng et al. 2019). Here, we first constructed a comprehensive landscape of AS alterations in ACC patients and identified PASEs and significant SFs on the basis of multiply platforms using RNA-seq data. A total of 3,614 PASEs and 2,261 corresponding genes were identified after matching SpliceSeq and clinical pathological data in TCGA-ACC patients. Notably, survival analyses in AA, AD, AP, AT, ES, ME, RI, and all types of AS event groups revealed that the higher risk scores were significant independent parameters for worse overall survival in ACC patients, which further confirms the stability of RNA-Seq technique to explore AS alterations in cancers.
To our knowledge, genome-wide analyses of tissue-specific AS events and SF expression profiles in ACC dataset have not been reported yet (Bates 2020). Hence, we integrated PSI-specific AS events and significant SF networks to investigate the splicing pathway mechanism. We found that alterations in SFs were markedly enriched in the mRNA splicing process. Interestingly, the vascular endothelial growth factor (VEGF) gene, which encodes a protein involved in angiogenesis, undergoes AS events that yield four different transcripts, and thus it is an important factor in the cascade that leads to the onset of carcinogenesis and metastasis .
The clinical implications of the tight connection between abnormal regulation of SFs and novel aspects of cancer biology have attracted considerable interest. In the present study, we conducted complex bioinformatic analyses to comprehensively elucidate the landscape of alternative SFs in ACC. Among the prognosis-related AS patterns, eight genes (BAG2, CXorf56, DDX21, HSPB1, MBNL3, MSI1, RBMXL2 and SEC31B) were identified as essential regulatory elements in splicing pathways and other biological processes. For example, DDX21 (RNA helicase II/Gu) is a crucial regulatory sensor for DNA polymerase I and II that modulates genome dynamics and safeguards genome integrity (Valdez et al. 2002). In rRNA processing, DDX21 unwinds double-stranded RNA (helicase) and enables the formation of secondary structures in single-stranded RNA (Treiber et al. 2017;Xing et al. 2017). In addition, elevated DDX21 expression has been found in many cancers and is associated with malignant progression (Cao et al. 2018;Zhang et al. 2018;Kim et al. 2019). We identified the differential expression of DDX21 in ACC samples compared with the controls and validated its prognostic value for patients with ACC. Comprehensively decoding the functions and potential mechanism of the PASEs landscape and significant SFs may help with the discovery of novel aspects of carcinogenesis that may have implications for therapeutic strategies. Therefore, these studies provide a structural basis for the RNA recognition and unwind the mechanism of DDX21 with a novel perspective for the virus-host interaction interface, which is conducive to the development of drugs targeting DDX21, and they brings more treatment options for patients with high-grade adrenocortical carcinoma.
This study has several limitations. First, in May 2020, ) reported survival-related alternative SFs based only on TCGA SpliceSeq data without external validation. Although there were some similarities in the initial data processing between their study and the present study, we used a more intuitive research flowchart and more accurate data processing methods, and in-depth investigations of hub SFs in the ACC samples from FUSCC were conducted. Second, our study had limited clinical relevance. To improve the clinical translational value, we performed a survival risk assessment and Cox regression analysis to evaluate the prognostic implications of AS events, and we constructed ROC models to show significant AUC values. We also performed IHC staining to identify the differential abundant hub SF DDX21 in ACC and normal tissues from FUSCC. To further improve the clinical value of our study, we enrolled 39 patients with ACC. The findings suggested that DDX21 expression was significantly correlated with aggressive progression and poor prognosis for the patients, which may shed light on the mechanism of ACC development. Third, this study lacks an in-depth mechanism investigation. In the future, we plan to conduct an in-depth study of the mechanism of ACC by complex bioinformatic screening and identification, and the real-world cohort data will be used for validation.

Conclusion
In conclusion, the implementation of strict standards ensured the comprehensive screening and identification of ubiquitous AS patterns relevant to ACC. Functional enrichment analysis identified eight SFs, including BAG2, CXorf56, DDX21, HSPB1, MBNL3, MSI1, RBMXL2 and SEC31B, as essential regulatory elements of splicing pathways in ACC. We found that PASEs not only play important roles in detecting potential mechanisms of AS in tumorigenesis, but also act as novel clinical signatures and targets for treatment strategies. Importantly, systematic profiles of AS events based on powerful large-scale genome-wide cohorts provide the opportunity to comprehensively elucidate the landscape and underlying mechanism of AS, providing valuable insights