Screening and Identification of Potential Prognostic Biomarkers in Adrenocortical Carcinoma

Abstract

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.

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/β-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 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. Log₂FC (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.

Figure 1

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 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, 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).

Figure 2

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.

Figure 3

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 (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).

Figure 5

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 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.

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/β-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 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.

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 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 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.

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).

References

• 1

  AsgharU.WitkiewiczA. K.TurnerN. C.KnudsenE. 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

  • CrossRef
  • Google Scholar

• 2

  AshburnerM.BallC. A.BlakeJ. A.BotsteinD.ButlerH.CherryJ. 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

  • CrossRef
  • Google Scholar

• 3

  BandettiniW. P.KellmanP.ManciniC.BookerO. J.VasuS.LeungS. 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

  • CrossRef
  • Google Scholar

• 4

  BapatS. A.KrishnanA.GhanateA. D.KusumbeA. P.KalraR. 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

  • CrossRef
  • Google Scholar

• 5

  BerthonA.Sahut-BarnolaI.Lambert-LanglaisS.de JoussineauC.Damon-SoubeyrandC.LouisetE.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

  • CrossRef
  • Google Scholar

• 6

  BerumenJ.EspinosaA. M.MedinaI. (2014). Targeting CDKN3 in cervical cancer. Expert Opin. Ther. Targets18 (10), 1149–1162. doi: 10.1517/14728222.2014.941808

  • CrossRef
  • Google Scholar

• 7

  BhattacharjeeA.RichardsW. G.StauntonJ.LiC.MontiS.VasaP.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

  • CrossRef
  • Google Scholar

• 8

  BindeaG.MlecnikB.HacklH.CharoentongP.TosoliniM.KirilovskyA.et al. (2009). ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics25 (8), 1091–1093. doi: 10.1093/bioinformatics/btp101

  • CrossRef
  • Google Scholar

• 9

  BindeaG.GalonJ.MlecnikB. (2013). CluePedia Cytoscape plugin: pathway insights using integrated experimental and in silico data. Bioinformatics29 (5), 661–663. doi: 10.1093/bioinformatics/btt019

  • CrossRef
  • Google Scholar

• 10

  CeramiE.GaoJ.DogrusozU.GrossB. E.SumerS. O.AksoyB. 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

  • CrossRef
  • Google Scholar

• 11

  ChenJ.RajasekaranM.XiaH.ZhangX.KongS. N.SekarK.et al. (2016). The microtubule-associated protein PRC1 promotes early recurrence of hepatocellular carcinoma in association with the Wnt/beta-catenin signalling pathway. Gut65 (9), 1522–1534. doi: 10.1136/gutjnl-2015-310625

  • CrossRef
  • Google Scholar

• 12

  de FraipontF.El AtifiM.CherradiN.Le MoigneG.DefayeG.HoulgatteR.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

  • CrossRef
  • Google Scholar

• 13

  EdgarR.DomrachevM.LashA. 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

  • CrossRef
  • Google Scholar

• 14

  ElseT.WilliamsA. R.SabolchA.JollyS.MillerB. S.HammerG. 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

  • CrossRef
  • Google Scholar

• 15

  EndohH.TomidaS.YatabeY.KonishiH.OsadaH.TajimaK.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

  • CrossRef
  • Google Scholar

• 16

  EnserinkJ. M.KolodnerR. D. (2010). An overview of Cdk1-controlled targets and processes. Cell Div.5, 11. doi: 10.1186/1747-1028-5-11

  • CrossRef
  • Google Scholar

• 17

  FanC.ChenL.HuangQ.ShenT.WelshE. A.TeerJ. K.et al. (2015). Overexpression of major CDKN3 transcripts is associated with poor survival in lung adenocarcinoma. Br. J. Cancer113 (12), 1735–1743. doi: 10.1038/bjc.2015.378

  • CrossRef
  • Google Scholar

• 18

  FassnachtM.KroissM.AllolioB. (2013). Update in adrenocortical carcinoma. J. Clin. Endocrinol. Metab.98 (12), 4551–4564. doi: 10.1210/jc.2013-3020

  • CrossRef
  • Google Scholar

• 19

  FranceschiniA.SzklarczykD.FrankildS.KuhnM.SimonovicM.RothA.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

  • CrossRef
  • Google Scholar

• 20

  FreedlandS. J.GerberL.ReidJ.WelbournW.TikishviliE.ParkJ.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

  • CrossRef
  • Google Scholar

• 21

  GaujouxS.GrabarS.FassnachtM.RagazzonB.LaunayP.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

  • CrossRef
  • Google Scholar

• 22

  GicquelC.BertagnaX.GastonV.CosteJ.LouvelA.BaudinE.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.

  • Google Scholar

• 23

  GiordanoT. J.ThomasD. G.KuickR.LizynessM.MisekD. E.SmithA. 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

  • CrossRef
  • Google Scholar

• 24

  GiordanoT. J.KuickR.ElseT.GaugerP. G.VincoM.BauersfeldJ.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

  • CrossRef
  • Google Scholar

• 25

  GloverA. R.ZhaoJ. T.GillA. J.WeissJ.MugridgeN.KimE.et al. (2015). MicroRNA-7 as a tumor suppressor and novel therapeutic for adrenocortical carcinoma. Oncotarget6 (34), 36675–36688. doi: 10.18632/oncotarget.5383

  • CrossRef
  • Google Scholar

• 26

  GopinathanL.TanS. L.PadmakumarV. C.CoppolaV.TessarolloL.KaldisP. (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

  • CrossRef
  • Google Scholar

• 27

  HannonG. J.CassoD.BeachD. (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

  • CrossRef
  • Google Scholar

• 28

  HeatonJ. H.WoodM. A.KimA. C.LimaL. O.BarlaskarF. M.AlmeidaM. 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

  • CrossRef
  • Google Scholar

• 29

  HuangD. W.ShermanB. T.TanQ.CollinsJ. R.AlvordW. G.RoayaeiJ.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

  • CrossRef
  • Google Scholar

• 30

  KanehisaM.SatoY.KawashimaM.FurumichiM.TanabeM. (2016). KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res.04, 44. doi: 10.1093/nar/gkv1070

  • CrossRef
  • Google Scholar

• 31

  KebebewE.ReiffE.DuhQ. Y.ClarkO. H.McMillanA. (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

  • CrossRef
  • Google Scholar

• 32

  KerkhofsT. M.VerhoevenR. H.Van der ZwanJ. M.DielemanJ.KerstensM. N.LinksT. P.et al. (2013). Adrenocortical carcinoma: a population-based study on incidence and survival in the Netherlands since 1993. Eur. J. Cancer49 (11), 2579–2586. doi: 10.1016/j.ejca.2013.02.034

  • CrossRef
  • Google Scholar

• 33

  KimD. H.ParkS. E.KimM.JiY. I.KangM. Y.JungE. 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. Cancer117 (17), 4080–4091. doi: 10.1002/cncr.25930

  • CrossRef
  • Google Scholar

• 34

  LiC.LinM.LiuJ. (2004). Identification of PRC1 as the p53 target gene uncovers a novel function of p53 in the regulation of cytokinesis. Oncogene23 (58), 9336–9347. doi: 10.1038/sj.onc.1208114

  • CrossRef
  • Google Scholar

• 35

  MeteO.GucerH.KefeliM.AsaS. L. (2018). Diagnostic and prognostic biomarkers of Adrenal cortical carcinoma. Am. J. Surg. Pathol.42 (2), 201–213. doi: 10.1097/PAS.0000000000000943

  • CrossRef
  • Google Scholar

• 36

  MohanD. R.LerarioA. M.HammerG. D. (2018). Therapeutic targets for adrenocortical carcinoma in the Genomics era. J. Endocr. Soc.2 (11), 1259–1274. doi: 10.1210/js.2018-00197

  • CrossRef
  • Google Scholar

• 37

  PaganoM.PepperkokR.VerdeF.AnsorgeW.DraettaG. (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

  • CrossRef
  • Google Scholar

• 38

  ParagliolaR. M.TorinoF.PapiG.LocantoreP.PontecorviA.CorselloS. 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

  • CrossRef
  • Google Scholar

• 39

  PoliG.GuastiD.RapizziE.FucciR.CanuL.BandiniA.et al. (2013). Morphofunctional effects of mitotane on mitochondria in human adrenocortical cancer cells. Endocr. Relat. Cancer20 (4), 537–550. doi: 10.1530/ERC-13-0150

  • CrossRef
  • Google Scholar

• 40

  SharanR.UlitskyI.ShamirR. (2007). Network-based prediction of protein function. Mol. Syst. Biol.3, 88. doi: 10.1038/msb4100129

  • CrossRef
  • Google Scholar

• 41

  SkotheimJ. M.Di TaliaS.SiggiaE. D.CrossF. R. (2008). Positive feedback of G1 cyclins ensures coherent cell cycle entry. Nature454 (7202), 291–296. doi: 10.1038/nature07118

  • CrossRef
  • Google Scholar

• 42

  SmootM. E.OnoK.RuscheinskiJ.WangP. L.IdekerT. (2011). Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics27 (3), 431–432. doi: 10.1093/bioinformatics/btq675

  • CrossRef
  • Google Scholar

• 43

  StarrD. A.SafferyR.LiZ.SimpsonA. E.ChooK. H.YenT. J.et al. (2000). HZwint-1, a novel human kinetochore component that interacts with HZW10. J. Cell Sci.113 (Pt 11), 1939–1950.

  • Google Scholar

• 44

  SubramanianA.TamayoP.MoothaV. K.MukherjeeS.EbertB. L.GilletteM. 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

  • CrossRef
  • Google Scholar

• 45

  SzabóP. M.TamásiV.MolnárV.AndrásfalvyM.TömbölZ.FarkasR.et al. (2010). Meta-analysis of adrenocortical tumour genomics data: novel pathogenic pathways revealed. Oncogene29 (21), 3163–3172. doi: 10.1038/onc.2010.80

  • CrossRef
  • Google Scholar

• 46

  UrbanucciA.SahuB.SeppäläJ.LarjoA.LatonenL. M.WalteringK. K.et al. (2012). Overexpression of androgen receptor enhances the binding of the receptor to the chromatin in prostate cancer. Oncogene31 (17), 2153–2163. doi: 10.1038/onc.2011.401

  • CrossRef
  • Google Scholar

• 47

  WajchenbergB. L.Albergaria PereiraM. A.MedoncaB. B.LatronicoA. C.Campos CarneiroP.AlvesV. A.et al. (2000). Adrenocortical carcinoma: clinical and laboratory observations. Cancer88 (4), 711–736. doi: 10.1002/(SICI)1097-0142(20000215)88:4<711::AID-CNCR1>3.0.CO;2-W

  • CrossRef
  • Google Scholar

• 48

  WuJ.ValleniusT.OvaskaK.WestermarckJ.MäkeläT. P.HautaniemiS. (2009). Integrated network analysis platform for protein-protein interactions. Nat. Methods6 (1), 75–77. doi: 10.1038/nmeth.1282

  • CrossRef
  • Google Scholar

• 49

  WuF.ShiX.ZhangR.TianY.WangX.WeiC.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

  • CrossRef
  • Google Scholar

• 50

  XiaoH.XuD.ChenP.ZengG.WangX.ZhangX. (2018). Identification of five genes as a potential biomarker for predicting progress and prognosis in adrenocortical carcinoma. J. Cancer9 (23), 4484–4495. doi: 10.7150/jca.26698

  • CrossRef
  • Google Scholar

• 51

  XuZ.ZhouY.CaoY.DinhT. L.WanJ.ZhaoM. (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

  • CrossRef
  • Google Scholar

• 52

  ZhanP.ZhangB.XiG. M.WuY.LiuH. B.LiuY. 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

  • CrossRef
  • Google Scholar