Screening and identification of key biomarkers in adrenocortical carcinoma based on bioinformatics analysis.

Adrenocortical carcinoma (ACC) is a rare malignancy with a poor prognosis. The presently available understanding of the pathogenesis of ACC is incomplete and the treatment options for patients with ACC are limited. Gene marker identification is required for accurate and timely diagnosis of the disease. In order to identify novel candidate genes associated with the occurrence and progression of ACC, the microarray datasets, GSE12368 and GSE19750, were obtained from Gene Expression Omnibus. Differentially expressed genes (DEGs) were identified, and functional enrichment analysis was performed. A protein-protein interaction network (PPI) was constructed to identify significantly altered modules, and module analysis was performed using Search Tool for the Retrieval of Interacting Genes and Cytoscape. A total of 228 DEGs were screened, consisting of 29 up and 199 downregulated genes. The enriched functions and pathways of the DEGs primarily included 'cell division', 'regulation of transcription involved in G1/S transition of mitotic cell cycle', 'G1/S transition of mitotic cell cycle', 'p53 signaling pathway' and 'oocyte meiosis'. A total of 14 hub genes were identified, and biological process analysis revealed that these genes were significantly enriched in cell division and mitotic cell cycle. Furthermore, survival analysis revealed that AURKA, TYMS, GINS1, RACGAP1, RRM2, EZH2, ZWINT, CDK1, CCNB1, NCAPG and TPX2 may be involved in the tumorigenesis, progression or prognosis of ACC. In conclusion, the 14 hub genes identified in the present study may aid researchers in elucidating the molecular mechanisms associated with the tumorigenesis and progression of ACC, and may be powerful and promising candidate biomarkers for the diagnosis and treatment of ACC. Copyright: © Xing et al.

Introduction

Adrenocortical carcinoma (ACC) is a rare endocrine malignancy that arises from the adrenal cortex, with an occurrence rate of 0.7–2.0 cases per 1,000,000 each year (1–3). The overall 5-year survival rate for patients in the Netherlands is ~32%, and patients with stage III and IV have a particularly lower survival rate, due to its highly aggressive biological behavior (2,4). ACC-associated mortality accounts for 0.02–0.20% of all cancer-associated mortalities (5,6). At present, complete surgical resection with or without mitotane adjuvant treatment is the only treatment option (7,8). However, 40–70% of patients with ACC present with metastasis at the time of diagnosis (9). Therefore, it is important to explore the molecular mechanisms underlying ACC and identify candidate biomarkers for its diagnosis and treatment.

At present, the diagnosis and classification of adrenocortical cancers relies on histological examination of tumor sections and immunohistochemical markers, such as Ki-67, IGF2 and SF-1, are used to support the diagnosis of ACC (10). The Weiss or modified Weiss score systems are most often used in diagnosis as the primary determinants of malignancy in adrenocortical tumors (10,11); however, the diagnosis of these tumors remains challenging, particularly for rare subtypes of ACC such as oncocytic, myxoid and sarcomatoid subtypes (12).

Following developments in microarray technology, several studies have demonstrated that abnormal expressed and mutated genes are involved in the tumorigenesis and progression of ACC (9,13–15). For example, using DNA microarray analysis, Giordano et al (13) demonstrated that several cell cycle and proliferation genes, such as Cyclin B2 (CCNB2), Abnormal Spindle Microtubule Assembly, Ribonucleotide Reductase Regulatory Subunit M2 (RRM2), DNA Topoisomerase II a and Cyclin Dependent Kinase Inhibitor 3, as well as genes associated with tumor invasion, such as Secreted Phosphoprotein 1 may serve as potential diagnostic biomarkers that could be developed into useful immunohistochemical tools. Kulshrestha and Suman (14) identified a total of 53 genes as common hubs of the disease system, which may exert important biological functions in pediatric adrenocortical adenoma and carcinoma. Cyclin Dependent Kinase 1 (CDK1), Cyclin B1 (CCNB1), Cell Division Cycle 20 and BUB1 Mitotic Checkpoint Serine/Threonine Kinase B may serve as potential biomarkers of pediatric ACC and as potential targets for its treatment (16). Yuan et al (15) analyzed 12 hub genes associated with the progression and prognosis of ACC by weighted gene co-expression network analysis. In addition, Duregon et al (17) assessed the expression of miRNAs associated with the regulation of the IGF2 gene and hypoxia induced microRNA in histological variants including 35 classical, 6 myxoid and 10 oncocytic cases of ACC and reported that miR-483-3p, miR-483-5p and miR-210, which were identified as candidates for tumor aggressiveness and poor prognosis in ACC, are differentially expressed in ACC variants. However, the latent molecular and pathway interactions of ACC have yet to be completely elucidated. There is a need for additional studies with comprehensive and integrated genomic characterization, combined with clinical data to explore the molecular mechanisms and identify candidate biomarkers for the diagnosis of ACC, as well as advance the presently available understanding of the tumorigenesis mechanism.

The aim of the present study was to analyze the differentially expressed genes (DEGs) in ACC by combining two mRNA microarray datasets from the Gene Expression Omnibus (GEO) database. Subsequently, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and protein-protein interaction (PPI) network analyses were performed to provide detailed insights into the biological mechanisms in ACC. GEO and Oncomine databases were subsequently combined to validate the importance of the hub genes. In conclusion, using bioinformatics methods, the present study identified 14 hub genes which provided significant diagnostic and prognostic value and may serve as candidate biomarkers for ACC.

Materials and methods

Data resources

The gene expression data was retrieved from the GEO database (ncbi.nlm.nih.gov/geo/). Two gene expression datasets, GSE12368 (18) and GSE19750 (19) were downloaded from GEO using GPL570 Affymetrix Human Genome U133 Plus 2.0 Array. The GSE12368 dataset contained 12 ACC and 6 normal samples. The GSE19750 dataset contained 44 ACC and 4 normal adrenal glands samples.

Identification of DEGs

The DEGs between ACC and normal adrenal gland samples were screened using GEO2R (ncbi.nlm.nih.gov/geo/geo2r). GEO2R is an interactive web tool using limma R packages (version R 3.2.3; limma 3.26.8) (20), which allows users to compare two or more datasets in a GEO series, to identify DEGs across experimental conditions (21,22). The adjusted P-values (adj. P) and Benjamini and Hochberg false discovery rate provided a balance between the discovery of statistically significant genes and the limitations of false-positive results. Probe sets without corresponding gene symbols or genes with >1 probe sets were removed or averaged, respectively. The cut-off criteria of |logFC(fold-change)|>1 and adj. P<0.05 were considered statistically significant.

KEGG and GO enrichment analysis of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov; version 6.8) is an online biological information database that integrates biological data and analysis tools, and provides a comprehensive set of functional annotation information of genes and proteins for users to extract biological information (23). KEGG (http://www.genome.ad.jp/kegg) is a knowledge base for systematic analysis of gene functions (24). GO analysis predicts the function of the target genes in three aspects, including biological processes (BPs), cellular components (CCs) and molecular functions (MFs) (25). To analyze the possible functions of DEGs, functional annotation was performed using the DAVID database; P<0.05 was considered statistically significant.

PPI network construction and module analysis

The PPI network was predicted using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org; version 10.5) database (26). Analyzing the functional interactions between proteins can predict the interaction relationship involved in the development and progression of ACC (27,28). In the present study, a PPI network of DEGs was constructed using the STRING database, and an interaction with a combined score of >0.4 was considered statistically significant. Cytoscape (version 3.6.1; http://cytoscape.org) is an open source bioinformatics platform for visualizing molecular interaction networks (29). The Molecular Complex Detection (MCODE; version 1.4.2) Cytoscape plugin allows for clustering a given network based on topology to identify densely connected regions (30). The PPI networks were constructed using Cytoscape and the most significant module in the PPI networks was identified using MCODE. The selection criteria were as follows: Degree cut-off=2, node score cut-off=0.2, max depth=100 and k-score=2. Subsequently, KEGG and GO analyses for genes in this module were performed using DAVID. Data analysis of biological processes in the hub genes was performed using R (http://www.r-project.org, version 3.2.4).

Hub gene selection and analysis

Hub genes with a degree of ≥10 were selected, and the network of the genes and their co-expression genes was analyzed using cBioPortal (http://www.cbioportal.org) (31,32). Hierarchical clustering of hub genes was constructed using UCSC Cancer Genomics Browser (http://genome-cancer.ucsc.edu) (33). To evaluate the prognostic value of the selected hub genes in ACC, overall and disease-free survival based on expression of the hub genes were performed using Kaplan-Meier curves in cBioPortal. In addition, the association between hub gene expression and tumor Weiss grade in patients with ACC was analyzed using the Oncomine database (www.oncomine.com). mRNA expression analyses of thymidylate synthetase (TYMS), GINS complex subunit 1 (GINS1), ribonucleotide reductase regulatory subunit M2 (RRM2), ZW10 interacting kinetochore protein (ZWINT) and structural maintenance of chromosomes 4 (SMC4) genes in ACC vs. normal tissues was performed in the Giordano Adrenal and Giordano Adrenal 2 datasets from the Oncomine database (13,34).

Results

Identification of DEGs in ACC

Following standardization of the microarray results by GEO2R, 970 and 998 DEGs were identified in the GSE12368 and GSE19750 gene expression profile datasets, respectively. The overlap between the 2 datasets contained 228 genes, as illustrated in the Venn diagram (Fig. 1A), consisting of 29 up and 199 downregulated genes in ACC tissue compared with normal tissues.

Figure 1.

Venn diagram, PPI network and the most significant module of DEGs. (A) DEGs were selected with a fold change >2 and an adjusted P<0.05 from GSE12368 and GSE1975 mRNA expression profiling sets. The 2 datasets revealed that there was an overlap of 228 genes. (B) The PPI network of DEGs was constructed using Cytoscape software and the 19 upregulated genes were marked in light red and the 94 downregulated genes are marked in light blue. (C) The most significant module was obtained from the PPI network which contained 14 nodes and 88 edges. PPI, protein-protein interaction network; DEGs, differentially expressed genes.

KEGG and GO enrichment analyses of DEGs

To investigate the functional annotation of the DEGs, GO terms and pathway enrichment analysis were performed using DAVID. The results indicated that the BPs of DEGs were significantly enriched in ‘Cell division’, ‘Regulation of transcription involved in G1/S transition of mitotic cell cycle’, ‘G1/S transition of mitotic cell cycle’, ‘aging’ and ‘signal transduction’ (Table I). Changes in MFs were primarily enriched in protein kinase binding, protein binding and calcium ion binding (Table I). Changes in the DEG CCs were primarily enriched in the nucleus and extracellular space (Table I). In addition, KEGG pathway analysis showed that the upregulated DEGs were mainly involved in the ‘p53 signaling pathway’, ‘Oocyte meiosis’ and ‘Progesterone-mediated oocyte maturation’, whereas the downregulated DEGs were mainly involved in ‘Tyrosine metabolism’, ‘Focal adhesion’ and the ‘Ras signaling pathway’ (Table I).

Table I.

GO and KEGG pathway enrichment analysis of DEGs in ACC samples.

A, Upregulated
ID	Description	Count	P-value
GO:0051301	Cell division	8	5.60×10−7
GO:0000083	Regulation of transcription involved in G1/S transition of mitotic cell cycle	4	5.06×10−6
GO:0000082	G1/S transition of mitotic cell cycle	5	1.47×10−5
GO:0005634	Nucleus	23	1.26×10−7
GO:0005654	Nucleoplasm	17	4.07×10−7
GO:0072686	Mitotic spindle	4	3.33×10−5
GO:0019901	Protein kinase binding	5	2.83×10−3
GO:0035173	Histone kinase activity	2	6.38×10−3
GO:0005515	Protein binding	21	1.64×10−2
Hsa04115	p53 signaling pathway	3	3.19×10−3
Hsa04914	Progesterone-mediated oocyte maturation	3	3.19×10−3
Hsa04114	Oocyte meiosis	3	5.33×10−3
B, Downregulated
ID	Description	Count	P-value
GO:0014068	Positive regulation of phosphatidylinositol 3-kinase signaling	6	4.30×10−4
GO:0006006	Glucose metabolic process	6	4.95×10−4
GO:0007165	Signal transduction	22	4.58×10−3
GO:0005615	Extracellular space	33	1.37×10−6
GO:0005578	Proteinaceous extracellular matrix	13	2.60×10−4
GO:0005576	Extracellular region	30	6.40×10−4
GO:0005509	Calcium ion binding	16	3.02×10−3
Hsa00350	Tyrosine metabolism	4	4.57×10−3
Hsa04510	Focal adhesion	7	1.44×10−2
Hsa04014	Ras signaling pathway	7	2.17×10−2

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes.

A total of 113 DEGs, consisting of 19 up and 94 downregulated genes, were filtered into the PPI network using the STRING database. The network contained 113 nodes and 272 edges, and were visualized using Cytoscape software (Fig. 1B). The most significant module was obtained using Cytoscape and it contained 14 nodes and 88 edges (Fig. 1C). Functional analysis of the genes in the most significant module was performed using DAVID (Fig. 2A; Table II). BP analysis revealed that genes in this module were mainly enriched in ‘Cell division’, ‘Anaphase-promoting complex-dependent catabolic process’, ‘Regulation of transcription involved in G1/S transition of mitotic cell cycle’, ‘G2/M transition of mitotic cell cycle’ and ‘Protein ubiquitination involved in ubiquitin-dependent protein catabolic process’ (Fig. 2A). The CC, MF and KEGG pathway analyses of hub genes are presented in Table II. KEGG pathway analysis revealed that those genes were predominantly enriched in the p53 signaling pathway, progesterone-mediated oocyte maturation, oocyte meiosis and cell cycle (Table II).

Figure 2.

Interaction network and BPs of the hub genes. (A) The most significant BPs of the hub genes (P<0.05). (B) The hub genes and their co-expression genes were analyzed using cBioPortal. Nodes with a bold black outline represent hub genes, whereas nodes with a thin black outline represent the co-expression genes. (C) Hierarchical clustering of hub genes was constructed using UCSC. The upregulated genes are marked in red and downregulated genes in blue. The samples were all ACC primary samples. The tumor stages are presented as follows: Light purple bar, stage I; blue bar, stage II; red bar, stage III; and orange bar, stage IV. BP, biological process; ACC, adrenocortical carcinoma.

Table II.

GO and KEGG pathway enrichment analysis of DEGs in the most significant module.

ID	Description	Count in gene set	P-value
GO:0005634	Nucleus	14	1.39×10−7
GO:0005829	Cytosol	12	3.86×10−7
GO:0005654	Nucleoplasm	11	1.25×10−6
GO:0019901	Protein kinase binding	5	1.48×10−4
GO:0035173	Histone kinase activity	2	3.08×10−3
GO:0005515	Protein binding	12	1.62×10−2
Hsa04115	P53 signaling pathway	3	1.36×10−3
Hsa04914	Progesterone-mediated oocyte maturation	3	2.28×10−3
Hsa04114	Oocyte meiosis	3	3.55×10−3
Hsa04110	Cell cycle	3	4.57×10−3

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes.

A total of 14 genes were identified as hub genes (degrees ≥10). A network of the hub and their co-expression genes was analyzed using the cBioPortal online platform and is illustrated in Fig. 2B. Hierarchical clustering showed that the hub genes were associated with a high clinical stage of ACC (Fig. 2C). In addition, the overall survival analysis of the hub genes performed using Kaplan-Meier analysis in cBioPortal demonstrated that patients with ACC with upregulation of AURKA, TYMS, GINS1, RACGAP1, RRM2, EZH2, ZWINT, CDK1, CCNB1, NCAPG and TPX2 presented with a decreased overall and disease-free survival (Fig. 3). Patients with ACC with upregulation of Mitotic arrest deficient 2 like 1 (MAD2L1) and PRC1 genes presented with a worse disease-free survival (Fig. 3B). Nonetheless, no change was observed in the patients with ACC with SMC4 alterations, according to the cBioPortal data (P=0.527 for overall survival and P=0.429 for disease-free survival; Fig. 3). In addition, Oncomine analysis of cancer vs. normal tissue showed that higher mRNA expression levels of those hub genes increased tumor Weiss grade in the Giordano Adrenal 2 dataset (Fig. 4). Of note, the mRNA expression levels of TYMS, GINS1, RRM2, ZWINT and SMC4 analyzed in the Oncomine database were observed to be significantly increased in ACC compared with normal tissues (Fig. 5).

Figure 3.

(A) Overall and (B) disease-free survival analysis of hub genes was performed using the cBioPortal online platform. P<0.05 was considered to indicate a statistically significant difference.

Figure 4.

Association between the expression of each hub gene and tumor Weiss grade in the Giordano Adrenal 2 dataset. 0, no value (n=32); 1, low (n=13); and 2, high (n=20).

Figure 5.

Oncomine analysis of the mRNA expression levels of (A) TYMS, (B) GINS1, (C) RRM2, (D) ZWINT and (E) SMC4 genes in clinical adrenocortical carcinoma samples vs. normal tissues. Red indicates high expression and blue indicates low expression. The rank for a gene is the median rank for that gene across each of the analyses. The P-value for a gene is its P-value for the median-ranked analysis. 1. Adrenal Cortex Carcinoma vs. Normal (Giordano Adrenal, Am J Pathol, 2003); 2. Adrenal Cortex Carcinoma vs. Normal (Giordano Adrenal 2, Clin Cancer Res, 2009).

Discussion

ACC is associated with a poor prognosis, limited treatment options and high tumor recurrence rates (1–3). The pathogenetic mechanisms of ACC includes alterations of the Insulin-like Growth Factor system (35), Wnt/β-catenin pathway activation (36), TP53 mutations and prognostic molecular markers involved in cancer cell invasion properties and angiogenesis, appear to be very promising in elucidating of tumorigenesis and progression of ACC (37,38). However, the molecular mechanisms of ACC remain poorly understood. The identification of biomarkers associated with ACC tumorigenesis, progression and prognosis are urgently required.

Microarray technology combined with bioinformatics analysis has enabled researchers to explore genetic alterations, and has been proven to be a useful approach in identifying novel biomarkers in several diseases, such as hepatocellular carcinoma and adrenocortical tumors (22,35). In the present study, a total of 228 DEGs were identified, 14 of which were selected as hub genes (degrees ≥10). BP analysis suggested that these hub genes were significantly enriched in cell division and the mitotic cell cycle, which indicated that the deregulation of the cell cycle may serve a key role in the tumorigenesis and development of ACC. The present study additionally combined various databases to identify and validate the diagnostic and prognostic value of hub genes in ACC. Kaplan-Meier analysis revealed that the expressions of AURKA, TYMS, GINS1, RACGAP1, RRM2, EZH2, ZWINT, CDK1, CCNB1, NCAPG and TPX2 were negatively associated with overall and disease-free survival, suggesting these genes may exert pivotal functions in the progression of ACC.

Some of these hub genes have previously been identified as biomarkers for ACC (18,39–41). For example, AURKA, which regulates cell cycle and meiotic division, was overexpressed in pediatric adrenocortical tumors, suggesting it may be associated with more aggressive disease and poor prognosis, and could help develop an interesting therapeutic approach against ACC (39,40). MAD2L1 and CCNB1 have also been reported as potential markers for differentiating ACCs from adenomas (18). In particular, overexpressed CCNB1 dysregulated the cell cycle in the G2-M phase transition, with poor survival in the majority of solid tumors (41). Similarly, in the present study, upregulation of MAD2L1 and CCNB1 in tumor tissues predicted a worse overall and disease-free survival in patients with ACC using the cBioPortal platform, which indicated a poor prognosis. CDK1 serves an important role in regulating cell cycle progression and mediating the phosphorylation of Bcl-2 by binding with cyclin B to form a complex called cyclin B-CDK1 (42). In adrenocortical tumors, CDK1 overexpression is associated with tumor suppressor miR-7 downregulation, which may serve as a target for inhibiting the progression of ACC (43,44). EZH2 was significantly associated with poorer outcomes in ACC (45). A recent study by Drelon et al (46) reported that EZH2, as a deregulated histone modifier, deregulated activity of the P53/RB/E2F pathway and WNT signaling modulation to promote cell proliferation, which may be a new therapeutic target for ACC. In addition, Yuan et al (15) underlined the potential of TPX2, PRC1 and RACGAP1 as markers for the diagnosis and prognosis of ACC, which was consistent with the hypothesis of the present study.

Although six other hub genes (TYMS, RRM2, ZWINT, GINS1, SMC4 and NCAPG) have not been extensively reported to participate in ACC progression, they were observed to be involved in various tumors. Using Oncomine evaluation, the present study identified that the mRNA expression levels of TYMS, GINS1, RRM2, ZWINT and SMC4 were higher in ACC compared with normal tissues. In addition, Oncomine analysis of ACC vs. normal tissues revealed that the upregulation of hub genes was significantly associated with a higher Weiss grade.

Among these six hub genes, the elevated expression of TYMS has also been reported in lung (47), gastric (48), colorectal (49,50), renal cell (51) and prostate cancer (52), suggesting it may serve as a valuable biomarker for the diagnosis, treatment and prognosis of tumors. In the present study, the PPI network demonstrated that TYMS directly interacted with other hub genes, such as CDK1, AURKA and PRC1, and that it may affect cell proliferation through modulation of the cell cycle and multiple signaling pathways. In addition, overexpression of TYMS was significantly associated with a shorter survival time and higher tumor Weiss grade, indicating a key role in the tumorigenesis or progression of ACC. RRM2 serves a key regulatory role in DNA synthesis and cell proliferation (53). RRM2 has been hypothesized to promote angiogenesis by producing reactive oxygen species to activate the ERK1/2 signaling pathway, and inducing HIF-1α and VEGF expression in human cervical cancer, which is associated with a poor outcome in certain types of cancer (54). ZWINT is essential for mitotic checkpoint signaling (55,56). It was recently demonstrated that ZWINT was significantly correlated with the expression of cell-cycle proteins such as PCNA, cyclin B1, Cdc25C and CDK1, and may be considered as potential therapeutic targets for hepatocellular carcinoma (55). GINS1, also known as PSF1 (57), was significantly associated with a worse overall survival, but not disease-free survival in the present study. However, PSF1 was highly expressed in several types of cancer (58–61). In addition, previous studies have demonstrated that the transcriptional activity of the PSF1 gene was associated with cancer cell malignancy by affecting the cell cycle and proliferation, highlighting it as a potentially useful biomarker for the identifying patients who may have unfavorable prognoses (59,61). SMC4 has been reported to be involved in tumor cell growth, migration and invasion (62,63). However, its role in ACC has yet to be completely elucidated. In the present study, although SMC4 alteration was not significantly associated with a worse overall and disease-free survival, hierarchical clustering for hub genes and data from Oncomine indicates that it serves a crucial role in ACC tumorigenesis. NCAPG organizes the coiling topology of individual chromatids during cell mitosis and meiosis, which is involved in the progression of liver carcinoma (64). Therefore, it was speculated that they may serve a critical role in the carcinogenesis and progression of ACC.

In conclusion, the present study identified and analyzed key biomarkers in ACC using bioinformatics analysis. Two databases were combined to screen 228 DEGs, and 14 hub genes, which may be regarded as powerful and promising biomarkers for predicting tumorigenesis and progression of ACC, were identified. These hub genes were associated with tumor cell proliferation and cell cycle regulation. Of note, candidate hub gene upregulation was associated with a worse survival rate and higher Weiss grade; and this may provide a basis for further clinical molecular target therapy experiments and diagnostic approaches for ACC, if these potential genes are developed as novel useful diagnostic as well as prognostic markers and the underlying pathological causative pathways or involved signaling targets are elucidated. However, further studies are required to confirm the biological functions and mechanisms of action of these genes in ACC.