PeerJ
Submitted 25 October 2018 Accepted 2 February 2019 Published 14 March 2019 Corresponding author Qing-Peng Kong, kongqp@mail.kiz.ac.cn
Academic editor Shi-Cong Tao Additional Information and Declarations can be found on page 10
DOI 10.7717/peerj.6555
cc Copyright 2019 Xia et al.
Distributed under Creative Commons CC-BY 4.0
OPEN ACCESS
Identification of four hub genes associated with adrenocortical carcinoma progression by WGCNA
Wang-Xiao Xia1,2,3,4, Qin Yu1,2,3,4, Gong-Hua Li1,2,3, Yao-Wen Liu1,2,3, Fu-Hui Xiao1,2,3, Li-Qin Yang1,2,3, Zia Ur Rahman1,2,3,4, Hao-Tian Wang1,2,3,4 and Qing-Peng Kong1,2,3
State Key Laboratory of Genetic Resources and Evolution/Key Laboratory of Healthy Aging Research of Yunnan Province, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming, China
2 Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming, China
3 Kunming Key Laboratory of Healthy Aging Study, Kunming, China
4 Kunming College of Life Science, University of Chinese Academy of Sciences, Beijing, China
ABSTRACT
Background: Adrenocortical carcinoma (ACC) is a rare and aggressive malignant cancer in the adrenal cortex with poor prognosis. Though previous research has attempted to elucidate the progression of ACC, its molecular mechanism remains poorly understood.
Methods: Gene transcripts per million (TPM) data were downloaded from the UCSC Xena database, which included ACC (The Cancer Genome Atlas, n = 77) and normal samples (Genotype Tissue Expression, n = 128). We used weighted gene co-expression network analysis to identify gene connections. Overall survival (OS) was determined using the univariate Cox model. A protein-protein interaction (PPI) network was constructed by the search tool for the retrieval of interacting genes.
Results: To determine the critical genes involved in ACC progression, we obtained 2,953 significantly differentially expressed genes and nine modules. Among them, the blue module demonstrated significant correlation with the “Stage” of ACC.
Enrichment analysis revealed that genes in the blue module were mainly enriched in cell division, cell cycle, and DNA replication. Combined with the PPI and co-expression networks, we identified four hub genes (i.e., TOP2A, TTK, CHEK1, and CENPA) that were highly expressed in ACC and negatively correlated with OS. Thus, these identified genes may play important roles in the progression of ACC and serve as potential biomarkers for future diagnosis.
Subjects Bioinformatics, Genetics, Oncology Keywords ACC, WGCNA, Hub genes, Progression
INTRODUCTION
Adrenocortical carcinoma (ACC) is a rare and aggressive malignant cancer found in the adrenal cortex (Fay et al., 2014). While this disease can occur at any age, it tends to show a bi-modal distribution with an initial peak in childhood (1-6 years old) and a second
How to cite this article Xia W-X, Yu Q, Li G-H, Liu Y-W, Xiao F-H, Yang L-Q, Rahman ZU, Wang H-T, Kong Q-P. 2019. Identification of four hub genes associated with adrenocortical carcinoma progression by WGCNA. PeerJ 7:e6555 DOI 10.7717/peerj.6555
peak in middle-age (40-50 years old) (Kiseljak-Vassiliades et al., 2018). As ACC has no obvious phenotypic traits at the early stage, almost 70% of patients are at stage III or IV when diagnosed (Bharwani et al., 2011; Fay et al., 2014). At these stages, ACC is invasive and metastatic, with patients at stage IV only having a 5-year survival of 6-13% (Else et al., 2014; Fassnacht et al., 2009; Fassnacht, Kroiss & Allolio, 2013). Unfortunately, current ACC therapies, such as surgery, chemotherapy, and radiotherapy, exhibit poor performance and outcomes (Allolio & Fassnacht, 2006). While next generation sequencing technology recently identified several genetic molecules associated with ACC (Soon et al., 2008; Assié et al., 2014; Greenhill, 2016; Zheng et al., 2016; Chortis et al., 2018), our understanding of ACC progression at each stage remains incomplete and treatment options are limited (Hoang, Ayala & Albores-Saavedra, 2002; Cherradi, 2016).
Thus, integrated analysis is required to further understand the molecular characterization of ACC gene expression, which may indicate stage and identify additional biomarkers for further research and clinical therapies.
Traditional methods of identifying the functional genes of ACC have focused on screening differentially expressed genes (DEGs) (Giordano et al., 2003; Slater et al., 2006; Lombardi et al., 2006), with limited attention paid to gene connections. Weighted gene co-expression network analysis (WGCNA) is a popular method in systems biology that can construct gene networks and detect gene modules (Clarke et al., 2013; Yang et al., 2014; Lee et al., 2015; Goldman et al., 2017; Sun et al., 2017). By analyzing the connectivity between modules and clinical traits, we can determine which modules are associated with which traits. Those genes found in the center of a regulation network usually exhibit more important functions. Thus, the degree of gene connectivity in one module can also be analyzed by the gene-gene interaction/regulation network, from which critical hub genes can be identified.
In this study, we identified genes involved in ACC progression via comprehensive transcriptome-wide analysis of ACC gene expression patterns. We systematically analyzed clusters of genes with similar expression patterns using WGCNA and found the MEblue module to be highly related to clinical stage. Further analysis identified four hub genes (i.e., TOP2A, TTK, CHEK1, and CENPA) from the module that were associated with ACC progression and prognosis. Thus, these hub genes may serve as candidate biomarkers of ACC in clinical treatment and contribute to a greater understanding of ACC progression.
MATERIALS AND METHODS
Data collection
We obtained gene expression transcripts per million (TPM) values (Table S1) from the UCSC Xena (http://xena.ucsc.edu/) database, which included 77 ACC samples from the cancer genome atlas (TCGA) (https://cancergenome.nih.gov/) and 128 normal samples from genotype tissue expression (GTEx) (https://www.gtexportal.org/home/). The two databases raw sequencing reads were recalculated with a unifying pipeline. Clinical data were downloaded from TCGA using the “cgdsr” package in R (v3.1.3) (R Core Team, 2015; Jacobsen, 2015).
DEG screening
Of the 60,498 genes in each sample, we removed genes with a mean TPM ≤ 2.5 (>1 is a common cutoff for determining if an isoform is expressed or not (Liu, Jing & Tu, 2016)) in the cancer and normal samples and thus retained 13,987 genes. For those genes in the samples that showed significant changes, we used analysis of variance (ANOVA) in R (R Core Team, 2013) to determine the variance in genes between the two groups. ANOVA is a collection of statistical models useful for DEG analysis (Alabi et al., 2018; Monterisi et al., 2015). We obtained 2,953 significant DEGs (Table S2) in ACC with a p < 0.001 and |log2 (fold-change)| > 1 cutoff.
Co-expression network construction by WGCNA
Weighted gene co-expression network analysis (v1.49) can be applied to identify global gene expression profiles as well as co-expressed genes. Therefore, we installed WGCNA package for co-expression analysis using Bioconductor (http://bioconductor.org/biocLite.R). We used the soft threshold method for Pearson correlation analysis of the expression profiles to determine the connection strengths between two transcripts to construct a weighted network. Average linkage hierarchical clustering was carried out to group transcripts based on topological overlap dissimilarity in network connection strengths. To obtain the correct module number and clarify gene interaction, we set the restricted minimum gene number to 30 for each module and used a threshold of 0.25 to merge the similar modules (see the detailed R script in File S1).
Identification of clinically significant modules
We used two methods to identify modules related to clinical progression traits. Module eigengenes (MEs) are the major component for principal component analysis of genes in a module with the same expression profile. Thus, we analyzed the relationship between MEs and clinical traits and identified the relevant modules. We used log10 to transform the p-value from the linear regression between gene expression and clinical stage, which was defined as gene significance. Average gene significance in a module was defined as module significance.
Functional and pathway enrichment analysis
The database for annotation visualization and integrated discovery (DAVID) (v6.8) (http://david.abcc.ncifcrf.gov/) was used for functional annotation of genes to better understand their biological functions. All genes in the blue module were uploaded for gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses, with cutoffs of p < 0.01 and p < 0.05 established for significant biological processes and pathways, respectively.
PPI and co-expression analysis
Genes were uploaded to the search tool for the retrieval of interacting genes (STRING) (v10.5) (https://string-db.org/) database. Confidence was set to more than 0.4 and other
parameters were set to default. We visualized the gene co-expression network with Cytoscape (v2.7.0) (Shannon et al., 2003).
Gene expression correlation with stage and survival analysis
The correlation between gene expression and stage was determined using GEPIA (http://gepia.cancer-pku.cn/index.html) (Tang et al., 2017). The correlation between gene expression and overall survival (OS) was established using the Cox model. A hazard ratio p-value of <0.01 was considered significant. Each gene with higher expression in the ACC samples had corresponding lower survival expectation. The “limma” (Ritchie et al., 2015) R package was used to test for the significantly expressed gene in GSE10927.
RESULTS
Construction and analysis of gene co-expression network with DEGs in ACC
Genes with a mean TPM ≤ 2.5 were removed from the two groups and the remaining 13,987 genes were used for differential expression analysis with ANOVA. In total, 2,953 significant DEGs were identified with a cutoff of p < 0.001 and |log2(fold-change)| > 1 (Fig. 1A), which included 1,181 up-regulated and 1,772 down-regulated genes (Fig. 1B). The 2,953 gene expression levels in ACC and normal samples are shown in the heatmap in Fig. 1C and Table S2.
Genes with similar expression patterns may participate in similar biological processes or networks (Mao et al., 2009). To better understand the gene expression network during ACC development, the co-expression network of the 2,953 DEGs was analyzed by WGCNA. First, to determine whether all 77 ACC samples were suitable for network analysis, the sample dendrogram and corresponding clinical traits were analyzed. We found that all samples were included in the clusters and passed the cutoff thresholds (Fig. 1D). The power value is a critical parameter that can affect the independence and average connectivity degree of the co-expression modules. Thus, network topology using different soft thresholding powers was screened, with $ = 6 (scale free R2 = 0.928) selected for later analysis (Figs. 1E, 1F and 1G). We then constructed the gene co- expression network using WGCNA based on the hierarchical clustering of the calculated dissimilarities, and nine modules were obtained (Fig. 1H; Table S3). We used eigengenes as representative profiles and quantified module similarity by eigengene correlation (Fig. 1I).
Correlation of blue module with clinical stage and progression
We investigated whether any module was correlated with clinical stage and tested the relevance between each module and ACC clinical traits. We found that module significance of the blue module was higher than that of any other, implying it had greater correlation with ACC stage (Fig. 2A). The blue module also displayed a positive correlation with ACC clinical stage (r = 0.5, p = 6e-06) and negative correlation with OS (r = - 0.56, p = 3e-07) (Fig. 2B). The eigengene dendrogram and heatmap indicated that the MEblue and MEyellow modules were highly correlated with clinical stage
A
B
C
Cancer
Normal
1,500
5
-Log10(P value)
· Down-regulated genes
. No change genes
· Up-regulated genes
Gene number
-5
10
1,000
500
0
0
-6
-3
0
3
6
Up-regulated gene
Down-regulated gene
Log2 (Fold change)
D
180
Height
TOGA.OR A5L9.01
TCGA.OR.ASK0.01
W
Scale Free Topology Model Fit, signed R
Scale independence
F
1.0
Mean connectivity
140
6 7 8 910
12
14 16
18
20
600
1
0.8
4
5
500
3
100
TEGA.PELASOG.61
AbLO.01
Mean Connectivity
0.6
400
TEGAOF
TEGA
300
TOGA
TESA
0.4
2
200
2
Age
0.2
Gender
100
Os status
3
Stage
0.0
4
T
1
0
6
7 8 910 12 14
16
18
20
M
N
5
10
15
20
5
10
15
20
Soft Threshold (power)
Soft Threshold (power)
G
H
scale R^2= 0.92 , slope =- 1.41
Network Heatmap Plot
1.0
Cluster Dendrogram
-0.5
0.8
Log10 (p(k))
-1.0
Height
0.6
-1.5
0.4
0.2
-2.0
Dynamic Tree Cut
0
Merged dynamic
0.6
0.8
1.0
1.2
1.4
1.6
1.8
Log10 (k)
A
C
E
Biological Process Molecular Function Cellular Component
Gene significance across modules, p-value = 2.9e-148
60
Gene significance
0.00 0.05 0.10 0.15 0.20 0.25
0.2 0.4 0.6 0.8 1.0
MEpink
40
MEblack
Stage
MEgreen
MEred
MEbrown
MEturquoise
MEblue
MEyellow
- Log2(FDR)
20
0
1
NÃOTOO
Cell division
Mitotic nuclear division Sister chromatid cohesion
G1/S transition of mitotic cell cycle
DNA replication
Protein binding
ATP binding
Chromosome, centromeric region Condensed chromosome kinetochore
midbody
Kinetochore
Spindle
0 8
00 )
Stage
Pink
Black
Green
Red
Brown
Turquoise
Blue
Yellow
Grey
Stage
B
D
TI
Module - traitrelationships
60
MEpink 76
-0.067
0.081
-0.14
0.074
0.13
(0.6)
(0.5)
(0.2)
(0.5)
(0.3)
0.0096 (0.9)
0.047 (0.7)
1
Module membership vs. gene significance cor = 0.76, p = 2.2e-123
0.5
MEblack
84
0.021
-0.23
-0.28
0.22
0.17
0.3
0.082
(0.9)
(0.05)
(0.01)
(0.05)
(0.1)
(0.009)
(0.5)
0
-0.5
Gene significance for stage
0.5
MEgreen 279
0.082
0.049
0.07
-0.088
-0.1
0.062
-0.11
-Log2 (P value)
(0.5)
(0.7)
(0.6)
(0.5)
(0.4)
(0.6)
(0.4)
40
-1
MEred
102
-0.12
-0.034
(0.3)
-0.11
0.0019 (1)
-0.017
0.11
-0.069
0.4
(0.4)
(0.8)
(0.9)
(0.3)
(0.6)
MEbrown 282
0.18
(0.1)
-0.052 (0.7)
-0.014 (0.9)
0.24
0.14
0.061
(0.04)
0.19 (0.1)
(0.2)
(0.6)
0.3
MEturquoise 994
-0.011
-0.25
-0.23
0.24
0.16
0.21
0.12
(0.9)
20
(0.03)
(0.05)
(0.04)
(0.2)
(0.07)
(0.3)
0.2
MEblue
650
-0.088
-0.037
-0.56
0.5
0.47
0.41
0.25
(0.5)
(0.8)
(3e-07)
(6e-06)
(3e-05)
(3e-04)
(0.03)
MEyellow
281
-0.11 (0.4)
-0.31
-0.31
0.25
0.1
0.2
0.28
(0.007)
(0.008)
(0.03)
0.16 (0.2)
(0.1)
(0.01)
MEgrey
205
0.00018
0.28
0.38
-0.11
-0.054 (0.7)
-0.2
-0.064
(1)
(0.02)
(0.001)
(0.4)
(0.09)
(0.6)
0.0
0
Age
Gender
Os status
Stage
1
M
N
0.2
0.4
0.6
0.8
Module membership in blue module
Cell cycle
DNA replication
Oocyte meiosis
Progesterone-mediated
oocyte maturation
HTLV-I infection
p53 signaling pathway
Fanconi anemia pathway
Base excision repair
Mismatch repair
To determine the function of the 650 genes in the blue module, GO and KEGG function and pathway enrichment analyses were performed by DAVID functional annotation (Huang, Sherman & Lempicki, 2009). For GO biological processes, genes in the module were significantly enriched in cell division (p = 1.05e-26) (Fig. 2E; Table S4), whereas for KEGG pathway analysis, the genes were mainly enriched in cell cycle (p = 2.7e-19) and DNA replication (p = 8.27e-8) (Fig. 2F; Table S5) pathways. These processes and pathways all play critical roles in cancer progression (Tachibana, Gonzalez & Coleman, 2005), implying that genes in this module may be involved in ACC progression.
PPI and co-expression networks to identify hub genes in ACC progression
To clarify high confidence hub genes, we entered the blue module genes into the STRING (Szklarczyk et al., 2015) database. The genes were ranked by the protein-protein interaction nodes and the top 5% of genes (16 genes) were chosen as candidate hub genes (Fig. 3A; Table S6). As highly connected hub genes in a module play important roles in biological processes (Liu, Jing & Tu, 2016), genes in the blue module were ranked by their degree of gene co-expression connectivity (Table S7). To identify genes that may play notable roles in ACC progression, the top 5% of genes (31 genes) (Fig. 3B) in the blue module with the highest connectivity were classified as candidate hub genes for further analysis. Finally, four common genes (i.e., TOP2A, TTK, CHEK1, and CENPA) in the two analysis were identified as hub genes in ACC (Fig. 3C). These four genes were highly expressed in ACC samples compared with normal samples (Figs. 3D-3G), indicating that they likely act as oncogenes in ACC. Further analysis of the GSE10927 dataset, which included microarray data of 10 normal samples and 33 ACC samples (Human Genome U133A 2.0 Plus; Affymetrix, Santa Clara, CA, USA) (Giordano et al., 2009), demonstrated that the four genes showed significant high expression in ACC (Figs. SIA-S1D). Furthermore, based on immunoreactivity experiments, TOP2A is reported to be highly expressed in ACC (Giordano et al., 2003).
Significant associations of hub genes with ACC stage and survival
We investigated the four hub genes to better understand their functions. We found that TOP2A, TTK, CHEK1, and CENPA play critical roles in biological processes that are highly correlated with cancer (Dominguez-Brauer et al., 2015), such as DNA topological structure, cell cycle progression, and mitosis (Hoffmann et al., 2011; Liu et al., 2000; De Resende et al., 2013; Thu et al., 2018), thereby suggesting their possible role in cancer development. Further exploration of their expression patterns during ACC clinical progression showed that the levels of these genes were significantly altered with clinical stage and markedly increased at stage III and IV (Figs. 4A-4D). This correlation between the expression levels of the four genes and ACC progression may be useful in ACC diagnosis.
Tumor prognosis is an important feature in cancer and has attracted considerable attention. To assess the utility of WGCNA at identifying hub genes indicative of ACC, we conducted survival analysis (Figs. 4E-4H). We separated the samples into two groups
A
POLD1
B
C
FAM
CEN
83D
PAK4
_1G
PPIA
RAD5
PK
AP1
RAN
PRIM1
PCNA
TK1
CDC
CDT1
NCA
20
PH
FEN1
TPX2
NDC
RAD51
TYMS
SMC4
80
UBE2T
NCAPG2
FANCI
Hub genes in co-expression network 27 genes
FEN1
EXO1
DTL
KIF
CDC6
TTK
RNAS
EH2A
MCME
GPS M2
SPC
4A
24
CDC 45
CHEK1
4 genes
RACGAPI
ZWINT
TOP2A
CEP55
CENPA
NCA
POL
PG2
DHFR
RRM2
NEK2
E2
NCA PG
CEN
TTK
TPX2
PH
Top2a Chek1 Ttk Cenpa
BUB1
RAD54L
HMMR
RRM2
476
CDK1
ASPM
GAS
PRC1
CHE
TOR
CDCAB
PLK1
RAD
MK1
MAD211
K1
2A
ANLN
54L
67
MXD3
2L3
12 genes
PTTG1
AURKB
HJURP
Hub genes in PPI network
KIF20A
CDC20
SPAG5
CDKN3
ASPM
STM
TAC
DLG
N1
C3
AP5
PLK1
HMG
TOP
CENPN
B2
2A
MK167
UBEZC KIF 2C
PBK
FOXMI
CDC25C
KIF4A
CCNF
REE
GTS
H2A
CEN
RACG
NEK2
CKS2
P4
E1
FZ
PA
AP1
RMI2
GTSE1
PBK
CONF
D
E
F
G
Top2a
Ttk
Chek1
5
Cenpa
5
4
+
Log2 (TPM+1)
4
3
in theint ..
Log2 (TPM+1)
Log2 (TPM+1)
Log2 (TPM+1)
3
+
3
V
N
2
2
A
-
L
0
0
0
0
Cancer
Normal
Cancer
Normal
Cancer
Normal-
Cancer
Normal
according to median gene expression levels and performed survival analysis using the Cox model. Survival analysis showed that the expression of all four genes was significantly correlated with OS (Figs. 4E-4H), with higher expression associated with lower patient survival time. The correlation between the hub genes and ACC prognosis suggests that these four genes likely contribute to the progression of ACC.
DISCUSSION
As ACC exhibits no obvious phenotypic traits during its early stages, diagnosis is often delayed in many patients (Bharwani et al., 2011; Fay et al., 2014). We systematically analyzed gene expression and found potential biomarker genes for ACC diagnosis. To identify genes that may play central roles in ACC progression, gene co-expression
A
Top2a
B
Ttk
C
Chek1
D
Cenpa
5
6
TPM (Transcripts Per Million)
F value = 7.28
Pr(>F) = 0.000251
TPM (Transcripts Per Million)
F value = 2.73
Pr(>F) = 0.0501
TPM (Transcripts Per Million)
F value = 5.17
Pr(>F) = 0.00274
5
F value = 8.49
Pr(>F) = 6.8e-05
1
TPM (Transcripts Per Million)
+
5
0
4
5
+
0
0
+
0
0
2
2
2
~
1
-
-
-
0
Stage I
Stage II
Stage III
Stage IV
0
Stage I
Stage II
Stage III
Stage IV
Stage I
Stage II
Stage III
Stage IV
Stage I
Stage II
Stage III
Stage IV
E
Top2a
F
Ttk
G
Chek1
H
Cenpa
1.0
Low (n=38)
1.0
Low (n=38)
1.0
Low (n=38)
1.0
Low (n=38)
High (n=38)
High (n=38)
High (n=38)
High (n=38)
0.8
0.8
0.8
0.8
Percent survival
0.6
Percent survival
0.6
Percent survival
0.6
Percent survival
0.6
0.4
0.4
0.4
0.4
0.2
Logrank p=0.00014
0.2
Logrank p=0.00029
0.2
Logrank p=0.0011
0.2
Logrank p=8.4e-07
HR(high)=4.7
p(HR)=0.00047
HR(high)=4.7
HR(high)=3.8
HR(high)=9.4
0.0
0.0
p(HR)=0.00095
0.0
p(HR)=0.0024
0.0
p(HR)=4.2e-05
0
50
100
150
0
50
100
150
0
50
100
150
0
50
100
150
Months
Months
Months
Months
Figure 4 Significant correlation between hub gene expression with pathological stage and survival. (A-D) Significant correlation between expression levels of TOP2A, TTK, CHEK1, and CENPA with ACC pathological stage. (E-H) Survival plot of OS in ACC. Higher expression (red line) of TOP2A, TTK, CHEK1, and CENPA indicates poorer prognosis. HR, hazard ratio. Full-size DOI: 10.7717/peerj.6555/fig-4
network analysis was conducted using WGCNA, which can describe correlation patterns among genes at the RNA level. Based on WGCNA, we obtained nine modules, with each module containing an average of 217 genes. Only 205 genes were unclassified in any module (in gray), accounting for 10.50% of DEGs. In comparison, previous studies have reported an average gene number in each module of 216 to 336 and percentage of genes not found in any module of 5.67-33.61% of DEGs (Liu et al., 2017, 2018; Yang et al., 2018; Zuo et al., 2018). In conclusion, our WGCNA results were comparable. We identified four hub genes (i.e., TOP2A, TTK, CHEK1, and CENPA) in the network center related to gene regulation and possible carcinogenesis.
Genes located in the central position of a gene-gene interaction network likely exhibit more important functions than other genes. Further investigation found that these four hub genes contribute to several tumor types indeed. For instance, TOP2A (topoisomerase II alpha), a specific marker of cell proliferation, is the primary molecular target of anthracyclines used for treating breast cancer (Villman et al., 2006; Wang et al., 2012). TTK, also known as monopolar spindle 1 (MPS1), plays a key role in cancer cell growth and proliferation, with its inhibition able to decrease tumor aggressiveness (Al-Ejeh et al., 2014; Maire et al., 2015; Zhu et al., 2018). CHEK1 (checkpoint kinase 1), a conserved serine/threonine kinase, plays a key role in tumor growth promotion (Zhang & Hunter, 2013). Furthermore, inhibition of CHEK1 expression by UCN-01, CEP-3891 (Zhu et al., 2018), AZD7762, or LY2606368 inhibitors (Manic et al., 2017) can prevent the proliferation of cancer cells (Bryant, Rawlinson & Massey, 2014; Schuler et al., 2017).
CENPA (centromere protein A), a histone H3 variant, is highly expressed in cancers, including breast, colorectal, liver, lung, ovarian, and osteosarcoma (Athwal et al., 2015; Sun et al., 2016; Filipescu et al., 2017). In addition, inhibition of CENPA expression in cancer cells can reduce sphere forming ability, proliferation, and cell viability (Behnan et al., 2016). Here, our study revealed that the expression levels of all four hub genes were significantly correlated with ACC progression (Figs. 4A-4D) and OS (Figs. 4E-4H), suggesting their critical function in ACC. Our results indicated that these four genes may play key roles in ACC tumorigenesis. However, the specific functions of these genes that contribute to ACC cell proliferation, differentiation, and metastasis need further study.
CONCLUSIONS
Based on gene co-expression network analysis, we identified four hub genes that likely contribute to the progression of ACC. The expressions of the four hub genes demonstrated significant correlation with ACC clinical stage and prognosis (Figs. 4A-4H). Thus, these four genes may act as potential biomarkers in predicting clinical outcomes and diagnosis of ACC. Furthermore, inhibitors of TOP2A, TTK, and CHEK1, which are already used for treating certain cancers, could potentially be used in ACC treatment. Further experimental and clinical studies are required to extend these findings.
ACKNOWLEDGEMENTS
We thank Qiong-Hua Gao for suggestions in modifying the paper and Christine Watts for help in honing the manuscript.
ADDITIONAL INFORMATION AND DECLARATIONS
Funding
This study was supported by grants from the National Natural Science Foundation of China (81602346), Applied Basic Research Projects of Yunnan Province (2018FB137), Chinese Academy of Sciences (QYZDB-SSW-SMC020 and KJZD-EW-L14). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Grant Disclosures
The following grant information was disclosed by the authors:
National Natural Science Foundation of China: 81602346. Applied Basic Research Projects of Yunnan Province: 2018FB137. Chinese Academy of Sciences: QYZDB-SSW-SMC020 and KJZD-EW-L14.
Competing Interests
The authors declare that they have no competing interests.
Author Contributions
· Wang-Xiao Xia conceived and designed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper.
· Qin Yu authored or reviewed drafts of the paper.
. Gong-Hua Li analyzed the data, authored or reviewed drafts of the paper.
· Yao-Wen Liu authored or reviewed drafts of the paper.
. Fu-Hui Xiao authored or reviewed drafts of the paper.
· Li-Qin Yang contributed reagents/materials/analysis tools.
· Zia Ur Rahman authored or reviewed drafts of the paper.
· Hao-Tian Wang prepared figures and/or tables.
· Qing-Peng Kong conceived and designed the experiments, authored or reviewed drafts of the paper, approved the final draft.
Data Availability
The following information was supplied regarding data availability: The raw measurements are available in the Supplemental Files.
Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.6555#supplemental-information.
REFERENCES
Alabi N, Sheka D, Gupta M, Kannappan S. 2018. Identification of a pathway-based 5-gene expression signature for predicting outcomes in gastric cancer. Journal of Proteomics & Bioinformatics 11(9):161-168 DOI 10.4172/jpb.1000482.
Al-Ejeh F, Simpson PT, Saunus JM, Klein K, Kalimutho M, Shi W, Miranda M, Kutasovic J, Raghavendra A, Madore J, Reid L, Krause L, Chenevix-Trench G, Lakhani SR, Khanna KK. 2014. Meta-analysis of the global gene expression profile of triple-negative breast cancer identifies genes for the prognostication and treatment of aggressive breast cancer. Oncogenesis 3(10):e124 DOI 10.1038/oncsis.2014.41.
Allolio B, Fassnacht M. 2006. Adrenocortical carcinoma: clinical update. Journal of Clinical Endocrinology & Metabolism 91(6):2027-2037 DOI 10.1210/jc.2005-2639.
Assié G, Letouzé E, Fassnacht M, Jouinot A, Luscap W, Barreau O, Omeiri H, Rodriguez S, Perlemoine K, René-Corail F, Elarouci N, Sbiera S, Kroiss M, Allolio B, Waldmann J, Quinkler M, Mannelli M, Mantero F, Papathomas T, De Krijger R, Tabarin A, Kerlan V, Baudin E, Tissier F, Dousset B, Groussin L, Amar L, Clauser E, Bertagna X, Ragazzon B, Beuschlein F, Libé R, De Reyniès A, Bertherat J. 2014. Integrated genomic characterization of adrenocortical carcinoma. Nature Genetics 46(6):607-612.
Athwal RK, Walkiewicz MP, Baek S, Fu S, Bui M, Camps J, Ried T, Sung MH, Dalal Y. 2015. CENP-A nucleosomes localize to transcription factor hotspots and subtelomeric sites in human cancer cells. Epigenetics & Chromatin 8(1):2 DOI 10.1186/1756-8935-8-2.
Behnan J, Grieg Z, Joel M, Ramsness I, Stangeland B. 2016. Gene knockdown of CENPA reduces sphere forming ability and stemness of glioblastoma initiating cells. Neuroepigenetics 7:6-18 DOI 10.1016/j.nepig.2016.08.002.
Bharwani N, Rockall AG, Sahdev A, Gueorguiev M, Drake W, Grossman AB, Reznek RH. 2011. Adrenocortical carcinoma: the range of appearances on CT and MRI. American Journal of Roentgenology 196(6):W706-W714 DOI 10.2214/ajr.10.5540.
Bryant C, Rawlinson R, Massey AJ. 2014. Chk1 inhibition as a novel therapeutic strategy for treating triple-negative breast and ovarian cancers. BMC Cancer 14(1):570 DOI 10.1186/1471-2407-14-570.
Cherradi N. 2016. microRNAs as potential biomarkers in adrenocortical cancer: progress and challenges. Frontiers in Endocrinology 6:195 DOI 10.3389/fendo.2015.00195.
Chortis V, Taylor AE, Doig CL, Walsh MD, Meimaridou E, Jenkinson C, Rodriguez-Blanco G, Ronchi CL, Jafri A, Metherell LA, Hebenstreit D, Dunn WB, Arlt W, Foster PA. 2018. Nicotinamide nucleotide transhydrogenase as a novel treatment target in adrenocortical carcinoma. Endocrinology 159(8):2836-2849 DOI 10.1210/en.2018-00014.
Clarke C, Madden SF, Doolan P, Aherne ST, Joyce H, O’Driscoll L, Gallagher WM, Hennessy BT, Moriarty M, Crown J, Kennedy S, Clynes M. 2013. Correlating transcriptional networks to breast cancer survival: a large-scale coexpression analysis. Carcinogenesis 34(10):2300-2308 DOI 10.1093/carcin/bgt208.
De Resende MF, Vieira S, Chinen LT, Chiappelli F, Da Fonseca FP, Guimarães GC, Soares FA, Neves I, Pagotty S, Pellionisz PA, Barkhordarian A, Brant X, Rocha RM. 2013. Prognostication of prostate cancer based on TOP2A protein and gene assessment: TOP2A in prostate cancer. Journal of Translational Medicine 11(1):36 DOI 10.1186/1479-5876-11-36.
Dominguez-Brauer C, Thu KL, Mason JM, Blaser H, Bray MR, Mak TW. 2015. Targeting mitosis in cancer: emerging strategies. Molecular Cell 60(4):524-536 DOI 10.1016/j.molcel.2015.11.006.
Else T, Williams AR, Sabolch A, Jolly S, Miller BS, Hammer GD. 2014. Adjuvant therapies and patient and tumor characteristics associated with survival of adult patients with adrenocortical carcinoma. Journal of Clinical Endocrinology & Metabolism 99(2):455-461 DOI 10.1210/jc.2013-2856.
Fassnacht M, Johanssen S, Quinkler M, Bucsky P, Willenberg HS, Beuschlein F, Terzolo M, Mueller HH, Hahner S, Allolio B. 2009. Limited prognostic value of the 2004 International union against cancer staging classification for adrenocortical carcinoma: proposal for a revised TNM classification. Cancer 115(2):243-250 DOI 10.1002/cncr.24030.
Fassnacht M, Kroiss M, Allolio B. 2013. Update in adrenocortical carcinoma. Journal of Clinical Endocrinology and Metabolism 98:4551-4564.
Fay AP, Elfiky A, Telo GH, Mckay RR, Kaymakcalan M, Nguyen PL, Vaidya A, Ruan DT, Bellmunt J, Choueiri TK. 2014. Adrenocortical carcinoma: the management of metastatic disease. Critical Reviews in Oncology/Hematology 92(2):123-132.
Filipescu D, Naughtin M, Podsypanina K, Lejour V, Wilson L, Gurard-Levin ZA, Orsi GA, Simeonova I, Toufektchan E, Attardi LD, Toledo F, Almouzni G. 2017. Essential role for centromeric factors following p53 loss and oncogenic transformation. Genes & Development 31(5):463-480 DOI 10.1101/gad.290924.116.
Giordano TJ, Kuick R, Else T, Gauger PG, Vinco M, Bauersfeld J, Sanders D, Thomas DG, Doherty G, Hammer G. 2009. Molecular classification and prognostication of adrenocortical tumors by transcriptome profiling. Clinical Cancer Research 15(2):668-676 DOI 10.1158/1078-0432.ccr-08-1067.
Giordano TJ, Thomas DG, Kuick R, Lizyness M, Misek DE, Smith AL, Sanders D, Aljundi RT, Gauger PG, Thompson NW, Taylor JM, Hanash SM. 2003. Distinct transcriptional profiles of adrenocortical tumors uncovered by DNA microarray analysis. American Journal of Pathology 162(2):521-531 DOI 10.1016/s0002-9440(10)63846-1.
Goldman M, Craft B, Brooks AN, Zhu J, Haussler D. 2017. Abstract 2584: the UCSC Xena system for cancer genomics data visualization and interpretation. Cancer Research 77(13 Supplement):2584 DOI 10.1158/1538-7445.am2017-2584.
Greenhill C. 2016. Adrenal gland: the genetics of adrenocortical carcinoma revealed. Nature Reviews Endocrinology 12(8):433.
Hoang MP, Ayala AG, Albores-Saavedra J. 2002. Oncocytic adrenocortical carcinoma: a morphologic, immunohistochemical and ultrastructural study of four cases. Modern Pathology 15(9):973-978 DOI 10.1038/modpathol.3880638.
Hoffmann S, Dumont M, Barra V, Ly P, Nechemia-Arbely Y, McMahon MA, Hervé S, Cleveland DW, Fachinetti D. 2011. CENP-A is dispensable for mitotic centromere function after initial centromere/kinetochore assembly. Cell Reports 17(9):2394-2404 DOI 10.1016/j.celrep.2016.10.084.
Huang DW, Sherman BT, Lempicki RA. 2009. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Research 37(1):1-13 DOI 10.1093/nar/gkn923.
Jacobsen A. 2015. cgdsr: R-based API for accessing the MSKCC cancer genomics data server (CGDS). Available at https://github.com/cBioPortal/cgdsr.
Kiseljak-Vassiliades K, Zhang Y, Bagby SM, Kar A, Pozdeyev N, Xu M, Gowan K, Sharma V, Raeburn CD, Albuja-Cruz M, Jones KL, Fishbein L, Schweppe RE, Somerset H, Pitts TM, Leong S, Wierman ME. 2018. Development of new preclinical models to advance adrenocortical carcinoma research. Endocrine-Related Cancer 25(4):437-451 DOI 10.1530/erc-17-0447.
Lee YS, Hwang SG, Kim JK, Park TH, Kim YR, Myeong HS, Kwon K, Jang CS, Noh YH, Kim SY. 2015. Topological network analysis of differentially expressed genes in cancer cells with acquired gefitinib resistance. Cancer Genomics & Proteomics 12(3):153.
Liu Q, Guntuku S, Cui XS, Matsuoka S, Cortez D, Tamai K, Luo G, Carattini-Rivera S, DeMayo F, Bradley A, Donehower LA, Elledge SJ. 2000. Chk1 is an essential kinase that is regulated by Atr and required for the G(2)/M DNA damage checkpoint. Genes & Development 14(12):1448-1459.
Liu X, Hu AX, Zhao JL, Chen FL. 2017. Identification of key gene modules in human osteosarcoma by co-expression analysis weighted gene co-expression network analysis (WGCNA). Journal of Cellular Biochemistry 118(11):3953-3959 DOI 10.1002/jcb.26050.
Liu J, Jing L, Tu X. 2016. Weighted gene co-expression network analysis identifies specific modules and hub genes related to coronary artery disease. BMC Cardiovascular Disorders 16(1):54 DOI 10.1186/s12872-016-0217-3.
Liu Z, Li M, Fang X, Shen L, Yao W, Fang Z, Chen J, Feng X, Hu, Zeng Z, Lin C, Weng J, Lai Y, Yi G. 2018. Identification of surrogate prognostic biomarkers for allergic asthma in nasal epithelial brushing samples by WGCNA. Journal of Cellular Biochemistry 120(4):5137-5150 DOI 10.1002/jcb.27790.
Lombardi CP, Raffaelli M, Pani G, Maffione A, Princi P, Traini E, Galeotti T, Rossi ED, Fadda G, Bellantone R. 2006. Gene expression profiling of adrenal cortical tumors by cDNA macroarray analysis. Results of a preliminary study. Biomedicine & Pharmacotherapy 60(4):186-190 DOI 10.1016/j.biopha.2006.03.006.
Maire V, Baldeyron C, Richardson M, Tesson B, Vincent-Salomon A, Gravier E, Marty- Prouvost B, De Koning L, Rigaill G, Dumont A, Gentien D, Barillot E, Roman-Roman S, Depil S, Cruzalegui F, Pierré A, Tucker GC, Dubois T. 2015. TTK/hMPS1 is an attractive therapeutic target for triple-negative breast cancer. PLOS ONE 8(5):e63712 DOI 10.1371/journal.pone.0063712.
Manic G, Signore M, Sistigu A, Russo G, Corradi F, Siteni S, Musella M, Vitale S, De Angelis ML, Pallocca M, Amoreo CA, Sperati F, Di Franco S, Barresi S, Policicchio E, De Luca G, De Nicola F, Mottolese M, Zeuner A, Fanciulli M, Stassi G, Maugeri-Saccà M, Baiocchi M, Tartaglia M, Vitale I, De Maria R. 2017. CHK1-targeted therapy to deplete DNA replication-
stressed, p53-deficient, hyperdiploid colorectal cancer stem cells. Gut 67(5):903-917 DOI 10.1136/gutjnl-2016-312623.
Mao L, Van Hemert JL, Dash S, Dickerson JA. 2009. Arabidopsis gene co-expression network and its functional modules. BMC Bioinformatics 10(1):346 DOI 10.1186/1471-2105-10-346.
Monterisi S, D’Ario G, Dama E, Rotmensz N, Confalonieri S, Tordonato C, Troglio F, Bertalot G, Maisonneuve P, Viale G, Nicassio F, Vecchi M, Di Fiore PP, Bianchi F. 2015. Mining cancer gene expression databases for latent information on intronic microRNAs. Molecular Oncology 9(2):473-487 DOI 10.1016/j.molonc.2014.10.001.
R Core Team. 2013. R: a language and environment for statistical computing. Version 3.0.2. Vienna: R Foundation for Statistical Computing. Available at https://www.R-project.org/.
R Core Team. 2015. R: a language and environment for statistical computing. Version 3.1.3. Vienna: R Foundation for Statistical Computing. Available at https://www.R-project.org/.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. 2015. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43(7):e47 DOI 10.1093/nar/gkv007.
Schuler F, Weiss JG, Lindner SE, Lohmüller M, Herzog S, Spiegl SF, Menke P, Geley S, Labi V, Villunger A. 2017. Checkpoint kinase 1 is essential for normal B cell development and lymphomagenesis. Nature Communications 8(1):1697 DOI 10.1038/s41467-017-01850-4.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. 2003. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Research 13(11):2498-2504.
Slater EP, Diehl SM, Langer P, Samans B, Ramaswamy A, Zielke A, Bartsch DK. 2006. Analysis by cDNA microarrays of gene expression patterns of human adrenocortical tumors. European Journal of Endocrinology 154(4):587-598 DOI 10.1530/eje.1.02116.
Soon PS, McDonald KL, Robinson BG, Sidhu SB. 2008. Molecular markers and the pathogenesis of adrenocortical cancer. Oncologist 13(5):548-561 DOI 10.1634/theoncologist.2007-0243.
Sun X, Clermont PL, Jiao W, Helgason CD, Gout PW, Wang Y, Qu S. 2016. Elevated expression of the centromere protein-A(CENP-A)-encoding gene as a prognostic and predictive biomarker in human cancers. International Journal of Cancer 139(4):899-907 DOI 10.1002/ijc.30133.
Sun Q, Zhao H, Zhang C, Hu T, Wu J, Lin X, Luo D, Wang C, Meng L, Xi L, Li K, Hu J, Ma D, Zhu T. 2017. Gene co-expression network reveals shared modules predictive of stage and grade in serous ovarian cancers. Oncotarget 8(26):42983-42996 DOI 10.18632/oncotarget.17785.
Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, Von Mering C. 2015. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Research 43(D1):D447-D452 DOI 10.1093/nar/gku1003.
Tachibana KE, Gonzalez MA, Coleman N. 2005. Cell cycle dependent regulation of DNA replication and its relevance to cancer pathology. Journal of Pathology 205(2):123-129 DOI 10.1002/path.1708.
Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. 2017. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Research 45(W1):W98-W102 DOI 10.1093/nar/gkx247.
Thu KL, Silvester J, Elliott MJ, Ba-Alawi W, Duncan MH, Elia AC, Mer AS, Smirnov P, Safikhani Z, Haibe-Kains B, Mak TW, Cescon DW. 2018. Disruption of the anaphase-promoting complex confers resistance to TTK inhibitors in triple-negative breast cancer. Proceedings of the
PeerJ
National Academy of Sciences of the United States of America 115(7):E1570-E1577 DOI 10.1073/pnas.1719577115.
Villman K, Sjöström J, Heikkilä R, Hultborn R, Malmström P, Bengtsson NO, Söderberg M, Saksela E, Blomqvist C. 2006. TOP2A and HER2 gene amplification as predictors of response to anthracycline treatment in breast cancer. Acta Oncologica 45(5):590-596 DOI 10.1080/02841860500543182.
Wang J, Xu B, Yuan P, Zhang P, Li Q, Ma F, Fan Y. 2012. TOP2A amplification in breast cancer is a predictive marker of anthracycline-based neoadjuvant chemotherapy efficacy. Breast Cancer Research & Treatment 135(2):531-537 DOI 10.1007/s10549-012-2167-5.
Yang Y, Han L, Yuan Y, Li J, Hei N, Liang H. 2014. Gene co-expression network analysis reveals common system-level properties of prognostic genes across cancer types. Nature Communications 5(1):3231 DOI 10.1038/ncomms4231.
Yang Q, Wang R, Wei B, Peng CG, Wang L, Hu GZ, Kong DL, Du C. 2018. Candidate biomarkers and molecular mechanism investigation for glioblastoma multiforme utilizing WGCNA. BioMed Research International 2018:1-10 DOI 10.1155/2018/4246703.
Zhang Y, Hunter T. 2013. Roles of Chk1 in cell biology and cancer therapy. International Journal of Cancer 134(5):1013-1023 DOI 10.1002/ijc.28226.
Zheng S, Cherniack AD, Dewal N, Moffitt RA, Danilova L, Murray BA, Lerario AM, Else T, Knijnenburg TA, Ciriello G, Kim S, Assie G, Morozova O, Akbani R, Shih J, Hoadley KA, Choueiri TK, Waldmann J, Mete O, Robertson AG, Wu HT, Raphael BJ, Shao L, Meyerson M, Demeure MJ, Beuschlein F, Gill AJ, Sidhu SB, Almeida MQ, Fragoso MCBV, Cope LM, Kebebew E, Habra MA, Whitsett TG, Bussey KJ, Rainey WE, Asa SL, Bertherat J, Fassnacht M, Wheeler DA. 2016. Comprehensive pan-genomic characterization of adrenocortical carcinoma. Cancer Cell 29(5):723-736.
Zhu D, Xu S, Deyanat-Yazdi G, Peng SX, Barnes LA, Narla RK, Tran T, Mikolon D, Ning Y, Shi T, Jiang N, Raymon HK, Riggs JR, Boylan JF. 2018. Synthetic lethal strategy identifies a potent and selective TTK and CLK1/2 inhibitor for treatment of triple-negative breast cancer with a compromised G1-S checkpoint. Molecular Cancer Therapeutics 17(8):1712-1738 DOI 10.1158/1535-7163.mct-17-1084.
Zuo Z, Shen JX, Pan Y, Pu J, Li YG, Shao XH, Wang WP. 2018. Weighted gene correlation network analysis (WGCNA) detected loss of MAGI2 promotes chronic kidney disease (CKD) by podocyte damage. Cellular Physiology and Biochemistry 51(1):244-261 DOI 10.1159/000495205.