This article is mentioned in:

Introduction

Cervical cancer (CC) is the fourth most common type of female malignancy with a high mortality rate worldwide, and almost all cases of CC are caused by human papillomavirus (HPV) infection (1). In China, HPV infection is the most common etiological factor underlying CC (2). Our previous study showed that the prevalence of HPV in the Jiangxi Province is 22.49% (3). However, although HPV infection is a prerequisite of CC, it is not sufficient to induce CC alone (4). The specific molecular mechanisms between the persistent high-risk HPV infection and the pathological process of CC are still controversial (5). Increasing evidence has shown that the abnormal expression of multiple genes are involved in the pathogenesis of CC (6). As tumorigenesis is a complex pathological process involving various genetic and epigenetic events, such as the inactivation of suppressor genes and/or the overexpression of oncogenes (7), the identification of dysregulated genes in cancer-associated pathways may highlight the molecular mechanisms underlying tumorigenesis, assisting in the development of novel strategies for treatment of CC. Therefore, it is important to elucidate the potential molecular mechanisms underlying CC, to provide novel therapeutic targets and prognostic biomarkers of CC (8,9).

To improve our understanding of the molecular mechanisms underlying CC, bioinformatics analysis has been widely used to identify differentially expressed genes (DEGs), functional pathways and hub genes involved in the carcinogenesis and progression of CC. Although several genes have been identified in predicting the clinical outcome of CC, there have been inconsistencies between previous studies (10,11), which may be due to small sample sizes, heterogeneous histological subtypes, different detection platforms and various data processing methods. In the present study, three mRNA microarray datasets were downloaded from the Gene Expression Omnibus (GEO) database to acquire DEGs between CC and normal tissues. Subsequently, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, protein-protein interaction (PPI) network analysis and the Cytoscape plugin, CytoHubba, were used to identify hub genes in CC. Furthermore, Kaplan-Meier plotter was used on the cBioPortal website and the DNA topoisomerase II α (TOP2A) gene was identified as an independent prognostic biomarker of disease-free survival in CC. A total of 188 DEGs and 8 hub genes were identified, which may be candidate biomarkers for CC.

Materials and methods

Microarray data information

GSE63514(12), GSE9750(13) and GSE7410(14) were obtained from the NCBI GEO database (ncbi.nlm.nih.gov/geo) (15,16). The GSE63514 dataset used the GPL570 [HG-U133 Plus 2] Affymetrix Human Genome U133 Plus 2.0 Array platform, which contained 28 CC tissue samples and 24 normal samples. The GSE9750 dataset used the GPL96 Affymetrix Human Genome U133A Array platform and included 33 CC tissue samples that were primarily marked by HPV16 or HPV18 and 21 normal cervical samples. GSE7410 contained 40 CC and 5 normal cervical samples, and used the GPL1708 Agilent-012391 Whole Human Genome Oligo Microarray G4112A platform.

Screening of DEGs in CC

The DEGs between CC and normal cervical samples were screened using GEO2R (ncbi.nlm.nih.gov/geo/geo2r). The adjusted P-values (adj.P) and Benjamini and Hochberg false discovery rates were used to provide a balance between the discovery of statistically significant genes and limit false-positives. Probe sets without corresponding gene symbols or genes with >1 probe set were removed or averaged, respectively. Only genes with a log fold-change >1 and adj.P-value of <0.01 were regarded as statistically significant DEGs.

Functional annotation and pathway enrichment analysis

The screened DEGs were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID; david.abcc.ncifcrf.gov/; version 6.7) (17), an online tool for functional annotation. KEGG is a database for understanding higher-level functions and biological information generated by high throughput experimental technologies (18). GO analysis is a useful bioinformatics tool for annotating genes and analyzing their biological processes (19). The significant enrichment analysis of DEGs was evaluated using GO and KEGG, with P<0.05 as the cutoff. If there were >5 terms enriched in the category, the top 10 were selected based on the P-value.

PPI network construction and module analysis

The identified DEGs were imported into the Search Tool for the Retrieval of Interacting Genes (STRING; string-db.org; version 10.0) (20) database to evaluate the functional interactions among them. Subsequently, Cytoscape (cytoscape.org/; version 3.6.1) was used to screen the modules from the PPI network (21), as a molecular interaction network tool. The plugins Cytoscape and CytoHubba were used to predict and explore the important nodes and subnetworks in the network with 12 topological algorithms, which included degree, edge percolated component, maximum neighborhood component and density of maximum neighborhood component (22). CytoHubba was used to rank nodes in a network by their network features and select the overlapping top 15 hub genes by 6 ranked methods. Finally, all identified hub genes in this module were mapped to DAVID to perform KEGG and GO enrichment analysis.

Key gene selection and analysis

Overall and disease-free survival analysis of hub genes using Kaplan-Meier survival curves, as well as correlation analysis of TOP2A and co-expression genes, were performed in cBioPortal (cbioportal.org; version 2.4.3) (23,24). The expression profile of TOP2A was analyzed and displayed using the online database UALCAN (ualcan.path.uab.edu). The relationship between TOP2A and expression patterns and HPV viral infection status were analyzed using the 201291_s_at and 201292_probes in the Zhai Cervix dataset available from Oncomine (oncomine.com) (13,14,25). The heat maps of TOP2A gene expression in clinical CC samples vs. normal tissues are respectively: Cervical Squamous Cell Carcinoma vs. Normal Biewenga Cervix, Gynecol Oncol, 2008(14); Cancer Type: CC Bittner Multi-cancer, unpublished data, 2006; Cervical Squamous Cell Carcinoma vs. Normal Scotto Cervix 2, Genes Chromosomes Cancer, 2008(15); and Cervical Squamous Cell Carcinoma Epithelia vs. Normal Zhai Cervix, Cancer Res, 2007(25).

Statistical analysis

All statistical analyses were performed using the aforementioned bioinformatics tools. Adj.P were corrected for in multiple comparisons using the Benjamini and Hochberg's false discovery rate. Co-expression and networks were calculated according to the cBioPortal online instructions. A log-rank test was performed to identify the significance of the Spearman's Correlation coefficient between the mRNA expression z-Scores (RNASeq V2 RSEM). Spearman's Correlation coefficient >0.5 and P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of DEGs

Following the standardization of the microarray results, 3,816 genes in the GSE63514 dataset, 2,625 genes in the GSE9750 dataset and 2,093 genes in the GSE7410 dataset were identified as DEGs. The 188 overlapping DEGs amongst the three datasets were determined using a Venn diagram (Fig. 1), and consisted of 56 downregulated genes and 132 upregulated genes between normal and cancerous tissues.

Figure 1 Differentially expressed genes were selected for using a fold-change of >2 and an adjusted P-value of <0.01 among the mRNA expression profiling datasets GSE63514, GSE9750 and GSE7410. An overlap of 188 genes was observed among the three datasets.

GO and KEGG enrichment analysis of DEGs

To determine the biological significance of the 188 DEGs, functional and pathway enrichment analyses were performed using DAVID. As shown by GO analysis, the DEGs were significantly enriched in cell cycle, cell cycle phase, mitotic cell cycle, cell cycle process and M phase (Fig. 2A). In terms of cellular components, the DEGs were primarily enriched in microtubule cytoskeleton, spindle, chromosome, nuclear lumen and chromosomal part (Fig. 2B). The molecular function terms were primarily enriched in ATP binding, adenyl ribonucleotide binding, adenyl nucleotide binding, purine nucleoside binding and nucleoside binding (Fig. 2C). These significantly enriched terms may assist in improving our understanding of the roles of the DEGs in the development and progression of CC. In the KEGG analysis, the DEGs were primarily enriched in cell cycle, DNA replication and mismatch repair (Fig. 2D).

Figure 2 Top 10 significant changes in the Gene Ontology and KEGG significant enrichment analyses of differentially expressed mRNAs. (A) Biological process. (B) Cellular components. (C) Molecular function. (D) KEGG. KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network construction and interaction network analysis of the hub genes

The screened DEGs were uploaded to the STRING website to reconstruct the PPI network. The resulting network contained 188 nodes and 2,854 edges (Fig. 3A). Subsequently, hub genes were filtered using CytoHubba. Cell division cycle 6 (CDC6), cyclin-dependent kinase 1 (CDK1), cell division control protein 45 (CDC45), budding uninhibited by benzimidazoles 1 (BUB1), TOP2A, minichromosome maintenance complex component 4 (MCM4), CCNB2 and CCNB1 were identified as the 8 hub genes, all of which were imported into STRING to reconstruct the PPI network of 8 nodes and 28 edges (Fig. 3B). A network of the hub genes and their co-expression genes was then analyzed using the cBioPortal online platform, which shows the gene-gene interaction network among the hub genes and the most frequently altered neighboring co-expression genes (Fig. 3C). These results suggest that the above hub genes may serve a critical role in the development of CC.

Figure 3 PPI network and the hub gene module of DEGs. (A) PPI network of DEGs was constructed with Cytoscape. Upregulated genes are marked in light red and downregulated genes in light blue. (B) Hub gene module obtained from the PPI network based on STRING with 8 nodes and 28 edges. (C) Hub genes and their co-expression genes were analyzed using cBioPortal. Nodes with bold black outlines represent the hub genes. Nodes with thin black outline represent the co-expression genes. PPI, protein-protein interaction; DEGs, differentially expressed genes; STRING, Search Tool for the Retrieval of Interacting Genes; CDC6, cell division cycle 6; CDK1, cyclin-dependent kinase 1; CDC45, cell division control protein 45; BUB1, budding uninhibited by benzimidazoles 1; MCM4, minichromosome maintenance complex component 4; TOP2A, DNA topoisomerase II α.

TOP2A mRNA expression and gene expression correlation analysis

cBioPortal was used to perform the survival analysis of these hub genes to evaluate their effects on CC. Briefly, only TOP2A was clearly associated with the prognosis of patients (Fig. 4A). Patients whose tissues exhibited upregulated expression levels of TOP2A had significantly shorter disease/progression-free survival times compared to those with lower expression levels (P<0.01). The expression profile of TOP2A in human tissues was displayed using UALCAN. It was found that the TOP2A mRNA expression in bladder, breast, cervical, esophageal, liver, kidney and lung cancer, as well as glioblastoma multiforme, was higher compared with the respective normal tissues (Fig. 4B). Oncomine analysis of cancer vs. normal tissue confirmed that TOP2A expression was significantly upregulated in CC (Fig. 4C). In the Zhai Cervix dataset, higher TOP2A mRNA expression levels were associated with cancer type and human papillomavirus (HPV) infection status (Fig. 5A-D).

Figure 4 Expression analysis of the TOP2A gene. (A) Disease-free survival analyses of hub genes were performed in the cBioPortal online platform. (B) Expression profiles for TOP2A in human cancer analyzed using UALCAN. (C) Oncomine analysis of TOP2A in cancer vs. normal tissue. TOP2A, DNA topoisomerase II α; TCGA, The Cancer Genome Atlas.

Figure 5 Association between the expression of TOP2A, cancer type and HPV Infection Status in the 201291_s_at and 201292_at probes in the Zhai Cervix dataset. (A and B) TOP2A mRNA expression in CC between normal CC and precursor tissues. 0, no value (3); 1, cervix squamous epithelium (10); 2, cervical squamous cell carcinoma (21); 3, high grade cervical squamous intraepithelial neoplasia (7). (C and D) TOP2A mRNA expression and HPV Infection Status. 0, no value (3); 1, HPV negative (10); 2, HPV type 16 positive (14); 3, HPV type 18 positive (weak) (1); 4, HPV type 18 positive (4); 5, HPV type 18 and 45 positive (1); 6, HPV type 33 positive (1); 7, HPV type 33, 52 and 58 positive (4); 8, HPV type 56 positive (1); 9, HPV type 58 positive (1); 10, HPV type 59 positive (1). TOP2A, DNA topoisomerase II α; CC, cervical cancer; HPV, human papillomavirus. The numbers in the brackets refers to the number of samples.

TOP2A gene expression correlation analysis results from the cBioPortal online platform showed that, of the hub genes, BUB1 (P<0.0001, rs=0.635) and CDK1 (P<0.0001, rs=0.511) had a close association with TOP2A (Fig. 6A and B). In addition, strong positive correlations were observed between TOP2A mRNA expression levels and centromere protein F (CENPF; P<0.0001, rs=0.677), Rac GTPase activating protein 1 (RACGAP1; P<0.0001, rs=0.612), F-box protein 5 (FBXO5; P<0.0001, rs=0.585) and BUB1 mitotic checkpoint serine/threonine kinase B (P<0.0001, rs=0.584; Fig. 6C-F).

Figure 6 Gene expression correlation analysis for BUB1, CDK1, CENPF, RACGAP1, FBXO5 and BUB1B with the TOP2A gene using cBioPortal. The scatter plot shows the Spearman's correlation of TOP2A expression with the expression of (A) BUB1, (B) CDK1, (C) CENPF, (D) RACGAP1, (E) FBXO5 and (F) BUB1B. BUB1, budding uninhibited by benzimidazoles 1; CDK1, cyclin-dependent kinase 1; CENPF, centromere protein F; RACGAP1, Rac GTPase activating protein 1; FBXO5, F-box protein 5; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase B; TOP2A, DNA topoisomerase II α.

Discussion

In the present study the key words ‘cervical cancer’ were used to search the GEO database and screen out the datasets which contained Homo sapiens mRNA expression profiles of both CC tissues and normal tissues. A total of three datasets (GSE63514, GSE9750 and GSE7410) were selected for further investigation. Among these, GSE63514 and GSE9750 have been analyzed together (26,27). Other studies have also mined the GSE63514(28) and GSE9750(29) datasets; only the GSE7410 dataset had not yet been thoroughly examined previously, to the best of our knowledge. Therefore, the GEO2R online tool was used to perform the analyses of the three datasets to determine any potentially relevant factors.

In the present study, 188 DEGs were identified in CC, including 132 upregulated and 56 downregulated genes. As shown by the results of KEGG and GO enrichment analysis, the 188 DEGs were primarily enriched in cell cycle. Previous studies have reported that the abnormal regulation of cell cycle process serves a vital role in the tumorigenesis and progression of several types of cancer (30-32). In addition, with regard to GO terms, the DEGs were primarily enriched in cell cycle, mitotic cell cycle, microtubule cytoskeleton, spindle, ATP binding and adenyl ribonucleotide binding, whereas changes in the KEGG pathway were largely enriched in cell cycle, DNA replication and mismatch repair. Previous studies have shown the antitumor effect of disrupting the microtubule cytoskeleton (33,34). Furthermore, HPV infection facilitating DNA replication and stimulating DNA damage response are potential important causes of a range of diseases, including CC (35). Therefore, these potential mechanisms support the results of the present study. The results of the present study may improve our understanding of the underlying mechanisms by which the identified DEGs may promote development or progression of CC. The other enriched functions and pathways may also be involved in CC carcinogenesis, and thus should be further studied.

Cytoscape is an important open source bioinformatics software platform for visualizing molecular interaction networks. The CytoHubba plugin was used to screen out core genes for subsequent prognosis analysis similar to previous studies (36,37). These hub genes may serve an important role in the development of CC. Through the PPI network, CDC6, CDK1, CDC45, BUB1, TOP2A, MCM4, CCNB2 and CCNB1 were further selected as hub genes and these hub gens were considered important. Subsequently, the hub genes were chosen for further survival analysis. Only TOP2A was associated with patient prognosis, and TOP2A dysregulation was significantly associated with reduced disease-free survival, but not overall survival (P<0.01). TOP2A is an essential nuclear enzyme for chromosome condensation, chromatid separation and the relief of supercoiled DNA during mitosis, and it is crucial for the segregation of daughter chromosomes at the end of cell division (38). Upregulated expression of TOP2A is significantly associated with increased CC cell division (39) and reduced survival periods (40-43). In addition, TOP2A may be used as an immunohistochemical biomarker for CC (39). Several studies have confirmed the use of TOP2A as a sensitive biomarker for the screening of carcinoma tissues (39,44,45). TOP2A was also frequently overexpressed in several types of cancer, including lung (46), prostate (47), breast (42,48), hepatic (49) and ovarian cancer (50), and glioma (51). Upregulated expression of TOP2A is significantly associated with the progression from cervical intraepithelial neoplasia grade 2 to more advanced cervical lesion division (38), and the Oncomine analysis performed in the present study confirmed these results. In addition, due to HPV infection being the primary etiological factor of CC, the data in the Zhai Cervix dataset suggested that high-risk HPV infection resulted in increased TOP2A expression, compared with no HPV infection. Santin et al (52) showed that TOP2A was coordinately dysregulated in primary HPV16 and HPV18-infected stage-IB-IIA CC, potentially representing a common signaling pathway initiated by HPV transformation. Thus, the encoding enzyme gene is regarded as the target of several anticancer agents (53-55) and a variety of mutations in this gene have been associated with the development of drug resistance (56-58).

To examine the underlying molecular mechanisms in CC, gene expression correlation analysis was performed in cBioPortal. CENPF, RACGAP1 FBXO5 and BUB1B had a strong positive correlation with TOP2A expression, as well as the hub genes BUB1 and CDK1 in CC. Nevertheless, previous studies have not shown interactions between TOP2A and co-expression genes BUB1, CDK1, CENPF, RACGAP1, FBXO5 and BUB1B in CC. However, Guo et al (59) suggested that BUB1, CDK1, RACGAP1 and TOP2A may be involved in the tumorigenesis of adrenocortical carcinoma. Li et al (60) showed that lower expression levels of TOP2A and CENPF were associated with improved overall survival in patients with bladder cancer. BUB1 is a mitotic checkpoint serine/threonine kinase that is crucial for physiological chromosomal segregation and is correlated with cancer stem cell potential; it may also be a target for anti-breast cancer stem cell therapies (61). CDK1 is a member of the Ser/Thr protein kinase family, which is associated with proliferation and has been demonstrated to eb a prognostic factor in the development of ovarian cancer (62,63). Schwermer et al (64) confirmed a pivotal role of the CDK1/CCNB1 complex in tumor cell survival. CENPF is a centromere-kinetochore complex and chromosomal segregation-associated protein during mitosis (65). Preclinical analysis of mouse models showed that FOXM1 and CENPF were master drug treatment-responsive genes in prostate cancer (65). RACGAP1 serves a key role in metastasis, regulating cell morphology, motility and establishment of cell polarity (66). Furthermore, RACGAP1 has been shown to function as a potential biomarker, due to its significant association with poor disease-free and overall survival in colorectal cancer (67). FBXO5, also known as EMI1, is an inhibitor of the anaphase promoting complex/cyclosome (APC/C) and is upregulated in HPV16 E7 oncoprotein-expressing mitotic cells, where it interferes with the degradation of APC/C substrates in CC (68). BUB1B encoded a kinase involved in spindle checkpoint function, and also serves a role in inhibiting the APC/C; the spindle checkpoint function has been shown to be impaired in several types of cancer (69-71). Further studies are required to confirm whether TOP2A also participates in any biological processes associated with cancer development or progression, due to its strong positive correlation with FBXO5 and BUB1B.

The present study has some limitations. Firstly, although upregulated TOP2A mRNA expression levels was a biomarker for the prognosis of CC, all the data in the present study are based on online databases; further studies with larger sample sizes are required to validate these results and to determine potential targets for diagnosis and treatment of CC. Gene Set Enrichment Analysis was not performed, which is a more powerful strategy for functional and pathway enrichment analysis, to assess the significant enrichment analysis of DEGs; thus, future studies are required to identify the significant molecular signaling pathways in CC. Finally, due to the lack of sufficient clinical sample information, multivariate Cox regression analysis could not be performed to further elucidate the significance of TOP2A. In future studies, the underlying mechanisms of TOP2A in CC should be clarified.

In conclusion, TOP2A may be used a biomarker which predicts poor prognosis and may serve as a therapeutic target for treatment of CC. Further studies are required to explore and demonstrate the potential use of the identified hub genes for the diagnosis, prognosis and treatment of CC.

Acknowledgements

The authors are grateful for the help provided by Dr Rui Hou (Harry Perkins Institute of Medical Research, The Western University of Australia).

Funding

This study was partially supported by the National Natural Science Foundation of China (grant no. 81702580) and the Natural Science Foundation of Jiangxi Province (grant nos. 2017ACB21066 and 20171BBG70051).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

QZ, HL and TZ designed the study. QZ, LZ, SH, XX, and TZ analyzed the data and prepared the figures. QZ and YL revised the manuscript. HL, YL and JL conceived the study and interpreted the data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References