NOP16 promotes hepatocellular carcinoma progression and triggers EMT through the Keap1-Nrf2 signaling pathway

2.1Public database

TCGA database: We used the TCGA (https://portal.gdc.cancer.gov/) website to download the LIHC transcriptomic data and the corresponding clinical information. The TCGA database contains Illumina RNA-seq gene expression data, including tumour grade and survival data, expression profile data and clinical data. A total of 424 LIHC patients were screened.

HPA database: We performed immunohistochemical experiments using a NOP16 histochemical antibody on normal hepatic tissue and liver cancer tissue. We entered the keyword “NOP16” in the Human Protein Atlas database (http://www.proteinatlas.org/) and selected “Pathology” and “liver cancer”.

2.2Analysis of the differentially expressed genes

NOP16 (cut-off value of 50%) was divided into low expression and high expression groups, and the RNAseq data in the level 3 HTSeq-Counts format in the project were integrated. The R package DESeq2 (version 1.26.0) was used to identify DEGs, and the heatmap was visualised using ggplot2 (version 3.3.3). We performed an unpaired analysis of NOP16 mRNA expression in 33 cancers in the TCGA and GTEx databases, including liver cancer.

2.3Functional enrichment analysis

Differentially expressed genes in the NOP16 high and low expression groups to explore the effect of differential NOP16 expression levels on cell signaling pathways. In this study, the KEGG.v7.4 data set was used for enrichment analysis according to the default reference settings, and q< 0.05 was selected as significantly enriched gene sets. The “ggplot2” package was used to analyze the biological processes, molecular functions, and cellular localization involved in the NOP16 co-expressed genes, and the KEGG pathway enrichment analysis to discover the signaling pathways that these genes are involved in COAD.

2.4Cell culture and transfection

Cells were seeded into 6-well cell culture plates 1 d before transfection at 6 × 105 cells per well. At 80% adherent cell growth, transfection was performed according to the Lipogectamine TM3000 transfection kit instructions. Negative control cells were transfected with negative control siRNA, siRNA-NOP16 group cells were transfected with NOP16 siRNA, and blank control cells were not transfected. After 8 h of transfection, the normal medium (containing 10% FBS) was replaced. Cell transfection efficiency was determined using qRT-PCR and Western blotting (NOP16 mRNA and protein expression levels in cells, respectively).

2.5Apoptosis and the cell cycle measurement using flow cytometry

For conventional culture-negative control and NOP16-siRNA cells at 70% confluence, a single cell suspension was collected without EDTA pancreatic enzyme digestion or centrifugation, and the cells were washed twice via PBS centrifugation. The cell apoptosis kit and cell cycle kit were used, and the cell apoptosis rate and cell cycle of each group were determined using flow cytometry.

2.6Cell wound scratch assay

Cells of the negative control and NOP16-siRNA groups were seeded in 6-well plates and cultured at 37∘C in a 5% CO2 incubator on the cell surface. Cell fragments were removed by washing with PBS, and serum-free culture medium was removed from the same field at 0 h, 24 h, and 48 h. Cell migration was calculated.

2.7Transwell assays

The treated cells were washed twice with PBS and treated with serum-free DMEM starvation for 12 h. The treated cells were collected using trypsin digestion, and a single-cell suspension was prepared and seeded in the Transwell upper chamber (8 μm, 24-well insert; Transwell; Corning Inc. Lowell, MA, USA) at 5 × 103 cells per well). Complete culture medium (600 μL) was added to the lower chamber for 24 h. The chambers were removed, and unmigrated cells were discarded, fixed in 4% paraformaldehyde for 10 min, washed twice in PBS, dried, stained with crystal violet and placed under a light microscope for counting.

2.8Western blotting

Total intracellular protein was extracted using RIPA lysis buffer. The total protein concentration was determined using BCA. SDS-PAGE loading buffer was added for 10 min using 10% PAGE gel electrophoresis. Isolated proteins were transferred onto PVDF membranes and blocked with 5% skim milk powder for 1 h. The corresponding primary antibody and β-actin antibody (1:1000 dilution) were added and incubated at 4∘C overnight. After incubation with the secondary antibodies, the relative expression levels of the target proteins were analysed and compared after three washes with TBST buffer.

2.9Real-time qRT-PCR

Total intracellular RNA was extracted using the TRIzol method, and first-strand cDNA was synthesised using reverse transcription and used as the template for qPCR. Each reaction system included components according to the PCR kit instructions. The qPCR procedure was 95∘C pre-denaturation for 2 min, 95∘C denaturation for 15 s and 60∘C annealing for 1 min for a total of 40 cycles. Ct values were calculated for each group and compared between groups using the 2-Δ⁢ΔCt method. Primer sequences and antibody information are shown in Table 1.

2.10Prognostic analysis

Clinical information of LIHC patients with effective prognostic information included in the TCGA database was screened. The independent prognosis and risk score and clinical characteristics (age, sex and disease classification) by R software to plot the forest map, observe the accuracy of the survival of patients by drawing ROC curve, verify the results of independent prognostic analysis, and draw the nomogram to predict the patient survival rate

2.11Construction and validation of the nomograms

To predict the overall survival probability, we constructed a nomogram based on the independent prognostic factors in the multivariate Cox analysis. The nomogram was drawn with 1000 equal repeated samples of Bootstrap for 1 year, 3 year and 5 year calibration curves.

2.12Single-cell RNA-seq

The scRNA-seq data (GSE112271) [42] of LIHC patients were obtained from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/, accessed November 10, 2022). Single-cell RNA sequencing (scRNA seq) data were processed using the R package “Seurat” then further standardised data and multiple books were integrated. “Louvain” clustering was performed with 23 principal components. Some classical cell subpopulation-defined markers were obtained from previous studies and manually annotated according to marker expression.

2.13Statistical analysis

The statistical significance of cell proliferation, migration, cycle, apoptosis, invasion and gene expression was determined using Student’s t test (for comparisons between two groups) or two-way ANOVA (for comparisons between multiple groups). The relevant data are shown as fold-changes (mean ± SEM or ± SD), and significant differences were calculated at p< 0.05.

3.1The expression level of NOP16 was increased in LIHC

Pan-cancer analysis showed that NOP16 was highly expressed in 20 cancer types, including LIHC, bladder and urothelial carcinoma, breast cancer, cervical squamous cell carcinoma, cholangiocarcinoma, colon cancer, glioblastoma, oesophageal cancer, kidney chromophobe, kidney renal clear cell carcinoma, lung adenocarcinoma, lung squamous cell carcinoma, pheochromocytoma and paraganglioma, prostate adenocarcinoma, rectum adenocarcinoma, stomach adenocarcinoma, thyroid carcinoma and uterine corpus endometrial carcinoma (Fig. 1A). NOP16 expression was significantly higher in LIHC samples than normal tissues, and NOP16 expression was more significant in paired sample tumours than normal samples (p< 0.001) (Fig. 1B and C). The ROC curve showed that NOP16 expression was also predictive, with an area under the curve (AUC) of 0.826 (95% CI [CI] = 0.770–0.882) (Fig. 1D).

Figure 1.

Expression levels of NOP16 in different types of tumours and LIHC. Expression of NOP16 (A) in different types of tumours compared to normal tissues in TCGA and GTEx databases, (B) in LIHC and non-matched normal tissues in the TCGA and GTEx databases, (C) in LIHC and matched normal tissues in the TCGA and GTEx databases. (D) ROC curves for classifying LIHC versus normal tissues in the TCGA database.

Figure 2.

Representative images of NOP16 expression in LIHC tissues and normal tissues. Images of immunohistochemistry staining for LIHC and normal tissues were collected with HPA. The greater the antigen content (representing the level of protein expression) and the higher the distribution density, the stronger the positive result color rendering. According to the degree of colorrendering of positive markers, they are classified as: blue, negative; light yellow, weaklypositive; brown, moderately positive; and dark brown, strongly positive.

To further identify the expression of NOP16 at the LIHC protein level, we used the HPA database. The NOP16 histochemical antibody (No. HPA036506) showed primary staining in the nucleus. The results showed that the positive expression area and positive intensity of NOP16 in LIHC tissues (n= 12) were significantly higher than normal tissues (n= 3) (Fig. 2). Case statistics for immunohistochemistry are shown in Table 1. This result was consistent with the pan-cancer analysis (Fig. 1).

3.2Association between NOP16 expression and clinicopathological variables

The 375 primary tumors collected from TCGA included clinical features and gene expression data. Patient characteristics included gender, race, differentiation, TNM stage, pathological stage, initial treatment, histological type, tumor anatomical site division, TP 53 and PIK3CA status. The analysis showed that high expression of NOP16 was significantly associated with gender (p= 0.002), histological grade (p= 0.033) and albumin (g/dl) (p= 0.035), as shown in Table 1. The expression of NOP16 in males was significantly higher than females (Fig. 1A), Nevertheless, NOP16 remains higher in female patients than in normal tissue (Fig. 1B). In addition, from the concise G grade reflecting tumor tissue atypia, the expression level of NOP 14 was significantly higher (P< 0.001) than G1 (Fig. 1C).

3.3Identification of differentially expressed genes (DEGs) in LIHC and functional enrichment analysis, including GO, KEGG, and GSEA

There were 1580 significantly different genes, which were used to analyse the differential expression of NOP16 (high expression group and low expression group). A total of 758 (47.9%) genes were upregulated, and 826 (52.1%) genes were downregulated (adjusted p value <0.05, |Log2-FC| > 1) (Fig. 3A). Immune-related DEGS (PDCD 1, CD274, CTLA 4, LAG 3, HAVCR2, TIGIT, CD48) are shown in Fig. 3B.

Figure 3.

NOP16-related differentially expressed genes (DEGs) and functional enrichment analysis of NOP16 in glioma using GO and KEGG. (A) Volcano plot of DEGs. Blue and red dots indicate the significantly downregulated and upregulated DEGs, respectively. (B) Heatmap of the coexpression of NOP16 and immune-related genes. (C) GO and KEGG analysis of DEGs. (GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes. p*< 0.05, p**< 0.01, and p*⁣**< 0.001.)

GO enrichment analysis, including biological processes, cellular compositions, and molecular functions, revealed that DEGs were enriched in different GO terms, such as extracellular matrix disassembly, chromatin assembly, nucleosome, protein-DNA complex, ion gated channel activity, and gated channel activity. KEGG pathway analysis showed that significantly DEG-enriched pathways included neuroactive ligand-receptor interaction and steroid hormone biosynthesis (Fig. 3C and Table 2). GSEA enrichment analysis between the NOP16 high and low expression groups revealed a significant enrichment of biological processes, such as the cell cycle and reactive oxygen species metabolism process, which suggest that NOP16 regulates LIHC occurrence and development via the oxidative stress pathway (Fig. 1B).

3.4Effect of NOP16-siRNA interference on the Keap-Nrf2 signalling pathway in HepG2 cells

To verify the role of NOP16 in the occurrence and development of HCC, we designed three siRNAs to inhibit NOP16 expression in HepG2 cells. The interference efficiency of NOP16 was verified using real-time q-PCR and Western blotting (WB). The results showed that interference with siRNA-3 increased the NOP16 interference efficiency greater than 50% in HepG2 cells (Fig. 4A and B). Therefore, we selected siRNA-3 for the subsequent cellular function experiments.

Figure 4.

The decrease in NOP16 leads to abnormal expression of the Keap1-Nrf2 signalling pathway. (A–B) siRNA interfered with the expression of NOP16 in HepG2 cells, and mRNA and protein levels were measured using q-PCR (A) and WB, respectively (B). (C–D) siRNA3 interfered with the expression of NOP16 and affected the Keap1-Nrf2 pathway. q-PCR (C) and WB (D) were used to detect its expression. The full name of the protein involved in Fig D is as follows: Keap1 (kelch like ECH associated protein 1), Nrf2 NRF2 (Nuclear Factor erythroid 2-Related Factor 2), GCLC (Recombinant Glutamate Cysteine Ligase, Catalytic), NQO1 (NAD (P) H: quinoneoxidoreductaseNQO1), GCLM (Recombinant Glutamate Cysteine Ligase, Modifier Subunit), HO-1 (NAD (P) H: quinoneoxidoreductaseNQO1) N⩾ 3, and “**” “$⁡$” indicates p⩽ 0.01.

Keapl-Nrf2 is the most important endogenous antioxidative stress pathway discovered in recent years [1]. Nrf2 is a key regulator of the oxidant/antioxidant balance and is primarily activated in response to oxidative stress. Because of the enrichment of NOP16 in the oxidative stress pathway, q-PCR and WB were used to detect key genes in the Keap1-Nrf2 signalling pathway. The results showed that Keap1 decreased significantly compared to negative results, and Nrf2, GCLC, NQO1, GCLM and H0-1 were significantly higher in the NOP16-siRNA-3 group than the negative controls (P< 0.05). Similar results were observed in the Western blot assays. These results indicate that the knockdown of the NOP16 gene activated the Keap1-Nrf2 pathway to exert its hepatoprotective effects [16] (Fig. 4C and D).

3.5Functional validation of NOP16 in HepG2 cells

The decreased expression of NOP16 caused abnormal expression of the Keap1-Nrf2 pathway (Fig. 4C and D). Abnormal expression of the Keap1-Nrf2 pathway leads to abnormalities in cell infection, migration, the cell cycle and apoptosis [15]. To define the biological function of NOP16 in HepG2 cells, Transwell and cell scratch experiments were used to examine the effect of NOP16 on cell migration. The Transwell experiments showed that decreased NOP16 expression reduced cell invasiveness (p< 0.001) (Fig. 5A and B). However, the cell scratch experiments showed that interfering with NOP16 inhibited the cell scratch healing rate in HepG2 cells (p< 0.01) (Fig. 5C–E). The effect of NOP16 on the HepG2 cell cycle and apoptosis was determined using flow cytometry, and the results showed that interfering with NOP16 arrested the cell cycle in the G1 phase (p< 0.001) (Fig. 5F and G). However, the apoptosis results showed that the silencing of NOP16 increased the apoptosis rate in HepG2 cells compared to the control group (p< 0.001) (Fig. 5H and I). Therefore, the results suggest that NOP16 may be involved in invasion, migration, the cell cycle and cell apoptosis via the Keap1-Nrf2 pathway.

Figure 5.

Analyses of cell migration, wound healing, cell cycle and apoptosis. (A) Transwell migration assays. The polycarbonate film stained with crystal violet was photographed on an inverted microscope. The number of cells that migrated across the membrane in NC and NOP16-siRNA-3 cells was calculated. (B) The number of infected cells in each NC and NOP16-siRNA-3 cell line was calculated according to the infected cell chart in Fig. A. (C–E) Wound healing assay. (C) NC and NOP16-siRNA-3 cells covering the bottom of the small dish were scratched with the tip of 10 μL pipette, and the healing of each cell line was recorded at 0 h, 24 h and 48 h using inverted microscopy. Pictures show representative examples. Scale bars, 750 μm. (D) According to the scratch width in Fig. C, the cell migration distance of each cell line was calculated at 24 h and 48 h. (E) According to the scratch, the cell mobility of each cell was calculated at 24 h and 48 h. (F) Cell cycle detected using flow cytometry. Cell number distribution maps of NC and NOP16-siRNA-3 cell lines at each stage of the cell cycle were detected using flow cytometry. (G) The percentage of NC and NOP16-siRNA-3 cell lines in G0/G1, S and G2/M phases was calculated. (H) Cell apoptosis was detected using flow cytometry. Cell apoptosis of NC and NOP16-siRNA-3 cell lines was detected using flow cytometry. The cells were divided into four zones: “Dead”, “Late Apop”, “Early Apop” and “Live”. (I) It was used to calculate the percentage of apoptotic cells in NC and NOP16-siRNA-3 cells. “LA” refers to Late Apop, and “EA” refers to Early Apop. N⩾ 3, and “**” “$⁡$” indicates p⩽ 0.01.

3.6NOP16-siRNA leads to the dysregulation of genes that control the cell cycle and apoptosis

The above results indicated that the silencing of NOP16 led to phenotypic changes in cell behaviour, such as the invasive capacity of the cells, migration, cell cycle, and apoptosis. Therefore, we examined key genes controlling cellular processes and behaviour using real-time q-PCR and WB.

The mRNA expression levels of key genes in EMT progression were determined using q-PCR. The results showed that the relative expression levels of CTNNA1, CTNNB1, and CDH1 were significantly higher in the NOP16-siRNA group than the negative control group. However, the relative expression levels of SNAl1, SNAl2, CDH2, FNT1, MMP17, MMP9, and vimentin were significantly higher than the negative controls (Fig. 6A) (P< 0.05). Similar results were observed using Western blot analysis compared to a negative control. The protein expression levels of SNA12, CDH2, FN1, MMP9 and vimentin were significantly reduced in the cells of the NOP16-siRNA group (Fig. 6B). NOP16 knockdown significantly reversed the EMT process in HepG2 cells.

Figure 6.

Decreased NOP16 results in abnormal expression of genes typical for EMT, the cell cycle and apoptosis. (A) Q-PCR was used to detect the mRNA expression of NOP16, CTNNA1, CTNNB1, SNAI1, SNAI2, CDH1, CDH2, FN1, MMP17, MMP9 and vimentin in NC and NOP16-siRNA-3 cells. (B) WB was used to detect the protein expression of NOP16 SNAI2, CDH2, FN1, MMP9 and vimentin in NC and NOP16-siRNA-3 cells. (C) Q-PCR was used to detect the mRNA expression of P53, P21, CDK2, CDK6 and cyclinD1 in NC and NOP16-siRNA-3 cells. (D) Quantitative real-time PCR was used to detect the mRNA expression of CASP3 CASP4 CASP6 BCL2 and BCL2L1 in NC and NOP16-siRNA-3 cells. N⩾ 3, and “**” “$⁡$” indicates p⩽ 0.01.

Figure 7.

Prognostic values of NOP16 expression in patients with LIHC evaluated using the Kaplan-Meier method. Overall survival (A) and disease-specific survival (B) for LIHC patients with high versus low NOP16. (C) Forest map based on multivariate Cox analysis for overall survival. HR, hazard ratio; CI, confidence interval.

Figure 8.

A nomogram and calibration curves for prediction of one-, three-, and five-year overall survival rates of patients with LIHC. (A) A nomogram for the prediction of one-, three-, and five-year overall survival rates of patients with LIHC. (B–D) Calibration curves of the nomogram prediction of one-, three-, and five-year overall survival rates of patients with LIHC.

Figure 9.

Immune cell type-specific expression of NOP16 in LIHC. (A) Proportion of 12 major cell lines; (B) sorted cells in the gland-labelled set; (C–D) expression of NOP16 in different cell subsets; (E) score of the corresponding pathway for each cell subset.

The effect of NOP16 on the mRNA expression levels of cell cycle regulators was further examined. The results showed that the relative cyclinD1 expression level of cyclin-dependent kinase 6 (CDK6) in the NOP16-siRNA group was significantly higher than the negative control group. However, the relative expression levels of p53, p21 and CDK2 were significantly lower than the negative controls (P< 0.05). These results indicate that knockdown of NOP16 caused cell cycle arrest and inhibited cell proliferation (Fig. 6C).

We also examined the mRNA expression levels of apoptosis-related genes, and the results suggested that the relative expression levels of Casp3, Casp4 and Casp6 in the NOP16-siRNA group were significantly higher than the negative control group. However, the expression levels of BCL-2 and BCL2L1 were significantly lower than the negative controls (all Ps < 0.05). The BCL-2 family of proteins control and regulate the apoptotic pathway of the mitochondria by promoting mitochondrial outer membrane permeability (MOMP). MOMP releases proapoptotic factors from the mitochondria into the cytosol to further activate the caspase protease cascade [17, 18]. Therefore, we combined our experimental results to show that knockdown of the NOP16 gene promoted the apoptosis of HepG2 cells. All of these results explain our study of the NOP16 phenotype (Fig. 6D).

3.7The prognostic value of NOP16 in LIHC therapy

The correlation between NOP16 expression and the prognosis of LIHC patients was calculated using the Kaplan-Meier method. Patients were divided into high and low NOP16 expression groups using the median NOP16 expression as the cut-off value. OS and DSS in NOP16 showed significantly poor prognosis (OS: hazard ratio [HR] = 1.46, 95% CI = 1.03–2.07, p= 0.032; DSS: HR = 2.75, 95% CI = 2.11–3.59, p< 0.001) (Fig. 7A and B). Univariate and multivariate Cox regression analyses were used to determine the prognostic measures. The results of the multivariate analysis showed that tumour status and NOP16 were independent factors of OS in patients with LIHC. The prognosis of patients with high NOP16 expression, regardless of OS or DSS, was significantly more unfavourable in several subgroups, including T stage 3, Tumour status: with tumour, Gender: male, Weight: ⩽ 70, Pathological stage: Stage III, BMI: ⩽ 25, Histological grade: G3, Prothrombin time: ⩾ 4 and Fibrosis ishak score: 1/2 (Fig. 7C and Table 2).

3.8Construction and verification of the nomogram based on independent factors

Since disease stage is a golden indicator of prognosis, more accurate results would be obtained if the expression status of NOP 16 and pathological stage were combined to jointly predict the patient outcome. Different clinicopathological factors were assigned to the corresponding cut-off values, adding up to obtain the total score, and finally the 1-, 3-, and 5-year survival rate of each HCC patient were evaluated based on the total score. A higher total score on the nomogram indicates a worse prognosis. The expression level of NOP 16 had relatively little effect on the total score (Fig. 8A). Nomograms were drawn with 1000 isoback replicates by Bootstrap (Fig. 8B and C). The results showed that the bootstrap-corrected c-index of the nomogram was 0.671 (95% CI = 0.638–0.7.4), which means that the survival outcome predicted by the model is consistent with the actual survival outcome.

3.9NOP16 expression correlated with T-lymphocyte infiltration in LIHC

Immune cells play a role in inducing EMT and promoting tumor metastasis. The various immune cells that accumulate in the tumor stroma interact with neighboring cancer cells to activate the previously dormant EMT progression [17]. Therefore, this study further analyzed the subsets of immune cells responsible for the upregulation of NOP 16 expression. To identify cell types expressing NOP16, we used single-cell RNA sequencing data, GSE112271, from all seven hepatocellular carcinoma samples. We used a canonical marker set to classify the cells into 12 major cell lines (Fig. 9A and B), HPLC-like (e.g., tumour-associated hepatic progenitor cells, epithelial cells, liver bud hepatic cells, enterocytes and other immune cells (e.g., NK cells, B cells, CD cells, monocytes, Treg cells, and giant cells) are uncommon cell types. NOP16 was particularly highly expressed in T cells (Fig. 9C and D).

We downloaded the LlHC scRNA sequencing data and analysed the pathways enriched by NOP16. The results suggested that EMT markers and genes upregulated by reactive oxygen species (ROS) were significantly enriched (Fig. 9E). This result is highly consistent with our previous experimental results.