Abstract Adrenocortical carcinoma (ACC) is a rare endocrine malignancy with poor prognosis. The disease originates from the cortex of adrenal gland and lacks effective treatment. Efforts have been made to elucidate the pathogenesis of ACC, but the molecular mechanisms remain elusive. To identify key genes and pathways in ACC, the expression profiles of [37]GSE12368, [38]GSE90713 and [39]GSE143383 were downloaded from the Gene Expression Omnibus (GEO) database. After screening differentially expressed genes (DEGs) in each microarray dataset on the basis of cut-off, we identified 206 DEGs, consisting of 72 up-regulated and 134 down-regulated genes in three datasets. Function enrichment analyses of DEGs were performed by DAVID online database and the results revealed that the DEGs were mainly enriched in cell cycle, cell cycle process, mitotic cell cycle, response to oxygen-containing compound, progesterone-mediated oocyte maturation, p53 signaling pathway. The STRING database was used to construct the protein–protein interaction (PPI) network, and modules analysis was performed using Cytoscape. Finally, we filtered out eight hub genes, including CDK1, CCNA2, CCNB1, TOP2A, MAD2L1, BIRC5, BUB1 and AURKA. Biological process analysis showed that these hub genes were significantly enriched in nuclear division, mitosis, M phase of mitotic cell cycle and cell cycle process. Violin plot, Kaplan-Meier curve and stage plot of these hub genes confirmed the reliability of the results. In conclusion, the results in this study provided reliable key genes and pathways for ACC, which will be useful for ACC mechanisms, diagnosis and candidate targeted treatment. Keywords: adrenocortical carcinoma, gene expression omnibus, differentially expressed genes, protein-protein interaction, Kaplan-Meier curve Introduction Adrenocortical carcinoma (ACC) is a rare endocrine malignancy with poor prognosis ([40]1). According to the new WHO criteria, ACC is subclassified according to its morphological characteristics, including conventional, oncocytic, myxoid and sarcomatoid subtypes ([41]2).The incidence of ACC is estimated to be 1 to 2 per million per year ([42]3). Women appear to have a slightly higher incidence rate compared to men, with a ratio of 1.5-2.5 to 1 ([43]1). ACC mainly presents as an incidental imaging finding ([44]4), and unfortunately, the 5-year survival rate is less than 40% in most literature reports ([45]5–[46]8). However, there are few treatments available for this disease with such low survival rates ([47]9). At present, surgical resection is considered the primary approach for treating ACC. Complete removal of the primary tumor is crucial for improving prognosis, while 5-year survival rates are as low as 20% and 10%, respectively, for patients with microscopic or macroscopic involvement at the tumor margin ([48]6). However, despite undergoing complete resection, many patients with ACC experience recurrence and disease progression, so adjuvant therapy is often recommended. The medication that is commonly used for adjuvant therapy is mitotane, which is the only approved treatment for advanced ACC, but it is less effective in the advanced stages of the disease and has toxicity and side effects ([49]10). Another effective adjuvant therapy is radiotherapy, which has been shown to reduce local recurrence rates but not improve clinical outcomes ([50]11–[51]14). Overall, regardless of the treatment, the prognosis for advanced ACC remains poor. In recent years, there have been ongoing efforts to gain a deeper understanding of ACC in terms of its tumorigenesis and progression mechanisms. For example, a genome-wide methylation analysis shed light on the prognostic implications of hypermethylated ACC. It was found that hypermethylation in ACC primarily leads to the silencing of specific tumor suppressor genes, thereby contributing to a poorer prognosis ([52]15). This research highlights the importance of epigenetic alterations, specifically DNA methylation, in ACC. Transcriptomic analysis has revealed that somatic inactivation mutations of the tumor suppressor gene TP53 and activation mutations of the proto-oncogene β-catenin (CTNNB1) are the most frequently observed mutations in ACC, these mutations have been associated with a poor prognosis in ACC patients ([53]16). Whole-exome sequencing has also been utilized to genetically characterize potential somatic mutations in ACC, in order to obtain genomic landscape of ACC and facilitate the identification of frequently altered pathways, which may guide targeted therapy research ([54]17). However, it is concerning that the various advancements in the international scientific community in the field of epigenetics, transcriptomics, and exome sequencing have not been reported. Consequently, there is an urgent need to identify new candidate targets for ACC. To address this gap, it is crucial to explore the potential of high-throughput techniques such as RNA-sequencing and microarrays, which have enabled significant insights into the molecular mechanisms of tumors over the past few decades. By analyzing DEGs and functional pathways between tumor and normal tissues, we can gain a deeper understanding of the pathogenesis and progression of tumor. In this study, we conducted an integrated bioinformatics analysis to identify hub genes and pathways in ACC, which can greatly contribute to our understanding of the underlying mechanisms and candidate targeted therapy options for ACC. Materials and methods Microarray data Gene expression profiles of ACC were obtained from the GEO database ([55]https://www.ncbi.nlm.nih.gov/geo/). The screening criteria were as follows: i) adrenocortical carcinoma, ii) samples containing tumor and normal adrenal gland tissues, iii) expression profiling by array and iiii) Homo sapiens. A total of three GEO datasets were selected, including [56]GSE12368 ([57]18), [58]GSE90713 ([59]19), and [60]GSE143383 ([61]20). The [62]GSE12368 dataset is based on the [63]GPL570 platform and contains 12 ACC samples (5 male samples) and 6 normal samples (0 male samples). [64]GSE90713 dataset is based on the [65]GPL15207 platform and contains 58 ACC samples and 5 normal samples. [66]GSE143383 dataset is based on the [67]GPL16043 platform and contains 57 ACC samples (17 male samples) and 5 normal samples (3 male samples). Data preprocessing and identification of differentially expressed genes (DEGs) The DEGs between ACC and normal samples in each GEO dataset were screened using GEO2R ([68]http://www.ncbi.nlm.nih.gov/geo/geo2r). GEO2R is an online web tool that compares the datasets in a GEO series to identify DEGs under different conditions. Genes with more than one probe set were averaged and probe sets without corresponding gene symbols were removed. |log[2]FC (fold change) |≥1 and adjusted P-values (adj. P) <0.05 were considered statistically significant. GO and KEGG enrichment analyses of DEGs GO is a main bioinformatics tool to analyze the biological process (BP), cellular component (CC) and molecular function (MF) of genes ([69]21, [70]22). KEGG is an online database for gene functional enrichment ([71]23). GO and KEGG enrichment analysis of DEGs were performed using DAVID (version 2021) ([72]24) database ([73]https://david.ncifcrf.gov/). P<0.05 was considered statistically significant. PPI network construction and module analysis Search Tool for the Retrieval of Interacting Genes (STRING) (version 11.5) ([74]25) online database ([75]http://www.string-db.org/) was used to construct the protein-protein interaction (PPI) network, thereby providing insights into the occurrence and development of diseases. Then, the PPI network was visualized by the Cytoscape (version 3.9.1) software ([76]26) ([77]http://www.cytoscape.org/) to further analyze the complex relationships among the DEGs. Significant modules in the PPI network were identified using plug-in Molecular Complex Detection (MCODE) (version 2.0.0) ([78]27) (Parameters were the default values) of Cytoscape. Subsequently, the GO and KEGG analyses were performed for genes in the most significant module using DAVID. Hub genes selection and analysis The plug-in cytoHubba ([79]28) of Cytoscape was used to screen hub genes based on the PPI network. Cytohubba ranks the importance of each gene according to different algorithms. Among 12 computing methods, we selected the maximal clique centrality (MCC) method to identify hub genes. Then, the BP analysis of hub genes was visualized using plug-in Biological Networks Gene Oncology tool (BiNGO)(version3.0.3) ([80]29) of Cytoscape. To verify the roles of hub genes, the FPKM gene expression profiles of GDC TCGA Adrenocortical Cancer (ACC) (14 datasets) and CTEX (11 datasets) were downloaded using the University of California, Santa Cruz (UCSC) Xena browser ([81]30) ([82]http://genome-cancer.ucsc.edu). The datasets contained 79 ACC samples and 128 normal samples. The violin plot was used to show the gene expression of hub genes in ACC and normal adrenal samples. The overall survival analyses and disease free survival analyses of hub genes were performed using Kaplan-Meier curve in cBioPortal (version 5.2.2) ([83]31, [84]32)([85]http://www.cbioportal.org/) online platform. The relationship between hub genes and tumor stage in ACC patients was presented using stage plot in Gene Expression Profiling Interactive Analysis (GEPIA) ([86]33) ([87]www.gepia.cancer-pku.cn). Results Identification of DEGs A total of 127 ACC and 16 normal tissues were obtained from the three GEO datasets. Using |log[2]FC|≥1 and adj. P <0.05 as cut-off criteria, DEGs (970 in [88]GSE12368, 581 in [89]GSE90713 and 784 in [90]GSE143383) were identified, including 397 up-regulated and 573 down-regulated genes in [91]GSE12368 dataset, 161 up-regulated and 420 down-regulated genes in [92]GSE90713 dataset, 250 up-regulated and 534 down-regulated genes in [93]GSE143383 dataset. The distributions of DEGs in three datasets are shown in the volcano plots ([94] Figures 1A–C ), respectively. The overlap among the three datasets contained 206 genes as shown in the Venn diagram ([95] Figure 1D ), including 72 up-regulated and 134 down-regulated DEGs ([96] Table 1 ) between ACC tissues and normal tissues. The expression heatmap of the top 10 up-regulated DEGs and top 10 down-regulated DEGs is shown in [97]Figure 1E . Figure 1. [98]Figure 1 [99]Open in a new tab Identification of DEGs between ACC and normal tissues. The volcano plots of DEGs for dataset [100]GSE12368 (A), [101]GSE90713 (B) and [102]GSE143383 (C). (D) The Venn diagram of the overlapping DEGs among the three datasets. (E) The expression heatmap of the top 10 up-regulated and top 10 down-regulated DEGs. Table 1. Two hundred and six differentially expressed genes (DEGs) were identified and confirmed from three profile datasets, including 72 up-regulated genes and 134 down-regulated genes in the ACC tissues, compared to normal tissues. Regulation DEGs (gene symbol) Up-regulated TPX2, DEPDC1B, CCNB1, GINS1, ANLN, BIRC5, PLEKHG4, CDK1, KIF11, RACGAP1, CDCA4, C12orf76, TMEM161B, CDC6, KIF2A, SMC4, PLXNC1, AURKA, DHFR, MAD2L1, TMPO, STC1, GAS2L3, SULF2, C5orf34, TYMS, TMEM106C, ZWINT, NDC80, H2AFZ, RFC4, CCNA2, TMEM132A, BUB1, PBK, SPATS2, TRIP13, CDCA5, UBE2T, ATF7IP, CRNDE, ANGPT2, MYBL1, TIGD7, DIAPH3, GGH, SHCBP1, ASPM, GADD45A, MDM2, UBE2C, MND1, CCNB2, PRC1, CEP55, TOP2A, FANCI, ZNF367, RAD51AP1, RNASEH2A, DTL, HMMR, GJC1, GMNN, ENC1, CENPK, ZNF404, IGFBP3, KIAA0101, CDKN3, NCAPG, ESM1 Down- regulated ANKS1A, MMP2, MUM1L1, IGF1, GOLM1, PHYHIP, NFIA, NR2F1, KIAA1210, SLCO2B1, TEK, CD55, SERPINB9, MRPL33, SLCO2A1, SERPING1, PARVA, SLC40A1, CRHBP, ALAD, ZNF331, EBPL, IGSF11, ZBED1, SLITRK4, CYP11B2, ITGA8, TLE2, STAT5A, DUSP26, SORBS2, USP13, STEAP4, KIAA1217, OMD, FBLN1, SLC16A9, PDGFD, SPON1, PLA2G4A, LMOD1, ABCC3, C1GALT1C1, ACADVL, TBXAS1, HOXA5, GALM, CTH, NPY1R, PLCXD3, ECHDC3, LUZP2, KLHL4, USP53, LIMS1, DNAJC12, GRAMD4, C7, IGFBP5, CTGF, EPHX2, CREM, KCNQ1, NR2F2, ANGPTL1, CREBL2, PLIN1, TAGLN, EMCN, KCNK2, FOSL2, TMEM150C, DLG2, AXL, SHE, GIPC2, OGN, ACADSB, MYLK, FZD5, IGFBP6, SLC37A2, FNDC5, ARMC9, HS6ST1, CXCL12, GLUL, EHD3, EFEMP2, ALDH1A1, ADORA3, AEBP1, NEXN, NIT1, WISP1, CBLN4, INMT, BRE, SEC62, TLR4, FMO2, ZNF185, CDC42EP4, RBMS3, CASP9, ADH1B, OLFML3, CAB39L, PDGFRA, CDH19, SIRT1, FAM65C, C11orf96, RSPO3, CPT1A, SREBF1, PRELP, WFDC1, NANOS1, PTH1R, CYBRD1, CRISPLD2, ATP1B3, C2orf40, AOX1, DAPL1, PHF11, HGF, APOC1, C11orf54, FLVCR2, PARM1, PTGIS, FBLN5 [103]Open in a new tab KEGG and GO enrichment analyses of DEGs The BP of the GO enrichment analysis showed that up-regulated DEGs were significantly enriched in cell cycle, cell cycle process, mitotic cell cycle, regulation of cell cycle and mitotic cell cycle process ([104] Figure 2A ) and down-regulated DEGs were enriched in response to oxygen-containing compound, lipid metabolic process, circulatory system development, vasculature development and cardiovascular system development ([105] Figure 2B ). In terms of CC, up-regulated DEGs were enriched in cytosol, nucleoplasm, chromosome, microtubule cytoskeleton and chromosomal part ([106] Figure 2C ), while the down-regulated DEGs were significantly enriched in extracellular region, extracellular region part, extracellular organelle, extracellular vesicle and extracellular exosome ([107] Figure 2D ). Up-regulated DEGs were significantly enriched in enzyme binding, adenyl nucleotide binding, adenyl ribonucleotide binding, ATP binding and kinase binding in the MF categories ([108] Figure 2E ), the most significantly enriched MF terms for down-regulated genes were identical protein binding, transition metal ion binding, protein complex binding, cofactor binding, coenzyme binding and extracellular matrix structural constituent ([109] Figure 2F ). The BP of the GO enrichment analysis showed that up-regulated and down-regulated DEGs were significantly enriched in spindle assembly involved in female meiosis I, intestinal epithelial cell maturation, positive regulation of exit from mitosis, mitotic sister chromatid segregation, mitotic spindle assembly checkpoint, animal organ regeneration, regulation of cyclin-dependent protein serine/threonine kinase activity, mitotic spindle organization, positive regulation of phosphatidylinositol 3-kinase signaling, positive regulation of fibroblast proliferation, mitotic cell cycle, cell division ([110] Figure 3A ). In terms of CC, up-regulated and down-regulated DEGs were enriched in insulin-like growth factor binding protein complex, insulin-like growth factor ternary complex, elastic fiber, spindle microtubule, chromosome, centromeric region, midbody, microtubule cytoskeleton, spindle, extracellular space, nucleoplasm, cytoplasm, nucleus ([111] Figure 3B ). Up-regulated and down-regulated DEGs were significantly enriched in insulin-like growth factor II binding, insulin-like growth factor I binding, ubiquitin-like protein ligase binding, cyclin-dependent protein serine/threonine kinase regulator activity, insulin-like growth factor binding, flavin adenine dinucleotide binding, monooxygenase activity, fibronectin binding, integrin binding, actin binding, protein kinase binding, identical protein binding in the MF categories ([112] Figure 3C ). According to KEGG pathway enrichment analysis, up-regulated DEGs were mainly enriched in cell cycle, progesterone-mediated oocyte maturation, p53 signaling pathway, cellular senescence, oocyte meiosis ([113] Figure 4A ), while the down-regulated DEGs were significantly enriched in PI3K-Akt signaling pathway, focal adhesion, Ras signaling pathway, Alcoholic liver disease, EGFR tyrosine kinase inhibitor resistance ([114] Figure 4B ). Up-regulated and down-regulated DEGs were significantly enriched in cell cycle, PI3K-Akt signaling pathway, progesterone-mediated oocyte maturation, p53 signaling pathway, cellular senescence, oocyte meiosis, melanoma, FoxO signaling pathway, AMPK signaling pathway, EGFR tyrosine kinase inhibitor resistance, fatty acid degradation, antifolate resistance ([115] Figure 3D ). Figure 2. [116]Figure 2 [117]Open in a new tab GO (Gene Ontology) pathway enrichment analyses of DEGs (differentially expressed genes). (A) Biological process for up-regulated DEGs. (B) Biological process for down-regulated DEGs. (C) Cellular component for up-regulated DEGs. (D) Cellular component for down-regulated DEGs. (E) Molecular function for up-regulated DEGs. (F) Molecular function for down-regulated DEGs. The X-axis represents the ratio of enriched pathways to total pathways, while the Y-axis displays the names of enriched pathways. “Count” indicates the number of enriched pathways, which corresponds to the size of the dot. The color indicates the P-value, with darker shades indicating smaller P-value. Figure 3. [118]Figure 3 [119]Open in a new tab GO (Gene Ontology) pathway enrichment analyses and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses of DEGs (differentially expressed genes). (A) Biological process for up-regulated and down-regulated DEGs. (B) Cellular component for up-regulated and down-regulated DEGs. (C) Molecular function for up-regulated and down-regulated DEGs. (D) Bubble plot of KEGG for up-regulated and down-regulated DEGs. Figure 4. [120]Figure 4 [121]Open in a new tab KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses of DEGs (differentially expressed genes). (A) Bubble plot of KEGG for up-regulated DEGs. (B) Bubble plot of KEGG for down-regulated DEGs. PPI network construction and module analysis The PPI network of the DEGs with 150 nodes and 1094 edges was constructed ([122] Figure 5A ). In addition, we obtain top three modules using MCODE in Cytoscape, and the genes of three modules were shown in [123]Table 2 . The most significant module is shown in [124]Figure 5B and the PPI network of this module consisted of 41 nodes and 780 edges. GO and KEGG pathway enrichment analysis of DEGs involved in the most significant module were also analyzed using DAVID and results were mainly enriched in cell cycle process, cell cycle, mitotic cell cycle process, nuclear division, progesterone-mediated oocyte maturation ([125] Table 3 ). Figure 5. [126]Figure 5 [127]Open in a new tab PPI networks of DEGs. (A) Up-regulated and down-regulated genes in PPI network (red nodes mean up-regulated genes, and blue nodes mean down-regulated genes). (B) The most significant module in PPI network. (C) Hub genes in PPI network. (D) Hub genes and other genes in PPI network (orange nodes mean hub genes). Table 2. The top three significant modules. Modules Nodes Edges Genes Module1 41 780 SMC4, CDCA5, ZWINT, TYMS, KIAA0101, HMMR, PBK, CENPK, DEPDC1B, CDC6, AURKA, MND1, FANCI, TOP2A, CCNA2, GINS1, CDK1, NDC80, CCNB1, CEP55, CCNB2, MAD2L1, TRIP13, BIRC5, CDKN3, UBE2T, ANLN, KIF11, SHCBP1, RFC4, UBE2C, DTL, RNASEH2A, RAD51AP1, TPX2, ASPM, PRC1, NCAPG, BUB1, GMNN, RACGAP1 Module2 6 12 ANGPT2, CTGF, MMP2, HGF, IGFBP3, TEK Module3 4 5 STAT5A, TLR4, CXCL12, PDGFRA [128]Open in a new tab Table 3. GO and KEGG pathway enrichment analysis of DEGs in the most significant module. pathway ID Pathway description Count in gene set FDR GO:0022402 Cell cycle process 32 3.72E-27 GO:0007049 Cell cycle 33 1.02E-24 GO:1903047 Mitotic cell cycle process 27 1.02E-24 GO:0000280 Nuclear division 24 3.76E-24 GO:0048285 Organelle fission 24 2.11E-23 GO:0000278 Mitotic cell cycle 27 5.60E-23 GO:0051301 Cell division 23 1.16E-20 GO:0051726 Regulation of cell cycle 27 2.02E-20 GO:0007067 Mitotic nuclear division 19 1.12E-19 GO:0000819 Sister chromatid segregation 17 6.54E-19 Hsa04914 Progesterone-mediated oocyte maturation 7 1.20E-06 Hsa04110 Cell cycle 7 2.13E-06 Hsa04114 Oocyte meiosis 6 5.75E-05 [129]Open in a new tab Annotation: FDR (false discovery rate, which corrects the P-value by using multiple hypothesis tests on the P-value to reduce the false positive rate). Hub gene selection and analysis A total of eight genes were identified as hub genes using the cytoHubba in Cytoscape. They were CDK1, CCNA2, CCNB1, TOP2A, MAD2L1, BIRC5, BUB1 and AURKA. The PPI networks of hub genes and its interaction and other genes were shown in [130]Figures 5C, D respectively. All eight hub genes were up-regulated and had high degree of connectivity in the network ([131] Table 4 ). The abbreviations, full names, aliases and functions for the hub genes are shown in [132]Table 5 . The BP analysis of hub genes was visualized using BiNGO and the result showed that these hub genes were significantly enriched in nuclear division, mitosis, M phase of mitotic cell cycle and cell cycle process ([133] Figure 6A ). To verify the roles of hub genes, the gene expression profiles of GDC TCGA and CTEX were downloaded. The datasets contained a total of 79 ACC samples and 128 normal samples. The violin plot was used to show the differential gene expression of hub genes in ACC and normal samples ([134] Figure 6B ). Subsequent survival analysis of these hub genes was performed using cBioPortal, revealing that ACC patients with alterations in CDK1, CCNA2, CCNB1, BIRC5, BUB1, and AURKA (depicted by the red solid line) demonstrated significantly lower overall survival rates and disease free survival rates compared to patients without alterations in these genes (depicted by the blue solid line). Similarly, overall survival was reduced in ACC patients with TOP2A alteration. The results of TOP2A progression free survival demonstrated that ACC patients with TOP2A alterations had lower progression free survival ([135] Figures 6C, D ). The relationship between the expression of hub genes and ACC tumor stage is showed in [136]Figure 6E , where higher ACC tumor stage is associated with higher hub gene expression. Table 4. Hub genes with high degree of connectivity. Gene symbol Degree Type CDK1 53 Up CCNA2 50 Up CCNB1 50 Up TOP2A 47 Up MAD2L1 47 Up BIRC5 47 Up BUB1 46 Up AURKA 46 Up [137]Open in a new tab Table 5. Functional roles of 8 hub genes. No. Gene symbol Full name Aliases Function References