Research Paper

Identification of survival-associated alternative splicing events and signatures in adrenocortical carcinoma based on TCGA SpliceSeq data

Ning Xu1,*, Zhi-Bin Ke1,*, Xiao-Dan Lin1,4, Fei Lin1,4, Shao-Hao Chen1, Yu-Peng Wu1, Ye-Hui Chen1, Yong Wei1, Qing-Shui Zheng1

1Departments of Urology, The First Affiliated Hospital of Fujian Medical University, Fuzhou 350005, China *Equal contribution

Correspondence to: Qing-Shui Zheng, Yong Wei; email: zhengqingshui@fjmu.edu.cn, weiyong2017@fjmu.edu.cn

Keywords: alternative splicing, adrenocortical carcinoma, TCGA, survival

Received: December 11, 2019 Accepted: March 2, 2020

Published: March 26, 2020

Copyright: Xu et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

ABSTRACT

Objective: To explore the correlations among alternative splicing (AS), splicing factors (SF) and survival outcome in adrenocortical carcinoma (ACC) patients.

Results: A total of 92 ACC patients were included. Univariate analysis identified 3919 AS events significantly associated with overall survival. Lasso method followed by multivariate analysis revealed that the prognostic capacity of these survival-related AS events is satisfactory. Interestingly, we found that the area under the curve (AUC) of AA, AD, AP and RI were more than 0.9, indicating that these four types of AS were of great significance. Independent prognostic analysis showed that only the risk score was the independent risk factor of ACC survival. Finally, we constructed an interesting interaction network between AS and SF. Conclusions: This is the first and most comprehensive study to explore the aberrant AS variants in ACC, which might provide novel insights into molecular mechanism of ACC.

Methods: The transcriptome data, clinical information and Percent Spliced In (PSI) values of the ACC were obtained from TCGA database and TCGA SpliceSeq data portal. Lasso method and uni/multivariate Cox regression analysis were used to identify survival-related AS events and develop multi-AS-based signatures. The relationship between AS events and SFs was also investigated.

INTRODUCTION

As a rare but greatly aggressive malignancy, the prognosis of adrenocortical carcinoma (ACC) is rather poor [1, 2]. Surgical excision currently remains the preferred treatment for ACC in patients with localized lesion [3]; however, the five-years overall survival is only 15-44% even if the tumor is completely removed [4]. Although there have been numerous studies exploring ACC, the specific molecular mechanism has not been fully elucidated.

As the most important post-transcriptional regulatory process, alternative splicing (AS) modifies more than 90% of human genes [5, 6], contributing greatly to

protein diversity and complexity [7]. Numerous physiological and pathological behaviors in the human body are related to AS, including angiogenesis, proliferation, hypoxia, etc [5, 8, 9]. It was reported that survival predictors of gene expression were consistently inferior to AS-based survival predictors [10]. There were emerging data revealing that abnormal AS events were involved in carcinogenesis, immunologic escape, cancer progression and metastasis [6, 10-12]. However, little work has been devoted to investigate the role of AS in ACC.

In this study, a total of 92 patients with ACC from the Cancer Genome Atlas (TCGA) database were enrolled to comprehensively profile genome-wide AS events,

which play a vital role in the development of ACC. Further, uni/multivariate Cox regression analysis and Lasso regression analysis were used to identify survival-related AS events and develop multi-AS-based signatures. We also explored the relationship between AS events and clinical feature using clinical data from TCGA database. Moreover, AS related splicing factors (SF) were identified and functional enrichment analyses were performed.

RESULTS

Clinical characteristics and integrated AS events

A total of 92 ACC patients from TCGA database were included in the present study. The clinicopathological characteristics of these 92 cases were summarized in Table 1. The AS events were comprehensively analyzed. Intersections among these seven types of AS events and quantitative analysis of these interactive sets were presented by the UpSet plot (Figure 1A). These data suggested that a single gene could possess several types of mRNA splicing events and ES was the predominant type in ACC cohort whereas ME was rare.

Survival associated AS events

Univariate Cox regression analysis between AS events and OS was conducted to determine the survival associated-AS events. A total of 3919 AS events were significantly associated with OS (P<0.05). Con- sequently, one gene might contain two or more AS events that were remarkably related to the OS, and ES was the most common survival associated event. The UpSet plot was generated to visualize the intersecting sets of genes and splicing types vividly (Figure 1B). The top 20 (if available) significant survival associated AS events of each type were selected and visualized in Figure 2. The distributions of AS events significantly associated with OS were displayed in Figure 2D.

Construction and evaluation of the prognostic AS signature

Lasso method followed by multivariate Cox regression analyses were performed to evaluate the prognostic capacity of these survival related AS events. The results of Lasso regression analysis including seven types of AS events and all AS events were presented in Figure 3. Using Lasso regression analysis, we selected the most highly correlated AS events. Next, the particular AS events were selected and risk scores were calculated based on multivariate Cox regression model. The detailed information of particular AS events in the eight prognostic models was shown in Table 2. Moreover, ACC patients were stratified into low-risk and high-risk

groups based on the median value of the risk score as cut off. The Kaplan-Meier curves were employed to demonstrate the survival probability variation of patients in the low-risk and high-risk groups. The results showed that survival time differed notably between these two groups in each type and the whole cohort (Figure 4). The ROC curves were generated to evaluate the predictive power of each prognostic signature. Finally, four prognostic signatures with AUC ≥ 0.9 were selected for subsequent analysis. The results revealed that risk score of AD performed the greatest prognostic power with an AUC of 0.983, followed by risk score of AA with an AUC of 0.969 and risk score of RI with an AUC of 0.928 (Figure 5). The detailed information of corresponding splicing pattern of the candidate AS events, living status as well as survival time ranked by the distribution of risk score was displayed (Figure 6).

Stratification survival analysis with the prognostic signatures

The following clinical parameters including gender, T stage, N stage and clinical stage were integrated into the uni/multivariable Cox regression analysis. Univariate Cox regression analysis indicated that several clinical features could predict poorer survival of ACC patients including high T stage, high N stage and high clinical stage and high risk-score. However, multivariate Cox regression showed that the risk score was the only independent risk factor of ACC survival (Figure 7). These results revealed that the risk scores with an AUC ≥ 0.9 were independent prognostic indicators after adjusted by the clinicopathological characteristics.

Functional enrichment analysis was performed for the sake of revealing the potential mechanisms of corresponding genes (known as SFs) with survival- associated AS events. The biological process terms of these genes were mainly related to “mRNA processing”, “RNA splicing”, “RNA splicing via transesterification reactions”. “RNA splicing via transesterification reactions with bulged adenosine as nucleophile”. “Spliceosomal complex”, “nuclear speck”, “catalytic step 2 spliceosome” were the three most significant cellular component terms. For molecular function, “mRNA binding”, “helicase activity” and “mRNA 3’- UTR binding” were three most enriched categories (Figure 8A). KEGG analysis revealed a total of 5 remarkably enriched pathways, including “spliceosome”, “RNA transport”, “mRNA surveillance pathway”, “RNA degradation”, and “Legionellosis” (Figure 8B).

Table 1. Clinicopathological characteristic of 92 patients with adrenocortical carcinoma from TCGA database.
Clinicopathological characteristicsValue
Gender, n(%)
Female60(65.2)
Male32(34.8)
TCGA stage, n(%)
Stage I9(9.7)
Stage II44(47.8)
Stage III19(20.7)
Stage IV18(19.6)
Unknown2(2.2)
T stage, n(%)
T19(9.7)
T249(53.3)
T311(12)
T421(22.8)
Unknown2(2.2)
N stage, n(%)
N080(87.0)
N110(10.8)
Unknown2(2.2)
Survival, n(%)
Yes59(64.1)
No33(35.9)

Construction of potential splicing regulatory network

Then, a correlation network between the expression of SFs and the PSI values of survival-associated AS events was conducted. As shown in Figure 8C, there were a total of 42 down-regulated AS events (purple triangles), 50 up-regulated AS events (yellow triangles) and 67 SFs in this network. In this network, we selected the five most significant nodes as the hub AS events or hub SFs according to the degree, including one upregulated AS event (UQCR11-46528-AT), two downregulated AS events (UQCR11-46527-AT and NASP-2754-ES) and two SF (EIF3A and HNRNPR).

DISCUSSION

It has been reported that there was an extremely high incidence of splicing defects in carcinogenesis and that RNA splicing modulators might be a promising target for cancer treatment [13-15]. Alternative splicing allows cells to generate diverse [13, 14] mRNA by modifying mRNA isoforms [16]. There have been numerous studies exploring the important role of AS in

various malignant tumors, including colorectal cancer [17], renal cell carcinoma [18], breast cancer [19] and so on. Gray et al. [20]. found a novel splicing event which showed both an allelic expression imbalance and preferential splicing for one of the alleles in ACC samples. Previous studies demonstrated that the activation of spliced X-box protein 1(XBP1) - mRNA splicing might be one of the important molecular mechanisms of mitotane in the treatment of ACC [21- 23]. Besides, Kroiss et al. [22]. found that combination of bortezomib and carfilzomib with low-dose mitotane could be used to treat ACC effectively by increasing XBP1-mRNA splicing. Freddi et al. [24]. explored the expression of splice variants of growth hormone- releasing hormone receptor (GHRH-R) in human primary adrenocortical tumors and found that splice variants might be a key pathogenic factor in adrenal carcinogenesis. However, limited studies focused on the AS events in occurrence and progression of ACC.

As far as we know, this is the first study to date lucubrating the aberrant alternative spliced variants in patients with ACC using high throughput data. In this study, we systematically explored the prognostic value

Figure 1. Upset plot of interactions among seven types of all AS events (A) and survival-associated AS events (B) in ACC.

A

50G

1500

403

Gene Intersections

1000

23

AS

500

51

aly19

1848473,64, 50, 0

40 \1988, Nº 14 10,00,00 g1 st 18 76 14 75 g8 60 68 gb 00 60 63 60 gb g0 at 15 a4 30 g8 35 33 3 3 31 25 0 20 21 9

0

.

ME

RI

.

AD

AA

AP

AT

.

ES

.

4000

2000

0

Set Size

199

B

750

ar

Gene Intersections

500

346

250

MAB,

04

0 10 20

0

ME

AA

RI

AD

AP

AT

ES

1000

750

500

250

o

Set Size

Figure 2. Bubble plots of top 20 significant prognostic AS events in AA (A), AD (B), ES (C), AP (E), AT (F), ME (G) and RI (H) type, respectively. Volcano plots of prognostic AS events (D). AA, alternate acceptor site; AD, alternate donor site; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron.

A

ZFAND6-32171-AA

TRAFD1-24580-AD

CIRBP-46443-ES

MED11-38579-AA

B

ZSCAN18-52407-AD

C

CIRBP-46445-ES

D

-

WASH4P-32782-AA

BEX2-89726-AD

HM13-58892-ES

Prognosis-related AS

No significant

FCF1-28423-AA

pvalue

UQCRQ-73323-AD

pvalue

MPND-46796-ES

ZNF682-10567-AA

60-04

NOP2-19895-AD

40-04

NDUFA12-23740-ES

pvalue

POURZH-67943-AA

GOLGA3-25307-AD

36-04

FLAD1-7858-ES

50-06

SRP19-72988-AA

40-04

NDUFAF2-72172-AD

20-04

STOML1-31622-ES

48-06

USMG5-13004-ES

-log10(pvalue)

PEX10-266-AA

PRKAG1-21510-AD

30-06

STRADA-42966-AA

20-04

RPS8-214603-AD

1e-04

CTRL-20076-ES

20-06

PHYKPL-74856-AA

UBL7-31724-AD

10-06

*

GGCX-54288-ES

GUK1-10188-AA

GTF3C1-35695-AD

CMC1-63790-ES

MYL6-22381-AA

-log10(pvalue)

HTATIP2-14714-AD

-log10(pvalue)

PCBP2-22056-ES

RHOC-4238-AA

-log10(pvalue)

3.5

SLC35F5-55071-AD

4

MYL6-22384-ES

2

CIRBP-46427-AA

4.0

ANAPC11-44212-AD

5

SEC24C-12227-ES

6

ACTRS-23907-AA

4.5

PPP1R8-1347-AD

6

SGSM3-62347-ES

7

CER55-21680-AA

5.0

ANAPC11-44222-AD

7

MBD1-45620-ES

B

FDXR-43335-AA

SRSF5-28161-AD

8

MMAA-70764-ES

O

RUSC1-8079-AA

FOXR-43338-AD

PSEN2-10033-ES

-4

-2

0

2

4

ASB8-21443-AA

DPY30-53146-AD

CDK5-82324-ES

IKBKB-83594-AA

ABCCS-67824-AD

.

TANK-55743-ES

2-score

-5.0

2.5

0.0

2.5

-6.0

-2.5

0.0

25

5.0

-3

Q

3

6

z-score

z-score

z-score

BLOC151-22229-AP

METTL15-14782-AT

DIF6-59080-P

E

UNG-24278-AP

F

KLHL7-78950-AT

THNSL2-54469-ME

CIREP-46441-P

UNG-24277-AP

KLHL7-78952-AT

G

H

FASTK-82334-FI

CMC2-37702-AP

pvalue

TECPR2-29416-AT

pvalue

ANGEL2-9778-ME

pvalue

TST-62070-P

TECPR2-29415-AT

DYNLL1-24763-RI

pvalue

CMC2-37703-AP

DUT-30484-AP

DNAJC12-11898-AT

1.50-06

0.03

PILRB-80936-RI

0.00020

DUT-30485-AP

20-05

DNAJC12-11899-ATT

RAB6A-17707-ME

0.00015

KIAA0391-27213-AP

SNX5-58744-AT

1.00-06

0.02

ZSWIM7-38393-FI

ARLGIP4-25024-RI

PORMC2-70575-AP

1e-05

0.00010

UCHLS-9240-AT

5.00-07

0.01

CLN3-35746-RI

0.00005

PIK2A-12729-AP

ABCA3-33266-AT

H2AFY-96831-ME

GSTK1-82082-RI

PHK2A-12728-AP

ABCA3-33265-AT

FASTK-82333-RI

MRPL51-19869-AP

-log10(pvalue)

USP4-64852-AT

-log10(pvalue)

GCDH-47884-ME

-log10(pvalue)

AKAPEL-48077-RI

-log10(pvalue)

MRPLS1-19868-AP

5

AIG1-77972-AT

6.0

2

Clorf50-2117-RI

4

GHDC-41016-AP

6

DFNB31-87326-AT

6.5

RADS1-30020-ME

3

C16orf91-33097-RI

5

GHDC-41017-AP

OKR1-84849-AT

CCDC107-86267-RI

6

PSMG3-78500-AP

7

7.0

8

UCHLS-9239-AT

4

NFE2L1-42158-AP

ARHOEF28-72492-AT

7.5

EMD1-61575-ME

5

ATPSD-46401-RI

PSMD11-40209-R

7

PSMA1-14471-AP

ARHGEF28-72493-AT

TMEM91-50046-RI

SUMO2-43379-AP

KLHL3-73481-AT

CCDC53-24021-ME

EPOR-47689-R

LDB1-12935-AP

.

KLHL3-73482-AT

.

USP39-54316-RI

·

-6

-3

0

3

6

-3

8

3

6

-2

0

2

-25

0.0

25

5.0

z-score

Z-score

z-score

z-score

Table 2. Multivariate Cox analysis of prognostic alternative splicing predicting overall survival.
gene symbolSpliceseq IDAS typecoefHRHR.95LHR.95Hpvalue
MEWASH4P32782AA-14.320720256.03E-074.69E-100.0007765078.85E-05
FCF128423AA10.5224208537138.9233926.6719393951713510.980.004385308
ZNF69210567AA-5.1565314630.005761650.0004815140.06894224.66E-05
STRADA42966AA-49.809837572.33E-224.63E-361.17E-080.001972983
PHYKPL74856AA9.0806547618783.7153770.268598401287245402.70.086875639
RHOC4238AA-23.487149556.30E-119.57E-210.4152814830.041735412
CIRBP46427AA-11.514328839.99E-061.55E-080.0064244440.000483316
FDXR43335AA-89.669821981.14E-391.11E-571.17E-212.26E-05
ASB821443AA-8.0542953930.0003177349.34E-070.108033020.006764598
IKBKB83594AA-12.554128813.53E-066.50E-110.1916602670.024010337
NOP219895AD5.78087704324.043264721.797217574817.3138182.69E-05
PRKAG121510AD7.9423085572813.84894835.46032682223284.62850.000372286
RPS6214603AD10.7965241848850.7093158.2105428840995869.170.001671678
UBL731724AD5.487061825241.54645764.52072576112906.045230.006867221
SLC35F555071AD9.40523433912151.821354.518077662708583.4140.000650922
DPY3053146AD13.831778341016400.8180.1197147778.62943E+120.089281102
ABCC567824AD-2.2977457890.1004851030.0147851870.6829305610.018773482
ZSCAN1852407AD12.12172017183821.44825.9565541913018031780.007364346
DUT30484AP-5.6200226040.0036245597.82E-050.1680120620.004088329
PGRMC270575AP6.70571319817.06053846.3319200514408.811944.66E-06
PSMG378590AP11.3083314381497.81036662.622971310023638.454.11E-06
PSMA114471AP-24.414348952.49E-115.66E-190.001099770.006556697
METTL1514782AT-29.333215871.82E-131.85E-201.79E-060.000356241
DNAJC1211898AT9.2648008210559.707130.507740663655053.1610.001898066
USP464852AT16.056096029398832.355334.09486812.6441E+110.002127869
KLHL373481AT-2.2593905660.1044140990.0059814251.8226933350.121496101
STOML131622ES-9.1600949290.0001051534.91E-070.0225078910.000820894
USMG513004ES35.971100894.18842E+152920100.9496.01E+240.000826184
C1RL20076ES-28.403539214.62E-136.34E-193.37E-073.72E-05
GGCX54288ES-9.1371412750.0001075945.07E-070.0228193390.000828759
MMAA70764ES-21.661726963.91E-105.37E-142.85E-061.81E-06
PSEN210033ES-23.240415318.07E-112.43E-200.2679747540.037737792
THNSL254469ME
FASTK82334RI15.249655354196054.9731668.459188105527767550.000134968
TST62070RI10.7925880848658.806260.097705043242329295470.106858488
PILRB80936RI-5.6254771120.0036048430.0001143490.1136424110.001397614
PSMD1140209RI-5.4256129370.0044023672.48E-067.815807430.155222621
AllCIRBP46443ES13.996959931198953.8320.0675357312.12849E+130.100277799
BLOC1S122229AP-25.727920576.71E-123.20E-200.0014043830.008491906
TRAFD124580AD6.050752744424.43239871.54910273116288.51830.034618673
METTL1514782AT-23.564122315.84E-116.33E-190.0053854260.011794019
HM1358892ES-1.9574543110.141217460.0112512751.7724543580.129386117
Figure 3. Lasso regression analysis of survival-associated AS events. AA cohort (A, B), AD cohort (C, D), ES cohort (E, F), the whole cohort (G, H), AP cohort (I, J), AT cohort (K, L), ME cohort (M, N) and RI cohort (O, P).

15

15

15

15

16

17

17

16

15

15

8

0

16

16

16

14

14

13

13

14

14

12

10

7

4

1

0

13

13

13

12

11

11

11

11

12

10

9

5

3

0

13

12

9

3

U

6

6

6

6

6

6

6

6

Đ

6

1 T 5 1

A

4.5

C

E

10.5

G

T

4

2

10.0

Partial Likelihood Deviance

9.0

Partial Likelihood Deviance

2

Partial Likelihood Deviance

Partial Likelihood Deviance

8

2

5

35

9.0

g

8.0

2

2

9

13

0

3

1

1

-3.0

-2.5

-2.0

-1.5

-3.5

-3.0

-2.5

-2.0

-1.5

-3.5

-3.0

-25

-2.0

-1.5

-1.0

-4.5

-4.0

-3.5

-3.0

-2.5

-2.0

-1.5

-1.0

log(Lambda)

log(Lambda)

log(Lambda)

log(Lambda)

15

17

16

8

16

14

13

12

5

13

13

11

11

9

0

13

10

8

6

6

6

7

0

B

D

1

F

H

8

0

10

5

F

2

0

Coefficients

Coefficients

0

Coefficients

Coefficients

#

0

,

-10

:

2

&

e

$

8

-3.0

-2.5

-2.0

-1.5

-3.5

-3.0

-2.5

-2.0

-1.5

-3.5

-3.0

-2.5

-2.0

-1.5

-1.0

=4.5

-4.0

=3.5

=3.0

-2.5

-2.0

=1.5

=1.0

Log Lambda

Log Lambda

Log Lambda

Log Lambda

13

12

13

12

13

13

11

11

11

11

12

8

4

0

14

14

12

12

10

9

9

8

9

9

9

8

5 0

8

8

8

7

7

7

1

7

7

7

7

7

7

7

4

3

1

0

16

15

15

16

16

15

15

15

14

14

14

12 9

5

3

0

10.0

K

M

¥

O

=

5

10.0

?

Partial Likelihood Deviance

Partial Likelihood Deviance

=

9.0

9.5

Partial Likelihood Deviance

Partial Likelihood Deviance

2

2

9.0

8.5

=

=

:

:

80

1

1

0

0

-3.5

-3.0

-2.5

-2.0

-1,5

-1.0

-5

-4

-3

-2

-1

-6

-5

-4

-3

-2

-4.0

-3.5

-3.0

-2.5

-2.0

-1.5

log(Lambda)

log(Lambda)

log(Lambda)

log(Lambda)

13

11

13

11

11

0

14

12

8

9

0

8

7

7

7

4

16

15

15

15

14

5

J

L

N

P

:

0

45

10

1

0

&

1

0

0

Coefficients

Coefficients

Coefficients

1

Coefficients

4

1

-40

=

=

8

10

£

2

&

-3.5

-3.0

-25

-2.0

-1.5

-1.0

-5

-4

-3

-2

-1

-6

-5

-4

-3

-2

Log Lambda

Log Lambda

Log Lambda

-4.0

-3.5

-3.0

-2.5

-2.0

-1.5

Log Lambda

of AS events in ACC patients for the first time. A total of 3919 AS events, which were significantly associated with survival outcome, were identified by univariate Cox regression analysis. Next, Lasso regression and multivariate Cox regression

analysis were preformed to constructed eight the predictive significances, including seven types of AS and all AS, respectively. These ACC patients from TCGA database were then divided into low- and high-risk groups according to the risk score.

Figure 4. The Kaplan-Meier curves analysis of overall survival between the low-risk patients and high-risk patients. AA cohort (A), AD cohort (B), ES cohort (C), the whole cohort (D), AP cohort (E), AT cohort (F), ME cohort (G) and RI cohort (H).

A

B

C

D

Risk

High risk

Low risk

Risk

High risk

Low risk

Risk

High risk

Low risk

Risk

High risk

Low risk

1.00-

1.00-

1.00

1.00-

A

Survival probability

0.75

Survival probability

0.75

Survival probability

0.75

Survival probability

0.75

0.50

0.50

0.50

0.50

0.25

p=5.121e-10

0.25

p=5.9088-08

0.25

p=6.684e-11

0.25

p=2.39e-08

0.00

0.00

0.00

0.00

0

1

2

3

4

5

6

7

B

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

0

1

2

3

4

5

6

7

8

9

10

11 1

12

Time(years)

Time(years)

Time(years)

Time(years)

AA

AD

ES

All

Risk

High risk Low risk

3

19

0

10

7

5

A

2

2

Risk

High risk-

38

34

7

1

A

ON

UN

Risk

High risk

38

34

33

16

3

0

10

2

4

2

g

Risk

High risk-

Low risk

39

34 39

13

Low risk

39

39

15

Low risk

39

25

22

19

13

27

6 22

4

La

3

18

11

1

3

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

0

1

2

3

4

5

6

7

8

9

10

11

12

Time(years)

Time(years)

Time(years)

Time(years)

E

F

G

H

Risk

High risk

Low risk,

Risk

High risk

Low risk

Risk

High risk

+

Low risk

Risk

High risk

Low risk

1.00

1.00

1.00

1.00-

Survival probability

0.75

Survival probability

0.75

Survival probability

0.75

Survival probability

0.75

0.50

0.50

0.50

0.50

0.25

P=6.622e-11

0.25

p=2.862e-07

0.25

p=1.1060-04

0.25

p=1.2460-07

0.00

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

5

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

1

8

9

Time(years)

10

11

12

Time(years)

Time(years)

Time(years)

AP

AT

ME

RI

Risk

High risk Low risk

34 39

33

7

0

9

g

8

2

Risk

High risk- Low risk

38

39 38

13

21

28

7

21 19

12

9

4

NO

NO

Risk

High risk- Low risk

1

Risk

1

39

0 4

High risk- Low risk

38 39

34 39

6 16

26

10

2

15

E

3

4

2

2

So

1

É

6

4

2

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

Time(years)

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

Time(years)

Time(years)

12

Time(years)

Figure 5. ROC curves to evaluate the predictive power of each prognostic signature. AA cohort (A), AD cohort (B), ES cohort (C), the whole cohort (D), AP cohort (E), AT cohort (F), ME cohort (G) and RI cohort (H).

A

ROC curve ( AUC = 0.969 )

B

ROC curve ( AUC = 0.983)

C

ROC curve ( AUC = 0.869 )

D

ROC curve ( AUC = 0.897)

0

1.0

1.0

0

0.8

0.8

0.8

0.8

True positive rate

0.6

True positive rate

0.6

True positive rate

True positive rate

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

AA

0.0

AD

0.0

ES

0.0

9

All

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

False positive rate

False positive rate

False positive rate

False positive rate

E

F

G

H

ROC curve ( AUC = 0.925 )

ROC curve ( AUC = 0.764 )

ROC curve ( AUC = 0.767 )

ROC curve ( AUC = 0.928 )

a

1.0

1.0

1.0

0.8

0.8

0.8

0.8

True positive rate

0.6

True positive rate

0.6

True positive rate

0.6

True positive rate

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

02

0.0

AP

0.0

AT

RI

0.0

ME

0.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

False positive rate

False positive rate

False positive rate

False positive rate

Figure 6. The distribution of risk score, the distribution of survival time, the expression Heatmap of the most four significant prognostic signatures with AUC > 0.9. AA cohort (A-C), AD cohort (D-F), AT cohort (G-I), and RI cohort (J-L).

A

AA

D

AD

1

2

00

-

Risk score

Risk score

6

.

.

2

2

0

0

0

20

40

60

80

0

20

40

60

80

Patients (increasing risk socre)

Patients (increasing risk socre)

B

E

Survival time (years)

2

Survival time (years)

:

-

8

-

=

2

0

0

..

0

20

40

60

80

0

20

40

60

80

Patients (increasing risk socre)

Patients (increasing risk socre)

C

F

type

type

type

type

FCF 1|28423|AA

0.8

high

low

ABCC5)67824|AD

0.8

high

low

IKBKB|83594JAA

0.6

PRKAG1|21510|AD

0.6

STRADAĮ42966JAA

0.4

0.4

FDXRĮ43335JAA

0.2

DPY30/53146|AD

0.2

PHYKPLĮ74856JAA

0

ZSCAN18|52407|AD

0

RHOC|4238|AA

RPS6|214603|AD

CIRBP|46427JAA

SLC35F5|55071|AD

AS88/21443JAA

WASH4P(32782JAA

NOP2|19895|AD

ZNF692|10567JAA

UBL7|31724|AD

AT

RI

G

J

10

1

-

8

Risk score

10

Risk score

4

«

2

2

0

0

0

20

40

60

80

0

20

40

60

80

Patients (increasing risk socre)

Patients (increasing risk socre)

H

K

~

~

Survival time (years)

2

Survival time (years)

0

00

®

0

6

.

.

2

0

0

0

20

40

60

80

0

20

40

60

80

Patients (increasing risk socre)

Patients (increasing risk socre)

I

type

type

L

type

type

0.8

high

high

low

0.6

low

DNAJC12|11896|AT

0.6

PILRB(80936|RI

0.4

0.4

USP464852|AT

0.2

TST|62070|RI

0.2

METTL15|14782|AT

FASTK(82334|RI

KLHL3|73481|AT

PSMD11|40209;R)

A

AA

pvalue

Hazard ratio

gender

0.928

0.963(0.425-2.180)

stage

<0.001

2.886(1.819-4.579)

T

<0.001

3.349(2.075-5.403)

N

<0.001

10.060(3.037-33.322)

riskScore

<0.001

1.001(1.001-1.002)

0.50

1.0

2.0

4.0

8.0

16.0

32.0

0.50

1.0

2.0

4.0

8.0

16.0

32.0

Hazard ratio

Hazard ratio

0.25

0.50

1.0

2.0

4.0

8.0

Hazard ratio

0.50

1.0

2.0

4.0

8.0

16.0

Hazard ratio

E

AT

G

RI

0.50

1.0

2.0

4.0

8.0

16.0

32.0

0.50

1.0

2.0

4.0

8.0

16.0

32.0

Hazard ratio

Hazard ratio

F

H

pvalue

Hazard ratio

gender

0.987

1.008(0.411-2.468)

stage

0.898

0.930(0.306-2.830)

T

0.048

2.965(1.008-8.726)

N

0.015

6.336(1.431-28.055)

riskScore

<0.001

1.068(1.031-1.107)

0.50

1.0

2.0

4.0

8.0

16.0

0.50

1.0

2.0

4.0

8.0

16.0

Hazard ratio

Hazard ratio

CAD
pvalueHazard ratio
gender0.9280.963(0.425-2.180)
stage<0.0012.886(1.819-4.579)
T<0.0013.349(2.075-5.403)
N<0.00110.060(3.037-33.322)
riskScore0.0011.002(1.001-1.003)
BpvalueHazard ratio
gender0.4320.695(0.280-1.723)
stage0.3691.490(0.625-3.552)
T0.1591.869(0.783-4.463)
N0.0683.637(0.910-14.540)
riskScore<0.0011.001(1.001-1.002)
DpvalueHazard ratio
gender0.8651.077(0.456-2.547)
stage0.8151.127(0.414-3.063)
T0.0412.753(1.043-7.265)
N0.0384.243(1.083-16.621)
riskScore<0.0011.003(1.001-1.004)
pvalueHazard ratiopvalueHazard ratio
gender0.9280.963(0.425-2.180)gender0.9280.963(0.425-2.180)
stage<0.0012.886(1.819-4.579)stage<0.0012.886(1.819-4.579)
T<0.0013.349(2.075-5.403)T<0.0013.349(2.075-5.403)
N<0.00110.060(3.037-33.322)N<0.00110.060(3.037-33.322)
riskScore<0.0011.095(1.060-1.130)riskScore<0.0011.081(1.050-1.113)
pvalueHazard ratio
gender0.8540.921(0.385-2.202)
stage0.9440.964(0.347-2.677)
T0.0422.846(1.039-7.794)
N0.0374.202(1.093-16.156)
riskScore<0.0011.080(1.044-1.117)

Figure 7. Independent prognostic analysis determining predictors of overall survival. AA cohort (A, B), AD cohort (C, D), AT cohort (E, F), and RI cohort (G, H).

Kaplan-Meier analysis demonstrated that the difference of overall survival between the low-risk patients and high-risk patients were significant. Besides, the results showed that AUC of AA, AD, AP and RI were more than 0.9, indicating that these four types of AS were of great significance in predicting survival outcome in ACC patients. As a consequence, these four types of AS were screened to further analysis. However, as is known to us, there were few studies related to AA, AD, AD and RI in ACC. Future study into these four AS types in ACC might have noticeable research value. Given the high prevalence of splicing defects in tumor, these prognostic AS events identified in this study could be novel therapeutic targets in ACC treatment.

As is known to us, splicing factors were one of the vital regulated factors of AS events [25]. Changes of expression or sequence mutations of SFs may affect splicing events [26]. To elucidate the underlying mechanism of AS events in the survival of ACC patients, we performed GO and KEGG enrichment analysis for the SFs related significantly to AS events. According to the results, “mRNA processing”, “RNA splicing” and “RNA splicing via transesterification reactions” were the three most significant biological

process. “Spliceosomal complex”, “nuclear speck”, “catalytic step 2 spliceosome” were the three most significant cellular component terms. For molecular function, “mRNA binding”, “helicase activity” and “mRNA 3’-UTR binding” were three most enriched categories. The GO results revealed that these SFs have a significant relationship with splicing. KEGG analysis demonstrated that spliceosome is the most significant pathway among these SFs. Besides, we also established an interaction network between AS events and SF. In this network, we selected the five most significant nodes as the hub AS events or hub SF according to the degree, including one upregulated AS event (UQCR11-46528-AT), two downregulated AS events (UQCR11-46527-AT and NASP-2754-ES) and two SF (EIF3A and HNRNPR). These new hub SF genes and AS events were not reported previously, which need further validation.

There were several limitations in this study. Firstly, we merely analyzed the survival-associated AS events of ACC using TCGA database due to no additional external cohort concerning splicing data. Secondly, this is a retrospective study with limited number of patients. Further study with larger sample size is warranted.

Figure 8. GO functional enrichment analysis (A) and KEGG pathway analysis (B) of AS events related SF genes. The interaction network between SF and AS events (C).

A

mRNA processing

C

RNA splicing

RNA splicing, via transesterification reactions

RNA splicing, via transesterification reactions with bulged adenosine as nucleophile

mRNA splicing, via spliceosome

regulation of mRNA metabolic process

PI4K24-1

regulation of RNA splicing

ribonucleoprotein complex assembly

regulation of mRNA processing

regulation of mRNA splicing, via spliceosome

Count

PRACA- 100

50

PRKONEY

spliceosomal complex

100

nuclear speck

150

PROFILO

catalytic step 2 spliceosome

200

U2-type spliceosomal complex

PSEN2

Sm-like protein family complex

250

PTART

U2-type precatalytic spliceosome

8

precatalytic spliceosome

p.adjust

small nuclear ribonucleoprotein complex

0.000000e+00

UQGR11-46528>A

spliceosomal snRNP complex

1.343717e-16

4.763-RJ

U2-type catalytic step 2 spliceosome

2.6874340-16

mRNA binding

4.0311520-16

helicase activity

5.3748804-16

mRNA 3’-UTR binding

EIF3A

ribonucleoprotein complex binding

pre-mRNA binding

single-stranded RNA binding

RNA helicase activity

snRNA binding

NASI 2754 ES

RNA-dependent ATPase activity

SLC41A

ATP-dependent RNA helicase activity

.

SLCA1ALLE

LØGRD-4602THAT

0.2

0.4 0.

GeneRatio

INFINE

B

SURF

INRNPR

SNXS-

Spliceosome

SPHKZ

SERY

p-adjust

RNA transport

0.004

0.00%

0.012

ST7-M

mRNA surveillance pathway

TAFTA

Count

35

TCEB2

50

RNA degradation

100

125

TRINEZES

Legionelosis

LAP

0.0

0.2

0.4

GeneRatio

0.6

CONCLUSIONS

In summary, this is the first and most comprehensive study to explore the aberrant alternative spliced variants in patients with ACC. These identified survival-related alternative spliced events and the interesting interaction network between AS and SF might provide novel insights into molecular mechanism of ACC. Further researches are needed to investigate the role of AS in ACC.

MATERIALS AND METHODS

Data collection and preprocessing

The transcriptome data and clinical information of the ACC cohort were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Percent Spliced In (PSI) values for AS events of 33 various tumor types and available adjacent normal samples, have been loaded into TCGA SpliceSeq data portal (https://bioinformatics.mdanderson.org/TCGA

SpliceSeg/) [13]. The AS-related data in TCGA SpliceSeq is of high quality [13]. The data, download procedure and operating environment in TCGA SpliceSeq have been widely used to explore survival-associated AS events and signatures of various malignant tumors [12, 13]. We obtained AS events with PSI value of ACC from the TCGA SpliceSeq data portal. There were seven types of AS events, including retained intron (RI), exon skip (ES), mutually exclusive exons (ME), alternate promoter (AP), alternate donor site (AD), alternate terminator (AT) and alternate acceptor site (AA).

Univariate Cox regression analysis was carried out to evaluate the association between AS events and overall survival (OS). UpSet plot and volcano Plot was used to present the survival-related AS events based on the seven types of AS events. Besides, the bubble chart was used to summarize the top 20 all AS events, top 20 RI events, top 20 ES events, top 8 ME events, top 20 AP events, top 20 AD events, top 20 AT events and top 20 AA events.

Prognostic risk scores calculation and survival analysis

Lasso regression analysis was further conducted to screen highly correlated AS events and avoid overfit the models constructed subsequently. Multivariate Cox regression analysis was applied to calculate the prognostic risk scores for OS prediction based on seven types of AS events. All patients were divided into high-risk group and low-risk group according to the median risk score. Kaplan-Meier curves were used to demonstrate the

survival probability variation between low-risk and high- risk subgroups. Receiver operating characteristic (ROC) curve was constructed and the area under the curve (AUC) was also calculated. Further, prognostic evaluation of the risk score was performed using risk curve. The expression heatmap, distribution of risk score and survival time of related signature were then visualized.

Independence of prognostic signature

The following clinicopathological characteristic from TCGA database including gender, T stage, N stage and clinical stage were subjected to subsequent analysis. Univariate and multivariate Cox regression analysis was performed to evaluate whether the final AS prognostic signature with an AUC ≥ 0.9 was an independent risk factor of OS.

Functional enrichment analyses

The list of splicing factor (SF) genes was collected from the SpliceAid 2 database (http://193.206.120.249/ splicing_tissue.html) [14] and was shown in Supplementary Table 1. The expression profiles of SF genes were downloaded from TCGA database. To further illustrated the underlying mechanisms of AS in ACC, we identified corresponding SF genes of AS events. Functional enrichment analyses, including Gene Ontology (GO) analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, for these SF genes was performed. The results of GO analysis were presented by three parts including biological processes (BP), molecular functions (MF), and cellular component (CC). Both GO analysis and KEGG analysis was conducted using R x64 3.6.1 software.

Construction of splicing regulatory network

Spearman test was conducted to analyze the correlation between the expression of SF genes and PSI values of survival-associated AS events. P value < 0.001 and correlation coefficient > 0.6 was considered statistically significant. Besides, we constructed an interaction network of AS and SF. The plug-in Molecular Complex Detection (MCODE) was used to select most significant modules from the interaction network. The most significant interaction network was then visualized by Cytoscape software (version 3.4.0). Furthermore, the top five hub nodes were selected according to the connection degree from this interaction network.

CONFLICTS OF INTEREST

The authors of this manuscript declare no conflicts of interest.

FUNDING

This study was supported by Foundation of Fujian Provincial Department of Finance (Grant number: 2018B011).

REFERENCES

1. Mpaili E, Moris D, Tsilimigras DI, Oikonomou D, Pawlik TM, Schizas D, Papalampros A, Felekouras E, Dimitroulis D. Laparoscopic Versus Open Adrenalectomy for Localized/Locally Advanced Primary Adrenocortical Carcinoma (ENSAT I-III) in Adults: Is Margin-Free Resection the Key Surgical Factor that Dictates Outcome? A Review of the Literature. J Laparoendosc Adv Surg Tech A. 2018; 28:408-14. https://doi.org/10.1089/lap.2017.0546 PMID: 29319399

2. Else T, Kim AC, Sabolch A, Raymond VM, Kandathil A, Caoili EM, Jolly S, Miller BS, Giordano TJ, Hammer GD. Adrenocortical carcinoma. Endocr Rev. 2014; 35:282-326. https://doi.org/10.1210/er.2013-1029 PMID: 24423978

3. Vanbrugghe C, Lowery AJ, Golffier C, Taieb D, Sebag F. Adrenocortical carcinoma surgery-surgical extent and approach. Langenbecks Arch Surg. 2016; 401:991-97. https://doi.org/10.1007/s00423-016-1462-8 PMID: 27412357

4. Paragliola RM, Torino F, Papi G, Locantore P, Pontecorvi A, Corsello SM. Role of Mitotane in Adrenocortical Carcinoma - Review and State of the art. Eur Endocrinol. 2018; 14:62-66. https://doi.org/10.17925/EE.2018.14.2.62 PMID: 30349596

5. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, Kingsmore SF, Schroth GP, Burge CB. Alternative isoform regulation in human tissue transcriptomes. Nature. 2008; 456:470-76. https://doi.org/10.1038/nature07509 PMID: 18978772

6. Zhang D, Duan Y, Cun J, Yang Q. Identification of Prognostic Alternative Splicing Signature in Breast Carcinoma. Front Genet. 2019; 10:278. https://doi.org/10.3389/fgene.2019.00278 PMID: 30984247

7. Baralle FE, Giudice J. Alternative splicing as a regulator of development and tissue identity. Nat Rev Mol Cell Biol. 2017; 18:437-51. https://doi.org/10.1038/nrm.2017.27 PMID: 28488700

8. Fu XD, Ares M Jr. Context-dependent control of alternative splicing by RNA-binding proteins. Nat Rev Genet. 2014; 15:689-701.

https://doi.org/10.1038/nrg3778 PMID: 25112293

9. Park E, Pan Z, Zhang Z, Lin L, Xing Y. The Expanding Landscape of Alternative Splicing Variation in Human Populations. Am J Hum Genet. 2018; 102:11-26. https://doi.org/10.1016/j.ajhg.2017.11.002 PMID: 29304370

10. Shen S, Wang Y, Wang C, Wu YN, Xing Y. SURVIV for survival analysis of mRNA isoform variation. Nat Commun. 2016; 7:11548. https://doi.org/10.1038/ncomms11548 PMID: 27279334

11. Climente-González H, Porta-Pardo E, Godzik A, Eyras E. The Functional Impact of Alternative Splicing in Cancer. Cell Rep. 2017; 20:2215-26. https://doi.org/10.1016/j.celrep.2017.08.012 PMID: 28854369

12. Hu YX, Zheng MJ, Zhang WC, Li X, Gou R, Nie X, Liu Q, Hao YY, Liu JJ, Lin B. Systematic profiling of alternative splicing signature reveals prognostic predictor for cervical cancer. J Transl Med. 2019; 17:379. https://doi.org/10.1186/s12967-019-02140-x PMID: 31744495

13. Ryan M, Wong WC, Brown R, Akbani R, Su X, Broom B, Melott J, Weinstein J. TCGASpliceSeq a compendium of alternative mRNA splicing in cancer. Nucleic Acids Res. 2016; 44:D1018-22. https://doi.org/10.1093/nar/gkv1288 PMID: 26602693

14. Piva F, Giulietti M, Burini AB, Principato G. SpliceAid 2: a database of human splicing factors expression data and RNA target motifs. Hum Mutat. 2012; 33:81-85. https://doi.org/10.1002/humu.21609 PMID: 21922594

15. Salton M, Misteli T. Small Molecule Modulators of Pre- mRNA Splicing in Cancer Therapy. Trends Mol Med. 2016; 22:28-37. https://doi.org/10.1016/j.molmed.2015.11.005 PMID: 26700537

16. Chen J, Weiss WA. Alternative splicing in cancer: implications for biology and therapy. Oncogene. 2015; 34:1-14. https://doi.org/10.1038/onc.2013.570 PMID: 24441040

17. Sakuma K, Sasaki E, Kimura K, Komori K, Shimizu Y, Yatabe Y, Aoki M. HNRNPLL, a newly identified colorectal cancer metastasis suppressor, modulates alternative splicing of CD44 during epithelial- mesenchymal transition. Gut. 2018; 67:1103-11. https://doi.org/10.1136/gutjnl-2016-312927 PMID: 28360095

18. Chen K, Xiao H, Zeng J, Yu G, Zhou H, Huang C, Yao W, Xiao W, Hu J, Guan W, Wu L, Huang J, Huang Q, et al. Alternative Splicing of EZH2 pre-mRNA by SF3B3

Contributes to the Tumorigenic Potential of Renal Cancer. Clin Cancer Res. 2017; 23:3428-41. https://doi.org/10.1158/1078-0432.CCR-16-2020 PMID: 27879367

19. Anczuków O, Akerman M, Cléry A, Wu J, Shen C, Shirole NH, Raimer A, Sun S, Jensen MA, Hua Y, Allain FH, Krainer AR. SRSF1-Regulated Alternative Splicing in Breast Cancer. Mol Cell. 2015; 60:105-17. https://doi.org/10.1016/j.molcel.2015.09.005 PMID: 26431027

20. Gray SG, Kjellman M, Larsson C, Ekström TJ. Novel splicing of an IGF2 polymorphic region in human adrenocortical carcinomas. Biochem Biophys Res Commun. 1997; 239:878-83. https://doi.org/10.1006/bbrc.1997.7554 PMID: 9367863

21. Sbiera S, Leich E, Liebisch G, Sbiera I, Schirbel A, Wiemer L, Matysik S, Eckhardt C, Gardill F, Gehl A, Kendl S, Weigand I, Bala M, et al. Mitotane Inhibits Sterol-O-Acyl Transferase 1 Triggering Lipid-Mediated Endoplasmic Reticulum Stress and Apoptosis in Adrenocortical Carcinoma Cells. Endocrinology. 2015; 156:3895-908. https://doi.org/10.1210/en.2015-1367 PMID: 26305886

22. Kroiss M, Sbiera S, Kendl S, Kurlbaum M, Fassnacht M. Drug Synergism of Proteasome Inhibitors and Mitotane by Complementary Activation of ER Stress in Adrenocortical Carcinoma Cells. Horm Cancer. 2016;

7:345-55. https://doi.org/10.1007/s12672-016-0273-2 PMID: 27631436

23. Sbiera S, Kendl S, Weigand I, Sbiera I, Fassnacht M, Kroiss M. Hsp90 inhibition in adrenocortical carcinoma: limited drug synergism with mitotane. Mol Cell Endocrinol. 2019; 480:36-41. https://doi.org/10.1016/j.mce.2018.10.009 PMID: 30315857

24. Freddi S, Arnaldi G, Fazioli F, Scarpelli M, Appolloni G, Mancini T, Kola B, Bertagna X, Mantero F, Collu R, Boscaro M. Expression of growth hormone-releasing hormone receptor splicing variants in human primary adrenocortical tumours. Clin Endocrinol (Oxf). 2005; 62:533-38.

https://doi.org/10.1111/j.1365-2265.2005.02253.x PMID: 15853821

25. Wahl MC, Will CL, Lührmann R. The spliceosome: design principles of a dynamic RNP machine. Cell. 2009; 136:701-18. https://doi.org/10.1016/j.cell.2009.02.009 PMID: 19239890

26. Dvinge H, Kim E, Abdel-Wahab O, Bradley RK. RNA splicing factors as oncoproteins and tumour suppressors. Nat Rev Cancer. 2016; 16:413-30. https://doi.org/10.1038/nrc.2016.51 PMID: 27282250

SUPPLEMENTARY MATERIALS

Supplementary Table

Please browse Full Text version to see the data of Supplementary Table 1. The list of splicing factor (SF) genes.