CDK1-PLK1/SGOL2/ANLN pathway mediating abnormal cell division in cell cycle may be a critical process in hepatocellular carcinoma

ABSTRACT

This study aims to investigate the potential mechanisms and identify core biomarkers of Hepatocellular carcinoma (HCC). The profile GSE113850 was downloaded to analyze the differentially expressed genes. Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and protein–protein interaction network analysis were used to reveal the main signal pathways of the differentially expressed genes (DEGs) and hub genes. The correlation between core gene expression and pathological stages, the disease-free survival analysis, the overall survival analysis were analyzed by Gene Expression Profiling Interactive Analysis. Furthermore, we reidentified the expression level of core genes of carcinoma tissues and para-carcinoma tissues from 14 HCC patients with real-time reverse transcription-polymerase chain reaction analysis (RT-PCR) and western blotting. After SK-Hep1 cell was treated with cyclin-dependent kinase 1 (CDK1) siRNA for 72 h, we detected the expression of the core genes and fluorescence-activated cell sorting analysis. A total of 378 DEGs were found. GO and KEGG analysis revealed that the DEGs were mainly enriched in the cell cycle. There were positive correlations among CDK1, polo-like kinase 1, shugoshin2 and anillin actin-binding protein. Moreover, the expression levels of four core genes were related to the HCC occurrence, pathological stages, and survivorship curve. The clinical HCC specimens verified the higher expression level of core genes by real-time RT-PCR. The transfection of siCDK1 in SK-Hep1 resulted in a disordered cell cycle. Furthermore, CDK1 knockdown suppressed the expression of PLK1, ANLN, and SGOL2. The CDK1-PLK1/SGOL2/ANLN pathway mediating abnormal cell division in the cell cycle might be a critical process in HCC.

KEYWORDS:

• Hepatocellular carcinoma
• CDK1
• PLK1
• SGOL2
• ANLN
• bioinformatics

1. Introduction

Hepatocellular carcinoma (HCC) is a representative of malignancies, which is characterized by poor prognosis and high mortality. Globally, it is the second leading cause of cancer-related mortality [1], presumptively more than 750,000 new HCC cases are diagnosed worldwide every year, 55% of which cases are in China [2] Despite advances in medical and surgical techniques, it is typically in an advanced stage at diagnosis and leaves very little hepatic functional reserve. As a highly heterogeneous tumor, HCC is a progressive accumulation of plenty of alterations of the genome and multifaceted molecular pathway [3]. The identification of target genes and related pathway is of great importance, not only for early diagnosis, optimal treatment, and better prognosis but also for a comprehensive understanding of the cellular and molecular mechanisms involved in carcinogenesis.

Microarray analysis is increasingly valued as an essential tool in the field of medical oncology with the superiority of simultaneously detecting the alterations of thousands of genes, which can bring in great clinical applications, such as molecular classification of cancers, molecular diagnosis, and prognosis prediction [4,5]. In the present study, original data (GSE113850) were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) and we systematically selected the differentially expressed genes (DEGs). We established protein–protein interaction (PPI) network of DEGs and top 20 genes were analyzed with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. Top four core genes, including CDK1, PLK1, SGOL2, and ANLN, were further analyzed in the perspective of disease-free survival and overall survival by the Gene Expression Profiling Interactive Analysis online database (http://gepia.cancer-pku.cn/) and the associations among the top genes were also identified there. To verify the data from bioinformatic analysis, we examined the expression of the four core genes in carcinoma tissues and para-carcinoma tissues of HCC patients. Furthermore, we investigated the mechanisms of CDK1 in HCC progression, and the data suggest that CDK1 promoted HCC invasion through PLK1/SGOL2/ANLN mediating cell cycle disorder, and the abnormal inaccurate chromosome segregation and cytokinesis were involved in the pathological process.

2. Materials and methods

2.1. Microarray data

Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database in the National Center for Biotechnology Information (NCBI) was used to store gene expression datasets, original series, and platform records. The gene expression profile of GSE113850 was downloaded from the GEO database, which was based on Invitrogen NCode Human lncRNA Array on the Platform of GPL16847 by Polycarpou-Schwarz M et al. A total of 39 samples, including 7 normal liver samples and 32 HCC samples, were included in the analysis.

2.2. Identification of DEGs

GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/), an interactive online analysis tool for the GEO database, was applied to identify the DEGs between HCC samples and normal liver samples. The adjusted P-values (adj P-values) were applied to correct for the occurrence of false-positive results using the Benjamini and Hochberg false discovery rate method by default. The logFC <−2 (downregulated genes) or logFC >2 (upregulated genes) and the adj P-value <0.05 were used as the inclusion criteria [6]. Based on that, totally 378 DEGs were found, including 142 upregulated genes and 236 downregulated genes, and we selected the top 20 genes with a high degree of connectivity as hub genes.

To show the results more directly and visually, we used visual hierarchical cluster analysis to depict the heat map and volcano plot of two groups by ImageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/index.html).

2.3. Gene Ontology and KEGG pathway analysis of DEGs and hub genes

Candidate DEGs` functions and pathway enrichments can be analyzed using multiple online databases. Among them, GO (Gene ontology) is a major bioinformatics initiative to unify the representation of gene and gene product and to interpret the characteristic biological attributes including cellular components, molecular function, and biological pathways [7]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism, and the ecosystem from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies [8]. In addition, the 20 hub genes were analyzed by GO and KEGG. The DEGs and 20 hub genes were mapped to the relevant biological annotation using the DAVID online database (https://david.ncifcrf.gov/). P < 0.05 was considered to indicate a statistically significant difference.

2.4. PPI network and module analysis

PPT information of DEGs was acquired from the Search Tool for the Retrieval of Interacting Genes (STRING, http://www.stringdb.org/) database. So far, STRING (version 10.0) contains 9.6 million proteins from 2,031 organisms, and 184 million interactions has been portrayed and integrated [9]. The cutoff criterion was set as a maximum number of interactors = 0 and a confidence score ≥ 0.4. The whole giant network was then uploaded into Cytoscape software (version 3.4) to evaluate primary modules of the PPI sub-network using the plug‑in Molecular Complex Detection (MCODE) according to node score cutoff = 0.2, degree cutoff = 2, max. Depth = 100, and k-core = 2. Also, the genes of three modules were analyzed by GO and KEGG to further investigate the biological function and potential pathway information, respectively.

2.5. Comparison of the hub gene expression level

Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) is a newly created web-based tool for cancer and normal gene expression profiling and interactive analyses. It integrates the tremendous amount of RNA sequencing data produced by large consortium projects such as TCGA and GTEx, in order to provide key interactive and customizable functions including differential expression analysis, profiling plotting, correlation analysis, patient survival analysis, similar gene detection, and dimensionality reduction analysis [10]. In our study, we utilized the box plots to depict the contract of hub gene expression between normal human hepatocytes and hepatoma cells. Meanwhile, the violin plots based on pathological stages were presented to further insight into correlation between gene expressions and disease phases. GEPIA also provides pairwise gene correlation analysis for given genes, and we input several suspicious genes for detailed exploration.

2.6. Survival analysis of core genes

The overall survival and disease-free survival analysis of core genes were carried out on GEPIA. The hazard ratio (HR) with 95% confidence intervals and logrank P value was measured and P < 0.05 was considered statistically significant.

2.7. Reidentification of core gene expression in HCC

2.7.1. Clinical specimens and cell line

The study protocol was approved by the Ethics Committee of Zhongnan Hospital, Wuhan University (No.2019016). We collected 14 patients who underwent curative surgery for HCC in Zhongnan Hospital of Wuhan University, for primary HCC and corresponding adjacent non-tumor (NT) liver specimens. The cell line SK-Hep1 was purchased from the Cell Bank of Chinese Academy of Sciences (No. TCHu109 and GNHu 6, Chinese academy of sciences, Shanghai, China) and cultured in Dulbecco’s modified Eagle’s medium (DMEM; Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, USA). All cells were maintained in a humidified incubator at 37°C with 5% CO₂.

2.7.2. RNA interference

The control nonspecific small interfering RNA and CDK1 siRNAs were synthesized by Guangzhou RiboBio Co., Ltd. Transfection of siRNA was carried out according to the procedure recommended by the manufacturer. The cell line SK-Hep1 was chosen for the RNA interference. The control group was disposed with nonspecific siRNA, and the CDK1 knockdown groups were disposed with CDK1 siRNA. The nucleotide sequences of CDK1 siRNAs were siCDK1 primer: GTACTGCAATTCGGGAAAT. Seventy-two hours after transfection with 50 nM siRNA, the cells were then collected and processed for quantitative real-time PCR and flow cytometry analysis.

2.7.3. RNA isolation, cDNA synthesis, and quantitative real-time PCR

Total RNA was extracted from human liver tissues and cells using TRIzol reagent (No.15596026, Thermo Fisher Scientific, Waltham, MA, USA). The reverse transcription of mRNA was performed using RevertAid RT Reverse Transcription Kit (No.K1691, Thermo Fisher Scientific, Waltham, MA, USA). Quantitative real-time PCR was performed using a standard protocol from the Hieff® qPCR SYBR Green Master Mix (No Rox) kits (No.11201ES08, shanghai, China). β-actin was used as an internal control. Each sample was analyzed in triplicate. Primers specific to target genes were synthesized by Tsing Ke (Beijing Qing Ke Xin Ye Biotechnology Co., Ltd, Wuhan, China). Data were analyzed according to the 2^−ΔCq (Comparative Cq) method. The primers used in this study are as follows: CDK1, forward primer: TTACAAAGATCAAGGGCTGTCCGCAACA, reverse primer: AAAGCCAAGATAAGCAACTCCTTCAGT; PLK1 forward primer: ACGGCAGCGTGCAGATCAACTT; reverse primer: TCCAGGAGACTCAGGCGGTATGT; SGOL2, forward primer: GAAAGCCATTCCCACTCAGACC; reverse primer: TTTGATTCCCAAGATGACCCAC; ANLN, forward primer: CACCGTTTCCATCGTCTC; reverse primer: CTGCTGGTTACTTGCTTCTG; β-actin, forward primer:5ʹ-TGGCACCCAGCACAATGAA-3ʹ, reverse primer: 5ʹ-CTAAGTCATAGTCCGCCTAGAAGCA-3ʹ

2.7.4. Western blotting analysis

Total protein was isolated from SK-Hep-1 cells using RIPA lysis buffer (Wuhan Servicebio Biotechnology, Wuhan, China). Protein concentrations were determined using bicinchoninic acid assay (Wuhan Servicebio Biotechnology). Equal amounts of protein for each sample were run on polyacrylamide gels and then transferred to PVDF membranes. Subsequently, PVDF membranes were incubated at 4°C overnight with primary antibodies at the indicated concentrations (Table 1). Protein bands were detected using a chemiluminescence reagent kit (Wuhan Servicebio Biotechnology). GAPDH content was used to standardize sample loading.

Table 1. Diluted concentrations of antibodies used for western blotting

Antibodies	Catalog number	Diluted concentrations
CDK1	10026-1-AP	1: 1000
PLK1	19553-1-AP	1: 1000
SGOL2	ab122258	1: 1000
ANLN	66643-1-lg	1: 1000
GAPDH	GB12002	1: 1000

2.7.5. Immunohistochemistry (IHC)

For IHC, the carcinoma tissues and para-carcinoma tissues of HCC patients were deparaffinized in xylene and rehydrated in ethanol, and submerged in distilled water. After treated in Hydrogen Peroxide, the slides were processed for antigen retrieval in Tris-EDTA buffer, and then cooled down slowly. After blocking with 10% goat serum, the slides were incubated in primary antibodies overnight and then incubated in HRP-conjugated secondary antibodies. Subsequent detection was performed using the standard substrate detection of HRP. Then, the sections were stained with hematoxylin and dehydration in graded alcohols and xylene. The CDK1 antibody was purchased from (No.10026-1-AP, Proteintech, Wuhan, China). The IHC images were observed under 400× magnification.

2.7.6. Flow cytometry analysis

Cells were transfected with CDK1 siRNAs or control sequence for 72 h before analysis. Briefly, cells were washed in PBS and fixed in 70% precooled ethanol at 4°C at least 4 h. Fixed cells were then washed twice in PBS and water-bath heated in 100ul RNase (50 µg/mL) at 37°C for 30 min, then stained away from light in 50 µg/mL propidium iodide for 30 min. They were analyzed by flow cytometry (FACS Calibur flow cytometer; BD Biosciences, Franklin Lakes, NJ, USA).

2.8. Statistical analysis

All values are presented as the mean ± SD. All statistical analyses were performed by SPSS 19.0 software. A difference of P < 0.05 was considered statistically significant.

3. Results

3.1. Identification of DEGs and hub genes

A total of 32 HCC samples and 7 normal samples were analyzed. By using GEO2R online tool with the inclusion criteria of logFC <−2 or logFC >2 and the adj P-value <0.05, totally 378 DEGs were found, including 142 up-regulated genes and 236 down-regulated genes. The results were shown in the heat map and volcano plot (Figure 1(a,b)). Meanwhile, we selected the top 20 genes with a high degree of connectivity as hub genes, including CDK1, PLK1, HIST3H3, CYP2B6, SYK, ACSL1, ANLN, EDN1, RAD21, CXCL12, IGF1, DMGDH, HDAC10, SHMT1, SGOL2, CYP2 C19, CYP3A5, HTR2 C, SYNE1, and MTR (Table 2).

Figure 1. Identification of DEGs of HCC in GSE113850. (a) Heat map of the expression levels of the top 60 DEGs in GSE113850 employed the GEO2R online analysis tool based on P-value < 0.05 and logFC ≤ −2 or logFC ≥ 2 criteria. Ordinate Abscissa represented sample types and abscissa referred to genes. Each block referred to the expression of a gene in a sample, and upper histogram showed color gradation. Red represented high expression and green represented low expression. (b) Volcano plot about total 378 DEGs from GSE113850. The red dots stand for 142 up-regulated DEGs, while the green dots represent 236 down-regulated DEGs. The blue dots denote the no-differentially expressed genes

Table 2. Top 20 hub genes with a higher degree of connectivity

Name	Degree	P -value	Log FC
CDK1	11	6.79E-04	2.69
PLK1	10	1.61E-04	2.07
HIST3H3	10	1.55E-06	2.35
CYP2B6	7	1.46E-03	−3.31
SYK	6	0.00325	−2.16
ACSL1	6	6.70E-04	−2.06
ANLN	6	2.62E-04	2.2
EDN1	6	0.0341	−2.38
RAD21	5	0.00195	−3.09
CXCL12	5	9.67E-05	−2.52
IGF1	5	8.14E-04	−2.49
DMGDH	5	2.48E-04	−2.69
HDAC10	5	3.79E-05	2.33
SHMT1	4	1.07E-05	−2.99
SGOL2	4	7.50E-03	−2.96
CYP2C19	4	6.51E-03	−2.13
CYP3A5	4	4.67E-03	−2.27
HTR2C	4	3.04E-04	2.01
SYNE1	4	1.81E-04	−2.24
MTR	3	5.58E-05	2.15

3.2. GO Function and KEGG pathway enrichment analysis

To investigate the functions of DEGs and further assess the signal pathway, we performed GO and KEGG pathway enrichment analysis using DAVID. Inputting all the DEGs into DAVID, we obtained up-regulated DEGs and down-regulated DEGs by GO analysis. To be more specific, these DEGs were mainly enriched in biological processes (BP), involving the regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle; the positive regulation of ubiquitin-protein ligase activity involved in regulation of mitotic cell cycle transition; the anaphase-promoting complex-dependent catabolic process and the single organismal cell-cell adhesion for upregulation. GO cell component (CC) analysis uncovered that the down-regulated DEGs were principally enriched in the plasma membrane, chromosome centromeric region, and extracellular region. To be more specific, we imported 20 hub genes into DAVID software. The KEGG result showed that the pathways mainly involved were metabolic pathways and cell cycle (Table 3).

Table 3. KEGG pathway analysis of 20 hub gene

Term	Count	%	P-value	Genes
bta00140: Steroid hormone biosynthesis	3	0.107953	0.005474	CYP3A5, CYP2B6
bta00830: Retinol metabolism	3	0.107953	0.005474	CYP3A5, CYP2B6
bta05204: Chemical carcinogenesis	3	0.107953	0.008162	CYP3A5, CYP2B6
bta04914: Progesterone-mediated oocyte maturation	3	0.107953	0.012674	CDK1, PLK1, IGF1
bta04114: Oocyte meiosis	3	0.107953	0.020026	CDK1, PLK1, IGF1
bta04110: Cell cycle	3	0.107953	0.024602	CDK1, RAD21, PLK1
bta01100: Metabolic pathways	7	0.251889	0.02564	SHMT1, CYP3A5, ACSL1, CYP2B6, MTR, DMGDH
bta00670: One carbon pool by folate	2	0.071968	0.035052	SHMT1, MTR
bta00591: Linoleic acid metabolism	2	0.071968	0.070805	CYP3A5
bta05203: Viral carcinogenesis	3	0.107953	0.076482	CDK1, HDAC10, SYK
bta00260: Glycine, serine and threonine metabolism	2	0.071968	0.078173	SHMT1, DMGDH

3.3. Module screening from the PPI network

Based on the information of the STRING protein query from public databases, we constructed the PPI network of all DEGs (Figure 2(a)). Furthermore, we applied the MCODE plug-in to focus on the most significant modules in this PPI network. The selected top three modules are shown in Figure 2(b–d).

Figure 2. PPI network and module analysis of HCC. (a) PPI network of all DEGs via the information of the STRING protein query with a maximum number of interactors ≤ 5, confidence score ≥ 0.4 which consists of 88 nodes and 112 edges based on STRING database. The top 10 DEGs with high degree of connectivity were selected as the hub genes of HCC. (b-d) PPI networks of Module 1–3. The top three gene modules were significantly enriched in metabolic pathways, mineral absorption, and cell cycle, respectively

What is more, we conducted GO and KEGG analysis of the three top modules, respectively, to reveal the predominant biological functions and pathways. KEGG analysis disclosed that the top three modules were mainly associated with metabolic pathways and mineral absorption (Tables 4 and 5).

Table 4. Gene ontology analysis of hub genes associated with hepatocellular carcinoma (Cluster1, 2, and 3)

	Category	Term	Count	%	P-value
Cluster1	GOTERM_BP_DIRECT	GO:0019373~epoxygenase P450 pathway	2	0.99	0.00
	GOTERM_MF_DIRECT	GO:0020037~heme binding	4	1.98	0.00
	GOTERM_MF_DIRECT	GO:0005506~iron ion binding	4	1.98	0.00
	GOTERM_MF_DIRECT	GO:0008392~arachidonic acid epoxygenase activity	2	0.99	0.00
	GOTERM_MF_DIRECT	GO:0016712~oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen	2	0.99	0.00
	GOTERM_MF_DIRECT	GO:0008395~steroid hydroxylase activity	2	0.99	0.00
Cluster2	GOTERM_BP_DIRECT	GO:0071276~cellular response to cadmium ion	4	19.05	0.00
	GOTERM_BP_DIRECT	GO:0045926~negative regulation of growth	4	19.05	0.00
	GOTERM_BP_DIRECT	GO:0071294~cellular response to zinc ion	4	19.05	0.00
	GOTERM_CC_DIRECT	GO:0048471~perinuclear region of cytoplasm	4	19.05	0.00
	GOTERM_CC_DIRECT	GO:0005737~cytoplasm	4	19.05	0.02
	GOTERM_CC_DIRECT	GO:0005634~nucleus	4	19.05	0.03
	GOTERM_MF_DIRECT	GO:0008270~zinc ion binding	4	19.05	0.00
	GOTERM_MF_DIRECT	GO:0046872~metal ion binding	4	19.05	0.00
Cluster3	GOTERM_BP_DIRECT	GO:0018105~peptidyl-serine phosphorylation	2	0.23	0.03
	GOTERM_BP_DIRECT	GO:0043066~negative regulation of apoptotic process	2	0,23	0.06
	GOTERM_CC_DIRECT	GO:0005813~centrosome	3	0.34	0.00
	GOTERM_CC_DIRECT	GO:0005876~spindle microtubule	2	0.23	0.00
	GOTERM_CC_DIRECT	GO:0000922~spindle pole	2	0.23	0.02
	GOTERM_CC_DIRECT	GO:0030496~midbody	2	0.23	0.02

Table 5. KEGG pathway analysis of hub genes associated with hepatocellular carcinoma (Cluster1, 2, 3)

Category	Term	Count	%	P-value	Genes
Cluster1	bta00830:Retinol metabolism	3	1.49	0.000	CYP3A5, CYP2B6
	bta00140:Steroid hormone biosynthesis	3	1.49	0.000	CYP3A5, CYP2B6
	bta05204:Chemical carcinogenesis	3	1.49	0.000	CYP3A5, CYP2B6
	bta01100:Metabolic pathways	4	1.98	0.004	CYP3A5, ACSL1, CYP2B6
	bta00260:Glycine, serine and threonine metabolism	2	0.99	0.015	CYP3A5
Cluster2	hsa04978:Mineral absorption	4	19.05	0.000	MT1E, MT1X, MT1G, MT1F
Cluster3	cfa04914:Progesterone-mediated oocyte maturation	2	0.226	0.025	CDK1, PLK1
	cfa04114:Oocyte meiosis	2	0.226	0.031	CDK1, PLK1
	cfa04110:Cell cycle	2	0.226	0.036	CDK1, PLK1

3.4. The Kaplan–Meier plotter and expression level of core genes

As a highly heterogeneous tumor, the pathogenesis of HCC is intricate, which is implicated with the progressive accumulation of plenty of alterations of genome and multifaceted molecular pathways. Based on the comprehensive but complicated results available of GO and KEGG analysis, we can hardly focalize on a pivotal pathway which is partly responsible for cancerization. GEPIA provides pairwise gene correlation analysis for given genes, and we input the hub genes for detailed exploration. Intriguingly, we found distinct positive correlations between CDK1, PLK1, SGOL2, and ANLN (Figure 3).

Figure 3. Correlation analysis of any two core genes in GEPIA. (a) CDK1 and PLK1 correlation analysis. CDK1 and PLK1 are obviously positively correlated (P-value = 0, R = 0.85). (b) CDK1 and SGOL2 correlation analysis. CDK1 and SGOL2 are obviously positively correlated (P-value = 0, R = 0.82). (c) CDK1 and ANLN correlation analysis. CDK1 and ANLN are obviously positively correlated (P-value = 0, R = 0.88). (d) PLK2 and SGOL2 correlation analysis. PLK2 and SGOL2 are obviously positively correlated (P-value = 0, R = 0.88). (e) PLK1 and ANLN correlation analysis. PLK1 and ANLN are obviously positively correlated (P-value  = 0, R = 0.88)

Cyclin-dependent kinase 1, also known as CDK1 or cell division cycle protein 2 homolog, a highly conserved protein that functions as a serine/threonine kinase, is a key player in cell cycle regulation [11]. Serine/threonine-protein kinase PLK1, also known as polo-like kinase 1 (PLK1) or serine/threonine-protein kinase 13 (STPK13), functions as an early trigger for G2/M transition. Often observed to overexpress in tumor cells, PLK1 is considered as a proto-oncogene [12]. Depletion of PLK1 by RNA interference results in metaphase arrest and formation of abnormal chromatin structures. Knockdown of PLK1 also reduces cell survival and elevates drug sensitivity of tumor, but not normal cells [13,14]. In more detail, PLK1 contributes to the tumorigenicity of BEL-7402 hepatoma cells via regulation of Survivin expression [15]. Shugoshin-like 2 (SGOL2) is vital for faithful chromosome segregation during mitosis and meiosis, which regulates sister chromatid cohesion by protecting the integrity of a complex named cohesion [16]. ANLN (Anillin) is an evolutionarily conserved actin-binding protein that serves as a key mediator of cytokinesis [17]. It binds to all three types of actin filaments, as well as other important cytoskeletal components, serving as a scaffolding protein to coordinate cytokinesis in space and time [18,19]. The high correlations hinted that they might participate jointly in a pro-tumorigenic or tumor-responsive pathway. According to the findings, we chose the four genes as core genes and tried to undercover the key pathways.

The prognostic information of the core genes was downloaded from http://gepia.cancer-pku.cn/. As shown in Figure 4(a–d), the results of disease-free survival verified the worse effect of CDK1 (HR (high) = 1.7, P (HR) = 0.00063, logrank P = 0.00057), PLK1 (HR (high) = 1.6, P (HR) = 0.0019, logrank P = 0.0019), ANLN (HR (high) = 1.7, P (HR) = 5e-04, logrank P = 0.00046) and SGOL2 (HR (high) = 1.7, P (HR) = 0.00058, logrank P = 0.00052). What is more, the high expression of CDK1 was associated with worse overall survival (OS) for HCC patients (HR (high) = 2, P (HR) = 0.00022, logrank P = 0.00017), as well as PLK1 (HR (high) = 1.8, P (HR) = 0.00097, logrank P = 0.00082), SGOL2(HR (high) = 1.9, P (HR) = 0.00031, logrank P = 0.00025) and ANLN (HR (high) = 1.8, P (HR) = 0.00099, logrank P = 0.00085). (Figure 4(e–h))

Figure 4. The Kaplan-Meier plotter of four core genes in HCC patients. The disease-free curves of CDK1 (a), PLK1 (b), SGOL2 (c), and ANLN (d). It was demonstrated that expression of CDK1 (HR (high) = 1.7, Logrank P =0.00057), PLK1 (HR (high) = 1.6, Logrank P = 0.0019), SGOL2 (HR (high) = 1.7, Logrank P = 0.00052) and ANLN (HR (high) = 1.7, Logrank P = 0.00046) were associated with worse disease-free survival (DFS) for HCC patients. The overall survival curves of CDK1 (e), PLK1 (f), SGOL2 (g), and ANLN (h) in HCC patients. It was demonstrated that expression of CDK1 (HR (high) = 2, Logrank P = 0.00017), PLK1 (HR (high) = 1.8, Logrank P = 0.00082), SGOL2 (HR (high) = 1.9, Logrank P = 0.00025) and ANLN (HR (high) = 1.8, Logrank P = 0.00085) were associated with worse overall survival (OS) for HCC patients. *P < 0.05 was regarded statistically different

In order to compare the expression level of hub genes in normal hepatocytes with that in hepatoma cells, we drew the boxplots which displayed that CDK1, PLK1, SGOL2, and ANLN appeared a higher expression level, which is consistent with the survival analysis (Figure 5(a–d)). The violin plots showed the relevance between gene expressions and pathological stages, indicating the influence of core genes on progression of disease (Figure 5(e–h)).

Figure 5. The boxplots and the violin plots of four core genes in HCC patients. The boxplots of CDK1 (a), PLK1 (b), SGOL2 (c), and ANLN (d). The expression levels of CDK1, PLK1, SGOL2, and ANLN were higher in hepatocellular carcinoma tissues than in liver tissues. The violin plots of CDK1 (e), PLK1 (f), SGOL2 (g), and ANLN (h) in HCC patients. The expression levels of CDK1, PLK1, SGOL2, and ANLN were significantly associated with tumor stage for hepatocellular carcinoma. *P < 0.05 was regarded statistically different

3.5. Reidentification of target gene expression in HCC

To reidentify the expression profiles of core genes shown in the integrated bioinformatic analysis, we examined the mRNA levels in 14 paired HCC samples. The results showed significantly higher CDK1, PLK1, ANLN, and SGOL2 expression in tumor tissues than those in adjacent nontumor tissues (Figure 6).

Figure 6. The relative mRNA level of four hub genes in HCC tissues and adjacent non-tumor tissues. (a) Up-regulation of CDK1 (log₂(HCC/NT)>0) was detected in 100% (14/14) of primary HCCs, and (b) up-regulation of PLK1 was detected in 85.71% (12/14), (c) up-regulation of SGOL2 was detected in 85.71% (12/14), (d) up-regulation of ANLN was detected in 100% (14/14). The results showed that the expression level of CDK1, PLK1, SGOL2, and ANLN was distinctly higher in HCC tissues

Immunohistochemistry (IHC) in three pairs of HCC tissues further confirmed the upregulation of CDK1 protein in HCC cells and indicated that multinuclear giant cells causing by failed cell division were involved in HCC pathogenesis (Figure 7).

Figure 7. Representative IHC images showing the protein level of CDK1 in HCC tissues and paired adjacent non-tumor tissues. The results shown that CDK1 protein was significant high expression in the HCC issues (d–e) compared with adjacent non-tumor tissues (a–c) and obvious delayed process of the mitosis was found in those HCC issues

3.6. CDK1 mediates disordered cell cycle

To investigate the potential biological function of CDK1, we knocked down CDK1 in SK-Hep1 cells. The siRNAs were transiently transfected into the cells and the knockdown efficiency was measured by real-time RT-PCR (Figure 8(a)). Previous studies suggest that CDK1 is implicated in cell cycle regulation, which is confirmed by our results. Fluorescence-activated cell sorting (FACS) analysis showed that the proportion of cells in G2/M phase increased while G0/G1 decreased accordingly following CDK1 knockdown using siRNAs (Figure 8(b–f)).

Figure 8. FACS analysis showing cell cycle arrest in G2/M phase after silencing of CDK1 in transiently transfected cell lines. (a) The siRNA knockdown efficiency of CDK1 was measured by quantitative real-time RT-PCR. (b) showed the FACS analysis of the control group; (c) showed the FACS analysis of the siCDK1 group. (d-f) showed that the proportion of cells in G2/M phase increased while G0/G1 decreased accordingly following CDK1 knockdown. *P < 0.05 was regarded statistically different

3.7. CDK1 knockdown suppresses the expression of downstream genes

PLK1, ANLN, and SGOL2 are implicated to involve in different aspects of cell division, and our results indicated that the effect of CDK1 on cell cycle was mediated by regulating downstream gene expressions including PLK1, ANLN, and SGOL2. The results from quantitative real-time RT-PCR showed the significantly decreased expression levels of PLK1, ANLN, and SGOL2 following CDK1 knockdown (Figure 10) and CDK1 knockdown suppressed the expression protein of PLK1, ANLN, and SGOL2 (Figure 9) in SK-Hep1 cell line.

Figure 9. CDK1 knockdown suppressed the expression protein of hub genes in transiently transfected cell lines. *P < 0.05 was regarded statistically different

Figure 10. CDK1 knockdown suppressed the relative mRNA level of hub genes in transiently transfected cell lines. *P < 0.05 was regarded statistically different

4. Discussion

In this study, a total of 378 DEGs, including 236 down-regulated and 142 up-regulated genes, were identified in HCC tissues compared to in normal livers. To have a deep exploration of these DEGs, we carried out GO function and KEGG pathway analysis of these DEGs. To confirm the clinical significance and practicability, we reexamined the expression level of core genes in clinical HCC specimens. Cell experiments worked out the specific pathway in which the core genes get involved in HCC occurrence and development.

4.1. CDK1/PLK1/SGOL2/ANLN genetic interaction mediates inaccurate chromosome segregation through cohesin and abnormal cytokinesis in mitosis

In mitosis, the genome is programmed to be duplicated and separated in the way that each of the daughter cells inherits a complete copy of the genome without alterations. Chromosome segregation and cytokinesis are essential for complete cell division and are processes shared by most mitotic cells [20]. The processes are controlled by numerous mechanisms ensuring correct cell division [21]. Any errors in the signal pathway are likely to initiate aberrant cell proliferation and facilitate the development of cancer [22]. A key feature for controlling genomic integrity is sister chromatid cohesion, which is mediated by the cohesin complex.

In this study, CDK1, PLK1, SGOL2, and ANLN, four of the selected hub genes, were confirmed to be associated with HCC pathological stages and disease-free survival. More importantly, the following correlation analysis exhibited high positive correlation of expression level between each pair, implying that potential genetic interactions between CDK1/PLK1/SGOL2/ANLN mediates disorganized cell cycle, exerting an influence on HCC progression and prognosis. The results of CDK1 knockdown experiments implied that PLK1, SGOL2, and ANLN were downstream genes regulated by CDK1.

Binding of CDK1 to cyclin B is essential for activation of the kinase, and the activation of the CDK1-cyclin B brings the onset of mitosis and is tightly regulated [23]. At the G2/M transition, CyclinB1 expression gradually increases and binds to CDK1 form a heterodimer. CDK1-cyclin B is activated through dephosphorylation on threonine 14 (Thr14) and tyrosine 15 (Tyr15) by Cdc25 and triggers the G2/M transition and initiates the biological processes and changes in cellular structure and organization necessary for mitotic progression [24]. Mouse knockout experiments have shown that CDK1 is prerequisite for mammalian cell proliferation; it is the only CDK that can initiate the onset of mitosis (that is, the M phase) [25], subsequently, CDK1-cyclin B phosphorylates multiple substrates, leading to nuclear envelope breakdown, chromosome condensation, and mitotic spindle assembly [26]. Analysis of human HCC tissue samples confirmed that CDK1 activity was detected in primary malignancies but not in healthy paraneoplastic tissues, what is more, small hairpin RNA (shRNA) knockdown of CDK1 significantly reduced the tumor-specific killing effects of apoptin [27].

PLK1 is a master mitotic regulator controlling a wide variety of processes during G2/M cell cycle progression, such as centrosome maturation, entry into mitosis, chromosome segregation, and cytokinesis. In prophase, PLK1 shuttles to the kinetochores, helping condensed chromosomes to align in the metaphase plate. At the onset of anaphase, PLK1 localizes at the spindle mid-zone and later to the cytokinetic bridge in order to coordinate cytokinesis and cell abscission [12,28]. As mentioned above, CDK1-cyclin B phosphorylates multiple substrates, including the conserved protein, Bora, which functions as a PLK1 activator. Phosphorylation of specific sites on Bora by CDK1 is essential for PLK1 function in mitotic entry in C. elegans embryos and during DNA damage checkpoint recovery in mammalian cells [29,30]. Cyclin-B accumulation in G2 triggers activation of Cdk1 and drives cells into mitosis; once all the chromosomes have been correctly aligned on the mitotic spindle, the activity of Cdk1 falls and the activation of CDK1 control the switch of PLK1 localization from centrosomes and kinetochores during metaphase to the central spindle during anaphase [31]. At the G2/M transition, cohesin is phosphorylated by the mitotic kinases CDK1, PLK1, and Aurora B. Then, Wapl will remove the bulk of cohesin from chromosome arms and the sister chromatids move to the opposite poles of the cell [24]. In accordance with this, the GEO analysis and quantitative real-time RT-PCR results showed the consistent uptrend of CDK1 and PLK1 expression in HCC tissues, and higher expression of both indicated a shorter overall survival and disease-free survival.

SGOL2 is vital to regulate sister chromatid cohesion by protecting the integrity of a complex named cohesin. Cohesion, a multi-subunit complex, consisting of two SMC (structural maintenance of chromosome) family proteins – a kleisin subunit Scc1/Rad21 and an accessory subunit Scc3/SA (stromal antigen) [16,32,33], is maintained until metaphase when sister kinetochores attach to microtubules emanating from the opposite spindle poles [34]. SGOL2 regulates chromosome alignment and retention of centromeric protection mainly by recruiting MCAK and PP2A to centromeres, and that the recruitment of these factors largely depends on SGOL2 phosphorylations by Aurora B in Hela cells [35]. There is a report that Aurora B phosphorylates PLK1 on Thr210 and activates PLK1 at centrosomes, active PLK1 and phosphorylates MCAK at Ser715 which promotes its microtubule depolymerase activity essential for faithful chromosome segregation [36]. In prophase, phosphorylation of the cohesin subunit SA by PLK1 promotes the sister chromatid dissociation [37,38], and phosphorylation of the Scc1 subunit by PLK1 enhances its cleavability by separase [39]. Tomoya S. Kitajima et al. revealed the mechanism of protection of cohesion against phosphorylation: SGOL1 associated with protein phosphatase 2A (PP2A) plays a crucial part in counteracting the phosphorylation by PLK1, while SGOL2 is required for the localization of PP2A at centromeres [40]. Hence, the expression of SGOL2 is crucial to the validity and accuracy of cell cycle, both in mitosis and meiosis. Further investigation into the activation mechanism of SGOL1 uncovered that SGOL1 is activated by CDK1-dependent phosphorylation at T346 during early mitosis in human cells, allowing SGOL1 binding to cohesion [41]. The interrelation between SGOL2 and CDK1 is unclear yet, but the high correlation of SGOL2/CDK1 expression indicates an underlying relation. What intrigued us is that disruption of the mouse SGOL2 gene leads to no alterations in sister chromatid cohesion in embryonic cultured fibroblasts or adult somatic tissues. What is more, SGOL2 knockout mice grow normally and survive to adulthood without any abnormality, while both male and female SGOL2-deficient mice are infertile [42]. However, in our study, the GEO analysis manifested that SGOL2 expresses at a statistically significant higher level in HCC tissues compared to normal liver tissues. We speculated that SGOL2-deficiency is unable to cause disordered cell cycle or induce carcinogenesis independently, but its disfunction or high expression is a significant part in the pathological process of HCC. The Kaplan-Meier plotter showed that SGOL2 high-expression is associated with shorter disease-free survival, indicating a worse prognosis. Inaccuracy of chromosome segregation in mitosis gives rise to aneuploidy, which may lead to tumorigenesis under certain contexts. Indeed, errors of cohesin and its regulators are proved to cause human diseases, including cancer [43]. A deep exploration into cohesion stability and dissociation will help to decipher the latent mechanism of carcinogenesis and identify new target genes for early diagnosis or treatment.

Multinuclear cells caused by cytokinesis failure are a common characteristic amongst cancers, and unscheduled polyploidization may cause genomic instability that might lead to neoplastic aneuploidy [44]. The actin-binding protein ANLN is a ubiquitously expressed protein required for cytokinesis. During the interphase of the cell cycle, ANLN is primarily located to the nucleus. At the onset of mitosis, ANLN protein relocates to the cytoplasm where it accumulates in the contractile ring and cleavage furrow during telophase [45]. ANLN knockdown experiment demonstrated that ANLN was required for cytokinesis, for that siANLN-treated cells failed to complete cytokinesis and became large and multinucleated. Further study showed that siANLN-treated cells still appeared to be frozen in anaphase and could not generate a cleavage furrow or contractile ring [46]. There is reported that Anln shRNA inhibits tumor suppressor but not oncogene by affecting the polyploidy in the liver [47]. Magnusson et al. revealed that transient knockdown of ANLN protein expression in breast cancer cell lines resulted in an increase of senescent cells and an accumulation of cells in the G2/M phase of the cell cycle with altered cell morphology including large, poly-nucleated cells [48]. The IHC figures exhibited distinct multinuclear giant cells in HCC tissues. The exact mechanisms by which CDK1 influences ANLN remain unknown, but its presence is important during multiple stages of cytokinesis. Previous research mentioned that inhibition of phosphoinositide 3-kinase/AKT activity in NSCLC (non-small-cell lung cancer) cells decreased the stability of ANLN and caused a reduction of the nuclear ANLN level, which might give us a clue that PI3 K/AKT pathway links up CDK1 and ANLN [49].

4.2. The disordered cell cycle is critical in the carcinogenesis of HCC

Although a multitude of mechanisms have been proposed by bioinformatic analysis, disordered cell cycle is the predominant contributing to the progression of HCC [50,51]. The inerrant cell cycle is precisely regulated throughout the global process, and ultimately determines the fate of organisms. In addition to intrinsic mechanisms, the cell cycle is very sensitive to the extracellular environment and can be altered by factors including growth factors, nutrients, and other stress signals [52]. The FACS analysis showed that the proportion of cells in G2/M phase increased while G0/G1 decreased accordingly following CDK1 knockdown induced by siRNAs.

As is stated above, PLK1/SGOL2 mediated inaccurate chromosome segregation and ANLN mediated abnormal cytokinesis regulated by CDK1 are likely to trigger disordered cell division, thus interfering inherent cell cycle.

What is more, KEGG analysis shows that the downregulated DEGs are implicated with oxidative stress. When the levels of reactive oxygen species (ROS) and reactive nitrogen species (RNS) increase beyond norm, intracellular oxidative stress response is triggered, thus activating a variety of mechanisms to allow them to cope with these changes, which may cause cell injury and lead to initiation and progression of cancer by regulating molecules such as DNA, enhancers, transcription factors, and cell cycle regulators [53]. Chul-Woong Pyo et al. [54] found that exposure of cells to excessive H₂O₂ induced ubiquitin-dependent proteasomal degradation of cyclin D1, subsequently followed by translational repression. The dual control mechanism was found to contribute to the induction of cell cycle arrest in G2 phase under oxidative stress [55].

The functional analysis of hub genes and DEGs approved metal homeostasis as an enrichment pathway involved in HCC. Correct metal homeostasis is essential for the control of physiology in eukaryotic cells, most of which act as a cofactor for many enzymes and excess of which may be toxic to cells [56]. For instance, via preventing the formation or activation of the cyclin and cyclin-dependent kinase complexes, iron acts as an indispensable regulator of cell cycle. Studies showed that intracellular iron disruption by chelators could cause cell cycle arrest, particularly in G1 and S phases [54,57]. On the other hand, excess intracellular iron may result in oxidative stress response, which can eventually lead to genome instability and disordered cell cycle [58,59]. Several reports have given evidence that alterations in the import, export, and storage of cellular iron may contribute to cancer development, behavior, and recurrence [56].

Taken all together, defect of the cell cycle may intermediate cell transformation and uncontrolled cell growth, dominating the pathological process of malignant tumors. Increased oxidative response and disequilibrium of metal homeostasis make an adverse effect on cell cycle directly or indirectly, subsequently facilitating the carcinogenesis. Therefore, therapies targeting the cell cycle of cancer cells may be a potential strategy for therapy of HCC. The new medical pattern, in which bioinformatic analysis provides clues of potential target genes and possible pathways related to diseases, and clinical specimens and cell experiments offer mechanistic investigation and verification, will facilitate future research on the molecular mechanisms contributing to HCC and promote the discovery of novel diagnosis biomarkers and therapeutic strategies.

Disclosure statement

The authors declare that they have no conflict of interests to state.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This project was supported by the Open Research Fund Program of the State Key Laboratory of Virology of China (No. 2018KF005), the Medical Science Advancement Program (Clinical Medicine) of Wuhan University (No. TFLC 2018003), the National Natural Science Foundation of China (No. 81673503 and 30973582) and the Innovative Experimental Project of Wuhan University Medical Department’s for College Students in 2018 (No. MS2018028 and MS2018029).