The gene expression and clinical data covering age, gender, stage, and survival data were rooted in the GEO (https://www.ncbi.nlm.nih.gov/geo/), ICGC (https://icgc.org/), and TCGA databases (https://tcga-data.nci.nih.gov/). Gene expression and clinical data of the GSE14520 series (Roessler et al., 2010) and the GSE76427 series (Grinchuk et al., 2018) originated from the GEO database. Gene expression and clinical data of GSE14520 based on the GPL3921 platform [HT_HG-U133A] (Affymetrix HT Human Genome U133A Array), including 241 non-tumor specimens and 247 tumor specimens (including 227 paired HCC specimens). Gene expression and clinical data of GSE76427 based on the GPL10558 platform (Illumina HumanHT-12 V4.0 expression bead chip), including 52 non-tumor specimens and 115 tumor specimens (including 52 paired HCC specimens). The gene expression files (ICGC-LIRI-JP) of the ICGC database based on the Illumina HiSeq RNA Seq platform, including 202 non-tumor specimens and 243 tumor specimens (including 202 paired HCC specimens), were obtained from the ICGC database. The RNA sequencing and clinical data (TCGA-LIHC-FPKM) based on the Illumina HiSeq RNA-Seq platform, including 50 non-tumor specimens and 374 tumor specimens (including 50 paired HCC specimens), were obtained from the TCGA database.

The DEGs of non-tumor specimens and tumor specimens were identified with R studio by dealing with the raw data of GSE14520, GSE76427, ICGC, and TCGA databases. In the four datasets, the RMA package was utilized to standardize the entire raw data, the Affy package was utilized to evaluate the quality of data, and the Limma package was utilized to filtrate DEGs whose standard was deemed as the | logFC| > 2 as well as the adjusted p < 0.05. In addition, the survival package of R studio was also utilized to filtrate survival-related genes from GSE76427 and ICGC, in which the cutoff criterion was p-values < 0.05. The common genes between DEGs and survival-related genes were selected by R studio, painting a Venn diagram.

The STRING database (http://string-db.org)(Szklarczyk et al., 2015) formed the PPI network where the MCODE plug-in chose the essential modules in the Cytoscape application. In addition, the standards of choosing the essential modules are k-core = 2, degree cutoff = 2, node score cutoff = 0.2, and maximum depth = 100. CytoHubba plug-in was another plug-in in Cytoscape software in which the top 10 hub genes were selected in the PPI network. What is more, the 12 calculating methods in CytoHubba presented as the following: DMNC, ClusteringCoefficient, EPC, Degree, BottleNeck, EcCentricity, MNC, Radiality, Betweenness, MCC, Stress, and Closeness, which generated 12 different outcomes; through comparing with them, the NAP1L1 was chosen as the objective of our study.

The Oncomine database (https://www.oncomine.org/resource/login.html)(Rhodes et al., 2007) consists of 86,733 specimens and 715 gene expression data. Moreover, it is a unitary integrated database aimed to accelerate data mining efforts effectively. Therefore, this database evaluated the NAP1L1 expression in human tumors. Accordingly, HPA (https://www.proteinatlas.org) is deemed a kind of online tool whose goal was to elucidate the level of NAP1L1 protein exiting in non-tumor tissues and HCC tissues (Pontén et al., 2011). In 2003, the Swedish-based program, the HPA, witnessed the world project of human proteins. Specifically, the cells, tissues, and organs applying the combination of diverse omics technologies containing antibody-based imaging were all taken in human proteins. HPA version 21.0, appearing in November 2021, was used in this study (Uhlén et al., 2015; Thul et al., 2017; Uhlen et al., 2017). The immunohistochemistry pictures originated from HPA. The expression of the NAP1L1 protein was explained in non-tumor tissues and HCC tissues among HCC patients.

R studio was utilized to achieve survival analysis, clinical feature analysis, and univariate/multivariate regression analysis. Patients’ gene expression datasets and clinical data were used in this study, including GSE14520, GSE76427, ICGC, and TCGA datasets. Nine packages are being taken in this research: the beeswarm, Limma, survival, survminer, Ggpubr, plyr, grid, ggplot2, and gridExtra packages. To be more specific, the beeswarm and Limma packages were used to conduct the differential expression analysis, the survival package and survminer package were applied to perform the survival analysis, the survival package also was utilized to do the univariate/multivariate regression analysis, the Ggpubr package was utilized to carry out the analysis of the clinical features, and the plyr, ggplot2, grid, and gridExtra packages were performed to analyze the multiple GSEA.

GSEA4.2.1 application was derived from the Broad Institute (http://www.gsea-msigdb.org/gsea/login.jsp) (Subramanian et al., 2005), which still played a significant role in generating enrichment results. Plus, R studio was visualized well to complete multiple GSEA results.

CIBERSORT was utilized to figure out the enrichment files of immune cells in the HCC sample ground on GSE14520, GSE76427, ICGC, and TCGA (Chen B. et al., 2018). TISIDB (http://cis.hku.hk/TISIDB/index.php) was utilized to detect the interactions between 28 types of TILs and 30 types of cancers in humans to check the tumor–immune system interactions (Ru et al., 2019). The association between NAP1L1 and TILs was tested through Spearman’s test, in which all hypothetical tests were deemed as statistically significant in that they were bilateral and p-value < 0.05. Then, to verify the enrichment of six immune cells (dendritic cells, B cells, CD4⁺T cells, CD8⁺T cells, macrophages, and neutrophils) in HCC patients, the TIMER online database was used to verify the relations between NAP1L1 expression and the level of immune infiltration (Li T. et al., 2017).

IPS is a terminology that comprises MHC molecules, immunomodulators, effector cells, and suppressor cells. According to the expression of representative genes, IPS could be figured out when applying it, and the TCIA (https://tcia.at/home) was applied to acquire the IPS about HCC patients. R studio was utilized to achieve immune cell infiltration analysis. The immune cell infiltration analysis needed Limma, ggExtra, vioplot, ggpubr, ggplot2, and Venn Diagram packages. Specifically, the Limma package was used to acquire the results of CIBERSORT, and the vioplot package was utilized to paint the violin plot. Moreover, ggpubr, ggplot2, ggpubr, and ggExtra were applied to perform the correlation analysis between NAP1L1’s expression and the immune cell, and to paint the correlation graph; the Venn Diagram package was used to paint the Venn graphs.

As the objective of our study, the normal liver cell lines (MIHA) and the HCC cell lines (MHCC-LM3, PRF-5, Huh-7, MHCC-97H, and Hep-G2) are originated from the cell bank of the Chinese Academy of Sciences (Shanghai, China). Furthermore, they were maintained in a DMEM cell culture (Gibco, New York, United States) supplemented with 10% FBS (Gemini, West Sacramento, United States) in the incubator at 37°C with 5% CO2.

NAP1L1 siRNA was synthesized from RiboBio (Guangzhou, China) and was used to reduce the NAP1L1 expression. MHCC-97H and Huh-7 cells were cultivated into each well of the six-well plates until they were 30–60% confluent. After that, cells were washed, placed in a serum-free medium, and transfected with siRNA using Lipofectamine™2000 according to the manufacturer’s instructions (Invitrogen, MA, United States). After 6 h, the medium was changed to a complete medium, and cells were cultured at 37°C in 5% CO₂. Four groups were generated for all experiments, a negative control group (negative control); a siRNA1 group (siRNA1); a siRNA2 group (siRNA2); a siRNA3 group (siRNA3).

The whole RNA was obtained at first by utilizing the TRIzol reagent and then going through the reverse transcription way to reach the cDNA under the guidance of the Thermo Fisher Scientific Reagent Kit. Finally, we did a qRT-PCR reaction of cDNA, taking a Bio-Rad qRT-PCR Checking System.

First of all, we centrifuged and gathered the cells when cells were transfected after 48 h. Then, the radioimmunoprecipitation assay buffer was applied to acquire the whole-cell protein lysates. Furthermore, detaching by SDS-PAGE, the protein lysates were blotted into PVF membranes (Bio-Rad). The principal antibodies, including NAP1L1 (14898-1-AP; 45kDa; 1:1000), CDK1 (19532-1-AP; 34kDa; 1:1000), β-catenin (19532-1-AP; 45kDa; 1:1000), and GAPDH (10494-1-AP; 36kDa; 1:1000), were attached to the membranes after blocking, which they stored one night at 4°C. Next, after the TBST was washed three times, the second antibody was hatched for 2 h at homeothermy, and then TBST was eluted three times once again. In the end, the Western blot was visualized by imaging systems.

The day before the transfection, MHCC-97H and Huh-7 cells were planted into the culture plates. When they were transfected successfully, those cells were eluted with PBS and fitted in 75% pre-cooling alcohol. Before applying the flow cytometer, MHCC-97H and Huh-7 cells were hatched in 500 μl specimen buffer and 0.25 mg/ml of RNase A, maintaining thirty minutes at homeothermy. Subsequently, BD FACSCalibur Flow Cytometer was utilized to count the quantities of the 2 cells in each stage, and FlowJo was applied to explain these data.

In terms of the statistical analysis, we considered the two software: GraphPad Prism (version: 8.2.1) and R studio (version:4.0.5). Data analysis was regarded as meaningful data, while the formula p < 0.05 showed the meaning of a statistically significant difference.

R studio demonstrated the DEGs in every database (3,641 in GSE14520, 13,299 in GSE76427, 5,172 in ICGC, and 3,642 in TCGA) and DEGs whose standard was deemed as the | logFC| > 2 as well as the adjusted p < 0.05. At the same time, R studio was also used to paint heatmaps and volcano plots about DEGs in every database, presented in figures (Figures 2A–H). The 2,145 overlapping DEGs were examined in the aforementioned databases, and the results were displayed in a Venn graph (Figure 3A). In addition, we continued to do a survival filter using R studio (1,710 in GSE76427 and 2,691 in ICGC), and the cutoff criterion was p-values < 0.05. Therefore, 101 survival-related genes were screened out. The outcomes were displayed in a Venn graph (Figure 3B). Furthermore, we intersected 33 common genes from DEGs and survival-related genes whose outcomes were displayed in a Venn graph (Figure 3C).

The STRING was applied to produce a PPI network including 33 common genes taken by the Cytoscape application. It covered 26 high-expression genes labeled with red color and seven low-expression genes labeled with green color. Here, we only showed 25 genes related to each other (Figure 3D). CytoHubba plug-in functions as a plug-in in Cytoscape application were selected the top 10 hub genes in the PPI network. The 12 calculating methods in cytoHubba presented as the following: DMNC, ClusteringCoefficient, EPC, Degree, BottleNeck, EcCentricity, MNC, Radiality, Betweenness, MCC, Stress, and Closeness, which generated 12 various results in the top 10 key genes through comparing with them. Moreover, we claimed that NAP1L1 obtained the highest score in DMNC (Figure 3E) and ClusteringCoefficient (Figure 3F) so that the NAP1L1 was chosen as the objective of our study.

Gene expression analyses using the ONCOMINE and the TIMER claimed that NAP1L1 mRNA levels were evidently higher in HCC than in non-tumor tissues (Figures 3G,H). The HPA database presented the expression levels of the NAP1L1 protein in HCC. The levels of the NAP1L1 protein were low in normal liver tissues, and the high expression levels of NAP1L1 were observed in HCC tissues (Figure 3I). Subsequently, the qRT-PCR measured the NAP1L1 mRNA level (Figure 3J). Thus, the expression of NAP1L1 was higher based on the normal liver tissue in HCC tissues.

First, the gene expression profile and the clinical information (the GSE14520, GSE76427 datasets, ICGC datasets, and TCGA datasets) were obtained from the GEO, ICGC, and TCGA databases. From those scatterplots (Figure 4A, Supplementary Figures S1A,2A, 3A), we concluded that the NAP1L1 expression of HCC specimens was evidently higher than that of the non-tumor specimens (p < 0.05). From those survival curves (Figure 4B, Supplementary Figures S1B,2B, 3B), we found that the high expression of NAP1L1 is closely associated with a bad OS (p < 0.05).

Furthermore, we also sought the correlations between the expression of NAP1L1 and the HCC patients’ clinical features from age (Figure 4C, Supplementary Figures S1C,S2C, S3C), gender (Figure 4D, Supplementary Figures S1D,S2D, S3D), and stage (Figure 4E, Supplementary Figures S1E,S2E, S3E). Those results showed a significant difference between the expression of NAP1L1 and the stage of HCC patients (p < 0.05). The Conclusion also showed that the significant difference did not exist in age (p > 0.05) and gender (p > 0.05). Then, the univariate (Figure 4F, Supplementary Figures S1F,S2F, S3F) and the multivariate regression analyses (Figure 4G, Supplementary Figures S1G,S2G, S3G) were performed with the clinical data of HCC patients, which reveals that NAP1L1 was statistically significant in the univariate regression analysis (p < 0.05) and the multivariate regression analysis (p < 0.05). This result claimed that NAP1L1 was an independent prognostic factor in HCC patients.

Finally, to deeply investigate how the NAP1L1 promotes HCC progression, the HCC patients were separated into the high and low expression groups ground on the expression of NAP1L1. Furthermore, we took multiple GSEA employing the GO and the KEGG, and the enrichment results of the top five genes of upregulation and the top five genes of downregulation were displayed in GO and KEGG gene sets (Figures 4H,I, Supplementary Figures S1H, S1I, S2H, S2I, S3H, S3I).

In Conclusion, the aforementioned results in GSE14520 datasets, GSE76427 datasets, ICGC datasets, and TCGA datasets showed that NAP1L1 is a high expression in HCC tissue and is closely associated with poorer prognosis of HCC patients. In addition, we found correlations between the expression of NAP1L1 and clinical features in HCC patients, and NAP1L1 is an independent prognostic factor in HCC patients. which could promote the progression of HCC by some complex mechanisms.

We closely understood the aforementioned GSEA results to explore complex mechanisms of NAP1L1 that promoted HCC progression. Hence, we further found that cell cycle and Wnt signaling pathway can be found in the GSE14520 dataset, GSE76427 dataset, ICGC dataset, and TCGA dataset, based on GSEA. The result demonstrated that NAP1L1 might accelerate the progression of HCC, which relies on the cell cycle (Figures 5A–D), especially the G2/M phase transition (Figures 5E–H) and Wnt signaling pathway (Figure 5I–L).

The aforementioned results show that the G2/M phase transition of the cell cycle and Wnt signaling pathway can be found in the GSE14520 dataset, GSE76427 dataset, ICGC dataset, and TCGA dataset through GSEA. Then, the NAP1L1 expression level was measured in NC, si-1, si-2, and si-3 groups in MHCC-97H and Huh7 cell lines (Figure 6A) via the Western blot. The outcomes demonstrated no changed features in NC and si-3 groups, but si-1 and si-2 groups were evidently low in terms of NAP1L1 expression levels in them.

β-catenin was regarded as a protein to regulate the Wnt signaling pathway. CDK1 was regarded as a protein to regulate G2 to M phase transition of the cell cycle. The expression of β-catenin and CDK1 was measured using the Western blot (Figure 6B, Supplementary Figure S4A). The outcomes showed that after NAP1L1 silencing, β-catenin and CDK1 both downregulate accordingly, which meant that the expression of NAP1L1 had a positive relationship between β-catenin and CDK1. In addition, after NAP1L1 silencing, cell numbers in the G2 phase (Supplementary Figures S4B, C) were significantly increased in both MHCC-97H and Huh7 cells, measured by the flow cytometry assay.

Therefore, we believed that NAP1L1 knockdown could inhibit HCC cell proliferation by inhibiting Wnt/β-catenin pathway activation and the G2/M phase of cell cycle transition.

For exploring the association between NAP1L1 and the immune cells in HCC patients, we conducted the following analysis:

The CIBERSORT algorithm was taken to acquire the landscape of tumor-infiltrating immune cells. Bar plots were painted through R studio to test the correlation between NAP1L1 expression and 22 immune cells. From those bar plots (Figure 7A, Supplementary Figures S5A, S6A, S7A), we can observe the percent of 22 immune cells in each HCC sample. The violin plot was utilized to recognize the variety of 22 immune cells between the high and the low expression groups in HCC samples ground on the median expression level of NAP1L1. Furthermore, from observing those violin plots (Figure 7B, Supplementary Figures S5B, S6B, S7B), we found a remarkable difference in activated memory CD4 T cells, activated dendritic cells, gamma delta T cells, regulatory T cells (Tregs), CD8 T cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, and M2 macrophages.

The correlation graph was utilized to recognize the correlation between NAP1L1 and 22 immune cells. Then these significant immune cells were selected and presented in the following figures. They activated Mast cells (Figure 7C), M0 macrophages (Figure 7D, Supplementary Figures S5C, S7D), M2 macrophages (Figure 7E, Supplementary Figures S6E), activated dendritic cells (Figure 7F), monocytes (Figure 7G), activated memory CD4 T cells (Figure 7H), CD8 T cells (Figure 7I), M1 macrophages (Supplementary Figures S6D, S7E), gamma delta T cells (Supplementary Figure S5D), naive B cells (Supplementary Figure S6C), activated NK cells (Supplementary Figure S6F), and memory B cells (Supplementary Figure S7C). From those Venn graphs of meaningful immune cells in the violin plot and correlation graph, common immune cells were selected in GSE14520 datasets (Figure 7J), GSE76427 datasets (Supplementary Figures S5E), ICGC datasets (Supplementary Figure S6G), and TCGA datasets (Supplementary Figure S7F). We found CD8 T cells, activated memory CD4 T cells, monocytes, M0 macrophages, M2 macrophages, and activated dendritic cells as the common immune cells in GSE14520 dataset. M0 macrophages and gamma delta T cells were the common immune cells in GSE76427 datasets. Moreover, activated NK cells and M1 macrophages were common immune cells in ICGC dataset. Finally, M0 macrophages and M1 macrophages were the common immune cells in TCGA dataset.

In Conclusion, the aforementioned analysis of immune infiltration in GSE14520 dataset, GSE76427 dataset, ICGC dataset, and TCGA dataset showed that remarkable differences in immune cells existed in HCC patients, and the correlations between the expression of NAP1L1 and some immune cells were statistically significant. Finally, we speculated that NAP1L1 expression is most associated with macrophages in HCC.

The association between the enrichment of TILs and the expression of NAP1L1 were inferred using the TISIDB database, and the correlation between the NAP1L1 expression and the TILs of people’s cancer is displayed in Figure 8A. It was a certain association between the enrichment of 28 TIL types and the expression of NAP1L1. Specifically speaking, NAP1L1 expression was closely positively related to the following immune cells: activated B cells (r = 0.104, p = 0.0441), activated CD4 T cells (r = 0.4, p = 0.6.9e-17), activated CD8 T cells (r = 0.266, p = 2e-07), activated dendritic cells (r = 0.104, p = 0.0292), immature B cells (r = 0.105, p = 0.0443), macrophage (r = 0.153, p = 0.00316), mast cells (r = 0.172, p = 0.000864), myeloid-derived suppressor cells (r = 0.163, p = 0.00162), regulatory T cells (r = 0.123, p = 0.0175), central memory CD4 T cells (r = 0.215, p = 3.07e-05), effector memory CD4 T cells (r = 0.148, p = 0.004), T follicular helper cells (r = 0.293, p = 9.33e-09), gamma delta T cells (r = 0.332, p = 6.02e-11), and Type-2 T helper cells (r = 0.186, p = 0.000309) (Figure 8B-O). In addition, the expression of NAP1L1 was deeply negatively relevant to the following immune cells: natural killer cells (r = −c0.207, p = 5.86e-05) and Type-17 helper cells (r = −0.147, p = 0.00455) (Figures 8P, Q). Last, the TIMER database was utilized to further verify the relations between the NAP1L1 expression and the immune cell infiltration in HCC (Figures 8R, S). The outcomes displayed that the NAP1L1 expression was closely relevant to the following immune cells: B cells (p = 3.58e-16), CD8⁺ T cells (p = 1.69e–14), CD4⁺ T cells (p = 2.59e–12), macrophages (p = 6.12e–29), neutrophils (p = 4.94e–16), and dendritic cells (p = 2.37e–24).

To summarize, through the comprehensive analysis of GSE14520, GSE76427, ICGC, TCGA, TISIDB, and TIMER databases, we had a hypothesis that NAP1L1 expression was related to immune cell and NAP1L1 may also accelerate the HCC progression by affecting macrophage through some underlying mechanisms.