Centromere protein N may be a novel malignant prognostic biomarker for hepatocellular carcinoma

Microarray data

The gene expression profiles of the GSE87630, GSE112790 and GSE25097 datasets were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (Edgar, Domrachev & Lash, 2002). GSE87630 contains 64 HCC samples and 30 normal samples and is based on the GLP6947 platform (Illumina HumanHT−12 V3.0 Expression BeadChip) (Woo et al., 2017); GSE112790 contains 183 HCC samples and 15 normal samples and is based upon the GLP570 platform ((HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array) (Shimada et al., 2019); and GSE25097 contains 268 HCC samples, 243 adjacent nontumor samples, and 40 cirrhotic liver samples, and also 6 healthy liver samples and is based upon the GPL10687 platform (Rosetta/Merck Human RSTA Affymetrix 1.0 microarray, Custom CDF) (Tung et al., 2011).

Identification of DEGs

DEGs between HCC samples and normal samples were identified utilizing GEO2R (Barrett et al., 2013), a tool for the online analysis of differential genes between datasets in a GEO series, with the following thresholds: | logFC | >1 and a false discovery rate (FDR) <  0.05. Volcano plots were constructed with GraphPad Prism 8 (GraphPad Software Inc.). Heatmaps were drawn with the heatmap R package (Galili et al., 2018). Overlapping DEGs in the two datasets were detected with the VennDiagram R package (Chen & Boutros, 2011).

Functional enrichment analysis

Gene Ontology (GO) as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were mainly utilized to explore the functions of these identified DEGs. GO consists three components: biological process (BP), cellular component (CC) as well as molecular function (MF) (Ashburner et al., 2000). KEGG mainly evaluates the pathways in which DEGs may be involved (Kanehisa & Goto, 2000). The results were visualized using clusterProfiler and ggplot2 R package (threshold: P <  0.05) (Ito & Murphy, 2013; Yu et al., 2012).

Protein-protein interaction (PPI) network construction

The PPI network comprising all the DEGs was predicted and analyzed by utilizing the Search Tool for Retrieval of Interacting Genes (STRING) database (http://string-db.org) (Szklarczyk et al., 2015) S. A combined score >0.7 was taken into consideration significant. The PPI network was used to clarify mechanisms relevant to the event and growth of HCC.

Hub gene selection

The public bioinformatics software Cytoscape is a system for visualizing complicated networks and incorporating related information (Smoot et al., 2011). The cytoHubba plugin application of Cytoscape was utilized to screen hub genes in a network through 11 topological analysis approaches (Chin et al., 2014). In this study, the intersecting sets of the top 30 genes revealed with the maximal clique centrality (MCC) and density of maximum neighborhood component (DMNC) methods were visualized as a Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/). Molecular Complex Detection (MCODE) was utilized to screen within the PPI network with the following parameters: degree cutoff >  2, node score cutoff >  0.2, K-core >  2 and max dep >  100 (Gary & Bader, 2003)). Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/index.html), an online tool that consists of information from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases, was utilized to analyze the expression of the genes as well as their association with patient survival and clinical stage (Tang et al., 2017). Then, after combining the statistical results of gene expression from the Oncomine (http://www.oncomine.com) (Rhodes et al., 2004), EMBL-EBI (https://www.ebi.ac.uk) and CCLE (https://www.broadinstitute.org/ccle) databases, the statistically significant hub genes were screened. EMBL-EBI, a web-based tool, provides a series of bioinformatics applications for sequence analysis (Li et al., 2015). CCLE contains genomic data, analytical data, and visualization data for approximately 1000 cell lines (Barretina et al., 2012). The levels of expression of the hub genes in clinical samples were compared with those in the GSE87630, GSE112790 and GSE25097 datasets using GraphPad Prism 8 (GraphPad Software Inc.). The area under the curve (AUC) assessed by the pROC R package was used to validate the prognostic efficacy of the hub genes in liver cancer (Robin et al., 2011).

Database analysis of the differential expression and prognosis of the CENPN gene

CENPN gene expression in the TCGA liver cancer database (TCGA-LIHC) was extracted. The limma R package was utilized to analyze the different scatter plots, and the survival R package was used to draw survival curves, and also the limma and ggpubr R packages were utilized to analyze CENPN expressing in the different clinical stages (Holleczek & Brenner, 2013; Ritchie et al., 2015).

Gene set enrichment analysis (GSEA)

The TCGA-LIHC gene expression dataset was downloaded and installed from the Xena web browser (https://xenabrowser.net/). Two groups from the TCGA-LIHC dataset were classified based on the median CENPN expression level (high and low CENPN). The GSEA 4.1.0 software program, which was downloaded from http://www.broad.mit.edu/gsea/, was utilized to carry out GSEA with the following predefined gene sets (Subramanian et al., 2005): hallmark gene sets, curated gene sets and gene ontology. The permutation number was established as 1000, and an FDR <0.25 was considered significant.

Cell culture and cell transfection

The human HCC cell lines HepG2 and Huh7 were purchased from Procell Life Science & Technology Company (Wuhan, China) and confirmed to contain no mycoplasma contamination via STR analysis. Cells were cultured in DMEM supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin at 37 °C and 5% CO₂. A small interfering RNA (siRNA) targeting CENPN and a negative control (NC) were purchased from GenePharma (Suzhou, China). We used Lipofectamine 3000 (L3000015, Invitrogen, Waltham, USA) as the transfection reagent according to the manufacturer’s instructions. Forty-eight hours following transfection, the cells were harvested for further experiments, including RNA or protein extraction.

RNA extraction, reverse transcription as well as qRT-PCR

Total RNA was extracted utilizing the traditional assay, that is to say using TRIzol reagent (Invitrogen, USA). The cDNA synthesis with total RNA (1 µg) was carried out with the Primescript™ RT reagent kit (Vazyme, Nanjing, China), and quantitative real-time PCR (qRT-PCR) was conducted utilizing 2x SYBR^® Green Supermix (Vazyme, Nanjing, China). Relative gene expression was quantified via the 2-^ΔΔCt approach and normalized to GAPDH expression. The primer sequences used for amplification were as follows: CENPN forward, 5′-CTGTGTGAGGAAAAGCGTGC-3′; CENPN reverse, 5′-TCACCTGGTCCTTTACTCATCTG-3′; GAPDH forward, 5′-GGAGCGAGATCCCTCCAAAAT-3′, and GAPDH reverse, 5′-GGCTGTTGTCATACTT- CTCATGG-3′.

X-irradiation

Cells were exposed to X-rays at doses of 0, 2 and 10 Gy according to the experimental design. The equipment used was a Siemens Primus Accelerator (6 Mv; Siemens AG, Munich, Germany)t at Zhongnan Hospital of Wuhan University.

CCK-8 and colony formation assays

For the CCK-8 assay, cells were seeded into plates with 96 wells at 3000 cells per well. Ten microliters of CCK-8 reagent (Dojindo Laboratories, Kumamoto, Japan) were added to each well, and the plates were incubated at 37 °C for 1–4 h. The absorbance of the wells at 450 nm was determined in a microplate reader (Molecular Devices, USA).

For the colony formation assay, 1000 HepG2 cells per 1 ml in an overall volume of 2 ml was seeded into a plate with six wells. For Huh7 cells, 1,500 cells per 1 ml in a total volume of 2 ml was seeded per well. After 24 h, the cells were irradiated with X-rays at doses of 0 Gy and 2 Gy. After 14 days, the cells were treated with 4% paraformaldehyde for 15 min as well as stained with 0.1% crystal violet for 1 h. Finally, the number of colonies was calculated after washing the wells with water three times and drying.

Flow cytometry analysis

HCC cells were harvested 48 h following transfection and washed in PBS. The cells were stained with 1 ml of DNA staining solution and 10 µl of propidium iodide (PI) at room temperature for 30 min protected from light and then were analyzed by flow cytometry (Cat. #FC500, Beckman, USA).

Transwell assay

Cells transfected with siRNA after 24 h were seeded in a Transwell chamber system (Corning, USA). Approximately 2 × 10⁵ HepG2 cells or 8 × 10⁴ Huh7 cells in 200 µl of serum-free medium were seeded in the upper chamber, while 600 µl of medium containing 20% serum was added to the lower chamber. After a 48-hour incubation, cells in the upper chamber were wiped off with cotton swabs, and those in the lower chamber were fixed and stained as described in the colony formation experiment before they were photographed under an optical microscope.

Western blotting

Protein lysates were obtained from HepG2 and Huh7 cells with RIPA lysis buffer (#P0013B, Beyotime Biotechnology, China) containing 1% cocktail (#HY-K0021, MCE, USA). The protein concentration was detected with a BCA assay kit (Beyotime Biotechnology, Shanghai, China), and 30 µg of protein sample was divided by SDS-PAGE through a 10% gel and then transferred to PVDF membranes, which were blocked in 5% skim milk in TBST before they were blotted with the appropriate primary antibody followed by the corresponding secondary antibody. Detailed information on the antibodies used is listed in Table S1. An enhanced chemiluminescence kit was utilized to develop the bands.

Immunofluorescence staining

The immunofluorescence method was based on a previously described method (Wang et al., 2017). Briefly, cells were seeded into a 6-well plate or on cell slides and cultured overnight. After fixation for 15 min, the cells were incubated with 0.2% Triton X-100 for 15 min at RT and blocked with 5% BSA for 30 min. Next, the cells were incubated successively with primary antibody targeting CENPN, E-cadherin, N-cadherin, vimentin or γ-H2AX and appropriate secondary antibody before they were stained with DAPI at RT for 5 min. Finally, the cells were photographed under a fluorescence or confocal microscope.

Statistical analysis

All prognostic information from the TCGA-LIHC dataset was collected, and samples without results were excluded. On the basis of optimal sample separation, Kaplan–Meier (K-M) survival curves were generated to calculate survival, and the difference in survival between groups was determined by the log-rank test. The relationships between variables and patient survival were analyzed with Cox models, and the results were evaluated with SPSS 22.0 software (SPSS, Chicago, IL, United States). All cell experimental results were individually replicated 3 times. The values are reported as the means ± standard deviations and were statistically analyzed with GraphPad Prism 8.0 software (GraphPad Software, La Jolla, CA, United States). Student’s t-test was utilized for comparisons between groups. P < 0.05 was thought about statistically significant.

Preliminary screening and enrichment analysis of DEGs in HCC

Figure 1 illustrates the study design. The DEGs in HCC vs. normal tissues in the GSE87630 and GSE112790 datasets were analyzed by processing and standardizing the gene expression profile data through GEO2R (Figs. 2A, 2B, Figs. S1 and S2). GSE87630 contains 1,162 DEGs, of which 394 are upregulated and 768 are downregulated, whereas GSE112790 contains 1,713 DEGs comprising 963 upregulated DEGs and 750 downregulated DEGs. Then, R was used to intersect the two datasets, and a total of 532 overlapping DEGs were acquired: 171 upregulated genes and 361 downregulated genes (Fig. 2C, Table S2).

To clarify the functions of these 532 DEGs, GO as well as KEGG enrichment were carried out, and the ggplot2 package of R was used to visualize the top 10 per classification of GO and the top 35 of KEGG. In the GO BP category, the identified DEGs were related to the organic acid catabolic process and epoxygenase P450 pathway (Fig. 2D, Table S3). In the GO CC category, the DEGs were enriched in MCM complex and blood microparticle extracellular matrix (Fig. 2D, Table S3). Finally, in the GO MF classification, the DEGs were enriched in monooxygenase activity and iron ion binding (Fig. 2D, Table S3). KEGG pathway enrichment result showed that these DEGs play important roles in multiple key signaling pathways, including DNA replication, the cell cycle as well as the p53 signaling pathway (Fig. 2E, Table S4).

Construction of the PPI network as well as screening of hub genes

The STRING database was utilized to analyze the potential interacting proteins among 532 DEGs with a combined score >  0.7. For the PPI network, which contained 347 nodes and 1560 edges, was built with Cytoscape software (Fig. 3A). Then, the hub genes were analyzed with the cytoHubba plugin. Two topological analysis approaches, MCC and DMNC, were used to rank the top 30 nodes in the created PPI network (Figs. 3C and 3D, Table 1). Then, the intersection was calculated online using a Venn diagram (Fig. 3E, Table S5). Module analysis was performed with MCODE, and the top 5 modules were listed (Table S6). CENPN was screened as the seed gene in module 1 (the most significant module), which consisted of 32 nodes as well as 479 edges (Fig. 3B, Table S7). Thus, CENPN was considered the hub gene.

Study of CENPN expression in HCC cells and tissues via online databases

CCLE and EMBL-EBI were used to explore the RNA expression level of CENPN in HCC cell lines. The results revealed that CENPN was highly expressed in liver cancer cells (Figs. 4A and 4B). With respect to tissue expression levels, Oncomine evaluation revealed that CENPN expression was significantly higher in tumor tissues than in normal tissues (Fig. 4F), and CENPN expression in the GSE87630 and GSE112790 datasets was consistent with the above results (Figs. 4C and 4D); furthermore, this result was verified in the GSE25097 dataset (Fig. 4E). GEPIA revealed that the overall survival (OS) and disease-free survival (DFS) of liver cancer patients with high CENPN expression were shorter than those of patients with low CENPN expression (Figs. 4G and 4H). To better understand the accuracy of CENPN in HCC tumorigenesis, receiver operating characteristic (ROC) curves were drawn, and the AUC values of GSE87630, GSE112790 and GSE25097 were 0.913, 0.904 and 0.787, respectively (Figs. 4I–4K). Our results revealed that CENPN was highly expressed in liver cancer as well as affected the prognosis of patients.

Verification of the correlation between CENPN and clinicopathological features of HCC patients through the TCGA database

The optimal cutoff value was calculated by the survminer package of R, and the patients were separated into 2 groups: high CENPN expression and low CENPN expression. Compared with the low expression group, the high expression group had dramatically shorter survival (Fig. 5C). When the clinical features of CENPN were combined in the analysis, the Wilcoxon rank-sum test revealed that CENPN expression in tumor tissues was dramatically higher than that in normal (Fig. 5A). Similar results were obtained in the paired analysis of normal and tumor tissues from the same patient (Fig. 5B). These results suggest that CENPN expression is negatively associated with HCC prognosis. In particular, with the progression of grade, stage and T stage, CENPN expression showed an upward trend (Figs. 5D–5F). Furthermore, univariate Cox analysis revealed that CENPN, stage, and T and M stages affected HCC prognosis (Figs. 6A–6C, Table 2). Multivariate Cox analysis showed that high CENPN expression was an independent predictor of poor prognosis in HCC patients (Table 2). In summary, CENPN acts as an oncogene in liver cancer.

GSEA of the CENPN gene

GSEA showed that high CENPN expression was mainly associated with liver cancer subclass proliferation, the cell cycle, the p53 signaling pathway, the Rb1 pathway, cell cycle checkpoints, positive regulation of the G1/S phase transition, reactome SCF SKP2-mediated degradation of p27 p21, E2F targets and more correlated gene sets (P < 0.05; Figs. 7A–7H). High CENPN expression was also associated with DNA damage, DNA damage detection, and DNA damage response signal transduction by p53 class mediators and to other DNA damage and repair-related gene sets (P < 0.05; Fig. 8A).

CENPN downregulation arrests the cell cycle in G1 phase and inhibits HCC cell proliferation but has no effect on migration and invasion

The above results suggest that CENPN functions as a biomarker for the diagnosis of HCC. To further explore the mechanism by which CENPN impacts the biological behavior of HCC, we established a cell model of CENPN deficiency (HepG2 and Huh7) by transfecting siRNA targeting CENPN and NC siRNA. qRT-PCR, western blot and immunofluorescence analyses were performed 48 h after transfection to confirm knockdown (Figs. 7I, 7J, 7M and 7O). CCK-8 assays revealed that CENPN deficiency considerably inhibited the proliferation and viability of HCC cells (Figs. 7K, 7L), and immunofluorescence analysis showed that when CENPN was knocked down, Ki67 expression decreased significantly (Figs. 7N, 7P). Cell cycle analysis showed that CENPN deficiency significantly inhibited the transition from G1 to S phase in HCC cells (Figs. 7Q–7X). However, Transwell assays revealed that there was no difference in the numbers of cells that migrated or invaded after CENPN knockdown compared with those subjected to NC siRNA (Figs. S3A, S3B). Immunofluorescence results indicated that there was no significant difference in E-cadherin, N-cadherin and vimentin expression between the groups (Figs. S3C–S3G), and the western blot results were consistent with the above data (Figs. S3H, SI). In summary, these findings indicate that CENPN knockdown inhibits HCC cell proliferation but does not affect invasion or migration.

CENPN affects the p21-CDK2/cyclin E, p27-CDK4/cyclin D and Rb/E2F1 signaling pathways in HCC

To elucidate the mechanism by which CENPN promotes HCC, we used western blotting to analyze the expression levels of related markers. The western blot results revealed that compared with the NC, the CENPN-knockdown group showed enhanced expression of P53, P27 and P21 and also lowered expression of CDK4, cyclin D1, CDK2, cyclin E as well as pRb, but there was no significant change in total Rb levels. In addition, the expression of E2F1 and c-Myc decreased (Fig. 8F). Taken together, these results indicate that CENPN downregulation prevents cell proliferation via regulating p21 and p27.

Decreased CENPN expression promotes radiation damage in HCC cells

GSEA indicated that CENPN is closely related to DNA damage and repair functions. Cells transfected with siRNA-CENPN and NC for 48 h were treated with X-rays at doses of 0 Gy and 10 Gy, and γ-H2AX foci were detected by immunofluorescence 24 h later. The number of γ-H2AX foci increased after CENPN knockdown and was further pronounced upon 10 Gy X-ray exposure (Figs. 8B, 8C). Moreover, colony formation assays suggested that the number of HCC cell colonies formed after CENPN knockdown combined with X-ray treatment was significantly lower than that after siRNA-CENPN treatment alone (Figs. 8D, 7E). These results suggest that interfering with CENPN expression enhances X-ray-induced radiation damage in HCC cells.