Skip to main content

ORIGINAL RESEARCH article

Front. Genet., 15 August 2018
Sec. Systems Endocrinology

Co-expression Network Analysis of Biomarkers for Adrenocortical Carcinoma

\r\nLushun YuanLushun Yuan1Guofeng QianGuofeng Qian2Liang ChenLiang Chen1Chin-Lee WuChin-Lee Wu3Han C. DanHan C. Dan4Yu Xiao,,*Yu Xiao1,5,6*Xinghuan Wang*Xinghuan Wang1*
  • 1Department of Urology, Zhongnan Hospital of Wuhan University, Wuhan, China
  • 2Department of Endocrinology, The First Affiliated Hospital of Zhejiang University, Hangzhou, China
  • 3Department of Urology, Massachusetts General Hospital, Harvard Medical School, Boston, MA, United States
  • 4Greenebaum Cancer Center, School of Medicine, University of Maryland, Baltimore, MD, United States
  • 5Laboratory of Precision Medicine, Zhongnan Hospital of Wuhan University, Wuhan, China
  • 6Department of Biological Repositories, Zhongnan Hospital of Wuhan University, Wuhan, China

Adrenocortical carcinoma (ACC) is a rare malignancy with a poor prognosis. And currently, there are no specific diagnostic biomarkers for ACC. In our study, we aimed to screen biomarkers for disease diagnosis, progression and prognosis. We firstly used the microarray data from public database Gene Expression Omnibus database to construct a weighted gene co-expression network, and then to identify gene modules associated with clinical features of ACC. Though this algorithm, a significant module with R2 = 0.64 (P = 9 × 10-5) was identified. Co-expression network and protein–protein interaction network were performed for screen the candidate hub genes. Checked by The Cancer Genome Atlas (TCGA) database, another independent dataset GSE19750, and GEPIA database, using one-way ANOVA, Pearson’s correlation, survival analysis, diagnostic capacity (ROC curve) and expression level revalidation, a total 12 real hub genes were identified. Gene ontology and KEGG pathway analysis of genes in the significant module revealed that the hub genes are significantly enriched in cell cycle regulation. Moreover, gene set enrichment analysis suggests that the samples with highly expressed hub genes are correlated with cell cycle. Taken together, our integrated analysis has identified 12 hub genes that are associated with the progression and prognosis of ACC; these hub genes might lead to poor outcomes by regulating the cell cycle.

Introduction

Adrenocortical carcinoma is a rare malignancy. Its incidence is approximately 0.5–2 new cases per million people per year (Fassnacht et al., 2011), with an increased occurrence during childhood and the fourth to fifth decades of life. Most ACCs occur sporadically, but some cases are associated with various genetic diseases, e.g., Li-Fraumeni syndrome (LFS), Beckwith–Wiedemann, and multiple endocrine neoplasia type I (MEN1). ACCs are often aggressive, cannot be completely resected and have a poor prognosis (Fassnacht et al., 2011; Else et al., 2014). In a study of 330 ACCs from M. D. Anderson, the median survival was 3.2 years and median overall survival for stage I, II, III, and IV were 24.1, 6.08, 3.47, and 0.89 years, respectively (Ayala-Ramirez et al., 2013). The 275 surgically resected tumors had a median local-recurrence free time of 1 year. In addition, a study of 3982 ACCs from 1985 to 2005 in United States exhibited that the survival did not improve (Beuschlein et al., 2015). Based on these facts, it is urgent to investigate the mechanisms that promote ACC progression and to discover novel molecular biomarkers for the prognosis.

With the development of high-throughput microarray technology, many studies have reported the genes involving in tumor progression (Chen L. et al., 2017; Robertson et al., 2017). Gene expression profiles have been used to identify genes associated with progression of ACCs in some previous studies (Pinto et al., 2015; Kulshrestha et al., 2016). However, these studies were mostly concerned with differentially expressed genes and did not consider the high degree of interconnection between genes, where genes with similar expression patterns may be functionally related. Merely screening the differentially expressed genes in normal and tumor samples had limitations and we should pay more attention to correlation between gene expression and clinical features. The algorithm, WGCNA can construct free-scale gene co-expression networks to explore the relationships between different gene sets or between gene sets and clinical features (Langfelder and Horvath, 2008). Currently, WGCNA is used to study various cancers and has identified biomarkers associated with critical features. However, it was not used in ACC samples and our group was the first to use this algorithm in ACC. Thus, we attempt to construct a co-expression network of relationships between genes through a systematic biology approach based on a weighted genome expression network and to identify network-centric genes associated with clinical features of ACC and then use various datasets and databases (GEO, TCGA, GEPIA, and STRING) to validate the value of the hub genes.

Materials and Methods

Data Collection

Expression profiles of mRNA and related clinical data of human ACCs were downloaded from Gene Expression Omnibus (GEO) database1. Dataset GSE10927 (Giordano et al., 2009) performed on Affymetrix Human Genome U133 Plus 2.0 Array was used as a training set to construct co-expression networks and identify hub genes in this study. This dataset included human samples of 33 ACCs, 22 adrenocortical adenomas, and 10 normal adrenal cortex samples, each from a different patient. The additional independent dataset GSE19750 (Demeure et al., 2013) based on the same platform as the training set was downloaded from GEO and used as a test set to verify our results. This dataset included 44 ACC tissues and four normal adrenal glands. Moreover, RNA-sequencing data of 79 ACC samples were also downloaded from The Cancer Genome Atlas (TCGA) database2 to further verify our results. The gene expression data were based on the RNA-sequencing technology IlluminaHiseq.

Data Preprocessing

Normalized data were first downloaded from GEO database. Microarray quality was assessed by sample clustering according to the distance between different samples in Pearson’s correlation matrices (Gautier et al., 2004). A brief study design is shown in Figure 1A.

FIGURE 1
www.frontiersin.org

FIGURE 1. Study design and clustering dendrogram of 32 tumor samples as well as the clinical traits. (A) Flow diagram of data preparation, processing, analysis, and validation in this study. (B) The clustering was based on the expression data of differentially expressed genes between tumor samples and normal samples in ACC. The color intensity was proportional to older age, gender, side, tumor diameter (cm), tumor weight (gm), weiss grade of tumor, mitotic rate of tumor, tumor stage, years to last follow-up and status. Color: Female: white, Male: red; Status: live: white, dead: red.

Differentially Expressed Genes (DEGs) Screening

The “limma” (linear models for microarray data) (Ritchie et al., 2015) R package was used to screen the DEGs between normal adrenal gland and ACC. The SAM (significance analysis of microarrays) with FDR < 0.05 and | log2 FC| > 0.263 were chosen as the cut-off criteria to select genes further considered in the network construction (Feng et al., 2016; Ling et al., 2017; Ma et al., 2017; Wang C.F. et al., 2017; Xiong et al., 2017; Yan et al., 2017; Yang et al., 2017). The differentially expressed genes between normal adrenal cortex and ACC in two dataset GSE10927 and dataset GSE19750 were screened.

Co-expression Network Construction

First, DEG expression data profiles were evaluated to determine whether the samples and genes were of sufficient quality. Then, we used the “WGCNA” package in R to construct scale-free co-expression network for the DEGs. To ensure that the results of network construction were reliable, outlier samples (which were distant from other samples in the clustering via average linkage method) were removed. An appropriate soft threshold power was selected in accordance with standard scale-free networks, with which adjacencies between all differential genes were calculated by a power function (Clarke et al., 2013; Botia et al., 2017). Then, the adjacency was transformed into a TOM, and the corresponding dissimilarity (1-TOM) was calculated. Module identification was accomplished with the dynamic tree cut method by hierarchically clustering genes using 1-TOM as the distance measure with a deepSplit value of two and a minimum size cutoff of 50 for the resulting dendrogram. To further analyze the module, we calculated the dissimilarity of module eigengenes (MEs), chose a cut line for module dendrogram and merged some modules.

Identification of Clinical Significant Modules

Two approaches were used to identify modules related to clinical traits of ACC. First, gene significance (GS) was defined as the log10 transformation of the P value (GS = lgP) in the linear regression between gene expression and the clinical traits. In addition, module significance (MS) was defined as the average GS for all the genes in a module. In general, the module with the absolute MS ranked first or second among all the selected modules was considered as the one related to a clinical trait. MEs were considered as the major component in the principal component analysis for each gene module and the expression patterns of all genes could be summarized into a single characteristic expression profile within a given module. In addition, we calculated the correlation between MEs and clinical traits to identify the relevant module. The module with the maximal absolute MS among all the selected modules was usually considered to be related to a clinical trait. Finally, the module highly correlated with certain clinical traits was selected for further analysis.

Identification of Candidate Hub Genes

Hub genes in the co-expression network were defined by module connectivity, measured by absolute value of the Pearson’s correlation and clinical trait relationship. We identified hub genes in module that were highly correlated with certain clinical traits. Furthermore, we uploaded all genes in the hub module to the Search Tool for the Retrieval of Interacting Genes (STRING) database3 (Szklarczyk et al., 2015, 2017), choosing a confidence > 0.4 to construct a PPI. In the PPI network, genes with a connectivity degree ≥ 4 (node/edge) were defined hub nodes. Hub genes common to both the co-expression network and PPI network were regarded as candidate hub genes for further analysis.

Real Hub Genes Identification

In the GSE19750 test set, linear regression analyses and survival analysis were performed to validate the role of hub genes in the progression and prognosis of ACC. Moreover, we used the Gene Expression Profiling Interactive Analysis (GEPIA) database4 (Tang et al., 2017) to perform survival analysis and stage plots of the candidate hub genes. The ANOVA test was performed using GSE19750 and TCGA-ACC data. Meanwhile, ROC curve analysis was performed and the AUC was calculated to distinguish the tumor tissues from the normal tissues as well as the malignant tumors from the benign tumors by using training set GSE10927 and test set GSE19750. Genes with a significant P value in survival analyses, linear regression analyses, stage plots, common DEGs in GSE19750 and GSE10927, ANOVA test and AUC ≥ 0.85 were defined as the “real” hub genes. Then the real hub genes were performed the Pearson’s correlation between their expression levels and MKi67 by using GSE10927, GSE19750, and TCGA-ACC.

Functional and Pathway Enrichment Analysis

The Database for Annotation, Visualization and Integrate Discovery5 (Dennis et al., 2003; Sherman et al., 2007) is an online program providing a comprehensive set of functional annotation tools for investigators to understand biological meaning behind large list of genes. Enriched biological themes of DEGs in hub modules, particularly GO terms, were identified, and those on KEGG pathway maps were visualized using DAVID. We chose the top 10 terms including the candidate hub genes to be the key BP and pathways. P < 0.05 was set as the cut-off criterion.

Gene Set Enrichment Analysis (GSEA)

In the TCGA-ACC data, 79 ACC samples were divided into two groups according to the median expression level of hub genes. To identify potential function of the hub gene, GSEA6 (Subramanian et al., 2005) was conducted to detect whether a series of a priori defined BPs were enriched in the gene rank derived from DEGs between the two groups.

For use with GSEA software, the collection of annotated gene sets of c2.cp.kegg.v6.0.symbols.gmt in Molecular Signatures Database (MSigDB)7 was chosen as the reference gene set. ES ≥ 0.6 and FDR < 0.05 were chosen as the cut-off criteria. We chose the terms enriched in all hub genes as the potential function of the hub genes.

Results

DEG Screening

After data preprocessing and quality assessment, the expression matrices were obtained from the 65 samples in training set GSE10927. Under the threshold of FDR < 0.05 and | log2 FC| > 0.263, a total of 1956 DEGs (666 up-regulated and 1290 down-regulated in ACC samples) were selected for subsequent analysis (Supplementary Table S1). In the test set GSE19750, a total of 3810 DEGs (1263 up-regulated and 2547 in ACC samples) were identified (Supplementary Table S2).

Weighted Co-expression Network Construction and Key Modules Identification

GSM277130 was an outlier sample and was removed from subsequent analysis in GSE10927 (Supplementary Figure S1). A total of 33 samples with clinical data were included in the co-expression analysis (Figure 1B). Using the “WGCNA” package in R, DEGs with similar expression patterns were grouped into modules via the average linkage hierarchical clustering. In this study, the power of β = 3 (scale free R2 = 0.89) was selected as the soft-thresholding to ensure a scale-free network (Figure 2). A total of four modules were identified (Figure 3A). Two methods were used to test the relevance between each module and the ACC clinical information; here, we focused on the tumor grade. First, modules with a greater MS were considered to have more connection with tumor grade. However, most of the correlations were low to moderate (R2 < 0.5), and we found that the MS of the turquoise module (P = 9 × 10-5, R2 = 0.64) was higher than that of any other module (Figure 3B). Afterward, the ME of the turquoise module showed a higher correlation with tumor grade than did other modules (Figure 3C). Based on the two methods, the turquoise module with tumor grade was identified as the clinical significant module, which was extracted for further analysis.

FIGURE 2
www.frontiersin.org

FIGURE 2. Determination of soft-thresholding power in the weighted gene co-expression network analysis (WGCNA). (A) Analysis of the scale-free fit index for various soft-thresholding powers (β). (B) Analysis of the mean connectivity for various soft-thresholding powers. (C) Histogram of connectivity distribution when β = 3. (D) Checking the scale free topology when β = 3.

FIGURE 3
www.frontiersin.org

FIGURE 3. Identification of modules associated with clinical information. (A) Dendrogram of all differentially expressed genes clustered based on a dissimilarity measure (1-TOM). (B) Heatmap of the correlation between module eigengenes (MEs) and different clinical information of ACC (age, gender, side, tumor diameter (cm), tumor weight (gm), weiss grade of tumor, mitotic rate of tumor, tumor stage, years to last follow-up and status). (C) Distribution of average gene significances and errors in the modules associated with the weiss grade of ACC.

Identification of Candidate Hub Genes for Tumor Grade

Highly connected hub genes in a module play important roles in the BPs. As to the PPI network, under the cutoff of confidence > 0.4 and connectivity degree of ≥ 4 (node/edge), 123 genes were identified as hub genes (Figure 4A) and the whole PPI network were showed in Supplementary Figure S2. Defined by module connectivity and measured by absolute value of the Pearson’s correlation and clinical trait relationship, 78 genes with high connectivity in the turquoise module were taken as candidate hub genes in the co-expression network (Figure 4B). A total of 60 genes were identified both in PPI network and co-expression network identified as candidate hub genes (Figure 4C).

FIGURE 4
www.frontiersin.org

FIGURE 4. Hub genes detection. (A) Protein–protein interaction network of the hub nodes. (B) Scatter plot of MEs in turquoise module. (C) Common hub genes in the co-expression network and PPI network. (C) Common genes with significant p value in survival analysis, ROC curve analysis, one-way ANOVA, Pearson’s correlation analysis and DEGs analysis.

Real Hub Gene Identification

Linear regression analyses and one-way ANOVA analysis were conducted to validate candidate hub genes in the test set GSE19750 and TCGA-ACC (Supplementary Table S3), and totally 46 and 41 genes were identified, respectively. As tumor progression always affects tumor prognosis, we validated the real hub genes by investigating their roles in ACC prognosis using GSE19750 and TCGA-ACC (Supplementary Table S4) and 30 genes showed significant P value in both test sets. To further test the diagnostic capacity of the candidate hub genes, ROC curve was performed and the AUC were calculated (Supplementary Table S5), and 34 genes were identified with AUC ≥ 0.85. In addition, we also investigated the expression levels of candidate hub genes in two datasets (GSE19750 and GSE10927), and 31 were differentially expressed in both datasets (GSE19750 and GSE10927). Common genes with significant P value in those five analyses were screened as the real hub genes, and 12 genes were eventually identified (ANLN, ASPM, CDCA5, CENPF, FOXM1, KIAA0101, MELK, NDC80, PRC1, RACGAP1, SPAG5, and TPX2) (Figure 4D). The grade or stage plot, Pearson’s correlation, ANOVA test results and survival analyses of the real hub genes were showed in Figures 57 and Supplementary Figures S3S5. MKi67 is a biomarker for tumor proliferation, the Pearson’s correlation of MKi67 and real hub genes were performed (Figure 8).

FIGURE 5
www.frontiersin.org

FIGURE 5. Stage and grade plot of the real hub genes. (A) Grade plot of the real hub genes using test set GSE19750. (B) Stage plot of the real hub genes using test set TCGA-ACC. (C) Pearson’s correlation analysis of the real hub genes using test set GSE19750 and TCGA-ACC. (D) One-way ANOVA analysis of the real hub genes using test set GSE19750 and TCGA-ACC.

FIGURE 6
www.frontiersin.org

FIGURE 6. Survival analysis of the association between the expression levels of real hub genes and overall survival time in ACC (based on TCGA data in GEPIA). (A) ANLN, (B) ASPM, (C) CCNA5, (D) CENPF, (E) FOXM1, (F) KIAA0101, (G) MELK, (H) NDC80, (I) PRC1, (J) RACGAP1, (K) SPAG5, (L) TPX2. Red line indicates the samples with gene highly expressed, and blue line shows the samples with gene lowly expressed. HR, hazard ratio.

FIGURE 7
www.frontiersin.org

FIGURE 7. Survival analysis of the association between the expression levels of real hub genes and disease free survival time in ACC (based on TCGA data in GEPIA). (A) ANLN, (B) ASPM, (C) CCNA5, (D) CENPF, (E) FOXM1, (F) KIAA0101, (G) MELK, (H) NDC80, (I) PRC1, (J) RACGAP1, (K) SPAG5, (L) TPX2. Red line indicates the samples with gene highly expressed, and blue line shows the samples with gene lowly expressed. HR, hazard ratio.

FIGURE 8
www.frontiersin.org

FIGURE 8. Pearson’s correlation analysis of the association between the Ki67 expression and the expression of real hub genes.

Functional and Pathway Enrichment Analysis

To obtain further insight into the function of DEGs in hub module, the DEGs were uploaded to the DAVID database. GO analysis results showed that hub genes were enriched in the top 10 BPs, including cell cycle phase, cell cycle, M phase, nuclear division, mitosis, organelle fission, M phase of mitotic cell cycle, cell cycle process, mitotic cell cycle and cell division. Moreover, hub genes were quite significantly enriched in cell cycle (Figure 9).

FIGURE 9
www.frontiersin.org

FIGURE 9. Bioinformatics analysis of genes in turquoise module. (A) GO analysis and (B) KEGG pathway enrichment of hub genes.

Gene Set Enrichment Analysis

To identify the potential function of those hub genes in ACC, GSEA was conducted to search BPs enriched in any hub gene highly expressed samples. Seven gene sets were enriched in the samples with all hub genes highly expressed, including “base excision repair,” “cell cycle,” “DNA replication,” “mismatch repair,” “nucleotide excision repair,” “spliceosome,” and “homologous recombination” (Supplementary Figure S6).

Discussion

Adrenocortical carcinoma is one of the most aggressive solid tumors in humans, with an overall poor prognosis (Fassnacht et al., 2011; Else et al., 2014). Successful surgical resection for early stage disease is the only curative treatment for ACC. However, the risk of recurrence is high. Therefore, it is particularly important to stratify patients with ACC into low- or high-risk groups to monitor disease recurrence and assign them to appropriate therapeutic interventions. Recently, a large multicentric study by the European Network for the Study of Adrenal Tumors has demonstrated that the Ki67 labeling index (Li) is the most powerful parameter predicting disease recurrence and survival in ACC patents after complete tumor resection (Beuschlein et al., 2015). Recent studies have attached values to biomarkers for outcomes of ACC; however, carcinogenesis does not occur due to a few altered genes. Thus, it is necessary to identify novel molecular biomarkers that can predict disease stage and clinical outcome of ACC patients; these biomarkers could play critical roles in pathogenesis and enable personalized treatment.

Weighted genes co-expression network analysis was used to construct gene co-expression networks on the basis of similarities of expression profiles among samples and then illustrate a global interpretation of genes. Many articles have been published using WGCNA to screen clinical-feature related biomarkers in several types of tumors (Chen P. et al., 2017; Wang F. et al., 2017; Yuan et al., 2017). This study is the first to use the WGCNA algorithm for the purpose of identifying potential biomarkers for ACC. Based on the clinical features (age, gender, tumor side, tumor diameter, tumor weight, Weiss grade of tumor, mitotic rate of tumor, tumor stage, years to last follow, and living status), we eventually screened the gene co-expression modules related to the tumor grade of ACC. The turquoise module was identified, and 78 hub genes were derived from the module in co-expression network; meanwhile, 123 hub genes were identified in PPI analysis. Relating the results of PPI network, 60 hub nodes were screened in both the co-expression module and PPI network. Further regression analysis, survival analysis, one-way ANOVA and ROC curve analysis were performed, and 12 real hub genes were finally selected that were associated with poor outcomes and tumor stage of ACC. Interestingly, we found that all real hub genes differed greatly between tumor, begin or normal tissue in the test set and GEPIA database. Additionally, among the 12 real hub genes, only few biomarkers were identified in previous studies. Jain et al. (2011) discovered that KIAA0101, which is a marker of cellular proliferation and promotes growth and invasion, was a good diagnostic marker for distinguishing benign from malignant adrenocortical neoplasm. Although other hub genes were not reported to participate in ACC progression, they were reported to be involved in various tumors. TPX2, overexpressed in cancers, is being established as marker for the diagnosis and prognosis of malignancies, and the functions of TPX2 are attributed to its Ran-regulated microtubule-associated protein properties and to its control of the Aurora A kinase (Neumayer et al., 2014). SPAG5 is reported as a poor prognostic biomarker in cervix cancer, breast cancer and bladder cancer (Yuan et al., 2014; Abdel-Fatah et al., 2016; Liu et al., 2018). RACGAP1 as a proliferation marker was indicated to involve in carcinogenesis of many cancers (Ke et al., 2013; Milde-Langosch et al., 2013; Imaoka et al., 2015; Bornschein et al., 2016; Sahin et al., 2016). In recent study, researchers found that abnormal PRC1 expression correlates with poor patient outcome in various malignancies, which may be caused by PRC1-mediated CIN and aneuploidy (Li et al., 2018). NDC80 was reported to associate with poor outcomes and tumor proliferation in multiple cancers via cell cycle regulation (Qu et al., 2014; Meng et al., 2015; Xing et al., 2016). In previous studies, MELK functioned as a modulator of intracellular signaling and affects various cellular and BPs, including cell cycle, cell proliferation, apoptosis, spliceosome assembly, gene expression, embryonic development, hematopoiesis, and oncogenesis (Jiang and Zhang, 2013). FOXM1 transcription factor is a regulator of myriad BPs, including cell proliferation, cell cycle progression, cell differentiation, DNA damage repair, tissue homeostasis, angiogenesis and apoptosis; meanwhile, elevated FOXM1 expression was found in most cancers, suggesting it has an integral role in tumorigenesis and recent research findings also placed FOXM1 at the center of cancer progression and drug sensitivity (Koo et al., 2012). Interestingly, co-expression of FOXM1 and CENPF is a robust prognostic indicator of poor survival and metastasis (Aytes et al., 2014). CDCA5 was identified as a negative prognostic marker in bladder cancer, breast cancer, gastric cancer and hepatocellular carcinoma (Chang et al., 2015; Phan et al., 2018; Shen et al., 2018; Zhang et al., 2018). ASPM played an important role in pancreatic tumor and glioma progression, in vascular invasion, early recurrence and poor prognosis of hepatocellular carcinoma and even represented a potential target for radiotherapy (Lin et al., 2008; Bikeye et al., 2010; Kato et al., 2011; Wang et al., 2013). According to the previous studies, we found all the hub genes participated in the tumorigenesis.

Enrichment analyses for the turquoise module indicated that BPs focused on cell cycle regulation, including cell cycle phase, cell cycle, M phase, nuclear division, mitosis, organelle fission, M phase of mitotic cell cycle, cell cycle process, mitotic cell cycle and cell division. Meanwhile, KEGG enrichment analysis revealed that genes in the turquoise module were significantly enriched in cell cycle. To further study the mechanism by which these genes regulate tumorigenesis, we performed GSEA analysis using ACC RNA Seq data from TCGA. All hub genes were enriched in base excision repair, cell cycle, DNA replication, homologous recombination, mismatch repair, nucleotide excision repair and spliceosome. The two enrichment analyses demonstrated that hub genes were closely related to the function of cell cycle regulation. Previous studies demonstrated that the presence of p53 was necessary for DNA-damaged cells to arrest, repair the damage, and reenter the cell cycle (Lukin et al., 2015). Cell cycle is a continuous and accurate process; meanwhile, monitoring and regulating the cell cycle are important (Xu, 2006). Deregulation of cell cycle is closely related to carcinogenesis and tumor progression (Kaistha et al., 2016). Regarding our candidate hub genes that were enriched in the function of regulating cell cycle, we hypothesize that they might influence ACC by affecting cell cycle.

Conclusion

In conclusion, we used weighted gene co-expression analysis to construct a gene co-expression network for ACC for the first time, as well as identified and validated network hub genes associated with tumor grade. Interestingly, candidate hub genes showed significant prognostic value. Meanwhile, we predicted the potential function of the stage- and prognosis-related hub genes, which participated in cell cycle regulation. These hub genes had potential roles in the translational medicine and might become therapeutic targets. However, our study only selected potential biomarkers of ACC that were associated with tumor progression and prognosis, and further molecular biological experiments are needed to confirm the function of candidate biomarkers in ACC.

Author Contributions

LY, YX, and XW conceived and designed the study. LY, GQ, LC, and YX performed the analysis procedures. LY, GQ, C-LW, HD, and YX analyzed the results. LY, GQ, and XW contributed analysis tools. LY, GQ, and YX contributed to the writing of the manuscript. All authors reviewed the manuscript.

Funding

This study was supported in part by grants from the Hubei Province Health and Family Planning Scientific Research Project (Grant No. WJ2017H0002) and National Natural Science Foundation of China (Grant No. 81772730). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The excellent technical assistance of Shanshan Zhang and Danni Shan is gratefully acknowledged. We would like to acknowledge the KEGG database developed by Kanehisa Laboratories. We would like to acknowledge the Oncomine, OncoLnc, and TCGA databases for free use. We would also like to acknowledge the KEGG database developed by Kanehisa Laboratories for the original source of KEGG pathway images.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2018.00328/full#supplementary-material

Abbreviations

ACC, adrenocortical carcinoma; BP, biological process; DAVID, the database for annotation, visualization, and integrated discovery; FC, fold change; FDR, false discovery rate; GO, gene ontology; GSEA, gene set enrichment analysis; PPI, protein–protein interaction; TOM, topological overlap matrix; WGCNA, weighted gene co-expression network analysis.

Footnotes

  1. ^http://www.ncbi.nlm.nih.gov/geo/
  2. ^https://genome-cancer.ucsc.edu/
  3. ^http://www.string-db.org/
  4. ^http://gepia.cancer-pku.cn/
  5. ^http://david.abcc.ncifcrf.gov/
  6. ^http://software.broadinstitute.org/gsea/index.jsp
  7. ^http://software.broadinstitute.org/gsea/msigdb/index.jsp

References

Abdel-Fatah, T. M. A., Agarwal, D., Liu, D. X., Russell, R., Rueda, O. M., Liu, K., et al. (2016). SPAG5 as a prognostic biomarker and chemotherapy sensitivity predictor in breast cancer: a retrospective, integrated genomic, transcriptomic, and protein analysis. Lancet Oncol. 17, 1004–1018. doi: 10.1016/S1470-2045(16)00174-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayala-Ramirez, M., Jasim, S., Feng, L., Ejaz, S., Deniz, F., Busaidy, N., et al. (2013). Adrenocortical carcinoma: clinical outcomes and prognosis of 330 patients at a tertiary care center. Eur. J. Endocrinol. 169, 891–899. doi: 10.1530/EJE-13-0519

PubMed Abstract | CrossRef Full Text | Google Scholar

Aytes, A., Mitrofanova, A., Lefebvre, C., Alvarez, M. J., Castillo-Martin, M., Zheng, T., et al. (2014). Cross-species regulatory network analysis identifies a synergistic interaction between FOXM1 and CENPF that drives prostate cancer malignancy. Cancer Cell 25, 638–651. doi: 10.1016/j.ccr.2014.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Beuschlein, F., Weigel, J., Saeger, W., Kroiss, M., Wild, V., Daffara, F., et al. (2015). Major prognostic role of Ki67 in localized adrenocortical carcinoma after complete resection. J. Clin. Endocrinol. Metab. 100, 841–849.

Google Scholar

Bikeye, S. N., Colin, C., Marie, Y., Vampouille, R., Ravassard, P., Rousseau, A., et al. (2010). ASPM-associated stem cell proliferation is involved in malignant progression of gliomas and constitutes an attractive therapeutic target. Cancer Cell Int. 10:1. doi: 10.1186/1475-2867-10-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Bornschein, J., Nielitz, J., Drozdov, I., Selgrad, M., Wex, T., Jechorek, D., et al. (2016). Expression of aurora kinase A correlates with the Wnt-modulator RACGAP1 in gastric cancer. Cancer Med. 5, 516–526. doi: 10.1002/cam4.610

PubMed Abstract | CrossRef Full Text | Google Scholar

Botia, J. A., Vandrovcova, J., Forabosco, P., Guelfi, S., D’Sa, K., United Kingdom, et al. (2017). An additional k-means clustering step improves the biological features of WGCNA gene co-expression networks. BMC Syst. Biol. 11:47. doi: 10.1186/s12918-017-0420-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, I. W., Lin, V. C., He, H. L., Hsu, C. T., Li, C. C., Wu, W. J., et al. (2015). CDCA5 overexpression is an indicator of poor prognosis in patients with urothelial carcinomas of the upper urinary tract and urinary bladder. Am. J. Transl. Res. 7, 710–722.

PubMed Abstract | Google Scholar

Chen, L., Yuan, L., Wang, G., Cao, R., Peng, J., Shu, B., et al. (2017). Identification and bioinformatics analysis of miRNAs associated with human muscle invasive bladder cancer. Mol. Med. Rep. 16, 8709–8720. doi: 10.3892/mmr.2017.7726

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, P., Wang, F., Feng, J., Zhou, R., Chang, Y., Liu, J., et al. (2017). Co-expression network analysis identified six hub genes in association with metastasis risk and prognosis in hepatocellular carcinoma. Oncotarget 8, 48948–48958.

PubMed Abstract | Google Scholar

Clarke, C., Madden, S. F., Doolan, P., Aherne, S. T., Joyce, H., O’Driscoll, L., et al. (2013). Correlating transcriptional networks to breast cancer survival: a large-scale coexpression analysis. Carcinogenesis 34, 2300–2308. doi: 10.1093/carcin/bgt208

PubMed Abstract | CrossRef Full Text | Google Scholar

Demeure, M. J., Coan, K. E., Grant, C. S., Komorowski, R. A., Stephan, E., Sinari, S., et al. (2013). PTTG1 overexpression in adrenocortical cancer is associated with poor survival and represents a potential therapeutic target. Surgery 154, 1405–1416. doi: 10.1016/j.surg.2013.06.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Dennis, G. Jr., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., et al. (2003). DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 4:P3.

Google Scholar

Else, T., Kim, A. C., Sabolch, A., Raymond, V. M., Kandathil, A., Caoili, E. M., et al. (2014). Adrenocortical carcinoma. Endocr. Rev. 35, 282–326. doi: 10.1210/er.2013-1029

PubMed Abstract | CrossRef Full Text | Google Scholar

Fassnacht, M., Libe, R., Kroiss, M., and Allolio, B. (2011). Adrenocortical carcinoma: a clinician’s update. Nat. Rev. Endocrinol. 7, 323–335. doi: 10.1038/nrendo.2010.235

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, A., Tu, Z., and Yin, B. (2016). The effect of HMGB1 on the clinicopathological and prognostic features of non-small cell lung cancer. Oncotarget 7, 20507–20519. doi: 10.18632/oncotarget.7050

PubMed Abstract | CrossRef Full Text | Google Scholar

Gautier, L., Cope, L., Bolstad, B. M., and Irizarry, R. A. (2004). affy–analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 20, 307–315. doi: 10.1093/bioinformatics/btg405

PubMed Abstract | CrossRef Full Text | Google Scholar

Giordano, T. J., Kuick, R., Else, T., Gauger, P. G., Vinco, M., Bauersfeld, J., et al. (2009). Molecular classification and prognostication of adrenocortical tumors by transcriptome profiling. Clin. Cancer Res. 15, 668–676. doi: 10.1158/1078-0432.CCR-08-1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Imaoka, H., Toiyama, Y., Saigusa, S., Kawamura, M., Kawamoto, A., Okugawa, Y., et al. (2015). RacGAP1 expression, increasing tumor malignant potential, as a predictive biomarker for lymph node metastasis and poor prognosis in colorectal cancer. Carcinogenesis 36, 346–354. doi: 10.1093/carcin/bgu327

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, M., Zhang, L., Patterson, E. E., and Kebebew, E. (2011). KIAA0101 is overexpressed, and promotes growth and invasion in adrenal cancer. PLoS One 6:e26866. doi: 10.1371/journal.pone.0026866

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., and Zhang, D. (2013). Maternal embryonic leucine zipper kinase (MELK): a novel regulator in cell cycle control, embryonic development, and cancer. Int. J. Mol. Sci. 14, 21551–21560. doi: 10.3390/ijms141121551

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaistha, B. P., Lorenz, H., Schmidt, H., Sipos, B., Pawlak, M., Gierke, B., et al. (2016). PLAC8 localizes to the inner plasma membrane of pancreatic cancer cells and regulates cell growth and disease progression through critical cell-cycle regulatory pathways. Cancer Res. 76, 96–107. doi: 10.1158/0008-5472.CAN-15-0216

PubMed Abstract | CrossRef Full Text | Google Scholar

Kato, T. A., Okayasu, R., Jeggo, P. A., and Fujimori, A. (2011). ASPM influences DNA double-strand break repair and represents a potential target for radiotherapy. Int. J. Radiat. Biol. 87, 1189–1195. doi: 10.3109/09553002.2011.624152

PubMed Abstract | CrossRef Full Text | Google Scholar

Ke, H. L., Ke, R. H., Li, S. T., Li, B., Lu, H. T., and Wang, X. Q. (2013). Expression of RACGAP1 in high grade meningiomas: a potential role in cancer progression. J. Neurooncol. 113, 327–332. doi: 10.1007/s11060-013-1121-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Koo, C. Y., Muir, K. W., and Lam, E. W. (2012). FOXM1: from cancer initiation to progression and treatment. Biochim. Biophys. Acta 1819, 28–37. doi: 10.1016/j.bbagrm.2011.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kulshrestha, A., Suman, S., and Ranjan, R. (2016). Network analysis reveals potential markers for pediatric adrenocortical carcinoma. Onco Targets Ther. 9, 4569–4581. doi: 10.2147/OTT.S108485

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Dallmayer, M., Kirchner, T., Musa, J., and Grunewald, T. G. P. (2018). PRC1: linking cytokinesis, chromosomal instability, and cancer evolution. Trends Cancer 4, 59–73. doi: 10.1016/j.trecan.2017.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, S. Y., Pan, H. W., Liu, S. H., Jeng, Y. M., Hu, F. C., Peng, S. Y., et al. (2008). ASPM is a novel marker for vascular invasion, early recurrence, and poor prognosis of hepatocellular carcinoma. Clin. Cancer Res. 14, 4814–4820. doi: 10.1158/1078-0432.CCR-07-5262

PubMed Abstract | CrossRef Full Text | Google Scholar

Ling, Q., Xie, H., Li, J., Liu, J., Cao, J., Yang, F., et al. (2017). Donor graft MicroRNAs: a newly identified player in the development of new-onset diabetes after liver transplantation. Am. J. Transplant. 17, 255–264. doi: 10.1111/ajt.13984

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J. Y., Zeng, Q. H., Cao, P. G., Xie, D., Yang, F., He, L. Y., et al. (2018). SPAG5 promotes proliferation and suppresses apoptosis in bladder urothelial carcinoma by upregulating Wnt3 via activating the AKT/mTOR pathway and predicts poorer survival. Oncogene. doi: 10.1038/s41388-018-0223-2 [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Lukin, D. J., Carvajal, L. A., Liu, W. J., Resnick-Silverman, L., and Manfredi, J. J. (2015). p53 Promotes cell survival due to the reversibility of its cell-cycle checkpoints. Mol. Cancer Res. 13, 16–28. doi: 10.1158/1541-7786.MCR-14-0177

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J., Kang, M., Zhang, Y., Guo, X., Tian, Z., Ding, C., et al. (2017). Self-defense of Escherichia coli against damages caused by nanoalumina. Environ. Toxicol. Pharmacol. 55, 110–117. doi: 10.1016/j.etap.2017.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, Q. C., Wang, H. C., Song, Z. L., Shan, Z. Z., Yuan, Z., Zheng, Q., et al. (2015). Overexpression of NDC80 is correlated with prognosis of pancreatic cancer and regulates cell proliferation. Am. J. Cancer Res. 5, 1730–1740.

PubMed Abstract | Google Scholar

Milde-Langosch, K., Karn, T., Muller, V., Witzel, I., Rody, A., Schmidt, M., et al. (2013). Validity of the proliferation markers Ki67, TOP2A, and RacGAP1 in molecular subgroups of breast cancer. Breast Cancer Res. Treat. 137, 57–67. doi: 10.1007/s10549-012-2296-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Neumayer, G., Belzil, C., Gruss, O. J., and Nguyen, M. D. (2014). TPX2: of spindle assembly, DNA damage response, and cancer. Cell. Mol. Life. Sci. 71, 3027–3047. doi: 10.1007/s00018-014-1582-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Phan, N. N., Wang, C. Y., Li, K. L., Chen, C. F., Chiao, C. C., Yu, H. G., et al. (2018). Distinct expression of CDCA 3, CDCA 5, and CDCA8 leads to shorter relapse free survival in breast cancer patient. Oncotarget 9, 6977–6992. doi: 10.18632/oncotarget.24059

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, E. M., Chen, X., Easton, J., Finkelstein, D., Liu, Z., Pounds, S., et al. (2015). Genomic landscape of paediatric adrenocortical tumours. Nat. Commun. 6:6302. doi: 10.1038/ncomms7302

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, Y., Li, J., Cai, Q., and Liu, B. (2014). Hec1/Ndc80 is overexpressed in human gastric cancer and regulates cell growth. J. Gastroenterol. 49, 408–418. doi: 10.1007/s00535-013-0809-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, A. G., Kim, J., Al-Ahmadie, H., Bellmunt, J., Guo, G., Cherniack, A. D., et al. (2017). Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell 171, 540.e25–556.e25. doi: 10.1016/j.cell.2017.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahin, S., Isik Gonul, I., Cakir, A., Seckin, S., and Uluoglu, O. (2016). Clinicopathological significance of the proliferation markers Ki67, RacGAP 1, and topoisomerase 2 alpha in breast cancer. Int. J. Surg. Pathol. 24, 607–613. doi: 10.1177/1066896916653211

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Z., Yu, X., Zheng, Y., Lai, X., Li, J., Hong, Y., et al. (2018). CDCA5 regulates proliferation in hepatocellular carcinoma and has potential as a negative prognostic marker. Onco Targets Ther. 11, 891–901. doi: 10.2147/OTT.S154754

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, B. T., Huang, da, W., Tan, Q., Guo, Y., Bour, S., et al. (2007). DAVID Knowledgebase: a gene-centered database integrating heterogeneous gene annotation resources to facilitate high-throughput gene functional analysis. BMC Bioinformatics 8:426. doi: 10.1186/1471-2105-8-426

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 45, D362–D368. doi: 10.1093/nar/gkw937

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017). GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 45, W98–W102.

PubMed Abstract | Google Scholar

Wang, C. F., Zhao, C. C., Weng, W. J., Lei, J., Lin, Y., Mao, Q., et al. (2017). Alteration in long non-coding RNA expression after traumatic brain injury in rats. J. Neurotrauma 34, 2100–2108. doi: 10.1089/neu.2016.4642

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Chang, Y., Li, J., Wang, H., Zhou, R., Qi, J., et al. (2017). Strong correlation between ASPM gene expression and HCV cirrhosis progression identified by co-expression analysis. Dig. Liver Dis. 49, 70–76. doi: 10.1016/j.dld.2016.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W. Y., Hsu, C. C., Wang, T. Y., Li, C. R., Hou, Y. C., Chu, J. M., et al. (2013). A gene expression signature of epithelial tubulogenesis and a role for ASPM in pancreatic tumor progression. Gastroenterology 145, 1110–1120. doi: 10.1053/j.gastro.2013.07.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing, X. K., Wu, H. Y., Chen, H. L., and Feng, H. G. (2016). NDC80 promotes proliferation and metastasis of colon cancer cells. Genet. Mol. Res. 15:gmr8312. doi: 10.4238/gmr.15028312

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, K., Long, L., Zhang, X., Qu, H., Deng, H., Ding, Y., et al. (2017). Overview of long non-coding RNA and mRNA expression in response to methamphetamine treatment in vitro. Toxicol. In Vitro 44, 1–10. doi: 10.1016/j.tiv.2017.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Y. (2006). DNA damage: a trigger of innate immunity but a requirement for adaptive immune homeostasis. Nat. Rev. Immunol. 6, 261–270. doi: 10.1038/nri1804

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, S., Zhang, H., Xie, W., Meng, F., Zhang, K., Jiang, Y., et al. (2017). Altered microRNA profiles in plasma exosomes from mesial temporal lobe epilepsy with hippocampal sclerosis. Oncotarget 8, 4136–4146. doi: 10.18632/oncotarget.13744

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, M., Chen, B. L., Huang, J. B., Meng, Y. N., Duan, X. J., Chen, L., et al. (2017). Angiogenesis-related genes may be a more important factor than matrix metalloproteinases in bronchopulmonary dysplasia development. Oncotarget 8, 18670–18679. doi: 10.18632/oncotarget.14722

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, L., Shu, B., Chen, L., Qian, K., Wang, Y., Qian, G., et al. (2017). Overexpression of COL3A1 confers a poor prognosis in human bladder cancer identified by co-expression analysis. Oncotarget 8, 70508–70520. doi: 10.18632/oncotarget.19733

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, L. J., Li, J. D., Zhang, L., Wang, J. H., Wan, T., Zhou, Y., et al. (2014). SPAG5 upregulation predicts poor prognosis in cervical cancer patients and alters sensitivity to taxol treatment via the mTOR signaling pathway. Cell Death Dis. 5:e1247. doi: 10.1038/cddis.2014.222

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Shen, M., and Zhou, G. (2018). Upregulation of CDCA5 promotes gastric cancer malignant progression via influencing cyclin E1. Biochem. Biophys. Res. Commun. 496, 482–489. doi: 10.1016/j.bbrc.2018.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: adrenocortical carcinoma (ACC), weighted gene co-expression network analysis (WGCNA), cell cycle, biomarker, progression and prognosis

Citation: Yuan L, Qian G, Chen L, Wu C-L, Dan HC, Xiao Y and Wang X (2018) Co-expression Network Analysis of Biomarkers for Adrenocortical Carcinoma. Front. Genet. 9:328. doi: 10.3389/fgene.2018.00328

Received: 05 February 2018; Accepted: 31 July 2018;
Published: 15 August 2018.

Edited by:

Louise Metherell, Queen Mary University of London, United Kingdom

Reviewed by:

Laura Ellestad, University of Georgia, United States
Vincenzo Pezzi, University of Calabria, Italy

Copyright © 2018 Yuan, Qian, Chen, Wu, Dan, Xiao and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Yu Xiao, yu.xiao@whu.edu.cn Xinghuan Wang, wangxinghuan@whu.edu.cn

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.