frontiers in Genetics
Check for updates
Screening and Identification of Potential Prognostic Biomarkers in Adrenocortical Carcinoma
Wen-Hao Xu1,21, Junlong Wu1,21, Jun Wang1,21, Fang-Ning Wan1,2, Hong-Kai Wang1,2, Da-Long Cao1,2, Yuan-Yuan Qu1,2*, Hai-Liang Zhang1,2* and Ding-Wei Ye1,2*
1 Department of Urology, Fudan University Shanghai Cancer Center, Shanghai, China, 2 Department of Oncology, Shanghai Medical College, Fudan University, Shanghai, China
Objective: Adrenocortical carcinoma (ACC) is a rare but aggressive malignant cancer that has been attracting growing attention over recent decades. This study aims to integrate protein interaction networks with gene expression profiles to identify potential biomarkers with prognostic value in silico.
Methods: Three microarray data sets were downloaded from the Gene Expression Omnibus (GEO) database to identify differentially expressed genes (DEGs) according to the normalization annotation information. Enrichment analyses were utilized to describe biological functions. A protein-protein interaction network (PPI) of the DEGs was developed, and the modules were analyzed using STRING and Cytoscape. LASSO Cox regression was used to identify independent prognostic factors. The Kaplan-Meier method for the integrated expression score was applied to analyze survival outcomes. A receiver operating characteristic (ROC) curve was constructed with area under curve (AUC) analysis to determine the diagnostic ability of the candidate biomarkers.
Results: A total of 150 DEGs and 24 significant hub genes with functional enrichment were identified as candidate prognostic biomarkers. LASSO Cox regression suggested that ZWINT, PRC1, CDKN3, CDK1 and CCNA2 were independent prognostic factors in ACC. In multivariate Cox analysis, the integrated expression scores of the modules showed statistical significance in predicting disease-free survival (DFS, P = 0.019) and overall survival (OS, P < 0.001). Meanwhile, ROC curves were generated to validate the ability of the Cox model to predict prognosis. The AUC index for the integrated genes scores was 0.861 (P < 0.0001).
Conclusion: In conclusion, the present study identifies DEGs and hub genes that may be involved in poor prognosis and early recurrence of ACC. The expression levels of ZWINT, PRC1, CDKN3, CDK1 and CCNA2 are of high prognostic value, and may help us understand better the underlying carcinogenesis or progression of ACC. Further studies are required to elucidate molecular pathogenesis and alteration in signaling pathways for these genes in ACC.
Keywords: adrenocortical carcinoma, bioinformatics analysis, biomarker, prognosis, network module
OPEN ACCESS Edited by: George Bebis, University of Nevada, United States
Reviewed by: Pawel Buczkowicz, Gene42, Inc., Canada Maria Candida Barisson Villares Fragoso, University of São Paulo, Brazil Reju Korah, Yale University, United States
*Correspondence: Yuan-Yuan Qu quyy1987@163.com Hai-Liang Zhang zhanghl918@163.com
Ding-Wei Ye dwyeli@163.com
+These authors have contributed equally to this work
Specialty section:
This article was submitted to Cancer Genetics, a section of the journal Frontiers in Genetics
Received: 21 December 2018 Accepted: 08 August 2019 Published: 11 September 2019
Citation:
Xu W-H, Wu J, Wang J, Wan F-N, Wang H-K, Cao D-L, Qu Y-Y, Zhang H-L and Ye D-W (2019) Screening and Identification of Potential Prognostic Biomarkers in Adrenocortical Carcinoma. Front. Genet. 10:821.
doi: 10.3389/fgene.2019.00821
INTRODUCTION
Adrenocortical carcinoma (ACC) is a rare endocrine malignancy with an annual incidence of 0.7-2.0 per million people, accounting for an estimated 0.02% of all cancers (Wajchenberg et al., 2000; Kebebew et al., 2006; Kerkhofs et al., 2013). Although comparatively uncommon, ACC patients often face aggressive progression, with merely less than 35% of patients surviving 5 years after initial diagnosis (Else et al., 2014). Currently, the preferred treatment regimen for ACC is surgical resection of the primary tumor (Fassnacht et al., 2013). However, almost half of ACC patients have disseminated metastasis, and approximately one-third of patients have locoregional metastases after surgery (Else et al., 2014). The first-line treatment, and the only ACC- specific medical therapy approved by the US Food and Drug Administration, is Mitotane, which is regularly used as an adjuvant agent in these patients (Else et al., 2014). Mitotane disrupts mitochondria and activates an apoptotic process (Poli et al., 2013). A major concern of the therapeutic management with Mitotane is the risk of toxicity, which may lead to severe adrenal insufficiency (Paragliola et al., 2018).
Accumulating evidence has demonstrated that gene expression levels and related pathways are involved in the carcinogenesis and progression of ACC. For example, the most frequent alterations observed in ACC are overexpression of insulin-like growth factor 2 (IGF-2) (Gicquel et al., 2001; Giordano et al., 2003; de Fraipont et al., 2005) and constitutive activation of the Wnt/B-Catenin pathway (Gaujoux et al., 2011). Despite these encouraging advances in ACC clinical strategies, only a minority of patients receive any significant survival benefit because of a lack of effective therapeutic strategies (Mohan et al., 2018). Therefore, it is crucial to understand the underlying molecular mechanisms involved in the carcinogenesis, proliferation and recurrence of ACC and thus develop effective diagnostic and therapeutic strategies.
Over the last decade, microarray technologies and bioinformatic analysis have been widely used to detect comprehensive mRNA expression levels, which have assisted in identifying the differentially expressed genes (DEGs) and functional pathways involved in the tumorigenesis and progression of ACC. However, because of the rarity of this tumor, there has been a problem in identifying potential markers to differentiate ACC from other renal neoplasms, and thus guiding potential treatment strategy. In the present study, three mRNA microarray datasets were downloaded from GEO database and analyzed to obtain DEGs between cancer tissues and adjacent normal tissues. Subsequently, functional pathway enrichment analyses were implemented to further understand the molecular mechanisms underlying carcinogenesis. The protein-protein interaction (PPI) network reveals the functions of all proteins and the importance
Abbreviations: ACC, adrenocortical carcinoma; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction; TCGA, the Cancer Genome Atlas; HPF, high power field; DFS, disease-free survival; OS, overall survival; HR, hazard ratio; CI, confidence interval; ROC, receiver operating characteristic curve; AUC, area under curve; GSEA, Gene set enrichment analysis.
of these interactions with regards to biological processes, molecular functions, and signal transduction (Sharan et al., 2007; Wu et al., 2009; Bapat et al., 2010). This may provide insights into the mechanisms of generation or development of diseases.
To investigate candidate biomarkers in tumor tissue and to define their value in ACC patients, this work focuses on analyzing the gene expression profiles, revealing the underlying biological interaction networks and assessing their prognostic value. We hypothesize that the oncogenic activity of significant hub genes correlates with poor prognosis, and might reveal potential prognostic markers and therapeutic targets for ACC.
MATERIALS AND METHODS
Raw Biological Microarray Data
The raw DNA microarray data were obtained from GEO (http:// www.ncbi.nlm.nih.gov/geo) (Edgar et al., 2002) for patients with ACC. Corresponding genes converted into the probes were converted into symbols according to the annotation information in the platform. Three chip data sets GSE14922, GSE19750 and GSE90713 (4 normal and 4 ACC samples in GSE14922, 4 normal and 44 ACC samples in GSE19750, and 5 normal and 58 ACC samples in GSE90713) were downloaded from GEO (Agilent GPL6480 platform, Affymetrix GPL570 platform and Affymetrix GPL15270 platform, respectively).
Normalization and Elucidation of DEGs
DNA microarray analysis begins with preprocessing and normalization of raw biological data. This process removes noise from the biological data and ensures its integrity. Next, background correction of probe data, normalization, and summarization were executed by robust multi-array average analysis algorithm17 in affy package of R.
The DEGs between ACC and non-cancerous samples were screened and identified across experimental conditions. Delineating parameters such as adjusted P-values (adj. P), Benjamini and Hochberg false discovery rate (FDR) and fold change were utilized for filtering of DEGs and applied to provide a balance between discovery of statistically significant genes and limitations of false- positives. Probe sets without corresponding gene symbols or genes with more than one probe set were removed or averaged. Log2FC (fold change) > 1 and adj. P-value <0.01 were considered statistically significant.
Functional Enrichment of DEGs
Discerning the role of DEGs in ACC, biological attributes including biological processes (BP), molecular functions (MF), and cellular components (CC) were extracted from Gene Ontology (GO) enrichment analysis (Ashburner et al., 2000). Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2016) is a database resource for understanding high-level functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies. The online Database for Annotation, Visualization
and Integrated Discovery (DAVID; https://david-d.ncifcrf. gov/summary.jsp Version 6.8) was used to explore the role of development-related signaling pathways in ACC (Huang et al., 2007). P-value < 0.05 was considered statistically significant. GO enrichment was analyzed and displayed using a bubble chart.
PPI Network Construction and Module Analysis
In the present study, the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 10.0) online database was used to predict PPI network of DEGs and analyze the functional interactions between proteins (Franceschini et al., 2013). An interaction with a combined score >0.4 was considered statistically significant.
Cytoscape (version 3.5), an open source bioinformatics software platform, was used to visualize molecular interaction networks (Smoot et al., 2011). Molecular Complex Detection (MCODE) (version 1.4.2) is a plug-in for Cytoscape used for clustering a given network based on topology to find densely connected regions (Bandettini et al., 2012). MCODE could identify the most significant module in the PPI networks with selection as follows: MCODE scores >5, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100 and k-score = 2. Subsequently, the KEGG and GO analyses for genes in this module were performed using DAVID.
Hub Genes Selection and Analysis
The hub nodes of network with connectivity degrees >10 were identified. A network of the 24 genes and their co-expression genes was analyzed using cBioPortal (http://www.cbioportal. org) online platform (Cerami et al., 2012). ClueGO is a Cytoscape plug-in that visualizes the non-redundant biological terms for large clusters of genes in a functionally grouped network (Bindea et al., 2009). The biological process from GO and KEGG pathway analysis of hub genes was performed and visualized using ClueGO (version 2.5.3) and CluePedia (version 1.5.3), a functional extension of ClueGO, plug-in of Cytoscape (Bindea et al., 2013). Potential coexpression relationship between the 24 hub genes and possible prognostic value are shown in a heat map.
Statistical Analysis
Phenotype and expression profiles of hub genes in 76 ACC patients from TCGA were analyzed and displayed. Clinical and pathological parameters of the cohort were summarized. Expression of hub genes was respectively identified as binary variables (high vs. low) referring to median expression of each hub gene in the TCGA cohort. Then, a LASSO Cox regression model was constructed to find independent prognostic factors. The significant hub gene expression profiles of common neoplasm were analyzed and displayed using Oncomine online database (http://www.oncomine.com) (Giordano et al., 2003; Giordano et al., 2009).
The Kaplan-Meier method was applied to analyze survival differences between groups. The primary end point was overall
survival (OS) for patients, which was evaluated from the date of first therapy to the date of death or last follow-up. Disease- free survival (DFS), as the secondary end point, was the length of time from the initiation of curative treatment to the date of progression or the start date of a second-line treatment or the date of death, whichever occurred first. The follow-up duration was estimated using the Kaplan-Meier method with 95% confidence intervals (95%CI) and log-rank test in separate curves. Univariate analyses were performed with Cox logistic regression models to find independent variables, including age at diagnosis, gender, laterality, TNM stage, pathologic stage, mitotic rate, invasion of tumor capsule, sinusoid invasion, necrosis, Weiss score, new tumor event after first treatment and integrated expression score. Parameters with P-value less than 0.1 were enrolled in multivariate Cox regression analyses of DFS and OS in “Back-LR” method. Integrated score was identified as the sum of the weight of each significant hub gene. X-tile software was utilized to take the cut-off value. All hypothetical tests were two- sided and P-values less than 0.05 were considered significant in all tests. The receiver operating characteristic curve (ROC) was constructed by predicting the probability of a diagnosis being of high or low integrated score of significant hub gene expression. Area under curve (AUC) analysis was performed to determine the diagnostic ability.
Sensitivity Analysis of Chip Datasets
In this study, GSE14922 only contains 4 ACC patients and 4 normal patients, while datasets 2 and 3 have 44 and 58 ACC samples respectively. To avoid penitential bias and see what other significant genes may have been missed, the analysis was re-run without GSE14922. Prognostic values of DEGs were then also tested against the TCGA validation cohort.
Data Processing of Gene Set Enrichment Analysis (GSEA)
TCGA database was implemented with the GSEA method using the Category version 2.10.1 package. For each separate analysis, Student’s-t-test statistical score was performed in consistent pathways, and the mean of the differential expression genes was calculated. A permutation test of 1000 times was used to identify the significantly changed pathways. The adjusted P values (adj. P) using Benjamini and Hochberg (BH) false discovery rate (FDR) method by default were applied to correct the occurrence of false positive results (Subramanian et al., 2005). The significant related genes were defined with an adj. P less than 0.01 and FDR less than 0.25. Statistical analysis and graphical plotting were conducted using R software (Version 3.3.2).
RESULTS
This study consisted of three stages. In the first stage, we assessed DEGs using three datasets hosted on the GEO platform. In the second stage, coexpression, functional annotation of hub genes and patient survival analysis were carried out. In the third stage,
the most significant hub genes were selected, evaluated and integrated to predict their prognostic value.
Identification of DEGs in ACC
After standardization and identification of the microarray results, the DEGs (1,804 probe samples with 1,539 DEGs in GSE14922, 2,454 probe samples with 2,040 DEGs in GSE19750 and 1,216 probe samples with 806 DEGs in GSE90713) were determined to be significant based on the analysis and the statistical parameters of the data processing steps. The overlap among the three datasets included 150 significant genes and is displayed in the Venn diagram in Figure 1A.
GO and KEGG Enrichment Assessment of DEGs
To analyze the biological classification of the DEGs, functional and pathway enrichment analyses were performed using DAVID. As shown in Supplementary Figure 1, gene ontology (GO) analysis indicated that changes in the biological processes of the DEGs were significantly associated with the mitotic cell cycle, cell cycle process, movement of cells or subcellular components and cell locomotion activity. Changes in molecular function were mostly enriched in growth factor binding, kinase activity, extracellular matrix structure constituents and insulin-like growth factor binding. Changes in
A
1055
GSE14922
195
150
139
GSE19750
GSE90713
1531
165
352
ORCS
SHID19
ERPINB
FANCI
EPHX2
KIFJA
SETBP!
STAB1
CIAA010
VSIG4
B
PARVA
AP151
NFIA
BIRC5
AURKA
TYMS
MAPK13
PRCI
NOCBO
GGH
MARCO
TNS1
EZH2
CDKN3
COK1
DHFR
AXL
TOP2A
MADEL1
CENPN
CIR
SPTBN1
DNAJB!
COLEC1
CCNA2
TPX2
CASP9
NASEH2
BUB1B
CENPF
TMPO
NCAPO
CONBI
IGFBPS
PAPPA
GFBP4
RSPO3
ANLN
MNDI
SLCMBAB
ZWINT
KONQ1
PLAT
RACGAP
GINS1
MYOTA
ERPING
PDGFRA
CTH
BRE
CONE1
COKNIC
TIMP3
GGTS
PTPRB
NEXN
TLR4
HGF
GPXB
ZNF367
NR4A2
PLN
TEK
CXCL12
TAGLN
MYLK
RRC32
NR2F1
FBLN1
GORBS1
DGATI
LMOD1
ANTXR1
CAB3OL
TLE2
PRELP
ADHIB
VASN
MFAPS
C
SRPX2
PTHIR
FANCI
FBLNS
AOX1
EFEMP2
BIRCS
CIAAD10
CDK1
TYMS
NDCAO
AURKA
PRC1
EZH2
TOP2A
CENPN
CDKN3
CCNA2
MAD2L
CENPF
TPX2
NASIEHQ
CAPG
ZWINT
80816
MND1
ANLN
CÒNB1
ACGAP
FIGURE 1 | Venn diagram, PPI network and the most significant module of DEGs. (A) DEGs were selected with a fold change >2 and P-value <0.01 among the mRNA expression profiling chip datasets GSE14922, GSE19750 and GSE90713. The 3 datasets show an overlap of 150 genes in the Venn diagram. (B) The PPI network of DEGs was constructed using Cytoscape. (C) The most significant module was obtained from PPI network with 24 nodes. Significant edges are marked in light blue with a K-score >0.800.
cellular components were mainly enriched in the chromosome, centromeric region, extracellular region, actin cytoskeleton and mitotic spindle. KEGG pathway analysis revealed that the DEGs were mainly enriched in cell cycle, progesterone-mediated oocyte maturation, oocyte meiosis, arachidonic acid metabolism and the p53 signaling pathway, summarized in Table 1.
PPI Network Establishment and Module Analysis
We constructed the PPI network of the DEGs (Figure 1B) and subsequently found the most significant module penal using a Cytoscape plugin (Figure 1C). The enrichment profiles from DAVID functional analyses of the 24 hub genes suggested that the hub genes in this module were primarily enriched in cell cycle phase, M phase, the mitotic cell cycle and mitosis (Table 2).
Hub Gene Selection and Analysis
After statistical selection, the significant hub nodes of the network included RACGAP1, AURKA, KIAA0101, MAD2L1, ZEH2, CCNB1, BIRC5, ZWINT, NDC80, NCAPG, TOP2A, PRC1,
| Term | Description | Count in gene set | P value |
|---|---|---|---|
| Has04110 | Cell cycle | 7 | 7.01E-04 |
| Has04914 | Progesterone-mediated oocyte maturation | 5 | 6.05E-03 |
| Has04114 | Oocyte meiosis | 5 | 0.01757 |
| Has00590 | Arachidonic acid metabolism | 4 | 0.01758 |
| Has04115 | p53 signaling pathway | 4 | 0.02312 |
| Has05133 | Pertussis | 4 | 0.02667 |
| Has06161 | Hepatitis B | 5 | 0.04010 |
| Has00380 | Tryptophan metabolism | 3 | 0.04228 |
KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; ACC, adrenocortical carcinoma.
| Term | Description | Count in gene set | P value |
|---|---|---|---|
| GO:0022403 | Cell cycle phase | 16 | 6.836E-19 |
| GO:0022402 | Cell cycle process | 17 | 1.154E-18 |
| GO:0000279 | M phase | 15 | 1.898E-18 |
| GO:0000278 | Mitotic cell cycle | 15 | 9.924E-18 |
| GO:0007067 | Mitosis | 13 | 6.478E-17 |
| GO:0000280 | Nuclear division | 13 | 6.477E-17 |
| GO:0000087 | M phase of mitotic cell cycle | 13 | 8.063E-17 |
| GO:0005819 | Spindle | 9 | 4.531E-11 |
| GO:0015630 | Microtubule cytoskeleton | 11 | 4.532E-9 |
| GO:0005694 | Chromosome | 10 | 3.922E-8 |
| hsa04110 | Cell cycle | 5 | 2.268E-5 |
| hsa04914 | Progesterone-mediated oocyte maturation | 4 | 2.461E-4 |
| hsa04114 | Oocyte meiosis | 4 | 5.097E-4 |
GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes.
CENPF, CENPN, FANCI, CDKN3, MND1, RNASEH2A, TYMS, CDK1, BUB1B, CCNA2, TPX2 and ANLN. A visual network of the 24 genes and their coexpressed genes was set up (Figure 2A). The GO biological processes and KEGG functional annotation analysis of the hub genes are shown in Figure 2B. The detailed functional notes and classification pie charts are provided in the Supplementary Figure 2. Of the GO biological processes, 66.67% of terms belonged to the mitotic cell cycle checkpoint, 15.79% to mitotic spindle organization, 12.28% to anaphase- promoting complex-dependent catabolic processes, 3.51% to protein localization to kinetochore, and 1.75% to chromosome condensation. A heat map shows that a potential coexpression relationship may exist between the 24 hub genes, which could suggest they have value for prognostic prediction (Figure 2C).
Clinicopathological Statistical Analysis
The clinical and pathological parameters from phenotype and expression profiles of the hub genes in 76 ACC patients from The Cancer Genome Atlas (TCGA) are summarized in Table 3. Each hub gene was classified into dichotomous variables according to the median expression in the analysis. Subsequently, the univariate survival analysis of the hub genes was performed using a Kaplan- Meier curve. Apart from MND1, ACC patients with elevated expression of the other 23 hub genes showed significantly worse OS and DFS (Supplementary Figure 3). LASSO Cox regression suggested that ZWINT, PRC1, CDKN3, CDK1 and CCNA2 are significant weighted prognostic factors, and that an integrated gene panel may serve as an independent penal in ACC samples. The five significant hub gene expression profiles showed significantly elevated expression in tumor tissues compared with the corresponding normal tissues (Figures 3A-E). In addition, differential analysis from the ONCOMINE online database of tumor and normal tissue in two cohorts indicated that ZWINT, PRC1, CDKN3, CDK1 and CCNA2 were highly expressed in ACC samples (Figure 3F). Elevated expression patterns were significantly associated with distant metastasis, necrosis, Weiss score and mitotic rate of >5 mitoses per 50 high power fields (HPF), plotted in Figure 4.
Cox Regression Analyses and Survival Outcomes of the Cohorts
In this study, the integrated expression score was identified as the sum of the weight of each binary gene expression. In univariate models, traditional prognostic factors, specifically T stage, M stage and pathologic stage, were significantly correlated with DFS (P < 0.001) and OS (P < 0.001) in ACC patients. Importantly, in univariate Cox regression analyses of DFS, subgroups of integrated expression score (High vs. Low) showed that integrated gene expression amplification significantly correlated with poor DFS (P < 0.001) for ACC patients. In addition, mitotic rate (≤5/50 HPF vs. >5/50 HPF) (P = 0.022), necrosis (Present vs. Absent) (P= 0.011), Weiss score(≤3 vs.>3) (P=0.008) and new tumor event (Present vs. Absent) (P <0.001) were correlated with poor DFS. In univariate Cox regression analyses of OS, invasion of tumor capsule (Present vs. Absent) (P = 0.013), sinusoid invasion (Present vs. Absent) (P = 0.011), necrosis (Present vs. Absent) (P = 0.026) and new tumor event
DCTN4
A
RACGAP1
B
Iver reomneration
CSNK1G3
AKT2
RAC1
aima orogn ergeneration
PPP2CA
CENPF
PIP5K1C
U2AF2
chromosome condensation
BIRC5
Progester firma dlaind oocyter
SKP
MADA
A
NP37
NUP107
negative regulation of 00M transition of miotic cell Cycle
histone phosphory/a
198260
RAD21
NCAPG
prchnoting eperkjent
89MG4
woce
NUPA
CSNK1A1
REDZ
ORVAR
CEUN
DNM1L
SUBSS
MAPK9
negative negu metaphaseiana
PSMOS
CONB
SORT
PRC1
protein localization lo chromosome,
centromenciongion
ODK1
mtotc cell cycle embly Checkpoint
ANAPCS
EZH2
COMEN
MAPSK1C
RAD17PSI 98
TOP2A
LIG1
CDKZ
ttott song
MBD3
CEK4
MND1
protein localization to kinetochore
SUPTEH
DYRK2
mitotic spindle organization
RAB5B
FANCI
RFC
AURKA
TAK2
segregation
water chromand singingation
kanelochore organization
organization involved in mitosis
Melachase svale congressos@lachment of soindie microtubules to kinelochare
MO
T4
C
MO
Stage II
MO
T2
T2
MO
MO
T4
T2
MO
Stage II
T2
MO
Stage II
T2
MO
Stage II
T2
RACGAP1
AURKA
KIAA0101
MAD2L1
EZH2
CCNB1
BIRC5
ZWINT
NDC80
NCAPG
TOP2A
PRC1
CENPF
CENPN
FANCI
CDKN3
MND3
RNASEH2A
TYMS
CDK1
BUB1B
CCNA2
TPX2
ANLN
Stage IV
M1
T4
Stage III Stage II Stage !
MO
T3 T2
T1
low
high
FIGURE 2 | Interaction network and biological process analysis of the hub genes. (A) Hub genes and their co-expression network were analyzed using cBioPortal. Nodes with bold black outline represent hub genes. Nodes with thin black outline represent the co-expression genes. (B) The biological process analysis of hub genes was constructed using ClueGO. Different color of nodes refers to the functional annotation of ontologies. Corrected P-value < 0.01 was considered statistically significant. (C) Hierarchical partitioning of 24 hub genes was obtained from DNA microarrays. It represent the level of expression of 24 genes across a number of comparable samples with high expression samples marked in red and low in blue.
(Present vs. Absent) (P < 0.001) were associated with shorter OS. However, in multivariate prognostic analysis, new tumor event after first treatment (P < 0.001) and integrated expression score (P = 0.019) were statistically significant parameters in predicting DFS (Table 4). Age at diagnosis (P = 0.003), M stage (P = 0.033), new tumor event after first treatment (P = 0.013) and integrated expression score (P < 0.001) were significantly associated with shorter OS (Table 5).
*Multivariate Cox regression analyses of OS in 76 enrolled ACC patients was run in “Back-LR” method. After integrating all the significant gene expression profiles in the Cox regression models, the Kaplan-Meier method was used to determine the significant survival outcomes (DFS: P < 0.0001; OS: P < 0.0001), shown in Figures 5A, B. Meanwhile, ROC curves were generated to validate the ability of the logistic model to predict prognosis. The AUC index for the integrated gene scores was 0.861 (P < 0.0001) (Figure 5C).
Sensitivity Analysis of Chip Datasets
To avoid penitential bias and see what other significant genes may have been missed, two datasets GSE19750 and GSE90713 were enrolled to re-run the analysis. The overlap among the two datasets, which includes 315 differential expressed genes (DEGs), is displayed in the Venn diagram (Supplementary Figure 4A). A PPI network of the new DEGs was constructed in Supplementary Figure 4B. Subsequently, we selected the most significant module penal using M-CODE, a plug-in of Cytoscape, and found 31 hub genes including NDC80, MND1, MAD2L1, UBE2C, NCAPG, GINS1, CENPN, CDKN3, CCNA2, ZWINT, BIRC5, KIAA0101, TOP2A, BUB1B, CCNB1, AURKA, SMC2, ATAD2, PRC1, TPX2, CDK1, RACGAP1, TYMS, ANLN, PRIM1, NUSAP1, CENPF, SPAG5, SMC4, EZH2, FANCI. Interestingly, five significant DEGs we have focused on (ZWINT, PRC1, CDKN3, CDK1, CCNA2) still consist of this new module penal (Supplementary Figure 4C), indicating a good stability of our molecular model. Eight different
| Characteristics | Entire cohort (N = 76) |
|---|---|
| N (%) | |
| Age, years | |
| ≤57 | 54 (71.1) |
| >58 | 22 (28.9) |
| Gender | |
| Male | 30 (39.5) |
| Female | 46 (60.5) |
| Germline testing performed | |
| Present | 12 (18.2) |
| Absent | 54 (81.1) |
| Laterality | |
| Left | 42 (55.3) |
| Right | 34 (44.7) |
| pT stage | |
| T1 -T2 | 50 (65.8) |
| T3 -T4 | 26 (34.2) |
| pN stage | |
| N0 | 68 (89.5) |
| N1 | 8 (10.5) |
| M stage | |
| M0 | 62 (81.6) |
| M1 | 14 (18.4) |
| Pathologic stage | |
| I- II | 46 (60.5) |
| III - IV | 30 (39.5) |
| Histological type | |
| Myxiod | 1 (1.3) |
| Oncocytic | 3 (3.9) |
| Usual | 72 (94.7) |
| Mitotic rate | |
| ≤5/50 HPF | 26 (38.8) |
| >5/50 HPF | 41 (61.2) |
| Invasion of tumor capsule | |
| Absent | 30 (42.9) |
| Present | 40 (57.1) |
| Sinusoid invasion | |
| Absent | 34 (58.6) |
| Present | 24 (41.4) |
| Necrosis | |
| Absent | 17 (23.9) |
| Present | 54 (76.1) |
| Weiss score | |
| ≤3 | 17 (22.0) |
| >4 | 47 (78.0) |
| Persistent distant metastasis | |
| Absent | 56 (73.7) |
| Present | 20 (26.3) |
ACC, adrenocortical carcinoma; TCGA, the Cancer Genome Atlas; HPF, high power field.
DEGs are found different from these in three-chipset study, including UBE2C, GINS1, SMC2, ATAD2, PRIM1, NUSAP1, SPAG5, SMC4. Kaplan-Meier method was used to analyze mRNA expression level of 8 hub genes in TCGA cohort, which also showed statistically significant correlation with progressive progression and poor prognosis (Supplementary Figure 5).
Significant Genes and Pathways Obtained by GSEA
A total of 100 significant genes were obtained by gene set enrichment analysis (GSEA) with positive and negative
correlation. Importantly, GSEA was used to perform hallmark analysis for ZWINT, PRC1, CDKN3, CDK1 and CCNA2. This suggested that the most involved significant pathways included mitotic spindle, G2M checkpoint and E2F targets. The details are shown in Figure 6.
DISCUSSION
Adrenocortical carcinoma (ACC) is a rare but aggressive cancer, with a typically high incidence in children with a TP53 germline mutation (Fassnacht et al., 2013). The Wnt/B-catenin pathway and IGF-2 signaling have been confirmed as altered signaling pathways in ACC patients, while increasing data indicate that the available evidence is inadequate for malignant phenotype and poor prognosis (Berthon et al., 2010; Heaton et al., 2012), especially for the diagnosis of low-grade ACC confined to the adrenal gland (Mete et al., 2018). Although there are diagnostic and prognostic molecular tests for ACC such as the IGF-2, Ki-67, p53, BUB1B, PBK, HURP, NEK2, DAX, Wnt/ß-catenin and PI3K signaling pathways, they remain largely unutilized in morphologic assessment coupled with ancillary diagnostic and prognostic modeling of ACC (Mete et al., 2018). Therefore, the major molecular mechanisms in the pathogenesis and progression are poorly understood. In 2003 and 2009, Giordano et al. performed unsupervised cluster analyses of transcriptome data to identify subgroups with different prognoses (Giordano et al., 2003; Giordano et al., 2009). These two studies laid the foundation for the molecular classification and prognostication of adrenocortical tumors and also provided a rich source of potential diagnostic and prognostic markers. Still, most cases of ACC were initially diagnosed with highly aggressive progression but were not candidates for curative therapies. Hence, potential biomarkers for diagnosis and treatment with high efficiency are urgently demanded.
Currently, microarray technology enables comprehensive mRNA expression profiling in ACC and can identify and investigate new biomarkers involved in tumorigenesis. A total of 150 DEGs and 24 hub genes were identified by microarray data analysis. GO and KEGG enrichment analysis showed association to the cell cycle, especially mitotic cycle checkpoint, mitotic spindle and oocyte meiosis, which was the most significant annotated function. Furthermore, among the 24 hub genes, the most significant molecular prognostic model integrated ZWINT, PRC1, CDKN3, CDK1 and CCNA2. Importantly, after reintegrating the weight of each gene, the new score was statistically the most significant parameter in both univariate and multivariate regression analysis. The gene set enrichment analysis (GSEA) method was used to visualize the significant signaling pathway analysis of ZWINT, PRC1, CDKN3, CDK1 and CCNA2.
ZW10 interactor (ZWINT), an interactor with ZW10, plays a vital role in rectifying incorrect centromere-microtubule attachment and regulating the miototic spindle checkpoint (Starr et al., 2000). Increased expression of ZWINT correlates with poor outcomes in human malignancies, including prostate, ovarian, bladder and lung cancers (Bhattacharjee et al., 2001; Endoh et al., 2004; Urbanucci et al., 2012; Xu et al., 2016). These new findings encourage further investigation of the potential clinical
A
ZWINT
B
PRC1
CESC LUSC
CESC
CHOL
CHOL
UCEC
LUSC
CESC
UCFC LUSC
CESC
CHOL
LIHC
LUSC
LUAD
LOAD
LUAD
SARC
LIHC
SARC
LHC
LUAD
GOM
CENA
LUAD
BRCA
BRCA
BRCA
BRCA
BLCA SARC
BLCA
PCPG
BEBE
ULEC
UCFC
ESCA
SARC
ESCA
GBM
ESCA
PCPG
ESCA
GBM
STAD
PCPG
STAD
STAD
HNSC
STAD
HNSC
BLCA
KIRC
BLCA
PRAD
PRAD
KICH
KIRC
KICH
PAAD
COAD
PAAD
COAD HNSC
COAD
COMO
HINSC
NIRP READ
CORD
-
KIRE
READ
DAAR
PARAD READ
PAAD
PARC
THCA
AIRE
SKCM THYM
THICA PRAD
BOAS
SKCM
THCA
-
PRAD
THICA
KICH
THYM SKCM
KICH THYM TGCT
SKCM
DLBC
DLBC
THYM
TGCT
UCs
TGCT UCS
TGCT UCS
ON
UCS
LAML
OV
LAML
LAML
OV
ov
ACC
LAML
MESO
ACE
DLRC
ACC
PLEC
DLBC
MESO
MESO
MEDO
UVM
MESO UVM
ACC
Toa
LGG
LGG
LGG
UVM
UVM
2000
1500
1000
500
0
500
1000
1500
2000
2500
2000
1500
1000
500
0
500
1000
1500
2000
2500
C
CDKN1
D
CDK1
LIHC
-
LIHC
GBM
CESC
CHOL LUSC
CESC
CESC
-
GBM
CHOL CHOL
CESC
LUSC LIHC
1
CHOL
LUAD
PCPG
LUAD LUSC
LIHC
PCPG SARC
LUSC SARC
UCEC
BRCA
UCEC
BRCA
ESCA
BRCA UCEC
BRÇA UCEC LUAD ESCA PCPG
ESCA
LUAD
GBM
KICH
GBM
KICH
ESCA
PCPG
SARC
KIRP
SARC MIED
BLCA
KIRC
STAD
BLCA
KIRC
COAD
STAD
BLCA
COAD
PICA
BICA
HINSC
BAAD
HNSC
STAD
PRAD
PAAD
STAD
PRAD
-
SIND SEAS
M
PRAD
PRAD
READ
PAAD
PARAD
DANS
KIRP
-
READ
READ
PHAD
READ
KIRC
MIDA
THYM
ICH
MIKE
-
THYM
THCA
MES
THICA
THCA
SKCM DLBC
SKCM
THICA SKCM
SKCM THYM
DLBC UCS
THYM
UCS
UCS
OV
UCS
OV
DLBC OV
DLBC
TGCT
ACC
TGCT TGCT
OV
TGCT
ACC
ACC
MESO
MESO LAML
ACC
MESO
LAML
LAML
UVM
LAML
UVM
LGG
MESO LGG
LGG
LGG
UVM
UVM
600
500
400
300
200
100
0
100
200
300
400
500
600
2000
1500
1000
500
0
500
1000
1500
2000
E
CCNA2
F
Median Rank
p-Value
Gene
1. P-value = 0.001
GBM
GBM
101.0
5.86E-4
ZWINT
CESC
CESC
2. P-value = 9.02E-12
CHO CHOL
1
2
CHOL
SABE
SARC
Median Rank
p-Value
Gene
LUAD
SARL
Tumor tissue
LUAD
1. P-value = null
PA
ABC
-
-
PRC1
BRCA
2. P-value = 5.97E-5
ESCA
BRCA
1
2
KICH
ESCA
KIRP KIRC
KICH
KIRP KIRC
Normal tissue
Median Rank
p-Value
Gene
BLCA
BLCA
9.0
4.45E-8
CDKN3
1. P-value = 8.89E-8
PCPG
STAD
PCPG
PRAD
STAD
HINSC
COAD PAAD READ
PRAD HNSC
1
2
2. P-value = 2.53E-14
COAD PAAD READ
Median Rank
p-Value
Gene
THCA
-
THICA
103.0
0.001
CDK1
1. P-value = 0.002
SKEM
2. P-value = 1.67E-13
TGCT
SKEM
TGCT
1
2
DLEC
DLBC
0
Median Rank
Gene
UCS
Uce UCS
p-Value
76.0
2.83E-4
CCNA2
1. P-value = 5.65E-4
LAML
LAML
ACC
ACC
2. P-value = 3.62E-12
MESO
LGG
MESO
1
2
UVM
LGG
UVM
1500
1200
900
600
300
0
300
600
900
1200
1500
1 5 10 25
25 10 5
1
M
L Not measured
%
A
A
New Tumor Event after First Treatment
B
Necrosis
…
.
-
.
12-
ZWINT
ZWINT
PRC1
12-
3
PRC1
10-
CDKN3
10-
CDKN3
8
CDK1
CCNA2
8
- CDK1
CCNA2
6-
6
4-
4
2
2
c
Weiss Score
D
Mitotic Rate
.
- ZWINT
.
.
**
*
*
12-
12-
- ZWINT
PRC1
- PRC1
10-
CDKN3
10
CDKN3
8-
CDK1
CCNA2
8
- CDK1
-
CCNA2
6-
6
4-
4
2
2
significance in human malignancies, yet the prognostic value of ZWINT in ACC has rarely been reported.
Protein Regulator of cytokinesis 1 (PRC1) protein is located in the nucleus. It is highly expressed in S and G2/M phases and shows an obvious drop in the G1 phase of the cell cycle (Freedland et al., 2013). During anaphase, it dynamically locates with the mitotic spindle and localizes to the cell midbody (Wu et al., 2018). Increasing evidence suggests that PRC1 may be involved in a cancer-specific manner, because of its negative correlation with p53 and overexpression in p53-defective cells in vitro (Li et al., 2004). In addition, Chen et al. and Zhan et al. demonstrated that PRC1 contributes to tumorigenesis by regulating the Wnt/ß-catenin signaling pathway in a positive feedback loop (Chen et al., 2016; Zhan et al., 2017), in which carcinogenesis and progression may feasibly be mediated in ACC.
| Covariates | Univariate analysis | Multivariate analysis | ||
|---|---|---|---|---|
| HR (95%CI) | P value | HR (95%CI) | P value | |
| Age at diagnosis (≤57 years vs. >58 years) | 1.703 (0.819 - 3.541) | 0.154 | ||
| Gender (male vs. female) | 0.977 (0.477 - 2.002) | 0.950 | ||
| Laterality (left vs. right) | 0.771 (0.381 - 1.563) | 0.471 | ||
| T stage (T1-T2 vs. T3-T4) | 3.846 (1.841 - 8.034) | <0.001 | ||
| N stage (N0 vs. N1) | 2.151 (0.820 - 5.641) | 0.119 | ||
| M stage (M0 vs. M1) | 3.104 (1.471 - 6.546) | 0.003 | 2.193 (0.977 - 4.921) | 0.057 |
| Pathologic stage (I - II vs. III - IV) | 3.937 (1.853 - 8.364) | <0.001 | ||
| Mitotic rate (≤5/50 HPF vs. >5/50 HPF) | 2.851 (1.164 - 6.984) | 0.022 | ||
| Invasion of tumor capsule (Present vs. Absent) | 2.074 (0.965 - 4.455) | 0.062 | ||
| Sinusoid invasion (Present vs. Absent) | 1.516 (0.678 - 3.389) | 0.311 | ||
| Necrosis (Present vs. Absent) | 6.501 (1.542 - 27.404) | 0.011 | ||
| Weiss score (≤3 vs. >3) | 2.816 (1.303 - 6.085) | 0.008 | ||
| New tumor event (Present vs. Absent) | 16.642 (5.673 - 48.822) | <0.001 | 9.041 (2.983 - 27.234) | <0.001 |
| Integrated expression score (High vs. Low) | 7.819 (3.569 - 17.114) | <0.001 | 2.767 (1.185 - 6.460) | 0.019 |
DFS, disease-free survival; ACC, clear cell renal cell carcinoma; HR, hazard ratio; CI, confidence interval; HPF, high power field.
*Multivariate Cox regression analyses of DFS in 76 enrolled ACC patients was run in “Back-LR” method. Statistically significant is considered as P value less than 0.05, indicated in bold.
| TABLE 5 | Univariate and multivariate Cox regression analyses of OS in 76 enrolled ACC patients. | ||||
|---|---|---|---|---|
| Covariates | Univariate analysis | Multivariate analysis* | ||
| HR (95%CI) | P value | HR (95%CI) | P value | |
| Age at diagnosis (≤57 years vs. >58 years) | 1.957 (0.913 - 4.196) | 0.085 | 4.959 (1.744 - 14.098) | 0.003 |
| Gender (male vs. female) | 0.996 (0.466 - 2.127) | 0.991 | ||
| Laterality (left vs. right) | 1.262 (0.591 - 2.695) | 0.548 | ||
| T stage (T1-T2 vs. T3-T4) | 10.693 (4.276 - 28.110) | <0.001 | ||
| N stage (N0 vs. N1) | 0.451 (0.171 - 1.191) | 0.108 | ||
| M stage (M0 vs. M1) | 7.340 (3.300 - 16.327) | <0.001 | 3.045 (1.094 - 8.477) | 0.033 |
| Pathologic stage (I - II vs. III - IV) | 7.157 (3.023 - 16.941) | <0.001 | ||
| Mitotic rate (≤5/50 HPF vs. >5/50 HPF) | 1.708 (0.743 - 3.926) | 0.208 | ||
| Invasion of tumor capsule (Present vs. Absent) | 3.015 (1.257 - 7.231) | 0.013 | ||
| Sinusoid invasion (Present vs. Absent) | 3.069 (1.297 - 7.262) | 0.011 | ||
| Necrosis (Present vs. Absent) | 5.135 (1.214 - 21.725) | 0.026 | ||
| Weiss score (≤3 vs. >3) | 3.467 (1.136 - 10.577) | 0.307 | ||
| New tumor event (Present vs. Absent) | 5.833 (2.351 - 14.473) | <0.001 | 3.609 (1.305 - 9.980) | 0.013 |
| Integrated expression score (High vs. Low) | 18.892 (6.830 - 52.255) | <0.001 | 18.719 (5.122 - 68.415) | <0.001 |
OS, overall survival; ACC, adrenocortical carcinoma; HR, hazard ratio; CI, confidence interval; HPF, high power field.
*Multivariate cox regression analyses of DFS in 76 enrolled ACC patients was run in “Back-LR” method. Statistically significant is considered as P value less than 0.05, indicated in bold.
Cyclin-dependent kinase inhibitor 3 (CDKN3) is part of the dual-specificity protein phosphatase family that dephosphorylates CDK2/CDK1 kinase and other cytokines (Hannon et al., 1994). Interestingly, a relationship between elevated CDKN3 expression and poor prognosis has been reported in many cancers by modulation of the cell cycle, mitotic spindle or p53 pathways (Berumen et al., 2014; Fan et al., 2015). A previous study has distinguished five genes modeling ACC using TOP2A, NDC80, CEP55, CDKN3 and CDK1, which may be utilized to form a board of progressive and predictive biomarkers for ACC for clinical purpose (Xiao et al., 2018). Thus, it is inferred that CDKN3 may be an oncogene in human ACC.
Cyclin dependent kinase 1 (CDK1) is a catalytic subunit of a highly conserved protein and is involved in many biological processes including cell cycle control, DNA damage repair, and checkpoint transcription (Skotheim et al., 2008; Enserink
and Kolodner, 2010). CDK1 plays an important regulatory role in the control of the eukaryotic cell cycle by modulating the centrosome cycle (Asghar et al., 2015). It has been previously reported that inhibition of CDK1 could serve as a therapeutic target via microRNA-7 for ACC samples in vivo (Glover et al., 2015). Meanwhile, CDC2, sharing approximately 63% amino- acid homology with CDK1, was found to be dysregulated in the cell cycle or retinoic acid signaling pathway by meta- analysis of genomic profiling data of adrenocortical tumors (Szabo et al., 2010).
Cyclin-A2 (CCNA2) belongs to a highly conserved cyclin family whose members function as regulators of the cell cycle. This protein interacts with CDK2 during G1/S and in G2/M phase, therefore promoting cell cycle transition (Pagano et al., 1992). There is accumulating evidence suggesting a role for CCNA2 in tumorigenesis of human malignancies. Kim et al. identified an SNP (rs769236) at the CCNA2 promoter
A
Disease Free Survival
C
Percent survival (%)
100
HIGH
100
LOW
50
P < 0.0001
80
Sensitivity
60
0
0
1000
2000
3000
4000
5000
Days
B
40
Overall Survival
Percent survival (%)
100
HIGH
20
LOW
AUC = 0.861
P < 0.0001
50
P < 0.0001
0
0
20
40
60
80
100
0
Specificity
0
1000
2000
3000
4000
5000
Days
that may be significantly associated with an increased risk of colon, liver and lung cancers (Kim et al., 2011). In addition, a significant delay in liver tumor formation was observed in mice with CCNA2-deficient hepatocytes (Gopinathan et al., 2014). As well as a prognostic value for CDK1 in ACC (Xiao et al., 2018), the mitotic checkpoint regulator CCNA2 may combine with other cell-cycle coding genes and be involved in aberrant regulation of the cell cycle network. It has not been evaluated whether this could be an effective approach to ACC treatment.
Our study represents the first attempt to construct a gene regulatory network incorporating DEGs and functional annotation of hub genes in ACC. An additional strength is that the alteration of ZWINT, PRC1, CDKN3, CDK1 and CCNA2 is significantly associated with worse OS and DFS, indicating that these genes may play important roles in the aggressive malignant phenotypes of ACC. At the same time, several limitations of this study are as follows. First, the data utilized in the study consisted of unbalanced ACC and normal control samples, which were restricted in quantity and downloaded from the GEO database, not generated by new
A
Enrichment plot: HALLMARK_MITOTIC_SPINDLE
Enrichment plot: HALLMARK_G2M_CHECKPOINT
C
Enrichment plot: HALLMARK_E2F_TARGETS
0.6
0.8
0.8
Enrichment score (ES)
0.5
ZWINT
Enrichment score (ES)
0.7
ZWINT
0.7
ZWINT
0.6
Enrichment score (ES)
0.6
0.4
0.5
0.5
0.3
0.4
0.4
0.2
0.3
03
0.2
0.2
0.1
0.1
0.1
0.0
0.0
0.0
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
1.0
HIGH (positively comelated)
1.0
HIGH (positively comrelated)
1.0
HIGH (positively conelated)
0.5
0.5
0.5
Zero cross at 9819
Zero cress at 9819
Zero cross at 9819
0.0
0.0
0.0
0.5
‘LOW (negatively comelated)
0.5
“LOW (negatively comelated)
0.5
“LOW (negatively cetrelated)
0
2,500
5.000
7,500
10,000
12,500
15,000
17,500
20,000
0
2,500
5,000
7,500
10,000
12.500
15,000
17,500
20,000
0
2,500
5,000
7,500
10,000
12,500
15,000
17,500
20,000
Rank in Ordered Dataset
Rank in Ordered Dataset
Rank in Ordered Dataset
Enrichment profile
Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
Enrichment profile -Hits
Ranking metric scores
D
Enrichment plot: HALLMARK_MITOTIC_SPINDLE
E
Enrichment plot: HALLMARK_G2M_CHECKPOINT
F
Enrichment plot: HALLMARK_E2F_TARGETS
0.40
0.5
Enrichment score (ES)
0.35
PRC1
PRC1
0.5
PRC1
0.30
Enrichment score (ES)
0,4
Enrichment score (ES)
0.25
0.4
0.20
0.3
0.15
0.3
0.10
0.2
0.2
0.05
0.1
0.00
0.1
0.05
0.0
0.0
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
1.0
HIGH (positively conelated)
1.0
‘HIGH (positively comrelated)
1.0
HIGH (positively correlated)
0.5
0.5
0.5
Zero cross at 9491
Zero cross at 9-491
Zero cross at 9491
0.0
0.0
0.0
0.5
‘LOW (negatively comelated)
0.5
“LOW (negatively comelated)
0.5
“LOW (negatively correlated)
0
2,500
5,000
7,500
10,000
12,500
15,000
17,500
0
2,500
5,000
7,500
10,000
12,500
15,000
17,500
0
2,500
5,000
7,500
10,000
12,500
15,000
17,500
Rank in Ordered Dataset
Rank in Ordered Dataset
Rank in Ordered Dataset
Enrichment profile - Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
G
Enrichment plot: HALLMARK_MITOTIC_SPINDLE
H
Enrichment plot: HALLMARK_G2M_CHECKPOINT
I
Enrichment plot: HALLMARK_E2F_TARGETS
0.6
0.8
0.0
0.5
CDKN3
0.7
Enrichment score (ES)
Enrichment score (ES)
CDKN3
Enrichment score (ES)
0.7
CDKN3
0.6
0.6
0.4
0.5
0.5
0.3
0.4
0.4
0.2
0.3
D
0.2
0.2
0.1
0.1
0.1
0.0
0.0
0.0
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
HIGH (positively comelated)
HIGIF (positively comrelated)
HIGH (positively conelated)
1.0
1.0
1.0
0.5
Zero cross at 0-481
O.
Zero cross at 0481
0.5
Zero cross at 9-481
0.0
0.0
0.0
0.5
‘LOW (negatively comelated)
0.5
‘LOW (negatively comelated)
0.5
LOW (negatively correlated)
0
2,500
5,000
7,500
10,000
12,500
15,000
17,500
20,000
0
2.500
5,000
7,500
10,000
12.500
15,000
17,500
20,000
0
2.500
5,000
7,500
10,000
12,500
15,000
17,500
20,000
Rank in Ordered Dataset
Rank in Ordered Dataset
Rank in Ordered Dataset
FIGURE 6 | Continued
Enrichment profile -Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
Enrichment profile -Hits
Ranking metric scores
J
Enrichment plot: HALLMARK_MITOTIC_SPINDLE
K
Enrichment plot: HALLMARK_G2M_CHECKPOINT
L
Enrichment plot: HALLMARK_E2F_TARGETS
0.6
0.0
0.8
Enrichment score (ES)
CDK1
Enrichment score (ES)
0.7
CDK1
0.7
CDK1
0.5
Enrichment score (ES)
0.6
0.6
0.4
0.5
0.5
0.3
0.4
0.4
0.2
0.3
0.3
0.1
0.2
0.2
0.1
0.1
0.0
0.0
0.0
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
1.0
‘HIGH’ (positively comelated)
1.0
HIGH (positively courelated)
1.0
HIGH (positively comrelated)
0.5
0.5
0.5
Zero cross at 9094
Zero cross at 9694
Zero cross at 9094
0.0
0.0
0.0
0.5
‘LOW (negatively correlated)
0.5
“LOW (negatively comelated)
0.5
‘LOW (negatively comrelated)
0
2,600
5,000
7,500
10,000
12,500
15,000
17,500
20,000
0
2,500
5,000
7,500
10,000
12,500
15,000
17,500
20,000
0
2,500
6,000
7,500
10,000
12,500
15,000
17,500
20,000
Rank in Ordered Dataset
Rank in Ordered Dataset
Rank in Ordered Dataset
Enrichment profile
Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
M
Enrichment plot: HALLMARK_MITOTIC_SPINDLE
N
Enrichment plot: HALLMARK_G2M_CHECKPOINT
Enrichment plot: HALLMARK_E2F_TARGETS
0.6
0.8
0.8
Enrichment score (ES)
CCNA2
0.7
CCNA2
0.7
CCNA2
0.5
Enrichment score (ES)
0.6
Enrichment score (ES)
0.6
0.4
0.5
0.5
0.3
0.4
0.4
0.2
0.3
0.3
0.2
0.2
0.1
0.1
0.1
0.0
0.0
0.0
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
1.0
‘HIGH (positively comelated)
1.0
‘HIGH’ (positively correlated)
1.0
HIGH (positively comelated)
0.5
0.6
0.5
Zero cross at 9804
Zero cross at 9804
Zero cross at 9804
0.0
0.0
0.0
0.5
‘LOW (negatively correlated)
0.5
“LOW (negatively comelated)
0.5
“LOW (negatively comrelated)
0
2,500
5.000
7.500
10,000
12,500
15,000
17,500
20,000
0
2,500
5.000
7.500
10,000
12,500
15,000
17,500
20,000
0
2,500
5,000
7.500
10,000
12,500
15,000
17,500
20,000
Rank in Ordered Dataset
Rank in Ordered Dataset
Rank in Ordered Dataset
Enrichment profile -Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
Enrichment profile - Hits
Ranking metric scores
DNA microarrays. Second, the microarray data contained relatively few ACC samples in the public database, and only 76 patients were enrolled from the TCGA cohort with corresponding transcriptome data. Third, prospective cohort was not used in this study. In addition, only the mRNA levels of hub genes are shown in this study, thus further functional works and validated cohorts are needed to verify these findings.
CONCLUSION
In conclusion, the present study identifies DEGs and hub genes that may be involved in poor prognosis and recurrence of ACC in silico. The transcriptional profiles of ZWINT, PRC1, CDKN3, CDK1 and CCNA2 are of prognostic value, and may assist in better understanding the underlying carcinogenesis or progression of ACC. Further studies are required to elucidate the molecular pathogenesis and alterations in signaling pathways of these genes in ACC.
ETHICS STATEMEMT
The ethics approval and consent to participate of the current study was approved and consented by the ethics committee of Fudan University Shanghai Cancer Center.
AUTHOR CONTRIBUTIONS
The work presented here was carried out in collaboration among all authors. D-WY, H-LZ, and Y-YQ defined the research theme, discussed analyses, interpretation and presentation. W-HX and JlW drafted the manuscript, analyzed the data, developed the algorithm and interpreted the results. JW co-worked on associated data collection, cohort validation and helped to draft the manuscript. F-NW, H-KW, and D-LC helped to perform the statistical analysis and reference collection. All authors read and approved the final manuscript.
FUNDING
This work is supported by Grants from the National Natural Science Foundation of China (No. 81202004, 81802525), Natural Science Foundation of Shanghai (No. 16ZR1406400), and Shanghai Sailing Program of China (No. 17YF1402700).
ACKNOWLEDGMENTS
We thank Catherine Perfect, MA (Cantab), from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.
SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00821/ full#supplementary-material
SUPPLEMENTARY FIGURE 1 | Functional and pathway enrichment analyses were performed using DAVID in bubble chart. (A) Changes in cellular components of DEGs were mainly enriched in the chromosome, centromeric region, extracellular region, actin cytoskeleton and mitotic spindle. (B) Changes in molecular funtions were mostly enriched in growth factor binding, kinase activity, extracellular matrix structure constituent and insulin-like grouth factor binding. (C) GO analysis results showed that changes in biological processes of DEGs were significantly enriched in mitotic cell cycle, cell cycle process, movement of cell or subcellular component and cell locomotion activity.
SUPPLEMENTARY FIGURE 2 | A network of the 24 genes and their co-expression genes was visualized and displayed in detail. (A) The biologic process and KEGG enrichment analysis of the hub genes were shown in different color. (B) The detailed functional notes and classification pie charts are listed as follows. 66.67% terms belong to mitotic cell cycle checkpoint, 15.79% to mitotic spindle organization, 12.28% to anaphase-promoting complex-dependent
catabolic process, 3.51% to protein localization to kinetochore and 1.75% to chromosome consederation.
SUPPLEMENTARY FIGURE 3 | Univariate survival analysis of the hub genes was performed using Kaplan-Meier curve. Besides MND1, each elevated expression in 24 hub gene showed markedly significant worse OS and DFS in ACC samples (P < 0.05).
SUPPLEMENTARY FIGURE 4 | Sensitivity analyze of GSE19750 and GSE90713 with Venn diagram, PPI network and the most significant module of DEGs. (A) DEGs were selected with a fold change >2 and P-value <0.01 among the mRNA expression profiling chip datasets GSE19750 and GSE90713. The 2 datasets showed an overlap of 315 genes in Venn diagram. (B) The PPI network of DEGs was constructed using Cytoscape. (C) The most significant module was obtained from PPI network with 31 nodes including ZWINT, PRC1, CDKN3, CDK1, CCNA2. Significant edges are marked in light blue with a K-score >0.800.
SUPPLEMENTARY FIGURE 5 | Univariate survival analysis of hub genes from sensitivity validated datasets was performed using Kaplan-Meier curve. Eight different DEGs are found different from these in three-chipset study, including UBE2C, GINS1, SMC2, ATAD2, PRIM1, NUSAP1, SPAG5, SMC4. Kaplan-Meier method was used to analysis mRNA expression level of 8 hub genes in TCGA cohort, which also showed significant correlation between elevated expression and progressive progression or poor prognosis (P < 0.05).
REFERENCES
Asghar U., Witkiewicz A. K., Turner, N. C., and Knudsen, E. S. (2015). The history and future of targeting cyclin-dependent kinases in cancer therapy. Nat. Rev. Drug Discov. 14 (2), 130-146. doi: 10.1038/nrd4504
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene ontology: tool for the unification of biology. Gene Ontol. Consort. Nat. Genet. 25 (1), 25-29. doi: 10.1038/75556
Bandettini, W. P., Kellman, P., Mancini, C., Booker, O. J., Vasu, S., Leung, S. W., et al. (2012). MultiContrast Delayed Enhancement (MCODE) improves detection of subendocardial myocardial infarction by late gadolinium enhancement cardiovascular magnetic resonance: a clinical validation study. J. Cardiovasc. Magn. Reson. 14, 83. doi: 10.1186/1532-429X-14-83
Bapat, S. A., Krishnan, A., Ghanate, A. D., Kusumbe, A. P., and Kalra, R. S. (2010). Gene expression: protein interaction systems network modeling identifies transformation-associated molecules and pathways in ovarian cancer. Cancer Res. 70 (12), 4809-4819. doi: 10.1158/0008-5472.CAN-10-0447
Berthon, A., Sahut-Barnola, I., Lambert-Langlais, S., de Joussineau, C., Damon- Soubeyrand, C., Louiset, E., et al. (2010). Constitutive beta-catenin activation induces adrenal hyperplasia and promotes adrenal cancer development. Hum. Mol. Genet. 19 (8), 1561-1576. doi: 10.1093/hmg/ddq029
Berumen, J., Espinosa, A. M., and Medina, I. (2014). Targeting CDKN3 in cervical cancer. Expert Opin. Ther. Targets 18 (10), 1149-1162. doi: 10.1517/14728222.2014.941808
Bhattacharjee, A., Richards, W. G., Staunton, J., Li, C., Monti, S., Vasa, P., et al. (2001). Classification of human lung carcinomas by mRNA expression profiling reveals distinct adenocarcinoma subclasses. Proc. Natl. Acad. Sci. U. S. A. 98 (24), 13790-13795. doi: 10.1073/pnas.191502998
Bindea, G., Mlecnik, B., Hackl, H., Charoentong, P., Tosolini, M., Kirilovsky, A., et al. (2009). ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25 (8), 1091- 1093. doi: 10.1093/bioinformatics/btp101
Bindea, G., Galon, J., and Mlecnik, B. (2013). CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics 29 (5), 661-663. doi: 10.1093/bioinformatics/btt019
Cerami, E., Gao, J., Dogrusoz, U., Gross, B. E., Sumer, S. O., Aksoy, B. A., et al. (2012). The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2 (5), 401-404. doi: 10.1158/2159-8290.CD-12-0095
Chen, J., Rajasekaran, M., Xia, H., Zhang, X., Kong, S. N., Sekar, K., et al. (2016). The microtubule-associated protein PRC1 promotes early recurrence of
hepatocellular carcinoma in association with the Wnt/beta-catenin signalling pathway. Gut 65 (9), 1522-1534. doi: 10.1136/gutjnl-2015-310625
de Fraipont, F., El Atifi, M., Cherradi, N., Le Moigne, G., Defaye, G., Houlgatte, R., et al. (2005). Gene expression profiling of human adrenocortical tumors using complementary deoxyribonucleic Acid microarrays identifies several candidate genes as markers of malignancy. J. Clin. Endocrinol. Metab. 90 (3), 1819-1829. doi: 10.1210/jc.2004-1075
Edgar, R., Domrachev, M., and Lash, A. E. (2002). Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 30 (1), 207-210. doi: 10.1093/nar/30.1.207
Else, T., Williams, A. R., Sabolch, A., Jolly, S., Miller, B. S., and Hammer, G. D. (2014). Adjuvant therapies and patient and tumor characteristics associated with survival of adult patients with adrenocortical carcinoma. J. Clin. Endocrinol. Metab. 99 (2), 455-461. doi: 10.1210/jc.2013-2856
Endoh, H., Tomida, S., Yatabe, Y., Konishi, H., Osada, H., Tajima, K., et al. (2004). Prognostic model of pulmonary adenocarcinoma by expression profiling of eight genes as determined by quantitative real-time reverse transcriptase polymerase chain reaction. J. Clin. Oncol. 22 (5), 811-819. doi: 10.1200/ JCO.2004.04.109
Enserink, J. M., and Kolodner, R. D. (2010). An overview of Cdk1-controlled targets and processes. Cell Div. 5, 11. doi: 10.1186/1747-1028-5-11
Fan, C., Chen, L., Huang, Q., Shen, T., Welsh, E. A., Teer, J. K., et al. (2015). Overexpression of major CDKN3 transcripts is associated with poor survival in lung adenocarcinoma. Br. J. Cancer 113 (12), 1735-1743. doi: 10.1038/ bjc.2015.378
Fassnacht, M., Kroiss, M., and Allolio, B. (2013). Update in adrenocortical carcinoma. J. Clin. Endocrinol. Metab. 98 (12), 4551-4564. doi: 10.1210/ jc.2013-3020
Franceschini, A., Szklarczyk, D., Frankild, S., Kuhn, M., Simonovic, M., Roth, A., et al. (2013). STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 41 (Database issue), D808-D815. doi: 10.1093/nar/gks1094
Freedland, S. J., Gerber, L., Reid, J., Welbourn, W., Tikishvili, E., Park, J., et al. (2013). Prognostic utility of cell cycle progression score in men with prostate cancer after primary external beam radiation therapy. Int. J. Radiat. Oncol. Biol. Phys. 86 (5), 848-853. doi: 10.1016/j.ijrobp.2013.04.043
Gaujoux, S., Grabar, S., Fassnacht, M., Ragazzon, B., Launay, P., Libé, R., et al. (2011). ß-catenin activation is associated with specific clinical and pathologic characteristics and a poor outcome in adrenocortical carcinoma. Clin. Cancer Res. 17 (2), 328-336. doi: 10.1158/1078-0432.CCR-10-2006
Gicquel, C., Bertagna, X., Gaston, V., Coste, J., Louvel, A., Baudin, E., et al. (2001). Molecular markers and long-term recurrences in a large cohort of patients with sporadic adrenocortical tumors. Cancer Res. 61 (18), 6762-6767.
Giordano, T. J., Thomas, D. G., Kuick, R., Lizyness, M., Misek, D. E., Smith, A. L., et al. (2003). Distinct transcriptional profiles of adrenocortical tumors uncovered by DNA microarray analysis. Am. J. Pathol. 162 (2), 521-531. doi: 10.1016/S0002-9440(10)63846-1
Giordano, T. J., Kuick, R., Else, T., Gauger, P. G., Vinco, M., Bauersfeld, J., et al. (2009). Molecular classification and prognostication of adrenocortical tumors by transcriptome profiling. Clin. Cancer Res. 15 (2), 668-676. doi: 10.1158/1078-0432.CCR-08-1067
Glover, A. R., Zhao, J. T., Gill, A. J., Weiss, J., Mugridge, N., Kim, E., et al. (2015). MicroRNA-7 as a tumor suppressor and novel therapeutic for adrenocortical carcinoma. Oncotarget 6 (34), 36675-36688. doi: 10.18632/oncotarget.5383
Gopinathan, L., Tan, S. L., Padmakumar, V. C., Coppola, V., Tessarollo, L., and Kaldis, P. (2014). Loss of Cdk2 and cyclin A2 impairs cell proliferation and tumorigenesis. Cancer Res. 74 (14), 3870-3879. doi: 10.1158/0008-5472. CAN-13-3440
Hannon, G. J., Casso, D., and Beach, D. (1994). KAP: a dual specificity phosphatase that interacts with cyclin-dependent kinases. Proc. Natl. Acad. Sci. U. S. A. 91 (5), 1731-1735. doi: 10.1073/pnas.91.5.1731
Heaton, J. H., Wood, M. A., Kim, A. C., Lima, L. O., Barlaskar, F. M., Almeida, M. Q., et al. (2012). Progression to adrenocortical tumorigenesis in mice and humans through insulin-like growth factor 2 and beta-catenin. Am. J. Pathol. 181 (3), 1017-1033. doi: 10.1016/j.ajpath.2012.05.026
Huang, D. W., Sherman, B. T., Tan, Q., Collins, J. R., Alvord, W. G., Roayaei, J., et al. (2007). The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 8 (9), R183. doi: 10.1186/gb-2007-8-9-r183
Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M. (2016). KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 04, 44. doi: 10.1093/nar/gkv1070
Kebebew, E., Reiff, E., Duh, Q. Y., Clark, O. H., and McMillan, A. (2006). Extent of disease at presentation and outcome for adrenocortical carcinoma: have we made progress? World J. Surg. 30 (5), 872-878. doi: 10.1007/s00268-005-0329-x
Kerkhofs, T. M., Verhoeven, R. H., Van der Zwan, J. M., Dieleman, J., Kerstens, M.N., Links, T. P., et al. (2013). Adrenocortical carcinoma: a population-based study on incidence and survival in the Netherlands since 1993. Eur. J. Cancer 49 (11), 2579-2586. doi: 10.1016/j.ejca.2013.02.034
Kim, D. H., Park, S. E., Kim, M., Ji, Y. I., Kang, M. Y., Jung, E. H., et al. (2011). A functional single nucleotide polymorphism at the promoter region of cyclin A2 is associated with increased risk of colon, liver, and lung cancers. Cancer 117 (17), 4080-4091. doi: 10.1002/cncr.25930
Li, C., Lin, M., and Liu, J. (2004). Identification of PRC1 as the p53 target gene uncovers a novel function of p53 in the regulation of cytokinesis. Oncogene 23 (58), 9336-9347. doi: 10.1038/sj.onc.1208114
Mete, O., Gucer, H., Kefeli, M., and Asa, S. L. (2018). Diagnostic and prognostic biomarkers of Adrenal cortical carcinoma. Am. J. Surg. Pathol. 42 (2), 201-213. doi: 10.1097/PAS.0000000000000943
Mohan, D. R., Lerario, A. M., and Hammer, G. D. (2018). Therapeutic targets for adrenocortical carcinoma in the Genomics era. J. Endocr. Soc. 2 (11), 1259- 1274. doi: 10.1210/js.2018-00197
Pagano, M., Pepperkok, R., Verde, F., Ansorge, W., and Draetta, G. (1992). Cyclin A is required at two points in the human cell cycle. EMBO J. 11 (3), 961-971. doi: 10.1002/j.1460-2075.1992.tb05135.x
Paragliola, R. M., Torino F., Papi, G., Locantore, P., Pontecorvi, A., and Corsello, S. M. (2018). Role of mitotane in Adrenocortical Carcinoma - Review and state of the art. Eur. Endocrinol. 14 (2), 62-66. doi: 10.17925/EE.2018.14.2.62
Poli, G., Guasti, D., Rapizzi, E., Fucci, R., Canu, L., Bandini, A., et al. (2013). Morphofunctional effects of mitotane on mitochondria in human
adrenocortical cancer cells. Endocr. Relat. Cancer 20 (4), 537-550. doi: 10.1530/ ERC-13-0150
Sharan, R., Ulitsky, I., and Shamir, R. (2007). Network-based prediction of protein function. Mol. Syst. Biol. 3, 88. doi: 10.1038/msb4100129
Skotheim, J. M., Di Talia, S., Siggia, E. D., and Cross, F. R. (2008). Positive feedback of G1 cyclins ensures coherent cell cycle entry. Nature 454 (7202), 291-296. doi: 10.1038/nature07118
Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P. L., and Ideker, T. (2011). Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27 (3), 431-432. doi: 10.1093/bioinformatics/btq675
Starr, D. A., Saffery, R., Li, Z., Simpson, A. E., Choo, K. H., Yen, T. J., et al. (2000). HZwint-1, a novel human kinetochore component that interacts with HZW10. J. Cell Sci. 113 (Pt 11), 1939-1950.
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102 (43), 15545-15550. doi: 10.1073/pnas.0506580102
Szabó, P. M., Tamási, V., Molnár, V., Andrásfalvy, M., Tömböl, Z., Farkas, R., et al. (2010). Meta-analysis of adrenocortical tumour genomics data: novel pathogenic pathways revealed. Oncogene 29 (21), 3163-3172. doi: 10.1038/ onc.2010.80
Urbanucci, A., Sahu, B., Seppälä, J., Larjo, A., Latonen, L. M., Waltering, K. K., et al. (2012). Overexpression of androgen receptor enhances the binding of the receptor to the chromatin in prostate cancer. Oncogene 31 (17), 2153-2163. doi: 10.1038/onc.2011.401
Wajchenberg, B. L., Albergaria Pereira, M. A., Medonca, B. B., Latronico, A. C., Campos Carneiro, P., Alves, V. A., et al. (2000). Adrenocortical carcinoma: clinical and laboratory observations. Cancer 88 (4), 711-736. doi: 10.1002/ (SICI) 1097-0142(20000215)88:4<711:AID-CNCR1>3.0.CO;2-W
Wu, J., Vallenius, T., Ovaska, K., Westermarck, J., Mäkelä, T. P., and Hautaniemi, S. (2009). Integrated network analysis platform for protein-protein interactions. Nat. Methods 6 (1), 75-77. doi: 10.1038/nmeth.1282
Wu, F., Shi, X., Zhang, R., Tian, Y., Wang, X., Wei, C., et al. (2018). Regulation of proliferation and cell cycle by protein regulator of cytokinesis 1 in oral squamous cell carcinoma. Cell Death Dis. 9 (5), 564. doi: 10.1038/ s41419-018-0618-6
Xiao, H., Xu, D., Chen, P., Zeng, G., Wang, X., and Zhang, X. (2018). Identification of five genes as a potential biomarker for predicting progress and prognosis in adrenocortical carcinoma. J. Cancer 9 (23), 4484-4495. doi: 10.7150/jca.26698
Xu, Z., Zhou, Y., Cao, Y., Dinh, T. L., Wan, J., and Zhao, M. (2016). Identification of candidate biomarkers and analysis of prognostic values in ovarian cancer by integrated bioinformatics analysis. Med. Oncol. 33 (11), 130. doi: 10.1007/ s12032-016-0840-y
Zhan, P., Zhang, B., Xi, G. M., Wu, Y., Liu, H. B., Liu, Y. F., et al. (2017). PRC1 contributes to tumorigenesis of lung adenocarcinoma in association with the Wnt/beta-catenin signaling pathway. Mol. Cancer. 16 (1), 108 doi: 10.1186/ s12943-017-0682-z
Conflict of Interest Statement: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Copyright @ 2019 Xu, Wu, Wang, Wan, Wang, Cao, Qu, Zhang and Ye. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.