Three NAFLD-related datasets GSE89632 (Arendt et al., 2015), GSE49541 (Moylan et al., 2014) and GSE164441 (Yang et al., 2021) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) including the gene expression profile data and related annotation file. The search keywords used were “Nonalcoholic fatty liver disease” AND “Homo sapiens.” The inclusion criteria for the dataset are as follows: (1) Each dataset included contained at least 10 samples; (2) Each sample contained pathological information; (3) The raw data or gene expression profiling from these datasets are accessible at the GEO.

The details of these datasets were shown in Table 1. The platform for GSE89632 was Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip, for GSE49541 was Affymetrix Human Genome U133 Plus 2.0 Array (HG-U133_Plus_2), and for GSE164441 was GPL20301, Illumina HiSeq 4000 (Homo sapiens). It is important to note that samples with no specific pathologic information were not analyzed.

The “limma” software package in R software was used to screen out DEGs based on the cut-off criteria (Ritchie et al., 2015). The conditions of |log2 FC| ≥ 0.5, adjust p-value < 0.05 (corrected by B and H method) was used for NAFLD without fibrosis vs. healthy liver sample and NAFLD with fibrosis vs. healthy liver sample group to obtain DEGs associated with NAFLD fibrosis. GSE49541 included 40 mild NAFLD liver tissue samples and 32 advanced NAFLD liver tissue samples. The cut-off criterion of adjust p-value < 0.05 (corrected by B and H method) was used to obtain DEGs associated with NAFLD fibrosis progression. The heatmap for the DEGs was created using the “ggplot2” package (https://ggplot2.tidyverse.org), and the Venn plot for co-expressed DEGs was integrated using the “ggvenn” (https://github.com/yanlinlin82/ggvenn) packages, and utilized for the following analysis.

To investigate the functional biological roles of the co-expressed DEGs, Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were automatically completed and visualized by the Metascape (Zhou et al., 2019). p-value < 0.05 was set as the significance threshold.

The common significant DEGs were mapped to the STRING database (http://www. string-db. org/), then set confidence >0.4 to conduct the PPI network analysis (Szklarczyk et al., 2019), which were visualized by Cytoscape software (v3.8.1) (Shannon et al., 2003). Module analysis was performed to identify significant modules and explore the hub genes using the Molecular Complex Detection (MCODE V2.0.0) (Bader and Hogue, 2003), which is a plug-in of Cytoscape (MCODE score ≥3.5).

First, we downloaded the gene expression matrix from GSE164441, composed of 10 cases with (NAFLD)-associated HCC tumor and paired adjacent non-tumor liver tissues. The “ggolot2” package visualized the expression of key genes in NAFLD-induced tumor tissues and adjacent non-tumor liver tissues. Clinical survival analyses were retrieved online using the online tools, Gene Expression Profiling Interactive Analysis (GEPIA2) (http://gepia.cancer-pku.cn/) and Kaplan-Meier Plotter (http://kmplot.com/analysis/). p-value = 0.05 and |log2FC| = 1 and were defined as the threshold for statistical significance (Tang et al., 2017). Finally, the HPA database (https://www.proteinatlas.org/) provides tissue and cellular distribution information of various human proteins using immunoassay technology (Uhlén et al., 2005; Berglund et al., 2008; Pontén et al., 2008; Uhlén et al., 2015; Thul et al., 2017; Uhlen et al., 2017).

The human hepatocyte cell lines HL-7702 and hepatoma cell line MHCC-97H (iCell Bioscience Inc., Shanghai, China) were cultured in a humidified incubator at 37 °C with 5% CO2 using RPMI-1640 with 10% FBS and 1% P/S solution (Procell). To establish the in vitro NAFLD cell model, HL-7702 and MHCC-97H cells were cultured in the presence or absence of 1 mmol/L free fatty acids (FFA, containing oleic acid and palmitic acid at a 2:1 volume ratio) for 24 h and then used for the indicated assays (Long et al., 2019).

Ten male C57/BL6 mice (7–8 weeks old) were randomly divided into two groups, one was administrated with a 60% high-fat diet (HFD) (New Brunswick, NJ08901, United States) (n = 5) and the other with a control diet (CD, 10% fat) for 12 weeks (n = 5) (Jiang et al., 2019). The animals were housed in a temperature-controlled environment, with a 12-h light/dark cycle and received food and water ad libitum, at the center of experimental animals of Benxi Laboratory, Shengjing Hospital of China Medical University. The animal experiments received ethical clearance from the Ethics Committee of Shengjing Hospital of China Medical University (2021PS570K).

After 12 weeks’ treatment, the mice were sacrificed under deeply anesthesia. The livers were then perfused with PBS via the portal vein to remove the blood and divided into two pieces. One was preserved in 10% neutral-buffered formalin, while the other was immediately frozen in liquid N₂ and stored at −80°C.

Trizol reagent (Invitrogen, CA, United States) was utilized to extract total RNA according to the manufacturer’s protocol. Total RNA was then reverse-transcribed into the complementary DNAs (cDNAs) with the PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China) and amplified by GoTaq qPCR Master Mix (Promega, Madison, WI, United States) using the ABI ViiA 7 Real-time PCR system (Applied Biosystems, United States). The relative mRNA levels were normalized to those of the housekeeping gene GAPDH. All samples were determined by the 2^–∆∆CT method. The primer pairs used in this study are listed in Supplementary Table S1 .

The NUSAP1 siRNA and negative controls were purchased from HANBIO (HANBIO, Shanghai, China). The transfection of Hl-7702 and MHCC-97H cells was conducted in accordance with manufacturer’s guidance.

Cells transfected with si-NUSAP1 were plated in 96-well plates and cultured for 24 h. Then, the transfected cells were add with 10 μL CCK-8 reagent. After 2 h incubation, the well’s absorbance at 450 nm was measured by microplate reader. Each group set was repeated for five times and each measurement was performed three times repeatedly.

Fresh or frozen tissues were fixed with 4% paraformaldehyde for 1 h at room temperature, and cell samples were washed with PBS and then fixed with 4% paraformaldehyde for 20min. Oil red O stain (Servicebio, G1015-100ML, Yikeshengwu, China) were operated in accordance with the manufacturer’s instructions.

Detailed steps of Western blotting are described (Zhou et al., 2010). Primary antibodies was NUSAP1 polyclonal antibody (ProteinTech, 12024-1-AP, United States) and GAPDH (ProteinTech Group, Rosemont, IL, United States) epitopes monoclonal antibodies. Then, the secondary antibody for Western blotting was the horseradish-peroxidase-labeled secondary antibody (ProteinTech Group, Rosemont, IL, United States).

All the contents of bioinformatics analysis were analyzed with R version 4.0.5 (https://www.r-project.org/). Student’s t-test was used to analyze the relative expression levels of hub genes. The p < 0.05 was considered statistically significant. *p < 0.05, **p < 0.01; ***p < 0.001; ns, not significant.

The bioinformatics analysis workflow is summarized in Figure 1. DEGs related to fibrosis, and advanced fibrosis were identified by integrated bioinformatics analysis. We first determined the DEGs among different groups of GSE89632, which consisted of 18 cases with fibrosis, 21 cases with non-fibrosis and 11 healthy controls (HC). Total number of 5510 DEGs were determined for the non-fibrosis vs. HC group and a total of 3913 DEGs were identified from the fibrosis vs. HC group. Then, we analyzed the datasets GSE49541 (32 cases with advanced fibrosis and 40 cases with mild fibrosis) and identified 739 DEGs associated with advanced fibrosis. The list of differential genes in each group were visualized by the heatmaps whose hierarchical clustering was based on Euclidean distance (Figure 2).

The common DEGs were showed using the “ggvenn” R package, which were significantly associated with progression of NAFLD from non-fibrosis to advanced fibrosis ( Figure 3A). More specific information about the DEGs was shown in the additional file (Supplementary Table S2). To gain insight into biological significance of these genes, we carried out enrichment analysis by using the Metascape online database. The results showed that pathways such as glucocorticoid receptor pathway, regulation of transmembrane transporter activity, peroxisome, proteoglycan biosynthetic process, regulation of wound healing, regulation of chemotaxis, etc., were the significantly enriched GO or pathways (Figure 3B).

To uncover the potential relationship between 112 common DEGs, a PPI network that included 112 nodes and 65 edges were constructed using the STRING tool (Figure 4A). Then, the most significant modules were recognized by the MCODE plug-in of cytoscape. The top six hub genes in the MCODE cluster1 were selected (Table 2), including Kinesin Family Member 22 (KIF22), ZW10 interacting kinetochore protein (ZWINT), nucleolar and spindle associated protein 1 (NUSAP1), proliferating cell nuclear antigen-binding factor (KIAA0101), ubiquitin-like with PHD and ring finger domains 1 (UHRF1), RAD51-associated protein 1 (RAD51AP1) (Figure 4B). ZWINT, NUSAP1 and RAD51AP1 expression showed an increasing trend in the progression of NAFLD (Figure 4C ), but KIF22, KIAA0101 and UHRF1 exhibits a different expression profile in different stages of the NAFLD ( Supplementary Figure S1 ).

To further confirm the role of the hub genes in the progression of NAFLD, experimental validation was conducted in vitro and vivo. The expression of ZWINT, NUSAP1, RAD51AP1 were highly expressed in HL-7702 cells treated with FFA for 24 h as mentioned above by RT-PCR screening (Figures 5A,B). Due to a high homology between mouse and human genes, we also used mouse models of NAFLD established by HFD feeding for 12 weeks to further examine the hub gene expression levels. H&E and Oil red O staining of liver tissues exhibited compressed liver sinusoids, more dispersed lipid vacuoles and increased hepatocyte volumes in NAFLD mice than those in control mice (Figure 5C). The expression of Zwint, and Nusap1 were significantly upregulated in NAFLD mice compared to the control group by RT-PCR, which was consistent with the results of our previous analysis. In contrast to our previous results, Kiaa0101 was significantly increased in NAFLD mouse. The expression levels of Rad51ap1 and Uhrf1 tended to increase in NAFLD mice but the difference was not statistically significant (p = 0.09 and p = 0.054) (Figure 5D). Therefore, NUSAP1 and ZWINT may be the key genes in the progression of NAFLD.

To further inquire the expression of NUSAP1 and ZWINT in the (NAFLD)-associated HCC patients, we downloaded the data set GSE164441 including 10 paired NAFLD-associated HCC tumor and adjacent normal tissues. It was found that the expression of NUSAP1 and ZWINT was significant higher (p < 0.05) in tumours specimens compared with para-cancer tissues ( Figure 6A). We then confirmed this result in in vitro cell line studies. NUSAP1 was highly expressed in MHCC-97H cell lines treated with 1 mmol/L of FFA for 24 h (Figures 6B,C). Then, we also investigated the protein expression of NUSAP1. Similarly, protein levels of NUSAP1 in liver normal cell HL-7702, liver cancer cell MHCC-97H and NAFLD mice liver significantly elevated (Figure 6D). What’s more, IHC results about the NUSAP1 expression, obtained from the HPA database, showed that NUSAP1 had a higher intensity in HCC tissues than that in normal liver tissues (Figure 6E). In order to further verify the functional characteristics of NUSAP1, we conducted the NUSAP1 silencing assay. We firstly confirmed the efficiency of si-NUSAP1 (Figure 7A). Then CCK-8 analysis and Wound healing assays indicated that knockdown of the NUSAP1 gene could result in a significant reduction in migration of MHCC-97H cell under high fat condition (Figures 7B,C). In addition, the lipid content of the cells reduced significantly after NUSAP1 silence (Figure 7D), suggesting that NUSAP1 may be associated with lipid accumulation in this study.

To further understand the role of NUSAP1 expression in HCC patients, the prognostic value of NUSAP1 was analyzed on the GEPIA and Kaplan-Meier Plotter website. The expression of NUSAP1 was significantly higher in liver cancer tissues than that in normal tissues (Figure 8A) and associated with tumor staging I-III but not stage IV (Figure 8B). In addition, the Kaplan-Meier curve and log-rank test analyses revealed that the increased NUSAP1 mRNA levels corresponded to the shortened OS (overall survival), RFS (relapse free survival), PFS (progression free survival) and DFS (distance free survival) in HCC (p < 0.05) (Figures 8C–G).