This article is mentioned in:

Introduction

Cervical Cancer (CC) is a highly aggressive tumor and is one of the leading causes of cancer-associated mortality in women, with an estimated 570,000 new cases and 311,000 deaths in 2018 worldwide (1). Women with CC are considered to have a lower quality of life (2). The progression of CC, from normal cervical mucosal epithelium to cervical intraepithelial neoplasia (CIN) grade 1, 2, and 3, to CC (3) is associated with persistent high-risk human papillomavirus (HPV) infection (4). Furthermore, a number of risk factors, including early sexual activity (5), multiple sexual partners (6), long-term use of oral contraceptives (7), genetic factors [active oncogenes, including PIK3CA (8), ATAD2 (9) and CRNDE (10); tumor suppressor genes, including p53 (11), Ras association domain family 1 isoform A (12) and NOL7 (13)], tobacco use [current smoker, started smoking age ≤15 years, smoking duration ≥30 years, ≥20 cigarettes/day (14)] and other viral infections (such as HIV, herpes simplex virus (HSV) type II and bacterial infections caused by Chlamydia trachomatis) (15) have been associated with CC progression.

HPV infection plays a leading role in CC (16). The DNA of HPV integrates into the host cell genome [HPV16: q21-q31 of chromosome no. 13; HPV18: q24 of chromosome no. 8 (17)], disrupts the open reading frame and causes overexpression of E6 and E7 genes (18). It has been verified that E6 and E7 exert carcinogenic effects by binding to the cell cycle regulators, p53 and retinoblastoma (Rb) (19). While E6 and E7 proteins are upregulated, E6 can interact with its associated protein [E6-associated protein, E6AP (20)] to form a complex and bind to p53. This binding hydrolyzes p53 and results in the loss of p53-induced negative regulation of cell proliferation, thereby leading to unchecked cellular proliferation and malignant transformation (21). E7 has a high affinity for Rb, which controls the cell cycle. Binding of E7 to Rb can dissociate the Rb-E2F complex, thus releasing E2F to exert its role as a transcription factor, which leads to an uncontrolled cell cycle and cellular immortalization (22,23). Furthermore, centrosomes are central regulators of mitosis that are often increased in numbers in cancer cells (24). A previous study indicated that an abnormally increased number of centrosomes is associated with structural chromosomal abnormality in cervical lesions with high risk of HPV infection (25). Duensing and Münger (26) reported that abnormal number of centrosomes and associated spindle mitotic abnormality can be found in cells infected by the high-risk HPV16 E6 and E7 proteins, but not in cells infected by the low-risk HPV6. Although numerous experimental studies on genes [HPV16 L1 protein (27), sonic hedgehog (28) and FGFR4 (29)], and signaling pathways [Wnt/β-catenin signaling pathway (30), adenosinergic pathway (31) and ERK signal transduction (32)], as well as bioinformatics analyses have focused on microRNAs (33) and genes (34) associated with CC, and have provided an understanding of the pathophysiological mechanisms of the disease over the last decade, the underlying molecular mechanisms remain unclear.

The development of CC occurs over a number of years and its complexity presents clinical challenges in patients screening and treatment. Currently, The Bethesda System (35), which is a tool that is used to report Pap smear results for cervical cytologic diagnoses, provides useful data that allows research into the epidemiology, biology and pathology of cervical lesions; however, its diagnostic value remains poor (36). Instead, direct biopsy remains the gold standard for diagnosis. Nevertheless, invasive examinations may cause adverse psychological effects, including anxiety, depression or distress (37). Surgery, chemotherapy and radiotherapy (38) are the three major therapeutic strategies in the treatment of CC; however, their uses may be limited for various reasons. Surgery may be limited by the status and stage of patients, including late stage or tolerance to anesthesia (39), whereas chemotherapy is limited due to the lack of sensitivity and the development of drug resistance (40). In addition, radiotherapy can be limited by the maximum tolerated dose to adjacent normal tissues (41). Thus, it is essential to understand the underlying molecular mechanisms in the initiation and development of CC, in order to develop methods for its accurate diagnosis and effective treatment. A number of studies have reported that multiple genes [CXCL12 (42), FGFR4 (29) and SHH (43)], proteins [cyclin D1 (44), FOXO1 (45) and BASP1 (46)] and pathways [Toll-like signaling pathway (47), VEGF signaling pathway (48) and Wnt signaling pathway (30)] are involved in the natural progression of CC; however, few studies have investigated the fundamental pathological molecular mechanisms in the progression of CC (from normal, to CIN1, CIN2, CIN3, to cancer). Thus, the specific pathological processes remain unclear.

The present study provided a systematic investigation of the development of CC and further understanding of the associations between the four phases of CC progression, and thus revealed additional targets for the detection and treatment of CC. A flow diagram of the present study is presented in Fig. 1.

Figure 1. Flow diagram of the present study. FC, fold change; DEGs, differentially expressed genes; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins; GO, Gene Ontology; FDR, false discovery rate; KEGG, Kyto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction.

Materials and methods

Identification of differentially expressed genes (DEGs)

The CC gene expression profile in the GSE63514 dataset, acquired using the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) provided by den Boon in 2015 (49), was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The profile contained 128 cervical specimens, including: Normal (n=24), CIN1 (n=14), CIN2 (n=22), CIN3 (n=40) and cancer (n=28) samples. All samples were divided into four phases as follows: Phase I, normal to CIN1; phase II, CIN1 to CIN2; phase III, CIN2 to CIN3 and phase IV, CIN3 to cancer, and GEO2R tools (https://www.ncbi.nlm.nih.gov/geo/geo2r/) (50) within the limma package version 3.26.8 (51) were used to screen the DEGs at the four phases. The criteria fold change (FC) of expression >2 and P<0.05 were used to identify DEGs.

Identification of hub genes

The Search Tool for the Retrieval of Interacting Genes (STRING) database (52) and Cytoscape software (version 3.5.1) (53) were used to identify the hub genes in the four phases. The PPI network was constructed by searching for gene symbols and the minimum required interaction score was set at 0.7, to ensure high confidence in the results. The nodes that not connect to the major network were removed to decrease the error detection rate. CytoHubba (54), a plug-in for Cytoscape software, was used to investigate notable nodes in the interactome network using 12 topological algorithms, including Degree, Edge Percolated Component, Maximum Neighborhood Component, Density of Maximum Neighborhood Component and Maximal Clique Centrality, and centralities based on shortest paths, such as Bottleneck, EcCentricity, Closeness, Radiality, Betweenness, Clustering Coefficient and Stress. The genes that ranked in the top 10 for each topological algorithm were extracted and the duplication of each gene was calculated. Genes duplicated <2 times were excluded, in order to guarantee that the genes were associated with CC. The remaining genes were considered as hub genes in the four phases.

Functional enrichment analyses of GO and pathways

The functional features of the genes associated with the four phases were examined using WebGestalt (55) and ToppGene (56). In WebGestalt, over-representation analysis was selected as the enrichment method, Biological Process in GO as the functional database, gene symbol as the gene ID type and genome as the reference set for enrichment analysis. In ToppGene, two frequently used databases, Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp/) and BioCarta (https://www.biocarta.com/), were utilized to perform pathway enrichment analysis, to improve the reliability of the results. Pathways with a false discovery rate of P<0.05 were considered to indicate significantly enriched pathways.

Pathway crosstalk analysis

Pathway crosstalk analysis was performed (57), to investigate the interactions among the significantly enriched pathways. The pathways with either a false discovery rate of P>0.05 or <3 genes were removed as selection criteria. The number of shared genes between pairwise pathways was calculated and pairwise pathways with <2 overlapping genes were removed. The Jaccard Coefficient (JC) and the Overlap Coefficient (OC) parameters were calculated, to measure the overlap between the pathways. Specifically, JC=|A∩BA∪B|, while OC=|A∩Bmin(|A|,|B|)|, where A and B represent the gene numbers of the tested pathways. The interrelationships between pathways were visualized using Cytoscape software.

To determine the overall progression and further detect an association between two pathways, the KEGG and BioCarta databases were used to identify the upstream or downstream associations between pathways. Furthermore, the nodal degree was calculated using Centiscape (58) to identify key nodes. According to Han et al (59), key nodes are considered as those with a nodal degree ≥5.

Integration of the PPI network

The PPI network (60) was used to identify key proteins for the four phases of CC. As Protein Interaction Network Analysis (PINA) (https://omics.bjcancer.org/pina/) (61) is an integrated platform for protein interaction network construction, analysis and visualization, it can identify the associations between the queried genes based on integration of data from six public PPI databases: IntAct (62), MINT (63), BioGRID (64), DIP (65), HPRD (66) and MIPS MPact (67). Thus, the PINA4MS plug-in for Cytoscape software was used to construct the PPI network, to identify CC progression-associated genes. As PINA4MS requires UniProt accession numbers, the UniProt Retrieve/ID mapping tool (https://www.uniprot.org/uploadlists/) was used to input gene symbols. The key nodes for the PPI network were also extracted using a criterion of nodal degree ≥5.

Comprehensive gene-pathway analysis

To determine the molecular mechanisms and associations between the key genes and pathways, the gene-pathway network was constructed by examining the key pathways, in order to determine which pathway contained at least one of the key genes.

Co-expression and survival analysis for key genes

To identify the co-expression of key genes and their impact on OS time, the LinkedOmics database (68) was used, which was based on TCGA (69). The co-expression analysis was performed using Pearson correlation and OS analysis was assessed with Cox regression method. For survival analysis, samples were divided by the median value of the investigated gene. P<0.05 was considered to indicate a statistically significant difference for both the co-expression correlation and OS time.

Results

Identification of DEGs

Analysis of the GSE63514 dataset using GEO2R, with a criteria of ≥2 FC and P<0.05, identified a total of 3,446 DEGs for the four phases as follows: 446 DEGs in phase I, of which 76 were upregulated and 370 were downregulaged; 382 DEGs in phase II, of which 146 were upregulated and 236 were downregulated; 756 DEGs in phase III, of which 435 were upregulated and 321 were downregulated; 1,862 DEGs in phase IV, of which 816 were upregulated and 1046 were downregulated.

Identification of hub genes

Following removal of 2,256 irrelevant genes (Phase I, 265; Phase II, 197; Phase III, 603; Phase IV, 1191), 12 topological algorithms were used and the top 10 genes for each method were extracted. A total of 107 genes that appeared at least twice were conserved as hub genes, as presented in Table I. A total of 29 genes were identified in phase I, among which five genes were members of the kinesin family (KIF11, KIF15, KIF23, KIF4A and KIF14), and five genes were associated with meiosis and the maturation of oocytes [BUB1B (70), BUB1 (71), CCNA2 (72), CCNB1 (72) and CDK1 (73)], as well as other genes associated with inflammation and innate immune responses [STAT1 (74), GBP1 (75) and RHOA (76)]. A total of 25 hub genes were verified in phase II, among which the involvement of seven interferon-induced genes was identified (IFI44L, IFIT3, IFIF1, IFIF5, IFI44, IFIT2 and IFI6), and several pattern recognition receptor-associated genes [IRF7 (77), STAT1 (78) and CXCL10 (79)], as well as some genes involved in invasion and metastasis of cancer cells [HGF (80), IGF1 (81), KIT (82), FN1 (83) and CXCL12 (84)]. A number of common cancer-associated signaling pathway genes were identified in phase III [CCND1 (85), STAT1 (86) and VEGFA (87)]. A total of three C-X-C motif chemokine ligands (CXCL8, CXCL11 and CXCL4), two integrin subunits (ITGB1 and ITGA1), and one mitogen-activated protein (MAPK12) were identified in phase IV. Furthermore, PIK3CA (88) and FOS (89) participated in cancer-associated pathways in phase IV. The diversity of genes within the four phases demonstrated that CC progression is a complex process and its molecular mechanisms are not constant.

Table I. Hub genes in phases I–IV.

Table I.

Hub genes in phases I–IV.

Phase	Gene	Regulation	Counts	LogFC	P-value
Phase I	CDK1	+	7	1.24	<0.01
	KIF11	+	7	1.27	0.01
	BUB1B	+	6	1.13	0.02
	BUB1	+	5	1.17	0.01
	CCNA2	+	5	1.01	<0.01
	HLA-DPA1	+	5	1.49	0.03
	CENPE	+	5	1.36	0.02
	RHOA	+	4	1.47	<0.01
	KIF15	+	4	1.11	<0.01
	CCNB1	+	4	1.23	0.01
	NDC80	+	4	1.47	0.01
	TTK	+	4	1.14	<0.01
	STAT1	–	4	−1.47	0.03
	CXCL10	+	4	2.47	<0.01
	KIF23	+	4	1.48	0.01
	KIF4A	+	3	1.00	0.05
	PSMB9	+	3	1.12	0.03
	GNG2	+	3	1.26	0.04
	SPAG5	+	2	1.15	0.01
	TRIP13	+	2	1.09	0.02
	ANLN	+	2	1.19	0.02
	CDKN3	+	2	1.83	<0.01
	KIF14	+	2	1.15	0.01
	MKI67	+	2	1.45	0.01
	NUSAP1	+	2	1.22	0.02
	NEK2	+	2	1.30	0.01
	NCAPG	+	2	1.18	<0.01
	DLGAP5	+	2	1.54	<0.01
	GBP1	+	2	1.45	<0.01
Phase II	STAT1	–	9	−1.03	0.02
	CXCL10	–	6	−2.16	0.02
	CXCL12	–	6	−1.48	0.04
	DCN	–	6	−1.37	0.01
	CCL2	–	5	−1.54	0.01
	KIT	–	5	−1.26	0.05
	IGF1	–	5	−1.86	<0.01
	OAS2	–	4	−1.37	0.01
	IRF7	–	4	−1.08	0.01
	ISG15	–	4	−1.83	0.02
	FN1	–	4	−1.34	0.04
	HGF	–	4	−1.22	0.02
	HERC6	–	3	−1.71	<0.01
	MX2	–	3	−1.86	0.01
	IFIT3	–	3	−1.79	<0.01
	IFIT1	–	3	−2.96	<0.01
	GBP1	–	3	−1.20	0.01
	CDC6	+	3	1.40	0.01
	IFIT5	–	2	−1.11	<0.01
	IFI6	–	2	−1.55	0.01
	SP110	–	2	−1.43	<0.01
	IFI44	–	2	−1.87	<0.01
	DDX60	–	2	−1.18	<0.01
	IFIT2	–	2	−1.67	<0.01
	RSAD2	–	2	−2.66	<0.01
Phase III	BIRC5	+	9	1.13	<0.01
	TOP2A	+	8	1.36	<0.01
	KIF2C	+	6	1.04	<0.01
	MCM10	+	6	1.16	0.01
	VEGFA	+	6	1.27	<0.01
	MAD2L1	+	5	1.03	<0.01
	KIF15	+	5	1.45	<0.01
	ASPM	+	5	1.70	<0.01
	FOXM1	+	5	1.18	0.01
	MX2	+	4	1.34	0.01
	STAT1	+	4	1.20	<0.01
	PLXNA4	–	4	−1.15	0.01
	AR	–	4	−1.69	<0.01
	CCND1	–	3	1.20	<0.01
	OAS2	+	3	1.05	<0.01
	ACLY	+	3	1.55	0.01
	GNG2	+	3	−1.21	<0.01
	RSAD2	+	2	1.85	<0.01
	ISG15	+	2	1.82	<0.01
	IFI35	+	2	1.35	0.01
	IRF5	+	2	1.46	<0.01
	SAMHD1	+	2	1.05	<0.01
	MKI67	+	2	1.20	<0.01
	PLK4	+	2	1.03	0.01
	AHCTF1	+	2	1.03	<0.01
	NUDC	+	2	1.07	<0.01
	EXO1	+	2	1.29	<0.01
	PLXNA3	+	2	1.01	0.01
	MMP9	+	2	2.55	<0.01
	PLAUR	+	2	1.01	<0.01
Phase IV	PIK3CA	+	9	1.42	<0.01
	CXCL8	+	8	1.28	0.05
	ITGB1	+	8	2.20	<0.01
	PTK2	+	8	1.33	<0.01
	GNG2	+	6	1.18	0.05
	ITGA1	+	6	1.33	<0.01
	GNG12	–	6	−1.14	<0.01
	FOS	–	5	−1.13	0.05
	EDN1	+	5	1.07	<0.01
	NMU	–	4	−2.23	<0.01
	LPAR5	–	4	−1.36	<0.01
	STAT1	+	3	1.75	<0.01
	FN1	+	3	3.61	<0.01
	GSTM1	–	3	−1.09	0.02
	PLA2G4A	–	3	−1.41	0.02
	CXCR4	+	2	1.58	<0.01
	HCAR3	–	2	−1.60	<0.01
	S1PR5	–	2	−1.56	<0.01
	CXCL5	–	2	−2.10	0.01
	NQO1	–	2	−1.38	0.01
	CXCL11	+	2	1.69	0.02
	COMP	+	2	1.62	0.01
	MAPK12	+	2	1.38	<0.01

[i] +, upregulated; -, downregulated.

GO enrichment analysis of hub genes

To further identify the biological functions and locations of hub genes, GO enrichment analysis (90) was performed (Fig. 2). Hub genes were notably enriched in ‘biological regulation’, ‘metabolic process’ and ‘cellular component organization’ in phase I and II, while ‘responses to stimulus’ and ‘biological regulation’ were predominantly enriched at phases III–IV in biological process. For the cellular components, ‘nucleus’, ‘membrane-enclosed lumen’ and ‘macromolecular complex’ was enriched at phases I–III, while ‘chromosome’ and ‘membrane’ was identified in phases II and IV, respectively. ‘Protein binding’ was enriched at all four phases for Molecular Function. Furthermore, ‘nucleic acid binding’ and ‘hydrolase activity’ were enriched at phase II and III, while ‘ion binding’ was enriched at phase III and IV.

Figure 2. Gene ontology enrichment analysis of hub genes for phase I, II, III and IV. The number in each phase represents the gene count.

Pathway enrichment analysis of hub genes

As presented in Table II, a total of 10 notably enriched pathways were identified at phase I, of which five pathways were associated with virus infections including, ‘influenza A’, ‘tuberculosis’, ‘herpes simplex infection’, ‘viral carcinogenesis’ and ‘Epstein-Barr virus infection’, and additional pathways involved in the ‘cell cycle’, ‘oocyte meiosis’ and ‘progesterone-mediated oocyte maturation’. Furthermore, the chemokine signaling pathway was also identified in phase I. The RIG-I-like receptor and Toll-like receptor signaling pathways were identified in phase II, and are associated with pattern-recognition receptors (91). In addition, several pathways, including ‘focal adhesion’, ‘Rap1 signaling pathway’, ‘Ras signaling pathway’, ‘PI3K-Akt signaling pathway’ and ‘Proteoglycans in cancer’ were associated with invasion and metastasis (92–96). The two common cancer-associated signaling pathways ‘Pathways in cancer’ and ‘Proteoglycans in cancer’, were enriched in phase III, while 70 pathways were significant enriched at phase IV (P<0.05). Apart from the common cancer-associated signaling pathways and virus infection pathways at phase IV, the ‘IL-17 signaling pathway’, ‘VEGF signaling pathway’ and ‘endocrine resistance’ also were also demonstrated to be associated with CC progression. Furthermore, the ‘AGE-RAGE signaling pathway in diabetic complications’ was also identified at phases II–IV.

Table II. Pathway enrichment analysis for phases I–IV.

Table II.

Pathway enrichment analysis for phases I–IV.

Phase	Pathway	FDR	Involved genes
Phase I	Cell cycle	<0.001	CDK1, TTK, CCNA2, BUB1, CCNB1, BUB1B
	Progesterone-mediated oocyte maturation	0.002	CDK1, CCNA2, BUB1, CCNB1
	Chemokine signaling pathway	0.006	RHOA, CXCL10, STAT1, GNG2
	Oocyte meiosis	0.018	CDK1, BUB1, CCNB1
	NOD-like receptor signaling pathway	0.036	RHOA, GBP1, STAT1
	Influenza A	0.036	HLA-DPA1, CXCL10, STAT1
	Tuberculosis	0.037	RHOA, HLA-DPA1, STAT1
	Herpes simplex infection	0.038	CDK1, HLA-DPA1, STAT1
	Viral carcinogenesis	0.044	RHOA, CDK1, CCNA2
	Epstein-Barr virus infection	0.044	CDK1, HLA-DPA1, CCNA2
Phase II	Influenza A	<0.001	CCL2, OAS2, IRF7, RSAD2, CXCL10, STAT1
	NOD-like receptor signaling pathway	<0.001	GBP1, CCL2, OAS2, IRF7, STAT1
	Herpes simplex infection	<0.001	CCL2, OAS2, IRF7, IFIT1, STAT1
	Pathways in cancer	0.001	HGF, IGF1, FN1, KIT, CXCL12, STAT1
	Hepatitis C	0.001	OAS2, IRF7, IFIT1, STAT1
	Cytokine-cytokine receptor interaction	0.001	HGF, CCL2, KIT, CXCL10, CXCL12
	RIG-I-like receptor signaling pathway	0.003	IRF7, ISG15, CXCL10
	Chemokine signaling pathway	0.003	CCL2, CXCL10, CXCL12, STAT1
	Genes encoding secreted soluble factors	0.003	HGF, CCL2, IGF1, CXCL10, CXCL12
	Proteoglycans in cancer	0.004	HGF, IGF1, FN1, DCN
	AGE-RAGE signaling pathway in diabetic complications	0.005	CCL2, FN1, STAT1
	Toll-like receptor signaling pathway	0.006	IRF7, CXCL10, STAT1
	Ensemble of genes encoding extracellular matrix and extracellular matrix-associated proteins	0.008	HGF, CCL2, IGF1, FN1, DCN, CXCL10, CXCL12
	Measles	0.009	OAS2, IRF7, STAT1
	PI3K-Akt signaling pathway	0.015	HGF, IGF1, FN1, KIT
	Focal adhesion	0.024	HGF, IGF1, FN1
	Rap1 signaling pathway	0.026	HGF, IGF1, KIT
	Ras signaling pathway	0.029	HGF, IGF1, KIT
	Ensemble of genes encoding ECM-associated proteins including ECM-affilaited proteins, ECM regulators and secreted factors	0.032	HGF, CCL2, IGF1, CXCL10, CXCL12
Phase III	Pathways in cancer	0.002	BIRC5, CCND1, MMP9, AR, STAT1, GNG2, VEGFA
	Bladder cancer	0.006	CCND1, MMP9, VEGFA
	Hepatitis B	0.011	BIRC5, CCND1, MMP9, STAT1
	Pancreatic cancer	0.012	CCND1, STAT1, VEGFA
	Proteoglycans in cancer	0.025	PLAUR, CCND1, MMP9, VEGFA
	AGE-RAGE signaling pathway in diabetic complications	0.027	CCND1, STAT1, VEGFA
Phase IV	Chemokine signaling pathway	<0.001	GNG12, CXCL11, CXCL5, PIK3CA, CXCR4, PTK2, STAT1, CXCL8, GNG2
	Pathways in cancer	<0.001	FN1, LPAR5, GNG12, ITGB1, PIK3CA, CXCR4, FOS, PTK2, STAT1, CXCL8, GNG2
	Fluid shear stress and atherosclerosis	<0.001	GSTM1, NQO1, MAPK12, PIK3CA, FOS, EDN1, PTK2
	Signaling of Hepatocyte Growth Factor Receptor	<0.001	ITGA1, ITGB1, PIK3CA, FOS, PTK2
	PI3K-Akt signaling pathway	<0.001	ITGA1, FN1, COMP, LPAR5, GNG12, ITGB1, PIK3CA, PTK2, GNG2
	B Cell Survival Pathway	<0.001	ITGA1, ITGB1, PIK3CA, FOS
	AGE-RAGE signaling pathway in diabetic complications	<0.001	MAPK12, FN1, PIK3CA, EDN1, STAT1, CXCL8
	Toll-like receptor signaling pathway	<0.001	MAPK12, CXCL11, PIK3CA, FOS, STAT1, CXCL8
	Aspirin Blocks Signaling Pathway Involved in Platelet Activation	<0.001	PLA2G4A, ITGA1, ITGB1, PTK2
	Erk and PI-3 Kinase Are Necessary for Collagen Binding in Corneal Epithelia	<0.001	ITGA1, ITGB1, PIK3CA, PTK2
	Pertussis	<0.001	MAPK12, CXCL5, ITGB1, FOS, CXCL8
	Focal adhesion	<0.001	ITGA1, FN1, COMP, ITGB1, PIK3CA, PTK2
	TNF signaling pathway	<0.001	MAPK12, CXCL5, PIK3CA, FOS, EDN1
	Leukocyte transendothelial migration	<0.001	MAPK12, ITGB1, PIK3CA, CXCR4, PTK2
	Regulation of actin cytoskeleton	<0.001	ITGA1, FN1, GNG12, TGB1, PIK3CA, PTK2
	VEGF signaling pathway	<0.001	PLA2G4A, MAPK12, PIK3CA, PTK2
	PTEN dependent cell cycle arrest and apoptosis	<0.001	ITGB1, PIK3CA, PTK2
	uCalpain and friends in Cell spread	<0.001	ITGA1, ITGB1, PTK2
	Trefoil Factors Initiate Mucosal Healing	<0.001	ITGB1, PIK3CA, PTK2
	Prolactin signaling pathway	<0.001	MAPK12, PIK3CA, FOS, STAT1
	Leishmaniasis	<0.001	MAPK12, ITGB1, FOS, STAT1
	Inhibition of Cellular Proliferation by Gleevec	<0.001	PIK3CA, FOS, STAT1
	CXCR4 Signaling Pathway	<0.001	PIK3CA, CXCR4, PTK2
	TPO Signaling Pathway	<0.001	PIK3CA, FOS, STAT1
	Bacterial invasion of epithelial cells	<0.001	FN1, ITGB1, PIK3CA, PTK2
	mCalpain and friends in Cell motility	<0.001	ITGA1, ITGB1, PTK2
	ECM-receptor interaction	<0.001	ITGA1, FN1, COMP, ITGB1
	Small cell lung cancer	<0.001	FN1, ITGB1, PIK3CA, PTK2
	EGF Signaling Pathway	<0.001	PIK3CA, FOS, STAT1
	IL-17 signaling pathway	<0.001	MAPK12, CXCL5, FOS, CXCL8
	PDGF Signaling Pathway	<0.001	PIK3CA, FOS, STAT1
	Amoebiasis	<0.001	FN1, PIK3CA, PTK2, CXCL8
	Endocrine resistance	<0.001	MAPK12, PIK3CA, FOS, PTK2
	Proteoglycans in cancer	<0.001	MAPK12, FN1, ITGB1, PIK3CA, PTK2
	Chagas disease (American trypanosomiasis)	<0.001	MAPK12, PIK3CA, FOS, CXCL8
	Agrin in Postsynaptic Differentiation	<0.001	ITGA1, ITGB1, PTK2
	Integrin Signaling Pathway	<0.001	ITGA1, ITGB1, PTK2
	Fc Epsilon Receptor I Signaling in Mast Cells	<0.001	PLA2G4A, PIK3CA, FOS
	Cholinergic synapse	<0.001	GNG12, PIK3CA, FOS, GNG2
	Platelet activation	0.001	PLA2G4A, MAPK12, ITGB1, PIK3CA
	Osteoclast differentiation	0.001	MAPK12, PIK3CA, FOS, STAT1
	Dopaminergic synapse	0.001	MAPK12, GNG12, FOS, GNG2
	Hepatitis C	0.001	MAPK12, PIK3CA, STAT1, CXCL8
	Cytokine-cytokine receptor interaction	0.001	CXCL11, CXCL5, CXCR4, ACKR3, CXCL8
	Hepatitis B	0.001	PIK3CA, FOS, STAT1, CXCL8
	Phospholipase D signaling pathway	0.001	PLA2G4A, LPAR5, PIK3CA, CXCL8
	Shigellosis	0.001	MAPK12, ITGB1, CXCL8
	Fc epsilon RI signaling pathway	0.002	PLA2G4A, MAPK12, PIK3CA
	Influenza A	0.002	MAPK12, PIK3CA, STAT1, CXCL8
	Axon guidance	0.002	ITGB1, PIK3CA, CXCR4, PTK2
	Integrin Signaling Pathway	0.002	ITGA1, PIK3CA, PTK2
	Salmonella infection	0.003	MAPK12, FOS, CXCL8
	MAPKinase Signaling Pathway	0.003	MAPK12, FOS, STAT1
	Rap1 signaling pathway	0.003	MAPK12, LPAR5, ITGB1, PIK3CA
	Rheumatoid arthritis	0.003	CXCL5, FOS, CXCL8
	Th1 and Th2 cell differentiation	0.003	MAPK12, FOS, STAT1
	Circadian entrainment	0.003	GNG12, FOS, GNG2
	Inflammatory mediator regulation of TRP channels	0.003	PLA2G4A, MAPK12, PIK3CA
	Ras signaling pathway	0.003	PLA2G4A, GNG12, PIK3CA, GNG2
	Choline metabolism in cancer	0.003	PLA2G4A, PIK3CA, FOS
	Retrograde endocannabinoid signaling	0.003	MAPK12,GNG12,GNG2
	T cell receptor signaling pathway	0.004	MAPK12, PIK3CA, FOS
	Th17 cell differentiation	0.004	MAPK12, FOS, STAT1
	Serotonergic synapse	0.004	PLA2G4A, GNG12, GNG2
	Toxoplasmosis	0.004	MAPK12, ITGB1, STAT1
	Glutamatergic synapse	0.004	PLA2G4A, GNG12, GNG2
	MAPK signaling pathway	0.004	PLA2G4A, MAPK12, GNG12, FOS
	Sphingolipid signaling pathway	0.005	S1PR5, MAPK12, PIK3CA
	NOD-like receptor signaling pathway	0.012	MAPK12, STAT1, CXCL8
	cAMP signaling pathway	0.018	HCAR3, PIK3CA, FOS

[i] FDR, false discovery rate.

Pathway crosstalk analysis

The pathway crosstalk analysis results are presented in Fig. S1. A total of two major modules were identified in phase I, one of which was predominantly associated with inflammatory responses to viral infections, such as the NOD-like receptor signaling pathway and the chemokine signaling pathway, while the other module was associated with cycle regulation of oocytes, including cell cycle, progesterone-mediated oocyte maturation and oocyte meiosis (Fig. S1A). Similarly, the pathways were grouped into two modules in phase II. One module consisted of immune responses (RIG-I-like receptor and Toll-like receptor signaling pathways, and NOD-like receptor signaling pathway), which may trigger rapid activation of innate immunity by inducing the production of proinflammatory cytokines (97–99). The other module was predominantly involved in the regulation of cell proliferation and invasion (Ras signaling pathway, PI3K-Akt signaling pathway, Rap1 signaling pathway, focal adhesion and other cancer-associated signaling pathways) (Fig. S1B). Notably, the two modules in phases I and II were not independent as they were demonstrated to connect with each other via several signaling pathways. All pathways formed a cluster and common cancer-associated signaling pathways were indicated to play a critical role in phase III (Fig. S1C), while the connection between the pathways became highly complex in phase IV (Fig. S1D). The possible molecular mechanisms, such as inflammation caused by virus infections, pathways associated with cell invasion, IL-17 and VEGF signaling pathways, and endocrine resistance, are associated with one another. These results further verified the complexity of CC. A comprehensive combination of pathway crosstalk analysis containing 47 nodes and 105 edges among the four phases is presented in Fig. 3A. By analyzing the nodes with degrees ≥5, a subnetwork containing 11 key pathways was extracted (Fig. 3B). It indicates that the MAPK signaling pathway (degree=33), PI3K-Akt signaling pathway (degree=21) and focal adhesion (degree=15), which are ranked as the top three nodes and most interactive, may play critical roles in the progression of CC.

Figure 3. Comprehensive pathways crosstalk analysis. (A) Combination of pathways crosstalk analysis for the four phases and the (B) subnetwork with nodal degree ≥5. Green, phase I; yellow, phase II; orange, phase III; red, phase IV. The arrow represents the up/downstream associations between the pathways. MAPK signaling pathway (degree=33), PI3K-Akt signaling pathway (degree=21) and Focal adhesion (degree=15) rank as the top three pathways, whereby the majority of other pathways transfer information with them.

PPI network analysis

A PPI network containing 51 nodes and 78 edges was constructed (Fig. 4A) by downloading the hub genes into the PINA database. Based on the description of a previous study (59), which defined the main nodes as nodes with degree >5, 14 key genes were identified from the PPI network (Fig. 4B). CDK1, FN1 and ITGB1 rank first, second and third, respectively, as the top three-degree levels (16, 12 and 8, respectively). STAT1 was the only gene demonstrated to be involved at all four phases. Furthermore, MMP9 presented self-regulating functions and was demonstrated to co-express with FN1 and ITGB1.

Figure 4. Protein-protein interaction network analysis. (A) Protein-protein interaction network downloaded from The Protein Interaction Network Analysis platform and (B) subnetwork with nodal degree ≥5. Green, phase I; blue, phase II; yellow, phase II; red, phase IV. The nodal size represents the degree of each node. CDK1 (degree=16), FN1 (degree=12) and ITGB1 (degree=8) rank as the top three proteins. MMP9 was self-regulated and co-expressed with FN1 and ITGB1.

Comprehensive gene-pathway analysis

After mapping the key genes onto the key pathways using the KEGG and BioCarta databases, a potential gene-pathway flowchart, including eight key pathways and six key genes was constructed (Fig. 5). The results demonstrated the following: For phase I, CDK1 and CCNB1 participated in the regulation of the cell cycle, while CDK1 was also involved in viral carcinogenesis; for phases II–IV, ‘pathways in cancer,’ ‘focal adhesion’ and ‘PI3K-Akt signaling pathway’ were ranked the top three pathways according to the number of genes involved; FN1, ITGB1 and MMP9 may be associated with metastasis of tumor cells, and STAT1 participated in ‘pathways in cancer’ and ‘Toll-like receptor signaling pathway’, which functioned at a phase IV.

Figure 5. Gene-pathway flowchart for the key genes and the pathways in phases I–IV. Circle, gene; rectangle, pathway. CDK1 and CCNB1 regulate the cell cycle and they are activated in phase I. For phases II–IV, ‘pathways in cancer’, ‘focal adhesion’ and ‘PI3K-Akt signaling pathway’ rank as the top three pathways according to the number of genes involved.

Co-expression network and survival analysis for key genes

By mining the data from LinkedOmics, the results of co-expression demonstrated that CDK1 had a significantly positive correlation with CCNB1 (P<0.0001), but negative correlation with FN1 (P=0.003) and MMP9 (P=0.001), respectively (Fig. 6). CCNB1 demonstrated a significantly negative correlation with ITGB1 (P=0.047); however, FN1, ITGB1 and MMP9 indicated a significantly positive correlation between each other (FN1 and ITGB1, P<0.001; FN1 and MMP9, P<0.0001; ITGB1 and MMP9, P=0.023). STAT1 was significantly positively correlated with MMP9 (P<0.0001).

Figure 6. Co-expression analyses for the key genes. Positive correlation was detected between CDK1 and CCNB1; between FN1 and ITGB1, and MMP9, and between STAT1 and MMP9. Color in each grid represents the correlation coefficient between two genes. The values in the color legend represent the correlation coefficient. *P<0.05, **P<0.01, ***P<0.001.

Survival analysis indicated that patients with higher FN1 and ITGB1 expression levels had a significantly worse OS time (FN1, P=0.00080; ITGB1, P=0.00005; Fig. 7). However, CDK1, CCNB1, MMP9 and STAT1 were not demonstrated to have a significant effect on OS.

Figure 7. Survival analysis for the key genes. (A) FN1, (B) ITGB1, (C) CDK1 and (D) CCNB1, (E) MMP9 and (F) STAT1. ITGB1 and FN1 have significant effect on overall survival.

Discussion

To date, the occurrence and development of CC is hypothesized to be linked with persistent HPV infection (100); however, the specific molecular mechanisms require further investigation. In addition, although a number of studies are examining the molecular mechanisms of CC (101–104), the detailed pathological process remains unclear.

The results of the present study indicated that CDK1 had the highest degree, and participated in the cell cycle with CCNB1 and viral carcinogenesis in phase I. CDK1 is one of the major cell cycle regulatory proteins and operates at the center of the cell cycle regulatory network (105). It regulates the G-S phase transition and initiates DNA replication (106). Furthermore, as a core molecule at the M phase checkpoint, CDK1 plays a role in the regulation of G2 phase, at M phase in the cell cycle (107). The HPV infection pathway is regulated within the viral carcinogenesis pathway (19). E6 proteins inactivate p53 by binding (108), while p53 negatively regulates CDK1 transcription under normal physiological conditions (109). In addition, the activation of cyclin B1-CDK1 is the key event that initiates the start of mitosis (110). Centrosome separation can be regulated by CDK1 (111), and cyclin B1-CDK1 remains activated following centrosome separation (110). Hence, overexpression of CDK1 can cause dysfunction in cell cycle progression, failure of normal proliferation and differentiation, and thereby lead to malignant proliferation of cancer cells and the formation of CC. CCNB1 is a notable member of the cyclin family, a key initiator and a stringent quality control step of mitosis (112). It also plays a key role in the regulation of CDK1, and its phosphorylated substrates can promote the transition of the cell cycle from G2 to mitosis (113,114). Amplification of the HPV genome depends on prolongation of the G2 phase in the cell cycle (115). CCNB1 is a downstream target of STAT3, which is a key gene that regulates the proliferation and differentiation of CC cells (116). In cells with inactivated STAT3, CCNB1 expression is downregulated and amplification of the HPV genome is also decreased, resulting in decreased activity of CC cells (116). As the results of the present study demonstrated that CDK1 and CCNB1 occurred in phase I and functioned as regulators of proliferation and differentiation, they may be potential promoters of CIN and CC.

In the process of tumor invasion and metastasis, cancer cells can bind to ligands of the extracellular matrix (ECM) via integrins and degrade the basement membrane (BM) by secreting proteases via the pathways of focal adhesion and the PI3K-Akt signaling pathway (117). This degradation is also the prerequisite for stromal infiltration and cancer cell migration (118). ITGB1 belongs to the integrin family and FN1 is the ligand. The binding of ITGB1 and FN1 induces the phosphorylation of tyrosine and directly affects cytoskeleton reconstruction and signal transduction activities of the Ras-MAPK signaling pathway via the RAP1 signaling pathway, which initiates the expression of MMP genes (119). MMPs are a family of calcium and zinc-dependent proteases that degrade a variety of components of the ECM (120). Collagen type IV is the main scaffold in the BM of the ECM and also the main substrate of MMP9 (121). MMP9 can decompose the nestin in the BM to destroy the cells integrity and promote the invasion and metastasis of cancer cells (122). MMP9 expression in HPV-positive patients with CC is higher than in HPV-negative patients (123). Cardeal et al (124) reported that MMP9 is upregulated in human keratinocytes expressing the HPV16 E7 protein. This may be due to TIMP2, an inhibitor of MMP9, which could be downregulated by HPV16 E7. It was also demonstrated that HPV can directly regulate the activity of MMP9 in lung cancer cells (125). There may be an association between HPV infection and the MMP family, which may be beneficial in the diagnosis of cervical precancerous lesions and CC as MMP9 may be considered as a novel biomarker. However, the specific molecular mechanisms require further investigation. As FN1 and ITGB1 were targets of miR-9-3p (126) and FN1 promoted migration and invasion by upregulating MMP9 in cancer (127), it is not surprising that these three genes are co-expressed as a reaction triplet. Furthermore, since higher levels of FN1 and ITGB1 are significantly associated with lower OS rate, these two genes may be developed as novel prognostic factors for CC.

STAT1, the only gene that participates in all four phases, in the present study, is involved in the cancer pathway at phases II, III and IV. It has been reported that STAT1 is upregulated in both CIN1 and CC (128), and the results of the present study that STAT1 is upregulated in phase I, III and IV confirmed this finding. A previous study demonstrated that activated STAT1 plays a tumor suppressive role in breast cancer cells (129). Nevertheless, STAT1 also exerts tumor promoter effects under specific conditions (130). In some malignant diseases, including breast and lung cancers, STAT1 can act as an oncoprotein or a tumor suppressor of the same cell type based on the specific genetic background (130). In CC, STAT1 may have a protective effect in the early stages of HPV infection but may act as a proto-oncogene during the invasive phase of the disease (128). STAT1 can promote cancer cell death by activating p53 expression, and it plays a role in immunosurveillance, and the inhibition of angiogenesis and metastasis in cancer cells (130); however, STAT1 can also promote tumor invasion and metastasis in chronic inflammation (131). The effect of STAT1 in CC still remains unclear; therefore, further verification is required as it may be a key target for the treatment of CC.

The current study presented several limitations. First, as this was an in-silico study, the identification of DEGs may change with additional data, thus the results of subsequent analyses may change accordingly. Secondly, some genes were excluded to decrease the false-positive rate; however, these genes may also have a vital effect on CC. Thirdly, although several genes associated with HPV infection (TP53TG1, RAC1, PAK2 and LTBP2) were identified in DEGs, HPV infection was not observed in the pathway analyses. This may be due to the fact that an insufficient amount of hub genes were identified, or the genes had a low or moderate effect on the HPV infection pathway.

In conclusion, the present study revealed that CDK1, CCNB1, ITGB1, FN1, MMP9 and STAT1 played different roles in the progression of CC through different signaling pathways. CDK1 and CCNB1 served as regulators of proliferation and differentiation via regulation of the cell cycle and viral tumorigenesis, and initiated CIN and CC, whereas FN1, ITGB1 and MMP9 were co-expressed as a reaction triplet to trigger metastasis via cancer pathways, PI3K-Akt signaling pathway and focal adhesion. FN1 and ITGB1 may be novel prognostic factors for CC. STAT1 may have a protective effect in the early stage of HPV infection, but may also act as a proto-oncogene during the invasive stage; however, the specific molecular mechanisms require further investigation.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

All data generated and/or analyzed during this study are included in this published article. The datasets generated and/or analyzed during the current study are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/) and LinkOmics (http://www.linkedomics.org/) repository.

Authors' contributions

YY, YF and WZ designed the present study and drafted the initial manuscript. KW and YL performed the literature review, acquired the data and performed the statistical analyses. All authors have read and approved the 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.

Glossary

Abbreviations

Abbreviations:

References