Identification of Flap endonuclease 1 as a potential core gene in hepatocellular carcinoma by integrated bioinformatics analysis

Data collection and DEG validation

Four liver cancer-related datasets (GSE14520, GSE29721, GSE45267 and GSE60502), were downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/). These datasets contained a total of 299 tumor samples and 289 non-tumor samples. GSE14520 contained 225 liver cancer samples and 220 adjacent controls; GSE29721 contained 10 pairs of liver cancer samples and adjacent control tissue; GSE45267 contained 46 liver cancer samples and 41 adjacent controls; and GSE60502 contained 18 pairs of liver cancer tissues and adjacent controls. The inclusive criteria for datasets were as follows: (1) the samples were from human HCC tissues and paired adjacent or non-tumor tissues; (2) gene expression profiling of mRNA; (3) each dataset included no less than ten paired samples; (4) the datasets were from 2009 to 2019. The exclusive criteria for datasets were as follows: (1) the HCC samples were induced by specific diseases, including chronic hepatitis B, hepatitis C, hepatitis D, alcoholic liver disease, nonalcoholic fatty liver disease and cholangiocarcinoma; (2) details of the samples were not available; (3) the datasets could not be analyzed by GEO2R; (4) no gene symbol. The specific platform information for these four datasets is compiled in Table 1. GEO2R was utilized to identify DEGs in these studies, using the screening criteria: |logFC(foldchange)| ≥ 1, P < 0.05, and adjusted P < 0.05. In addition, genes with multiple probe set or probe sets lacking matched gene symbols were removed or averaged, respectively. Then, the overall DEGs, as well as those that were up- and down-regulated in the four datasets were intersected and visualized using Funrich (v 3.0, http://funrich.org/index.html).

Functional enrichment analyses

Gene ontology analyses focuses on three domains: biological processes (BP), cellular components (CC), and molecular functions (MF), and such analyses are commonly used to understand the biological functions, pathways, or localization of DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis database surveys as a valuable resource for assessing how particular DEGs may be involved in or influenced by specific signaling pathways and disease states. The DAVID (https://david.ncifcrf.gov/) was utilized for functional enrichment analyses, with a P < 0.05 cutoff for significance.

PPI network analysis

The STRING database (v10.5; https://string-db.org/) was used to generate a DEG PPI, with a minimum interaction score cutoff of 0.4. Cytoscape (v 3.4.0, https://cytoscape.org/) was used for network visualization, and the CytoHubba plug-in was used to identify hub genes with the criteria of filtering degree ≥10. The MCODE plug-in was used to construct key clustering modules (MCODE score >10, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100 and k-score = 2).

Validation of core genes

Hub genes among overall DEGs and the most critical clustering module were identified through the CytoHubba plug-in, and intersecting core genes were identified. The UCSC Cancer Genomics Browser (https://genome-cancer.ucsc.edu/) was then used for hierarchical clustering of these core genes. In addition, the expression profiles of these core genes in 421 TCGA liver cancer tissues, including 50 solid normal tissues and 371 primary tumors, were determined by analyzing available datasets. A core gene PPI network was constructed with the cBioportal online database (http://www.cbioportal.org/). Furthermore, correlations between the expression of FEN1 and MCM2, RFC4, or BIRC5 in TCGA liver cancer tissues were investigated. Finally, we used the ONCOMINE (https://www.oncomine.org/resource/login.html) and FIREBROWSE online database (http://firebrowse.org/) to investigate the FEN1 expression in various cancers including HCC.

Survival analysis

The GEPIA database (http://gepia.cancer-pku.cn/) was used to conduct survival analyses based on core gene expression, with hazard ratios (HRs) and 95% confidence intervals being calculated, and logrank P value <0.05 being the threshold of statistical significance.

Clinical samples

A total of 34 paired HCC tumor tissues as well as corresponding adjacent non-cancerous tissues were obtained from our Hospital’s Department of Hepatobiliary Surgery, with all patients providing informed consent. The content of the informed consent includes research purposes, risks and discomfort, benefits, and privacy issues. The form of the informed consent is in written. The study was examined and approved by the Second affiliated Hospital of Chongqing Medical University Ethics Committee (approval number: 2018-60).

Cell culture

Seven human liver cancer cell lines (SMMC-7721, BEL-7404, HCCLM3, HepG2, MHCC97-H, SK-HEP-1, and Huh-7) as well as normal human liver HL7702 cells were donated by the Institute for Viral Hepatitis, Chongqing Medical University. All cells were grown in high glucose DMEM (Gibco, Waltham, MA, USA) containing 10% FBS (Corning, Corning, NY, USA) at 37 °C in a 5% CO₂ incubator.

Hematoxylin and eosin (H & E) staining

Paraffin-embedded tissues were dewaxed in xylene I, II, and III for 20 min each, and then dehydrated in an ethanol gradient (100%, 95%, 90%, 80%, and 70%), 3 min per step. Sections were then rinsed with distilled water for 5 min, and the nucleus was counterstained with hematoxylin for 3 min. Sections were washed again in water, followed by differentiation for 30 s in a 75% hydrochloric acid alcohol solution, and blue color was returned by washing with distilled water for 5 min. A red dye was then used for counterstaining for 5 min, after which samples were dehydrated in 70%, 80%, 95%, and 100% ethanol, 1 min per concentration. Sections were then cleared with xylene and sealed using neutral gum.

IHC analysis

IHC staining was conducted as previously described (Zhang et al., 2016). Briefly, paraffin-embedded sections were incubated at 56 °C for 2 h, and then 3% hydrogen peroxide was used for antigen retrieval. Afterwards, the sections were incubated with rabbit anti-human FEN1 (1:100, A1175, ABclonal, China) at four °C overnight. Then, sections were probed for 1 h using HRP-conjugated secondary antibodies at 37 °C, after which a DAB substrate kit was utilized, and hematoxylin was used for nuclear staining. The Image-Pro Plus (IPP) software (Media Cybernetics, Rockville, MD, USA) was used to quantify staining intensity. The frequency of positively-stained cells was determined on a 0–100 scale, while staining intensity was scored as follows: 0 = negative; 1 = weak; 2 = moderate; 3 = strong. These two scores were then multiplied together to yield an IHC score between 0 and 300. The final scores were assigned by two independent pathologists. The mean IHC score was used as a cutoff value to separate patients into low- and high-expression groups.

RT-qPCR analysis

TRIzol (ThermoFisher Scientific, Waltham, MA, USA) was used for total RNA extraction, and RNA was then reverse transcribed with the PrimeScript RT-PCR kit (Takara Bio, Dalian, China) based on provided protocols. The SYBR Premix Ex Taq II (Takara, Japan) kit was used to conduct RT-qPCR analysis on a Bio-Rad CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA). The reaction conditions were as follows: pre-denaturation at 95 °C for 30 s, followed by 40 cycles of amplification at 95 °C for 5 s and 60 °C for 30 s, and a final cycle along the melting curve from 65 °C to 95 °C in increments of 0.5 °C for 5 s. mRNA levels were determined via the 2^−ΔΔCt method, with GAPDH used for normalization. Primers used were: GAPDH F: 5′-GGTGGTCTCCTCTGACTTCAACA-3′ and R: 5′-GTTGCTGTAGCCAAATTCGTTGT-3′, FEN1 F: 5′-CTGTGGACCTCATCCAGAAGCA-3′ and R: 5′-CCAGCACCTCAGGTTCCAAGA-3′.

Statistical analysis

Statistical analyses and graphing were performed with SPSS v19.0 (SPSS Inc., Chicago, IL, USA) and Graph Pad Prism v8.0 (Graph Pad Software, San Diego, CA, USA), respectively. Data are means ± standard deviation (SD). Student’s t-tests were used to compare groups. Fisher’s exact test was used to assess correlations between the expression of FEN1 and HCC patient clinicopathological features. Spearman’s correlation analyses were used to compare the expression of pairs of genes in TCGA liver cancer tissues. P < 0.05 was the significance threshold (*P < 0.05, **P < 0.01).

HCC-associated DEG identification

In this study, 1,088 total DEGs (505 and 583 up- and down-regulated, respectively) in GSE14520, 1,449 total DEGs (837 and 612 up- and down-regulated, respectively) in GSE29721, 1,604 total DEGs (713 and 891 up- and down-regulated, respectively) in GSE45267, and 1,533 total DEGs (792 and 741 up- and down-regulated, respectively) in GSE60502 were screened. Based on these datasets, a total of 409 overlapping DEGs were identified among these four datasets (142 and 267 up- and down-regulated, respectively), as visualized with the Funrich software (Fig. 1). DEGs are listed in Table S1.

DEG functional enrichment analyses

To explore the biological activities of these DEGs, the DAVID database was used to conduct GO and KEGG enrichment analysis. With respect to BPs, up-regulated DEGs were primarily enriched in processes such as mitotic nuclear division, cell division, cell cycle, DNA replication, and mitotic sister chromatid segregation (Fig. S1A), while down-regulated DEGs were primarily enriched in processes such as redox process, the cytochrome 450 pathway, drug metabolism, and negative regulation of growth (Fig. S2A). With respect to CCs, up-regulated DEGs were primarily enriched in the nucleoplasm, nucleus, cytoplasm, spindle, and cellular intermediates (Fig. S1B), while down-regulated DEGs were mostly enriched in extracellular exosomes, organelle membranes, blood microparticles, extracellular regions, and the mitochondrial matrix (Fig. S2B). With respect to MFs, up-regulated DEGs were primarily enriched in functions such as protein binding, ATP binding, DNA helicase activity, protein kinase binding, single-stranded DNA binding, and chromatin binding (Fig. S1C), whereas down-regulated DEGs were primarily associated with iron ion binding, oxidoreductase activity, heme binding, monooxygenase activity, and oxygen binding (Fig. S2C).

A KEGG analysis revealed that up-regulated DEGs were particularly enriched in pathways such as the cell cycle, DNA replication, P53 signaling, and tumor pathways (Fig. S1D). However, down-regulated DEGs were mostly associated with metabolic pathways, fatty acid degradation, chemical carcinogenesis, and PPAR signaling (Fig. S2D).

PPI network and module analyses

To better understand interactions among DEGs, the STRING online database was used to generate a PPI network consisting of 403 nodes and 3,502 edges, which was visualized using Cytoscape. Six of the 409 DEGs were not included in this network (Fig. 2A). This network was then analyzed using the MCODE plug-in, and three clustering modules were filtered out according to the chosen screening conditions. Clustering module 1 scored 58.492 with 62 nodes and 1,784 edges (Fig. 2B), clustering module 2 scored 11.529 with 18 nodes and 98 edges (Fig. 2C), and clustering module 3 scored 10.917 with 25 nodes and 131 edges (Fig. 2D). The genes in clustering module 1 were up-regulated DEGs, whereas those in the other two modules were primarily down-regulated DEGs.

Functional enrichment analysis of key clustering modules

The DAVID database was next used to explore the biological functions of genes in these key clustering modules (Tables S2–S4). With respect to BPs, clustering module 1 was primarily enriched in cell differentiation, mitotic nuclear division, DNA replication, and DNA helicase activity, while clustering module 2 was primarily enriched in plasminogen activation, coagulation, cytolysis, and complement activation regulation, and clustering module 3 was primarily enriched in steroid metabolism, heterogeneous biomass metabolism, and exogenous drug catabolism. With respect to CCs, clustering module 1 was primarily enriched for the nucleoplasm, nucleus, intermediate, spindle, cytoplasm, nuclear chromosome, while clustering module 2 was primarily enriched for exosomes, extracellular regions, membrane attack complexes, and extracellular vesicles, and clustering module 3 was primarily enriched for organelle membranes, the endoplasmic reticulum membrane, and high-density lipoproteins. With respect to MFs, clustering module 1 was primarily enriched for protein binding, ATP binding, DNA helicase activity, protein kinase binding, and DNA binding, while clustering module 2 was primarily enriched for endopeptidase activity, transcription factor binding, steroid binding, RNA polymerase II transcription factor activity, and enzyme binding activity, and clustering module 3 was primarily enriched for aerobic binding, iron ion binding and heme binding.

A KEGG analysis revealed that clustering module 1 was primarily enriched in the cell cycle, oocyte meiosis, DNA replication, and p53 signaling, while clustering module 2 was primarily enriched in the complement system, prion disease, and systemic lupus erythematosus, and clustering module 3 was primarily enriched in chemical carcinogenesis, retinol metabolism, P450 drug metabolism, and metabolic pathways.

Identification of core genes and analysis of their clinical significance

Next, core genes involved in HCC were identified based on their levels of interaction via analyzing our PPI network using the Cytoscape program. Based on our clustering module analysis, clustering module 1 included 62 genes was found closely related to the progression of HCC. Then, nine total genes (BIRC5, DLGAP5, DTL, FEN1, KIAA0101, KIF4A, MCM2, MKI67, and RFC4) were identified based on the intersecting genes among the top 40 genes derived from 12 different algorithms by the cytoHubba plug-in (Table S5). In addition, 185 hub genes were identified across 403 nodes in our PPI network based on the filtering degree ≥10 criteria (Table S6). The nine identified core genes were found belonged to this larger subset of hub genes (Fig. 3A).

To investigate core gene expression in HCC, a hierarchical clustering analysis was performed using the UCSC Cancer Genomics Browser, revealing that these nine core genes were highly expressed in most liver cancer samples (Fig. 3C). Next, the expression profiles of these nine core genes in 421 TCGA liver cancer tissues, including 50 solid normal tissues and 371 primary tumors, were downloaded and analyzed, revealing that the expression of these core genes was significantly elevated in HCC (Fig. 3B). The correlations between core gene expression levels and patient prognosis in 182 total HCC samples were further assessed, revealing that BIRC5 expression (HR = 2, logrank P = 6.7e−05) was correlated with worse overall survival (OS) for HCC patients, as was that of DLGAP5 (HR = 1.9, logrank P = 0.00039), DTL (HR = 1.7, logrank P = 0.0049), FEN1 (HR = 1.5, logrank P = 0.022), KIAA0101 (HR = 1.7, logrank P = 0.002), KIF4A (HR = 1.8, logrank P = 0.001), MCM2 (HR = 1.7, logrank P = 0.0022), MKI67 (HR = 1.9, logrank P = 0.00045), and RFC4 (HR = 1.7, logrank P = 0.004) (Fig. 4). Furthermore, the expression of BIRC5 (HR = 1.6, logrank P = 0.002) was correlated with decreased disease-free survival (DFS) for HCC patients, as was that of DLGAP5 (HR = 1.6, logrank P = 0.0033), DTL (HR = 1.6, logrank P = 0.0016), FEN1 (HR = 1.5, logrank P = 0.0075), KIAA0101 (HR = 1.6, logrank P = 0.0022), KIF4A (HR = 1.6, logrank P = 0.0011), MCM2 (HR = 1.6, logrank P = 0.0034), MKI67 (HR = 1.9, logrank P = 4.2e−05), and RFC4 (HR = 1.5, logrank P = 0.011) (Fig. 5). High expression of these nine core genes was associated with significantly reduced survival among HCC patients.

FEN1 may be a key candidate gene in HCC

To clarify the PPI network for these nine core genes, the cBioportal online database were explored and identified an interaction between FEN1, MCM2, BIRC5 and RFC4 (Fig. 6A). Previous studies have confirmed that MCM2, BIRC5, and RFC4 are abnormally highly expressed in HCC, and that they participate in the regulation of HCC tumor biology. However, the expression of FEN1 and its clinical significance in HCC is unclear. Therefore, the correlations between the expression of FEN1 and these three other genes were further analyzed. The expression of FEN1 in TCGA liver cancer tissues was positively correlated with that of MCM2 (r = 0.853, P = 0.000), BIRC5 (r = 0.809, P = 0.000), and RFC4 (r = 0.852, P = 0.000) (Figs. 6B–6D), suggesting that FEN1 may play as important a role in the progression of liver cancer as do MCM2, BIRC5, and RFC4. The ONCOMINE and FIREBROWSE databases were further utilized to investigate the expression of FEN1 in various cancers, including HCC. The results revealed that FEN1 expression was clearly elevated in most cancers, including bladder, breast, colorectal, esophageal, lung, and liver cancer (Figs. 7A and 7B). In addition, FEN1 was found overexpressed in three HCC-related datasets (Fig. 7C).

Experimental validation

In this study, 34 paired HCC and adjacent control tissues were used to verify the expression of FEN1 in HCC. The pathologic diagnosis of HCC and matched adjacent samples were confirmed by H&E staining (Figs. 8A and 8C). FEN1 staining was localized to the nucleus and cytoplasm, with stronger expression in the HCC samples relative to the adjacent controls (Figs. 8B and 8D). Furthermore, IHC analysis indicated that FEN1 levels were significantly higher in HCC relative to adjacent tissues (Fig. 8E), with IHC scores being significantly higher in HCC samples (Fig. 8F).

The correlation between FEN1 expression and HCC patient clinicopathological features was investigated via Fisher’s exact test. As shown in Table 2, there were significant correlations between FEN1 expression and tumor size (P = 0.047 < 0.05), metastasis (P = 0.013 < 0.05) and vascular invasion (P = 0.024 < 0.05). FEN1 expression did not significantly correlated with gender, age, tumor multiplicity, TNM stage, pathological grade, HBsAg, liver cirrhosis or serum alpha-fetoprotein (AFP) (P > 0.05). Furthermore, FEN1 mRNA levels were found significantly elevated in six human hepatoma cell lines relative to that in the normal human liver cell line, HL7702 (Fig. 8G).