Identification of hub genes and pathways in adrenocortical carcinoma by integrated bioinformatic analysis

Jinshuai Guo1 İD Yinzhong Gu 2

İD Xiaoyu Ma1 1

İD Lu Zhang1 1 Huimin Li1 D İD |

Zhongyi Yan1 1

İD Yali Han’ 1

İD Longxiang Xie 1

İD

Xiangqian Guo1 İD

1Department of Predictive Medicine, Institute of Biomedical Informatics, Cell Signal Transduction Laboratory, Bioinformatics Center, Henan Provincial Engineering Center for Tumor Molecular Medicine, School of Basic Medical Sciences, Henan University, Kaifeng, China 2Education Bureau of Longgang District, Shenzhen, China

Correspondence

Longxiang Xie and Xiangqian Guo, Department of Predictive Medicine, Institute of Biomedical Informatics, Cell Signal Transduction Laboratory, Bioinformatics Center, Henan Provincial Engineering Center for Tumor Molecular Medicine, School of Basic Medical Sciences, Henan University, Kaifeng 475004, China. Emails: xielongxiang123@126.com (LX); xqguo@henu.edu.cn (XG)

Funding information

Kaifeng Science and Technology Major Project, Grant/Award Number: 18ZD008; Program for Science and Technology Development in Henan Province, Grant/Award Number: 162102310391; 172102210187; China Postdoctoral Science Foundation, Grant/Award Number: 2017M62237; Yellow River Scholar Program, Grant/Award Number: H2016012; Projects for College Students in Henan University, Grant/Award Number: 201819002; Supporting grants of Henan University, Grant/Award Number: 2015YBZR048; B2015151; Key Project of Science and Technology Research of Henan Provincial Department of Education, Grant/Award Number: 20A180001; National Natural

Abstract

Adrenocortical carcinoma (ACC), a rare malignant neoplasm originating from adrenal cortical cells, has high malignancy and few treatments. Therefore, it is necessary to explore the molecular mechanism of tumorigenesis, screen and verify potential bio- markers, which will provide new clues for the treatment and diagnosis of ACC. In this paper, three gene expression profiles (GSE10927, GSE12368 and GSE90713) were downloaded from the Gene Expression Omnibus (GEO) database. Differentially ex- pressed genes (DEGs) were obtained using the Limma package. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched by DAVID. Protein-protein interaction (PPI) network was evaluated by STRING database, and PPI network was constructed by Cytoscape. Finally, GEPIA was used to validate hub genes’ expression. Compared with normal adrenal tissues, 74 up-regulated DEGs and 126 down-regulated DEGs were found in ACC samples; GO analysis showed that up-regulated DEGs were enriched in organelle fission, nuclear division, spindle, et al, while down-regulated DEGs were enriched in angiogenesis, proteinaceous ex- tracellular matrix and growth factor activity; KEGG pathway analysis showed that up-regulated DEGs were significantly enriched in cell cycle, cellular senescence and progesterone-mediated oocyte maturation; Nine hub genes (CCNB1, CDK1, TOP2A, CCNA2, CDKN3, MAD2L1, RACGAP1, BUB1 and CCNB2) were identified by PPI net- work; ACC patients with high expression of 9 hub genes were all associated with worse overall survival (OS). These hub genes and pathways might be involved in the tumorigenesis, which will offer the opportunities to develop the new therapeutic tar- gets of ACC.

Jinshuai Guo and Yingzhong Gu are Co-first authors.

@ 2020 The Authors. Journal of Cellular and Molecular Medicine published by Foundation for Cellular and Molecular Medicine and John Wiley & Sons Ltd

Science Foundation of China, Grant/ Award Number: 81602362; Program for Scientific and Technological Research of Henan Education Department, Grant/ Award Number: 14B520022; Program for Innovative Talents of Science and Technology in Henan Province, Grant/ Award Number: 18HASTIT048; Henan Postdoctoral Foundation, Grant/Award Number: 001702052; Supporting grant of Bioinformatics Center of Henan University, Grant/Award Number: 2018YLJC01

KEYWORDS

adrenocortical carcinoma, Gene Expression Omnibus, GEPIA, overall survival, prognostic

1 INTRODUCTION |

Adrenal gland is an important endocrine organ of the human body and is also one of the most common organs with high tumour me- tastasis rate. According to the latest edition of pathological and genetic classification criteria of adrenal tumours published by World Health Organization (WHO), adrenal tumours can be di- vided into five categories: adrenocortical tumour, adrenal medul- lary tumour, extra-adrenal paraganglioma, secondary tumour and other adrenal tumour.1 Adrenocortical tumour mainly includes adrenocortical adenoma (ACA) and adrenocortical carcinoma (ACC).2 ACC is a rare malignant tumour originating from adre- nal cortical cells,3 and its incidence ranges from 0.7/10 00 000 to 2.0/10 00 000, which is under a high degree of malignancy, aggressiveness, high recurrence rate and poor prognosis. Most of the patients are found to have metastases and relapse easily after treatment, and the overall 5-year survival rate was <35%.4 Early and accurate diagnosis is particularly important for the treatment and prognosis of ACC.5 At present, surgical resection is the only feasible method to cure ACC, but it is difficult to con- trol its quality.6 Therefore, identifying new therapeutic targets or biomarkers for prognosis, diagnosis or prediction of ACC is urgently needed.

In recently years, many microarray profiling studies have been performed in ACC,7,8 and hundreds of differentially expressed genes (DEGs) have been obtained. However, the results are lim- ited or inconsistent due to molecular heterogeneity, and the re- sults are usually generated from a single cohort study. Until now, no reliable biomarkers have been used in ACC clinics. Hence, the bioinformatics methods integrating multi-cohorts analysed by gene microarray or RNAseq will be innovative and valuable for ACC research.

In this work, we downloaded three different Gene Expression Omnibus (GEO) datasets (GSE10927, GSE12368 and GSE90713) and screened differentially DEGs using the Limma package. Then, the PPI network in STRING database was constructed to screen the hub genes and pathways, and GEPIA database was used to verify hub genes and potential pathways. This study will offer the opportunities to develop the new therapeutic targets of ACC.

2 MATERIALS AND METHODS |

The flow diagram of this study was shown in Figure 1. The raw expression data were operated through a series of databases and software. This study has been proved by the Henan University insti- tutional committee.

FIGURE 1 The flow diagram of this study

GEO

Download of GSE10927 GSE12368 GSE90713

Limma

Identification of DEGs

DAVID

STRING

GO analysis

PPI network

KEGG

Cytoscape

Pathway analysis

Hub genes analysis

TCGA

GEPIA

Validation

2.1 Data collection

Gene expression profiles of GSE10927,9 GSE123681º and GSE9071311 were obtained from GEO database. The GSE10927 dataset includes 55 neoplastic samples and 10 non-neoplastic sam- ples (33 cases of ACC, 22 cases of ACA and 10 cases of normal ad- renal cortex). The GSE12368 dataset is composed of 28 neoplastic samples and 6 non-neoplastic samples (12 cases of ACC, 16 cases of ACA and 6 cases of normal adrenal cortex). The GSE90713 dataset includes 58 neoplastic samples and 5 non-neoplastic samples (58 cases of ACC and 5 cases of normal adrenal cortex).

2.2 Screening differentially expressed genes

DEGs between ACC samples and non-neoplastic samples were screened by using Limma package based on R language. DEGs were defined as representing differences with |log2FC| > 1, P < . 05.12

2.3 Gene ontology and KEGG pathway enrichment analysis of DEGs

Gene ontology (GO) analysis annotates genes and gene products with functions of molecular function (MF), biological process (BP) and cellular component (CC).13,14 Kyoto Encyclopedia of Genes and Genomes (KEGG) includes a series of genomics and enzymol- ogy methods and an online database of biochemical energy.15 The Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) is an online program providing a comprehensive set of functional annotation tools for investigators to understand biological meaning behind large list of genes.16 We performed GO terms and KEGG pathway analysis of DEGs by using DAVID database.

2.4 PPI network construction and hub module selection

Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/) was used to evaluate protein-protein in- teraction (PPI).17 In addition, the database was used to quantify the relationships among the DEGs. Then, we used Cytoscape software to construct PPI network.18 The genes with the highest node score and the strongest connectivity were selected. P < . 05 was consid- ered to have statistical significance.

2.5 Hub genes validation

GEPIA (http://gepia.cancer-pku.cn/) is a powerful interactive web server that can analyse the RNA sequencing expression data of 9736 tumours and 8587 normal samples from the TCGA and the GTEx

projects.19 GEPIA can be directly used for tumour/normal differential expression analysis according to cancer types, and the box plot will be shown to visualize the relationship. GEPIA was used to verify the hub genes and perform validation, P < . 05 showed statistical significance.19

TCGA-ACC RNA sequencing data with patient survival data were downloaded from the University of California, Santa Cruz (UCSC) Xena browser.20 Clinicopathological parameters of ACC patients with primary tumours, including age at diagnosis, gender, pathologic stage, living status, and overall survival (OS), were used for surviv- al-curve analysis.

2.6 Statistical analyses

Clinicopathologic parameter association analysis, and univariate and multivariate Cox regression analysis of 9 hub genes were performed with SPSS 22.0. Statistical significance was set at probability values of P < . 05.

3 RESULTS

|

3.1 Identification of differentially expressed genes

Using P < . 05 and |log2FC| > 1 as cut-off, we identified 1040 up- regulated and 1821 down-regulated genes in ACCs compared with normal tissues from GSE10927 dataset (Figure S1A), 2295 up-regulated and 3300 down-regulated genes from GSE12368 dataset (Figure S1B), 487 up-regulated and 1129 down-regulated genes from GSE90713 dataset (Figure S1C). By intersecting the DEGs across three datasets, a total of 200 consistently differen- tially expressed genes including 74 up-regulated and 126 down- regulated genes were identified to be significant in all the three above gene expression profiles (Figure 2A and B). The heat map of top 20 down-regulated genes and top 20 up-regulated genes’ expression was shown in Figure 2C.

3.2 | Gene ontology analysis of differentially expressed genes

Subsequently, GO analysis of up-regulated and down-regulated DEGs was carried out by using DAVID online analysis tool.16 In terms of biological process (BP), up-regulated DEGs were significantly en- riched organelle fission, nuclear division, mitotic nuclear division, chromosome segregation, regulation of cell cycle phase transition, regulation of mitotic cell cycle phase transition, nuclear chromo- some segregation and sister chromatid segregation (Figure 3A). Down-regulated DEGs were involved in angiogenesis, regulation of cell growth, fatty acid metabolic process, ageing, developmental maturation and protein kinase B signalling (Figure 3B). In terms of cellular component (CC), the up-regulated DEGs were significantly enriched in spindle, condensed chromosome, chromosomal region,

|

FIGURE 2 Venn diagram from intersection of differentially expressed genes (DEGs) in the three Gene Expression Omnibus (GEO) datasets and the heat map of top 40 DEGs. A, Up-regulated genes in ACCs. B, Down- regulated genes in ACCs. C, Heat map of top 40 representative DEGs. Blue represents a lower expression level, and red represents higher expression level

A

GSE10927

GSE12368

B

GSE10927

GSE12368

439

126

196

654

236

200

74

126

33

2

144

11

53

141

GSE90713

GSE90713

C

4.37

1.95

2.47

DEPDC1B

4.81

3.48

1.75

ASPM

4

4.32

3.97

2.97

1.89

DTL

2.50

1.79

BUB1

4.09

2.29

1.52

CENPK

2

3.64

3.20

1.50

TPX2

3.57

2.70

BIRC5

3.77

3.04

2.72

1.70

1.39

ANGPT2

0

3.57

2.52

PRC1

2.62

3.27

2.56

UBE2C

2.76

2.95

2.54

ENC1

-2

3.83

3.56

3.51

3.47

2.15

CCNB2

1.89

PBK

3.83

3.13

2.39

CDKN3

-4

3.61

3.24

2.29

ANLN

3.75

3.63

2.70

2.18

STC1

3.02

1.91

HMMR

-6

4.36

3.12

2.76

TOP2A

3.72

3.46

2.98

CDK1

3.51

3.65

2.90

KIAA0101

-6.03

-4.79

-4.10

ADH1B

-7.36

-4.20

-1.90

DAPL1

-6.34

-3.24

-1.70

C2orf40

-4.99

-1.31

-3.97

CYP11B2

-5.01

-2.31

-3.77

-2.62

MUM1L1

-3.17

-2.89

WFDC1

-3.48

-3.45

-3.52

INMT

-3.78

-3.62

-3.09

CXCL12

-3.64

-3.44

-3.12

SLC16A9

-3.71

-3.43

-2.93

KCNQ1

-3.93

-2.86

-3.25

C7

-3.16

-2.93

-3.75

STEAP4

-4.33

-4.06

-2.59

FMO2

-4.62

-4.36

-3.32

-2.66

ANGPTL1

-3.58

-2.48

OGN

-4.32

-3.57

-2.28

FBLN1

-4.97

-3.25

-3.33

NPY1R

-4.72

-2.86

-3.49

SLITRK4

-5.31

-3.81

-3.25

RSPO3

-4.81

-3.91

-3.33

PDGFRA

0

G

G

00

5

m

m

10

8

2

~

w

microtubule, chromosome and centromeric region (Figure 3C), while the down-regulated DEGs were significantly enriched in proteina- ceous extracellular matrix (Figure 3D). In terms of molecular func- tion (MF), up-regulated DEGs were significantly enriched in tubulin binding, microtubule binding, drug binding and cyclin-dependent protein kinase activity (Figure 3E), while the down-regulated DEGs were significantly enriched in growth factor activity, growth factor binding and hormone binding (Figure 3F).

3.3 | KEGG pathway enrichment analysis of differentially expressed genes

KEGG pathway analysis of all DEGs showed that most of up- regulated DEGs were enriched in cell cycle, cell senescence,

progesterone-mediated oocyte maturation, oocyte, p53 signalling pathway and folic acid resistance (Figure 4), while down-regulation of DEGs did not significantly enrich KEGG pathway.

3.4 | PPI network construction and hub gene selection

Through analysing STRING database17 and constructing PPI network by Cytoscape software,18 as shown in Figure 5, the PPI network of DEGs consists of 120 nodes and 1032 edges with the highest degree of 53, including 56 up-regulated genes and 64 down-regulated genes. It was considered that top nine DEGs with high degree of connectivity as the hub genes of ACC: CCNB1 (Cyclin B1), CDK1 (cyclin-dependent kinase 1), TOP2A

|

FIGURE 3 Gene ontology analysis of differentially expressed genes (DEGs) with |log_FC| > 1, P < . 05. A, Biological process terms for up-regulated DEGs. B, Biological process terms for down-regulated DEGs. C, Cellular component terms for up-regulated DEGs. D, Cellular component terms for down-regulated DEGs. E, Molecular function terms for up-regulated DEGs. F, Molecular function terms for down- regulated DEGs

A

B

Organelle fission

Angiogenesis

Nuclear division

Regulation of cell growth

Mitotic nuclear division

Fatty acid metabolic process

Chromosome segregation

Aging

Regulation of cell cycle phase transition

Developmental maturation

Regulation of mitotic cell cycle phase transition

p.adjust

Protein kinase B signaling

p.adjust

Nuclear chromosome segregation

1e-08

Respiratory system development-

0.04

Sister chromatid segregation

Respiratory tube development-

5e-09

0.03

Cell cycle checkpoint-

Lung development-

Mitotic sister chromatid segregation-

Phosphatidylinositol 3-kinase signaling

0.02

Mitotic cell cycle checkpoint-

Regulation of phosphatidylinositol 3-kinase signaling

Regulation of nuclear division

Count

Positive regulation of phosphatidylinositol 3-kinase signaling

Count

Regulation of mitotic nuclear division

8

Response to CAMP

12

5.0

Positive regulation of mitotic cell cycle

Substrate adhesion-dependent cell spreading

16

7.5

Regulation of chromosome segregation

Positive regulation of fibroblast proliferation

Chromosome separation

20

10.0

Sulfur compound catabolic process

Regulation of sister chromatid segregation

Icosanoid biosynthetic process

Metaphase/anaphase transition of mitotic cell cycle

Prostaglandin biosynthetic process

Regulation of metaphase/anaphase transition of cell cycle

Export across plasma membrane

Regulation of mitotic metaphase/anaphase transition

Keratan sulfate catabolic process

0.10

0.1

.15

0.20

0.25

0.30

0.35

0.04

0.06

0.08

0.10

GeneRatio

GeneRatio

C

Spindle

D

Condensed chromosome

Chromosomal region

Microtubule

Chromosome, centromeric region

p.adjust

Spindle pole

Condensed chromosome, centromeric region

0.006

Count

Midbody

0.004

11

Kinetochore

0.002

Condensed chromosome kinetochore

p.adjust

Mitotic spindle

Proteinaceous extracellular matrix

Spindle microtubule

Count

Condensed nuclear chromosome

2.5

0.004897994

5.0

Condensed nuclear chromosome, centromeric region

7.5

Cell division site

10.0

Cyclin-dependent protein kinase holoenzyme complex

12.5

Spindle midzone

Condensed chromosome outer kinetochore

Condensed nuclear chromosome kinetochore

Contractile ring

0.05

0.10

0.15

0.06

0.08

0.10

0.12

0.14

GeneRatio

GeneRatio

E

Tubulin binding

F

Growth factor activity

Microtubule binding

Drug binding

Count

Growth factor binding

Count

2

3

3

Cyclin-dependent protein kinase activity

4

Hormone binding

4

5

5

Cyclin-dependent protein serine/threonine kinase activity

6

6

Monocarboxylic acid transmembrane transporter activity

Histone kinase activity

p.adjust

p.adjust

0.04

DNA binding, bending

0.03

Insulin-like growth factor binding

0.03

0.02

0.02

Translation repressor activity, nucleic acid binding

0.01

Peptide hormone binding

0.01

Folic acid binding

Sequence-specific mRNA binding

Fibronectin binding

0.03 0.04 0.05 0.06 0.07 0.08

0.0300.0350.0400.0450.0500.055 GeneRatio

GeneRatio

FIGURE 4 Significantly enriched pathway terms of up-regulated differentially expressed genes in ACC

Cell cycle

Cellular senescence

Count

Progesterone-mediated oocyte maturation

2

☒ 3

Oocyte meiosis

☒ 4

☒ 5

☒ 6

p53 signalling pathway

p.adjust

Antifolate resistance

0.04

0.03

DNA replication

0.02

0.01

Folate biosynthesis

One carbon pool by folate

0.10

0.15

0.20

GeneRatio

FIGURE 5 Protein-protein interaction (PPI) network, module analysis and hub gene identification. Red nodes represent up-regulated genes. Green nodes represent down-regulated genes. A, PPI network of differentially expressed genes was constructed in STRING database. B, Top nine hub genes were selected by Cytoscape software based on the degree of each node

A

B

YBRD

MUM1

5orf3

RIP13

STOV

SMC4

TPX2

TOIS

DHAIC 2

MND1

GMNN

SC40A

ASPM

OSL

ARV

PY IR

NFIA

HF 1

JBE2

CDK1

BXAS

TYMS

AA010

FANC

THIR

PLINI

GADD45

YBL

ASP

JSPI

R2F

TLR4

TGA

ASEH

2A

12AF2

R

D51A

1

EXCL1

STATS

KIF2

MGB

CCNA2

PBK

A

NGPT

SIRT

SLC16A

CONB2

CDC6

VISP

IDC80

FEMP?

DGF

GF

VDM

Cr

PRC 1

BUB

NEXN

CONB

AURKA

LAPHE

1AD2L

CDKN3

IMP

GGH

BLN

DORA

R

R

IBE2

HOF

DOFF

KIF11

FBP

ENPK

ASPM

CDCA

CAPG

HMMR

ZWINT

TEK

OHFR

ICAPG

THRE

CONB

GFBP

SERPING

ANLN

BRE

PRC1

CTGE

RN

MPO

SMCA

REL

AAIZ

GFUP

DLG2

RFC4

BLN

D

POC

B

CEP55

AXL

RSPO

OMD

ANLN

DTL

KIF11

R

CGAP1

NF 36

OGN

MOD

TPRE

CTH

BCC

AGLI

AB39

JBE20

GLUIL

BIRC

A

DHIA

OP2A

S CO2E

1

URKA

IDC80

CADS

MYLK

ALAD

CENP

DH1

ADX1

S

HICBP1

(topoisomerase lla), CCNA2 (CyclinA2), CDKN3 (cyclin-depend- ent kinase inhibitor 3), MAD2L1 (mitosis arrest-deficient 2 like 1), RACGAP1 (Rac GTPase activating protein 1), BUB1 (benzimida- zole 1 homolog beta) and CCNB2 (Cyclin B2). In terms of biologi- cal process, these hub genes are significantly enriched in mitotic spindle assembly checkpoint (Table S1). In terms of molecular function, these hub genes are significantly enriched in ATP bind- ing (Table S1). KEGG pathway enrichment analysis showed that these hub genes are associated with progesterone-mediated

oocyte maturation, cell cycle, oocyte meiosis, p53 signal- ling pathway and progesterone-mediated oocyte maturation (Table S1).

3.5 Evaluate the prognostic value of hub genes

TCGA-ACC dataset was used to evaluate the prognostic value of nine hub genes by GEPIA. All patients with high hub gene expression

FIGURE 6 Prognostic value of 9 hub genes in ACCs by GEPIA. A, BUB1. B, CCNA2. C, CCNB1. D, CCNB2. E, CDK1. F, CDKN3. G, MAD2L1. H, RACGAP1. I, TOP2A P < . 05 was as statistically significant

A

Overall Survival

B

Overall Survival

C

Overall Survival

1.0

Low BUB1 TPM

1.0

High BUB1 TPM

Low CCNA2 TPM

1.0

Low CCNB1 TPM

Logrank P = 2.6e-06

High CCNA2 TPM

High CCNB1 TPM

HR(high) = 7.7

Logrank P = 3e-06

Logrank P = 1e-04

0.8

p(HR) = 5e-05

0.8

HR(high) = 7.5

p(HR) = 5 8e-05

0.8

HR(high) = 4.9

Per cent survival

n(high) = 38

Per cent survival

p(HR) = 0.00038

n(low) = 37

n(high) = 38

Per cent survival

n(high) = 38

0.6

0.6

n(low) = 38

0.6

n(low) = 38

0.4

0.4

0.4

0.2

0.2

0.2

0.0

0.0

0.0

0

50

100

150

0

50

100

150

0

50

100

150

Months

Months

Months

D

Overall Survival

E

Overall Survival

F

Overall Survival

1.0

Low CCNB2 TPM

1.0

High CCNB2 TPM

Low CDK1 TPM

1.0

Low CDKN3 TPM

Logrank P = 9.3e-06

High CDK1 TPM

High CDKN3 TPM

0.8

HR(high) = 6.9

Logrank P = 7e-08

(HR) = 0.0001

0.8

HR(high) = 11

Logrank P = 2.3e-06

p(HR) = 1.20-05

0.8

HR(high) = 7.6

Per cent survival

p(HR)= 4.9e-05

n(high) = 38

Per cent survival

n(high) = 38

Per cent survival

n(high) = 38

0.6

n(low) = 38

0.6

n(low) = 38

0.6

n(low) = 38

0.4

0.4

0.4

0.2

0.2

0.2

0.0

0.0

0.0

0

50

100

150

0

50

100

150

0

50

100

150

Months

Months

Months

G

Overall Survival

H

Overall Survival

Overall Survival

1.0

Low MAD2L 1 TPM

1.0

High MAD2L1 TPM

Low RACGAP1 TPM

1.0

High

Low TOP2A

Logrank P = . 00038

HR(high) = 4.4

Logrank P = 7.3e-06

High TOP2A TPM

0.8

p(HR) = 0.001

0.8

HR(high) = 7.1

Logrank P = . 00014

HR(high) = 4.7

Per cent survival

p (HR) = 9.8e-05

0.8

p(HR) = 0.00047

n(high) =38

n(high) = 38

0.6

n(low) =38

Per cent survival

n(high) = 38

Per cent survival

0.6

n(low) = 38

0.6

n(low) = 38

0.4

0.4

0.4

0.2

0.2

0.2

0.0

0.0

0.0

0

50

100

150

0

50

100

150

0

50

100

150

Months

Months

Months

were associated with worse OS (Figure 6). The additional univariate and multivariate Cox regression analysis showed that the hub gene BUB1 (budding uninhibited by benzimidazole 1) was an independent prognostic factor for ACC patients, and BUB1 was significantly as- sociated with living status and clinical stage in TCGA data (Table 1 and Table 2). The analysis results of the remaining genes were shown in Tables S2-S17. In addition, all nine hub genes were validated to be significantly up-regulated in ACCs as they were in above three GEO datasets (Figure 7).

|

4 DISCUSSION

In this study, by integrated three expression profiling datasets from GEO, we identified 200 commonly changed DEGs (74 up- regulated and 126 down-regulated) in ACCs. These DEGs were further analysed by GO analysis (molecular function, biological process and cellular component). In terms of biological process, the up-regulated DEGs were mainly enriched in organelle fission and nuclear division, which were typical malignant indicators of

TABLE 1 Clinicopathological parameters and BUB1 expression according to the TCGA database
ParametersBUB1 mRNA expression
GroupLow(n = 38)High(n = 39)×2P value
Age (Mean ± SD)46.63 ± 15.75646.59 ± 16.106
GenderFemale24240.0211.000
Male1415
Clinical stageI/II311514.877.000
III/IV724
Recurrence statusNo29112.363.188
Yes77
Null221
Living statusLiving341619.841.000
Dead423
TABLE 2 Univariate and multivariate Cox regression analysis of BUB1 clinical pathologic features according to the TCGA database
Parameters OSUnivariate analysisMultivariate analysis
HR95% CIPHR95% CIP
Age ≥60 vs < 601.5490.6773.548.3000.4940.2071.180.112
Gender Female vs Male0.9860.4512.154.971
Clinical stage I/II vs III/IV6.4672.70215.481.0000.2080.0770.558.002
BUB1 expression Low vs High9.0243.09426.320.0005.9071.92018.176.002

histopathological examination. Down-regulated DEGs were sig- nificantly enriched in angiogenesis. For example, Plk1 is closely related to a series of mitotic events such as centrosome replica- tion, spindle formation, chromosome segregation and cytokinesis and is also related to chromosome stability.21 In addition, down- regulated DEGs were significantly enriched in proteinaceous ex- tracellular matrix. In terms of molecular function, up-regulated DEGs were significantly enriched in tubulin binding and microtu- bule binding. It had been reported that STMN1 regulated micro- tubule dynamics and participates in the malignant phenotype of cancer cells.22 Down-regulated DEGs were significantly enriched in growth factor activity and growth factor binding. Some studies have shown that the stimulation of growth factors such as insulin- like growth factors may promote tumour proliferation.23 These results can help us to further understand the role of DEGs in the development and progress of ACC. The additional KEGG path- way analysis showed that up-regulated DEGs were significantly enriched in cell cycle, cell senescence, progesterone-mediated oocyte maturation, oocyte, p53 signalling pathway and folic acid resistance, further confirmed the important roles of p53 signalling pathway in ACC.24

By DEGs PPI network analysis, the hub genes with highest degree of communication were identified, and they are CCNB1,

CDK1, TOP2A, CCNA2, CDKN3, MAD2L1, RACGAP1, BUB1 and CCNB2. CCNB1, a key molecule to initiate mitosis, was associated with tumorigenesis and development.25 CCNB2 can promote cell proliferation, migration and invasion in lung adenocarcinoma 26 and hepatocellular carcinoma.27 In addition, CCNB1 and CCNB2 promoted gastric cancer cell proliferation and tumour growth.28 CDK1 and TOP2A played an important role in the regulation of cell cycle and regulated the proliferation of tumours.29,30 Glover et al found that CDK1 was up-regulated in ACC tissues com- pared with normal tissues.31 TOP2A was a potential biomarker for the progression and prognosis of various tumours.32 CCNA2 was found to promote the proliferation of breast cancer.33 Xing et al found that the expression of CDKN3 was generally increased in hepatocellular carcinoma tissues and was positively correlated with the pathological stage and differentiation of the tumours.34 MAD2L1 and BUB1 were important components of mitotic check- point complex proteins. High expression of these two genes was related to the poor disease-free survival of invasive tumours.35 RACGAP1 was found to be highly expressed in colorectal cancer36 and breast cancer. 37

Finally, we evaluated the prognostic value of 9 hub genes by using GEPIA database and found ACC patients with high expression of 9 hub genes (CCNB1, CDK1, TOP2A, CCNA2, CDKN3, MAD2L1,

WILEY

FIGURE 7 Nine hub genes are highly expressed in ACC tissues compared with normal tissues in GEPIA. The red and grey boxes represent cancer and normal tissues, respectively. A, BUB1. B, CCNA2. C, CCNB1. D, CCNB2. E, CDK1. F, CDKN3. G, MAD2L1. H, RACGAP1. I, TOP2A P < . 05 was as statistically significant

A

B

C

5

8

BUB1

*

1

CCNA2

*

CCNB1

*

4

6

6

5

0

+

4

2

0

2

2

-

-

0

0

0

ACC (num(T)=77; num(N)=128)

ACC (num(T)=77; num(N)=128)

ACC (num(T)=77; num(N)=128)

D

E

F

CCNB2

CDK1

CDKN3

6

6

6

5

5

5

4

+

4

0

3

0

~

2

~

-

-

1

0

0

0

ACC (num(T)=77; num(N)=128)

ACC (num(T)=77; num(N)=128)

ACC (num(T)=77; num(N)=128)

G

H

I

7

MAD2L1 -*

6

RACGAPI

TOP2A

*

6

5

6

5

+

+

4

3

0

2

2

2

1

1

:

0

0

0

ACC

ACC (num(T)=77; num(N)=128)

ACC (num(T)=77; num(N)=128)

(num(T)=77; num(N)=128)

RACGAP1, BUB1 and CCNB2) were significantly associated with worse OS.

In conclusion, using multiple cohort profiling datasets and inte- grated bioinformatics analysis, we identified hub genes and potential pathways that may be involved in the progress of ACC.

This study has been proved by the Henan University institutional committee.

ACKNOWLEDGEMENTS

This study was supported by Key Project of Science and Technology Research of Henan Provincial Department of Education (No.20A180001), National Natural Science Foundation of China (No.81602362), Supporting grants of Henan University (No.2015YBZR048; No.B2015151), Yellow River Scholar Program (No.H2016012) and Program for Innovative Talents of Science and Technology in Henan Province (No. 18HASTIT048), Projects for College Students in Henan University (No. 201819002), China Postdoctoral Science Foundation (No.2017M62237) and Henan Postdoctoral Foundation (No. 001702052), Program for Science and Technology Development in Henan Province (No.162102310391, No.172102210187), Program for Scientific and Technological Research of Henan Education Department (No.14B520022), Kaifeng Science and Technology Major Project (18ZD008) and Supporting grant of Bioinformatics Center of Henan University (No.2018YLJC01).

CONFLICT OF INTEREST

The authors declare that there is no conflict of interests. No animal or human studies were carried out by the authors for this article.

AUTHOR CONTRIBUTIONS

XG and LX involved in study concept and design; LX, JG, YG and XG involved in acquisition of data; JG, YG, XM, LZ, HL, ZY and YH in- volved in analysis and interpretation of data. LX, JG and XG drafted the manuscript; LX, JG and XG involved in critical revision of the manuscript for intellectual content.

DATA AVAILABILITY STATEMENT

All data generated or analysed during this study are included in this article.

ORCID

Jinshuai Guo İD https://orcid.org/0000-0002-2475-7847

Yinzhong Gu ID https://orcid.org/0000-0002-9335-0055

Xiaoyu Ma ID https://orcid.org/0000-0002-8279-6035

Lu Zhang ID https://orcid.org/0000-0001-8219-466X

Huimin Li İD https://orcid.org/0000-0001-8751-7501

Zhongyi Yan İD https://orcid.org/0000-0002-9835-5875

Yali Han İD https://orcid.org/0000-0002-5492-8117

Longxiang Xie İD https://orcid.org/0000-0001-6825-9690 Xiangqian Guo https://orcid.org/0000-0001-6645-5786 İD

REFERENCES

1. Delellis RA. Pathology and genetics of tumours of endocrine or- gans. 2004.

2. NIH state-of-the-science statement on management of the clin- ically inapparent adrenal mass (“incidentaloma”). NIH Consensus State Sci Statements. 2002;19:1.

3. Nakamura Y, Yamazaki Y, Felizola S, et al. Adrenocortical carcinoma: review of the pathologic features, production of adrenal steroids, and molecular pathogenesis. Endocrinol Metabo Clin North Am. 2015;44:399-410.

4. Bruno A, Martin F. Clinical review: adrenocortical carcinoma: clini- cal update. J Clin Endocrinol Metab. 2006;91:2027-2037.

5. Guillaume A, Eric L, Martin F, et al. Integrated genomic characteri- zation of adrenocortical carcinoma. Nat Genet. 2014;46:607-612.

6. Tobias E, Kim AC, Aaron S, et al. Adrenocortical carcinoma. Endocr Rev. 2014;35:282-326.

7. Gara SK, Lack J, Zhang L, Harris E, Cam M, Kebebew E. Metastatic adrenocortical carcinoma displays higher mutation rate and tumor heterogeneity than primary tumors. Nat Commun. 2018;9:4172.

8. Xie L, Wang Q, Nan F, et al. OSacc: Gene Expression-Based Survival Analysis Web Tool For Adrenocortical Carcinoma. Cancer Manag Res. 2019;11:9145.

9. Giordano TJ, Kuick R, Else T, et al. Molecular classification and prog- nostication of adrenocortical tumors by transcriptome profiling. Clin Cancer Res. 2009;15(2):668-676.

10. Soon PSH, Gill AJ, Benn DE, et al. Microarray gene expression and immunohistochemistry analyses of adrenocortical tumors identify IGF2 and Ki-67 as useful in differentiating carcinomas from adeno- mas. Endocr Relat Cancer. 2009;16(2):573-583.

11. Weiss ID, Huff LM, Evbuomwan MO, et al. Correction: Screening of cancer tissue arrays identifies CXCR4 on adrenocortical carcinoma: correlates with expression and quantification on metastases using (64)Cu-plerixafor PET. Oncotarget. 2018;8:73387-73406.

12. Ward A, Balwierz A, Zhang JD, et al. Re-expression of microRNA-375 reverses both tamoxifen resistance and accompanying EMT-like properties in breast cancer. Oncogene. 2013;32:1173-1182.

13. Consortium GO. The gene ontology (GO) project in 2006. Nucleic Acids Res. 2006;34:322-326.

14. Han YJ, Wang XH, Chen WC, et al. Differential expression of ca- rotenoid-related genes determines diversified carotenoid color- ation in flower petal of Osmanthus fragrans. Tree Genet Genomes. 2014;10:329-338.

15. Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;27:29-34.

16. Dennis G, Sherman BT, Hosack DA, et al. DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 2003;4:P3.

17. Damian S, Andrea F, Stefan W, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43:D447.

18. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environ- ment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498.

19. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analy- ses. Nucleic Acids Res. 2017;45:W98-W102.

20. Goldman M, Craft B, Zhu J, Haussler D. The UCSC Xena system for cancer genomics data visualization and interpretation. AACR. 2017;77:2584.

21. Klaus S, Axel U. Targeting polo-like kinase 1 for cancer therapy. Nat Rev Cancer. 2006;6:321-330.

22. JoãO Agostinho MN, Saad STO, Fabiola T. Stathmin 1 in normal and malignant hematopoiesis. BMB Rep. 2014;47:660-665.

23. Martino MCD, Koetsveld PMV, Pivonello R, Hofland LJ. Role of the mTOR Pathway in normal and tumoral adrenal cells. Neuroendocrinology. 2010;92:28-34.

24. Faillot S, Assie G. ENDOCRINE TUMOURS: The genomics of adre- nocortical tumors. Eur J Endocrinol. 2016;174(6):R249-R265.

25. Ding K, Li W, Zou Z, Zou X, Wang C. CCNB1 is a prognostic bio- marker for ER+ breast cancer. Med Hypotheses. 2014;83:359-364.

26. Stav D, Bar I, Sandbank J. Usefulness of CDK5RAP3, CCNB2, and RAGE genes for the diagnosis of lung adenocarcinoma. Int J Biol Markers. 2007;22:108-113.

27. Gao CL, Wang GW, Yang GQ, Yang H, Zhuang L. Karyopherin sub- unit-a 2 expression accelerates cell cycle progression by upregu- lating CCNB2 and CDK1 in hepatocellular carcinoma. Oncol Lett. 2017;15:2815-2820.

28. Shi Q, Wang W, Jia Z, Chen P, Ma K, Zhou C. ISL1, a novel reg- ulator of CCNB1, CCNB2 and c-MYC genes, promotes gastric cancer cell proliferation and tumor growth. Oncotarget. 2016;7: 36489-36500.

29. Xi Q, Huang M, Wang Y, et al. The expression of CDK1 is associated with proliferation and can be a prognostic factor in epithelial ovar- ian cancer. Tumour Biol. 2015;36:4939-4948.

30. Kang X, Song C, Du X. PTEN stabilizes TOP2A and regulates the DNA decatenation. Sci Rep. 2015;5:17873.

31. Glover AR, Ting ZJ, Gill AJ, et al. microRNA-7 as a tumor suppres- sor and novel therapeutic for adrenocortical carcinoma. Oncotarget. 2015;6:36675-36688.

32. Lee YE, He HL, Lee SW, et al. AMACR overexpression as a poor prognostic factor in patients with nasopharyngeal carcinoma. Tumour Biol. 2015;35:7983-7991.

33. Tian G, Yong H, Ling Y, Sheng A, Ziyu L, Jiafu J. CCNA2 is a prognos- tic biomarker for ER+ breast cancer and tamoxifen resistance. PLoS ONE. 2014;9:e91771.

34. Xing C, Xie H, Lin Z, et al. Cyclin-dependent kinase inhibitor 3 is overexpressed in hepatocellular carcinoma and promotes tumor cell proliferation. Biochem Biophys Res Commun. 2012;420:29-35.

35. Zhanwei W, Dionyssios K, Yi S, et al. Biological and clinical sig- nificance of MAD2L1 and BUB1, genes frequently appearing in expression signatures for breast cancer prognosis. PLoS ONE. 2015;10:e0136246.

36. Hiroki I, Yuji T, Susumu S, et al. RacGAP1 expression, increasing tumor malignant potential, as a predictive biomarker for lymph node metastasis and poor prognosis in colorectal cancer. Carcinogenesis. 2015;36:346-354.

37. Pliarchopoulou K, Kalogeras KT, Kronenwett R, et al. Prognostic significance of RACGAP1 mRNA expression in high-risk early breast cancer: a study in primary tumors of breast cancer patients partici- pating in a randomized Hellenic Cooperative Oncology Group trial. Cancer Chemother Pharmacol. 2013;71:245-255.

SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

How to cite this article: Guo J, Gu Y, Ma X, et al. Identification of hub genes and pathways in adrenocortical carcinoma by integrated bioinformatic analysis. J Cell Mol Med. 2020;24:4428-4438. https://doi.org/10.1111/jcmm.15102