Check for updates
Systematic Genome-Wide Profiles Reveal Alternative Splicing Landscape and Implications of Splicing Regulator DExD-Box Helicase 21 in Aggressive Progression of Adrenocortical Carcinoma
Wenhao Xu1,2 . Aihetaimujiang Anwaier1,2 . Wangrui Liu3 . Xi Tian1,2 . Wen-Kai Zhu1,2 . Jian Wang3 . Yuanyuan Qu1,2 . Hailiang Zhang1,2 . Dingwei Ye1,2(D
Received: 1 June 2021 / Revised: 14 September 2021 / Accepted: 18 September 2021 / Published online: 29 October 2021 @ International Human Phenome Institutes (Shanghai) 2021
Abstract
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 predict- ing prognosis for patients with ACC.
Keywords Adrenocortical carcinoma · Genome-wide analysis · Prognosis · Alternative splicing · Bioinformatics · DExD- Box Helicase 21 (DDX21)
Wenhao Xu, Aihetaimujiang Anwaier and Wangrui Liu are co- first authors and contribute equally to this study. Corresponding authors: Yuan-Yuan Qu, Hai-Liang Zhang and Ding-Wei Ye.
☒ Yuanyuan Qu quyy1987@163.com
☒ Hailiang Zhang zhanghl918@alu.fudan.edu.cn
☒ Dingwei Ye
1 Department of Urology, Fudan University Shanghai Cancer Center, No. 270 Dong’an Road, Shanghai 200032, People’s Republic of China
2 Department of Oncology, Shanghai Medical College, Fudan University, Shanghai 200032, People’s Republic of China
3 Department of Transplantation, Xinhua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, Shanghai 200092, People’s Republic of China
Abbreviations
| ACC | Adrenocortical carcinoma |
|---|---|
| AS | Alternative splicing |
| AS | Alternative splicing |
| AA | Alternate acceptor site |
| AD | Alternate donor site |
| AP | Alternate promotor |
| AT | Alternate terminator |
| ES | Exon skip |
| ME | Mutually exclusive exons |
| RI | Retained intron |
| GO | Gene ontology |
| BP | Biological process |
| CC | Cellular components |
| MF | Molecular functions |
| KEGG | Kyoto Encyclopedia of Genes and Genomes |
| PPI | Protein-protein interaction |
| TCGA | The Cancer Genome Atlas |
| SFs | Splicing factors |
| OS | Overall survival |
| DFS | Disease-free survival |
| CI | Confidence interval |
| PSI | Percent-spliced-in |
| PASEs | Prognosis-related AS events |
| GSEA | Gene set enrichment analysis |
| IHC | Immunohistochemical |
Introduction
Adrenocortical carcinoma (ACC) is a rare endocrine malig- nancy 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 rela- tively difficult, so most patients have metastasized tumors at initial diagnosis, resulting in few chances of surgical inter- vention 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 Administra- tion for ACC treatment for patients with metastasis (Schtein- gart 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 stand- ard 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 pro- cess, 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 bioin- formatics 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 accep- tor 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 and Rio 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 break- throughs in many cancers (Chen et al. 2019; 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.
Materials and Methods
Acquisition, Normalization, and Processing of Raw Data Based on Multiple Databases
The RNA sequencing profiles and matched clinicopatho- logical data of 76 ACC patients were obtained from TCGA database (Tomczak et al. 2015). University of North Caro- lina 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://bioinformatics.mdand erson.org/TCGASpliceSeq/) (Ryan et al. 2012). The per- cent-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 (Li et al. 2019). The PSI value can be estimated as the ratio of the density of inclusion reads (i.e., reads per position in regions sup- porting 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” pack- age 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 soft- ware was used for LASSO regression analysis among AA, AD, AP, AT, ES, ME, RI, and all AS events. A multivari- ate Cox regression analysis was used to calculate the risk score for each ACC patient on the basis of the correspond- ing PSI value of highly important prognostic AS (Li et al. 2019).
Survival Analysis and ROC Construction of PASEs
Clinical and pathological parameters, including age, gen- der, 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 fur- ther investigate significant independent clinicopathological factors, including age (ref. Low), gender (ref. Female), T stage (ref. T1), N stage (ref. N0), M stage (ref. M0), stage (ref. stage1), and risk score (ref. Low) of ACC, univariate and multivariate Cox regression analyses were performed using “Survival” package of eight groups. The receiver- operating 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.introni.it/spliceaid.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 Spear- man’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.genemania.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 path- ways related to mRNA splicing was developed with cor- relation 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” pack- age 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 Univer- sity Shanghai Cancer Center (FUSCC; Shanghai, China) from April 2015 to May 2019, were enrolled in this study. Tissue samples were collected during surgery and avail- able from FUSCC tissue bank. IHC staining of DExD-Box Helicase 21 (DDX21) was performed using a mouse mon- oclonal 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 inten- sity 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 microen- vironment. 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 Collec- tion (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% CO2 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 rea- gent (Invitrogen) according to the manufacturer’s instruc- tions. 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, includ- ing 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-AACt 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, includ- ing Upset plot, Cox regression analyses, LASSO regres- sion 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 iden- tified after matching SpliceSeq sequences and clinical
pathological data in TCGA-ACC patients. Second, sur- vival 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 pat- tern 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 (Sup- plementary 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 prog- nosis of ACC patients (p<0.001, Fig. 3A-H) with high- risk 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).
TCGA SpliceSeq
TCGA database
SpliceAid2
1
Genomic Data Commons Data Portal
Average PSI value ≥ 0.05 PSI(min) standard deviation ≥ 0.01
-
-
-
—
Delete duplicated variables External validation
1
[305,465
3.142.246
M
23984 AS events, 8266 genes (No. Samples = 76)
Clinicopathological information and mRNA expression profiles of ACC (No. Samples 76)
Splicing factors (SF) (n=8)
=1
L
!
7
V
1
/
y
4
MatchedACC samples (No. Samples = 76)
SW-13
NCI-H295R
SPRY1
SPRY1
Control
sh-DDX21
MSI
Tel3
Tel3
MSPB:
DOK21
PCBP2
PCBP2
BAGE
MBNL3
NFATc3
NFATc3
.
KDM5A
KDMSA
CNOT4
CNOT4
-
Prognosis-related AS Eevents
DDX21
DDX21
-
0
1
2
3
Fold change
0
3
Fold change
4
4
.
2
-
LASSO and Survival analysis
SF and PPI networks
Functional enrichments
Clustering
Real-world validation
Risk
High risk
9
1.00
og
Overall Survival
Survival probability
0.75
True positive rate
₹
M
0.50
0
1:
0.25
p=1.370-09
AA AUCH0.972
AD AUC-0.972
-
Tine (montha)
0.00
0
1
2
3
1
6
9
10 11 12
0
AT AUG-9912
Progression-free Survival
Time(years)
ES AUG-9.972
ME AUC-9.972
AUCH9.972
₴
-18
0021™
Risk
All
High risk Low
12
?
AUC-9.972
13
&
3
0
2
3
10 1
12
0.0
0.2
0.4
0.6
0.8
1.0
.
Time(years)
False positive rate
-0.020, Hiljhigh)=2.781.
20
Time |montha)
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). Univari- ate Cox analyses showed that age (ref. Low), N stage (ref. N0), M stage (ref. M0), and risk score (ref. Low) of eight groups were independent prognostic parameters in 76 ACC patients (p<0.05). Multivariate Cox analysis
indicated that poorer OS was significantly correlated with risk score (ref. Low) of eight groups, T stage (ref. T1-T2) of AA, AD, AP, AT, ES, ME, RI, and all of AS groups (p < 0.05). N stage (ref. N0) was significantly associated with prognosis in AA, AT, ES, RI, and all of AS groups (p<0.05).
Interaction and PPI Networks of SFs
A total of eight SFs, including BAG2, CXorf56, DDX21, HSPB1, MBNL3, MSI1, RBMXL2 and SEC31B, were iden- tified 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
A
B
10000
1
2
3
4
5
6
2
3.1
3.2
Alternate acceptor site ( AA )
Numbers of cases
7000
4.1
4.2
5
Alternate Donor site ( AD )
1.1
1.2
2
Alternate Promoter ( AP )
4000-
5
6.1
6.2
Alternate Terminator ( AT )
2
3
4
Exon skip ( ES)
1000
“
100
2
3
4
5
Mutually exclusive exions ( ME )
50
0
T
AA Gene
T
T
T
1
RI
Retained intron ( RI )
AD
Gene
AP
Gene
AT
Gene
ES
Gene
ME
Gene
Gene
3
4
C
D
.595
651
31
1500
600
Gene Intersections
Gene Intersections
1000
400
02
DO
500
៛
200
2Tyr
140 3%
3
2go
0
U
ME
V
0
ME
ra
AA
AD
RI
AA
AD
AP
AP
AT
AT
ES
ES
4000
2000
0
800
600
400
200
0
Set Size
E
Set Size
F
AA
G
AD
H
AP
Prognosis-related AS
FCF1-28423-AA
TRAFD1-24580-AD
BLOCIS1-22229-AP
No significant
ZNF692-10567-AA
ZSCAN18-52407-AD
UNG-24277-AP-
MED11-38579-AA
BEX2-89726-AD-
UNG-24278-AP
CERSS-21680-AA
UQCRQ-73323-AD-
pvalue
WASHAP-32782-AA-
pvalue
RPS6-214603-AD
0.00025
CMC2-37703-AP
DUT-30485-AP
pvalue
POLR2H-67943-AA
NDUFAF2-72172-AD
0.00020
CMC2-37702-AP
30-05
-log10(pvalue)
RHOC-4238-AA
60-04
PPP1R8-1347-AD
0.00015
0.00010
DUT-30484-AP
20-05
CIRBP-46427-AA
40-04
NOP2-19895-AD
PGRMC2-70575-AP
GUK1-10188-AA
20-04
GTF3C1-35655-AD
0.00005
MRPL51-19869-AP
10-05
៛
CDK10-38114-AA
ANAPC11-44212-AD
MRPL51-19868-AP
PEX10-266-AA
GOLGA3-25307-AD-
RUSC1-8079-AA
HTATIP2-14714-AD
-log 10(pvalue)
P14K2A-12728-AP
-log10(pvalue)
4
KIAA0391-27213-AP
-log 10[pvalue)
ABHD148-65145-AA
GCDH-47885-AD
LNPEP-72874-AP
2
RABGGTA-26951-AA
3.5
UBL7-31724-AD
5
PIK2A-12729-AP
5
HDDC3-32524-AA
4.0
FOXR-43338-AD
6
SUMO2-43379-AP
6
MECR-1425-AA
4.5
MIF4GD-43427-AD
7
TFDP1-26386-AP
7
RNF181-54295-AA
DPY30-53146-AD
8
ZKSCAN5-80654-AP
8
ZNF577-51379-AA
SRSF5-28151-AD
CCAR2-83035-AP
0
HMGA2-22883-AA
ABCC5-67824-AD
TCEB2-33298-AP
-4
-2
0
2
4
SEC16A-88181-AA
.
UBL5-47437-AD
.
GHDC-41016-AP
-2.5
0.0
2.5
-5.0
-25
0.0
25
5.0
-6
-3
0
3
6
z-score
z-score
z-score
I
AT
J
ES
K
ME
L
z-score RI
KIHL7-78950-AT
CIRBP-46443-ES
EIF6-59080-RI
KLHL7-78952-AT-
CIRBP-46445-ES-
THNSL2-54459-ME
TST-62070-RI-
TECPR2-29416-AT-
HM13-58892-ES
CIRBP-45441-RI-
TECPR2-29415-AT
pvalue
USMG5-13004-ES-
pvalue
FASTK-82334-RI
DNAJC12-11899-AT-
1.20-06
FLAD1-7868-ES-
1.25e-05
H2AFY-95931-ME
pvalue
DYNLL1-24763-RI
pvalue
DNAJC12-11898-AT-
9.00-07
MPND-46796-ES-
1.000-05
0.04
PILRB-80936-RI
30-04
SNX5-58744-AT-
0.03
6.00-07
PCBP2-22056-ES-
ZSWIM7-39393-RI
UCHL5-9240-AT-
MYL6-22384-ES-
7.500-06
RAD51-30020-ME
CMC1-63790-ES-
0.02
FASTK-82333-RI-
2e-04
METTL15-14782-AT-
5.000-06
3.00-07
AIG1-77972-AT
PSEN2-10033-ES-
2.50e-06
0.01
ARLEIP4-25024-RI
10-04
CTNNB1-64249-RI-
USP4-64852-AT
GOCX-54288-ES-
RABA-17707-ME
EPOR-47689-RI-
DFNB31-87326-AT
-log10(pvalue)
INO80E-36016-ES
-log10(pvalue)
-log10(pvalue)
C160091-33097-RI
-log10(pvalue)
ABCA3-33266-AT
6.0
PCBP2-22055-ES
· 5
GCOH-47884-ME
· 2
ATP5D-45401-RI-
ABCA3-33265-AT
6.4
GTF3C5-88006-ES
6
CLN3-35746-RI
4
UCHLS-9239-AT
5
ARHGEF28-72493-AT
6.8
EIF4H-80063-ES-
7
3
CCDC107-86267-RI-
COA4-17737-ES
6
ARHGEF28-72492-AT
72
ISCU-24244-ES
8
ACADS-24779-ME
4
Cfort50-2117-RI-
NAA38-81579-RI
7
RNF180-72196-AT
SEC24C-12227-ES
GSTK1-82082-RI
RNF180-72195-AT
PPHLN1-21221-ES
EMID1-61575-ME
TRAPPC1-39078-RI-
KLHL3-73481-AT
.
PTRH2-42793-ES
.
PSMD11-40209-RI
.
-3
0
3
-6
-3
0
3
6
-4
-2
0
2
-5.0
-25
0.0
25
5.0
z-score
z-score
z-score
z-score
A
Risk
High risk
+ Low risk
B
Risk
High risk
Low risk
C
Risk
High risk
Low risk
1.00-
1.00
1.00-
Survival probability
0.75
Survival probability
0.75
Survival probability
0.75
0.50
0.50
0.50
0.25
p=1.37e-09
0.25
p=1.983e-10
0.25
p=2.031e-09
0.00
0.00
0.00
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
Time(years)
Time(years)
Time(years)
Risk
High risk
38
34 20
16
8
5
3
2
1
1
1
1
1
Risk
High risk
1
38
34
21
13
3
2
Low risk
0
0
0
0
0
0
0
Risk
High risk
38
34
19
17
8
6
3
1
2
2
1
3
2
1
1
38
38
35
26
21
19
13
9
6
Low risk
38
38
34
29
26
22
16
11
8
7
4
2
2
Low risk
38
38
36
25
21
18
13
9
6
5
3
1
1
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
Time(years)
Time(years)
Time(years)
D
Risk
High risk +
Low risk
E
F
+
Risk
+
High risk + Low risk
Risk
L
High risk
+
Low risk
1.00
1.00-
1.00-
Survival probability
0.75
Survival probability
0.75
Survival probability
0.75
0.50
0.50
0.50
0.25
p=6.229e-09
0.25
p=5.009e-10
0.25
p=1.142e-09
0.00
0.00
0.00
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
Time(years)
Time(years)
Time(years)
Risk
High risk
38
38
34 20
38
12 7
5
4
2
Low risk
35
30
22
19
12
CO
0
CO
0
7
0
0
NO
Risk
4
High risk
38
34
19
17
8
5
2
1
Risk
34
20
Low risk
38
38
1
1
0
0
0
High risk
38
6
4
2
36
25
21
19
14
10
7
6
4
2
2
Low risk
38
38
35
15
0
0
8
0
7
0
27
23 20
14
11
4
0 2
0 2
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
Time(years)
Time(years)
Time(years)
G
Risk
H
ROC curve
+
High risk +
Low risk
Risk
+
High risk
++
Low risk
1.0
1.00
1.00
Survival probability
0.75
Survival probability
0.75
0.8
True positive rate
0.50
0.50
0.6
0.25
AA
AUC=0.972
p=9.991e-06
0.25
p=1.316e-08
0.4
AD
AUC=0.992
AP
AUC=0.859
0.00
0.00
AT
AUC=0.797
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
0.2
ES
AUC=0.921
Time(years)
Time(years)
ME
AUC=0.817
RI
0.0
AUC=0.944
Risk
High risk Low risk
38
35
21
18
11
7
4
3
3
0
Risk
High risk-
38
34
19
16
9
6
4
0
All
38 37
0
NO
AUC=0.907
34
24
3
y
1
18
17
12
8
4
3
NO
Low risk
38
38
2
36
26
20
18
12
NO
NO
2
9
8
4
0
1
2
3
4
5
6
7
8
9
10
11
12
0
1
2
3
4
5
6
7
8
9
10
11
12
0.0
0.2
0.4
0.6
0.8
1.0
Time(years)
Time(years)
False positive rate
and transfected in human ACC cells. The results suggested that expression levels of DDX21 and the targeted genes (CNOT4, KDM51, NFATc3, PCBP2, Tcf3, and SPRY1) were significantly decreased in sh-DDX21 group compared with normal control in SW-13 and NCI-H295R cells, as shown in Fig. 5C.
cantly correlated with prognosis of ACC patients (p<0.001, Fig. 4A- H) with high-risk group in red and low-risk group in blue. I ROC analysis of each group
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, cellu- lar response to VEGF stimulus, mRNA processing and RNA splicing; alterations of SFs in CC were significantly enriched in chaperone complex, endoplasmic reticulum exit site, poly- some 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
u
2ª
M
Hazard 1200
4
ka
Hazard ratio
| A | Univariate Cox Regression | Multivariate Cox Regression | B | Univariate Cox Regression | Multivariate Cox Regression | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| All | pvalue | Hazard ratio | All | pvalue | Hazard ratio | AA | pvalue | Hazard ratio | AA | pvalue | Hazard ratio | |||||
| age | 0.334 | 1.012(0.988-1.037) | age | 0.450 | 1.012(0.982-1.043) | age | 0.334 | 1.012(0.988-1.037) | age | 0.794 | 1.004(0.974-1.035) | |||||
| gender | 0.991 | 1.004(0.470-2.145) | gender | 0.994 | 1.003(0.441-2.281) | gender | 0.991 | 1.004(0.470-2.145) | gender | 0.777 | 0.884(0.376-2.080) | |||||
| T | <0.001 | 10.980(4.283-28.150) | T | 0.014 | 31.521(2.039-487.407) | T | <0.001 | 10.980(4.283-28.150) | T | 0.005 | 58.123(3.365-1003.999) | |||||
| N | 0.107 | 2.221(0.842-5.857) | N | 0.026 | 5.882(1.234-28.027) | N | 0.107 | 2.221(0.842-5.857) | N | 0.007 | 11.778(1.962-70.709) | |||||
| M | <0.001 | 7.351(3.305-16.351) | M | 0.988 | 1.009(0.326-3.120) | M | <0.001 | 7.351(3.305-16.351) | M | 0,471 | 1.514(0.490-4.681) | |||||
| stage | <0.001 | 7.166(3.027-16.960) | stage | 0.338 | 0.267(0.018-3.968) | stage | <0.001 | 7.166(3.027-16.960) | stage | 0.139 | 0.115(0.007-2.011) | |||||
| riskScore | <0.001 | 1.019(1.012-1.026) | riskScore | <0.001 | 1.018(1.010-1.027) | riskScore | <0.001 | 1.011(1.007-1.015) | riskScore | <0.001 | 1.010(1.005-1.016) | |||||
| 2.0 - wa | 12 = | - 45 | ||||||||||||||
| C AD | pvalue | Hazard ratio | AD | pvalue | Hazard ratio | D AP | pvalue | Hazard ratio | Hatand ratio | AP | pvalue | Hazard ratio | ||||
| age | 0.334 | 1.012(0.988-1.037) | age | 0.334 | 1.016(0.983-1.050) | age | 0.334 | 1.012(0.988-1.037) | age | 0.079 | 1.028(0.997-1.061) | |||||
| gender | 0.991 | 1.004(0.470-2.145) | gender | 0.355 | 0.663(0.278-1.584) | gender | 0.991 | 1.004(0.470-2.145) | gender | 0.380 | 0.682(0.290-1.604) | |||||
| T | <0.001 | 10.980(4.283-28.150) | T | 0.017 | 26.456(1.784-392.328) | T | <0.001 | 10.980(4.283-28.150) | T | 0.165 | 6.699(0.458-97.905) | |||||
| N | 0.107 | 2.221(0.842-5.857) | N | 0.051 | 4.855(0.995-23.702) | N | 0.107 | 2.221(0.842-5.857) | N | 0.705 | 1.341(0.294-6.112) | |||||
| M | <0.001 | 7.351(3.305-16.351) | M | 0.283 | 1.897(0.589-6.107) | M | <0.001 | 7.351(3.305-16.351) | M | 0.282 | 1.853(0.603-5.693) | |||||
| stage | <0.001 | 7.166(3.027-16.960) | stage | 0.373 | 0.287(0.018-4.468) | stage | <0.001 | 7.166(3.027-16.960) | stage | 0.829 | 1.332(0.098-18.049) | |||||
| riskScore | <0.001 | 1.026(1.012-1.040) | riskScore | <0.001 | 1.026(1.013-1.039) | riskScore | <0.001 | 1.035(1.020-1.049) | riskScore | <0.001 | 1.031(1.015-1.048) | |||||
| 3. . - | Hazard | 24 | 4g | |||||||||||||
| E AT | pvalue | Hazard ratio | Harand rado | AT | pvalue | Hazard ratio | ratio | F ES | pvalue | Hazard ratio | Macard tado | ES | pvalue | Hazard ratio | Harand radio | |
| age | 0.334 | 1.012(0.988-1.037) | age | 0.803 | 1.004(0.974-1.034) | age | 0.334 | 1.012(0.988-1.037) | age | 0.528 | 1.010(0.980-1.040) | |||||
| gender | 0.991 | 1.004(0.470-2.145) | gender | 0.792 | 0.897(0.400-2.013) | gender | 0.991 | 1.004(0.470-2.145) | gender | 0.932 | 0.964(0.414-2.242) | |||||
| T | <0.001 | 10.980(4.283-28.150) | T | 0.017 | 25.446(1.768-366.224) | T | <0.001 | 10.980(4.283-28.150) | T | 0.019 | 25.491(1.696-383.040) | |||||
| N | 0.107 | 2.221(0.842-5.857) | N | 0.032 | 5.194(1.148-23.494) | N | 0.107 | 2.221(0.842-5.857) | N | 0.022 | 6.438(1.301-31.845) | |||||
| M | <0.001 | 7.351(3.305-16.351) | M | 0.625 | 1.306(0.449-3.799) | M | <0.001 | 7.351(3.305-16.351) | M | 0.170 | 2.148(0.721-6.399) | |||||
| stage | <0.001 | 7.166(3.027-16.960) | stage | 0.299 | 0.246(0.017-3.475) | stage | <0.001 | 7.166(3.027-16.960) | stage | 0.305 | 0.240(0.016-3.660) | |||||
| riskScore | <0.001 | 1.078(1.049-1.107) | riskScore | <0.001 | 1.062(1.028-1.097) | riskScore | <0.001 | 1.038(1.024-1.052) | riskScore | <0.001 | 1.037(1.018-1.057) | |||||
| G ME | pvalue | 0.50 20 48 Hazard ratio . 16.0 | H | as La Ł 4 | 8 4ª | Hazard ratio | ||||||||||
| Hazard ratio | ME | pvalue | Hazard ratio | RI | pvalue | Hazard ratio | RI | pvalue | Hazard ratio | |||||||
| age | 0.334 | 1.012(0.988-1.037) | age | 0.336 | 1.014(0.986-1.043) | age | 0.334 | 1.012(0.988-1.037) | age | 0.254 | 1.017(0.988-1.046) | |||||
| gender | 0.991 | 1.004(0.470-2.145) | gender | 0.677 | 0.840(0.371-1.903) | gender | 0.991 | 1.004(0.470-2.145) | gender | 0.626 | 0.807(0.340-1.916) | |||||
| T | <0.001 | 10.980(4.283-28.150) | T | 0.023 | 18.042(1.485-219.197) | T | <0.001 | 10.980(4.283-28.150) | T | 0.006 | 56.137(3.250-969.593) | |||||
| N | 0.107 | 2.221(0.842-5.857) | N | 0.156 | 2.743(0.681-11.057) | N | 0.107 | 2.221(0.842-5.857) | N | 0.008 | 10.590(1.836-61.088) | |||||
| M | <0.001 | 7.351(3.305-16.351) | M | 0.224 | 1.922(0.671-5.509) | M | <0.001 | 7.351(3.305-16.351) | M | 0.903 | 1.073(0.348-3.305) | |||||
| stage | <0.001 | 7.166(3.027-16.960) | stage | 0.528 | 0.447(0.037-5.446) | stage | <0.001 | 7.166(3.027-16.960) | stage | 0.200 | 0.159(0.009-2.647) | |||||
| riskScore | <0.001 | 1.391(1.227-1.578) | riskScore | <0.001 | 1.515(1.268-1.812) | riskScore | <0.001 | 1.073(1.049-1.099) | 1 | riskScore | <0.001 | 1.074(1.043-1.106) | ||||
Fig. 4 Univariate and multivariate Cox regression analyses of PASEs in eight groups were performed in forest plots of A AA, B AD, C AP, D AT, E ES, F ME, G RI, and H all types of AS events, respectively
(Fig. 5D). Pathway analysis suggested that alterations in KEGG are mainly enriched in protein processing in endo- plasmic reticulum, VEGF signaling pathway, mRNA sur- veillance pathway, amoebiasis and spliceosome; alterations in Reactome pathway enrichment analysis existed in AUF1 that bind and destabilize mRNA, regulation of HSF1-medi- ated 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 enrich- ment in blue, and other top 20 related enrichment biologi- cal 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 exam- ple, 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 sig- nificantly associated with shorter OS (p<0.0001, HR 4.6, Fig. 6B) and DFS (p<0.0001, HR 4.3, Fig. 6C). The expres- sion 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 accord- ing to differential DDX21 expression groups (Fig. 6H).
A
E
€
B
a
STUB1
TRAF6
TUBA1A TUBB3
PPI
Activation
RPL 10
SIK2
SIM2
SRRM2
STRN
TP53
TYK2
UBC
PSMD6
VCP
PPP2R2B
WDR83
Phosphorylation
RELA
XRCC5
CXorf56
Inhibition
OTUB1
PAN2
APP
JUN
SNW1
NPM1
DDX21
ESR1
MSI1
NFKB2
NOS2
GABARAP
MYC
GABARAPL1
HSPB1
DDX21
MED19
GABARAPL2
MAPK6
MAP3K14
HSCB
LRRK2
HUWE1
BAG2
MBNL3
ITGA4
LIMK2
KPNA3
IRAK1
MDM2
ILK
IKBKE
NFKB1
CXorf56
SEC31B
HSPAS
SAP18
RBMXL2
HSPA4
HSP90AA1
HNRNPU
SUMO1
SRPK1
HNRNPA1
YBX1
VCAM1
FN1
EGFR
DHX9
CUL3
CFTR
BUB1B
MAPKAPKS
MAPKAPK3
MAPKAPK2
MAP1LC3B
MAP1LC3A
HSPB1
YWHAE
BAG3
ARAF
AICDA
AARSD1
BAG2
SEC31B
RBMXL2
MSI1
MBNL3
C
D
ubiquitin-like protein binding-
single-stranded RNA binding
SW-13
NCI-H295R
Control
double-stranded RNA binding
rRNA binding
sh-DDX21
protein kinase C binding
SPRY1
H ***
SPRY1
tau protein binding
snRNA binding
Enrichment_Ratio
regulatory RNA binding
Tcf3.
30
Tcf3
snoRNA binding
protein binding involved in protein folding
60
PCBP2
PCBP2
4 **
coated membrane
90
peptidase complex
NFATc3.
NFATc3
polysome
P_Value
endoplasmic reticulum exit site
0.04
KDM5A
H **
KDM5A
H **
chaperone complex
cellular response to VEGF stimulus
0.03
mRNA processing
CNOT4.
0.02
CNOT4.
positive regulation of cytokine production
RNA splicing
0.01
DDX21
H ***
DDX21
H ***
cargo loading into vesicle
response to virus
0
1
2
3
0
1
2
3
regulation of mRNA metabolic process
Fold change
Fold change
negative regulation of transferase activity
I-kappaB kinase/NF-kappaB signaling
protein folding
Biologic Process
Cellular Component
Molecular Function
E
F
ESPONSE
TRANSCRIPTION FROM
ENDOGENOUS
MERASE
Signaling by VEGF
EPOR
Positive epigenetic regulation of rRNA expression
VEGFA-VEGFR2 Pathway
P_Value
RRMA
B-WICH complex positively regulates rRNA expression
0.06
MAPK6/MAPK4 signaling-
ETABOL
Regulation of mRNA stability by proteins that bind AU-rich
0.04
elements
Cellular response to heat stress
0.02
Regulation of HSF1-mediated heat shock response
000639
TION
Metabolism of RNA
Enrichment_Ratio
AUF1 (hnRNP DO) binds and destabilizes mRNA
NUCLEOBASENU
HE OSIDENUCLEOTID
10
POLTHIS ETABOL
Spliceosome
20
PROCE
Amoebiasis
30
ROTEIN
mRNA surveillance pathway
40
COMPTEK ASSEMBLY
CESSING
VEGF signaling pathway
LOTEIN
Protein processing in endoplasmic reticulum
COMPTE ASSENDLY
KEGG
Reactome
STIMULUS
TABOL
AND ASSEMBLY
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 con- structed using all significant AS event and DDX21 mRNA
A
12.00
SEC31B
RBMXL2
MSI1
HSPB1
MBNL3
DDX21
CXorf56
BAG2
PSI_ALL
PSI_ES
2
₹
7
PSI_ME
B
D
AA
Overall Survival
11.50
Spearman rho=0.283
PSI
PSI
82
PS
PS
9
Low DDX21 TPM
p=0.0127
High DDX21.TPM
11.00
1
Logrank p=0.00017
SEC31B
log2 RSEM value
10.50
1
0.22
0.28
0.24
-0.01
-0.07
-0.08
-0.18
0.14
0.05
0.02
0.04
0.03
0.05
0.06
-0.14
3
HR(high)=4.6
p(HR)=0.00056
Percent survival
n(high)=38
10.00
RBMXL2
0.22
1
0.3
0.25
-0.13
-0.27
-0.16
-0.16
0.13
0.05
-0.04
-0.01
-0.1
0.01
0
-0.12
0.8
0.6
n|low)=38
9.50
MSI1
0.28
0.3
1
0.61
0.16
-0.41
-0.02
-0.18
-0.13
0.01
0.08
0.02
0.1
0.11
-0.31
-0.2
9.00
HSPB1
0
0.24
0.25
0.61
1
-0.1
-0.6
-0.14
8.50
0.04
-0.05
-0.01
0.04
-0.01
0.06
0.05
-0.14
=0.24
0.6
8.00
MBNL3
-0.01
-0,13
0.16
-0.1
1
0,38
0.26
-0.26
0
0,07
0.09
0.09
0.09
0.04 -0.09
-0.06
0
0,4
7.50
DDX21
-0.07
-0.27
-0.41
-0.6
0.38
1
0.56
-0.15
0.25
0.16
0.14
0.18
0.17
0.13
0.24
0.17
0
7.00
STAGE I
STAGE II
STAGE III
STAGE IV
CXorf56
-0.08
-0.16
-0.02
-0.14
0.26
0.56
1
-0.16
-0.03
-0.08
0
-0.03
0
-0.03
-0.1
-0.17
0.2
0
50
100
150
Months
BAG2
-0.18
-0.16
-0.18
0.04
-0.26
-0.15
-0.16 1
0.23
0.15
0.13
0.15
0.11
0.1
0.26
0.17
0
E
PSI_ALL
0.14
0.13
-0.13
-0.05
0
0.25
-0.03
0.23
1
0.85
0.75
0.82
0.73
0.76
0.68
0.28
C
PSI_ES
0.05
0.06
0.01
-0.01
0.07
0.16
-0.08
0.15
0.85
1
0.92
0.96
0.89
0.92
0.41
0.24
-0.2
Disease Free Survival
9
Low TPM
PSI_AD
0.02
-0.04
0.08
0.04
0.09
0.14
0
0.13
0.75
0.92
1
0.98
0.98
0.98
0.24
0.3
High DDX21 TPM
Logrank p=3.6e-05
Normal
HR(high)=4.3
PSI_RI
0.04
-0.01 0.02
-0.01
0.09
0.18
-0.03
0.15
0.82
0.96
0.98
1
0.95
0.97
0.34
0.32
-0.4
3
p(HR)=0.00013
Percent survival
talhigh) 36
PSI_AA
0.03
-0.1
0.1
0.06
0.09
0.17
0
0.11
0.73
0.89
0.98
0.96
1
0.97
0.29
0.34
0.6
n(low)=38
-0.6
PSI_AP
0.05
0.01
0.11
0.06
0.04
0.13
-0.03
0.1
0.76
0.92
0.98
0.97
0.97
1
0.29
0.36
0.4
PSI_AT
0.06
0
-0.31
-0.14
-0.09
0.24
-0.1
0.26
0.68
0,41
0.24
0.34
0.29
0.29
1
0.3
-0.8
PSI_ME -0.14
-0.12
-0.2
-0.24
-0.06 0.17
-0.17
0.17
0.28
0.24
0.3
0,32
0.34
0,35
0.3
1
8
Tumor
-1
F
G
8
Enrichment plot: HALLMARK_MITOTIC_SPINDLE
0
50
100
150
Enrichment score (ES)
NES=0.183
4
Months
La
NOM p=0.018
12
L1
Data & Annotation: Junction Junction usage
pathological_stage:
Ranked list metric (Signal)tvoise)
Transcription pattern Gene expression
STAGE I
STAGE II
STAGE III
MA
6
STAGE IV
14
UNDEFINED
·
5.000
Rank in Ordered Dataset
Enrichment proflig . Hits
Ranking metric scores
pathological_stage
H
DDX21
-
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 supple- mented 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 aggres- sive 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 aggres- sive 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 characteri- zation have limited clinical prediction abilities for patients with ACC, who have high mortality rates and complex treat- ment 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 regula- tion in diseases such as cancers (Rahman et al. 2020). How- ever, 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
A
Overall Survival
C
Overall Survival
DDX21lov
DDX21low
Percent survival (%)
100
n=16
DDX21high
Percent survival (%)
100
n=8
DDX21high
50-
50-
n=17
n=23
p=0.005, HR(high)=4.171
p=0.011, HR(high)=8.556
All patients
Necrosispos patients
0
0
0
20
40
60
80
100
0
20
40
60
80
100
Time (months)
Time (months)
B
Progression-free Survival
D
Progression-free Survival
DDX21lov
DDX21low
Percent survival (%)
100
n=8
n=16
DDX21 high
Percent survival (%)
100
DDX21high
50.
50-
n=23
n=17
p=0.020, HR(high)=2.781
p=0.020, HR(high)=4.748
0
0
0
20
40
60
80
0
20
40
60
80
Time (months)
Time (months)
diverse cancers, while the lack of current sequencing tech- nologies and analysis processes hinders systematic analy- sis. 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). Accu- mulating 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 cor- responding genes were identified after matching SpliceSeq and clinical pathological data in TCGA-ACC patients. Nota- bly, 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-spe- cific AS events and SF expression profiles in ACC data- set have not been reported yet (Bates 2020). Hence, we integrated PSI-specific AS events and significant SF net- works to investigate the splicing pathway mechanism. We found that alterations in SFs were markedly enriched in the mRNA splicing process. Interestingly, the vascular endothe- lial 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 (Zhao and Zhang 2018).
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 (heli- case) 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. Comprehen- sively 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. There- fore, 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 target- ing DDX21, and they brings more treatment options for patients with high-grade adrenocortical carcinoma.
This study has several limitations. First, in May 2020, (Xu et al. 2020) reported survival-related alternative SFs based only on TCGA SpliceSeq data without external validation. Although there were some similarities in the initial data pro- cessing 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 clin- ical translational value, we performed a survival risk assess- ment 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 analy- sis identified eight SFs, including BAG2, CXorf56, DDX21, HSPB1, MBNL3, MSI1, RBMXL2 and SEC31B, as essen- tial 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 strate- gies. 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
into the potential of DDX21 for predicting the prognosis of ACC patients.
Supplementary Information The online version contains supplemen- tary material available at https://doi.org/10.1007/s43657-021-00026-x.
Acknowledgements We are grateful to all patients for their dedicated participation in the current study. We expressed our sincere gratitude to Ms. Zoo for editing figures and drawing the schematic diagram for this study.
Authors’ contributions Conceptualization: WX, AA, WL, and HZ. Data curation and formal analysis: WX, AA, WL, JW, XT, and WZ. Funding acquisition: YQ, HZ, and DY. Investigation and methodology: WX, AA, WZ, and WL. Resources and software: WL, WZ, YQ, HZ, and DY. Supervision: YQ, JW, HZ and DY. Validation and visualiza- tion: WX, WL, XT, and AA. Original draft: WX, AA, and WL. Edit- ing: YQ, HZ, and DY.
Funding This work is supported by National Key Research and Devel- opment Program of China (No. 2019YFC1316000), Fuqing Scholar Student Scientific Research Program of Shanghai Medical College, Fudan University (No. FQXZ202112B), Natural Science Foundation of Shanghai (No. 20ZR1413100), and Shanghai Municipal Health Bureau (No. 2020CXJQ03).
Availability of data and materials The datasets analyzed in this study were obtained from the corresponding author upon reasonable request or open-access online databases.
Code availability The R code is publicly available, and it can be found by name on the GitHub site listed in the Data and Material Availability (https://github.com/) section.
Declarations
Conflicts of interest The authors have no conflicts of interest to declare that are relevant to the content of this article.
Ethics approval and consent to participate Study procedures were approved by Fudan University Shanghai Cancer Center (Shang- hai, China) (ID: 050432-4-1805C). Written informed consents were acquired from open-access online databases.
Consent for publication Not applicable.
References
Bates SE (2020) Epigenetic therapies for cancer. N Engl J Med 383(7):650-663. https://doi.org/10.1056/NEJMra1805035
Cao J, Wu N, Han Y, Hou Q, Zhao Y, Pan Y, Xie X, Chen F (2018) DDX21 promotes gastric cancer proliferation by regulating cell cycle. Biochem Biophys Res Commun 505(4):1189-1194. https:// doi.org/10.1016/j.bbrc.2018.10.060
Chen T, Zheng W, Chen J, Lin S, Zou Z, Li X, Tan Z (2019) Systematic analysis of survival-associated alternative splicing signatures in clear cell renal cell carcinoma. J Cell Biochem. https://doi.org/ 10.1002/jcb.29590
Climente-González H, Porta-Pardo E, Godzik A, Eyras E (2017) The functional impact of alternative splicing in cancer. Cell Rep 20(9):2215-2226. https://doi.org/10.1016/j.celrep.2017.08.012
Consortium SEQC/MAQC-III (2014) A comprehensive assessment of RNA-seq accuracy reproducibility and information content by the Sequencing Quality Control Consortium. Nat Biotechnol 32(9):903-914. https://doi.org/10.1038/nbt.2957
Feng H, Qin Z, Zhang X (2013) Opportunities and methods for study- ing alternative splicing in cancer with RNA-Seq. Cancer Lett 340(2):179-191. https://doi.org/10.1016/j.canlet.2012.11.010
Giulietti M, Piva F, D’Antonio M, D’Onorio De Meo P, Paoletti D, Castrignanò T, D’Erchia AM, Picardi E, Zambelli F, Principato G, Pavesi G, Pesole G (2013) SpliceAid-F: a database of human splicing factors and their RNA-binding sites. Nucleic Acids Res 41(Database issue):D125-131. https://doi.org/10.1093/nar/gks997 Jasim S, Habra MA (2019) Management of adrenocortical car- cinoma. Curr Oncol Rep 21(3):20. https://doi.org/10.1007/ s11912-019-0773-7
Jonasch E, Walker CL, Rathmell WK (2021) Clear cell renal cell car- cinoma ontogeny and mechanisms of lethality. Nat Rev Nephrol 17(4):245-261. https://doi.org/10.1038/s41581-020-00359-2
Kim DS, Camacho CV, Nagari A, Malladi VS, Challa S, Kraus WL (2019) Activation of PARP-1 by snoRNAs controls ribosome biogenesis and cell growth via the RNA helicase DDX21. Mol Cell 75(6):1270-1285 e1214. https://doi.org/10.1016/j.molcel. 2019.06.020
Kostiainen I, Hakaste L, Kejo P, Parviainen H, Laine T, Löyttyniemi E, Pennanen M, Arola J, Haglund C, Heiskanen I, Schalin-Jäntti C (2019) Adrenocortical carcinoma: presentation and outcome of a contemporary patient series. Endocrine 65(1):166-174. https:// doi.org/10.1007/s12020-019-01918-9
Lee Y, Rio DC (2015) Mechanisms and regulation of alternative Pre- mRNA splicing. Annu Rev Biochem 84:291-323. https://doi.org/ 10.1146/annurev-biochem-060614-034316
Li ZX, Zheng ZQ, Wei ZH, Zhang LL, Li F, Lin L, Liu RQ, Huang XD, Lv JW, Chen FP, He XJ, Guan JL, Kou J, Ma J, Zhou GQ, Sun Y (2019) Comprehensive characterization of the alterna- tive splicing landscape in head and neck squamous cell carci- noma reveals novel events associated with tumorigenesis and the immune microenvironment. Theranostics 9(25):7648-7665. https://doi.org/10.7150/thno.36585
Libe R (2015) Adrenocortical carcinoma (ACC): diagnosis progno- sis and treatment. Front Cell Dev Biol 3:45. https://doi.org/10. 3389/fcell.2015.00045
Liu C, Guo T, Xu G, Sakai A, Ren S, Fukusumi T, Ando M, Sadat S, Saito Y, Khan Z, Fisch KM, Califano J (2018) Characterization of alternative splicing events in HPV-negative head and neck squamous cell carcinoma identifies an oncogenic DOCK5 vari- ant. Clin Cancer Res 24(20):5123-5132
Maniatis T, Tasic B (2002) Alternative pre-mRNA splicing and proteome expansion in metazoans. Nature 418(6894):236-243. https://doi.org/10.1158/1078-0432.CCR-18-0752
Meng T, Huang R, Zeng Z, Huang Z, Yin H, Jiao C, Yan P, Hu P, Zhu X, Li Z, Song D, Zhang J, Cheng L (2019) Identification of prognostic and metastatic alternative splicing signatures in kid- ney renal clear cell carcinoma. Front Bioeng Biotechnol 7:270. https://doi.org/10.3389/fbioe.2019.00270
Nilsen TW, Graveley BR (2010) Expansion of the eukaryotic pro- teome by alternative splicing. Nature 463(7280):457-463. https://doi.org/10.1038/nature08909
Rahman MA, Krainer AR, Abdel-Wahab O (2020) SnapShot: splic- ing alterations in cancer. Cell 180(1):208-208 e201. https://doi. org/10.1016/j.cell.2019.12.011
Ryan MC, Cleland J, Kim R, Wong WC, Weinstein JN (2012) SpliceSeq: a resource for analysis and visualization of RNA- Seq data on alternative splicing and its functional impacts.
Bioinformatics 28(18):2385-2387. https://doi.org/10.1093/bioin formatics/bts452
Ryan M, Wong WC, Brown R, Akbani R, Su X, Broom B, Melott J, Weinstein J (2016) TCGASpliceSeq a compendium of alterna- tive mRNA splicing in cancer. Nucleic Acids Res 44(D1):D1018- 1022. https://doi.org/10.1093/nar/gkv1288
Schteingart DE, Doherty GM, Gauger PG, Giordano TJ, Hammer GD, Korobkin M, Worden FP (2005) Management of patients with adrenal cancer: recommendations of an international consensus conference. Endocr Relat Cancer 12(3):667-680. https://doi.org/ 10.1677/erc.1.01029
Sun W, Duan T, Ye P, Chen K, Zhang G, Lai M, Zhang H (2018) TSVdb: a web-tool for TCGA splicing variants analysis. BMC Genom 19(1):405. https://doi.org/10.1186/s12864-018-4775-x
Tierney JF, Chivukula SV, Poirier J, Pappas SG, Schadde E, Hertl M, Kebebew E, Keutgen X (2019) National treatment practice for adrenocortical carcinoma: have they changed and have we made any progress? J Clin Endocrinol Metab 104(12):5948-5956. https://doi.org/10.1210/jc.2019-00915
Tomczak K, Czerwińska P, Wiznerowicz M (2015) The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (pozn) 19(1A):A68-77. https://doi.org/10.5114/ wo.2014.47136
Treiber T, Treiber N, Plessmann U, Harlander S, Daiß JL, Eichner N, Lehmann G, Schall K, Urlaub H, Meister G (2017) A com- pendium of RNA-binding proteins that regulate microRNA bio- genesis. Mol Cell 66(2):270-284 e213. https://doi.org/10.1016/j. molcel.2017.03.014
Urbanski LM, Leclair N, Anczuków O (2018) Alternative-splicing defects in cancer: splicing regulators and their downstream targets guiding the way to novel cancer therapeutics. Wiley Interdiscip Rev RNA 9(4):e1476. https://doi.org/10.1002/wrna.1476
Valdez BC, Yang H, Hong E, Sequitin AM (2002) Genomic structure of newly identified paralogue of RNA helicase II/Gu: detection of pseudogenes and multiple alternatively spliced mRNAs. Gene 284(1-2):53-61. https://doi.org/10.1016/s0378-1119(01)00888-5
Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, King- smore SF, Schroth GP, Burge CB (2008) Alternative isoform regu- lation in human tissue transcriptomes. Nature 456(7221):470-476. https://doi.org/10.1038/nature07509
Wang J, Xu W, Wang B, Lin G, Wei Y, Abudurexiti M, Zhu W, Liu C, Qin X, Dai B, Wan F, Zhang H, Zhu Y, Ye D (2020) GLUT1 is an AR target contributing to tumor growth and glycolysis in
castration-resistant and enzalutamide-resistant prostate cancers. Cancer Lett 485:45-55. https://doi.org/10.1016/j.canlet.2020.05. 007
Xing YH, Yao RW, Zhang Y, Guo CJ, Jiang S, Xu G, Dong R, Yang L, Chen LL (2017) SLERT regulates DDX21 rings associated with pol I transcription. Cell 169(4):664-678 e616. https://doi.org/10. 1016/j.cell.2017.04.011
Xu WH, Shi SN, Xu Y, Wang J, Wang HK, Cao DL, Shi GH, Qu YY, Zhang HL, Ye DW (2019a) Prognostic implications of Aqua- porin 9 expression in clear cell renal cell carcinoma. J Transl Med 17(1):363. https://doi.org/10.1186/s12967-019-2113-y
Xu WH, Wu J, Wang J, Wan FN, Wang HK, Cao DL, Qu YY, Zhang HL, Ye DW (2019b) Screening and identification of potential prognostic biomarkers in adrenocortical carcinoma. Front Genet 10:821. https://doi.org/10.3389/fgene.2019.00821
Xu WH, Xu Y, Wang J, Wan FN, Wang HK, Cao DL, Shi GH, Qu YY, Zhang HL, Ye DW (2019c) Prognostic value and immune infiltration of novel signatures in clear cell renal cell carcinoma microenvironment. Aging (albany NY) 11(17):6999-7020. https:// doi.org/10.18632/aging.102233
Xu N, Ke ZB, Lin XD, Lin F, Chen SH, Wu YP, Chen YH, Wei Y, Zheng QS (2020) Identification of survival-associated alternative- splicing events and signatures in adrenocortical carcinoma based on TCGA SpliceSeq data. Aging (Albany NY) 12(6):4996-5009. https://doi.org/10.18632/aging.102924
Xu W, Tian X, Liu W, Anwaier A, Su J, Zhu W, Wan F, Shi G, Wei G, Qu Y, Zhang H, Ye D (2021) m6A regulator-mediated methylation modification model predicts prognosis tumor microenvironment characterizations and response to immunotherapies of clear cell renal cell carcinoma. Front Oncol. https://doi.org/10.3389/fonc. 2021.709579
Zhang H, Zhang Y, Chen C, Zhu X, Zhang C, Xia Y, Zhao Y, Andrisani OM, Kong L (2018) A double-negative feedback loop between DEAD-box protein DDX21 and Snail regulates epithelial-mes- enchymal transition and metastasis in breast cancer. Cancer Lett 437:67-78. https://doi.org/10.1016/j.canlet.2018.08.021
Zhao N, Zhang J (2018) Role of alternative splicing of VEGF-A in the development of atherosclerosis. Aging (albany NY) 10(10):2695- 2708. https://doi.org/10.18632/aging.101580
Zuberi K, Franz M, Rodriguez H, Montojo J, Lopes CT, Bader GD, Morris Q (2013) GeneMANIA prediction server 2013 update. Nucleic Acids Res 41(Web Server issue): W115-122. https://doi. org/10.1093/nar/gkt533